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 }
