OpenCV之getOptimalNewCameraMatrix

去畸变后的图像四周会出现黑色区域，如果畸变较大，如鱼眼镜头，四周黑色区域将会更大。
opencv中给我们提供了一个函数getOptimalNewCameraMatrix()，用于去除畸变矫正后图像四周黑色的区域。

该函数根据给定参数alpha计算最优的新相机内参矩阵。alpha=0，则去除所有黑色区域，alpha=1，则保留所有原始图像像素，其他值则得到介于两者之间的效果。

通过该函数，我们得到新的相机内参矩阵，然后利用initUndistortRectifyMap()即可获得用于remap()的映射矩阵。

主要步骤

获取内切和外切矩形

通过将整个图像划分9*9网格，然后对每个网格点去畸变。
通过比较这些去畸变后的网格点坐标，确定内切和外切矩形。

源码中提到：Get inscribed and circumscribed rectangles in normalized (independent of camera matrix) coordinates.
这是因为，在调用icvGetRectangles()获取内切和外切矩形时，R和newCameraMatrix为空。从cvUndistortPoints()源码可以知道，当R和newCameraMatrix为空，尤其是newCameraMatrix为空时，像素点去畸变后得到的是相机坐标系下的坐标。

static void
icvGetRectangles(const CvMat* cameraMatrix, const CvMat* distCoeffs,const CvMat* R, const CvMat* newCameraMatrix, CvSize imgSize,cv::Rect_<float>& inner, cv::Rect_<float>& outer)
{const int N = 9;int x, y, k;cv::Ptr<CvMat> _pts(cvCreateMat(1, N * N, CV_32FC2));CvPoint2D32f* pts = (CvPoint2D32f*)(_pts->data.ptr);for (y = k = 0; y < N; y++)for (x = 0; x < N; x++)pts[k++] = cvPoint2D32f((float)x * imgSize.width / (N - 1),(float)y * imgSize.height / (N - 1));// 由于R和newCameraMatrix为0，因此，去畸变后的坐标为相机坐标系下的坐标，（如 mm）cvUndistortPoints(_pts, _pts, cameraMatrix, distCoeffs, R, newCameraMatrix);float iX0 = -FLT_MAX, iX1 = FLT_MAX, iY0 = -FLT_MAX, iY1 = FLT_MAX;float oX0 = FLT_MAX, oX1 = -FLT_MAX, oY0 = FLT_MAX, oY1 = -FLT_MAX;// find the inscribed rectangle.// the code will likely not work with extreme rotation matrices (R) (>45%)for (y = k = 0; y < N; y++)for (x = 0; x < N; x++){CvPoint2D32f p = pts[k++];oX0 = MIN(oX0, p.x);oX1 = MAX(oX1, p.x);oY0 = MIN(oY0, p.y);oY1 = MAX(oY1, p.y);if (x == 0)iX0 = MAX(iX0, p.x);if (x == N - 1)iX1 = MIN(iX1, p.x);if (y == 0)iY0 = MAX(iY0, p.y);if (y == N - 1)iY1 = MIN(iY1, p.y);}inner = cv::Rect_<float>(iX0, iY0, iX1 - iX0, iY1 - iY0);outer = cv::Rect_<float>(oX0, oY0, oX1 - oX0, oY1 - oY0);
}

新的相机内参矩阵

通过icvGetRectangles()我们得到的内切和外切矩形（相机坐标系下）。接下来我们需要一个新的内参矩阵，将指定区域（内切矩形、外切矩形或介于两者之间的区域）重新投影到整个图像区域。

分别计算出内切矩形和外切矩形对应的相机内参矩阵，然后根据alpha，加权得到最终的相机内参矩阵。

double fx0 = (newImgSize.width - 1) / inner.width; // 可以认为inner所在坐标系下单位长度对应的像素个数
double fy0 = (newImgSize.height - 1) / inner.height;
double cx0 = -fx0 * inner.x; // 将矩形的左上角，移动到图像坐标系的原点
double cy0 = -fy0 * inner.y;// Projection mapping outer rectangle to viewport
double fx1 = (newImgSize.width - 1) / outer.width;
double fy1 = (newImgSize.height - 1) / outer.height;
double cx1 = -fx1 * outer.x;
double cy1 = -fy1 * outer.y;// Interpolate between the two optimal projections
M[0][0] = fx0 * (1 - alpha) + fx1 * alpha;
M[1][1] = fy0 * (1 - alpha) + fy1 * alpha;
M[0][2] = cx0 * (1 - alpha) + cx1 * alpha;
M[1][2] = cy0 * (1 - alpha) + cy1 * alpha;

有效的像素区域

上述我们获得了新的相机内参矩阵，利用这个矩阵，我们可以得到像素坐标系下的内切和外切矩形。

