＜图像处理＞ 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 accuracyCV_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值elsecornerMinEigenVal( image, eig, blockSize, gradientSize ); //计算每个像素最小特征值double maxVal = 0;minMaxLoc( eig, 0, &maxVal, 0, 0, _mask ); //计算图像中所有角点的最大值threshold( eig, eig, maxVal*qualityLevel, 0, THRESH_TOZERO ); //步骤3操作，阈值为maxVal*qualityLeveldilate( eig, tmp, Mat()); //对角点的特征值图像矩阵进行灰度形态学的膨胀操作，结构元素3x3，膨胀运算的目的是在3x3邻域内选择最大值Size imgsize = image.size();std::vector<const float*> tmpCorners;// collect list of pointers to features - put them into temporary imageMat 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 checkx1 = 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);
}