Opencv2.4.9源码分析——SURF



SURF (Speeded Up Robust Features)是一种具有鲁棒性的局部特征检测算法，它首先由Herbert Bay等人于2006年提出，并在2008年进行了完善。其实该算法是Herbert Bay在博士期间的研究内容，并作为博士毕业论文的一部分发表。

SURF算法的部分灵感来自于SIFT算法，但正如它的名字一样，该算法除了具有重复性高的检测器和可区分性好的描述符特点外，还具有很强的鲁棒性以及更高的运算速度，如Bay所述，SURF至少比SIFT快3倍以上，综合性能要优于SIFT算法。与SIFT算法一样，SURF算法也在美国申请了专利。

与我的上一篇SIFT文章一样，该文分为原理解释、源码分析和应用实例。但该文的内容也较多，对公式和图表的粘贴复制太麻烦，所以我也把这篇文章上传到了百度文库和csdn，这样更方便大家阅读，网址为：

http://download.csdn.net/detail/zhaocj/8347849

http://wenku.baidu.com/view/2f1e4d8ef705cc1754270945.html

在这里我仅把源码分析部分复制到这里，建议大家还是把原理和源码结合着看，这样更容易理解，所以最好还是看详细的全文。

实现SURF算法的文件为sources/modules/nonfree/src/surf.cpp

SURF类的构造函数为：

SURF::SURF()

SURF::SURF(doublehessianThreshold, int nOctaves=4, int nOctaveLayers=2, bool extended=true,

bool upright=false )

参数的含义为：

hessianThreshold为Hessian矩阵行列式响应值的阈值

nOctaves为图像堆的组数

nOctaveLayers为图像堆中每组中的中间层数，该值加2等于每组图像中所包含的层数

extended表示是128维描述符，还是64维描述符，为true时，表示128维描述符

upright表示是否采用U-SURF算法，为true时，采用U-SURF算法

SURF类的重载运算符( )为：

void SURF::operator()(InputArray _img, InputArray _mask,
                      CV_OUT vector<KeyPoint>& keypoints,
                      OutputArray _descriptors,
                      bool useProvidedKeypoints) const
//_img表示输入8位灰度图像
//_mask为掩码矩阵
//keypoints为特征点矢量数组
//_descriptors为特征点描述符
//useProvidedKeypoints表示是否运行特征点的检测，该值为true表示不进行特征点的检测，而只是利用输入的特征点进行特征点描述符的运算。
{
    //定义、初始化各种矩阵，sum为输入图像的积分图像矩阵，mask1和msum为掩码所需矩阵
    Mat img = _img.getMat(), mask = _mask.getMat(), mask1, sum, msum;
    //变量doDescriptors表示是否进行特征点描述符的运算
    bool doDescriptors = _descriptors.needed();
    //确保输入图像的正确
    CV_Assert(!img.empty() && img.depth() == CV_8U);
    if( img.channels() > 1 )
        cvtColor(img, img, COLOR_BGR2GRAY);
    //确保各类输入参数的正确
    CV_Assert(mask.empty() || (mask.type() == CV_8U && mask.size() == img.size()));
    CV_Assert(hessianThreshold >= 0);
    CV_Assert(nOctaves > 0);
    CV_Assert(nOctaveLayers > 0);
    //计算输入图像的积分图像
    integral(img, sum, CV_32S);

    // Compute keypoints only if we are not asked for evaluating the descriptors are some given locations:
    //是否进行特征点检测运算，useProvidedKeypoints=false表示进行该运算
    if( !useProvidedKeypoints )
    {
        if( !mask.empty() )    //掩码
        {
            cv::min(mask, 1, mask1);
            integral(mask1, msum, CV_32S);    //msum为掩码图像的积分图像
        }
        //进行Fast-Hessian特征点检测，fastHessianDetector函数在后面给出详细解释
        fastHessianDetector( sum, msum, keypoints, nOctaves, nOctaveLayers, (float)hessianThreshold );
    }
    //N为所检测到的特征点数量
    int i, j, N = (int)keypoints.size();
    if( N > 0 )    //如果检测到了特征点
    {
        Mat descriptors;    //描述符变量矩阵
        bool _1d = false;
        //由extended变量得到描述符是128维的还是64维的
        int dcols = extended ? 128 : 64;
        //一个描述符的存储空间大小
        size_t dsize = dcols*sizeof(float);

        if( doDescriptors )    //需要得到描述符
        {
            //定义所有描述符的排列形式
            _1d = _descriptors.kind() == _InputArray::STD_VECTOR && _descriptors.type() == CV_32F;
            if( _1d )
            {
                //所有特征点的描述符连在一起，组成一个矢量
                _descriptors.create(N*dcols, 1, CV_32F);
                descriptors = _descriptors.getMat().reshape(1, N);
            }
            else
            {
                //一个特征点就是一个描述符，一共N个描述符
                _descriptors.create(N, dcols, CV_32F);
                descriptors = _descriptors.getMat();
            }
        }

        // we call SURFInvoker in any case, even if we do not need descriptors,
        // since it computes orientation of each feature.
        //方向角度的分配和描述符的生成，SURFInvoker在后面给出了详细的解释
        parallel_for_(Range(0, N), SURFInvoker(img, sum, keypoints, descriptors, extended, upright) );

        // remove keypoints that were marked for deletion
        //删除那些被标注为无效的特征点
        // i表示特征点删除之前的特征点索引，j表示删除之后的特征点索引
        for( i = j = 0; i < N; i++ )
        {
            if( keypoints[i].size > 0 )    //size大于0，表示该特征点为有效的特征点
            {
                if( i > j )    // i > j表示在i之前有被删除的特征点
                {
                    //用后面的特征点依次填补被删除的那些特征点的位置
                    keypoints[j] = keypoints[i];    
                    //相应的描述符的空间位置也要调整
                    if( doDescriptors )
                        memcpy( descriptors.ptr(j), descriptors.ptr(i), dsize);
                }
                j++;    //索引值计数
            }
        }
        //N表示特征点删除之前的总数，j表示删除之后的总数，N > j表明有一些特征点被删除      
        if( N > j )
        {
            //重新赋值特征点的总数
            N = j; 
            keypoints.resize(N);
            if( doDescriptors )
            {
                //根据描述符的排列形式，重新赋值
                Mat d = descriptors.rowRange(0, N);
                if( _1d )
                    d = d.reshape(1, N*dcols);
                d.copyTo(_descriptors);
            }
        }
    }
}

