目录

1.RANSAC原理  

2. RANSAC算法步骤: 

3. RANSAC源码解析

step one niters最初的值为2000,这就是初始时的RANSAC算法的循环次数,getSubset()函数是从一组对应的序列中随机的选出4组(因为要想计算出一个3X3的矩阵,至少需要4组对应的坐标),m1和m2是我们输入序列,ms1和ms2是随机选出的对应的4组匹配。

参考https://blog.csdn.net/laobai1015/article/details/51683076

cv::findFundamentalMat

4. 针孔相机模型和变形

5.照相机定标

5.1 ProjectPoints2

5.2 FindHomographyFindHomography

5.3 CalibrateCamera2

5.4 FindExtrinsicCameraParams2

计算指定视图的摄像机外参数

5.5 Rodrigues2

5.6 Undistort2

5.7 InitUndistortMap

5. 8 FindChessboardCorners

DrawChessBoardCorners

6. 姿态估计

6.1 CreatePOSITObject

6.2 POSIT

6.3 ReleasePOSITObject

6.4 CalcImageHomography

6.5 FindFundamentalMat

6.6 ComputeCorrespondEpilines

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

1.RANSAC原理  

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

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

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

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

2. 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函数调用,下面是个例程:

[cpp]  view plain  copy

#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:掩码

[cpp]  view plain  copy

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;  
}  

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

[cpp]  view plain  copy

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;  
}  

step one niters最初的值为2000,这就是初始时的RANSAC算法的循环次数,getSubset()函数是从一组对应的序列中随机的选出4组(因为要想计算出一个3X3的矩阵,至少需要4组对应的坐标),m1和m2是我们输入序列,ms1和ms2是随机选出的对应的4组匹配。
 随机的选出4组匹配后,就应该根据这4个匹配计算相应的矩阵,所以函数runKernel()就是根据4组匹配计算矩阵,参数里的models就是得到的矩阵。

stept wo 这个矩阵只是一个假设,为了验证这个假设,需要用其他的点去计算,看看其他的点是内点还是外点。

findInliers()函数就是用来计算内点的。利用前面得到的矩阵,把所有的序列带入,计算得出哪些是内点,哪些是外点,函数的返回值为goodCount,就是此次计算的内点的个数。

step three  函数中还有一个值为maxGoodCount,每次循环的内点个数的最大值保存在这个值中,一个估计的矩阵如果有越多的内点,那么这个矩阵就越有可能是正确的。

step four   所以计算内点个数以后,紧接着判断一下goodCount和maxGoodCount的大小关系,如果goodCount>maxGoodCount,则把goodCount赋值给maxGoodCount。赋值之后的一行代码非常关键,我们单独拿出来说一下,代码如下:

niters = cvRANSACUpdateNumIters( confidence,  
                   (double)(count - goodCount)/count, modelPoints, niters );

 
niters本来是迭代的次数,也就是循环的次数。但是通过这行代码我们发现,每次循环后,都会对niters这个值进行更新,也就是每次循环后都会改变循环的总次数。

cvRANSACUpdateNumIters()函数利用confidence(置信度)count(总匹配个数)goodCount(当前内点个数)niters(当前的总迭代次数)这几个参数,来动态的改变总迭代次数的大小。该函数的中心思想就是当内点占的比例较多时,那么很有可能已经找到了正确的估计,所以就适当的减少迭代次数来节省时间。这个迭代次数的减少是以指数形式减少的,所以节省的时间开销也是非常的可观。因此最初设计的2000的迭代次数,可能最终的迭代次数只有几十。同样的,如果你自己一开始把迭代次数设置成10000或者更大,进过几次迭代后,niters又会变得非常小了。所以初始时的niters设置的再大,其实对最终的运行时间也没什么影响。

所以,们现在应该清楚为什么输入数据增加,而算法运行时间不会增加了。openCV的RANSAC算法首先把迭代的次数设置为2000,然后再迭代的过程中,动态的改变总迭代次数,无论输入数据有多少,总的迭代次数不会增加,并且通过4个匹配计算出估计的矩阵这个时间是不变的,通过估计矩阵来计算内点,这方面的增加的时间开销基本上可以忽略。所以导致的最终结果就是,无论输入点有多少,运算时间基本不会有太大变化。

以上就是RANSAC算法的核心代码,其中用到的一些函数,下面一一给出。

(3)getSubset

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

maxAttempts:最大迭代次数,300

[cpp]  view plain  copy
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

[cpp]  view plain  copy

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:外点比例

[cpp]  view plain  copy
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);  
}  

参考https://blog.csdn.net/laobai1015/article/details/51683076

cv::findFundamentalMat

在处理立体图像对的时候经常会用到对极几何的知识,计算基础矩阵也是很常见的事。OpenCV实现了基本矩阵的算法。对于老版本的C代码,计算基本矩阵的RANSAC方法中有一个迭代次数不收敛的bug,可能导致每次计算的采样次数都达到最大限制的次数,其根源在于计算采样可信度的公式有误,新版本的代码还没有仔细看过,不敢确定是否已经修正了这个bug。但是这个bug并不会对最后的结果造成多大影响,只是会影响计算的效率——原本几次采样即可结束迭代,在这个bug的影响下可能要采样数百次。

    新版本的计算基本矩阵的函数的使用也有一些问题,下面来看看cv::findFundamentalMat函数:

