图像算法中会经常用到摄像机的畸变校正，有必要总结分析OpenCV中畸变校正方法，其中包括普通针孔相机模型和鱼眼相机模型fisheye两种畸变校正方法。
普通相机模型畸变校正函数针对OpenCV中的cv::initUndistortRectifyMap()，鱼眼相机模型畸变校正函数对应OpenCV中的cv::fisheye::initUndistortRectifyMap()。两种方法算出映射Mapx和Mapy后，统一用cv::Remap()函数进行插值得到校正后的图像。
1. FishEye模型的畸变校正。
方便起见，直接贴出OpenCV源码，我在里面加了注释说明。建议参考OpenCV官方文档看畸变模型原理会更清楚：https://docs.opencv.org/3.0-beta/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html#fisheye-initundistortrectifymap
简要流程就是：
1. 求内参矩阵的逆，由于摄像机坐标系的三维点到二维图像平面，需要乘以旋转矩阵R和内参矩阵K。那么反向投影回去则是二维图像坐标乘以 K*R的逆矩阵。
2. 将目标图像中的每一个像素点坐标(j,i)，乘以1中求出的逆矩阵iR，转换到摄像机坐标系（_x,_y,_w）,并归一化得到z=1平面下的三维坐标(x,y,1)；
3.求出平面模型下像素点对应鱼眼半球模型下的极坐标(r, theta)。
4.利用鱼眼畸变模型求出拥有畸变时像素点对应的theta_d。
5.利用求出的theta_d值将三维坐标点重投影到二维图像平面得到(u,v)，(u,v)即为目标图像对应的畸变图像中像素点坐标
6.使用cv::Remap（）函数，根据mapx,mapy取出对应坐标位置的像素值赋值给目标图像，一般采用双线性插值法，得到畸变校正后的目标图像。
#include <opencv2\opencv.hpp> void cv::fisheye::initUndistortRectifyMap( InputArray K, InputArray D, InputArray R, InputArray P,
const cv::Size& size, int m1type, OutputArray map1, OutputArray map2 )
{
CV_Assert( m1type == CV_16SC2 || m1type == CV_32F || m1type <= );
map1.create( size, m1type <= ? CV_16SC2 : m1type );
map2.create( size, map1.type() == CV_16SC2 ? CV_16UC1 : CV_32F ); CV_Assert((K.depth() == CV_32F || K.depth() == CV_64F) && (D.depth() == CV_32F || D.depth() == CV_64F));
CV_Assert((P.empty() || P.depth() == CV_32F || P.depth() == CV_64F) && (R.empty() || R.depth() == CV_32F || R.depth() == CV_64F));
CV_Assert(K.size() == Size(, ) && (D.empty() || D.total() == ));
CV_Assert(R.empty() || R.size() == Size(, ) || R.total() * R.channels() == );
CV_Assert(P.empty() || P.size() == Size(, ) || P.size() == Size(, )); //从内参矩阵K中取出归一化焦距fx,fy; cx,cy
cv::Vec2d f, c;
if (K.depth() == CV_32F)
{
Matx33f camMat = K.getMat();
f = Vec2f(camMat(, ), camMat(, ));
c = Vec2f(camMat(, ), camMat(, ));
}
else
{
Matx33d camMat = K.getMat();
f = Vec2d(camMat(, ), camMat(, ));
c = Vec2d(camMat(, ), camMat(, ));
}
//从畸变系数矩阵D中取出畸变系数k1,k2,k3,k4
Vec4d k = Vec4d::all();
if (!D.empty())
k = D.depth() == CV_32F ? (Vec4d)*D.getMat().ptr<Vec4f>(): *D.getMat().ptr<Vec4d>(); //旋转矩阵RR转换数据类型为CV_64F，如果不需要旋转，则RR为单位阵
cv::Matx33d RR = cv::Matx33d::eye();
if (!R.empty() && R.total() * R.channels() == )
{
cv::Vec3d rvec;
R.getMat().convertTo(rvec, CV_64F);
RR = Affine3d(rvec).rotation();
}
else if (!R.empty() && R.size() == Size(, ))
R.getMat().convertTo(RR, CV_64F); //新的内参矩阵PP转换数据类型为CV_64F
cv::Matx33d PP = cv::Matx33d::eye();
if (!P.empty())
P.getMat().colRange(, ).convertTo(PP, CV_64F); //关键一步：新的内参矩阵*旋转矩阵，然后利用SVD分解求出逆矩阵iR，后面用到
cv::Matx33d iR = (PP * RR).inv(cv::DECOMP_SVD); //反向映射，遍历目标图像所有像素位置，找到畸变图像中对应位置坐标(u,v)，并分别保存坐标(u,v)到mapx和mapy中
for( int i = ; i < size.height; ++i)
{
float* m1f = map1.getMat().ptr<float>(i);
float* m2f = map2.getMat().ptr<float>(i);
short* m1 = (short*)m1f;
ushort* m2 = (ushort*)m2f; //二维图像平面坐标系->摄像机坐标系
double _x = i*iR(, ) + iR(, ),
_y = i*iR(, ) + iR(, ),
_w = i*iR(, ) + iR(, ); for( int j = ; j < size.width; ++j)
{
//归一化摄像机坐标系，相当于假定在Z=1平面上
double x = _x/_w, y = _y/_w; //求鱼眼半球体截面半径r
double r = sqrt(x*x + y*y);
//求鱼眼半球面上一点与光心的连线和光轴的夹角Theta
double theta = atan(r);
//畸变模型求出theta_d，相当于有畸变的角度值
double theta2 = theta*theta, theta4 = theta2*theta2, theta6 = theta4*theta2, theta8 = theta4*theta4;
double theta_d = theta * ( + k[]*theta2 + k[]*theta4 + k[]*theta6 + k[]*theta8);
//利用有畸变的Theta值，将摄像机坐标系下的归一化三维坐标，重投影到二维图像平面，得到(j,i)对应畸变图像中的(u,v)
double scale = (r == ) ? 1.0 : theta_d / r;
double u = f[]*x*scale + c[];
double v = f[]*y*scale + c[]; //保存(u,v)坐标到mapx,mapy
if( m1type == CV_16SC2 )
{
int iu = cv::saturate_cast<int>(u*cv::INTER_TAB_SIZE);
int iv = cv::saturate_cast<int>(v*cv::INTER_TAB_SIZE);
m1[j*+] = (short)(iu >> cv::INTER_BITS);
m1[j*+] = (short)(iv >> cv::INTER_BITS);
m2[j] = (ushort)((iv & (cv::INTER_TAB_SIZE-))*cv::INTER_TAB_SIZE + (iu & (cv::INTER_TAB_SIZE-)));
}
else if( m1type == CV_32FC1 )
{
m1f[j] = (float)u;
m2f[j] = (float)v;
} //这三条语句是上面 ”//二维图像平面坐标系->摄像机坐标系“的一部分，是矩阵iR的第一列，这样写能够简化计算
_x += iR(, );
_y += iR(, );
_w += iR(, );
}
}
}
2.普通相机模型的畸变校正
同样建议参考OpenCV官方文档阅读代码 https://docs.opencv.org/3.0-beta/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html#。
主要流程和上面Fisheye模型差不多，只有第4部分的畸变模型不一样，普通相机的畸变模型如下：
同样把源代码贴上，并加上注解：
#include <opencv2\opencv.hpp> void cv::initUndistortRectifyMap( InputArray _cameraMatrix, InputArray _distCoeffs,
InputArray _matR, InputArray _newCameraMatrix,
Size size, int m1type, OutputArray _map1, OutputArray _map2 )
{
Mat cameraMatrix = _cameraMatrix.getMat(), distCoeffs = _distCoeffs.getMat();
Mat matR = _matR.getMat(), newCameraMatrix = _newCameraMatrix.getMat(); if( m1type <= )
m1type = CV_16SC2;
CV_Assert( m1type == CV_16SC2 || m1type == CV_32FC1 || m1type == CV_32FC2 );
_map1.create( size, m1type );
Mat map1 = _map1.getMat(), map2;
if( m1type != CV_32FC2 )
{
_map2.create( size, m1type == CV_16SC2 ? CV_16UC1 : CV_32FC1 );
map2 = _map2.getMat();
}
else
_map2.release(); Mat_<double> R = Mat_<double>::eye(, );
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(, , CV_64F);
distCoeffs = .;
} CV_Assert( A.size() == Size(,) && A.size() == R.size() );
CV_Assert( Ar.size() == Size(,) || Ar.size() == Size(, )); //LU分解求新的内参矩阵Ar与旋转矩阵R乘积的逆矩阵iR
Mat_<double> iR = (Ar.colRange(,)*R).inv(DECOMP_LU);
const double* ir = &iR(,); //从旧的内参矩阵中取出光心位置u0,v0,和归一化焦距fx,fy
double u0 = A(, ), v0 = A(, );
double fx = A(, ), fy = A(, ); //尼玛14个畸变系数，不过大多用到的只有(k1,k2,p1,p2)，最多加一个k3，用不到的置为0
CV_Assert( distCoeffs.size() == Size(, ) || distCoeffs.size() == Size(, ) ||
distCoeffs.size() == Size(, ) || distCoeffs.size() == Size(, ) ||
distCoeffs.size() == Size(, ) || distCoeffs.size() == Size(, ) ||
distCoeffs.size() == Size(, ) || distCoeffs.size() == Size(, ) ||
distCoeffs.size() == Size(, ) || distCoeffs.size() == Size(, )); if( distCoeffs.rows != && !distCoeffs.isContinuous() )
distCoeffs = distCoeffs.t(); const double* const distPtr = distCoeffs.ptr<double>();
double k1 = distPtr[];
double k2 = distPtr[];
double p1 = distPtr[];
double p2 = distPtr[];
double k3 = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .;
double k4 = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .;
double k5 = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .;
double k6 = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .;
double s1 = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .;
double s2 = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .;
double s3 = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .;
double s4 = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .;
double tauX = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .;
double tauY = distCoeffs.cols + distCoeffs.rows - >= ? distPtr[] : .; //tauX,tauY这个是什么梯形畸变，用不到的话matTilt为单位阵
// Matrix for trapezoidal distortion of tilted image sensor
cv::Matx33d matTilt = cv::Matx33d::eye();
cv::detail::computeTiltProjectionMatrix(tauX, tauY, &matTilt); for( int i = ; i < size.height; i++ )
{
float* m1f = map1.ptr<float>(i);
float* m2f = map2.empty() ? : map2.ptr<float>(i);
short* m1 = (short*)m1f;
ushort* m2 = (ushort*)m2f; //利用逆矩阵iR将二维图像坐标(j,i)转换到摄像机坐标系(_x,_y,_w)
double _x = i*ir[] + ir[], _y = i*ir[] + ir[], _w = i*ir[] + ir[]; for( int j = ; j < size.width; j++, _x += ir[], _y += ir[], _w += ir[] )
{
//摄像机坐标系归一化，令Z=1平面
double w = ./_w, x = _x*w, y = _y*w;
//这一部分请看OpenCV官方文档，畸变模型部分
double x2 = x*x, y2 = y*y;
double r2 = x2 + y2, _2xy = *x*y;
double kr = ( + ((k3*r2 + k2)*r2 + k1)*r2)/( + ((k6*r2 + k5)*r2 + k4)*r2);
double xd = (x*kr + p1*_2xy + p2*(r2 + *x2) + s1*r2+s2*r2*r2);
double yd = (y*kr + p1*(r2 + *y2) + p2*_2xy + s3*r2+s4*r2*r2);
//根据求取的xd,yd将三维坐标重投影到二维畸变图像坐标(u,v)
cv::Vec3d vecTilt = matTilt*cv::Vec3d(xd, yd, );
double invProj = vecTilt() ? ./vecTilt() : ;
double u = fx*invProj*vecTilt() + u0;
double v = fy*invProj*vecTilt() + v0;
//保存u,v的值到Mapx,Mapy中
if( m1type == CV_16SC2 )
{
int iu = saturate_cast<int>(u*INTER_TAB_SIZE);
int iv = saturate_cast<int>(v*INTER_TAB_SIZE);
m1[j*] = (short)(iu >> INTER_BITS);
m1[j*+] = (short)(iv >> INTER_BITS);
m2[j] = (ushort)((iv & (INTER_TAB_SIZE-))*INTER_TAB_SIZE + (iu & (INTER_TAB_SIZE-)));
}
else if( m1type == CV_32FC1 )
{
m1f[j] = (float)u;
m2f[j] = (float)v;
}
else
{
m1f[j*] = (float)u;
m1f[j*+] = (float)v;
}
}
}
}
如有错误，望不吝赐教！
另附上CUDA实现两种畸变校正方法的代码，放在我的码云上：https://gitee.com/rxdj/camera-calibration.git。见cudaUndistort中的两个.cu文件