【特征匹配】RANSAC算法原理与源码解析

转载请注明出处：http://blog.csdn.net/luoshixian099/article/details/50217655

随机抽样一致性（RANSAC）算法，可以在一组包含“外点”的数据集中，采用不断迭代的方法，寻找最优参数模型，不符合最优模型的点，被定义为“外点”。在图像配准以及拼接上得到广泛的应用，本文将对RANSAC算法在OpenCV中角点误匹配对的检测中进行解析。

1.RANSAC原理

OpenCV中滤除误匹配对采用RANSAC算法寻找一个最佳单应性矩阵H，矩阵大小为3×3。RANSAC目的是找到最优的参数矩阵使得满足该矩阵的数据点个数最多，通常令h33=1来归一化矩阵。由于单应性矩阵有8个未知参数，至少需要8个线性方程求解，对应到点位置信息上，一组点对可以列出两个方程，则至少包含4组匹配点对。

                                                                   
   其中(x,y)表示目标图像角点位置,(x‘,y‘)为场景图像角点位置，s为尺度参数。

RANSAC算法从匹配数据集中随机抽出4个样本并保证这4个样本之间不共线，计算出单应性矩阵，然后利用这个模型测试所有数据，并计算满足这个模型数据点的个数与投影误差(即代价函数)，若此模型为最优模型，则对应的代价函数最小。

-----------------------------------------------------------------------------------------------------------------

RANSAC算法步骤：

1. 随机从数据集中随机抽出4个样本数据 (此4个样本之间不能共线)，计算出变换矩阵H，记为模型M；

2. 计算数据集中所有数据与模型M的投影误差，若误差小于阈值，加入内点集 I ；

3. 如果当前内点集 I 元素个数大于最优内点集 I_best , 则更新 I_best = I，同时更新迭代次数k ;

4. 如果迭代次数大于k,则退出 ; 否则迭代次数加1，并重复上述步骤；

注：迭代次数k在不大于最大迭代次数的情况下，是在不断更新而不是固定的；

其中，p为置信度，一般取0.995；w为"内点"的比例 ; m为计算模型所需要的最少样本数=4；

-----------------------------------------------------------------------------------------------------------------

2.例程

OpenCV中此功能通过调用findHomography函数调用，下面是个例程：

#include <iostream>
#include "opencv2/opencv.hpp"
#include "opencv2/core/core.hpp"
#include "opencv2/features2d/features2d.hpp"
#include "opencv2/highgui/highgui.hpp"
using namespace cv;
using namespace std;
int main(int argc, char** argv)
{
	Mat obj=imread("F:\\Picture\\obj.jpg");   //载入目标图像
	Mat scene=imread("F:\\Picture\\scene.jpg"); //载入场景图像
	if (obj.empty() || scene.empty() )
	{
		cout<<"Can't open the picture!\n";
		return 0;
	}
	vector<KeyPoint> obj_keypoints,scene_keypoints;
	Mat obj_descriptors,scene_descriptors;
   	ORB detector;     //采用ORB算法提取特征点
	detector.detect(obj,obj_keypoints);
	detector.detect(scene,scene_keypoints);
	detector.compute(obj,obj_keypoints,obj_descriptors);
	detector.compute(scene,scene_keypoints,scene_descriptors);
	BFMatcher matcher(NORM_HAMMING,true); //汉明距离做为相似度度量
	vector<DMatch> matches;
	matcher.match(obj_descriptors, scene_descriptors, matches);
	Mat match_img;
	drawMatches(obj,obj_keypoints,scene,scene_keypoints,matches,match_img);
	imshow("滤除误匹配前",match_img);

	//保存匹配对序号
	vector<int> queryIdxs( matches.size() ), trainIdxs( matches.size() );
	for( size_t i = 0; i < matches.size(); i++ )
	{
		queryIdxs[i] = matches[i].queryIdx;
		trainIdxs[i] = matches[i].trainIdx;
	}	

	Mat H12;   //变换矩阵

	vector<Point2f> points1; KeyPoint::convert(obj_keypoints, points1, queryIdxs);
	vector<Point2f> points2; KeyPoint::convert(scene_keypoints, points2, trainIdxs);
    int ransacReprojThreshold = 5;  //拒绝阈值

    H12 = findHomography( Mat(points1), Mat(points2), CV_RANSAC, ransacReprojThreshold );
    vector<char> matchesMask( matches.size(), 0 );
	Mat points1t;
	perspectiveTransform(Mat(points1), points1t, H12);
	for( size_t i1 = 0; i1 < points1.size(); i1++ )  //保存‘内点’
	{
		if( norm(points2[i1] - points1t.at<Point2f>((int)i1,0)) <= ransacReprojThreshold ) //给内点做标记
		{
			matchesMask[i1] = 1;
		}
	}
	Mat match_img2;   //滤除‘外点’后
	drawMatches(obj,obj_keypoints,scene,scene_keypoints,matches,match_img2,Scalar(0,0,255),Scalar::all(-1),matchesMask);

	//画出目标位置
	std::vector<Point2f> obj_corners(4);
	obj_corners[0] = cvPoint(0,0); obj_corners[1] = cvPoint( obj.cols, 0 );
	obj_corners[2] = cvPoint( obj.cols, obj.rows ); obj_corners[3] = cvPoint( 0, obj.rows );
	std::vector<Point2f> scene_corners(4);
	perspectiveTransform( obj_corners, scene_corners, H12);
    line( match_img2, scene_corners[0] + Point2f(static_cast<float>(obj.cols), 0),
		scene_corners[1] + Point2f(static_cast<float>(obj.cols), 0),Scalar(0,0,255),2);
	line( match_img2, scene_corners[1] + Point2f(static_cast<float>(obj.cols), 0),
		scene_corners[2] + Point2f(static_cast<float>(obj.cols), 0),Scalar(0,0,255),2);
	line( match_img2, scene_corners[2] + Point2f(static_cast<float>(obj.cols), 0),
		scene_corners[3] + Point2f(static_cast<float>(obj.cols), 0),Scalar(0,0,255),2);
	line( match_img2, scene_corners[3] + Point2f(static_cast<float>(obj.cols), 0),
		scene_corners[0] + Point2f(static_cast<float>(obj.cols), 0),Scalar(0,0,255),2);

    imshow("滤除误匹配后",match_img2);
	waitKey(0);

    return 0;
}