//! the algorithm for finding fundamental matrix
enum
{
    FM_7POINT = CV_FM_7POINT, //!< 7-point algorithm
    FM_8POINT = CV_FM_8POINT, //!< 8-point algorithm
    FM_LMEDS = CV_FM_LMEDS,  //!< least-median algorithm
    FM_RANSAC = CV_FM_RANSAC  //!< RANSAC algorithm
};

//! finds fundamental matrix from a set of corresponding 2D points
CV_EXPORTS Mat findFundamentalMat( const Mat& points1, const Mat& points2,
                                     CV_OUT vector<uchar>& mask, int method=FM_RANSAC,
                                     double param1=3., double param2=0.99 );

//! finds fundamental matrix from a set of corresponding 2D points
CV_EXPORTS_W Mat findFundamentalMat( const Mat& points1, const Mat& points2,
                                     int method=FM_RANSAC,
                                     double param1=3., double param2=0.99 );

上面是OpenCV计算基本矩阵的C++接口,其内部实现仍然是调用的C代码函数,仅仅是在上层进行了一次封装。网上有些文档中说const Mat& points1和points2这两个参数可以直接由vector<Point2f>类型作为传入参数,其实这是错误的。直接用vector<Point2f>传递,在编译时可能不会报错,但是运行时会直接崩溃。因为cv::Mat的构造函数并不会把Point2f转化为两个浮点数存于Mat的两个元素中,而是仍然以Point2f存于Mat的一个元素中,于是findFundamentalMat一读取这个Mat就会出错。

    因此我们最好老老实实地去构建Mat points1和Mat points2。该矩阵的维度可以是2xN,也可以是Nx2的,其中N是特征点的数目。另一个要注意的地方就是该矩阵的类型不能是CV_64F,也就是说findFundamentalMat内部解析points1和points2时是按照float类型去解析的,设为CV_64F将导致读取数据失败,程序崩溃。最好设为CV_32F。下面是从匹配好的特征点计算F的代码示例:

// vector<KeyPoint> m_LeftKey;
// vector<KeyPoint> m_RightKey;

// vector<DMatch> m_Matches;

// 以上三个变量已经被计算出来,分别是提取的关键点及其匹配,下面直接计算F

// 分配空间

int ptCount = (int)m_Matches.size();
Mat p1(ptCount, 2, CV_32F);
Mat p2(ptCount, 2, CV_32F);

// 把Keypoint转换为Mat

Point2f pt;
for (int i=0; i<ptCount; i++)
{
     pt = m_LeftKey[m_Matches[i].queryIdx].pt;
     p1.at<float>(i, 0) = pt.x;
     p1.at<float>(i, 1) = pt.y;
  
     pt = m_RightKey[m_Matches[i].trainIdx].pt;
     p2.at<float>(i, 0) = pt.x;
     p2.at<float>(i, 1) = pt.y;
}

// 用RANSAC方法计算F

// Mat m_Fundamental;

// 上面这个变量是基本矩阵

// vector<uchar> m_RANSACStatus;

// 上面这个变量已经定义过,用于存储RANSAC后每个点的状态

m_Fundamental = findFundamentalMat(p1, p2, m_RANSACStatus, FM_RANSAC);

// 计算野点个数

int OutlinerCount = 0;
for (int i=0; i<ptCount; i++)
{
     if (m_RANSACStatus[i] == 0) // 状态为0表示野点
     {
          OutlinerCount++;
     }
}

// 计算内点

// vector<Point2f> m_LeftInlier;
// vector<Point2f> m_RightInlier;
// vector<DMatch> m_InlierMatches;

// 上面三个变量用于保存内点和匹配关系

int InlinerCount = ptCount - OutlinerCount;
m_InlierMatches.resize(InlinerCount);
m_LeftInlier.resize(InlinerCount);
m_RightInlier.resize(InlinerCount);
InlinerCount = 0;
for (int i=0; i<ptCount; i++)
{
     if (m_RANSACStatus[i] != 0)
     {
          m_LeftInlier[InlinerCount].x = p1.at<float>(i, 0);
          m_LeftInlier[InlinerCount].y = p1.at<float>(i, 1);
          m_RightInlier[InlinerCount].x = p2.at<float>(i, 0);
          m_RightInlier[InlinerCount].y = p2.at<float>(i, 1);
          m_InlierMatches[InlinerCount].queryIdx = InlinerCount;
          m_InlierMatches[InlinerCount].trainIdx = InlinerCount;
          InlinerCount++;
     }
}

// 把内点转换为drawMatches可以使用的格式

vector<KeyPoint> key1(InlinerCount);
vector<KeyPoint> key2(InlinerCount);
KeyPoint::convert(m_LeftInlier, key1);
KeyPoint::convert(m_RightInlier, key2);

// 显示计算F过后的内点匹配

// Mat m_matLeftImage;
// Mat m_matRightImage;

// 以上两个变量保存的是左右两幅图像

Mat OutImage;
drawMatches(m_matLeftImage, key1, m_matRightImage, key2, m_InlierMatches, OutImage);
cvNamedWindow( "Match features", 1);
cvShowImage("Match features", &(IplImage(OutImage)));
cvWaitKey( 0 );
cvDestroyWindow( "Match features" );

4. 针孔相机模型和变形

这一节里的函数都使用针孔摄像机模型,这就是说,一幅视图是通过透视变换将三维空间中的点投影到图像平面。投影公式如下:

