＜图像处理＞ 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之所以把方法称为，是源于提出该方法的那篇论文题目；

检测步骤

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

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)	
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

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(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 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 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 > 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  &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);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173

参考

1. 《图像局部特征检测与描述》