opencv 去畸变
opencv迭代去畸变算法
- 函数简介
opencv中函数undistortPoints()用于对图像点坐标进行去畸变,以下为该函数解释:
void undistortPoints(InputArray src, OutputArray dst, InputArray cameraMatrix, InputArray distCoeffs, InputArray R=noArray(), InputArray P=noArray())src-原图像坐标;dst-输出图像坐标;cameraMatrix-相机内参矩阵;distCoeffs-畸变系数,有四种畸变模型,分别含有4,5,8个元素,通常使用具有4/5个参数的模型,如果该向量为NULL,那么设定该图像没有畸变;R-相机坐标系的矫正矩阵(即对相机坐标系的位姿调整,见stereoRectify函数中的Rl,Rr),如果矩阵为空,那么默认使用单位矩阵;P-新的相机矩阵(3x3)或者新的投影矩阵(3x4,包含相机坐标系相对世界坐标系的相对位姿,见stereoRectify函数中的Pl,Pr),如该矩阵为空的话,将设置该矩阵为单位阵。
- 源码分析
void cvUndistortPointsInternal( const CvMat* _src, CvMat* _dst, const CvMat* _cameraMatrix, const CvMat* _distCoeffs, const CvMat* matR, const CvMat* matP, cv::TermCriteria criteria) { // 判断迭代条件是否有效 CV_Assert(criteria.isValid()); // 定义中间变量--A相机内参数组,和matA共享内存;RR-矫正变换数组,和_RR共享内存 // k-畸变系数数组 double A[3][3], RR[3][3], k[14]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; CvMat matA=cvMat(3, 3, CV_64F, A), _Dk; CvMat _RR=cvMat(3, 3, CV_64F, RR); cv::Matx33d invMatTilt = cv::Matx33d::eye(); cv::Matx33d matTilt = cv::Matx33d::eye(); // 检查输入变量是否有效 CV_Assert( CV_IS_MAT(_src) && CV_IS_MAT(_dst) && (_src->rows == 1 || _src->cols == 1) && (_dst->rows == 1 || _dst->cols == 1) && _src->cols + _src->rows - 1 == _dst->rows + _dst->cols - 1 && (CV_MAT_TYPE(_src->type) == CV_32FC2 || CV_MAT_TYPE(_src->type) == CV_64FC2) && (CV_MAT_TYPE(_dst->type) == CV_32FC2 || CV_MAT_TYPE(_dst->type) == CV_64FC2)); CV_Assert( CV_IS_MAT(_cameraMatrix) && _cameraMatrix->rows == 3 && _cameraMatrix->cols == 3 ); cvConvert( _cameraMatrix, &matA );// _cameraMatrix <--> matA / A // 判断输入的畸变系数是否有效 if( _distCoeffs ) { CV_Assert( CV_IS_MAT(_distCoeffs) && (_distCoeffs->rows == 1 || _distCoeffs->cols == 1) && (_distCoeffs->rows*_distCoeffs->cols == 4 || _distCoeffs->rows*_distCoeffs->cols == 5 || _distCoeffs->rows*_distCoeffs->cols == 8 || _distCoeffs->rows*_distCoeffs->cols == 12 || _distCoeffs->rows*_distCoeffs->cols == 14)); _Dk = cvMat( _distCoeffs->rows, _distCoeffs->cols, CV_MAKETYPE(CV_64F,CV_MAT_CN(_distCoeffs->type)), k);// _Dk和数组k共享内存指针 cvConvert( _distCoeffs, &_Dk ); if (k[12] != 0 || k[13] != 0) { cv::detail::computeTiltProjectionMatrix<double>(k[12], k[13], NULL, NULL, NULL, &invMatTilt); cv::detail::computeTiltProjectionMatrix<double>(k[12], k[13], &matTilt, NULL, NULL); } } if( matR ) { CV_Assert( CV_IS_MAT(matR) && matR->rows == 3 && matR->cols == 3 ); cvConvert( matR, &_RR );// matR和_RR共享内存指针 } else cvSetIdentity(&_RR); if( matP ) { double PP[3][3]; CvMat _P3x3, _PP=cvMat(3, 3, CV_64F, PP); CV_Assert( CV_IS_MAT(matP) && matP->rows == 3 && (matP->cols == 3 || matP->cols == 4)); cvConvert( cvGetCols(matP, &_P3x3, 0, 3), &_PP );// _PP和数组PP共享内存指针 cvMatMul( &_PP, &_RR, &_RR );// _RR=_PP*_RR 放在一起计算比较高效 } const CvPoint2D32f* srcf = (const CvPoint2D32f*)_src->data.ptr; const CvPoint2D64f* srcd = (const CvPoint2D64f*)_src->data.ptr; CvPoint2D32f* dstf = (CvPoint2D32f*)_dst->data.ptr; CvPoint2D64f* dstd = (CvPoint2D64f*)_dst->data.ptr; int stype = CV_MAT_TYPE(_src->type); int dtype = CV_MAT_TYPE(_dst->type); int sstep = _src->rows == 1 ? 1 : _src->step/CV_ELEM_SIZE(stype); int dstep = _dst->rows == 1 ? 1 : _dst->step/CV_ELEM_SIZE(dtype); double fx = A[0][0]; double fy = A[1][1]; double ifx = 1./fx; double ify = 1./fy; double cx = A[0][2]; double cy = A[1][2]; int n = _src->rows + _src->cols - 1; // 开始对所有点开始遍历 for( int i = 0; i < n; i++ ) { double x, y, x0 = 0, y0 = 0, u, v; if( stype == CV_32FC2 ) { x = srcf[i*sstep].x; y = srcf[i*sstep].y; } else { x = srcd[i*sstep].x; y = srcd[i*sstep].y; } u = x; v = y; x = (x - cx)*ifx;//转换到归一化图像坐标系(含有畸变) y = (y - cy)*ify; //进行畸变矫正 if( _distCoeffs ) { // compensate tilt distortion--该部分系数用来弥补沙氏镜头畸变?? // 如果不懂也没管,因为普通镜头中没有这些畸变系数 cv::Vec3d vecUntilt = invMatTilt * cv::Vec3d(x, y, 1); double invProj = vecUntilt(2) ? 1./vecUntilt(2) : 1; x0 = x = invProj * vecUntilt(0); y0 = y = invProj * vecUntilt(1); double error = std::numeric_limits<double>::max();// error设定为系统最大值 // compensate distortion iteratively // 迭代去除镜头畸变 // 迭代公式 x′= (x−2p1 xy−p2 (r^2 + 2x^2))∕( 1 + k1*r^2 + k2*r^4 + k3*r^6) // y′= (y−2p2 xy−p1 (r^2 + 2y^2))∕( 1 + k1*r^2 + k2*r^4 + k3*r^6) for( int j = 0; ; j++ ) { if ((criteria.type & cv::TermCriteria::COUNT) && j >= criteria.maxCount)// 迭代最大次数为5次 break; if ((criteria.type & cv::TermCriteria::EPS) && error < criteria.epsilon)// 迭代误差阈值为0.01 break; double r2 = x*x + y*y; double icdist = (1 + ((k[7]*r2 + k[6])*r2 + k[5])*r2)/(1 + ((k[4]*r2 + k[1])*r2 + k[0])*r2); double deltaX = 2*k[2]*x*y + k[3]*(r2 + 2*x*x)+ k[8]*r2+k[9]*r2*r2; double deltaY = k[2]*(r2 + 2*y*y) + 2*k[3]*x*y+ k[10]*r2+k[11]*r2*r2; x = (x0 - deltaX)*icdist; y = (y0 - deltaY)*icdist; // 对当前迭代的坐标加畸变,计算误差error用于判断迭代条件 if(criteria.type & cv::TermCriteria::EPS) { double r4, r6, a1, a2, a3, cdist, icdist2; double xd, yd, xd0, yd0; cv::Vec3d vecTilt; r2 = x*x + y*y; r4 = r2*r2; r6 = r4*r2; a1 = 2*x*y; a2 = r2 + 2*x*x; a3 = r2 + 2*y*y; cdist = 1 + k[0]*r2 + k[1]*r4 + k[4]*r6; icdist2 = 1./(1 + k[5]*r2 + k[6]*r4 + k[7]*r6); xd0 = x*cdist*icdist2 + k[2]*a1 + k[3]*a2 + k[8]*r2+k[9]*r4; yd0 = y*cdist*icdist2 + k[2]*a3 + k[3]*a1 + k[10]*r2+k[11]*r4; vecTilt = matTilt*cv::Vec3d(xd0, yd0, 1); invProj = vecTilt(2) ? 1./vecTilt(2) : 1; xd = invProj * vecTilt(0); yd = invProj * vecTilt(1); double x_proj = xd*fx + cx; double y_proj = yd*fy + cy; error = sqrt( pow(x_proj - u, 2) + pow(y_proj - v, 2) ); } } } // 将坐标从归一化图像坐标系转换到成像平面坐标系 double xx = RR[0][0]*x + RR[0][1]*y + RR[0][2]; double yy = RR[1][0]*x + RR[1][1]*y + RR[1][2]; double ww = 1./(RR[2][0]*x + RR[2][1]*y + RR[2][2]); x = xx*ww; y = yy*ww; if( dtype == CV_32FC2 ) { dstf[i*dstep].x = (float)x; dstf[i*dstep].y = (float)y; } else { dstd[i*dstep].x = x; dstd[i*dstep].y = y; } } }该函数所实现的迭代算法很简单,只要了解相机的畸变模型就可以很容易理解函数的实现过程了。
