图像桶形畸变矫正

转载自：_寒潭雁影的 C/C++ 图像处理（9）------图像の桶形畸变矫正

原文：https://blog.csdn.net/weixinhum/article/details/50719031

        广角镜头的摄像设备拍摄出来的图像经常会有桶形畸变的问题。原因在于广角镜头使用的是凸透镜，初中物理知识告诉我们凸透镜会对光线起汇聚作用，这是光的折射造成的。而离镜头中心越远，折射效果越强，因而其拍出来的照片会以镜头中心为圆心，呈圆形向外扩展失真，如下图所示：

        像上面这样的图像，如果用在一些还原性要求较高的场景是不行的，需要对图像做畸变矫正。而由于很多时候我们并不知道摄像头的物理参数和其他一些信息，只是拿到一个可以输出画面的摄像头，因而比较常采取的桶形畸变矫正算法为多项式修正算法。下面本文将展开介绍该算法并给出矫正后的结果，如果感觉本文的描述不够清楚，也可以参考这个网址中的文章（本文也是参考该文章的）。

    下面进入主题。

    首先，我们需要制作或者购买一个同心圆标定板，如下图所示：

         然后将该同心圆标定板的圆心和镜头中心（一般镜头中心就是图像的中心点，当然也有像本文章实验时这样比较奇葩的，镜头中心点在图像中心点下方的情况）对准获得一张待处理图像（在采集图像的时候要将同心圆标定板靠近镜头，使得图像中的同心圆覆盖到图像边界，不然图像的边界位置矫正效果会不大好），如下图所示：

       由上面图像我们可以看到，越外围的同心圆凑的越紧，也就是说偏离镜头中心越远，图像畸变越厉害，与上面的论述一致。

        这里我们认为离镜头中心最近的同心圆没有畸变，然后测出其与镜头中心的距离为143，则在没有畸变的图像中从近到远的同心圆（1-7）半径应该为R1（143，286，429，572，715，858，1001）。我们在上面图像中测出来的同心圆半径为R（143,260.5,341.5,395,430.5,454,469）。

         对图像做畸变矫正，实际上就是将图像中的像素点放到其理论上应该在的位置而已，因此我们可以由上面的两组值推算出计算出畸变后图像像素与镜头中心的距离和理论距离的比值K，再根据这一组K值来拟合出一道R与K的关系式，从而得到图像中像素与镜头中心距离与实际距离的映射关系，以此矫正图像。

    我们将理论值与实际畸变后的值相除，得到畸变矫正系数K（1，1.0979，1.2562，1.4481，1.6609，1.8899，2.1343），将R-K绘制成折线图如下所示：

        由于上面我们以第一个圆半径作为标准，半径之前的数据为空，不利于后面计算，因此我们将多取一个（0,1）点，折线图修正如下：

        粗略地，有了上图，我们已经可以求出图像像素点距离镜头中心各个距离所对应的K值。但由于该图只是折线图，拟合的并不好，且要分段编程比较麻烦，所以我们用四次多项式重新拟合这些点（MATLAB实现），得到四项多项式公式：

        拟合结果如下图所示：

        有了这个公式，便有了畸变图像像素点与正常图像像素点的映射关系了。借由这个映射关系编写下面的代码我们便可得到初步的矫正图像。

char *filename = "2542.bmp";//图像路径  
img = cvLoadImage(filename);
paintObj.drawpicinit(img, IDC_STATIC_SHOW, this);
CPoint lenscenter(480,295);//镜头中心在图像中的位置
uchar *data = (uchar*)img->imageData;//原图像数据
uchar *newimgdata = new uchar[img->widthStep*img->height];//新图像数据
for (int row = 0; row<img->height; row++)//操作数据区,要注意OpenCV的RGB的存储顺序为GBR  
	for (int cols = 0; cols<img->width; cols++)//示例为亮度调节  
	{
		int pointsrc = row*img->widthStep + cols * 3;
		double r = sqrt((row- lenscenter.y)*(row - lenscenter.y) + (cols - lenscenter.x)*(cols - lenscenter.x));	
		double s = 1.0008 - 0.0031*r + 3.7101*pow(10, -5)*pow(r, 2) - 1.3491*pow(10, -7)*pow(r, 3)+1.7184*pow(10, -10)*pow(r, 4);//比例
		int x = (cols - lenscenter.x)*s + lenscenter.x;
		int y = (row - lenscenter.y)*s + lenscenter.y;
		
		if (x >= 0 && x < img->width && y >= 0 && y < img->height)
		{
			int pointdrc = y*img->widthStep + x * 3;
			newimgdata[pointdrc + 0] = data[pointsrc + 0];
			newimgdata[pointdrc + 1] = data[pointsrc + 1];
			newimgdata[pointdrc + 2] = data[pointsrc + 2];
		}
	}
