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個特征向量， 該函數如下：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角點檢測完成！
OpenCV角點檢測源代碼分析（Harris和ShiTomasi角點）