标签：
图像矫正的本质，其实就是重投影的过程，即【像素坐标→物理坐标→像素坐标】的过程。只不过在重投影过程中我们可以改变投影矩阵(修改后的投影矩阵我把它称为扩展投影矩阵)从而模拟镜头缩放和平移的效果。
图像矫正可通过两种方式执行，我称之为正向矫正和逆向矫正。正向矫正是通过畸变坐标算出标准坐标，而逆向矫正是通过标准坐标算出畸变坐标。Opencv中UndistortPoints就是执行的正向矫正过程，而initUndistortRectifyMap执行的是逆向矫正过程。
正向矫正的流程为：畸变像素坐标→畸变物理坐标→标准物理坐标→标准像素坐标
逆向矫正的流程为：标准像素坐标→标准物理坐标→畸变物理坐标→畸变像素坐标
从流程可以看出，不管正向矫正还是逆向矫正都要先过渡到物理坐标，然后再转换为像素坐标，这也是为什么说图像矫正的本质其实就是重投影的原因。
畸变物理坐标和标准物理坐标之间的桥梁是畸变向量。而物理坐标和像素坐标之间的桥梁是投影矩阵。其中，畸变像素坐标和畸变物理坐标之间的桥梁是相机原生的投影矩阵，而标准像素坐标和标准物理坐标之间的桥梁既可以是相机的原生的投影矩阵也可以是自定义的扩展投影矩阵，扩展投影矩阵的目的是模拟镜头缩放和平移的效果，从而满足个性化需求。
OpenCV中单图像矫正的主要函数有：
undistortPoints
icvGetRectangles
getOptimalNewCameraMatrix
initUndistortRectifyMap
为了代码和公式更好地对应，这四函数我都移除对双目图像矫正的支持，并对代码进行了整理和简化。
对于undistortPoints和initUndistortRectifyMap，查看源查源码及源码中的注释即可与理论对应上。这里对讲解一下icvGetRectangles和getOptimalNewCameraMatrix的功能和原理
icvGetRectangles功能是获取标准图像的有效像素(有效像素指在畸变图像中有对应像素的像素)所构成的区域的最大内接矩阵和最小外接矩阵，其实现方式是：取畸变图像中的一些特殊点(四条边上的点和中间区域的一些点，详见源码)，调用UndistortPoints算出这些点的标准坐标，通过这些标准坐标即可算出标准图像的有效像素所构成的区域的最大内接矩阵和最小外接矩阵。
icvGetOptimalNewCameraMatrix功能是设计扩展投影矩阵，其中个性化参数包括：是否使主点是矫正后的图像的中点，以及是否使矫正后的图像包含畸变图像映射过来的所有点。这里要注意的是，若矫正后的图像包含畸变图像映射过来的所有点，则也包含一些在畸变图像中图像中无对应的点；若矫正后的图像只包含畸变图像映射过来的点，则也会遗漏一些畸变图像映射过来的点
1 void cv_UndistortPoints_monocular(Point2f *srcPoints, Point2f *dstPoints, int n, double *A, double *D, double *AA) 2 { 3 //1.获取投影参数 4 double ifx = 1. / A[0]; 5 double ify = 1. / A[4]; 6 double icx = -ifx * A[2]; 7 double icy = -ify * A[5]; 8 9 //2.获取扩展投影参数 10 double ffx = 1; 11 double ffy = 1; 12 double ccx = 0; 13 double ccy = 0; 14 if (AA != 0) 15 {//若AA为空则返回标准物理坐标 16 ffx = AA[0]; 17 ffy = AA[4]; 18 ccx = AA[2]; 19 ccy = AA[5]; 20 } 21 22 //3.获取畸变参数 23 double k1 = D[0]; 24 double k2 = D[1]; 25 double p1 = D[2]; 26 double p2 = D[3]; 27 double k3 = D[4]; 28 double k4 = D[5]; 29 double k5 = D[6]; 30 double k6 = D[7]; 31 double s1 = D[8]; 32 double s2 = D[9]; 33 double s3 = D[10]; 34 double s4 = D[11]; 35 //double a = D[12]; 36 //double b = D[13]; 37 38 //4.根据畸变像素坐标计算标准像素坐标 39 for (int i = 0; i < n; i++) 40 { 41 //4.1根据畸变像素坐标(u,v)和投影矩阵A反解畸变物理坐标(xxx,yyy) 42 double xxx = srcPoints[i].x * ifx + icx; 43 double yyy = srcPoints[i].y * ify + icy; 44 //4.2根据畸变物理坐标(xxx,yyy)和畸变向量D反解标准物理坐标(xx, yy) 45 double xx, xxx0 = xxx; 46 double yy, yyy0 = yyy; 47 for (int j = 0; j < 5; j++) 48 { 49 double tm1 = xxx * xxx; 50 double tm2 = yyy * yyy; 51 double tm3 = tm1 + tm2; 52 double tm4 = 2 * xxx * yyy; 53 double tm5 = (1 + ((k6 * tm3 + k5)*tm3 + k4)*tm3) / (1 + ((k3 * tm3 + k2)*tm3 + k1)*tm3); 54 double deltaX = p1 * tm4 + p2 * (tm3 + 2 * tm1) + s1 * tm3 + s2 * tm3 * tm3; 55 double deltaY = p1 * (tm3 + 2 * tm2) + p2 * tm4 + s3 * tm3 + s4 * tm3 * tm3; 56 xxx = xx = (xxx0 - deltaX)*tm5; 57 yyy = yy = (yyy0 - deltaY)*tm5; 58 } 59 //4.3根据标准物理坐标(xx,yy)和扩展投影矩阵矩阵AA计算标准像素坐标(uu,vv) 60 dstPoints[i].x = (float)(ffx * xx + ccx); 61 dstPoints[i].y = (float)(ffy * yy + ccy); 62 } 63 }
1 void cv_GetRectangles_monocular(double *A, double *D, double *AA, Size srcSize, Rect2f *inner, Rect2f *outer) 2 { 3 //1.取得一定数量的畸变像素点 4 int N = 9; 5 Point2f points[81]; 6 for (int y = 0, k = 0; y < N; y++) 7 for (int x = 0; x < N; x++) 8 points[k++] = Point2f((float)x*srcSize.width / (N - 1), (float)y*srcSize.height / (N - 1)); 9 10 11 //2.矫正这些畸变像素点 12 cv_UndistortPoints_monocular(points, points, 81, A, D, AA); 13 14 //3.根据标准像素点计算最大内接矩形和最小外接矩形 15 float inX1 = -FLT_MAX, inX2 = FLT_MAX, inY1 = -FLT_MAX, inY2 = FLT_MAX; 16 float outX1 = FLT_MAX, outX2 = -FLT_MAX, outY1 = FLT_MAX, outY2 = -FLT_MAX; 17 for (int y = 0, k = 0; y < N; y++) 18 for (int x = 0; x < N; x++) 19 { 20 Point2f p = points[k++]; 21 //对于最小外接矩形 22 outX1 = __min(outX1, p.x); 23 outY1 = __min(outY1, p.y); 24 outX2 = __max(outX2, p.x); 25 outY2 = __max(outY2, p.y); 26 //对于最大内接矩形 27 if (x == 0) inX1 = __max(inX1, p.x); 28 if (y == 0) inY1 = __max(inY1, p.y); 29 if (x == N - 1) inX2 = __min(inX2, p.x); 30 if (y == N - 1) inY2 = __min(inY2, p.y); 31 } 32 *inner = Rect2f(inX1, inY1, inX2 - inX1, inY2 - inY1); 33 *outer = Rect2f(outX1, outY1, outX2 - outX1, outY2 - outY1); 34 }
1 void cv_GetOptimalNewCameraMatrix(Mat_<double> &A, Mat_<double> &D, Size srcSize, Size newSize, double alpha, bool centerPrincipalPoint, Mat_<double> &AA, Rect &validROI) 2 { 3 if (centerPrincipalPoint == true) 4 { 5 Rect2f inner; 6 Rect2f outer; 7 cv_GetRectangles_monocular((double*)A.data, (double*)D.data, (double*)A.data, srcSize, &inner, &outer); 8 9 double cx = A(0, 2); 10 double cy = A(1, 2); 11 double ccx = (newSize.width - 1)*0.5; 12 double ccy = (newSize.height - 1)*0.5; 13 14 double s0 = max(max(max((double)ccx / (cx - inner.x), (double)ccy / (cy - inner.y)), (double)ccx / (inner.x + inner.width - cx)), (double)ccy / (inner.y + inner.height - cy)); 15 double s1 = min(min(min((double)ccx / (cx - outer.x), (double)ccy / (cy - outer.y)), (double)ccx / (outer.x + outer.width - cx)), (double)ccy / (outer.y + outer.height - cy)); 16 double s = s0*(1 - alpha) + s1*alpha; //alpha==1则矫正后的图像包含畸变图像映射过来的所有点 17 18 AA(0, 0) = A(0, 0) * s; 19 AA(0, 1) = A(0, 1); 20 AA(0, 2) = ccx; 21 AA(1, 0) = A(1, 0); 22 AA(1, 1) = A(1, 1) * s; 23 AA(1, 2) = ccy; 24 AA(2, 0) = A(2, 0); 25 AA(2, 1) = A(2, 1); 26 AA(2, 2) = A(2, 2); 27 28 inner = Rect2f((float)((inner.x - cx)*s + ccx), (float)((inner.y - cy)*s + ccy), (float)(inner.width*s), (float)(inner.height*s)); 29 validROI = Rect(cvCeil(inner.x), cvCeil(inner.y), cvFloor(inner.width), cvFloor(inner.height)) & Rect(0, 0, newSize.width, newSize.height); 30 } 31 else 32 { 33 //1.获取基于物理坐标的最大内接矩形和最小外接矩阵 34 Rect2f inner; 35 Rect2f outer; 36 cv_GetRectangles_monocular((double*)A.data, (double*)D.data, 0, srcSize, &inner, &outer); 37 38 //2.基于最大内接矩阵计算扩展投影矩阵 39 double ffx0 = (newSize.width - 1) / inner.width; 40 double ffy0 = (newSize.height - 1) / inner.height; 41 double ccx0 = -ffx0 * inner.x; 42 double ccy0 = -ffy0 * inner.y; 43 44 //3.基于最小外接矩阵计算扩展投影矩阵 45 double ffx1 = (newSize.width - 1) / outer.width; 46 double ffy1 = (newSize.height - 1) / outer.height; 47 double ccx1 = -ffx1 * outer.x; 48 double ccy1 = -ffy1 * outer.y; 49 50 //4.获得最佳的扩展投影矩阵 51 AA(0, 0) = ffx0*(1 - alpha) + ffx1*alpha; 52 AA(0, 1) = A(0, 1); 53 AA(0, 2) = ccx0*(1 - alpha) + ccx1*alpha; 54 AA(1, 0) = A(1, 0); 55 AA(1, 1) = ffy0*(1 - alpha) + ffy1*alpha; 56 AA(1, 2) = ccy0*(1 - alpha) + ccy1*alpha; 57 AA(2, 0) = A(2, 0); 58 AA(2, 1) = A(2, 1); 59 AA(2, 2) = A(2, 2); 60 61 cv_GetRectangles_monocular((double*)A.data, (double*)D.data, (double*)AA.data, srcSize, &inner, &outer); 62 validROI = (Rect)inner & Rect(0, 0, newSize.width, newSize.height); 63 } 64 }
1 void cv_InitUndistortRectifyMap_monocular(Mat_<double> &A, Mat_<double> &D, Mat_<double> &AA, Mat_<float> &map1, Mat_<float> &map2) 2 { 3 //1.获取投影参数 4 double cx = A(0, 2); 5 double cy = A(1, 2); 6 double fx = A(0, 0); 7 double fy = A(1, 1); 8 9 //2.获取扩展投影参数 10 double iffx = 1.0 / AA(0, 0); 11 double iffy = 1.0 / AA(1, 1); 12 double iccx = -iffx * AA(0, 2); 13 double iccy = -iffy * AA(1, 2); 14 15 //3.获取畸变参数 16 double k1 = D(0); 17 double k2 = D(1); 18 double p1 = D(2); 19 double p2 = D(3); 20 double k3 = D.total() >= 5 ? D(4) : 0.; 21 double k4 = D.total() >= 8 ? D(5) : 0.; 22 double k5 = D.total() >= 8 ? D(6) : 0.; 23 double k6 = D.total() >= 8 ? D(7) : 0.; 24 double s1 = D.total() >= 12 ? D(8) : 0.; 25 double s2 = D.total() >= 12 ? D(9) : 0.; 26 double s3 = D.total() >= 12 ? D(10) : 0.; 27 double s4 = D.total() >= 12 ? D(11) : 0.; 28 //double a = D.total() >= 14 ? D(12) : 0.; 29 //double b = D.total() >= 14 ? D(13) : 0.; 30 31 //4.根据标准像素坐标计算畸变像素坐标 32 for (int vv = 0; vv < map1.size().height; vv++) 33 for (int uu = 0; uu < map2.size().width; uu++) 34 { 35 //4.1根据标准像素坐标(uu,vv)和扩展投影矩阵AA计算标准物理坐标(xx,yy) 36 double xx = uu * iffx + iccx; 37 double yy = vv * iffy + iccy; 38 //4.2根据标准物理坐标(xx,yy)和畸变向量D计算畸变物理坐标(xxx, yyy) 39 double tm1 = xx * xx; 40 double tm2 = yy * yy; 41 double tm3 = tm1 + tm2; 42 double tm4 = 2 * xx * yy; 43 double tm5 = (1 + ((k3*tm3 + k2)*tm3 + k1)*tm3) / (1 + ((k6*tm3 + k5)*tm3 + k4)*tm3); 44 double xxx = (xx*tm5 + p1*tm4 + p2*(tm3 + 2 * tm1) + s1*tm3 + s2*tm3*tm3); 45 double yyy = (yy*tm5 + p1*(tm3 + 2 * tm2) + p2*tm4 + s3*tm3 + s4*tm3*tm3); 46 //4.3根据畸变物理坐标(xxx,yyy)和投影矩阵矩阵A计算畸变像素坐标(u,v) 47 double u = fx*xxx + cx; 48 double v = fy*yyy + cy; 49 map1(vv, uu) = (float)u; 50 map2(vv, uu) = (float)v; 51 } 52 }
标签：
原文地址：http://www.cnblogs.com/dzyBK/p/5579206.html
