網上有很多Harris角點檢測的資料，但是或多或少有一些缺陷，有的講清楚了原理，卻沒有給出程式碼的實現，有的給出了程式碼的實現，卻沒給出原理性的解釋。為此，我找出了幾篇部落格，將他們綜合整理在一起，有一套理論到程式碼的完整實現過程，這才是學一門演算法或運算元的正確方法。

本文主要參考部落格：

原理部分正文：

Harris角點檢測運算元是於1988年由CHris Harris & Mike Stephens提出來的。在具體展開之前，不得不提一下Moravec早在1981就提出來的Moravec角點檢測運算元。

1.Moravec角點檢測運算元

Moravec角點檢測的思想其實特別簡單，在影象上取一個W*W的滑動視窗，不斷的移動這個視窗並檢測視窗的畫素變化情況E。

A. 視窗影象平坦 ---------------E的變化不大

B.視窗影象是一條邊 ------------1.沿邊滑動E的變化不大      2.垂直於邊滑動，E的變化會很大

C.視窗影象為一個角點 ------------視窗影象沿任何方向移動，E的值變化都會很大

上圖很形象的描述了上述過程。

數學表示式也很簡單：

其中(x,y)就表示四個移動方向（1,0）（0,1）（1,1）（-1,1）

E --------------------畫素的變化值

Moravec運算元對四個移動方向加權求和來確定變化的大小

然後設定閾值,來確定到底是邊還是角點

2.Harris角點檢測運算元（Moravec 運算元的改進版）

在原文中，作者提出了三點Moravec運算元的缺陷並提出了改良方法

1.Moravec運算元方向依賴性太強（只移動了4個45度角的離散方向）

其中：           

可以推出：

2.由於Moravec運算元採用的是正方形Windows, 因此E的響應比較容易受到干擾（？？？？這個地方有點不明白）

                    Harris採用了一個較為平滑的視窗--------------高斯函式：

3.Moravec運算元對邊緣響應過於靈敏。為此，Harris提出了對E進行變形：

變成了二次型（大家對這個概念肯定不陌生，但就是想不起來，不要緊，我看到這裡也是這樣，後面會有詳細的講解），其中，

α，β表示矩陣M的特徵值,這樣會產生三種情況：

A.α，β都很小,說明高斯Windows中的影象接近平坦

B.如果一個大一個小,則表示檢測到邊

C.α，β都很大,那麼表示檢測到角點。

矩陣理論不說了，個人覺得非常難學，更不用說把那些概念牢記於心，一般都是用到了再來翻閱，為了防止今後再花時間去翻閱

幾個概念特地copy如下，博主的地址已經在樓上給出，相信看完下面的補充大家的不解之謎都會得到解決。

1. 關於特徵值和特徵向量：

特徵值的特徵向量的概念忘了就自己查吧，這裡只說關鍵的。對於實對稱矩陣M（設階數為n），則一定有n個實特徵值，每個特徵值對應一組特徵向量（這組向量中所有向量共線），不同特徵值對應的特徵向量間相互正交；（注意這裡說的是實對稱矩陣，不是所有的矩陣都滿足這些條件）


2. 關於對角化：
對角化是指存在一個正交矩陣Q，使得  Q’MQ 能成為一個對角陣（只有對角元素非0），其中Q’是Q的轉置（同時也是Q的逆，因為正交矩陣的轉置就是其逆）。一個矩陣對角化後得到新矩陣的行列式和矩陣的跡（對角元素之和）均與原矩陣相同。如果M是n階實對稱矩陣，則Q中的第 j 列就是第 j 個特徵值對應的一個特徵向量（不同列的特徵向量兩兩正交）。


3. 關於二次型：
對於一個n元二次多項式，f(x1,x2....xn) = ∑ ( aij*xi*xj ) ，其中 i 和 j 的求和區間均為 [1,n] ，可將其各次的係數 aij 寫成一個n*n矩陣M，由於 aij 和 aji 的對稱等價關係，一般將 aij 和 aji 設為一樣的值，均為 xi*xj 的係數的二分之一。這樣，矩陣M就是實對稱矩陣了。即二次型的矩陣預設都是實對稱矩陣


4. 關於二次型的標準化（正交變換法）：
二次型的標準化是指通過構造一個n階可逆矩陣 C，使得向量 ( x1,x2...xn ) = C * (y1,y2...yn)，把n維向量 x 變換成n維向量 y ，並代入f(x1,x2....xn) 後得到 g(y1,y2...yn)，而後者的表示式中的二次項中不包含任何交叉二次項 yi*yj（全部都是平方項 yi^2），也即表示式g的二次型矩陣N是對角陣。用公式表示一下 f 和 g ，（下面的表示式中 x 和 y都代表向量，x' 和 y' 代表轉置）


f = x' * M * x ;


g = f = x' * M * x = (Cy)' * M * (Cy) = y' * (C'MC) * y = y' * N * y  ;

