Harris Corner Detector

Introduction

本次作業目的在於介紹Harris Corner Detector的原理，總共分四大部分說明：

Moravec Detector: 說明Moravec Detector的操作方法與限制
Mathematics: 說明Harris Corner Detector的數學原理
Algorithmic Procedure: 說明實踐Harris Corner Detector的流程
Invariance Property: 說明scale invariance, rotational invariance and affine invariance三項特性

Moravec Detector

points of interest

Hans P. Moravec 在1977年在他的研究中提出Moravec operator，在研究中定義了"points of interest " 為在每個方向中都有很大的強度變化的點，Moravec operator被認為是一種 corner detector。

flat region

若移動參考區域，每次求得的梯度值不會有太多的變化

edge region

任意移動參考區域，會發現在垂直移動時梯度值沒有變化，其他方向會有變化

corner region

往任意方向移動，梯度值均會發生變化

Algorithm

Moravec operator 的輸入影像為grayscale image，有以下四步驟：

計算window 的強度變化，一共需計算8個位移量，分別為: (1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (-1, -1)
求出V(x, y) 裡的最小值
V(x, y) 裡的最小值令為Threshold 值，並將影像中小於Threshold 值得pixel另為0
使用non-maximal suppression找區區域內的最大值

Problem

Moravec operator有一些限制：

window funtion 為binary window，意思是window裡的數值為1 ，window外的數值為0，所以計算時binary window 內的權重都一樣，當影像的雜訊很大時，會將雜訊誤判為corner。
考慮的偏移數量太少
因為令Threshold 值為V(x, y) 裡的最小值，所以找到的點可能是edge。

因此，Harris Corner Detector將會解決以上三點限制。

Moravec operator：

Harris Corner Detector：

Mathematics

Corner

首先，需先介紹何謂 Corner及Corner的特性。

Corner 意指兩條邊緣相交處，在一個參考區內會有兩個不同方向的顏色變化，欲了解何處有兩條邊緣相交，就需求得影像的梯度(gradient)，找出梯度變化大的部分。

對影像作梯度計算後從可以分為三類，如右圖：

flat region: 計算梯度後，不論是對x方向做計算還是對y方向做計算，參考區域內均沒有突出的線條，只有由雜訊造成的細小線條。
edge region: 計算梯度後，兩張結果圖中有一條strong positive value 的對角線
corner region: 因為梯度計算有方向性，所以在對x方向和y方向計算後，會產生各一條strong positive value和strong negative value。

了解以上三種區域間不同的差異後，接下來就依照這個原理求得每個參考區域間的差異，需求sum of square difference，若sum of square difference值越大，代表每個參考區域的差異度越高，代表越接近corner。

Sum of square difference

以下開始求參考區域的sum of square difference，列式如下：

Taylor Series for 2D Functions

在計算過程中會使用到Taylor Series，用來求得近似函式，將式子以Taylor Series展開如下：

接下來將以上的算式整理一下，以Taylor Series簡化，再將變化後的式子以矩陣形式表示：

將整理完的式子帶回原式中：

Fitting Ellipse to each Set of Points

以下要來分析E(u, v)，在前面有提到sum of square difference值越大越接近corner，因此以下欲藉由E(u, v)找出參考區式位在flat region、edge region、corner region哪一個部分。

首先將上述E(u, v)計算出來，並繪製出圖形如下圖，可以觀察到，flat region像數間差異最小，所以集中在原點，edge region則呈現像是線性的分布，可以大致判斷出邊緣的走向，corner region有大的Ix 值和大的Iy值，和在CORNER 裡分析原圖的結果一致。

flat region edge region corner region

接下來，要分析數值的分布情形。E(u, v)可以近似於quadratic functions，取出E(u, v)矩陣的部分如右圖，因為E(u, v)是求sum of square difference結果為正數，矩陣M裡的數字同為正數，而將矩陣M取轉至矩陣後不變，M矩陣為對稱矩陣，因此矩陣M為一positive define matrix，若為positive define matrix函數圖形為鐘形。

我將問題簡化，以便繪製出立體的鐘形函數，假設梯度在像素上只能是垂直或水平的，代表IxIy = 0。我使用Maltab將flat region(Ix、Iy 數值小)、edge region (Ix、Iy 數值一大一小)、corner region Ix、Iy 數值大)各繪製一張圖，並繪製出函數的等高線圖。

驗證取近似於quadratic functions的圖形確實為橢圓，flat region為軸長短且接近圓形，edge region為細長的橢圓，corner region為軸長長的橢圓且接近圓形。

