利用RANSAC算法筛选SIFT特征匹配

关于RANSAC算法的基本思想，可从网上搜索找到，这里只是RANSAC用于SIFT特征匹配筛选时的一些说明。

RANSAC算法在SIFT特征筛选中的主要流程是：

(1) 从样本集中随机抽选一个RANSAC样本，即4个匹配点对

(2) 根据这4个匹配点对计算变换矩阵M

(3) 根据样本集，变换矩阵M，和误差度量函数计算满足当前变换矩阵的一致集consensus，并返回一致集中元素个数

(4) 根据当前一致集中元素个数判断是否最优(最大)一致集，若是则更新当前最优一致集

(5) 更新当前错误概率p，若p大于允许的最小错误概率则重复(1)至(4)继续迭代，直到当前错误概率p小于最小错误概率

下面结合RobHess的源码说明一下RANSAC算法在SIFT特征匹配筛选中的实现，

具体的源码分析见此系列文章：RobHess的SIFT源码分析：综述

在RobHess的源码中，RANSAC算法的声明和实现在xform.h和xform.c文件中

实现RANSAC算法的主函数是ransac_xform，如下：

CvMat* ransac_xform( struct feature* features, int n, int mtype,ransac_xform_fn xform_fn, int m, double p_badxform,ransac_err_fn err_fn, double err_tol,struct feature*** inliers, int* n_in )

函数说明：利用RANSAC算法进行特征点筛选，计算出最佳匹配的变换矩阵
参数：
features：特征点数组，只有当mtype类型的匹配点存在时才被用来进行单应性计算
n：特征点个数
mtype：决定使用每个特征点的哪个匹配域进行变换矩阵的计算，应该是FEATURE_MDL_MATCH，
FEATURE_BCK_MATCH，FEATURE_MDL_MATCH中的一个。若是FEATURE_MDL_MATCH，
对应的匹配点对坐标是每个特征点的img_pt域和其匹配点的mdl_pt域，
否则，对应的匹配点对是每个特征点的img_pt域和其匹配点的img_pt域。
xform_fn：函数指针，指向根据输入的点对进行变换矩阵计算的函数，一般传入lsq_homog()函数
m：在函数xform_fn中计算变换矩阵需要的最小特征点对个数
p_badxform：允许的错误概率，即允许RANSAC算法计算出的变换矩阵错误的概率，当前计算出的模型的错误概率小于p_badxform时迭代停止
err_fn：函数指针，对于给定的变换矩阵，计算推定的匹配点对之间的变换误差，一般传入homog_xfer_err()函数
err_tol：容错度，对于给定的变换矩阵，在此范围内的点对被认为是内点
inliers：输出参数，指针数组，指向计算出的最终的内点集合，若为空，表示没计算出符合要求的一致集
此数组的内存将在本函数中被分配，使用完后必须在调用出释放：free(*inliers)
n_in：输出参数，最终计算出的内点的数目
返回值：RANSAC算法计算出的变换矩阵，若为空，表示出错或无法计算出可接受矩阵

注释如下：

