＜图像处理＞ Shi-Tomasi角点检测

Shi-Tomasi角点检测

参考Harris角点检测方法，它是通过下式来判断是否为角点的：
$R=det(C(x,y))-k(trace(C(x,y)))^2$
式中，C(x,y)为自相关矩阵。Shi和Tomasi对此提出一种改进方法，通过两个特征值 $\lambda_1$ 和 $\lambda_2$ 的最小值 $\lambda_{min}$ 来判断是否为角点，因为角点检测的不确定性完全取决于这个最小值。 $\lambda_{min}$ 很大，则说明是角点； $\lambda_{min}$ 很小，说明要么是边缘（最大值很大），要么是平坦区（两个值都很小）。
该方法称为Shi-Tomasi角点检测方法，而OpenCV之所以把方法称为，是源于提出该方法的那篇论文题目；

检测步骤

计算图像中的每个像素的自相关矩阵的两个特征值中的最小特征值，但有时为了加快运行速度，最小特征值也用上式 $R$ 代替；
设定一个阈值，抛弃那些小于阈值的特征值所对应的角点；
在 $3 * 3$ 的邻域内对特征值进行非极大值抑制；
对剩下的特征值进行由大到小的顺序排序；
设定一个最小距离，计算特征值与其排序靠前的特征值的两点之间的坐标距离，抛弃那些小于最小距离的特征值，也就是在一定范围内只保留那些特征值最小的角点。

OpenCV函数

函数位于sources\modules\imgproc\src\featureselect.cpp

void cv::goodFeaturesToTrack(InputArray image, OutputArray corners, int maxCorners, double qualityLevel, double minDistance, InputArray mask = noArray(), int blockSize = 3, bool useHarrisDetector = false, double k = 0.04)	

Parameters
image				输入图像，单通道8位整型或32位浮点数；
corners				输出检测到角点的坐标；
maxCorners			最大输出角点数；如果角比找到的更多，则返回其中最强的角。maxCorners <= 0表示没有设置最大值的限制，并返回所有检测到的角点。
qualityLevel		决定步骤2的阈值大小；描述图像角点最小可接受质量的参数。参数值乘以最佳角点质量度量，即最小特征值(见cornerMinEigenVal)或Harris函数响应(见cornerHarris)。质量指标低于产品的边角被拒收。例如，如果最好的角点的质量度量值为1500，且qualityLevel=0.01，那么所有质量度量值小于15的角点都被拒绝。
minDistance			决定步骤5中的最小距离（欧氏距离）；
mask				掩膜区域，决定着图像中那些区域需要计算角点；
blockSize			Harris检测中步骤3中的加权窗口尺寸，可参见cornerEigenValsAndVecs。
useHarrisDetector	是使用cornerHarris函数还是使用cornerMinEigenVal函数，cornerMinEigenVal函数可以直接得到最小的特征值，cornerHarris函数使用上式中的R代替最小特征值，默认false使用cornerMinEigenVal函数；
k					上式中Harris计算中的k值；

下面重载的方法：
新增参数gradientSize，参数意指Sobel算子用于导数计算的孔径参数，参见cornerEigenValsAndVecs。
void cv::goodFeaturesToTrack(InputArray image, OutputArray corners, int maxCorners, double qualityLevel, double minDistance, InputArray mask, int blockSize, int gradientSize, bool useHarrisDetector = false, double k = 0.04)

OpenCV源码分析

void cv::goodFeaturesToTrack( InputArray _image, OutputArray _corners,
                              int maxCorners, double qualityLevel, double minDistance,
                              InputArray _mask, int blockSize, int gradientSize,
                              bool useHarrisDetector, double harrisK )
{
    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))

    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))

    //步骤1，eig为输出特征值矩阵，存储着最小特征值或Harris中的R值；
    if( useHarrisDetector )
        cornerHarris( image, eig, blockSize, gradientSize, harrisK ); //计算每个像素R值
    else
        cornerMinEigenVal( image, eig, blockSize, gradientSize ); //计算每个像素最小特征值

    double maxVal = 0;
    minMaxLoc( eig, 0, &maxVal, 0, 0, _mask ); //计算图像中所有角点的最大值
    threshold( eig, eig, maxVal*qualityLevel, 0, THRESH_TOZERO ); //步骤3操作，阈值为maxVal*qualityLevel
    dilate( eig, tmp, Mat()); //对角点的特征值图像矩阵进行灰度形态学的膨胀操作，结构元素3x3，膨胀运算的目的是在3x3邻域内选择最大值

    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]; //当前像素最小特征值
            //判断3x3邻域内，进行非极大值抑制，val == tmp_data[x]为判断是否为3x3邻域内最大值
            if( val != 0 && val == tmp_data[x] && (!mask_data || mask_data[x]) )
                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;
    }

    //步骤4，排序
    std::sort( tmpCorners.begin(), tmpCorners.end(), greaterThanPtr() );

    //步骤5，距离筛选
    if (minDistance >= 1)
    {
         // Partition the image into larger grids
        //把图像分割成大小相同的块，在块的8邻域范围内（即9个块）计算特征值之间的L2范数
        int w = image.cols;
        int h = image.rows;

        const int cell_size = cvRound(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);

        //距离平方，节省开根号的操作
        minDistance *= minDistance;

        //从大到小遍历所有角点
        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)); //横坐标

            bool good = true; //标记，true为距离大于最小距离，false为舍弃该角点

            //得到该角点所在的块位置
            int x_cell = x / cell_size;
            int y_cell = y / cell_size;

            //得到该角点所在块的8邻域范围，（x1,y1）左上角，（x2,y2）右下角
            int x1 = x_cell - 1;
            int y1 = y_cell - 1;
            int x2 = x_cell + 1;
            int y2 = y_cell + 1;

            // 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);

            //遍历8邻域内的所有块（即3x3=9块）
            for( int yy = y1; yy <= y2; yy++ )
            {
                for( int xx = x1; xx <= x2; xx++ )
                {
                    //得到块中所有角点的坐标，这些角点是在以前遍历角点时存储到向量grid中，所以向量grid内的角点的特征值一定大于当前角点的特征值，只需考虑距离即可
                    std::vector <Point2f> &m = grid[yy*grid_width + xx];

                    if( m.size() ) //块内有角点
                    {
                        for(j = 0; j < m.size(); j++)
                        {
                            //计算当前角点与块内角点的距离
                            float dx = x - m[j].x;
                            float dy = y - m[j].y;

                            //只要有一个距离小于minDistance 就抛弃该角点
                            if( dx*dx + dy*dy < minDistance )
                            {
                                good = false;
                                goto break_out;
                            }
                        }
                    }
                }
            }

            break_out:

            if (good)
            {
                //角点坐标存入向量grid中
                grid[y_cell*grid_width + x_cell].push_back(Point2f((float)x, (float)y));

                //再次把该角点坐标存入向量corners中
                corners.push_back(Point2f((float)x, (float)y));
                ++ncorners;

                //当前角点数等于maxcorners时，停止检测角点
                if( maxCorners > 0 && (int)ncorners == maxCorners )
                    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);
}