初学OpenCV之摄像机标定
最近方向定下来是双目立体视觉，主要是做重建这块的研究。大致过程是图像获取->摄像机标定->特征提取->匹配->三维重建，当然开始可以进行图像预处理，矫正，后期可以进行点云的进一步处理，如渲染表面使其更接近于现实物体。
图像获取相对来说比较简单，用相机拍摄目标物（大型场景或特定小型的室内物体）。但有两点需要注意：
1、双目重建所需的图像一般为两张，角度相差不应过大，否则公共部分太少以至于重建效果不佳；整个过程简便，成本也不高，但缺陷是只有两张图像的点云所表示的物体信息不会很全面；
2、标定所需的图像又是另外拍摄的，用张正友标定法的话，把印有黑白棋盘格的图像粘至硬纸板上，然后左右摄像机各自进行拍摄，理论上获得角度（图像）越多，最终标定结果越精确；标定板见下图：
这里主要结合OpenCV对左右摄像机标定做一个简单的介绍，望朋友们指正，一起交流、进步。
摄像机的标定问题是机器视觉领域的入门问题，可以分为传统的摄像机定标方法和摄像机自定标方法。定标的方法有很多中常见的有：Tsai（传统）和张正友（介于传统和自定标）等，
摄像机成像模型和四个坐标系（通用原理）。
摄像机模型采用经典的小孔模型，如图中Oc（光心），像面π表示的是视野平面，其到光心的距离为f（镜头焦距）。
四个坐标系分别为：世界坐标系（Ow），摄像机坐标系（Oc），图像物理坐标系（O1，单位mm），图像像素坐标系（O，位于视野平面的左上角，单位pix）。
空间某点P到其像点p的坐标转换过程主要是通过这四套坐标系的三次转换实现的，首先将世界坐标系进行平移和转换得到摄像机坐标系，然后根据三角几何变换得到图像物理坐标系，最后根据像素和公制单位的比率得到图像像素坐标系。（实际的应用过程是这个的逆过程，即由像素长度获知实际的长度）。
ps：通过摄像头的标定，可以得到视野平面上的mm/pix分辨率，对于视野平面以外的物体还是需要通过坐标转换得到视野平面上。
转化的过程和公式参见：摄像机标定原理(关键是三个坐标系).ppt
2 张正友算法的原理
zhang法通过对一定标板在不同方向多次（三次以上）完整拍照，不需要知道定标板的运动方式。直接获得相机的内参（参考文献上矩阵A）和畸变系数。该标定方法精度高于自定标法，且不需要高精度的定位仪器。
ZHANG的算法包含两个模型：一.经典针孔模型，包含四个坐标系，二畸变模型（这个来源未知）
公式三项依次表示，径向畸变，切线畸变，薄棱镜畸变。OPENCV中函数只能给出k1,k2,p1,p2。
还存在另外一种畸变模型，见《摄像机标定算法库的设计和试验验证》一文26 page。
张正友标定有matlab的工具箱TOOLBOX_CAL,以及OpenCV库，下面是C++结合OpenCV的代码:
1 #include "cvut.h" 2 #include <iostream> 3 #include <fstream> 4 #include <string> 5 using namespace cvut; 6 using namespace std; 7 void main() 8 { 9 ifstream fin("calibdata.txt"); 10 ofstream fout("calibration_result.txt"); 11 //****************开始提取角点***********************// 12 cout<<"开始提取角点………………"; 13 int image_count=0; 14 CvSize image_size; 15 CvSize board_size = cvSize(5,7); 16 CvPoint2D32f * image_points_buf = 17 new CvPoint2D32f[board_size.width*board_size.height]; 18 Seq<CvPoint2D32f> image_points_seq; 19 string filename; 20 while (std::getline(fin,filename)) 21 { 22 cout<<"\n 将鼠标焦点移到标定图像所在窗口" 23 <<"并输入回车进行下一幅图像的角点提取 \n"; 24 image_count++; 25 int count; 26 Image<uchar> view(filename); 27 if (image_count == 1) 28 { 29 image_size.width = view.size().width; 30 image_size.height = view.size().height; 31 } 32 if (0 == cvFindChessboardCorners( view.cvimage, board_size, 33 image_points_buf, &count, CV_CALIB_CB_ADAPTIVE_THRESH )) 34 { 35 cout<<"can not find chessboard corners!\n"; 36 exit(1); 37 } 38 else 39 { 40 Image<uchar> view_gray(view.size(),8,1); 41 rgb2gray(view,view_gray); 42 cvFindCornerSubPix( view_gray.cvimage, 43 image_points_buf, count, cvSize(11,11), 44 cvSize(-1,-1), 45 cvTermCriteria( CV_TERMCRIT_EPS+CV_TERMCRIT_ITER, 30, 0.1 )); 46 image_points_seq.push_back(image_points_buf,count); 47 cvDrawChessboardCorners( view.cvimage, board_size, 48 image_points_buf, count, 1); 49 view.show("calib"); 50 cvWaitKey(); 51 view.close(); 52 } 53 }//角点提取循环 54 delete []image_points_buf; 55 cout<<"角点提取完成！\n"<<endl; 56 cout<<"开始定标………………"<<"\n"<<endl; 57 CvSize square_size = cvSize(10,10); 58 Matrix<double> object_points(1, 59 board_size.width*board_size.height*image_count,3); 60 Matrix<double> image_points(1,image_points_seq.cvseq->total,2); 61 Matrix<int> point_counts(1,image_count,1); 62 Matrix<double> intrinsic_matrix(3,3,1); 63 Matrix<double> distortion_coeffs(1,4,1); 64 Matrix<double> rotation_vectors(1,image_count,3); 65 Matrix<double> translation_vectors(1,image_count,3); 66 int i,j,t; 67 for (t=0;t<image_count;t++) { 68 for (i=0;i<board_size.height;i++) { 69 for (j=0;j<board_size.width;j++) { 70 object_points(0,t*board_size.height*board_size.width 71 + i*board_size.width + j,0) = i*square_size.width; 72 object_points(0,t*board_size.height*board_size.width 73 + i*board_size.width + j,1) = j*square_size.height; 74 object_points(0,t*board_size.height*board_size.width 75 + i*board_size.width + j,2) = 0; 76 } 77 } 78 } 79 char str[10]; 80 itoa(image_points_seq.cvseq->total,str,10); 81 cout<<str<<"\n"<<endl; 82 for (i=0;i<image_points_seq.cvseq->total;i++) 83 { 84 image_points(0,i,0) = image_points_seq[i].x; 85 image_points(0,i,1) = image_points_seq[i].y; 86 } 87 for (i=0;i<image_count;i++) 88 { 89 point_counts(0,i) = board_size.width*board_size.height; 90 } 91 cvCalibrateCamera2(object_points.cvmat, 92 image_points.cvmat, 93 point_counts.cvmat, 94 image_size, 95 intrinsic_matrix.cvmat, 96 distortion_coeffs.cvmat, 97 rotation_vectors.cvmat, 98 translation_vectors.cvmat, 99 0); 100 cout<<"定标完成！\n"; 101 cout<<"标定结果显示\n"; 102 cout<<"***************************************\n"; 103 cout<<"相机内参intrinsic_matrix\n"; 104 for(int h=0;h<3;h++) 105 { 106 cout<<"X:"<<intrinsic_matrix(h,0,0)<<"\tY:"<< 107 intrinsic_matrix(h,1,0)<<"\tZ:"<<intrinsic_matrix(h,2,0) 108 <<"\n"; 109 } 110 cout<<"\n畸变系数：distortion_coeffs\n"; 111 for(int ndis=0;ndis<4;ndis++) 112 { 113 cout<<distortion_coeffs(0,ndis,0)<<"\\"; 114 } 115 cout<<"\n"; 116 cout<<"\nrotation_vectors\n"; 117 for(int rot=0;rot<7;rot++) 118 { 119 cout<<"X:"<<rotation_vectors(0,rot,0)<<"\tY:" 120 <<rotation_vectors(0,rot,1)<<"\tZ:"<<rotation_vectors(0,rot,2) 121 <<"\n"; 122 } 123 cout<<"\ntranslation_vectors\n"; 124 for(i=0;i<7;i++) 125 { 126 cout<<"第"<<i+1<<"张图"<<"\tX:"<<translation_vectors(0,i,0) 127 <<"\tY:"<<translation_vectors(0,i,1) 128 <<"\tZ:"<<translation_vectors(0,i,2)<<"\n"; 129 } 130 cout<<"*********************\n"; 131 cout<<"开始评价定标结果………………\n"; 132 double total_err = 0.0; 133 double err = 0.0; 134 Matrix<double> image_points2(1,point_counts(0,0,0),2); 135 int temp_num = point_counts(0,0,0); 136 cout<<"\t每幅图像的定标误差：\n"; 137 fout<<"每幅图像的定标误差：\n"; 138 for (i=0;i<image_count;i++) 139 { 140 cvProjectPoints2(object_points.get_cols(i * 141 point_counts(0,0,0),(i+1)*point_counts(0,0,0)-1).cvmat, 142 rotation_vectors.get_col(i).cvmat, 143 translation_vectors.get_col(i).cvmat, 144 intrinsic_matrix.cvmat, 145 distortion_coeffs.cvmat, 146 image_points2.cvmat, 147 0,0,0,0); 148 err = cvNorm(image_points.get_cols(i*point_counts(0,0,0),(i+1) 149 *point_counts(0,0,0)-1).cvmat, 150 image_points2.cvmat, 151 CV_L1); 152 total_err += err/=point_counts(0,0,0); 153 cout<<"****************************\n"; 154 cout<<"\t\t第"<<i+1<<"幅图像的平均误差："<<err<<"像素"<<'\n'; 155 fout<<"\t第"<<i+1<<"幅图像的平均误差："<<err<<"像素"<<'\n'; 156 cout<<"显示image_point2\n"; 157 for(int ih=0;ih<7;ih++) 158 { 159 cout<<"X:"<<image_points2(0,ih,0)<<"\tY:" 160 <<image_points2(0,ih,1)<<"\n"; 161 } 162 cout<<"显示object_Points\n"; 163 for(int iw=0;iw<7;iw++) 164 { 165 cout<<"X:"<<image_points.get_cols(i*point_counts(0,0,0),(i+1) 166 *point_counts(0,0,0)-1)(0,iw,0) 167 <<"\tY:"<<image_points.get_cols(i*point_counts(0,0,0),(i+1) 168 *point_counts(0,0,0)-1)(0,iw,1)<<"\n"; 169 } 170 } 171 cout<<"\t总体平均误差："<<total_err/image_count<<"像素"<<'\n'; 172 fout<<"总体平均误差："<<total_err/image_count<<"像素"<<'\n'<<'\n'; 173 cout<<"评价完成！\n"; 174 cout<<"开始保存定标结果………………"; 175 Matrix<double> rotation_vector(3,1); 176 Matrix<double> rotation_matrix(3,3); 177 fout<<"相机内参数矩阵：\n"; 178 fout<<intrinsic_matrix<<'\n'; 179 fout<<"畸变系数：\n"; 180 fout<<distortion_coeffs<<'\n'; 181 for (i=0;i<image_count;i++) { 182 fout<<"第"<<i+1<<"幅图像的旋转向量：\n"; 183 fout<<rotation_vectors.get_col(i); 184 for (j=0;j<3;j++) 185 { 186 rotation_vector(j,0,0) = rotation_vectors(0,i,j); 187 } 188 cvRodrigues2(rotation_vector.cvmat,rotation_matrix.cvmat); 189 fout<<"第"<<i+1<<"幅图像的旋转矩阵：\n"; 190 fout<<rotation_matrix; 191 fout<<"第"<<i+1<<"幅图像的平移向量：\n"; 192 fout<<translation_vectors.get_col(i)<<'\n'; 193 } 194 cout<<"完成保存\n"; 195 }
这段程序可以直接运行，本人是在vs2010+OpenCV2.4.0环境下，配置在前面文章中有介绍。标定板图像序列名称放在文件"calibdata.txt"中，每行一幅图像名，路劲问题这里就不论述了。如
image/chess1.jpg
image/chess2.jpg
……
最后标定结果存储在"calibration_result.txt"中，命令行窗口也会显示一些信息。观察结果会发现，摄像机内参数矩阵M1和畸变系数k1,k2,p1,p2是直接调用OpenCV中标定函数求解出的，以及每幅图像相对于对应摄像机成像平面的旋转矩阵R'(或旋转向量)和平移向量t'，而不是世界坐标系相对于摄像机坐标系，也不是两部摄像机之间的空间关系。而对于外参[R t]，一直都是很抽象的问题。
在重建的过程中，主要需要两个元素：
1、摄像机内M1、外参数M2（或者说是投影矩阵M=M1*M2）
2、左右成像平面上的匹配特征点对
最后通过视差原理计算出物体上特征点对应的空间点坐标，得到最终点云。
当然，有兴趣的可以继续学习、研究下去，外参可以通过图像上的点和预先得到的世界坐标系下的物体上的空间点坐标，利用一系列数学关系计算。
在实际操作中，可以在重建过程进行自标定，即不用分别进行标定求得摄像机参数，而是直接通过特征点对计算投影矩阵，后面步骤同上。当然自标定得到的结果精度不是很高，目前较多的是上面所述的张正友标定法，介于传统与自标定之间。