CvMat* ransac_xform( struct feature* features, int n, int mtype,ransac_xform_fn xform_fn, int m, double p_badxform,ransac_err_fn err_fn, double err_tol,struct feature*** inliers, int* n_in )
{//matched：所有具有mtype类型匹配点的特征点的数组，也就是样本集//sample：单个样本，即4个特征点的数组//consensus：当前一致集；//consensus_max：当前最大一致集(即当前的最好结果的一致集)struct feature** matched, ** sample, ** consensus, ** consensus_max = NULL;struct ransac_data* rdata;//每个特征点的feature_data域的ransac数据指针CvPoint2D64f* pts, * mpts;//每个样本对应的两个坐标数组：特征点坐标数组pts和匹配点坐标数组mptsCvMat* M = NULL;//当前变换矩阵//p：当前计算出的模型的错误概率，当p小于p_badxform时迭代停止//in_frac：内点数目占样本总数目的百分比double p, in_frac = RANSAC_INLIER_FRAC_EST;//nm：输入的特征点数组中具有mtype类型匹配点的特征点个数//in：当前一致集中元素个数//in_min：一致集中元素个数允许的最小值，保证RANSAC最终计算出的转换矩阵错误的概率小于p_badxform所需的最小内点数目//in_max：当前最优一致集(最大一致集)中元素的个数//k：迭代次数，与计算当前模型的错误概率有关int i, nm, in, in_min, in_max = 0, k = 0;//找到特征点数组features中所有具有mtype类型匹配点的特征点，放到matched数组(样本集)中，返回值nm是matched数组的元素个数nm = get_matched_features( features, n, mtype, &matched );//若找到的具有匹配点的特征点个数小于计算变换矩阵需要的最小特征点对个数，出错if( nm < m ){   //出错处理，特征点对个数不足fprintf( stderr, "Warning: not enough matches to compute xform, %s" \" line %d\n", __FILE__, __LINE__ );goto end;}/* initialize random number generator */srand( time(NULL) );//初始化随机数生成器//计算保证RANSAC最终计算出的转换矩阵错误的概率小于p_badxform所需的最小内点数目in_min = calc_min_inliers( nm, m, RANSAC_PROB_BAD_SUPP, p_badxform );//当前计算出的模型的错误概率,内点所占比例in_frac越大，错误概率越小；迭代次数k越大，错误概率越小p = pow( 1.0 - pow( in_frac, m ), k );i = 0;//当前错误概率大于输入的允许错误概率p_badxform，继续迭代while( p > p_badxform ){//从样本集matched中随机抽选一个RANSAC样本(即一个4个特征点的数组)，放到样本变量sample中sample = draw_ransac_sample( matched, nm, m );//从样本中获取特征点和其对应匹配点的二维坐标，分别放到输出参数pts和mpts中extract_corresp_pts( sample, m, mtype, &pts, &mpts );//调用参数中传入的函数xform_fn，计算将m个点的数组pts变换为mpts的矩阵，返回变换矩阵给MM = xform_fn( pts, mpts, m );//一般传入lsq_homog()函数if( ! M )//出错判断goto iteration_end;//给定特征点集，变换矩阵，误差函数，计算出当前一致集consensus，返回一致集中元素个数给inin = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);//若当前一致集大于历史最优一致集，即当前一致集为最优，则更新最优一致集consensus_maxif( in > in_max ){if( consensus_max )//若之前有最优值，释放其空间free( consensus_max );consensus_max = consensus;//更新最优一致集in_max = in;//更新最优一致集中元素个数in_frac = (double)in_max / nm;//最优一致集中元素个数占样本总个数的百分比}else//若当前一致集小于历史最优一致集，释放当前一致集free( consensus );cvReleaseMat( &M );iteration_end:release_mem( pts, mpts, sample );p = pow( 1.0 - pow( in_frac, m ), ++k );}//根据最优一致集计算最终的变换矩阵/* calculate final transform based on best consensus set *///若最优一致集中元素个数大于最低标准，即符合要求if( in_max >= in_min ){//从最优一致集中获取特征点和其对应匹配点的二维坐标，分别放到输出参数pts和mpts中extract_corresp_pts( consensus_max, in_max, mtype, &pts, &mpts );//调用参数中传入的函数xform_fn，计算将in_max个点的数组pts变换为mpts的矩阵，返回变换矩阵给MM = xform_fn( pts, mpts, in_max );/***********下面会再进行一次迭代**********///根据变换矩阵M从样本集matched中计算出一致集consensus，返回一致集中元素个数给inin = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);cvReleaseMat( &M );release_mem( pts, mpts, consensus_max );//从一致集中获取特征点和其对应匹配点的二维坐标，分别放到输出参数pts和mpts中extract_corresp_pts( consensus, in, mtype, &pts, &mpts );//调用参数中传入的函数xform_fn，计算将in个点的数组pts变换为mpts的矩阵，返回变换矩阵给MM = xform_fn( pts, mpts, in );if( inliers ){*inliers = consensus;//将最优一致集赋值给输出参数：inliers，即内点集合consensus = NULL;}if( n_in )*n_in = in;//将最优一致集中元素个数赋值给输出参数：n_in，即内点个数release_mem( pts, mpts, consensus );}else if( consensus_max ){   //没有计算出符合要求的一致集if( inliers )*inliers = NULL;if( n_in )*n_in = 0;free( consensus_max );}//RANSAC算法结束：恢复特征点中被更改的数据域feature_data，并返回变换矩阵M
end:for( i = 0; i < nm; i++ ){//利用宏feat_ransac_data来提取matched[i]中的feature_data成员并转换为ransac_data格式的指针rdata = feat_ransac_data( matched[i] );//恢复feature_data成员的以前的值matched[i]->feature_data = rdata->orig_feat_data;free( rdata );//释放内存}free( matched );return M;//返回求出的变换矩阵M
}

实验测试：

//特征提取和匹配
void match(IplImage *img1,IplImage *img2)
{IplImage *img1_Feat = cvCloneImage(img1);//复制图1，深拷贝，用来画特征点IplImage *img2_Feat = cvCloneImage(img2);//复制图2，深拷贝，用来画特征点struct feature *feat1, *feat2;//feat1：图1的特征点数组，feat2：图2的特征点数组int n1, n2;//n1:图1中的特征点个数，n2：图2中的特征点个数struct feature *feat;//每个特征点struct kd_node *kd_root;//k-d树的树根struct feature **nbrs;//当前特征点的最近邻点数组int matchNum;//经距离比值法筛选后的匹配点对的个数struct feature **inliers;//精RANSAC筛选后的内点数组int n_inliers;//经RANSAC算法筛选后的内点个数,即feat1中具有符合要求的特征点的个数//默认提取的是LOWE格式的SIFT特征点//提取并显示第1幅图片上的特征点n1 = sift_features( img1, &feat1 );//检测图1中的SIFT特征点,n1是图1的特征点个数export_features("feature1.txt",feat1,n1);//将特征向量数据写入到文件draw_features( img1_Feat, feat1, n1 );//画出特征点cvShowImage("img1_Feat",img1_Feat);//显示//提取并显示第2幅图片上的特征点n2 = sift_features( img2, &feat2 );//检测图2中的SIFT特征点，n2是图2的特征点个数export_features("feature2.txt",feat2,n2);//将特征向量数据写入到文件draw_features( img2_Feat, feat2, n2 );//画出特征点cvShowImage("img2_Feat",img2_Feat);//显示Point pt1,pt2;//连线的两个端点double d0,d1;//feat1中每个特征点到最近邻和次近邻的距离matchNum = 0;//经距离比值法筛选后的匹配点对的个数//将2幅图片合成1幅图片,上下排列stacked = stack_imgs( img1, img2 );//合成图像，显示经距离比值法筛选后的匹配结果stacked_ransac = stack_imgs( img1, img2 );//合成图像，显示经RANSAC算法筛选后的匹配结果//根据图2的特征点集feat2建立k-d树，返回k-d树根给kd_rootkd_root = kdtree_build( feat2, n2 );//遍历特征点集feat1，针对feat1中每个特征点feat，选取符合距离比值条件的匹配点，放到feat的fwd_match域中for(int i = 0; i < n1; i++ ){feat = feat1+i;//第i个特征点的指针//在kd_root中搜索目标点feat的2个最近邻点，存放在nbrs中，返回实际找到的近邻点个数int k = kdtree_bbf_knn( kd_root, feat, 2, &nbrs, KDTREE_BBF_MAX_NN_CHKS );if( k == 2 ){d0 = descr_dist_sq( feat, nbrs[0] );//feat与最近邻点的距离的平方d1 = descr_dist_sq( feat, nbrs[1] );//feat与次近邻点的距离的平方//若d0和d1的比值小于阈值NN_SQ_DIST_RATIO_THR，则接受此匹配，否则剔除if( d0 < d1 * NN_SQ_DIST_RATIO_THR ){   //将目标点feat和最近邻点作为匹配点对pt1 = Point( cvRound( feat->x ), cvRound( feat->y ) );//图1中点的坐标pt2 = Point( cvRound( nbrs[0]->x ), cvRound( nbrs[0]->y ) );//图2中点的坐标(feat的最近邻点)pt2.y += img1->height;//由于两幅图是上下排列的，pt2的纵坐标加上图1的高度，作为连线的终点cvLine( stacked, pt1, pt2, CV_RGB(255,0,255), 1, 8, 0 );//画出连线matchNum++;//统计匹配点对的个数feat1[i].fwd_match = nbrs[0];//使点feat的fwd_match域指向其对应的匹配点}}free( nbrs );//释放近邻数组}qDebug()<<tr("经距离比值法筛选后的匹配点对个数：")<<matchNum<<endl;//显示并保存经距离比值法筛选后的匹配图cvNamedWindow(IMG_MATCH1);//创建窗口cvShowImage(IMG_MATCH1,stacked);//显示//利用RANSAC算法筛选匹配点,计算变换矩阵HCvMat * H = ransac_xform(feat1,n1,FEATURE_FWD_MATCH,lsq_homog,4,0.01,homog_xfer_err,3.0,&inliers,&n_inliers);qDebug()<<tr("经RANSAC算法筛选后的匹配点对个数：")<<n_inliers<<endl;//遍历经RANSAC算法筛选后的特征点集合inliers，找到每个特征点的匹配点，画出连线for(int i=0; i<n_inliers; i++){feat = inliers[i];//第i个特征点pt1 = Point(cvRound(feat->x), cvRound(feat->y));//图1中点的坐标pt2 = Point(cvRound(feat->fwd_match->x), cvRound(feat->fwd_match->y));//图2中点的坐标(feat的匹配点)qDebug()<<"("<<pt1.x<<","<<pt1.y<<")--->("<<pt2.x<<","<<pt2.y<<")"<<endl;pt2.y += img1->height;//由于两幅图是上下排列的，pt2的纵坐标加上图1的高度，作为连线的终点cvLine(stacked_ransac,pt1,pt2,CV_RGB(255,0,255),1,8,0);//画出连线}cvNamedWindow(IMG_MATCH2);//创建窗口cvShowImage(IMG_MATCH2,stacked_ransac);//显示
}

结果：