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个特征向量， 该函数如下：2×2矩阵特征值和特征向量的计算，有线性代数基础的都学过，就不再赘述
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角点检测完成！