s \cdot m' = A\cdot[R|t] \cdot M'或者s\cdot  \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = \begin{bmatrix} fx & 0 & cx \\ 0 & fy & cy \\ 0 & 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} r_{11} & r_{12} & r_{13} & t_{1} \\ r_{21} & r_{22} & r_{23} & t_{2} \\ r_{31} & r_{32} & r_{33} & t_{3} \end{bmatrix} \cdot \begin{bmatrix} X \\ Y \\ Z \\ 1 \end{bmatrix}

这里(X, Y, Z)是一个点的世界坐标,(u, v)是点投影在图像平面的坐标,以像素为单位。A被称作摄像机矩阵,或者内参数矩阵。(cx, cy)是基准点(通常在图像的中心),fx, fy是以像素为单位的焦距。所以如果因为某些因素对来自于摄像机的一幅图像升采样或者降采样,所有这些参数(fx, fy, cx和cy)都将被缩放(乘或者除)同样的尺度。内参数矩阵不依赖场景的视图,一旦计算出,可以被重复使用(只要焦距固定)。旋转-平移矩阵[R|t]被称作外参数矩阵,它用来描述相机相对于一个固定场景的运动,或者相反,物体围绕相机的的刚性运动。也就是[R|t]将点(X, Y, Z)的坐标变换到某个坐标系,这个坐标系相对于摄像机来说是固定不变的。上面的变换等价与下面的形式(z≠0):
\begin{bmatrix}x \\ y \\z \end{bmatrix} = R \cdot \begin{bmatrix}X \\ Y \\Z \end{bmatrix} + t

x' = x / z

y' = y / z

u=fx \cdot x' + cx

v=fy \cdot y' + cy

真正的镜头通常有一些形变,主要的变形为径向形变,也会有轻微的切向形变。所以上面的模型可以扩展为:

\begin{bmatrix}x \\ y \\z \end{bmatrix} = R \cdot \begin{bmatrix}X \\ Y \\Z \end{bmatrix} + t

x' = x / z

y' = y / z

x'' = x' \cdot (1 + k_1 \cdot r^2 + k_2 \cdot r^4) + 2 \cdot p_1 \cdot x'\cdot y' + p_2 \cdot (r^2+2x'^2)

y'' = y' \cdot (1 + k_1 \cdot r^2 + k_2 \cdot r^4) + p_1 \cdot (r^2+2 \cdot y'^2) + 2 \cdot p_2 \cdot x'\cdot y'

这里 r2 = x'2 + y'2

u = fx \cdot x'' + cx

v = fy \cdot y'' + cy

k1和k2是径向形变系数,p1和p1是切向形变系数。OpenCV中没有考虑高阶系数。形变系数跟拍摄的场景无关,因此它们是内参数,而且与拍摄图像的分辨率无关。

后面的函数使用上面提到的模型来做如下事情:

  • 给定内参数和外参数,投影三维点到图像平面。
  • 给定内参数、几个三维点坐标和其对应的图像坐标,来计算外参数。
  • 根据已知的定标模式,从几个角度(每个角度都有几个对应好的3D-2D点对)的照片来计算相机的外参数和内参数。

5.照相机定标

5.1 ProjectPoints2

投影三维点到图像平面

void cvProjectPoints2( const CvMat* object_points, const CvMat* rotation_vector,
                       const CvMat* translation_vector, const CvMat* intrinsic_matrix,
                       const CvMat* distortion_coeffs, CvMat* image_points,
                       CvMat* dpdrot=NULL, CvMat* dpdt=NULL, CvMat* dpdf=NULL,
                       CvMat* dpdc=NULL, CvMat* dpddist=NULL );

object_points

物体点的坐标,为3xN或者Nx3的矩阵,这儿N是视图中的所有所有点的数目。

rotation_vector

旋转向量,1x3或者3x1。

translation_vector

平移向量,1x3或者3x1。

intrinsic_matrix

摄像机内参数矩阵A:\begin{bmatrix}fx & 0 & cx\\ 0 & fy & cy\\ 0&0&1\end{bmatrix}

istortion_coeffs

形变参数向量,4x1或者1x4,为[k1,k2,p1,p2]。如果是NULL,所有形变系数都设为0。

image_points

输出数组,存储图像点坐标。大小为2xN或者Nx2,这儿N是视图中的所有点的数目。

dpdrot

可选参数,关于旋转向量部分的图像上点的导数,Nx3矩阵。

dpdt

可选参数,关于平移向量部分的图像上点的导数,Nx3矩阵。

dpdf

可选参数,关于fx和fy的图像上点的导数,Nx2矩阵。

dpdc

可选参数,关于cx和cy的图像上点的导数,Nx2矩阵。

dpddist

可选参数,关于形变系数的图像上点的导数,Nx4矩阵。

函数cvProjectPoints2通过给定的内参数和外参数计算三维点投影到二维图像平面上的坐标。另外,这个函数可以计算关于投影参数的图像点偏导数的雅可比矩阵。雅可比矩阵可以用在cvCalibrateCamera2和cvFindExtrinsicCameraParams2函数的全局优化中。这个函数也可以用来计算内参数和外参数的反投影误差。注意,将内参数和(或)外参数设置为特定值,这个函数可以用来计算外变换(或内变换)。

5.2 FindHomographyFindHomography

计算两个平面之间的透视变换

void cvFindHomography( const CvMat* src_points,
                       const CvMat* dst_points,
                       CvMat* homography );
src_points

原始平面的点坐标,大小为2xN,Nx2,3xN或者 Nx3矩阵(后两个表示齐次坐标),这儿N表示点的数目。

