RANSAC算法实现去除误匹配并计算拼接矩阵-随笔

1.RANSAC算法实现去除误匹配并计算拼接矩阵流程

(1) 从样本集中随机抽选一个RANSAC样本，即4个匹配点对（至少4个匹配点对，才能计算出3*3变换矩阵）；

(2) 计算当错误概率为0.1（即90%正确率）的情况下所需最少匹配点对数，赋值给in_min；

(3) 计算初始的错误概率p = pow( 1.0 - pow( in_frac, m ), k )（初始状态k=0，m为随机抽样匹配点对数通常为4）；

(4) srand( time(NULL) );//初始化随机数生成器

(5) while（p > p_badxform ）循环执行迭代，终止条件是实际计算错误概率小于预设p_badxform 错误概率，通常p_badxform =0.01（即99%正确率）；

(5-1) 分别计算当前点和匹配点的坐标（x,y）分别存入pts和mpts的Point类型点集中；

(5-2) 利用传入的lsq_homog函数根据pts、mpts和m来计算根据随机选取的四个匹配对所得变换矩阵；

(5-2-1) lsq_homog函数通过调用OpenCv的cvSolve函数来计算pts到mpts的变换矩阵，cvSolve详解请参考http://blog.csdn.net/u011028345/article/details/78402625；

(5-3) 若能成功计算出变换矩阵，执行(5-4)，若不能成功计算出变换矩阵，更新当前错误概率p（因为即使没有计算出变换矩阵，但迭代k值增加了，依旧要重新计算错误概率），执行(5-1),继续迭代随机选择4对匹配点；

(5-4) 给定特征点集，变换矩阵，误差函数，计算出当前一致集consensus，返回一致集中元素个数给in；

(5-5) 若当前一致集大于历史最优一致集，即当前一致集为最优，则更新最优一致集consensus_max；

(5-6) 更新当前错误概率p；

 重复执行(5)迭代，直到条件终止。

(6) 若最优一致集中元素个数大于最低标准，即符合要求，根据最优一致集中所有符合要求的匹配点对，利用lsq_homog函数中的cvSolve的最小二乘法求得变换矩阵。

RANSAC算法的迭代终止条件是当前错误概率p小于允许的最小错误概率p_badxform，一般传入的值为0.01，即要求正确率达到99%

当前错误概率p的计算公式为：p=( 1 - in_fracm)k，

其中in_frac是当前最优一致集中元素(内点)个数占样本总数的百分比，表征了当前最优一致集的优劣；

m是计算变换矩阵需要的最小特征点对个数，是固定值，一般是4；k是迭代次数。

所以，内点所占比例in_frac越大，错误概率越小；迭代次数k越大，错误概率越小，用来保证收敛。

2.RobHess的RANSAC源码分析

/*利用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和匹配点坐标数组mpts  
    CvMat* M = NULL;//当前变换矩阵  
    //p：当前计算出的模型的错误概率，当p小于p_badxform时迭代停止  
    //in_frac：内点数目占样本总数目的百分比  
    double p, in_frac = RANSAC_INLIER_FRAC_EST;  
    //nm：输入的特征点数组中具有mtype类型匹配点的特征点个数  
    //in：当前一致集中元素个数  
    //in_min：一致集中元素个数允许的最小值，保证RANSAC最终计算出的转换矩阵错误的概率大于0.1小于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的矩阵，返回变换矩阵给M  
        M = xform_fn( pts, mpts, m );//一般传入lsq_homog()函数  
        if( ! M )//出错判断  
            goto iteration_end;  
        //给定特征点集，变换矩阵，误差函数，计算出当前一致集consensus，返回一致集中元素个数给in  
        in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);  
  
        //若当前一致集大于历史最优一致集，即当前一致集为最优，则更新最优一致集consensus_max  
        if( 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的矩阵，返回变换矩阵给M  
        M = xform_fn( pts, mpts, in_max );  
        /***********下面会再进行一次迭代**********/  
        //根据变换矩阵M从样本集matched中计算出一致集consensus，返回一致集中元素个数给in  
        in = 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的矩阵，返回变换矩阵给M  
        M = 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  
}