3. RANSAC源码解析

(1)findHomography内部调用cvFindHomography函数

srcPoints:目标图像点集

dstPoints:场景图像点集

method: 最小中值法、RANSAC方法、最小二乘法

ransacReprojTheshold:投影误差阈值

mask:掩码

cvFindHomography( const CvMat* objectPoints, const CvMat* imagePoints,
                  CvMat* __H, int method, double ransacReprojThreshold,
                  CvMat* mask )
{
    const double confidence = 0.995;  //置信度
    const int maxIters = 2000;    //迭代最大次数
    const double defaultRANSACReprojThreshold = 3; //默认拒绝阈值
    bool result = false;
    Ptr<CvMat> m, M, tempMask;

    double H[9];
    CvMat matH = cvMat( 3, 3, CV_64FC1, H );  //变换矩阵
    int count;

    CV_Assert( CV_IS_MAT(imagePoints) && CV_IS_MAT(objectPoints) );

    count = MAX(imagePoints->cols, imagePoints->rows);
    CV_Assert( count >= 4 );           //至少有4个样本
    if( ransacReprojThreshold <= 0 )
        ransacReprojThreshold = defaultRANSACReprojThreshold;

    m = cvCreateMat( 1, count, CV_64FC2 );  //转换为齐次坐标
    cvConvertPointsHomogeneous( imagePoints, m );

    M = cvCreateMat( 1, count, CV_64FC2 );//转换为齐次坐标
    cvConvertPointsHomogeneous( objectPoints, M );

    if( mask )
    {
        CV_Assert( CV_IS_MASK_ARR(mask) && CV_IS_MAT_CONT(mask->type) &&
            (mask->rows == 1 || mask->cols == 1) &&
            mask->rows*mask->cols == count );
    }
    if( mask || count > 4 )
        tempMask = cvCreateMat( 1, count, CV_8U );
    if( !tempMask.empty() )
        cvSet( tempMask, cvScalarAll(1.) );

    CvHomographyEstimator estimator(4);
    if( count == 4 )
        method = 0;
    if( method == CV_LMEDS )  //最小中值法
        result = estimator.runLMeDS( M, m, &matH, tempMask, confidence, maxIters );
    else if( method == CV_RANSAC )    //采用RANSAC算法
        result = estimator.runRANSAC( M, m, &matH, tempMask, ransacReprojThreshold, confidence, maxIters);//(2)解析
    else
        result = estimator.runKernel( M, m, &matH ) > 0; //最小二乘法

    if( result && count > 4 )
    {
        icvCompressPoints( (CvPoint2D64f*)M->data.ptr, tempMask->data.ptr, 1, count );  //保存内点集
        count = icvCompressPoints( (CvPoint2D64f*)m->data.ptr, tempMask->data.ptr, 1, count );
        M->cols = m->cols = count;
        if( method == CV_RANSAC )  //
            estimator.runKernel( M, m, &matH );  //内点集上采用最小二乘法重新估算变换矩阵
        estimator.refine( M, m, &matH, 10 );//
    }

    if( result )
        cvConvert( &matH, __H );  //保存变换矩阵

    if( mask && tempMask )
    {
        if( CV_ARE_SIZES_EQ(mask, tempMask) )
           cvCopy( tempMask, mask );
        else
           cvTranspose( tempMask, mask );
    }

    return (int)result;
}