dst_points

目标平面的点坐标大小为2xN,Nx2,3xN或者 Nx3矩阵(后两个表示齐次坐标)。

homography

输出的3x3的homography矩阵。

函数cvFindHomography计算源平面和目标平面之间的透视变换H=\begin{bmatrix}h_{ij}\end{bmatrix}_{i,j}

s_i \begin{bmatrix}x'_i \\ y'_i \\ 1\end{bmatrix}  \approx  H  \begin{bmatrix}x_i \\ y_i \\ 1\end{bmatrix}

使得反投影错误最小:

\sum_i((x'_i-\frac{h_{11}x_i + h_{12}y_i + h_{13}}{h_{31}x_i + h_{32}y_i + h_{33}})^2+          (y'_i-\frac{h_{21}x_i + h_{22}y_i + h_{23}}{h_{31}x_i + h_{32}y_i + h_{33}})^2)

这个函数可以用来计算初始的内参数和外参数矩阵。由于Homography矩阵的尺度可变,所以它被规一化使得h33= 1

5.3 CalibrateCamera2

利用定标来计算摄像机的内参数和外参数

void cvCalibrateCamera2( const CvMat* object_points, const CvMat* image_points,
                         const CvMat* point_counts, CvSize image_size,
                         CvMat* intrinsic_matrix, CvMat* distortion_coeffs,
                         CvMat* rotation_vectors=NULL,
                         CvMat* translation_vectors=NULL,
                         int flags=0 );

object_points

定标点的世界坐标,为3xN或者Nx3的矩阵,这里N是所有视图中点的总数。

image_points

定标点的图像坐标,为2xN或者Nx2的矩阵,这里N是所有视图中点的总数。

point_counts

向量,指定不同视图里点的数目,1xM或者Mx1向量,M是视图数目。

image_size

图像大小,只用在初始化内参数时。

intrinsic_matrix
输出内参矩阵(A)\begin{bmatrix}fx & 0 & cx\\ 0 & fy & cy \\ 0&0&1\end{bmatrix}

,如果指定CV_CALIB_USE_INTRINSIC_GUESS和(或)CV_CALIB_FIX_ASPECT_RATION,fx、 fy、 cx和cy部分或者全部必须被初始化。

distortion_coeffs

输出大小为4x1或者1x4的向量,里面为形变参数[k1, k2, p1, p2]。

rotation_vectors

输出大小为3xM或者Mx3的矩阵,里面为旋转向量(旋转矩阵的紧凑表示方式,具体参考函数cvRodrigues2)

translation_vectors

输出大小为3xM或Mx3的矩阵,里面为平移向量。

flags

不同的标志,可以是0,或者下面值的组合:

CV_CALIB_USE_INTRINSIC_GUESS - 内参数矩阵包含fx,fy,cx和cy的初始值。否则,(cx, cy)被初始化到图像中心(这儿用到图像大小),焦距用最小平方差方式计算得到。注意,如果内部参数已知,没有必要使用这个函数,使用cvFindExtrinsicCameraParams2则可。
CV_CALIB_FIX_PRINCIPAL_POINT - 主点在全局优化过程中不变,一直在中心位置或者在其他指定的位置(当CV_CALIB_USE_INTRINSIC_GUESS设置的时候)。
CV_CALIB_FIX_ASPECT_RATIO - 优化过程中认为fx和fy中只有一个独立变量,保持比例fx/fy不变,fx/fy的值跟内参数矩阵初始化时的值一样。在这种情况下, (fx, fy)的实际初始值或者从输入内存矩阵中读取(当CV_CALIB_USE_INTRINSIC_GUESS被指定时),或者采用估计值(后者情况中fx和fy可能被设置为任意值,只有比值被使用)。
CV_CALIB_ZERO_TANGENT_DIST – 切向形变参数(p1, p2)被设置为0,其值在优化过程中保持为0。
函数cvCalibrateCamera2从每个视图中估计相机的内参数和外参数。3维物体上的点和它们对应的在每个视图的2维投影必须被指定。这些可以通过使用一个已知几何形状且具有容易检测的特征点的物体来实现。这样的一个物体被称作定标设备或者定标模式,OpenCV有内建的把棋盘当作定标设备方法(参考cvFindChessboardCorners)。目前,传入初始化的内参数(当CV_CALIB_USE_INTRINSIC_GUESS不被设置时)只支持平面定标设备(物体点的Z坐标必须为全0或者全1)。不过3维定标设备依然可以用在提供初始内参数矩阵情况。在内参数和外参数矩阵的初始值都计算出之后,它们会被优化用来减小反投影误差(图像上的实际坐标跟cvProjectPoints2计算出的图像坐标的差的平方和)。

5.4 FindExtrinsicCameraParams2

计算指定视图的摄像机外参数

void cvFindExtrinsicCameraParams2( const CvMat* object_points,
                                   const CvMat* image_points,
                                   const CvMat* intrinsic_matrix,
                                   const CvMat* distortion_coeffs,
                                   CvMat* rotation_vector,
                                   CvMat* translation_vector );

object_points

定标点的坐标,为3xN或者Nx3的矩阵,这里N是视图中的个数。

image_points

定标点在图像内的坐标,为2xN或者Nx2的矩阵,这里N是视图中的个数。

intrinsic_matrix

内参矩阵(A) \begin{bmatrix}fx & 0 & cx\\ 0 & fy & cy \\ 0&0&1\end{bmatrix}

distortion_coeffs