memset(data, 0, img->widthStep*img->height);
memcpy(data, newimgdata, img->widthStep*img->height);
cvSaveImage("save33.bmp", img);
cvReleaseImage(&img);

       由上图可以看到图像的桶形畸变基本被矫正了，但是多出了很多波纹，这个是没有用差值算法造成的。通过改变R-K成R1-K，然后再进行拟合，得到如下的四项多项式公式：

       拟合结果如下图所示：

        有了这个公式便可以算出畸变矫正后的图像像素点与畸变图像的像素点的映射关系，与之前不同的是这个公式以畸变矫正后的图像像素点距离镜头中心的距离作为输入值，方便我们后面对图像进行双线性插值算法的处理。

        那么我们依据这个公式编写代码如下：

//畸变矫正
char *filename = "2542.bmp";//图像路径  
img = cvLoadImage(filename);
paintObj.drawpicinit(img, IDC_STATIC_SHOW, this);
CPoint lenscenter(480,295);//镜头中心在图像中的位置
uchar *data = (uchar*)img->imageData;//原图像数据
uchar *newimgdata = new uchar[img->widthStep*img->height];//新图像数据
for (int row = 0; row<img->height; row++)//操作数据区,要注意OpenCV的RGB的存储顺序为GBR  
	for (int cols = 0; cols<img->width; cols++)//示例为亮度调节  
	{
		int pointsrc = row*img->widthStep + cols * 3;
		double r = sqrt((row- lenscenter.y)*(row - lenscenter.y) + (cols - lenscenter.x)*(cols - lenscenter.x));	
		double s = 0.9998 - 4.2932*pow(10, -4)*r + 3.4327*pow(10, -6)*pow(r, 2) -2.8526*pow(10, -9)*pow(r, 3)+9.8223*pow(10, -13)*pow(r, 4);//比例
		double d_original_img_wnum = (cols - lenscenter.x)/s + lenscenter.x;
		double d_original_img_hnum = (row - lenscenter.y)/s + lenscenter.y;
		
		int i_original_img_hnum = d_original_img_hnum;
		int i_original_img_wnum = d_original_img_wnum;
		double distance_to_a_x = d_original_img_wnum - i_original_img_wnum;//在原图像中与a点的水平距离  
		double distance_to_a_y = d_original_img_hnum - i_original_img_hnum;//在原图像中与a点的垂直距离  
 
		int original_point_a = i_original_img_hnum*img->widthStep + i_original_img_wnum * 3;//数组位置偏移量，对应于图像的各像素点RGB的起点,相当于点A    
		int original_point_b = i_original_img_hnum*img->widthStep + (i_original_img_wnum + 1) * 3;//数组位置偏移量，对应于图像的各像素点RGB的起点,相当于点B  
		int original_point_c = (i_original_img_hnum + 1)*img->widthStep + i_original_img_wnum * 3;//数组位置偏移量，对应于图像的各像素点RGB的起点,相当于点C   
		int original_point_d = (i_original_img_hnum + 1)*img->widthStep + (i_original_img_wnum + 1) * 3;//数组位置偏移量，对应于图像的各像素点RGB的起点,相当于点D  
		if (row == img->height - 1)
		{
			original_point_c = original_point_a;
			original_point_d = original_point_b;
		}
		if (cols == img->width - 1)
		{
			original_point_a = original_point_b;
			original_point_c = original_point_d;
		}
		newimgdata[pointsrc] =
			data[original_point_a] * (1 - distance_to_a_x)*(1 - distance_to_a_y) +
			data[original_point_b] * distance_to_a_x*(1 - distance_to_a_y) +
			data[original_point_c] * distance_to_a_y*(1 - distance_to_a_x) +
			data[original_point_c] * distance_to_a_y*distance_to_a_x;
		newimgdata[pointsrc + 1] =
			data[original_point_a + 1] * (1 - distance_to_a_x)*(1 - distance_to_a_y) +
			data[original_point_b + 1] * distance_to_a_x*(1 - distance_to_a_y) +
			data[original_point_c + 1] * distance_to_a_y*(1 - distance_to_a_x) +
			data[original_point_c + 1] * distance_to_a_y*distance_to_a_x;
		newimgdata[pointsrc + 2] =
			data[original_point_a + 2] * (1 - distance_to_a_x)*(1 - distance_to_a_y) +
			data[original_point_b + 2] * distance_to_a_x*(1 - distance_to_a_y) +
			data[original_point_c + 2] * distance_to_a_y*(1 - distance_to_a_x) +
			data[original_point_c + 2] * distance_to_a_y*distance_to_a_x;
 
	}