(2) runRANSAC

maxIters:最大迭代次数

m1,m2 :数据样本

model:变换矩阵

reprojThreshold:投影误差阈值

confidence:置信度  0.995

bool CvModelEstimator2::runRANSAC( const CvMat* m1, const CvMat* m2, CvMat* model,
                                    CvMat* mask0, double reprojThreshold,
                                    double confidence, int maxIters )
{
    bool result = false;
    cv::Ptr<CvMat> mask = cvCloneMat(mask0);
    cv::Ptr<CvMat> models, err, tmask;
    cv::Ptr<CvMat> ms1, ms2;

    int iter, niters = maxIters;
    int count = m1->rows*m1->cols, maxGoodCount = 0;
    CV_Assert( CV_ARE_SIZES_EQ(m1, m2) && CV_ARE_SIZES_EQ(m1, mask) );

    if( count < modelPoints )
        return false;

    models = cvCreateMat( modelSize.height*maxBasicSolutions, modelSize.width, CV_64FC1 );
    err = cvCreateMat( 1, count, CV_32FC1 );//保存每组点对应的投影误差
    tmask = cvCreateMat( 1, count, CV_8UC1 );

    if( count > modelPoints )  //多于4个点
    {
        ms1 = cvCreateMat( 1, modelPoints, m1->type );
        ms2 = cvCreateMat( 1, modelPoints, m2->type );
    }
    else  //等于4个点
    {
        niters = 1; //迭代一次
        ms1 = cvCloneMat(m1);  //保存每次随机找到的样本点
        ms2 = cvCloneMat(m2);
    }

    for( iter = 0; iter < niters; iter++ )  //不断迭代
    {
        int i, goodCount, nmodels;
        if( count > modelPoints )
        {
           //在(3)解析getSubset
            bool found = getSubset( m1, m2, ms1, ms2, 300 ); //随机选择4组点，且三点之间不共线，(3)解析
            if( !found )
            {
                if( iter == 0 )
                    return false;
                break;
            }
        }

        nmodels = runKernel( ms1, ms2, models );  //估算近似变换矩阵，返回1
        if( nmodels <= 0 )
            continue;
        for( i = 0; i < nmodels; i++ )//执行一次
        {
            CvMat model_i;
            cvGetRows( models, &model_i, i*modelSize.height, (i+1)*modelSize.height );//获取3×3矩阵元素
            goodCount = findInliers( m1, m2, &model_i, err, tmask, reprojThreshold );  //找出内点，(4)解析

            if( goodCount > MAX(maxGoodCount, modelPoints-1) )  //当前内点集元素个数大于最优内点集元素个数
            {
                std::swap(tmask, mask);
                cvCopy( &model_i, model );  //更新最优模型
                maxGoodCount = goodCount;
				//confidence 为置信度默认0.995,modelPoints为最少所需样本数=4,niters迭代次数
                niters = cvRANSACUpdateNumIters( confidence,  //更新迭代次数，(5)解析
                    (double)(count - goodCount)/count, modelPoints, niters );
            }
        }
    }

    if( maxGoodCount > 0 )
    {
        if( mask != mask0 )
            cvCopy( mask, mask0 );
        result = true;
    }

    return result;
}

(3)getSubset

ms1,ms2:保存随机找到的4个样本

maxAttempts:最大迭代次数，300