大小为4x1或者1x4的向量,里面为形变参数[k1,k2,p1,p2]。如果是NULL,所有的形变系数都为0。

rotation_vector

输出大小为3x1或者1x3的矩阵,里面为旋转向量(旋转矩阵的紧凑表示方式,具体参考函数cvRodrigues2)。

translation_vector

大小为3x1或1x3的矩阵,里面为平移向量。

函数cvFindExtrinsicCameraParams2使用已知的内参数和某个视图的外参数来估计相机的外参数。3维物体上的点坐标和相应的2维投影必须被指定。这个函数也可以用来最小化反投影误差。

5.5 Rodrigues2

进行旋转矩阵和旋转向量间的转换

int  cvRodrigues2( const CvMat* src, CvMat* dst, CvMat* jacobian=0 );

src

输入的旋转向量(3x1或者1x3)或者旋转矩阵(3x3)。

dst

输出的旋转矩阵(3x3)或者旋转向量(3x1或者1x3)

jacobian

可选的输出雅可比矩阵(3x9或者9x3),关于输入部分的输出数组的偏导数。

函数转换旋转向量到旋转矩阵,或者相反。旋转向量是旋转矩阵的紧凑表示形式。旋转向量的方向是旋转轴,向量的长度是围绕旋转轴的旋转角。旋转矩阵R,与其对应的旋转向量r,通过下面公式转换:
\theta \leftarrow norm(r)

r \leftarrow r/\theta

R = \cos(\theta)I + (1-\cos(\theta))rr^T + \sin(\theta) \begin{bmatrix}0&-r_z&r_y\\ r_z&0&-r_x\\ -r_y&r_x&0\end{bmatrix}

反变换也可以很容易的通过如下公式实现:

\sin(\theta) \begin{bmatrix}0&-r_z&r_y\\ r_z&0&-r_x\\ -r_y&r_x&0\end{bmatrix} = \frac{R-R^T}{2}

旋转向量是只有3个自由度的旋转矩阵一个方便的表示,这种表示方式被用在函数cvFindExtrinsicCameraParams2和cvCalibrateCamera2内部的全局最优化中。

5.6 Undistort2

校正图像因相机镜头引起的变形

void cvUndistort2( const CvArr* src, CvArr* dst,
                   const CvMat* intrinsic_matrix,
                   const CvMat* distortion_coeffs );

src

原始图像(已经变形的图像)。只能变换32fC1的图像。

dst

结果图像(已经校正的图像)。

intrinsic_matrix

相机内参数矩阵,格式为 \begin{bmatrix}fx & 0 & cx\\ 0 & fy & cy\\ 0&0&1\end{bmatrix}

distortion_coeffs

四个变形系数组成的向量,大小为4x1或者1x4,格式为[k1,k2,p1,p2]。

函数cvUndistort2对图像进行变换来抵消径向和切向镜头变形。相机参数和变形参数可以通过函数cvCalibrateCamera2取得。使用本节开始时提到的公式,对每个输出图像像素计算其在输入图像中的位置,然后输出图像的像素值通过双线性插值来计算。如果图像得分辨率跟定标时用得图像分辨率不一样,fx、fy、cx和cy需要相应调整,因为形变并没有变化。

5.7 InitUndistortMap

计算形变和非形变图像的对应(map)

void cvInitUndistortMap( const CvMat* intrinsic_matrix,
                         const CvMat* distortion_coeffs,
                         CvArr* mapx, CvArr* mapy );

intrinsic_matrix

摄像机内参数矩阵(A) [fx 0 cx; 0 fy cy; 0 0 1].

distortion_coeffs

形变系数向量[k1, k2, p1, p2],大小为4x1或者1x4。

mapx

x坐标的对应矩阵。

mapy

y坐标的对应矩阵。

函数cvInitUndistortMap预先计算非形变对应-正确图像的每个像素在形变图像里的坐标。这个对应可以传递给cvRemap函数(跟输入和输出图像一起)。

5. 8 FindChessboardCorners

寻找棋盘图的内角点位置

int cvFindChessboardCorners( const void* image, CvSize pattern_size,
                             CvPoint2D32f* corners, int* corner_count=NULL,
                             int flags=CV_CALIB_CB_ADAPTIVE_THRESH );

image

输入的棋盘图,必须是8位的灰度或者彩色图像。

pattern_size

棋盘图中每行和每列角点的个数。

corners

检测到的角点

corner_count

输出,角点的个数。如果不是NULL,函数将检测到的角点的个数存储于此变量。

flags

各种操作标志,可以是0或者下面值的组合:

CV_CALIB_CB_ADAPTIVE_THRESH - 使用自适应阈值(通过平均图像亮度计算得到)将图像转换为黑白图,而不是一个固定的阈值。
CV_CALIB_CB_NORMALIZE_IMAGE - 在利用固定阈值或者自适应的阈值进行二值化之前,先使用cvNormalizeHist来均衡化图像亮度。
CV_CALIB_CB_FILTER_QUADS - 使用其他的准则(如轮廓面积,周长,方形形状)来去除在轮廓检测阶段检测到的错误方块。
函数cvFindChessboardCorners试图确定输入图像是否是棋盘模式,并确定角点的位置。如果所有角点都被检测到且它们都被以一定顺序排布(一行一行地,每行从左到右),函数返回非零值,否则在函数不能发现所有角点或者记录它们地情况下,函数返回0。例如一个正常地棋盘图右8x8个方块和7x7个内角点,内角点是黑色方块相互联通地位置。这个函数检测到地坐标只是一个大约地值,如果要精确地确定它们的位置,可以使用函数cvFindCornerSubPix。