memset(data, 0, img->widthStep*img->height);
memcpy(data, newimgdata, img->widthStep*img->height);
cvSaveImage("save33.bmp", img);
cvReleaseImage(&img);  

         最近发现用OpenCV2.0的变量操作快感提升了一个层次，因此贴出用2.0变量所写的代码以供参考

#include <opencv2/opencv.hpp>    
#include <time.h> 
using namespace cv;
 
void main()
{
	long time = clock();  
	char *filename = "4.bmp";//图像路径    
	Mat img = imread(filename);
	Mat drcimg(img.rows, img.cols, CV_8UC3);
	cv::imshow("矫正前", img);	
	CvPoint lenscenter(img.cols / 2, img.rows / 2);//镜头中心在图像中的位置  
	CvPoint src_a, src_b, src_c, src_d;//a、b、c、d四个顶点
	//矫正参数
	double r;//矫正前像素点跟镜头中心的距离
	double s;//矫正后像素点跟镜头中心的距离
	CvPoint2D32f mCorrectPoint;//矫正后点坐标
	double distance_to_a_x, distance_to_a_y;//求得中心点和边界的距离
 
	for (int row = 0; row< img.rows; row++)//操作数据区,要注意OpenCV的RGB的存储顺序为GBR    
		for (int cols = 0; cols<img.cols; cols++)//示例为亮度调节    
		{	
			r = sqrt((row - lenscenter.y)*(row - lenscenter.y) + (cols - lenscenter.x)*(cols - lenscenter.x))*0.8;
			s = 0.9998 - 4.2932*pow(10, -4)*r + 3.4327*pow(10, -6)*pow(r, 2) - 2.8526*pow(10, -9)*pow(r, 3) + 9.8223*pow(10, -13)*pow(r, 4);//比例  
			mCorrectPoint = cvPoint2D32f((cols - lenscenter.x) / s *1.35 + lenscenter.x, (row - lenscenter.y) / s*1.35 + lenscenter.y);
			//越界判断
			if (mCorrectPoint.y < 0 || mCorrectPoint.y >= img.rows-1)
			{
				continue;
			}
			if (mCorrectPoint.x < 0 || mCorrectPoint.x >= img.cols-1)
			{
				continue;
			}
			src_a = cvPoint((int)mCorrectPoint.x, (int)mCorrectPoint.y);
			src_b = cvPoint(src_a.x + 1, src_a.y);
			src_c = cvPoint(src_a.x, src_a.y + 1);
			src_d = cvPoint(src_a.x + 1, src_a.y + 1);
			distance_to_a_x = mCorrectPoint.x - src_a.x;//在原图像中与a点的水平距离    
			distance_to_a_y = mCorrectPoint.y - src_a.y;//在原图像中与a点的垂直距离    	
	
			drcimg.at<Vec3b>(row, cols)[0] =
				img.at<Vec3b>(src_a.y, src_a.x)[0] * (1 - distance_to_a_x)*(1 - distance_to_a_y) +
				img.at<Vec3b>(src_b.y, src_b.x)[0] * distance_to_a_x*(1 - distance_to_a_y) +
				img.at<Vec3b>(src_c.y, src_c.x)[0] * distance_to_a_y*(1 - distance_to_a_x) +
				img.at<Vec3b>(src_d.y, src_d.x)[0] * distance_to_a_y*distance_to_a_x;
			drcimg.at<Vec3b>(row, cols)[1] =
				img.at<Vec3b>(src_a.y, src_a.x)[1] * (1 - distance_to_a_x)*(1 - distance_to_a_y) +
				img.at<Vec3b>(src_b.y, src_b.x)[1] * distance_to_a_x*(1 - distance_to_a_y) +
				img.at<Vec3b>(src_c.y, src_c.x)[1] * distance_to_a_y*(1 - distance_to_a_x) +
				img.at<Vec3b>(src_d.y, src_d.x)[1] * distance_to_a_y*distance_to_a_x;
			drcimg.at<Vec3b>(row, cols)[2] =
				img.at<Vec3b>(src_a.y, src_a.x)[2] * (1 - distance_to_a_x)*(1 - distance_to_a_y) +
				img.at<Vec3b>(src_b.y, src_b.x)[2] * distance_to_a_x*(1 - distance_to_a_y) +
				img.at<Vec3b>(src_c.y, src_c.x)[2] * distance_to_a_y*(1 - distance_to_a_x) +
				img.at<Vec3b>(src_d.y, src_d.x)[2] * distance_to_a_y*distance_to_a_x;
		}
	printf("时间花费%dms\n",clock()-time);//修改前200左右,修改成2.0版本的变量之后在210左右
	cv::imwrite("矫正完成.bmp", drcimg);
	cv::imshow("矫正完成", drcimg);
	cv::waitKey(0);
	img.release();
	drcimg.release();
}

       处理结果如下图所示：

        由上图看到，图像的桶形畸变矫正已经完成了。由于桶形畸变越往外畸变的越严重，因此矫正过后的图像会比原图像大，而我们以图像的原尺寸来放置矫正后的图像，因此视角会比原图像小。这个问题如果需要处理的话可以通过缩小矫正后图像的大小等方法让图像有更大的视角，本文没做处理。
       另外，如果觉得这个方法矫正的不是很准确，可以利用Matlab的矫正工具箱对图像进行标定，得到图像的校准参数，然后导出xml文件，再用OpenCV读取这些导出的文件进行矫正。OpenCV的实现代码如下：

