OpenCV goodFeaturesToTrack特征提取解析笔记

1. goodFeaturesToTrack算法描述

$\text{goodFeaturesToTrack}$算法来自于$\text{Shi et al.}$的一篇论文，论文名称就叫$\text{Good Features to track}$。算法流程大致分为以下几个步骤：

• 生成特征矩阵：计算各个像素点用于判断$\text{Harris}$角点或$\text{Shi-Tomasi}$角点的相关特征值。
• 对特征矩阵进行分割：保留不小于特征矩阵中最大值某个百分比的所有特征值，并将其余特征值置零。
• 对特征矩阵进行非极大值抑制：如果某个特征值的$8$领域内存在有比其更大的特征值，那么该特征值将会被抑制掉。
• 对特征矩阵进行掩码抑制：根据掩码矩阵，对被掩中的特征值进行抑制。
• 剔除距离相近的特征值：对特征值进行排序，并由大到小进行筛选。如果待筛选特征值所对应的位置距离已筛选特征值所对应的位置相近，那么该特征值将会被抑制掉。

2. goodFeaturesToTrack代码解析

在$\text{OpenCV}$中，关于$\text{goodFeaturesToTrack}$函数的核心代码位于$\text{featureselect.cpp}$文件中：

void cv::goodFeaturesToTrack( InputArray _image, OutputArray _corners,
                              int maxCorners, double qualityLevel, double minDistance,
                              InputArray _mask, int blockSize, int gradientSize,
                              bool useHarrisDetector, double harrisK )
                                                                // _image：输入图像
                                                                // _corners：输出角点
                                                                // _maxCorners：最大角点数
                                                                // qualityLevel：角点质量等级（百分比）
                                                                // minDistance：角点最小间距
                                                                // _mask：像素掩码矩阵
                                                                // blockSize = 3：角点观察窗大小
                                                                // gradientSize = 3：梯度算子阶数（ Sobel 算子阶数）
                                                                // useHarrisDetector = false：是否检测 Harris 角点
                                                                // harrisK = 0.04： Harris 角点阈值
{
    CV_INSTRUMENT_REGION();

    CV_Assert( qualityLevel > 0 && minDistance >= 0 && maxCorners >= 0 );
                                                                // 参数检查
    CV_Assert( _mask.empty() || (_mask.type() == CV_8UC1 && _mask.sameSize(_image)) );
                                                                // 参数检查

    CV_OCL_RUN(_image.dims() <= 2 && _image.isUMat(),
               ocl_goodFeaturesToTrack(_image, _corners, maxCorners, qualityLevel, minDistance,
                                    _mask, blockSize, gradientSize, useHarrisDetector, harrisK))
                                                                // OpenCL 加速

    Mat image = _image.getMat(), eig, tmp;
    if (image.empty())
    {
        _corners.release();
        return;
    }

    // Disabled due to bad accuracy
    CV_OVX_RUN(false && useHarrisDetector && _mask.empty() &&
               !ovx::skipSmallImages<VX_KERNEL_HARRIS_CORNERS>(image.cols, image.rows),
               openvx_harris(image, _corners, maxCorners, qualityLevel, minDistance, blockSize, gradientSize, harrisK))
                                                                // OpenVX 加速

    if( useHarrisDetector )
        cornerHarris( image, eig, blockSize, gradientSize, harrisK );
    else
        cornerMinEigenVal( image, eig, blockSize, gradientSize );
                                                                // 计算各像素点的最小特征值，生成特征值矩阵

    double maxVal = 0;
    minMaxLoc( eig, 0, &maxVal, 0, 0, _mask );                  // 根据 _mask 矩阵，查找特征值矩阵中的最大值
    threshold( eig, eig, maxVal*qualityLevel, 0, THRESH_TOZERO );
                                                                // 对特征值矩阵进行阈值分割，将小于 maxVal * qualityLevel 的特征值清零。
    dilate( eig, tmp, Mat());                                   // 对特征值矩阵进行膨胀操作，使用 3 x 3 结构元素，进行非极大值抑制

    Size imgsize = image.size();
    std::vector<const float*> tmpCorners;

    // collect list of pointers to features - put them into temporary image
    Mat mask = _mask.getMat();
    for( int y = 1; y < imgsize.height - 1; y++ )
    {
        const float* eig_data = (const float*)eig.ptr(y);
        const float* tmp_data = (const float*)tmp.ptr(y);
        const uchar* mask_data = mask.data ? mask.ptr(y) : 0;

        for( int x = 1; x < imgsize.width - 1; x++ )
        {
            float val = eig_data[x];
            if( val != 0 && val == tmp_data[x] && (!mask_data || mask_data[x]) )
                                                                // 将没有被 mask 掩中的极大值特征值地址放入 tmpCorners 向量
                tmpCorners.push_back(eig_data + x);
        }
    }

    std::vector<Point2f> corners;
    size_t i, j, total = tmpCorners.size(), ncorners = 0;

    if (total == 0)
    {
        _corners.release();
        return;
    }

    std::sort( tmpCorners.begin(), tmpCorners.end(), greaterThanPtr() );
                                                                // 对 tmpCorners 进行排序

    if (minDistance >= 1)                                       // 如果对 tmpCorners 有 minDistance 约束
    {
         // Partition the image into larger grids
        int w = image.cols;
        int h = image.rows;

        const int cell_size = cvRound(minDistance);             // 对 minDistance 进行四舍五入，确定网格大小
        const int grid_width = (w + cell_size - 1) / cell_size; // 判断图像宽度可以容纳多少网格
        const int grid_height = (h + cell_size - 1) / cell_size;// 判断图像高度可以容纳多少网格

        std::vector<std::vector<Point2f> > grid(grid_width*grid_height);
                                                                // 根据 grid_width 以及 grid_height 生成网格

        minDistance *= minDistance;                             // 对 minDistance 进行平方

        /* 对 tmpCorners 进行遍历  */
        for( i = 0; i < total; i++ )
        {
            int ofs = (int)((const uchar*)tmpCorners[i] - eig.ptr());
                                                                // tmpCorners 在特征值矩阵中的偏移量
            int y = (int)(ofs / eig.step);                      // tmpCorners 对应像素点的x坐标
            int x = (int)((ofs - y*eig.step)/sizeof(float));    // tmpCorners 对应像素点的y坐标

            bool good = true;                                   // good 标志置为 true

            int x_cell = x / cell_size;                         // 判断 tmpCorners 落在哪个网格中（ x 方向编号）
            int y_cell = y / cell_size;                         // 判断 tmpCorners 落在哪个网格中（ y 方向编号）

            int x1 = x_cell - 1;                                // 确定 tmpCorners 网格左上方网格 x 编号
            int y1 = y_cell - 1;                                // 确定 tmpCorners 网格左上方网格 y 编号
            int x2 = x_cell + 1;                                // 确定 tmpCorners 网格右下方网格 x 编号
            int y2 = y_cell + 1;                                // 确定 tmpCorners 网格右下方网格 y 编号

            // boundary check
            x1 = std::max(0, x1);                               // 对网格编号进行修正，防止溢出
            y1 = std::max(0, y1);                               // 对网格编号进行修正，防止溢出
            x2 = std::min(grid_width-1, x2);                    // 对网格编号进行修正，防止溢出
            y2 = std::min(grid_height-1, y2);                   // 对网格编号进行修正，防止溢出

            /* 对九个网格进行遍历 */
            for( int yy = y1; yy <= y2; yy++ )
            {
                for( int xx = x1; xx <= x2; xx++ )
                {
                    std::vector <Point2f> &m = grid[yy*grid_width + xx];

                    /* 如果当前网格存在特征点，则遍历所有特征点，并判断 tmpCorners 与该特征点的距离是否满足 minDistance 约束，如果不满足，则将 good 标志置为 false 并跳出网格遍历 */
                    if( m.size() )
                    {
                        for(j = 0; j < m.size(); j++)
                        {
                            float dx = x - m[j].x;
                            float dy = y - m[j].y;

                            if( dx*dx + dy*dy < minDistance )
                            {
                                good = false;
                                goto break_out;
                            }
                        }
                    }
                }
            }

            break_out:

            if (good)                                           // 如果 tmpCorners 满足 minDistance 约束
            {
                grid[y_cell*grid_width + x_cell].push_back(Point2f((float)x, (float)y));
                                                                // 将 tmpCorners 的坐标记录到网格中

                corners.push_back(Point2f((float)x, (float)y)); // 将 tmpCorners 记录到 corners 向量
                ++ncorners;                                     // ncorners 计数器++

                if( maxCorners > 0 && (int)ncorners == maxCorners )
                                                                // 如果当记录的角点数已经达到 maxCorners ，则跳出 tmpCorners 遍历
                    break;
            }
        }
    }
    else
    {
        for( i = 0; i < total; i++ )
        {
            int ofs = (int)((const uchar*)tmpCorners[i] - eig.ptr());
            int y = (int)(ofs / eig.step);
            int x = (int)((ofs - y*eig.step)/sizeof(float));

            corners.push_back(Point2f((float)x, (float)y));
            ++ncorners;
            if( maxCorners > 0 && (int)ncorners == maxCorners )
                break;
        }
    }

    Mat(corners).convertTo(_corners, _corners.fixedType() ? _corners.type() : CV_32F);
                                                                // 将角点拷贝至 _corners
}