DrawChessBoardCorners
绘制检测到的棋盘角点

void cvDrawChessboardCorners( CvArr* image, CvSize pattern_size,
                              CvPoint2D32f* corners, int count,
                              int pattern_was_found );

image

结果图像,必须是8位彩色图像。

pattern_size

每行和每列地内角点数目。

corners

检测到地角点数组。

count

角点数目。

pattern_was_found

指示完整地棋盘被发现(≠0)还是没有发现(=0)。可以传输cvFindChessboardCorners函数的返回值。

当棋盘没有完全检测出时,函数cvDrawChessboardCorners以红色圆圈绘制检测到的棋盘角点;如果整个棋盘都检测到,则用直线连接所有的角点。

6. 姿态估计

6.1 CreatePOSITObject

初始化包含对象信息的结构

CvPOSITObject* cvCreatePOSITObject( CvPoint3D32f* points, int point_count );

points

指向三维对象模型的指针

point_count

对象的点数

函数 cvCreatePOSITObject 为对象结构分配内存并计算对象的逆矩阵。

预处理的对象数据存储在结构CvPOSITObject中,只能在OpenCV内部被调用,即用户不能直接读写数据结构。用户只可以创建这个结构并将指针传递给函数。

对象是在某坐标系内的一系列点的集合,函数 cvPOSIT计算从照相机坐标系中心到目标点points[0] 之间的向量。

一旦完成对给定对象的所有操作,必须使用函数cvReleasePOSITObject释放内存。

6.2 POSIT

执行POSIT算法

void cvPOSIT( CvPOSITObject* posit_object, CvPoint2D32f* image_points, 
              double focal_length,
              CvTermCriteria criteria, CvMatr32f rotation_matrix, 
              CvVect32f translation_vector );

posit_object

指向对象结构的指针

image_points

指针,指向目标像素点在二维平面图上的投影。

focal_length

使用的摄像机的焦距

criteria

POSIT迭代算法程序终止的条件

rotation_matrix

旋转矩阵

translation_vector

平移矩阵.

函数 cvPOSIT 执行POSIT算法。图像坐标在摄像机坐标系统中给出。焦距可以通过摄像机标定得到。算法每一次迭代都会重新计算在估计位置的透视投影。

两次投影之间的范式差值是对应点中的最大距离。如果差值过小,参数criteria.epsilon就会终止程序。

6.3 ReleasePOSITObject

释放3D对象结构

void cvReleasePOSITObject( CvPOSITObject** posit_object );

posit_object

指向 CvPOSIT 结构指针的指针。

函数 cvReleasePOSITObject 释放函数 cvCreatePOSITObject分配的内存。

6.4 CalcImageHomography

计算长方形或椭圆形平面对象(例如胳膊)的Homography矩阵

void cvCalcImageHomography( float* line, CvPoint3D32f* center,
                            float* intrinsic, float* homography );

line

对象的主要轴方向,为向量(dx,dy,dz).

center

对象坐标中心 ((cx,cy,cz)).

intrinsic

摄像机内参数 (3x3 matrix).

homography

输出的Homography矩阵(3x3).

函数 cvCalcImageHomography 为从图像平面到图像平面的初始图像变化(defined by 3D oblong object line)计算Homography矩阵。


对极几何(双视几何)

6.5 FindFundamentalMat

由两幅图像中对应点计算出基本矩阵

int cvFindFundamentalMat( const CvMat* points1,
                          const CvMat* points2,
                          CvMat* fundamental_matrix,
                          int    method=CV_FM_RANSAC,
                          double param1=1.,
                          double param2=0.99,
                          CvMat* status=NULL);

points1

第一幅图像点的数组,大小为2xN/Nx2 或 3xN/Nx3 (N 点的个数),多通道的1xN或Nx1也可以。点坐标应该是浮点数(双精度或单精度)。:

points2

第二副图像的点的数组,格式、大小与第一幅图像相同。

fundamental_matrix

输出的基本矩阵。大小是 3x3 或者 9x3 ,(7-点法最多可返回三个矩阵).

method

计算基本矩阵的方法

CV_FM_7POINT – 7-点算法,点数目= 7
CV_FM_8POINT – 8-点算法,点数目 >= 8
CV_FM_RANSAC – RANSAC 算法,点数目 >= 8
CV_FM_LMEDS - LMedS 算法,点数目 >= 8
param1

这个参数只用于方法RANSAC 或 LMedS 。它是点到对极线的最大距离,超过这个值的点将被舍弃,不用于后面的计算。通常这个值的设定是0.5 or 1.0 。

param2

这个参数只用于方法RANSAC 或 LMedS 。 它表示矩阵正确的可信度。例如可以被设为0.99 。

status

具有N个元素的输出数组,在计算过程中没有被舍弃的点,元素被被置为1;否则置为0。这个数组只可以在方法RANSAC and LMedS 情况下使用;在其它方法的情况下,status一律被置为1。这个参数是可选参数。

对极几何可以用下面的等式描述:
p_2^T \cdot F \cdot p_1=0

其中 F 是基本矩阵,p1 和 p2 分别是两幅图上的对应点。

函数 FindFundamentalMat 利用上面列出的四种方法之一计算基本矩阵,并返回基本矩阵的值:没有找到矩阵,返回0,找到一个矩阵返回1,多个矩阵返回3。 计算出的基本矩阵可以传递给函数cvComputeCorrespondEpilines来计算指定点的对极线。