#include <opencv2/opencv.hpp>
using namespace cv;
 
int main()
{
	//从xml文件中读入校准参数，注意这里的xml文件可以是用matlab的校准工具箱生成的
	CvMat *Intrinsics_Camera = (CvMat*)cvLoad("Intrinsics_Camera.xml");
	CvMat *Distortion_Camera = (CvMat*)cvLoad("Distortion_Camera.xml");
	//读入图像
	IplImage *img = cvLoadImage("待带矫正图像.jpg");
	IplImage *img_Rect = cvCloneImage(img);
	//在校准参数的基础上创建映射矩阵
	CvMat *Mapx = cvCreateMat(img->height, img->width, CV_32F);
	CvMat *Mapy = cvCreateMat(img->height, img->width, CV_32F);
	cvInitUndistortMap(Intrinsics_Camera, Distortion_Camera, Mapx, Mapy);
	//矫正图像
	cvRemap(img, img_Rect, Mapx, Mapy);
	//显示
	cvNamedWindow("Image");
	cvShowImage("Image", img_Rect);
	cvSaveImage("矫正后图像.jpg", img_Rect);
	//等待
	cvWaitKey();
	//释放对象内存
	cvDestroyAllWindows();
	cvReleaseImage(&img);
	cvReleaseImage(&img_Rect);
	cvReleaseMat(&Intrinsics_Camera);
	cvReleaseMat(&Distortion_Camera);
	cvReleaseMat(&Mapx);
	cvReleaseMat(&Mapy);
	return 0;
}