added GC_COLOR_GRAD cost function type into opencv_stitching

This commit is contained in:
Alexey Spizhevoy 2011-05-27 11:41:35 +00:00
parent 8e3777676c
commit 56f7e54cce
7 changed files with 218 additions and 63 deletions

View File

@ -45,6 +45,7 @@
using namespace std;
using namespace cv;
using namespace cv::gpu;
Ptr<ExposureCompensator> ExposureCompensator::createDefault(int type)
@ -53,6 +54,8 @@ Ptr<ExposureCompensator> ExposureCompensator::createDefault(int type)
return new NoExposureCompensator();
if (type == OVERLAP)
return new OverlapExposureCompensator();
if (type == SEGMENT)
return new SegmentExposureCompensator();
CV_Error(CV_StsBadArg, "unsupported exposure compensation method");
return NULL;
}
@ -61,13 +64,15 @@ Ptr<ExposureCompensator> ExposureCompensator::createDefault(int type)
void OverlapExposureCompensator::feed(const vector<Point> &corners, const vector<Mat> &images,
const vector<Mat> &masks)
{
CV_Assert(corners.size() == images.size() && images.size() == masks.size());
const int num_images = static_cast<int>(images.size());
Mat_<int> N(num_images, num_images); N.setTo(0);
Mat_<double> I(num_images, num_images); I.setTo(0);
Rect dst_roi = resultRoi(corners, images);
Mat subimg1, subimg2;
Mat_<uchar> submask1, submask2, overlap;
Mat_<uchar> submask1, submask2, intersect;
for (int i = 0; i < num_images; ++i)
{
@ -81,9 +86,9 @@ void OverlapExposureCompensator::feed(const vector<Point> &corners, const vector
submask1 = masks[i](Rect(roi.tl() - corners[i], roi.br() - corners[i]));
submask2 = masks[j](Rect(roi.tl() - corners[j], roi.br() - corners[j]));
overlap = submask1 & submask2;
intersect = submask1 & submask2;
N(i, j) = N(j, i) = countNonZero(overlap);
N(i, j) = N(j, i) = countNonZero(intersect);
double Isum1 = 0, Isum2 = 0;
for (int y = 0; y < roi.height; ++y)
@ -92,7 +97,7 @@ void OverlapExposureCompensator::feed(const vector<Point> &corners, const vector
const Point3_<uchar>* r2 = subimg2.ptr<Point3_<uchar> >(y);
for (int x = 0; x < roi.width; ++x)
{
if (overlap(y, x))
if (intersect(y, x))
{
Isum1 += sqrt(static_cast<double>(sqr(r1[x].x) + sqr(r1[x].y) + sqr(r1[x].z)));
Isum2 += sqrt(static_cast<double>(sqr(r2[x].x) + sqr(r2[x].y) + sqr(r2[x].z)));
@ -122,7 +127,7 @@ void OverlapExposureCompensator::feed(const vector<Point> &corners, const vector
}
}
solve(A, b, gains_, DECOMP_SVD);
solve(A, b, gains_);
}
@ -130,3 +135,14 @@ void OverlapExposureCompensator::apply(int index, Point /*corner*/, Mat &image,
{
image *= gains_(index, 0);
}
void SegmentExposureCompensator::feed(const vector<Point> &/*corners*/, const vector<Mat> &/*images*/,
const vector<Mat> &/*masks*/)
{
}
void SegmentExposureCompensator::apply(int /*index*/, Point /*corner*/, Mat &/*image*/, const Mat &/*mask*/)
{
}

View File

@ -48,7 +48,7 @@
class ExposureCompensator
{
public:
enum { NO, OVERLAP };
enum { NO, OVERLAP, SEGMENT };
static cv::Ptr<ExposureCompensator> createDefault(int type);
virtual void feed(const std::vector<cv::Point> &corners, const std::vector<cv::Mat> &images,
@ -78,4 +78,12 @@ private:
};
class SegmentExposureCompensator : public ExposureCompensator
{
public:
void feed(const std::vector<cv::Point> &corners, const std::vector<cv::Mat> &images,
const std::vector<cv::Mat> &masks);
void apply(int index, cv::Point corner, cv::Mat &image, const cv::Mat &mask);
};
#endif // __OPENCV_EXPOSURE_COMPENSATE_HPP__

View File

@ -71,7 +71,7 @@ void printUsage()
<< "\t[--wavecorrect (no|yes)]\n"
<< "\t[--warp (plane|cylindrical|spherical)]\n"
<< "\t[--exposcomp (no|overlap)]\n"
<< "\t[--seam (no|voronoi|graphcut)]\n"
<< "\t[--seam (no|voronoi|gc_color|gc_colorgrad)]\n"
<< "\t[--blend (no|feather|multiband)]\n"
<< "\t[--numbands <int>]\n"
<< "\t[--output <result_img>]\n\n";
@ -95,7 +95,7 @@ int warp_type = Warper::SPHERICAL;
int expos_comp_type = ExposureCompensator::OVERLAP;
bool user_match_conf = false;
float match_conf = 0.6f;
int seam_find_type = SeamFinder::GRAPH_CUT;
int seam_find_type = SeamFinder::GC_COLOR;
int blend_type = Blender::MULTI_BAND;
int numbands = 5;
string result_name = "result.png";
@ -201,6 +201,8 @@ int parseCmdArgs(int argc, char** argv)
expos_comp_type = ExposureCompensator::NO;
else if (string(argv[i + 1]) == "overlap")
expos_comp_type = ExposureCompensator::OVERLAP;
else if (string(argv[i + 1]) == "segment")
expos_comp_type = ExposureCompensator::SEGMENT;
else
{
cout << "Bad exposure compensation method\n";
@ -214,8 +216,10 @@ int parseCmdArgs(int argc, char** argv)
seam_find_type = SeamFinder::NO;
else if (string(argv[i + 1]) == "voronoi")
seam_find_type = SeamFinder::VORONOI;
else if (string(argv[i + 1]) == "graphcut")
seam_find_type = SeamFinder::GRAPH_CUT;
else if (string(argv[i + 1]) == "gc_color")
seam_find_type = SeamFinder::GC_COLOR;
else if (string(argv[i + 1]) == "gc_colorgrad")
seam_find_type = SeamFinder::GC_COLOR_GRAD;
else
{
cout << "Bad seam finding method\n";

View File

@ -65,8 +65,8 @@ namespace
class CpuSurfFeaturesFinder : public FeaturesFinder
{
public:
inline CpuSurfFeaturesFinder(double hess_thresh, int num_octaves, int num_layers,
int num_octaves_descr, int num_layers_descr)
CpuSurfFeaturesFinder(double hess_thresh, int num_octaves, int num_layers,
int num_octaves_descr, int num_layers_descr)
{
detector_ = new SurfFeatureDetector(hess_thresh, num_octaves, num_layers);
extractor_ = new SurfDescriptorExtractor(num_octaves_descr, num_layers_descr);
@ -80,20 +80,12 @@ namespace
Ptr<DescriptorExtractor> extractor_;
};
void CpuSurfFeaturesFinder::find(const Mat &image, ImageFeatures &features)
{
Mat gray_image;
CV_Assert(image.depth() == CV_8U);
cvtColor(image, gray_image, CV_BGR2GRAY);
detector_->detect(gray_image, features.keypoints);
extractor_->compute(gray_image, features.keypoints, features.descriptors);
}
class GpuSurfFeaturesFinder : public FeaturesFinder
{
public:
inline GpuSurfFeaturesFinder(double hess_thresh, int num_octaves, int num_layers,
int num_octaves_descr, int num_layers_descr)
GpuSurfFeaturesFinder(double hess_thresh, int num_octaves, int num_layers,
int num_octaves_descr, int num_layers_descr)
{
surf_.keypointsRatio = 0.1f;
surf_.hessianThreshold = hess_thresh;
@ -113,6 +105,17 @@ namespace
int num_octaves_descr_, num_layers_descr_;
};
void CpuSurfFeaturesFinder::find(const Mat &image, ImageFeatures &features)
{
Mat gray_image;
CV_Assert(image.depth() == CV_8U);
cvtColor(image, gray_image, CV_BGR2GRAY);
detector_->detect(gray_image, features.keypoints);
extractor_->compute(gray_image, features.keypoints, features.descriptors);
}
void GpuSurfFeaturesFinder::find(const Mat &image, ImageFeatures &features)
{
GpuMat gray_image;
@ -132,7 +135,8 @@ namespace
d_descriptors.download(features.descriptors);
}
}
} // anonymous namespace
SurfFeaturesFinder::SurfFeaturesFinder(bool try_use_gpu, double hess_thresh, int num_octaves, int num_layers,
int num_octaves_descr, int num_layers_descr)

View File

@ -52,8 +52,10 @@ Ptr<SeamFinder> SeamFinder::createDefault(int type)
return new NoSeamFinder();
if (type == VORONOI)
return new VoronoiSeamFinder();
if (type == GRAPH_CUT)
return new GraphCutSeamFinder();
if (type == GC_COLOR)
return new GraphCutSeamFinder(GraphCutSeamFinder::COST_COLOR);
if (type == GC_COLOR_GRAD)
return new GraphCutSeamFinder(GraphCutSeamFinder::COST_COLOR_GRAD);
CV_Error(CV_StsBadArg, "unsupported seam finding method");
return NULL;
}
@ -64,25 +66,33 @@ void PairwiseSeamFinder::find(const vector<Mat> &src, const vector<Point> &corne
{
if (src.size() == 0)
return;
images_ = src;
corners_ = corners;
masks_ = masks;
for (size_t i = 0; i < src.size() - 1; ++i)
{
for (size_t j = i + 1; j < src.size(); ++j)
{
Rect roi;
if (overlapRoi(corners[i], corners[j], src[i].size(), src[j].size(), roi))
findInPair(src[i], src[j], corners[i], corners[j], roi, masks[i], masks[j]);
findInPair(i, j, roi);
}
}
}
void VoronoiSeamFinder::findInPair(const Mat &img1, const Mat &img2, Point tl1, Point tl2,
Rect roi, Mat &mask1, Mat &mask2)
void VoronoiSeamFinder::findInPair(size_t first, size_t second, Rect roi)
{
const int gap = 10;
Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
Mat img1 = images_[first], img2 = images_[second];
Mat mask1 = masks_[first], mask2 = masks_[second];
Point tl1 = corners_[first], tl2 = corners_[second];
// Cut submasks with some gap
for (int y = -gap; y < roi.height + gap; ++y)
{
@ -127,27 +137,62 @@ void VoronoiSeamFinder::findInPair(const Mat &img1, const Mat &img2, Point tl1,
}
class GraphCutSeamFinder::Impl
class GraphCutSeamFinder::Impl : public PairwiseSeamFinder
{
public:
Impl(int cost_type, float terminal_cost, float bad_region_penalty)
: cost_type_(cost_type), terminal_cost_(terminal_cost), bad_region_penalty_(bad_region_penalty) {}
void findInPair(const Mat &img1, const Mat &img2, Point tl1, Point tl2,
Rect roi, Mat &mask1, Mat &mask2);
void find(const vector<Mat> &src, const vector<Point> &corners, vector<Mat> &masks);
void findInPair(size_t first, size_t second, Rect roi);
private:
void setGraphWeightsColor(const Mat &img1, const Mat &img2, const Mat &mask1, const Mat &mask2,
GCGraph<float> &graph);
void setGraphWeightsColor(const Mat &img1, const Mat &img2,
const Mat &mask1, const Mat &mask2, GCGraph<float> &graph);
void setGraphWeightsColorGrad(const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2,
const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2,
GCGraph<float> &graph);
vector<Mat> dx_, dy_;
int cost_type_;
float terminal_cost_;
float bad_region_penalty_;
};
void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &img2, const Mat &mask1, const Mat &mask2,
GCGraph<float> &graph)
void GraphCutSeamFinder::Impl::find(const vector<Mat> &src, const vector<Point> &corners,
vector<Mat> &masks)
{
// Compute gradients
dx_.resize(src.size());
dy_.resize(src.size());
Mat dx, dy;
for (size_t i = 0; i < src.size(); ++i)
{
CV_Assert(src[i].channels() == 3);
Sobel(src[i], dx, CV_32F, 1, 0);
Sobel(src[i], dy, CV_32F, 0, 1);
dx_[i].create(src[i].size(), CV_32F);
dy_[i].create(src[i].size(), CV_32F);
for (int y = 0; y < src[i].rows; ++y)
{
const Point3f* dx_row = dx.ptr<Point3f>(y);
const Point3f* dy_row = dy.ptr<Point3f>(y);
float* dx_row_ = dx_[i].ptr<float>(y);
float* dy_row_ = dy_[i].ptr<float>(y);
for (int x = 0; x < src[i].cols; ++x)
{
dx_row_[x] = normL2(dx_row[x]);
dy_row_[x] = normL2(dy_row[x]);
}
}
}
PairwiseSeamFinder::find(src, corners, masks);
}
void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &img2,
const Mat &mask1, const Mat &mask2, GCGraph<float> &graph)
{
const Size img_size = img1.size();
@ -162,9 +207,8 @@ void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &
}
}
const float weight_eps = 1e-3f;
// Set regular edge weights
const float weight_eps = 1e-3f;
for (int y = 0; y < img_size.height; ++y)
{
for (int x = 0; x < img_size.width; ++x)
@ -175,11 +219,9 @@ void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &
float weight = normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) +
normL2(img1.at<Point3f>(y, x + 1), img2.at<Point3f>(y, x + 1)) +
weight_eps;
if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y, x + 1) ||
!mask2.at<uchar>(y, x) || !mask2.at<uchar>(y, x + 1))
weight += bad_region_penalty_;
graph.addEdges(v, v + 1, weight, weight);
}
if (y < img_size.height - 1)
@ -187,11 +229,9 @@ void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &
float weight = normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) +
normL2(img1.at<Point3f>(y + 1, x), img2.at<Point3f>(y + 1, x)) +
weight_eps;
if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y + 1, x) ||
!mask2.at<uchar>(y, x) || !mask2.at<uchar>(y + 1, x))
weight += bad_region_penalty_;
graph.addEdges(v, v + img_size.width, weight, weight);
}
}
@ -199,15 +239,77 @@ void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &
}
void GraphCutSeamFinder::Impl::findInPair(const Mat &img1, const Mat &img2, Point tl1, Point tl2,
Rect roi, Mat &mask1, Mat &mask2)
void GraphCutSeamFinder::Impl::setGraphWeightsColorGrad(
const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2,
const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2,
GCGraph<float> &graph)
{
const int gap = 10;
const Size img_size = img1.size();
// Set terminal weights
for (int y = 0; y < img_size.height; ++y)
{
for (int x = 0; x < img_size.width; ++x)
{
int v = graph.addVtx();
graph.addTermWeights(v, mask1.at<uchar>(y, x) ? terminal_cost_ : 0.f,
mask2.at<uchar>(y, x) ? terminal_cost_ : 0.f);
}
}
// Set regular edge weights
const float weight_eps = 1e-3f;
for (int y = 0; y < img_size.height; ++y)
{
for (int x = 0; x < img_size.width; ++x)
{
int v = y * img_size.width + x;
if (x < img_size.width - 1)
{
float grad = dx1.at<float>(y, x) + dx1.at<float>(y, x + 1) +
dx2.at<float>(y, x) + dx2.at<float>(y, x + 1) + weight_eps;
float weight = (normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) +
normL2(img1.at<Point3f>(y, x + 1), img2.at<Point3f>(y, x + 1))) / grad +
weight_eps;
if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y, x + 1) ||
!mask2.at<uchar>(y, x) || !mask2.at<uchar>(y, x + 1))
weight += bad_region_penalty_;
graph.addEdges(v, v + 1, weight, weight);
}
if (y < img_size.height - 1)
{
float grad = dy1.at<float>(y, x) + dy1.at<float>(y + 1, x) +
dy2.at<float>(y, x) + dy2.at<float>(y + 1, x) + weight_eps;
float weight = (normL2(img1.at<Point3f>(y, x), img2.at<Point3f>(y, x)) +
normL2(img1.at<Point3f>(y + 1, x), img2.at<Point3f>(y + 1, x))) / grad +
weight_eps;
if (!mask1.at<uchar>(y, x) || !mask1.at<uchar>(y + 1, x) ||
!mask2.at<uchar>(y, x) || !mask2.at<uchar>(y + 1, x))
weight += bad_region_penalty_;
graph.addEdges(v, v + img_size.width, weight, weight);
}
}
}
}
void GraphCutSeamFinder::Impl::findInPair(size_t first, size_t second, Rect roi)
{
Mat img1 = images_[first], img2 = images_[second];
Mat dx1 = dx_[first], dx2 = dx_[second];
Mat dy1 = dy_[first], dy2 = dy_[second];
Mat mask1 = masks_[first], mask2 = masks_[second];
Point tl1 = corners_[first], tl2 = corners_[second];
const int gap = 10;
Mat subimg1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3);
Mat subimg2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3);
Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
Mat subdx1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
Mat subdy1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
Mat subdx2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
Mat subdy2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F);
// Cut subimages and submasks with some gap
for (int y = -gap; y < roi.height + gap; ++y)
@ -220,11 +322,15 @@ void GraphCutSeamFinder::Impl::findInPair(const Mat &img1, const Mat &img2, Poin
{
subimg1.at<Point3f>(y + gap, x + gap) = img1.at<Point3f>(y1, x1);
submask1.at<uchar>(y + gap, x + gap) = mask1.at<uchar>(y1, x1);
subdx1.at<float>(y + gap, x + gap) = dx1.at<float>(y1, x1);
subdy1.at<float>(y + gap, x + gap) = dy1.at<float>(y1, x1);
}
else
{
subimg1.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0);
submask1.at<uchar>(y + gap, x + gap) = 0;
subdx1.at<float>(y + gap, x + gap) = 0.f;
subdy1.at<float>(y + gap, x + gap) = 0.f;
}
int y2 = roi.y - tl2.y + y;
@ -233,11 +339,15 @@ void GraphCutSeamFinder::Impl::findInPair(const Mat &img1, const Mat &img2, Poin
{
subimg2.at<Point3f>(y + gap, x + gap) = img2.at<Point3f>(y2, x2);
submask2.at<uchar>(y + gap, x + gap) = mask2.at<uchar>(y2, x2);
subdx2.at<float>(y + gap, x + gap) = dx2.at<float>(y2, x2);
subdy2.at<float>(y + gap, x + gap) = dy2.at<float>(y2, x2);
}
else
{
subimg2.at<Point3f>(y + gap, x + gap) = Point3f(0, 0, 0);
submask2.at<uchar>(y + gap, x + gap) = 0;
subdx2.at<float>(y + gap, x + gap) = 0.f;
subdy2.at<float>(y + gap, x + gap) = 0.f;
}
}
}
@ -252,6 +362,10 @@ void GraphCutSeamFinder::Impl::findInPair(const Mat &img1, const Mat &img2, Poin
case GraphCutSeamFinder::COST_COLOR:
setGraphWeightsColor(subimg1, subimg2, submask1, submask2, graph);
break;
case GraphCutSeamFinder::COST_COLOR_GRAD:
setGraphWeightsColorGrad(subimg1, subimg2, subdx1, subdx2, subdy1, subdy2,
submask1, submask2, graph);
break;
default:
CV_Error(CV_StsBadArg, "unsupported pixel similarity measure");
}
@ -281,8 +395,8 @@ GraphCutSeamFinder::GraphCutSeamFinder(int cost_type, float terminal_cost, float
: impl_(new Impl(cost_type, terminal_cost, bad_region_penalty)) {}
void GraphCutSeamFinder::findInPair(const Mat &img1, const Mat &img2, Point tl1, Point tl2,
Rect roi, Mat &mask1, Mat &mask2)
void GraphCutSeamFinder::find(const vector<Mat> &src, const vector<Point> &corners,
vector<Mat> &masks)
{
impl_->findInPair(img1, img2, tl1, tl2, roi, mask1, mask2);
impl_->find(src, corners, masks);
}

View File

@ -47,7 +47,7 @@
class SeamFinder
{
public:
enum { NO, VORONOI, GRAPH_CUT };
enum { NO, VORONOI, GC_COLOR, GC_COLOR_GRAD };
static cv::Ptr<SeamFinder> createDefault(int type);
virtual ~SeamFinder() {}
@ -66,34 +66,36 @@ public:
class PairwiseSeamFinder : public SeamFinder
{
public:
void find(const std::vector<cv::Mat> &src, const std::vector<cv::Point> &corners,
std::vector<cv::Mat> &masks);
virtual void find(const std::vector<cv::Mat> &src, const std::vector<cv::Point> &corners,
std::vector<cv::Mat> &masks);
protected:
virtual void findInPair(const cv::Mat &img1, const cv::Mat &img2, cv::Point tl1, cv::Point tl2,
cv::Rect roi, cv::Mat &mask1, cv::Mat &mask2) = 0;
virtual void findInPair(size_t first, size_t second, cv::Rect roi) = 0;
std::vector<cv::Mat> images_;
std::vector<cv::Point> corners_;
std::vector<cv::Mat> masks_;
};
class VoronoiSeamFinder : public PairwiseSeamFinder
{
private:
void findInPair(const cv::Mat &img1, const cv::Mat &img2, cv::Point tl1, cv::Point tl2,
cv::Rect roi, cv::Mat &mask1, cv::Mat &mask2);
void findInPair(size_t first, size_t second, cv::Rect roi);
};
class GraphCutSeamFinder : public PairwiseSeamFinder
class GraphCutSeamFinder : public SeamFinder
{
public:
// TODO add COST_COLOR_GRAD support
enum { COST_COLOR };
GraphCutSeamFinder(int cost_type = COST_COLOR, float terminal_cost = 10000.f,
enum { COST_COLOR, COST_COLOR_GRAD };
GraphCutSeamFinder(int cost_type = COST_COLOR_GRAD, float terminal_cost = 10000.f,
float bad_region_penalty = 1000.f);
private:
void findInPair(const cv::Mat &img1, const cv::Mat &img2, cv::Point tl1, cv::Point tl2,
cv::Rect roi, cv::Mat &mask1, cv::Mat &mask2);
void find(const std::vector<cv::Mat> &src, const std::vector<cv::Point> &corners,
std::vector<cv::Mat> &masks);
private:
class Impl;
cv::Ptr<Impl> impl_;
};

View File

@ -91,10 +91,17 @@ B Graph::walkBreadthFirst(int from, B body) const
//////////////////////////////////////////////////////////////////////////////
// Some auxiliary math functions
static inline
float normL2(const cv::Point3f& a)
{
return a.x * a.x + a.y * a.y + a.z * a.z;
}
static inline
float normL2(const cv::Point3f& a, const cv::Point3f& b)
{
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z);
return normL2(a - b);
}