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 }