opencv mat 修改_OpenCV中initUndistortRectifyMap函数存在bug原因探究
點擊上方“3D視覺工坊”，選擇“星標(biāo)”
干貨第一時間送達
3D視覺工坊的第52篇文章
最近在運行如下一段代碼時，生成的mapx和mapy有點異常。
代碼片段如下：
#include
#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(0, 0) = 274.135594f;
mK.at(1, 1) = 274.733768f;
mK.at(0, 2) = 324.330264f;
mK.at(1, 2) = 231.508914f;
cv::Mat mDistCoef = (Mat_(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更進一步分析，如果統(tǒng)計相鄰兩元素差值的絕對值對于10或者該位置處的像素值低于兩邊或者高于兩邊，得到的mapx和mapy的異常點位置處如下圖：
mapx存在的異常位置分布(白色區(qū)域為異常)
mapy存在的異常位置分布如下圖(白色區(qū)域為異常)
放大了細看，如下圖：
上述中的27.786低于兩邊，上述的最后一行29.453低于左邊同時也低于右邊。
根據(jù)以上，想與大家探討兩個問題。
1)為什么mapx和mapy矩陣會發(fā)生突變？
2) ?如何有效地消除上述產(chǎn)生的突變？
相信學(xué)習(xí)視覺的小伙伴們對于畸變矯正時用到的函數(shù)initUndistortRectifyMap并不陌生，此處翻開OpenCV Documentation官網(wǎng)，查詢了一下對于該函數(shù)的解釋，截圖如下：
將上述的介紹簡單翻譯成中文，如下：
函數(shù)原型：
CV_EXPORTS_W void initUndistortRectifyMap( InputArray cameraMatrix, InputArray distCoeffs,
?????????????????????????? InputArray R, InputArray newCameraMatrix,
?????????????????????????? Size size, int m1type, OutputArray map1, OutputArray map2 );
函數(shù)功能：計算無畸變和修正轉(zhuǎn)換映射。
函數(shù)說明：
?? 這個函數(shù)主要用于計算無畸變和修正轉(zhuǎn)換關(guān)系，為了重映射，將結(jié)果以映射的形式表達。無畸變的圖像看起來就像原始的圖像，就好比這個圖像是用內(nèi)參為newCameraMatrix且無畸變的相機采集得到的。
?? 在單目相機的情況下，newCameraMatrix通常等于cameraMatrix，或者可以通過cv::getOptimalNewCameraMatrix來計算，以便更好地控制縮放。
在雙目立體相機的情況下，newCameraMatrix通常設(shè)置為由cv::stereoRectify計算的P1或P2。
?? 此外，根據(jù)矩陣R，新相機在坐標(biāo)空間中的取向也是不同的。例如，它幫助配準(zhǔn)雙目相機的兩個相機方向，從而使得圖像的極線是水平的，且y坐標(biāo)相同。
?? 該函數(shù)實際上為反向映射算法構(gòu)建映射，供反向映射使用。也就是說，對于在已經(jīng)修正畸變的圖像中的每個像素(u，v)，該函數(shù)計算原來圖像(從相機中獲得的原始圖像)中對應(yīng)的坐標(biāo)系。這個過程是這樣的，見上述OpenCV Documentation中的計算公式。
在雙目相機的例子中，這個函數(shù)被調(diào)用兩次：一次是為了確定每個相機的朝向，經(jīng)過stereoRectify之后，依次調(diào)用cv::stereoCalibrate。但是如果雙目立體相機沒有被標(biāo)定，依然可以使用cv::stereoRectifyUncalibrated直接從單應(yīng)性矩陣H中計算修正變換。對每個相機，函數(shù)計算像素域中的單應(yīng)性矩陣H作為修正變換，而不是3D空間中的旋轉(zhuǎn)矩陣R。R可以通過H 矩陣計算得來：
我們翻出OpenCV3.2.0中關(guān)于OpenCV中的initUndistortRectifyMap函數(shù)源碼，重新命名為一個函數(shù)，代入原工程中，分析存在異常的原因。
首先，我們先看一下initUndistortRectifyMap函數(shù)在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_ R = Mat_::eye(3, 3);
Mat_ A = Mat_(cameraMatrix), Ar;
if (!newCameraMatrix.empty())
Ar = Mat_(newCameraMatrix);
else
Ar = getDefaultNewCameraMatrix(A, size, true);
if (!matR.empty())
R = Mat_(matR);
if (!distCoeffs.empty())
distCoeffs = Mat_(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_ 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 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(i);//指向第i+1行第一個元素指針
float* m2f = map2.empty() ? 0 : map2.ptr(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(u*INTER_TAB_SIZE);
int iv = saturate_cast(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;
}
將上述函數(shù)替換掉OpenCV中的函數(shù)，目的是分析A、B以及r2_A，r2_B，kr_res等變量為何會引起異常。其中，kr_res=(1+r2_A)/(1+r2_B)，分析了四處明顯產(chǎn)生異常的位置，數(shù)據(jù)如下：
經(jīng)過上述分析，留給大家思考：
1)為何會產(chǎn)生跳變呢？
2)如何有效解決跳變呢？
3) 源代碼如何優(yōu)化便可以解決呢？
PS：經(jīng)過測試，OpenCV最新版本4.1.0仍然會出現(xiàn)此bug。
歡迎大家留言積極討論，同時我在【3D視覺工坊】星球里也發(fā)起了作業(yè)，感興趣的小伙伴歡迎積極參與回答。
上述內(nèi)容，如有侵犯版權(quán)，請聯(lián)系作者，會自行刪文。
薦讀
2D、3D視覺技術(shù)干貨之雜談
再談「相機標(biāo)定」
計算機視覺基本原理——RANSAC
一分鐘詳解本質(zhì)矩陣的推導(dǎo)過程
一分鐘詳解OpenCV之相機標(biāo)定函數(shù)calibrateCamera()
藏在標(biāo)定板背后的秘密
回復(fù)關(guān)鍵詞——1，前往知識星球
總結(jié)
以上是生活随笔為你收集整理的opencv mat 修改_OpenCV中initUndistortRectifyMap函数存在bug原因探究的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 本息转存是什么意思
- 下一篇： 农行纪念币预约怎么操作