OpenCV角点检测源代码分析（Harris和ShiTomasi角点）
OpenCV中常用的角点检测为Harris角点和ShiTomasi角点。
以OpenCV源代码文件 .\opencv\sources\samples\cpp\tutorial_code\TrackingMotion\cornerDetector_Demo.cpp为例，主要分析其中的这两种角点检测源代码。角点检测数学原理请参考我之前转载的一篇博客 http://www.cnblogs.com/riddick/p/7645904.html，分析的很详细，不再赘述。本文主要分析其源代码：
1. Harris角点检测
根据数学上的推导，可以根据图像中某一像素点邻域内构建的协方差矩阵获取特征值和特征向量，根据特征值建立特征表达式，如下：
(αβ) - k(α+β)^2
可以根据上式的值得大小来判断该像素点是平坦区域内点、边界点还是角点。下面说一下怎么在原图像中建立协方差矩阵并求取特征值α和β和特征向量t1, t2。
该例程代码中调用cornerEigenValsAndVecs（）函数计算特征值和特征向量。函数原型如下：
void cv::cornerEigenValsAndVecs( InputArray _src, OutputArray _dst, int blockSize, int ksize, int borderType )
src为输入灰度图像，dst为输出（6通道 CV_32FC(6)，依次保存的是α， t1, β， t2），blockSize为邻域大小，ksize为sobel求取微分时的窗口大小。
该函数内部调用cornerEigenValsVecs（）函数，原型如下：
static void cornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,int aperture_size, int op_type, double k=0.,int borderType=BORDER_DEFAULT )
主要介绍一下op_type这个参数，该参数是一个枚举值，有三个值可以选择（MINEIGENVAL， HARRIS， EIGENVALSVECS）
①MINEIGENVAL用于ShiTomasi角点检测中获取两个特征值中较小的那个值，用以获取强角点，随后介绍；
②HARRIS在cornerHarris（）函数中用到，用于直接利用协方差矩阵获取特征表达式值的大小，k值在此时会被设置，通常为0.04，其他情况下设置为0；
③EIGENVALSVECS就是本例程中设置的，求取两个特征值和特征向量。
在cornerEigenValsVecs（）函数中，先利用sobel算子求水平方向和竖直方向的微分，窗口大小为前述，如下代码：
Mat Dx, Dy; if( aperture_size > 0 ) { Sobel( src, Dx, CV_32F, 1, 0, aperture_size, scale, 0, borderType ); Sobel( src, Dy, CV_32F, 0, 1, aperture_size, scale, 0, borderType ); } else { Scharr( src, Dx, CV_32F, 1, 0, scale, 0, borderType ); Scharr( src, Dy, CV_32F, 0, 1, scale, 0, borderType ); }
然后初始化协方差矩阵cov（三通道，依次保存dx*dx, dx*dy, dy*dy），如下：
for( ; j < size.width; j++ ) { float dx = dxdata[j]; float dy = dydata[j]; cov_data[j*3] = dx*dx; cov_data[j*3+1] = dx*dy; cov_data[j*3+2] = dy*dy; }
接下来对协方差矩阵进行在前述设定窗口内进行均值（盒式）滤波：
boxFilter(cov, cov, cov.depth(), Size(block_size, block_size), Point(-1,-1), false, borderType ); if( op_type == MINEIGENVAL ) calcMinEigenVal( cov, eigenv ); else if( op_type == HARRIS ) calcHarris( cov, eigenv, k ); else if( op_type == EIGENVALSVECS ) calcEigenValsVecs( cov, eigenv );
然后就是利用滤波后的协方差矩阵求取特征值和特征向量了，根据设定不同的op_type调用不同的函数计算，本例程中为调用最后一个calcEigenValsVecs（）函数，该函数如下：
static void calcEigenValsVecs( const Mat& _cov, Mat& _dst ) { Size size = _cov.size(); if( _cov.isContinuous() && _dst.isContinuous() ) { size.width *= size.height; size.height = 1; } for( int i = 0; i < size.height; i++ ) { const float* cov = _cov.ptr<float>(i); float* dst = _dst.ptr<float>(i);
//调用该函数计算2x2协方差矩阵的特征值和特征向量 eigen2x2(cov, dst, size.width); } }
该函数中调用eigen2x2（）函数计算每个像素点处协方差矩阵的2个特征值和2个特征向量，协方差矩阵为如下形式，数据都保存在cov的三个通道中：
eigen2x2（）函数如下：2x2矩阵特征值和特征向量的计算，有线性代数基础的都学过，就不再赘述
static void eigen2x2( const float* cov, float* dst, int n ) { for( int j = 0; j < n; j++ ) { double a = cov[j*3]; double b = cov[j*3+1]; double c = cov[j*3+2]; double u = (a + c)*0.5; double v = std::sqrt((a - c)*(a - c)*0.25 + b*b);
//计算两个特征值l1,l2 double l1 = u + v; double l2 = u - v;
//计算特征值l1对应的特征向量 double x = b; double y = l1 - a; double e = fabs(x); if( e + fabs(y) < 1e-4 ) { y = b; x = l1 - c; e = fabs(x); if( e + fabs(y) < 1e-4 ) { e = 1./(e + fabs(y) + FLT_EPSILON); x *= e, y *= e; } } double d = 1./std::sqrt(x*x + y*y + DBL_EPSILON);
//保存特征值l1及其对应的特征向量 dst[6*j] = (float)l1; dst[6*j + 2] = (float)(x*d); dst[6*j + 3] = (float)(y*d);
//计算特征值l2对应的特征向量 x = b; y = l2 - a; e = fabs(x); if( e + fabs(y) < 1e-4 ) { y = b; x = l2 - c; e = fabs(x); if( e + fabs(y) < 1e-4 ) { e = 1./(e + fabs(y) + FLT_EPSILON); x *= e, y *= e; } } d = 1./std::sqrt(x*x + y*y + DBL_EPSILON);
//保存特征值l2及其对应的特征向量 dst[6*j + 1] = (float)l2; dst[6*j + 4] = (float)(x*d); dst[6*j + 5] = (float)(y*d); } }
求得2个特征值α、β和2个特征向量之后，就是要利用特征值构建特征表达式，通过表达式的值（ (αβ) - k(α+β)^2 ）来区分角点，k的值通常设置为0.04：
/* calculate Mc */ for( int j = 0; j < src_gray.rows; j++ ) { for( int i = 0; i < src_gray.cols; i++ ) { float lambda_1 = myHarris_dst.at<Vec6f>(j, i)[0]; float lambda_2 = myHarris_dst.at<Vec6f>(j, i)[1]; Mc.at<float>(j,i) = lambda_1*lambda_2 - 0.04f*pow( ( lambda_1 + lambda_2 ), 2 ); } }
代码中利用 minMaxLoc( Mc, &myHarris_minVal, &myHarris_maxVal, 0, 0, Mat() ); 函数获取特征表达式的最大值min和最小值max，通过选取不同的阈值min<=thresh<=max，来指定大于阈值thresh的表达式值对应的点为检测出的角点。并利用circle（）函数显示出来。
circle( myHarris_copy, Point(i,j), 4, Scalar( rng.uniform(0,255), rng.uniform(0,255), rng.uniform(0,255) ), -1, 8, 0 );
至此，Harris角点检测完成！
2. ShiTomasi角点检测
ShiTomasi角点提取是获取harris角点中的强角点，怎么获取强角点呢，那就是只选取两个特征值中较小的那个特征值构建特征表达式，如果较小的特征值都能够满足设定的阈值条件，那么该角点就视为强角点。
调用 cornerMinEigenVal( src_gray, myShiTomasi_dst, blockSize, apertureSize, BORDER_DEFAULT ); 函数来获取较小的特征值，其实该函数内部依然调用上面所述的函数 cornerEigenValsVecs( src, dst, blockSize, ksize, MINEIGENVAL, 0, borderType ); ，然后将op_type设置为MINEIGENVAL枚举值，进而调用 static void calcMinEigenVal( const Mat& _cov, Mat& _dst ) 函数计算较小的特征值。该函数代码如下：
static void calcMinEigenVal( const Mat& _cov, Mat& _dst ) { int i, j; Size size = _cov.size(); #if CV_TRY_AVX bool haveAvx = CV_CPU_HAS_SUPPORT_AVX; #endif #if CV_SIMD128 bool haveSimd = hasSIMD128(); #endif if( _cov.isContinuous() && _dst.isContinuous() ) { size.width *= size.height; size.height = 1; } for( i = 0; i < size.height; i++ ) { const float* cov = _cov.ptr<float>(i); float* dst = _dst.ptr<float>(i); #if CV_TRY_AVX if( haveAvx ) j = calcMinEigenValLine_AVX(cov, dst, size.width); else #endif // CV_TRY_AVX j = 0; #if CV_SIMD128 if( haveSimd ) { v_float32x4 half = v_setall_f32(0.5f); for( ; j <= size.width - v_float32x4::nlanes; j += v_float32x4::nlanes ) { v_float32x4 v_a, v_b, v_c, v_t; v_load_deinterleave(cov + j*3, v_a, v_b, v_c); v_a *= half; v_c *= half; v_t = v_a - v_c; v_t = v_muladd(v_b, v_b, (v_t * v_t)); v_store(dst + j, (v_a + v_c) - v_sqrt(v_t)); } } #endif // CV_SIMD128 for( ; j < size.width; j++ ) { float a = cov[j*3]*0.5f; float b = cov[j*3+1]; float c = cov[j*3+2]*0.5f;
//求根公式计算较小的根，即为较小的特征值 dst[j] = (float)((a + c) - std::sqrt((a - c)*(a - c) + b*b)); } } }
所有像素点处较小的特征值求出后，直接将该特征值作为特征表达式的值。利用 minMaxLoc( Mc, &myHarris_minVal, &myHarris_maxVal, 0, 0, Mat() ); 函数选取最小的min和最大的max，通过调整阈值thresh来设定大于阈值thresh的为显示出来的强角点。
至此，ShiTomasi角点检测完成！
#自己写了一个简单的Harris和ShiTomasi角点检测的代码，如下，仅供参考：
1 #include <opencv2\opencv.hpp> 2 #include <iostream> 3 #include <string> 4 5 using namespace std; 6 7 #define HARRIS 8 9 cv::RNG rng(12345); 10 11 void calEigen2x2(cv::Mat cov,cv::Mat &eigenValue) 12 { 13 int height = cov.rows; 14 int width = cov.cols; 15 16 float *pCov = (float*)cov.data; 17 float *pEigenValue = (float*)eigenValue.data; 18 for (int i = 0; i < height; i++) 19 { 20 for (int j = 0; j < width; j++) 21 { 22 double a = pCov[(i*width + j) * 3 + 0]; 23 double b = pCov[(i*width + j) * 3 + 1]; 24 double c = pCov[(i*width + j) * 3 + 2]; 25 26 double tmp1 = (a + c) / 2.; 27 double tmp2 = sqrtf(b*b + (a - c)*(a - c) / 4.); 28 29 double alpha = tmp1 - tmp2; 30 double beta = tmp1 + tmp2; 31 32 pEigenValue[(i*width + j) * 2 + 0] =(float) alpha; 33 pEigenValue[(i*width + j) * 2 + 1] =(float) beta; 34 } 35 } 36 } 37 38 void myCalEigenValues(cv::Mat srcImg, cv::Mat &eigenValue, int covWin, int sobelWin) 39 { 40 //求微分 41 cv::Mat sobelx, sobely; 42 cv::Sobel(srcImg, sobelx, CV_32FC1, 1, 0, sobelWin, 1. / (255*12), 0, 4); 43 cv::Sobel(srcImg, sobely, CV_32FC1, 0, 1, sobelWin, 1. / (255*12), 0, 4); 44 45 cv::Mat cov = cv::Mat::zeros(srcImg.size(), CV_32FC3); 46 int height = srcImg.rows; 47 int width = srcImg.cols; 48 float *pSobelX = (float*)sobelx.data; 49 float *pSobelY = (float*)sobely.data; 50 float *pCov = (float*)cov.data; 51 for (int i = 0; i < height; i++) 52 { 53 for (int j = 0; j < width; j++) 54 { 55 float dx = pSobelX[i*width + j]; 56 float dy = pSobelY[i*width + j]; 57 58 pCov[(i*width + j) * 3 + 0] = dx*dx; 59 pCov[(i*width + j) * 3 + 1] = dx*dy; 60 pCov[(i*width + j) * 3 + 2] = dy*dy; 61 } 62 } 63 64 cv::boxFilter(cov, cov, cov.depth(), cv::Size(covWin, covWin), cv::Point(-1, -1), false, 4); 65 66 calEigen2x2(cov, eigenValue); 67 } 68 69 void main() 70 { 71 string imgPath = "data/srcImg/0.png"; 72 cv::Mat srcImg = cv::imread(imgPath, 1); 73 cv::Mat grayImg; 74 cv::cvtColor(srcImg, grayImg, CV_BGR2GRAY); 75 int height = srcImg.rows; 76 int width = srcImg.cols; 77 78 cv::Mat eigenValue = cv::Mat::zeros(grayImg.size(), CV_32FC2); 79 int covWin = 3, sobelWin = 3; 80 myCalEigenValues(grayImg, eigenValue, covWin, sobelWin); 81 82 cv::Mat Mc = cv::Mat::zeros(grayImg.size(), CV_32FC1); 83 #ifndef HARRIS 84 //计算特征表达式值 85 for (int i = 0; i<height; i++) 86 { 87 for (int j = 0; j < width; j++) 88 { 89 float alpha = eigenValue.at<float>(i, j*2+0); 90 float beta = eigenValue.at<float>(i, j*2+1); 91 92 Mc.at<float>(i, j) = alpha*beta - 0.04f*pow((alpha + beta), 2); 93 } 94 } 95 #else //ShiTomasi 96 for (int i = 0; i<height; i++) 97 { 98 for (int j = 0; j < width; j++) 99 { 100 float alpha = eigenValue.at<float>(i, j * 2 + 0); 101 float beta = eigenValue.at<float>(i, j * 2 + 1); 102 103 float minEigenValue = (alpha > beta) ? (beta) : (alpha); 104 105 Mc.at<float>(i, j) = minEigenValue; 106 } 107 } 108 #endif 109 110 double minVal, maxVal; 111 cv::minMaxLoc(Mc, &minVal, &maxVal, 0, 0, cv::Mat()); 112 113 double thresh = (maxVal + minVal) / 2.; 114 for (int i = 0; i < height; i++) 115 { 116 for (int j = 0; j < width; j++) 117 { 118 double value = (double)Mc.at<float>(i, j); 119 if (value > thresh) 120 { 121 cv::circle(srcImg, cv::Point(j, i), 4, cv::Scalar(rng.uniform(0, 255), rng.uniform(0, 255), rng.uniform(0, 255)), -1, 8, 0); 122 } 123 } 124 } 125 126 cv::namedWindow("show", 0); 127 cv::imshow("show", srcImg); 128 cv::waitKey(0); 129 }
3. cornerHarris（）函数详解
前面讲述cornerEigenValsVecs（）这个函数是提到op_type这个枚举类型，有三个枚举值可以设置。其中MINEIGENVAL 和 EIGENVALSVECS都在前面介绍过。而 HARRIS则在cornerHarris（）函数中使用，这是一个公开的OpenCV API函数，函数原型如下：
void cv::cornerHarris( InputArray _src, OutputArray _dst, int blockSize, int ksize, double k, int borderType )
k值为上述特征表达式中的常数项。该函数内部其实还是调用cornerEigenValsVecs（）函数，只不过调用时将op_type设置为枚举值HARRIS。意思就是提取HARRIS角点，然后调用内部的 static void calcHarris( const Mat& _cov, Mat& _dst, double k ) 函数。只不过还内部函数不再计算特征值和特征向量，而是直接计算特征表达式的值。而特征表达式用下式表示：
其中矩阵M就是前面说的协方差矩阵，det(M)为M的行列式，Tr(M)为M的迹。在程序中代码如下：
for( ; j < size.width; j++ ) { float a = cov[j*3]; float b = cov[j*3+1]; float c = cov[j*3+2];
//求特征表达式的值 dst[j] = (float)(a*c - b*b - k*(a + c)*(a + c)); }
通常算出特征表达式的值后将其归一化到（0-255），然后可以直接设置阈值Thresh，特征表达式的值>Thresh对应的点视为角点。具体可以参见OpenCV例程源代码：.\opencv\sources\samples\cpp\tutorial_code\TrackingMotion\cornerHarris_Demo.cpp