例子1:使用 RANSAC 算法估算基本矩阵。
int    numPoints = 100;
CvMat* points1;
CvMat* points2;
CvMat* status;
CvMat* fundMatr;
points1 = cvCreateMat(2,numPoints,CV_32F);
points2 = cvCreateMat(2,numPoints,CV_32F);
status  = cvCreateMat(1,numPoints,CV_32F);
 
/* 在这里装入对应点的数据... */
 
fundMatr = cvCreateMat(3,3,CV_32F);
int num = cvFindFundamentalMat(points1,points2,fundMatr,CV_FM_RANSAC,1.0,0.99,status);
if( num == 1 )
     printf("Fundamental matrix was found\n");
else
     printf("Fundamental matrix was not found\n");
 
 
例子2:7点算法(3个矩阵)的情况。
CvMat* points1;
CvMat* points2;
CvMat* fundMatr;
points1 = cvCreateMat(2,7,CV_32F);
points2 = cvCreateMat(2,7,CV_32F);
 
/* 在这里装入对应点的数据... */
 
fundMatr = cvCreateMat(9,3,CV_32F);
int num = cvFindFundamentalMat(points1,points2,fundMatr,CV_FM_7POINT,0,0,0);
printf("Found %d matrixes\n",num); 

6.6 ComputeCorrespondEpilines

为一幅图像中的点计算其在另一幅图像中对应的对极线

void cvComputeCorrespondEpilines( const CvMat* points,
                                  int which_image,
                                  const CvMat* fundamental_matrix,
                                  CvMat* correspondent_lines);

points

输入点,是2xN 或者 3xN 数组 (N为点的个数)

which_image

包含点的图像指数(1 or 2)

fundamental_matrix

基本矩阵

correspondent_lines

计算对极点, 3xN数组

函数 ComputeCorrespondEpilines 根据外级线几何的基本方程计算每个输入点的对应外级线。如果点位于第一幅图像(which_image=1),对应的对极线可以如下计算 :
l_2=F \cdot p_1

其中F是基本矩阵,p1 是第一幅图像中的点, l2 - 是与第二幅对应的对极线。如果点位于第二副图像中 which_image=2),计算如下:

l_1=F^T \cdot p_2

其中p2 是第二幅图像中的点,l1 是对应于第一幅图像的对极线,每条对极线都可以用三个系数表示 a, b, c:

a\cdot x + b\cdot y + c = 0

归一化后的对极线系数存储在correspondent_lines 中。 

Freak特征描述+BruteForceMatcher匹配+RANSAC剔除误匹配

时间一晃飞快,前几日的广州之行让人难忘,置身炎热的户外仿佛蒸桑拿一样。这一段时间一直在研究一种从一幅图片中找到给定目标的方法,而基于局部特征提取并且进行特征的匹配可以实现。 
目前局部特征描述子有以下几种,