A = [1 0; 0 1]

A = [1 0; 0 4]

A = [4 0; 0 1]

A = [4 0; 0 4]

因為矩陣M為對稱矩陣，可以使用spectral decomposition，再將式子整理成橢圓標準式來表示，求得橢圓的長軸和短軸的對大值。

將計算完的橢圓疊加在sum of square difference數值分布圖上，可以清楚的知道橢圓跟書值分布的關係，若sum of square difference越大像數間的差異越大，橢圓的軸長越大。

flat region

edge region

corner region

conclusion

經過一連串的推倒知道可以由eigenvalue來判斷參考區域位於plat region、edge region、corner region哪個區域。若兩個eigenvalue數值很小則位於plat region，若其中一個eigenvalue遠大於另一個則位於edge region，若兩個eigenvalue均很大則位於corner region。

Algorithmic Procedure

實踐Harris Corner Detector的方法可分為以下五步驟：

此圖為以下步驟的範例圖片。

step1: Compute image gradients over small region

對window裡的每個pixel計算x 方向和y方向的梯度(gradient)，為解決noisy response，window functuon為gaussan function。原理上是將window 和一階微分的Gaussian做convolution，但在實際操作上運算過於龐大，所以在實作時，是將window 和 Sobel 做 convolution，求得Ix和Iy ，再將兩數相乘求Ix2、Iy2、IxIy。

step2: Compute the Autocorrelation Matrix

這步驟目的在於求出在 Mathematics 推倒出的M矩陣，這步驟會很耗時，因為需將Ix2、Iy2、IxIy分別對Gaussian Weight Function做convolution。對Gaussian Weight Function做convolution的會得到總和，而Gaussian Weight Function的特色是中心權重較重，周圍權重較低。

step3: Compute corner response function

依照原理，需求出eigenvalue，但求出eigenvalue的計算量過於龐大，所以改為求出 corner response function。將 Mathematics 中原先是使用eigenvalue來判斷window 是位於plat region、edge region、corner region哪個區域，改由corner response function R 來判斷的話圖形如下，區域相對位置沒有大改變，若R很小則位於plat region，若R小於0則位於edge region，若R很大則位於corner region。

step4: Threshold R

設定Threshold 值，將R值大於Threshold 值改為白色，R值小於Threshold 值改為黑色繪製成二元影像。

step5: Compute Non-Maximum Suppression

這個步驟的目的在找出區域最好的一個corner點，因為做convolution所以corner四周的數值都會很高，但僅需找出最佳的點作為corner。以為範例結果圖，左圖為最後求出的corner點，右圖為將corner 點和原圖片疊加後驗證結果。

Invariance Property

Classification

Geometry
- 1. Rotation
  2. scale
  3. Affine (scale dependent on direction)

Photometry
- 1. Affine intensity change (I --> a I + b)

scale invariance

scale invariance 代表說影像的強度(intensity)不會根據影像的大小改變，而是根據選擇區域改變，也就是說選擇同一個區域做縮放，區域內的平均強度(average intensity)不會改變。如下圖，參考區域選擇白色房子所在的位置，並將參考區域縮小，會發現強度分布圖上升至最高點、由最高點下降，兩面的坡度變得陡峭，根據這個特性，對圖片縮放使用Harris Corner Detector仍可以有效地找到corner，甚至對圖片縮放可以讓鋒值突出。

使用Harris Corner Detector時會出現一些問題。首先，若將參考區域選擇太小，容易受到雜訊的影響而選擇到過多的corner。一個良好的patch應該具有突出的峰值，且在patch裡只有一個峰值。再來，當影像放大後，原本的一個pixel會變成多個pixel，並且原本的corner部分可能會被偵測為edge。因此，選擇適當的縮放值對Harris Corner Detector的結果有很大的影響。

因此需在Harris Corner Detector裡加入選擇scale 的機制，有兩種方法：

