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角點檢測的代碼，如下，僅供參考：

#include <opencv2\opencv.hpp>
#include <iostream>
#include <string>

using namespace std;

#define HARRIS

cv::RNG rng(12345);

void calEigen2x2(cv::Mat cov,cv::Mat &eigenValue)
{
    int height = cov.rows;
    int width = cov.cols;

    float *pCov = (float*)cov.data;
    float *pEigenValue = (float*)eigenValue.data;
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            double a = pCov[(i*width + j) * 3 + 0];
            double b = pCov[(i*width + j) * 3 + 1];
            double c = pCov[(i*width + j) * 3 + 2];

            double tmp1 = (a + c) / 2.;
            double tmp2 = sqrtf(b*b + (a - c)*(a - c) / 4.);

            double alpha = tmp1 - tmp2;
            double beta = tmp1 + tmp2;

            pEigenValue[(i*width + j) * 2 + 0] =(float) alpha;
            pEigenValue[(i*width + j) * 2 + 1] =(float) beta;
        }
    }
}

void myCalEigenValues(cv::Mat srcImg, cv::Mat &eigenValue, int covWin, int sobelWin)
{
    //求微分
    cv::Mat sobelx, sobely;
    cv::Sobel(srcImg, sobelx, CV_32FC1, 1, 0, sobelWin, 1. / (255*12), 0, 4);
    cv::Sobel(srcImg, sobely, CV_32FC1, 0, 1, sobelWin, 1. / (255*12), 0, 4);

    cv::Mat cov = cv::Mat::zeros(srcImg.size(), CV_32FC3);
    int height = srcImg.rows;
    int width = srcImg.cols;
    float *pSobelX = (float*)sobelx.data;
    float *pSobelY = (float*)sobely.data;
    float *pCov = (float*)cov.data;
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            float dx = pSobelX[i*width + j];
            float dy = pSobelY[i*width + j];

            pCov[(i*width + j) * 3 + 0] = dx*dx;
            pCov[(i*width + j) * 3 + 1] = dx*dy;
            pCov[(i*width + j) * 3 + 2] = dy*dy;
        }
    }

    cv::boxFilter(cov, cov, cov.depth(), cv::Size(covWin, covWin), cv::Point(-1, -1), false, 4);

    calEigen2x2(cov, eigenValue);
}

void main()
{
    string imgPath = "data/srcImg/0.png";
    cv::Mat srcImg = cv::imread(imgPath, 1);
    cv::Mat grayImg;
    cv::cvtColor(srcImg, grayImg, CV_BGR2GRAY);
    int height = srcImg.rows;
    int width = srcImg.cols;

    cv::Mat eigenValue = cv::Mat::zeros(grayImg.size(), CV_32FC2);
    int covWin = 3, sobelWin = 3;
    myCalEigenValues(grayImg, eigenValue, covWin, sobelWin);

    cv::Mat Mc = cv::Mat::zeros(grayImg.size(), CV_32FC1);
#ifndef HARRIS
    //計算特征表達式值
    for (int i = 0; i<height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            float alpha = eigenValue.at<float>(i, j*2+0);
            float beta = eigenValue.at<float>(i, j*2+1);

            Mc.at<float>(i, j) = alpha*beta - 0.04f*pow((alpha + beta), 2);
        }
    }
#else //ShiTomasi
    for (int i = 0; i<height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            float alpha = eigenValue.at<float>(i, j * 2 + 0);
            float beta = eigenValue.at<float>(i, j * 2 + 1);

            float minEigenValue = (alpha > beta) ? (beta) : (alpha);

            Mc.at<float>(i, j) = minEigenValue;
        }
    }    
#endif

    double minVal, maxVal;
    cv::minMaxLoc(Mc, &minVal, &maxVal, 0, 0, cv::Mat());

    double thresh = (maxVal + minVal) / 2.;
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            double value = (double)Mc.at<float>(i, j);
            if (value > thresh)
            {
                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);
            }
        }
    }

    cv::namedWindow("show", 0);
    cv::imshow("show", srcImg);
    cv::waitKey(0);
}

 

---

 

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

 

---