bool CvModelEstimator2::getSubset( const CvMat* m1, const CvMat* m2,
                                   CvMat* ms1, CvMat* ms2, int maxAttempts )
{
    cv::AutoBuffer<int> _idx(modelPoints); //modelPoints 所需要最少的样本点个数
    int* idx = _idx;
    int i = 0, j, k, idx_i, iters = 0;
    int type = CV_MAT_TYPE(m1->type), elemSize = CV_ELEM_SIZE(type);
    const int *m1ptr = m1->data.i, *m2ptr = m2->data.i;
    int *ms1ptr = ms1->data.i, *ms2ptr = ms2->data.i;
    int count = m1->cols*m1->rows;

    assert( CV_IS_MAT_CONT(m1->type & m2->type) && (elemSize % sizeof(int) == 0) );
    elemSize /= sizeof(int); //每个数据占用字节数

    for(; iters < maxAttempts; iters++)
    {
        for( i = 0; i < modelPoints && iters < maxAttempts; )
        {
            idx[i] = idx_i = cvRandInt(&rng) % count;  // 随机选取1组点
            for( j = 0; j < i; j++ )  // 检测是否重复选择
                if( idx_i == idx[j] )
                    break;
            if( j < i )
                continue;  //重新选择
            for( k = 0; k < elemSize; k++ )  //拷贝点数据
            {
                ms1ptr[i*elemSize + k] = m1ptr[idx_i*elemSize + k];
                ms2ptr[i*elemSize + k] = m2ptr[idx_i*elemSize + k];
            }
            if( checkPartialSubsets && (!checkSubset( ms1, i+1 ) || !checkSubset( ms2, i+1 )))//检测点之间是否共线
            {
                iters++;  //若共线，重新选择一组
                continue;
            }
            i++;
        }
        if( !checkPartialSubsets && i == modelPoints &&
            (!checkSubset( ms1, i ) || !checkSubset( ms2, i )))
            continue;
        break;
    }

    return i == modelPoints && iters < maxAttempts;  //返回找到的样本点个数
}

(4) findInliers & computerReprojError

int CvModelEstimator2::findInliers( const CvMat* m1, const CvMat* m2,
                                    const CvMat* model, CvMat* _err,
                                    CvMat* _mask, double threshold )
{
    int i, count = _err->rows*_err->cols, goodCount = 0;
    const float* err = _err->data.fl;
    uchar* mask = _mask->data.ptr;

    computeReprojError( m1, m2, model, _err );  //计算每组点的投影误差
    threshold *= threshold;
    for( i = 0; i < count; i++ )
        goodCount += mask[i] = err[i] <= threshold;//误差在限定范围内，加入‘内点集’
    return goodCount;
}
//--------------------------------------
void CvHomographyEstimator::computeReprojError( const CvMat* m1, const CvMat* m2,const CvMat* model, CvMat* _err )
{
	int i, count = m1->rows*m1->cols;
	const CvPoint2D64f* M = (const CvPoint2D64f*)m1->data.ptr;
	const CvPoint2D64f* m = (const CvPoint2D64f*)m2->data.ptr;
	const double* H = model->data.db;
	float* err = _err->data.fl;

	for( i = 0; i < count; i++ )        //保存每组点的投影误差，对应上述变换公式
	{
		double ww = 1./(H[6]*M[i].x + H[7]*M[i].y + 1.);
		double dx = (H[0]*M[i].x + H[1]*M[i].y + H[2])*ww - m[i].x;
		double dy = (H[3]*M[i].x + H[4]*M[i].y + H[5])*ww - m[i].y;
		err[i] = (float)(dx*dx + dy*dy);
	}
}

(5)cvRANSACUpdateNumIters

对应上述k的计算公式

p：置信度

ep:外点比例

cvRANSACUpdateNumIters( double p, double ep,
                        int model_points, int max_iters )
{
    if( model_points <= 0 )
        CV_Error( CV_StsOutOfRange, "the number of model points should be positive" );

    p = MAX(p, 0.);
    p = MIN(p, 1.);
    ep = MAX(ep, 0.);
    ep = MIN(ep, 1.);

    // avoid inf's & nan's
    double num = MAX(1. - p, DBL_MIN);  //num=1-p,做分子
    double denom = 1. - pow(1. - ep,model_points);//做分母
    if( denom < DBL_MIN )
        return 0;

    num = log(num);
    denom = log(denom);

    return denom >= 0 || -num >= max_iters*(-denom) ?
        max_iters : cvRound(num/denom);
}
时间： 2024-10-27 11:56:04

【特征匹配】RANSAC算法原理与源码解析的相关文章

【特征匹配】SURF原理与源码解析

SURF (Speed Up Robust Features)是SIFT改进版也是加速版,提高了检测特征点的速度,综合性能要优于SIFT. 下面先逐次介绍SURF的原理,最后解析opencv上SURF源码. 转载请注明出处:http://blog.csdn.net/luoshixian099/article/details/47778143 1.积分图像 SURF是对积分图像进行操作,从而实现了加速,采用盒子滤波器计算每个像素点的Hessian矩阵行列式时,只需要几次加减法运算,而且运算量与盒子

【数据压缩】LZW算法原理与源码解析