Harris-Laplacian : 使用 Laplacia，將輸入影像和Laplacian做convolution
SIFT (Lowe's Scale Invariant Feature Transfrom): 使用Difference of Gaussians (DoG)，將輸入影像和DoG做convolution

因為求autocorrelation matrix 對noise的敏感度(sensitivity)很高，而這兩種方法裡都含有Gaussian可以濾除雜訊，調整不同的sigma 以選擇到適當的scale。

rotational invariance

對影像做旋轉後，計算gradient後得到的量值不變但方向(orientation)會改變，當量值不變計算出來的eigenvalue仍不會改變，Corner response 同樣不變，所以Harris Corner Detector有rotational invariance的特性。

在應用方面，若要兩張影像做image matching，但拍攝角度不同，物件被旋轉造成特徵點的方向改變，在偵測時就沒辦法被徵側為同一個特徵點，因為Harris Corner Detector有rotational invariance的特性，因此為了提高image matching的精準度和特徵點的辨識度，需將特徵點做orientation estimation。實踐orientation estimation的方法是平均特徵點和周圍區域的gradient，將特徵點和Gaussian weighting function做convolution。

affine invariance

affine transform 分為兩種，一種是幾何(geometry)上的轉變，意思是將影像做旋轉和非等比例縮放，而這兩種在前面已經討論過；另一種是光學(photometry)上的轉變，改變影像的intensity，對於Harris corner detector，將影像調亮並不會改變corner被偵測出的機率和位置，因此將影像調亮是invariance，但調整影像的對比(contrast)，就非invariance。

將影像的intensity增強是屬於affine invariance，但如左圖所示，若將intensity增強但threshold 維持，則會出現更多特徵點，因此threshold 需隨著intensity增強需做增加。threshold的最小值則被定義為最大穩定度(maximally stable )，還需求出maximally stable extremal regions(MSER)，只要在這個範圍內，無論幾何上的轉變還是光學上的轉換均屬於affine invariance。

Conclusion

這次作業遇到幾項困難。首先是在收集資料的部分，因為我想確保資料的正確性，所以參考資料很多是學校的教學投影片，或是一些網路上的公開課，但無論是上課投影片還是網路上的教學影片，都是以較精簡的方式呈現，所以在看數學式的時候會遇到很多不理解的地方，覺得算式前後不連貫，甚至不同老師的切入點不同，會使用不同的方式去講解，再加上大一線性代數也已經忘了很多基礎概念，所以這次作業以我理解的思考邏輯來呈現。

我這次也嘗試了使用ChatGPT來解決問題，會使用是因為有些資訊在網路上找不到完全解決我的疑問的答案，我之前有聽過一個關於ChatGPT的演講，講師有提到ChatGPT應該是補助工具並非學習工具，所以我將問題先給ChatGPT等真的釐清我的問題，或是找到一個解決方法時，再將關鍵字重新在google上查詢。一開始ChatGPT回答的答案離我的疑問很遠，但在對答的過程，發現我也越來越清楚自己的問題點，往往回到google上查詢時能更快能找到答案。

Reference

Corner Detection | Edge Detection , Shree K. Nayar, First Principles of Computer Vision, Mar 3, 2021, YouTube

Geometric Properties | Binary Images , Shree K. Nayar, First Principles of Computer Vision, Mar 3, 2021, YouTube

A Comparative Study between Moravec and Harris Corner Detection of Noisy Images Using Adaptive Wavelet Thresholding Technique , Nilanjan Dey, Pradipti Nandi, Nilanjana Barman , Debolina Das, Subhabrata Chakraborty /International Journal of Engineering Research and Applications (IJERA) ISSN: 2248-9622 www.ijera.com Vol. 2, Issue 1, Jan-Feb 2012, pp.599-606

Corner detection , Basics of Image Processing , Vincent Mazet (Université de Strasbourg), 2020-2023

7.7. Application: Conics and Quadric Surfaces , 13 June, 2019

Expressing a quadratic form with a matrix , Khan Academy , 2017

Harris corner detector , OpenCV 4.6.0, Open Source Computer Vision

Lecture 06:Harris Corner Detector , Robert Collins , CSE486, Penn State

Master MATLAB: visualize the matrix quadratic form , Mike X Cohen , 2019 , YouTube

Features Digital Visual Effects , Yung-Yu Chuang

Lecture 9: Edge + Corner Detection , Justin Johnson , February 6, 2020

Harris Corners , Kris Kitani , Carnegie Mellon University

二次型的意义是什么？有什么应用？

An Analysis and Implementation of the Harris Corner Detector , Javier Sánchez, Nelson Monzón , Agustín Salgado , Published in Image Processing On Line on 2018–10–03

Gaussian Weight Function , Marcus Brubaker , Feb 3, 2006

Harris Corner Detection Explained , baeldung , August 19, 2022

Feature detection and matching , Computer Vision: Algorithms and Applications, by R. Szeliski, Springer, 2011

Scale Invariant Detection 3 , Udacity , 2015 , YouTube

Scale Invariant Detection 4 , Udacity , 2015 , YouTube

MSER (Maximally Stable Extremal Regions) , P. Forss´en, D. Lowe, S-H Wang

Invariant Features , By Estrada, Jepson, and Fleet.