2004年发展成熟的SIFT
后来改进的SURF,它们的描述符是高维度整型的数字。
2010年出现的BRIEF是一种二进制的描述符,这种描述符在提高计算速度方面给出了不错的思路,但是不具备旋转不变,尺度不变,对噪声也比较敏感(虽然它采用了高斯平滑来抗噪声)
2011年出现的ORB算法同样是二进制,但是解决了旋转不变和噪声敏感的问题,但是不具备尺度不变性。
2011年同时出现了BRISK算法,解决旋转不变,尺度不变和抗噪声问题,网上有视频,实时性也不错,先有个直观认识,爽死。链接[(http://v.youku.com/v_show/id_XMTI5MzI3Mzk0OA==.html)]
2012年提出FREAK算法,其实是在BRISK上的改进。具备BRISK的优点。这篇论文中作者表示FREAK算法可以outperform以前所有的描述子。(注意:这里是特征描述子,因为论文中没有给出特征点检测方法,只有特征点的描述子) 
看完了综述,我决定用FREAK来实现。足足耗了我15天,今天我这一段所得写下,勉励自己,服务他人。 
先贴上代码

/************************************************************************
* Copyright(c) 2016  唯疯
* All rights reserved.
*
* Brief: FAST特征点提取以及FREAK描述子的图像匹配,基于OpenCV2.4.8
* Version: 1.0
* Author: 唯疯
* Date: 2016/07/21
* Address: 广州&北京
************************************************************************/
#include <opencv2/core/core.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/nonfree/features2d.hpp>
#include <opencv2/legacy/legacy.hpp>
#include <iostream>
#include <vector>
 
using namespace std;
using namespace cv;
 
int main()
{
    Mat img1_src = imread("im5.jpg",0);
    Mat img2_src = imread("im6.jpg",0);
    //FastFeatureDetector fast(40);
    SurfFeatureDetector fast(2000,4);
    FREAK extractor;
    vector<KeyPoint> keypoints1,keypoints2;
    Mat descriptor1,descriptor2;
    vector<DMatch> final_matches; 
    vector<DMatch> matches;
 
    double t = (double)getTickCount();
    fast.detect(img1_src,keypoints1);
    fast.detect(img2_src,keypoints2);
    //drawKeypoints(img1_src,keypoints1,img1_src,Scalar(0,255,0));
    //drawKeypoints(img2_src,keypoints2,img2_src,Scalar(0,255,0));//问题在这里!!!醉了,这里的问题,浪费了我5天,欧耶,就是整整5天,由于我想在这里看一看检测出来的特征点是啥样的,就跟父母想第一眼就立刻看到刚出生的婴儿是一个心情的。结果,后来特征匹配就几乎没有什么正确匹配,害我还以为是FREAK这个描述子有问题呢。于是就各种找问题,看论文,翻墙逛论坛。最后一个不经意间发现了。问题就是:这俩个语句是把特征点又原封不动的画到了img1_src中,也就是原图像里面,而后来我进行特征点描述的时候,就直接在画满了特征点的图片下进行描述,而不是原图!不是原图啊!是充满了特征点的图片!所以后期再进行匹配的时候,显然,各种乱匹配,就跟隔壁家小狗似的,见了猫都想干坏事。于是乎,我直接注释掉了这两句!
    extractor.compute(img1_src,keypoints1,descriptor1);
    extractor.compute(img2_src,keypoints2,descriptor2);
    BFMatcher matcher(NORM_HAMMING,true);//暴力匹配,并且进行crosscheck,就是说第二个参数选择true。
    matcher.match(descriptor1,descriptor2,matches);
    final_matches=matches;
    cout<<"number of total_matches : "<<final_matches.size()<<endl;
 
//接下来是RANSAC剔除误匹配
    vector<Point2f> querymatches, trainmatches;
    vector<KeyPoint> p1,p2;
    for(int i=0;i<final_matches.size();i++)
    {
        p1.push_back(keypoints1[final_matches[i].queryIdx]);
        p2.push_back(keypoints2[final_matches[i].trainIdx]);
    }
    for(int i=0;i<p1.size();i++)
    {
        querymatches.push_back(p1[i].pt);
        trainmatches.push_back(p2[i].pt);
    }
    cout<<querymatches[1]<<" and "<<trainmatches[1]<<endl;
    vector<uchar> status;
    Mat h = findHomography(querymatches,trainmatches,status,CV_FM_RANSAC,10);
    int index=0;
    vector<DMatch> super_final_matches;
    for (int i=0;i<final_matches.size();i++)
    {
        cout<<status[i];
        if (status[i] != 0)
        {
        super_final_matches.push_back(final_matches[i]);
            index++;
        }
    }
    cout<<"number of inlier_matches : "<<index<<endl;
 
    Mat imgMatch;
    drawMatches(img1_src,keypoints1,img2_src,keypoints2,super_final_matches,imgMatch);
    imshow("imgMatch",imgMatch);
 
    t = ((double)getTickCount()-t)/getTickFrequency();
    cout<<" total time [s] : "<<t<<endl;
    waitKey(0);
    return 0;
}

看看简单的效果

一: 
首先说说特征点检测,为什么我用surf,而不用速度很快的fast或者Brisk,我自己做了测试,发现用fast和brisk的话,再用freak描述,最终匹配效果并不好,比sift要差…我其实不太明白为什么,因为有一些后期的论文对比效果说freak综合效果最好,然而我却得不到这样的结果。要继续钻研了,大家有想法欢迎交流,哦不,欢迎指点! 
二: 
RANSAC的这个算法,请用opencv库中的cvFindHomography,不要用cvFindFundamentalMat,因为这两个是不一样的,字面上简单看一个是找到单应性矩阵,一个是找到基础矩阵,这两个到底有神马区别??? 
《计算机视觉中的多视图几何》中这样描述:

2D单应:一幅2D图像中点集到另一幅图像中点集的射影变换。
基本矩阵:它使一个队所有的点都是的x`Fx=0成立的奇异矩阵F。这是一个3*3的奇异矩阵,而且秩为2,不可逆。它是2维图像上的点通过对极线束约束的映射,是2维到1维的映射。
具体的大家可以看书或者百度,或者等以后我有空了,再具体描述。其实我也没有太懂哈哈哈!还要再看看!慎言慎言! 
之前,我曾经自己编写过一个RANSAC的算法啊,但是效果不是特别理想,由于RANSAC算法本身也比较耗时,各路大神们也都提出了很多改进方法,不知道opencv跟进这些改进算法没有?推荐综述型文章《A Universal Framework for Random 
Sample Consensus》,里面提到的prosac,lo-ransac速度提升很多!大家自行汲取,我想以后再写一些关于RANSAC的博文。 
三: 
代码中用到BFMatcher这个类,对象的参数有两个BFMatcher::BFMatcher(int normType=NORM_L2, bool crossCheck=false )。第一个指的是距离,如果是sift或者surf描述的话,就选择NORM_L2和NORM_L2;ORB,BRISK,FREAK这种二进制的描述子就用汉明距离,也就是NORM_HAMMING;NORM_HAMMING2是用作ORB的。第二个参数是是否crossCheck,最好选择true,这个也经过我实验验证,选择true相当于knn匹配中将k置1,也就是只返回最好的那一组匹配,而不是多个近邻。(注意:如果要用knnmatch的话,这个crossCheck的参数就选择false)。这种只返回最好的一组匹配的方式可以作为ratio test的一种替代选择。实验证明,如果选择false,那么将会返回一对多匹配,和多对一的匹配,后期用ransac的效果可能非常不理想。crossCheck就是交叉匹配,剔除不好的匹配。有论文支撑,大家感兴趣可以去搜。

参考博文:
https://blog.csdn.net/u011867581/article/details/43818183
https://blog.csdn.net/luoshixian099/article/details/50217655

https://blog.csdn.net/lcb_coconut/article/details/52145293?locationNum=4