转载请注明出处:http://blog.csdn.net/luoshixian099/article/details/50331883 <勿在浮沙筑高台> LZW压缩算法原理非常简单,因而被广泛地采用,已经被引入主流图像文件格式中.该算法由Lempel-Ziv-Welch三人发明,这种技术将定长码字分配给变长信源符号序列,它不需要知道被压缩文件的符号出现概率的先验知识,只需要动态地建立和维护一个字典,和其他压缩算法相比既是缺点也是优点. 1. LZW原理 1.1 概念的理解 LZW通过建立一个

【特征匹配】BRIEF特征描述子原理及源码解析

相关:Fast原理及源码解析 Harris原理及源码解析 SIFT原理及源码解析 SURF原理及源码解析 转载请注明出处: http://blog.csdn.net/luoshixian099/article/details/48338273 传统的特征点描述子如SIFT,SURF描述子,每个特征点采用128维(SIFT)或者64维(SURF)向量去描述,每个维度上占用4字节,SIFT需要128×4=512字节内存,SURF则需要256字节.如果对于内存资源有限的情况下,这种描述子方法显然不适应

OpenCV学习笔记（27）KAZE 算法原理与源码分析（一）非线性扩散滤波

http://blog.csdn.net/chenyusiyuan/article/details/8710462 OpenCV学习笔记(27)KAZE 算法原理与源码分析(一)非线性扩散滤波 2013-03-23 17:44 16963人阅读 评论(28) 收藏 举报 分类: 机器视觉(34) 版权声明:本文为博主原创文章,未经博主允许不得转载. 目录(?)[+] KAZE系列笔记: OpenCV学习笔记(27)KAZE 算法原理与源码分析(一)非线性扩散滤波 OpenCV学习笔记(28)KA

Sprig AOP原理及源码解析

在介绍AOP之前,想必很多人都听说AOP是基于动态代理和反射来实现的,那么在看AOP之前,你需要弄懂什么是动态代理和反射及它们又是如何实现的. 想了解JDK的动态代理及反射的实现和源码分析,请参见下面三篇文章 JDK的动态代理源码分析之一 (http://blog.csdn.net/weililansehudiefei/article/details/73655925) JDK的动态代理源码分析之二(http://blog.csdn.net/weililansehudiefei/article/

LinkedList原理及源码解析

简介 LinkedList是一个双向线性链表,但是并不会按线性的顺序存储数据,而是在每一个节点里存到下一个节点的指针(Pointer).由于不必须按顺序存储,链表在插入的时候可以达到O(1)的复杂度,比另一种线性表顺序表快得多,但是查找一个节点或者访问特定编号的节点则需要O(n)的时间,而顺序表相应的时间复杂度分别是O(logn)和O(1). UML关系图 使用示例 LinkedList list = new LinkedList<>(); //新增 list.add("a"

Spring核心框架 - AOP的原理及源码解析

一.AOP的体系结构 如下图所示:(引自AOP联盟) 层次3语言和开发环境:基础是指待增加对象或者目标对象:切面通常包括对于基础的增加应用:配置是指AOP体系中提供的配置环境或者编织配置,通过该配置AOP将基础和切面结合起来,从而完成切面对目标对象的编织实现. 层次2面向方面系统:配置模型,逻辑配置和AOP模型是为上策的语言和开发环境提供支持的,主要功能是将需要增强的目标对象.切面和配置使用AOP的API转换.抽象.封装成面向方面中的逻辑模型. 层次1底层编织实现模块:主要是将面向方面系统抽象封

【Spring】Spring IOC原理及源码解析之scope=request、session

一.容器 1. 容器 抛出一个议点:BeanFactory是IOC容器,而ApplicationContex则是Spring容器. 什么是容器?Collection和Container这两个单词都有存放什么东西的意思,但是放在程序猿的世界,却注定是千差万别.Collection,集合,存放obj instanceof Class为true的一类对象,重点在于存放:Container,容器,可以存放各种各样的obj,但不仅仅是存放,他被称为容器,更重要的是他能管理存放对象的生命周期和依赖. 容器:

可阻塞队列-原理及源码解析

阻塞原理:比如,一个队列中有8个格子,代表可放入8条数据,当一条信息到来就放入一个格子中,然后就进行处理.但是这个时候一次性来了8条数据,格子满了,数据还没有处理完,就来个一条数据.这个时候就把这条数据进行阻塞. 示例:假定有一个绑定的缓冲区,它支持 put 和 take 方法.如果试图在空的缓冲区上执行 take 操作,则在某一个项变得可用之前,线程将一直阻塞:如果试图在满的缓冲区上执行 put 操作,则在有空间变得可用之前,线程将一直阻塞.我们喜欢在单独的等待 set 中保存 put 线程和