因此 C‘MC = N。正交變換法，就是直接將M對角化得到N，而N中對角線的元素就是M的特徵值。正交變換法中得到的 C 正好是一個正交矩陣，其每一列都是兩兩正交的單位向量，因此 C 的作用僅僅是將座標軸旋轉（不會有放縮）。

OK，基礎知識補充完了，再來說說Harris角點檢測中的特徵值是怎麼回事。這裡的 M 是

將M對角化後得到矩陣N，他們都是2階矩陣，且N的對角線元素就是本文中提到的 α 和 β。

本來 E(x,y) = A*x^2 + 2*C*x*y + B*y^2 ，而將其標準後得到新的座標 xp和yp，這時表示式中就不再含有交叉二次項，新表示式如下：

E(x,y) = Ep (xp,yp) = α*xp^2 + β*yp^2，

我們不妨畫出 Ep (xp,yp) = 1 的等高線L ，即

α*xp^2 + β*yp^2 = 1 ，

可見這正好是(xp,yp)空間的一個橢圓，而α 和 β則分別是該橢圓長、短軸平方的倒數（或者反過來），且長短軸的方向也正好是α 和 β對應的特徵向量的方向。由於(x,y)空間只是 (xp,yp)空間的旋轉，沒有放縮，因此等高線L在(x,y)空間也是一個全等的橢圓，只不過可能是傾斜的。

現在就能理解下面的圖片中出現的幾個橢圓是怎麼回事了，圖(a)中畫的正是高度為 1 的等高線，（其他”高度“處的等高線也是橢圓，只不過長短軸的長度還要乘以一個係數）。其他的幾幅圖片中可以看到，“平坦”區域由於（高度）變化很慢，等高線（橢圓）就比較大；而”邊緣“區域則是在一個軸向上高度變化很快，另一個與之垂直的軸向上高度變化很慢，因此一個軸很長一個軸很短；“角點”區域各個方向高度都變化劇烈，因此橢圓很小。我們人眼可以直觀地看到橢圓的大小胖瘦，但如何讓計算機識別這三種不同的幾何特徵呢？為了能區分出角點、邊緣和平坦區域我們現在需要用α 和 β構造一個特徵表示式，使得這個特徵式在三種不同的區域有明顯不同的值。一個表現還不錯的特徵表示式就是：

(αβ) - k(α+β)^2

表示式中的 k 的值怎麼選取呢？它一般是一個遠小於 1 的係數，OpenCV的預設推薦值是 0.04(=0.2的平方)，它近似地表達了一個閾值：當橢圓短、長軸的平方之比（亦即α 和 β兩個特徵值之比）小於這個閾值時，認為該橢圓屬於“一個軸很長一個軸很短”，即對應的點會被認為是邊緣區域。

對於邊緣部分，（假設較大的特徵值為β）由於 β>>α且α<kβ，因此特徵式 ：

(αβ) - k(α+β)^2 ≈ αβ - kβ^2 < (kβ)β - kβ^2 = 0

即邊緣部分的特徵值小於0 ；

對於非邊緣部分，α 和 β相差不大，可認為 (α+β)^2 ≈ 4αβ，因此特徵式：


(αβ) - k(α+β)^2 ≈ αβ - 4kαβ = ( 1 - 4k ) * αβ

由於 k 遠小於1，因此 1 - 4k ≈ 1，這樣特徵式進一步近似為：

(αβ) - k(α+β)^2 ≈ αβ

在角點區域，由於α 和 β都較大，對應的特徵式的值也就很大；而在平坦區域，特徵式的值則很小。


因此，三種不同區域的判別依據就是： 如果特徵表示式的值為負，則屬於邊緣區域；如果特徵表示式的值較大，則屬於角點區域；如果特徵表示式的值很小，則是平坦區域。

最後，由於αβ和(α+β)正好是M對角化後行列式和跡，再結合上面補充的基礎知識第2條中提到的行列式和跡在對角化前後不變，就可以得到 (αβ) - k(α+β)^2 = det(M) - k*Tr(M)^2，這就是Harris檢測的表示式。

上述橢圓的來源上述註解中說的很清楚（見紅色部分）

接著有人又要問了，你怎麼知道我檢測到α，β算大還是算小？對此天才哈里斯定義了一個角點響應函式：

我們驚喜的發現，R為正值時，檢測到的是角點，R為負時檢測到的是邊，R很小時檢測到的是平坦區域。至於他怎麼想出來的，

天才的思維誰知道呢?哈哈哈哈

Harris角點檢測演算法由諸多優點：A旋轉不變性，橢圓轉過一定角度但是其形狀保持不變（特徵值保持不變）

B對於影象灰度的仿射變化具有部分的不變性,由於僅僅使用了影象的一階導數,對於影象灰度平移變化不變；對於影象灰度尺度變化不變

當然Harris也有許多不完善的地方： A它對尺度很明感，不具備幾何尺度不變性。

B提取的角點是畫素級的。

好的，理論部分就到此為止，如有任何問題，歡迎留言交流，我們再來看看OpenCV Harris角點的實現。

