OpenCV中initUndistortRectifyMap函数存在bug原因探究

原文首发于公众号「3D视觉工坊」：OpenCV中initUndistortRectifyMap函数存在bug原因探究。

最近在运行如下一段代码时，生成的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 = iir[1] + ir[2], _y = iir[4] + ir[5], _w = iir[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 = _xw, y = _yw;

double x2 = xx, y2 = yy;

double r2 = x2 + y2, _2xy = 2 * xy;

double kr = (1 + ((k3r2 + k2)r2 + k1)r2) / (1 + ((k6r2 + k5)r2 + k4)r2);

double xd = (xkr + p1_2xy + p2(r2 + 2 * x2) + s1r2 + s2r2r2);

double yd = (ykr + p1(r2 + 2 * y2) + p2*_2xy + s3r2 + s4r2r2);

cv::Vec3d vecTilt = matTiltcv::Vec3d(xd, yd, 1);

double invProj = vecTilt(2) ? 1. / vecTilt(2) : 1;

double u = fxinvProjvecTilt(0) + u0;

double v = fyinvProjvecTilt(1) + v0;

yong.qi added

//double A = ((k3r2 + k2)r2 + k1);

//double B = ((k6r2 + k5)r2 + k4);

//double r2_A = r2A;

//double r2_B = r2B;

//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>(uINTER_TAB_SIZE);

int iv = saturate_cast<int>(vINTER_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。

欢迎大家留言积极讨论，同时我在【3D视觉工坊】星球里也发起了作业，感兴趣的小伙伴欢迎积极参与回答。

上述内容，如有侵犯版权，请联系作者，会自行删文。

荐读

2D、3D视觉技术干货之杂谈

再谈「相机标定」

计算机视觉基本原理——RANSAC

一分钟详解本质矩阵的推导过程

一分钟详解OpenCV之相机标定函数calibrateCamera（）

藏在标定板背后的秘密

回复关键词——1，前往知识星球

OpenCV中initUndistortRectifyMap函数存在bug原因探究相关推荐

opencv mat 修改_OpenCV中initUndistortRectifyMap函数存在bug原因探究
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达 3D视觉工坊的第52篇文章最近在运行如下一段代码时,生成的mapx和mapy有点异常. 代码片段如下 ...
imfilter c语言,opencv中cvFilter2D( ) 函数filter2D（）函数与MATLAB中imfilter（）函数的差异...
出处: 1:cvFilter2D() 函数为opencv中c语言函数 2:filter2D()函数为opencv中c++函数 3:imfilter()函数为matlab版本函数计算结果的异同: 2( ...
Opencv中直方图函数calcHist
Opencv中直方图函数calcHist calcHist函数在Opencv中是极难理解的一个函数,一方面是参数说明晦涩难懂,另一方面,说明书给出的实例也不足以令人完全搞清楚该函数的使用方式.最难理解 ...
转载：opencv中imshow函数运行中断的解决办法
OpenCV中 imshow函数运行中断的解决方法羊和咩咩 2017-03-10 16:00:49 5919 收藏最后发布:2017-03-10 16:00:49首发:2017-03-10 16: ...
JS中Promise函数then的奥秘探究
JS中Promise函数then的奥秘探究 Promise概述 Promise对象是CommonJS工作组提出的一种规范,目的是为异步操作提供统一接口. 那么,什么是Promises? 首先,它是一个 ...
关于opencv中cvtcolor函数的code参数用CV_BGR2YUV与CV_BGR2YCrCb的了解
根据Opencv中cvtcolor函数的源码描述,,参数CV_BGR2YCrCb用的系数为[1868, 9617, 4899, 11682, 9241], 而参数CV_BGR2YUV对应的系数只是将C ...
关于Opencv中Filter2D函数的补全方式
目录关于Opencv中Filter2D函数的补全方式环境验证 C++举例 Python举例关于Opencv中Filter2D函数的补全方式环境 OpenCV3.4.16(C++) openc ...
OpenCV中flip函数实现
一.flip()函数原型介绍 void cv::flip(InputArray src,OutputArray dst, int flipCode) 各参数含义 src:输入图像. dst:输出图像. ...
分段线性插值c语言程序_【短道速滑】OpenCV中cvResize函数使用双线性插值缩小图像长宽大小一半时速度飞快（比最近邻还快）之异象解析和自我实现。...
点击上方↑↑↑"OpenCV学堂"关注我作者网名:laviewpbt 是图像处理,算法实现与加速优化方面的大神!其开发的imageshop软件大小只有1MB,却实现了非常丰富与复 ...

OpenCV中initUndistortRectifyMap函数存在bug原因探究

荐读

OpenCV中initUndistortRectifyMap函数存在bug原因探究相关推荐

最新文章

热门文章