cvEstimateRigidTransform是opencv中求取仿射變換的函數，定義在lkpyramid.cpp文件中，該函數先利用ransac算法從所有特征點中選取一定數目的特征點，選取出的這些特征點性質都較好，然后利用icvGetRTMatrix函數求取仿射變換系數，下面是cvEstimateRigidTransform函數的詳細注解。
1 CV_IMPL int
2 cvEstimateRigidTransform( const CvArr* matA, const CvArr* matB, CvMat* matM, int full_affine ) 3 { 4 const int COUNT = 15; 5 const int WIDTH = 160, HEIGHT = 120; 6 const int RANSAC_MAX_ITERS = 500; 7 const int RANSAC_SIZE0 = 3; 8 const double RANSAC_GOOD_RATIO = 0.5; 9
10 cv::Ptr<CvMat> sA, sB; //智能指針，相當於c++中的shared_ptr
11 cv::AutoBuffer<CvPoint2D32f> pA, pB; 12 cv::AutoBuffer<int> good_idx; 13 cv::AutoBuffer<char> status; 14 cv::Ptr<CvMat> gray; 15
16 CvMat stubA, *A = cvGetMat( matA, &stubA ); //將CvArr*類型的matA轉化為CvMat類型的stubA，A是1*192
17 CvMat stubB, *B = cvGetMat( matB, &stubB ); 18 CvSize sz0, sz1; 19 int cn, equal_sizes; 20 int i, j, k, k1; 21 int count_x, count_y, count = 0; 22 double scale = 1; 23 CvRNG rng = cvRNG(-1);//初始化隨機數發生器
24 double m[6]={0}; 25 CvMat M = cvMat( 2, 3, CV_64F, m ); 26 int good_count = 0; 27 CvRect brect; 28
29 if( !CV_IS_MAT(matM) ) 30 CV_Error( matM ? CV_StsBadArg : CV_StsNullPtr, "Output parameter M is not a valid matrix" ); 31
32 if( !CV_ARE_SIZES_EQ( A, B ) ) 33 CV_Error( CV_StsUnmatchedSizes, "Both input images must have the same size" ); 34
35 if( !CV_ARE_TYPES_EQ( A, B ) ) 36 CV_Error( CV_StsUnmatchedFormats, "Both input images must have the same data type" ); 37
38 if( CV_MAT_TYPE(A->type) == CV_8UC1 || CV_MAT_TYPE(A->type) == CV_8UC3 ) //8位無符號
39 { 40 cn = CV_MAT_CN(A->type); //返回通道數
41 sz0 = cvGetSize(A); 42 sz1 = cvSize(WIDTH, HEIGHT); //160,120
43
44 scale = MAX( (double)sz1.width/sz0.width, (double)sz1.height/sz0.height ); 45 scale = MIN( scale, 1. ); //scale需小於1
46 sz1.width = cvRound( sz0.width * scale ); //sz1的寬高比與原圖像的寬高比變得一致
47 sz1.height = cvRound( sz0.height * scale ); 48
49 equal_sizes = sz1.width == sz0.width && sz1.height == sz0.height; //如果equal_sizes=1，說明窗口sz1與原圖像sz0一樣大
50
51 if( !equal_sizes || cn != 1 ) //sz1與圖像大小不等或者通道數不為1
52 { 53 sA = cvCreateMat( sz1.height, sz1.width, CV_8UC1 ); 54 sB = cvCreateMat( sz1.height, sz1.width, CV_8UC1 ); 55
56 if( cn != 1 ) //通道數不為1
57 { 58 gray = cvCreateMat( sz0.height, sz0.width, CV_8UC1 ); 59 cvCvtColor( A, gray, CV_BGR2GRAY ); //先轉化成灰度圖
60 cvResize( gray, sA, CV_INTER_AREA ); //再改變圖像大小為160*120
61 cvCvtColor( B, gray, CV_BGR2GRAY ); 62 cvResize( gray, sB, CV_INTER_AREA ); 63 gray.release(); 64 } 65 else
66 { 67 cvResize( A, sA, CV_INTER_AREA ); //不管輸入圖像多大，進來之后都會被改成160*120大小
68 cvResize( B, sB, CV_INTER_AREA ); 69 } 70
71 A = sA; 72 B = sB; 73 } 74
75 count_y = COUNT; //15
76 count_x = cvRound((double)COUNT*sz1.width/sz1.height); 77 count = count_x * count_y; 78
79 pA.allocate(count); 80 pB.allocate(count); 81 status.allocate(count); 82
83 for( i = 0, k = 0; i < count_y; i++ ) 84 for( j = 0; j < count_x; j++, k++ ) 85 { 86 pA[k].x = (j+0.5f)*sz1.width/count_x; //初始化
87 pA[k].y = (i+0.5f)*sz1.height/count_y; 88 } 89
90 // find the corresponding points in B
91 cvCalcOpticalFlowPyrLK( A, B, 0, 0, pA, pB, count, cvSize(10,10), 3, 92 status, 0, cvTermCriteria(CV_TERMCRIT_ITER,40,0.1), 0 ); 93
94 // repack the remained points
95 for( i = 0, k = 0; i < count; i++ ) 96 if( status[i] ) // 需要保留的點
97 { 98 if( i > k ) 99 { 100 pA[k] = pA[i]; 101 pB[k] = pB[i]; 102 } 103 k++; 104 } 105
106 count = k; 107 } 108 else if( CV_MAT_TYPE(A->type) == CV_32FC2 || CV_MAT_TYPE(A->type) == CV_32SC2 ) 109 { 110 count = A->cols*A->rows; //A是CvMat*類型,上面有A = cvGetMat( matA, &stubA );
111 CvMat _pA, _pB; 112 pA.allocate(count); // pA, pB是AutoBuffer<CvPoint2D32f> 類型
113 pB.allocate(count); 114 _pA = cvMat( A->rows, A->cols, CV_32FC2, pA ); //注意這里CV_32FC2是兩個通道
115 _pB = cvMat( B->rows, B->cols, CV_32FC2, pB ); 116 cvConvert( A, &_pA ); //#define cvConvert(src, dst ) cvConvertScale((src), (dst), 1, 0 )
117 cvConvert( B, &_pB ); 118 } 119 else
120 CV_Error( CV_StsUnsupportedFormat, "Both input images must have either 8uC1 or 8uC3 type" ); 121
122 good_idx.allocate(count); 123
124 if( count < RANSAC_SIZE0 ) 125 return 0; 126
127 CvMat _pB = cvMat(1, count, CV_32FC2, pB); 128 brect = cvBoundingRect(&_pB, 1); 129
130 // RANSAC stuff: 131 // 1. find the consensus
132 for( k = 0; k < RANSAC_MAX_ITERS; k++ ) //如果中途出現無法選到足夠的點等情況，則重新開始新一輪選點過程，因此這里有個循環
133 { 134 int idx[RANSAC_SIZE0]; 135 CvPoint2D32f a[3]; 136 CvPoint2D32f b[3]; 137
138 memset( a, 0, sizeof(a) ); // 將a所指向的某一塊內存中的每個字節的內容全部設置為0， 塊的大小由第三個參數指定，這個函數通常為新申請的內存做初始化工作， 其返回值為指向S的指針。
139 memset( b, 0, sizeof(b) ); 140
141 // choose random 3 non-complanar points from A & B
142 for( i = 0; i < RANSAC_SIZE0; i++ ) //每個點
143 { 144 for( k1 = 0; k1 < RANSAC_MAX_ITERS; k1++ ) //每次選取當前點的迭代次數
145 { 146 idx[i] = cvRandInt(&rng) % count; //從所有特征點中隨機抽一個點的索引
147
148 for( j = 0; j < i; j++ ) //前面已經抽好的點
149 { 150 if( idx[j] == idx[i] ) 151 break; 152 // check that the points are not very close one each other
153 if( fabs(pA[idx[i]].x - pA[idx[j]].x) +
154 fabs(pA[idx[i]].y - pA[idx[j]].y) < FLT_EPSILON ) 155 break; 156 if( fabs(pB[idx[i]].x - pB[idx[j]].x) +
157 fabs(pB[idx[i]].y - pB[idx[j]].y) < FLT_EPSILON ) 158 break; 159 } 160
161 if( j < i ) //是從上面的break跳出來的
162 continue;//當前選取的點不行，結束當前點此次的迭代
163
164 if( i+1 == RANSAC_SIZE0 ) //最后一個點
165 { 166 // additional check for non-complanar vectors不共線
167 a[0] = pA[idx[0]]; 168 a[1] = pA[idx[1]]; 169 a[2] = pA[idx[2]]; 170
171 b[0] = pB[idx[0]]; 172 b[1] = pB[idx[1]]; 173 b[2] = pB[idx[2]]; 174
175 double dax1 = a[1].x - a[0].x, day1 = a[1].y - a[0].y; 176 double dax2 = a[2].x - a[0].x, day2 = a[2].y - a[0].y; 177 double dbx1 = b[1].x - b[0].x, dby1 = b[1].y - b[0].y; 178 double dbx2 = b[2].x - b[0].x, dby2 = b[2].y - b[0].y; 179 const double eps = 0.01; 180
181 if( fabs(dax1*day2 - day1*dax2) < eps*sqrt(dax1*dax1+day1*day1)*sqrt(dax2*dax2+day2*day2) ||
182 fabs(dbx1*dby2 - dby1*dbx2) < eps*sqrt(dbx1*dbx1+dby1*dby1)*sqrt(dbx2*dbx2+dby2*dby2) ) 183 continue; 184 } 185 break; //程序能運行到這里說明上面對當前點的要求均滿足，因此當前點可用，不需再迭代尋找當前點
186 } //當前點的一次迭代結束
187
188 if( k1 >= RANSAC_MAX_ITERS ) //說明迭代了RANSAC_MAX_ITERS次都沒找到合適的第i個點
189 break; //不再繼續往后找第i+1,i+2,i+3個點，而是准備新一輪的找點，即重新找第0,1,2,3....個點
190 } //當前第i個點結束
191
192 if( i < RANSAC_SIZE0 ) //如果從if( k1 >= RANSAC_MAX_ITERS )跳出循環，即沒有找到足夠多的點，則會執行此句
193 continue; //跳出當前的第k次迭代，准備第k+1輪迭代，即重新找第0,1,2,3....個點 194
195 // estimate the transformation using 3 points
196 icvGetRTMatrix( a, b, 3, &M, full_affine ); //函數定義在lkpyramid.cpp中，如果能執行到這里，說明找到了足夠多的符合條件的點
197
198 for( i = 0, good_count = 0; i < count; i++ ) //count是所有角點的總個數
199 { 200 if( fabs( m[0]*pA[i].x + m[1]*pA[i].y + m[2] - pB[i].x ) +
201 fabs( m[3]*pA[i].x + m[4]*pA[i].y + m[5] - pB[i].y ) < MAX(brect.width,brect.height)*0.05 ) 202 good_idx[good_count++] = i; 203 } 204
205 if( good_count >= count*RANSAC_GOOD_RATIO ) //如果第k次迭代找到的點能很好的代表所有點，則break不再迭代
206 break; 207 } //第k次迭代結束
208
209 if( k >= RANSAC_MAX_ITERS ) //所有的迭代結束都沒找到合適的一組的點
210 return 0; //此時直接返回，M中保留的是最后一次改寫后的結果或者為全0(如果最外層的RANSAC_MAX_ITERS次迭代每次都從if( i < RANSAC_SIZE0 )行跳出循環的話)
211
212 if( good_count < count ) //如果執行這句，則說明k < RANSAC_MAX_ITERS
213 { 214 for( i = 0; i < good_count; i++ ) 215 { 216 j = good_idx[i]; 217 pA[i] = pA[j]; 218 pB[i] = pB[j]; 219 } 220 } 221
222 icvGetRTMatrix( pA, pB, good_count, &M, full_affine ); 223 m[2] /= scale; 224 m[5] /= scale; 225 cvConvert( &M, matM ); 226
227 return 1; 228 }