OpenCV cvEstimateRigidTransform函数详细注解

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 }

 

上一篇:C++入门Step03【字符串与指针】


下一篇:射频PA测试