if (validPixROI)
{// 利用上述得到的投影矩阵（“新相机内参矩阵”），得到像素坐标系下的内切和外切矩形icvGetRectangles(cameraMatrix, distCoeffs, 0, &matM, imgSize, inner, outer);// 这里我们只关注内切矩形cv::Rect r = inner;r &= cv::Rect(0, 0, newImgSize.width, newImgSize.height);*validPixROI = cvRect(r);
}

完整代码

void cvGetOptimalNewCameraMatrix(const CvMat* cameraMatrix, const CvMat* distCoeffs,CvSize imgSize, double alpha,CvMat* newCameraMatrix, CvSize newImgSize,CvRect* validPixROI, int centerPrincipalPoint)
{cv::Rect_<float> inner, outer;newImgSize = newImgSize.width * newImgSize.height != 0 ? newImgSize : imgSize;double M[3][3];CvMat matM = cvMat(3, 3, CV_64F, M);cvConvert(cameraMatrix, &matM);if (centerPrincipalPoint){double cx0 = M[0][2];double cy0 = M[1][2];double cx = (newImgSize.width - 1) * 0.5;double cy = (newImgSize.height - 1) * 0.5;icvGetRectangles(cameraMatrix, distCoeffs, 0, cameraMatrix, imgSize, inner, outer);double s0 = std::max(std::max(std::max((double)cx / (cx0 - inner.x), (double)cy / (cy0 - inner.y)),  (double)cx / (inner.x + inner.width - cx0)), (double)cy / (inner.y + inner.height - cy0));double s1 = std::min(std::min(std::min((double)cx / (cx0 - outer.x), (double)cy / (cy0 - outer.y)), (double)cx / (outer.x + outer.width - cx0)), (double)cy / (outer.y + outer.height - cy0));double s = s0 * (1 - alpha) + s1 * alpha;M[0][0] *= s;M[1][1] *= s;M[0][2] = cx;M[1][2] = cy;if (validPixROI){inner = cv::Rect_<float>((float)((inner.x - cx0) * s + cx),(float)((inner.y - cy0) * s + cy),(float)(inner.width * s),(float)(inner.height * s));cv::Rect r(cvCeil(inner.x), cvCeil(inner.y), cvFloor(inner.width), cvFloor(inner.height));r &= cv::Rect(0, 0, newImgSize.width, newImgSize.height);*validPixROI = cvRect(r);}}else{// Get inscribed and circumscribed rectangles in normalized// (independent of camera matrix) coordinatesicvGetRectangles(cameraMatrix, distCoeffs, 0, 0, imgSize, inner, outer);// Projection mapping inner rectangle to viewportdouble fx0 = (newImgSize.width - 1) / inner.width; // 可以认为inner所在坐标系下单位长度对应的像素个数double fy0 = (newImgSize.height - 1) / inner.height;double cx0 = -fx0 * inner.x; // 将矩形的左上角，移动到图像坐标系的原点double cy0 = -fy0 * inner.y;// Projection mapping outer rectangle to viewportdouble fx1 = (newImgSize.width - 1) / outer.width;double fy1 = (newImgSize.height - 1) / outer.height;double cx1 = -fx1 * outer.x;double cy1 = -fy1 * outer.y;// Interpolate between the two optimal projectionsM[0][0] = fx0 * (1 - alpha) + fx1 * alpha;M[1][1] = fy0 * (1 - alpha) + fy1 * alpha;M[0][2] = cx0 * (1 - alpha) + cx1 * alpha;M[1][2] = cy0 * (1 - alpha) + cy1 * alpha;if (validPixROI){// 利用上述得到的投影矩阵（新相机内参矩阵），得到像素坐标系下的内切和外切矩形icvGetRectangles(cameraMatrix, distCoeffs, 0, &matM, imgSize, inner, outer);// 这里我们只关注内切矩形cv::Rect r = inner;r &= cv::Rect(0, 0, newImgSize.width, newImgSize.height);*validPixROI = cvRect(r);}}cvConvert(&matM, newCameraMatrix);
}

最后

该函数提供了一种思路，通过划分网格来确定图像有效区域的矩形边界。
在去畸变函数cvUndistortPoints()中，通过内参矩阵将像素点转化到相机坐标系下，然后根据畸变系数进行畸变矫正。如果给定旋转矩阵matR，可以对去畸变后的点进行旋转变换（可用于旋转补偿，如视频防抖）。如果给定matP，则可以对去畸变后的点重新投影到像素坐标系（如投影到指定图像区域）。

void cvUndistortPoints(const CvMat* _src, CvMat* _dst, const CvMat* _cameraMatrix,const CvMat* _distCoeffs, const CvMat* matR, const CvMat* matP)