OpenCV中cornerHarris函式可用於檢測影象的Harris角點。 cornerHarris函式名引數的說明：

void cornerHarris( InputArray src,  //輸入8bit 單通道灰度Mat矩陣
                    OutputArray dst, //用於儲存Harris角點檢測結果,32位單通道,大小與src相同
                    int blockSize,   //滑塊視窗的尺寸
                    int ksize,       //Sobel邊緣檢測濾波器大小
                    double k,        //Harries中間引數,經驗值0.04~0.06
                    int borderType=BORDER_DEFAULT   //插值型別
                    );

Samples:

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <iostream>

using namespace cv;
using namespace std;


Mat image;
Mat imageGray;
int thresh=200;
int MaxThresh=255;

void Trackbar(int,void*);  //閾值控制

int main()  
{  
	image=imread("Test.bmp");
	cvtColor(image,imageGray,CV_RGB2GRAY);
	GaussianBlur(imageGray,imageGray,Size(5,5),1); // 濾波
	namedWindow("Corner Detected");
	createTrackbar("threshold：","Corner Detected",&thresh,MaxThresh,Trackbar);
	imshow("Corner Detected",image);
	Trackbar(0,0);
	waitKey();
	return 0;
}  

void Trackbar(int,void*)
{
	Mat dst,dst8u,dstshow,imageSource;
	dst=Mat::zeros(image.size(),CV_32FC1);  
	imageSource=image.clone();
	cornerHarris(imageGray,dst,3,3,0.04,BORDER_DEFAULT);
	normalize(dst,dst8u,0,255,CV_MINMAX);  //歸一化
	convertScaleAbs(dst8u,dstshow);
	imshow("dst",dstshow);  //dst顯示
	for(int i=0;i<image.rows;i++)
	{
		for(int j=0;j<image.cols;j++)
		{
			if(dstshow.at<uchar>(i,j)>thresh)  //閾值判斷
			{
				circle(imageSource,Point(j,i),2,Scalar(0,0,255),2); //標註角點
			}
		}
	}
	imshow("Corner Detected",imageSource);
}

兩個可能不懂的函式註解：

1.convertScaleAbs函式是OpenCV中的函式，使用線性變換轉換輸入陣列元素成8位無符號整型。

執行效果如下：

原始碼分析環節

//--------------------【cornerHarris原始碼分析】------------------------------------
//Harries角點實現函式,擷取cornerHarries中的關鍵程式碼做了簡化
void myHarries( const Mat& src, Mat& eigenv, int block_size, int aperture_size, double k=0.)
{
    eigenv.create(src.size(), CV_32F);
    Mat Dx, Dy;
    //soble operation get Ix, Iy
    Sobel(src, Dx, CV_32F, 1, 0, aperture_size);
    Sobel(src, Dy, CV_32F, 0, 1, aperture_size);

    //get covariance matrix
    Size size = src.size();
    Mat cov(size, CV_32FC3);     //建立一個三通道cov矩陣分別儲存[Ix*Ix, Ix*Iy, Iy*Iy];

    for( int i=0; i < size.height;i++)
    {
        float* cov_data = cov.ptr<float>(i);
        const float* dxdata = Dx.ptr<float>(i);
        const float* dydata = Dy.ptr<float>(i);

        for( int j=0; j < size.width; j++ )
        {
            float dx = dxdata[j];
            float dy = dydata[j];

            cov_data[j*3] = dx * dx;    //即  Ix*Ix
            cov_data[j*3 + 1] = dx*dy;  //即  Ix*Iy
            cov_data[j*3 + 2] = dy*dy;  //即  Iy*Iy
        }
    }


    //方框濾波W(x,y)卷積, 也可用高斯核加權...
    //W(Y,Y)與矩陣cov卷積運算得到 H 矩陣,後面通過H矩陣的特徵值決定是否是角點
    boxFilter(cov, cov, cov.depth(), Size(block_size,block_size),Point(-1,-1),false);

    //cale Harris
    size = cov.size();
    if( cov.isContinuous() && eigenv.isContinuous())
    {
        size.width *= size.width;
        size.height = 1;
        //cout << "yes"<<endl;
    }

    //此處計算響應R= det(H) - k*trace(H)*trace(H);
    for (int i = 0; i < size.height; i++)
    {
        const float* covPtr = cov.ptr<float>(i);
        float* dstPtr = eigenv.ptr<float>(i);

        for( int j = 0; j < size.width; j++)
        {
            float a = covPtr[j*3];
            float b = covPtr[j*3 + 1];
            float c = covPtr[j*3 + 2];

            //根據公式 R = det(H) - k*trace(H)*trace(H);
            dstPtr[j] = (float)(a*c - b*b - k*(a + c)(a + c));
        }
    }

    double max, min;
    minMaxLoc(eigenv, &min, &max);
    //cout<< max << endl;
    double threshold = 0.1*max;
    cv::threshold(eigenv, eigenv,threshold, 1, CV_THRESH_BINARY);   //eigenv的型別是CV_32F,
    imshow("eigenv", eigenv);
}