SURF算法中特征点检测函数fastHessianDetector：

static void fastHessianDetector( const Mat& sum, const Mat& mask_sum, vector<KeyPoint>& keypoints, int nOctaves, int nOctaveLayers, float hessianThreshold )
{
    /* Sampling step along image x and y axes at first octave. This is doubled
       for each additional octave. WARNING: Increasing this improves speed,
       however keypoint extraction becomes unreliable. */
    //为了减少运算量，提高计算速度，程序设置了采样间隔。即在第一组的各层图像的每个像素都通过盒状滤波器模块进行滤波处理，而在第二组的各层图像中采用隔点（采样间隔为21）滤波器处理的方式，第三组的采样间隔为22，其他组以此类推。
    //定义初始采样间隔，即第一组各层图像的采样间隔
    const int SAMPLE_STEP0 = 1;
    //总的层数，它等于组数乘以每组的层数
    int nTotalLayers = (nOctaveLayers+2)*nOctaves;
    //中间层数，即进行非最大值抑制的层数
    int nMiddleLayers = nOctaveLayers*nOctaves;
    //各类矢量变量，dets为Hessian矩阵行列式矢量矩阵，traces为Hessian矩阵迹矢量矩阵，sizes为盒状滤波器的尺寸大小矢量，sampleSteps为每一层的采样间隔矢量，middleIndices为中间层（即进行非最大值抑制的层）的索引矢量
    vector<Mat> dets(nTotalLayers);
    vector<Mat> traces(nTotalLayers);
    vector<int> sizes(nTotalLayers);
    vector<int> sampleSteps(nTotalLayers);
    vector<int> middleIndices(nMiddleLayers);
    //清空特征点变量
    keypoints.clear();

    // Allocate space and calculate properties of each layer
    //index为变量nTotalLayers的索引，middleIndex为变量nMiddleLayers的索引，step为采样间隔
    int index = 0, middleIndex = 0, step = SAMPLE_STEP0;
    //分配内存空间，计算各层的数据
    //由于采用了采样间隔，因此各层图像所检测到的行列式、迹等数量会不同，这样就需要为每一层分配不同的空间大小，保存不同数量的数据
    for( int octave = 0; octave < nOctaves; octave++ )
    {
        for( int layer = 0; layer < nOctaveLayers&#