3D视觉工坊的第52篇文章
最近在运行如下一段代码时，生成的mapx和mapy有点异常。
代码片段如下：
#include<opencv2/imgproc/detail/distortion_model.hpp> #include"opencv.hpp" using namespace std; using namespace cv; int main(int argc, char ** argv) { if (argc < 2) { cout << " if(argc < 2)" << endl; return -1; } string img_name = argv[1]; cv::Mat img = imread(img_name, CV_LOAD_IMAGE_UNCHANGED); if (img.cols != 640) { resize(img, img, cv::Size(640, 480)); } ////存在bug ***有明显突变位置 cv::Mat mK = cv::Mat::eye(3, 3, CV_32F); mK.at<float>(0, 0) = 274.135594f; mK.at<float>(1, 1) = 274.733768f; mK.at<float>(0, 2) = 324.330264f; mK.at<float>(1, 2) = 231.508914f; cv::Mat mDistCoef = (Mat_<float>(8, 1) << 0.195393, -0.113167, -0.000077, -0.000027, -0.003565, 0.529943, -0.133162, -0.022701); Rect validPixROI; double alpha = 0; Size image_size(int(img.cols), int(img.rows)); Size new_image_size(image_size); cout << "image_size" << endl << image_size << endl; Mat mK2 = getOptimalNewCameraMatrix(mK, mDistCoef, image_size, alpha, new_image_size, &validPixROI); Mat mapx, mapy; initUndistortRectifyMap(mK, mDistCoef, Mat(), mK2,new_image_size, CV_32F, mapx, mapy); return 0; }
程序中用到的图片示例如下图所示：
运行程序之后，生成的mapx和mapy，存在的较为明显的异常点位置bug如下图所示。
如果我们对mapx和mapy更进一步分析，如果统计相邻两元素差值的绝对值对于10或者该位置处的像素值低于两边或者高于两边，得到的mapx和mapy的异常点位置处如下图：
mapx存在的异常位置分布（白色区域为异常）
mapy存在的异常位置分布如下图（白色区域为异常）
放大了细看，如下图：
上述中的27.786低于两边，上述的最后一行29.453低于左边同时也低于右边。
根据以上，想与大家探讨两个问题。
1）为什么mapx和mapy矩阵会发生突变？
2) 如何有效地消除上述产生的突变？
相信学习视觉的小伙伴们对于畸变矫正时用到的函数initUndistortRectifyMap并不陌生，此处翻开OpenCV Documentation官网，查询了一下对于该函数的解释，截图如下：
将上述的介绍简单翻译成中文，如下：
函数原型：
CV_EXPORTS_W void initUndistortRectifyMap( InputArray cameraMatrix, InputArray distCoeffs,
InputArray R, InputArray newCameraMatrix,
Size size, int m1type, OutputArray map1, OutputArray map2 );
函数功能：计算无畸变和修正转换映射。
函数说明：
这个函数主要用于计算无畸变和修正转换关系，为了重映射，将结果以映射的形式表达。无畸变的图像看起来就像原始的图像，就好比这个图像是用内参为newCameraMatrix且无畸变的相机采集得到的。
在单目相机的情况下，newCameraMatrix通常等于cameraMatrix，或者可以通过cv::getOptimalNewCameraMatrix来计算，以便更好地控制缩放。
在双目立体相机的情况下，newCameraMatrix通常设置为由cv::stereoRectify计算的P1或P2。
此外，根据矩阵R，新相机在坐标空间中的取向也是不同的。例如，它帮助配准双目相机的两个相机方向，从而使得图像的极线是水平的，且y坐标相同。
该函数实际上为反向映射算法构建映射，供反向映射使用。也就是说，对于在已经修正畸变的图像中的每个像素（u，v），该函数计算原来图像（从相机中获得的原始图像）中对应的坐标系。这个过程是这样的，见上述OpenCV Documentation中的计算公式。
在双目相机的例子中，这个函数被调用两次：一次是为了确定每个相机的朝向，经过stereoRectify之后，依次调用cv::stereoCalibrate。但是如果双目立体相机没有被标定，依然可以使用cv::stereoRectifyUncalibrated直接从单应性矩阵H中计算修正变换。对每个相机，函数计算像素域中的单应性矩阵H作为修正变换，而不是3D空间中的旋转矩阵R。R可以通过H 矩阵计算得来：
我们翻出OpenCV3.2.0中关于OpenCV中的initUndistortRectifyMap函数源码，重新命名为一个函数，代入原工程中，分析存在异常的原因。
首先，我们先看一下initUndistortRectifyMap函数在OpenCV3.2.0版本中的源码（稍作了修改，并添加了一点注释），如下：
void initUndistortRectifyMap(cv::Mat _cameraMatrix, cv::Mat _distCoeffs, cv::Mat _matR, cv::Mat _newCameraMatrix,Size size, int m1type, cv::Mat &_map1, cv::Mat &_map2) { Mat cameraMatrix = _cameraMatrix.clone(), distCoeffs = _distCoeffs.clone(); Mat matR = _matR.clone(), newCameraMatrix = _newCameraMatrix.clone(); if (m1type <= 0) m1type = CV_16SC2; CV_Assert(m1type == CV_16SC2 || m1type == CV_32FC1 || m1type == CV_32FC2); _map1.create(size, m1type); Mat map1 = _map1.clone(), map2; if (m1type != CV_32FC2) { _map2.create(size, m1type == CV_16SC2 ? CV_16UC1 : CV_32FC1); map2 = _map2.clone(); } else _map2.release(); Mat_<double> R = Mat_<double>::eye(3, 3); Mat_<double> A = Mat_<double>(cameraMatrix), Ar; if (!newCameraMatrix.empty()) Ar = Mat_<double>(newCameraMatrix); else Ar = getDefaultNewCameraMatrix(A, size, true); if (!matR.empty()) R = Mat_<double>(matR); if (!distCoeffs.empty()) distCoeffs = Mat_<double>(distCoeffs); else { distCoeffs.create(14, 1, CV_64F); distCoeffs = 0.; } CV_Assert(A.size() == Size(3, 3) && A.size() == R.size()); CV_Assert(Ar.size() == Size(3, 3) || Ar.size() == Size(4, 3)); Mat_<double> iR = (Ar.colRange(0, 3)*R).inv(DECOMP_LU); //LU分解求逆，矩阵求逆共有「LU，cholesky,eig以及SVD」 const double* ir = &iR(0, 0); double u0 = A(0, 2), v0 = A(1, 2); double fx = A(0, 0), fy = A(1, 1); CV_Assert(distCoeffs.size() == Size(1, 4) || distCoeffs.size() == Size(4, 1) || distCoeffs.size() == Size(1, 5) || distCoeffs.size() == Size(5, 1) || distCoeffs.size() == Size(1, 8) || distCoeffs.size() == Size(8, 1) || distCoeffs.size() == Size(1, 12) || distCoeffs.size() == Size(12, 1) || distCoeffs.size() == Size(1, 14) || distCoeffs.size() == Size(14, 1)); if (distCoeffs.rows != 1 && !distCoeffs.isContinuous()) distCoeffs = distCoeffs.t(); const double* const distPtr = distCoeffs.ptr<double>(); double k1 = distPtr[0]; double k2 = distPtr[1]; double p1 = distPtr[2]; double p2 = distPtr[3]; double k3 = distCoeffs.cols distCoeffs.rows - 1 >= 5 ? distPtr[4] : 0.; double k4 = distCoeffs.cols distCoeffs.rows - 1 >= 8 ? distPtr[5] : 0.; double k5 = distCoeffs.cols distCoeffs.rows - 1 >= 8 ? distPtr[6] : 0.; double k6 = distCoeffs.cols distCoeffs.rows - 1 >= 8 ? distPtr[7] : 0.; double s1 = distCoeffs.cols distCoeffs.rows - 1 >= 12 ? distPtr[8] : 0.; double s2 = distCoeffs.cols distCoeffs.rows - 1 >= 12 ? distPtr[9] : 0.; double s3 = distCoeffs.cols distCoeffs.rows - 1 >= 12 ? distPtr[10] : 0.; double s4 = distCoeffs.cols distCoeffs.rows - 1 >= 12 ? distPtr[11] : 0.; double tauX = distCoeffs.cols distCoeffs.rows - 1 >= 14 ? distPtr[12] : 0.; double tauY = distCoeffs.cols distCoeffs.rows - 1 >= 14 ? distPtr[13] : 0.; // Matrix for trapezoidal distortion of tilted image sensor //倾斜图像传感器的梯形畸变矩阵 cv::Matx33d matTilt = cv::Matx33d::eye(); cv::detail::computeTiltProjectionMatrix(tauX, tauY, &matTilt); for (int i = 0; i < size.height; i ) { float* m1f = map1.ptr<float>(i);//指向第i 1行第一个元素指针 float* m2f = map2.empty() ? 0 : map2.ptr<float>(i); short* m1 = (short*)m1f; ushort* m2 = (ushort*)m2f; double _x = i*ir[1] ir[2], _y = i*ir[4] ir[5], _w = i*ir[7] ir[8]; for (int j = 0; j < size.width; j , _x = ir[0], _y = ir[3], _w = ir[6]) { double w = 1. / _w, x = _x*w, y = _y*w; double x2 = x*x, y2 = y*y; double r2 = x2 y2, _2xy = 2 * x*y; double kr = (1 ((k3*r2 k2)*r2 k1)*r2) / (1 ((k6*r2 k5)*r2 k4)*r2); double xd = (x*kr p1*_2xy p2*(r2 2 * x2) s1*r2 s2*r2*r2); double yd = (y*kr p1*(r2 2 * y2) p2*_2xy s3*r2 s4*r2*r2); cv::Vec3d vecTilt = matTilt*cv::Vec3d(xd, yd, 1); double invProj = vecTilt(2) ? 1. / vecTilt(2) : 1; double u = fx*invProj*vecTilt(0) u0; double v = fy*invProj*vecTilt(1) v0; ////yong.qi added //double A = ((k3*r2 k2)*r2 k1); //double B = ((k6*r2 k5)*r2 k4); //double r2_A = r2*A; //double r2_B = r2*B; //double kr_res = (1 r2_A) / (1 r2_B); //fout << "A: " << A << endl // << "B: " << B << endl // << "r2_A: " << r2_A << endl // << "r2_B: " << r2_B << endl // << "kr_res: " << kr_res << endl // << "w:" << w << endl // << "x2:" << x2 << endl // << "y2:" << y2 << endl // << "p1:" << p1 << endl // << "p2:" << p2 << endl // << "s1:" << s1 << endl // << "s2:" << s2 << endl // << "s3:" << s3 << endl // << "s4:" << s4 << endl // << "_2xy:" << _2xy << endl // << "r2:" << r2 << endl // << "kr:" << kr << endl // << "xd:" << xd << endl // << "yd:" << yd << endl // << "fx:" << fx << endl // << "fy:" << fy << endl // << "u0:" << u0 << endl // << "v0:" << v0 << endl // << "matTilt:" << matTilt << endl // << "vecTilt:" << vecTilt << endl // << "invProj:" << invProj << endl // << "vecTilt(0):" << vecTilt(0) << endl // << "vecTilt(1):" << vecTilt(1) << endl << endl // << "u:" << u << endl // << "v:" << v << endl; //fout.close(); ////yong.qi end if (m1type == CV_16SC2) { int iu = saturate_cast<int>(u*INTER_TAB_SIZE); int iv = saturate_cast<int>(v*INTER_TAB_SIZE); m1[j * 2] = (short)(iu >> INTER_BITS); m1[j * 2 1] = (short)(iv >> INTER_BITS); m2[j] = (ushort)((iv & (INTER_TAB_SIZE - 1))*INTER_TAB_SIZE (iu & (INTER_TAB_SIZE - 1))); } else if (m1type == CV_32FC1) { m1f[j] = (float)u; m2f[j] = (float)v; } else { m1f[j * 2] = (float)u; m1f[j * 2 1] = (float)v; } } } _map1 = map1; _map2 = map2; }
将上述函数替换掉OpenCV中的函数，目的是分析A、B以及r2_A，r2_B，kr_res等变量为何会引起异常。其中，kr_res=（1 r2_A）/(1 r2_B)，分析了四处明显产生异常的位置，数据如下：
经过上述分析，留给大家思考：
1）为何会产生跳变呢？
2）如何有效解决跳变呢？
3) 源代码如何优化便可以解决呢？
PS：经过测试，OpenCV最新版本4.1.0仍然会出现此bug。
上述内容，如有侵犯版权，请联系作者，会自行删文。
