1. 不同類型的角點
在現實世界中，角點對應於物體的拐角，道路的十字路口、丁字路口等。從圖像分析的角度來定義角點可以有以下兩種定義：
- 角點可以是兩個邊緣的角點；
- 角點是鄰域內具有兩個主方向的特征點；
前者往往需要對圖像邊緣進行編碼，這在很大程度上依賴於圖像的分割與邊緣提取，具有相當大的難度和計算量，且一旦待檢測目標局部發生變化，很可能導致操作的失敗。早期主要有Rosenfeld和Freeman等人的方法，后期有CSS等方法。
基於圖像灰度的方法通過計算點的曲率及梯度來檢測角點，避免了第一類方法存在的缺陷，此類方法主要有Moravec算子、Forstner算子、Harris算子、SUSAN算子等。
這篇文章主要介紹的Harris角點檢測的算法原理，比較著名的角點檢測方法還有jianbo Shi和Carlo Tomasi提出的Shi-Tomasi算法，這個算法開始主要是為了解決跟蹤問題，用來衡量兩幅圖像的相似度，我們也可以把它看為Harris算法的改進。OpenCV中已經對它進行了實現，接口函數名為GoodFeaturesToTrack()。另外還有一個著名的角點檢測算子即SUSAN算子，SUSAN是Smallest Univalue Segment Assimilating Nucleus（最小核值相似區）的縮寫。SUSAN使用一個圓形模板和一個圓的中心點，通過圓中心點像素與模板圓內其他像素值的比較，統計出與圓中心像素近似的像元數量，當這樣的像元數量小於某一個閾值時，就被認為是要檢測的角點。我覺得可以把SUSAN算子看為Harris算法的一個簡化。這個算法原理非常簡單，算法效率也高，所以在OpenCV中，它的接口函數名稱為：FAST() 。
2. Harris角點
2.1 基本原理
人眼對角點的識別通常是在一個局部的小區域或小窗口完成的。如果在各個方向上移動這個特征的小窗口，窗口內區域的灰度發生了較大的變化，那么就認為在窗口內遇到了角點。如果這個特定的窗口在圖像各個方向上移動時，窗口內圖像的灰度沒有發生變化，那么窗口內就不存在角點；如果窗口在某一個方向移動時，窗口內圖像的灰度發生了較大的變化，而在另一些方向上沒有發生變化，那么，窗口內的圖像可能就是一條直線的線段。
對於圖像$I(x,y)$，當在點$(x,y)$處平移$(\Delta x,\Delta y)$后的自相似性，可以通過自相關函數給出：
$$c(x,y;\Delta x,\Delta y) = \sum_{(u,v)\in W(x,y)}w(u,v)(I(u,v) – I(u+\Delta x,v+\Delta y))^2$$
其中，$W(x,y)$是以點$(x,y)$為中心的窗口，$w(u,v)$為加權函數，它既可是常數，也可以是高斯加權函數。
根據泰勒展開，對圖像$I(x,y)$在平移$(\Delta x,\Delta y)$后進行一階近似：
$$I(u+\Delta x,v+\Delta y) = I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y+O(\Delta x^2,\Delta y^2)\approx I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y$$
其中，$I_x,I_y$是圖像$I(x,y)$的偏導數，這樣的話，自相關函數則可以簡化為：
$$c(x,y;\Delta x,\Delta y)\approx \sum_w (I_x(u,v)\Delta x+I_y(u,v)\Delta y)^2=[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix}$$
其中
$$M(x,y)=\sum_w \begin{bmatrix}I_x(x,y)^2&I_x(x,y)I_y(x,y) \\ I_x(x,y)I_y(x,y)&I_y(x,y)^2\end{bmatrix} = \begin{bmatrix}\sum_w I_x(x,y)^2&\sum_w I_x(x,y)I_y(x,y) \\\sum_w I_x(x,y)I_y(x,y)&\sum_w I_y(x,y)^2\end{bmatrix}=\begin{bmatrix}A&C\\C&B\end{bmatrix} $$
也就是說圖像$I(x,y)$在點$(x,y)$處平移$(\Delta x,\Delta y)$后的自相關函數可以近似為二項函數：
$$c(x,y;\Delta x,\Delta y)\approx A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2$$
其中
$$A=\sum_w I_x^2, B=\sum_w I_y^2,C=\sum_w I_x I_y$$
二次項函數本質上就是一個橢圓函數。橢圓的扁率和尺寸是由$M(x,y)$的特征值$\lambda_1、\lambda_2$決定的，橢賀的方向是由$M(x,y)$的特征矢量決定的，如下圖所示，橢圓方程為：
$$[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix} = 1$$
橢圓函數特征值與圖像中的角點、直線（邊緣）和平面之間的關系如下圖所示。共可分為三種情況：
- 圖像中的直線。一個特征值大，另一個特征值小，$\lambda_1\gg \lambda_2$或$\lambda_2\gg \lambda_1$。自相關函數值在某一方向上大，在其他方向上小。
- 圖像中的平面。兩個特征值都小，且近似相等；自相關函數數值在各個方向上都小。
- 圖像中的角點。兩個特征值都大，且近似相等，自相關函數在所有方向都增大。
根據二次項函數特征值的計算公式，我們可以求$M(x,y)$矩陣的特征值。但是Harris給出的角點差別方法並不需要計算具體的特征值，而是計算一個角點響應值$R$來判斷角點。$R$的計算公式為：
$$R=det \boldsymbol{M} - \alpha(trace\boldsymbol{M})^2$$
式中，$det\boldsymbol{M}$為矩陣$\boldsymbol{M}=\begin{bmatrix}A&B\\B&C\end{bmatrix}$的行列式；$trace\boldsymbol{M}$為矩陣$\boldsymbol{M}$的直跡；$\alpha$為經常常數，取值范圍為0.04~0.06。事實上，特征是隱含在$det\boldsymbol{M}$和$trace\boldsymbol{M}$中，因為，
$$det\boldsymbol{M} = \lambda_1\lambda_2=AC-B^2$$
$$trace\boldsymbol{M}=\lambda_2+\lambda_2=A+C$$
2.2 Harris角點算法實現
根據上述討論，可以將Harris圖像角點檢測算法歸納如下，共分以下五步：
1. 計算圖像$I(x,y)$在$X$和$Y$兩個方向的梯度$I_x、I_y$。
$$I_x=\frac{\partial I}{\partial x}=I\otimes(-1\ 0\ 1)，I_y =\frac{\partial I}{\partial x}=I\otimes(-1\ 0\ 1)^T $$
2. 計算圖像兩個方向梯度的乘積。
$$I_x^2=I_x\cdot I_y，I_y^2=I_y\cdot I_y，I_{xy}=I_x\cdot I_y$$
3. 使用高斯函數對$I_x^2、I_y^2和I_{xy}$進行高斯加權（取$\sigma=1$），生成矩陣$M$的元素$A、B$和$C$。
$$A=g(I_x^2)=I_x^2\otimes w，C=g(I_y^2)=I_y^2\otimes w，B=g(I_{x,y})=I_{xy}\otimes w$$
4. 計算每個像素的Harris響應值$R$，並對小於某一閾值$t$的$R$置為零。
$$R=\{R:det \boldsymbol{M} - \alpha(trace\boldsymbol{M})^2<t\}$$
5. 在$3\times3$或$5\times5$的鄰域內進行非最大值抑制，局部最大值點即為圖像中的角點。
Harris角點檢測的C++實現代碼：https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harris.cpp
2.3 Harris角點的性質
1. 參數$\alpha$對角點檢測的影響
假設已經得到了矩陣$\boldsymbol{M}$的特征值$\lambda_1\ge\lambda_2\ge0$，令$\lambda_2=k\lambda_1,0\le k\le 1$。由特征值與矩陣$\boldsymbol{M}$的直跡和行列式的關系可得：
$$det\boldsymbol{M}=\prod_i\lambda_i \ \ \ \ \ \ trace\boldsymbol{M}=\sum_i\lambda_i$$
從而可以得到角點的響應
$$R=\lambda_2\lambda_2=\alpha(\lambda_2+\lambda_2)^2=\lambda^2(k-\alpha(1+k)^2)$$
假設$R\ge0$，則有：
$$0\le \alpha \le\frac{k}{(1+k)^2}\le0.25$$
對於較小的$k$值，$R\approx\lambda^2(k-\alpha),\alpha<k$。
由此，可以得出這樣的結論：增大$\alpha$的值，將減小角點響應值$R$，降低角點檢測的靈性，減少被檢測角點的數量；減小$\alpha$值，將增大角點響應值$R$，增加角點檢測的靈敏性，增加被檢測角點的數量。
2. Harris角點檢測算子對亮度和對比度的變化不敏感
這是因為在進行Harris角點檢測時，使用了微分算子對圖像進行微分運算，而微分運算對圖像密度的拉升或收縮和對亮度的抬高或下降不敏感。換言之，對亮度和對比度的仿射變換並不改變Harris響應的極值點出現的位置，但是，由於閾值的選擇，可能會影響角點檢測的數量。
3. Harris角點檢測算子具有旋轉不變性
Harris角點檢測算子使用的是角點附近的區域灰度二階矩矩陣。而二階矩矩陣可以表示成一個橢圓，橢圓的長短軸正是二階矩矩陣特征值平方根的倒數。當特征橢圓轉動時，特征值並不發生變化，所以判斷角點響應值$R$也不發生變化，由此說明Harris角點檢測算子具有旋轉不變性。
4. Harris角點檢測算子不具有尺度不變性
如下圖所示，當右圖被縮小時，在檢測窗口尺寸不變的前提下，在窗口內所包含圖像的內容是完全不同的。左側的圖像可能被檢測為邊緣或曲線，而右側的圖像則可能被檢測為一個角點。
2.4 Harris的OpenCV接口
OpenCV的Hairrs角點檢測的函數為cornerHairrs()，但是它的輸出是一幅浮點值圖像，浮點值越高，表明越可能是特征角點，我們需要對圖像進行閾值化。
C++: void cornerHarris(InputArray src, OutputArray dst, int blockSize, int apertureSize, double k, int borderType = BORDER_DEFAULT);
- src – 輸入的單通道8-bit或浮點圖像。
- dst – 存儲着Harris角點響應的圖像矩陣，大小與輸入圖像大小相同，是一個浮點型矩陣。
- blockSize – 鄰域大小。
- apertureSize – 擴展的微分算子大。
- k – 響應公式中的，參數$\alpha$。
- boderType – 邊界處理的類型。
int main() { Mat image = imread("../buliding.png"); Mat gray; cvtColor(image, gray, CV_BGR2GRAY); Mat cornerStrength; cornerHarris(gray, cornerStrength, 3, 3, 0.01); threshold(cornerStrength, cornerStrength, 0.001, 255, THRESH_BINARY); return 0; }
從上面上間一幅圖像我們可以看到，有很多角點都是粘連在一起的，我們下面通過加入非極大值抑制來進一步去除一些粘在一起的角點。
非極大值抑制原理是，在一個窗口內，如果有多個角點則用值最大的那個角點，其他的角點都刪除，窗口大小這里我們用3*3，程序中通過圖像的膨脹運算來達到檢測極大值的目的，因為默認參數的膨脹運算就是用窗口內的最大值替代當前的灰度值。
int main() { Mat image = imread("buliding.png"); Mat gray; cvtColor(image, gray, CV_BGR2GRAY); Mat cornerStrength; cornerHarris(gray, cornerStrength, 3, 3, 0.01); double maxStrength; double minStrength; // 找到圖像中的最大、最小值 minMaxLoc(cornerStrength, &minStrength, &maxStrength); Mat dilated; Mat locaMax; // 膨脹圖像，最找出圖像中全部的局部最大值點 dilate(cornerStrength, dilated, Mat()); // compare是一個邏輯比較函數，返回兩幅圖像中對應點相同的二值圖像 compare(cornerStrength, dilated, locaMax, CMP_EQ); Mat cornerMap; double qualityLevel = 0.01; double th = qualityLevel*maxStrength; // 閾值計算 threshold(cornerStrength, cornerMap, th, 255, THRESH_BINARY); cornerMap.convertTo(cornerMap, CV_8U); // 逐點的位運算 bitwise_and(cornerMap, locaMax, cornerMap); drawCornerOnImage(image, cornerMap); namedWindow("result"); imshow("result", image); waitKey(); return 0; } void drawCornerOnImage(Mat& image, const Mat&binary) { Mat_<uchar>::const_iterator it = binary.begin<uchar>(); Mat_<uchar>::const_iterator itd = binary.end<uchar>(); for (int i = 0; it != itd; it++, i++) { if (*it) circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1); } }
現在我們得到的效果就比默認的函數得到的結果有相當的改善，如上面最右邊效果圖。
3. 多尺度Harris角點
3.1 多尺度Harris角點的原理
雖然Harris角點檢測算子具有部分圖像灰度變化的不變性和旋轉不變性，但它不具有尺度不變性。但是尺度不變性對圖像特征來說至關重要。人們在使用肉眼識別物體時，不管物體遠近，尺寸的變化都能認識物體，這是因為人的眼睛在辨識物體時具有較強的尺度不變性。在圖像特征提取：尺度空間理論這篇文章里就已經講到了高斯尺度空間的概念。下面將Harris角點檢測算子與高斯尺度空間表示相結合，使用Harris角點檢測算子具有尺度的不變性。
仿照Harris角點檢測中二階矩的表示方法，使用$M=\mu(x,\sigma_I,\sigma_D)$為尺度自適應的二階矩：
$$\boldsymbol{M}=\mu(x,\sigma_I,\sigma_D)=\sigma_D^2g(\sigma_I)\otimes\begin{bmatrix}L_x^2(x,\sigma_D)&L_xL_y(x,\sigma_D)\\L_xL_y(x,\sigma_D)&L_y^2(x,\sigma_D)\end{bmatrix} $$
其中，$g(\sigma_I)$表示尺度為$sigma_I$的高斯卷積核，$x$表示圖像的位置。與高斯測度空間類似，使用$L(x)$表示經過高斯平滑后的圖像，符號$\otimes$表示卷積，$L_x(x,\sigma_D)$和$L_y(x,\sigma_D)$表示對圖像使用高斯$g(\sigma_D)$函數進行平滑后，在$x$或$y$方向取其微分的結果，即$L_x=\partial_xL$和$L_y=\partial_yL$。通常將$\sigma_I$稱為積分尺度，它是決定Harris角點當前尺度的變量，$\sigma_D$為微分尺度或局部尺度，它是決定角點附近微分值變化的變量。顯然，積分尺度$\sigma_I$應該大於微分尺度$\sigma_D$。
3.2 多尺度Harris角點實現
首先，檢測算法從預先定義的一組尺度中進行積分尺度搜索，這一組尺度定義為
$$\sigma_1\dots\sigma_n = \sigma_0\dots k^n\sigma_0$$
一般情況下使用k=1.4。為了減少搜索的復雜性，對於微分尺度$\sigma_D$的選擇，我們采用在積分尺度的基礎上，乘以一個比例常數，即$\sigma_D=s\sigma_I$，一般取s=0.7。這樣，通常使用積分和微分的尺度，便可以生成$\mu(x,\sigma_I,\sigma_D)$，再利用Harris角點判斷准則，對角點進行搜索，具體可以分兩步進行。
1. 與Harris角點搜索類似，對於給定的尺度空間值$\sigma_D$，進行如下角點響應值計算和判斷：
$$cornerness = det(\mu(x,\sigma_n)-\alpha trace^2(\mu(x,\sigma_n)))>threshold_H$$
2. 對於滿足1中條件的點，在點的8鄰域內進行角點響應最大值搜索（即非最大值抑制）出在8鄰域內角點響應最大值的點。對於每個尺度$\sigma_n(1,2,\dots,n)$都進行如上搜索。
由於位置空間的候選點並不一定在尺度空間上也能成為候選點，所以，我們還要在尺度空間上進行搜索，找到該點的所謂特征尺度值。搜索特征尺度值也分兩步。
1. 對於位置空間搜索到的每個候選點，進行拉普拉斯響應計算，並滿足其絕對值大於給定的閾值條件：
$$F(x,\sigma_n) = \sigma_n^2|L_{xx}(x,\sigma_n)+L_{yy}(x,\sigma_n)| \ge threshold_L$$
2. 與鄰近的兩個尺度空間的拉普拉斯響應值進行比較，使其滿足：
$$F(x,\sigma_n)>F(x,\sigma_l), \ \ \ l\in \{n-1.n+1\}$$
滿足上述條件的尺度值就是該點的特征尺度值。這樣，我們就找到了在位置空間和尺度空間都滿足條件的Harris角點。
多尺度Harris角點檢測C++實現：https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harrisLaplace.cpp
4. 更多的討論
在上面描述的Harris角點具有光照不變性、旋轉不變性、尺度不變性，但是嚴格意義上來說並不具備仿射不變性。Harris-Affine是一種新穎的檢測仿射不變特征點的方法，可以處理明顯的仿射變換，包括大尺度變化和明顯的視角變化。Harris-Affine主要是依據了以下三個思路：
- 特征點周圍的二階矩的計算對區域進行的歸一化，具有仿射不變性；
- 通過在尺度空間上歸一化微分的局部極大值求解來精化對應尺度；
- 自適應仿射Harris檢測器能夠精確定位牲點；
這篇文章不對Harris-Affine作進一步的描述，有時間會對這一算法做專門的分析，有興趣的可以參考原論文：Scale & Affine Invariant Interest Point Detectors.
5. 參考資料
[1] 《圖像局部不變特征與描述》王永明，王貴錦。
[3] 圖像特征提取PPT
[4] Harris角點檢測算法 1