角檢測

角檢測（Corner detection）或興趣點檢測（interest point detection），是電腦視覺系統中用來提取特徵以及推測圖像內容的一種方法。角檢測的應用很廣，經常用在運動檢測、跟蹤、馬賽克（image mosaicing）、全景圖縫合（panorama stiching）、三維建模以及物體識別中。

## 問題定義

兩條邊的交點形成一個角（點）。而圖像的要點（也稱為受關注點）是指圖像中具有代表性以及穩健性（即指該點能夠在有噪聲干擾的情況下也能穩定的被定位，在大陸亦被稱為：魯棒性）的點。也就是說，要點可以是角（點），也可以不是，例如局部亮點或暗點，線段終點，或者曲線上的曲率最大值點。在實際應用中，很多所謂的（角）點檢測演算法其實是檢測要點，而不僅僅是角（點）。所以，如果我們只想檢測角的話，還需要對檢測出的要點進一步分析。

例如也可以先經過邊檢測，之後在做一些後處理來檢測角，像是Kirsch operator和Frei-Chen masking set兩種方法。為了分辨區辨識時圖像含有一部分目標圖像的資訊是否正確，角檢測通常不是一個非常確定且同時需要很多額外的確認才可以確定，角檢測在圖像上最簡單的做法是利用相關性，但這樣的作法需要耗費大量的運算資源以及得到的結果可能只是局部最大值，因此有另一個常見的作法是使用Harris和Stephens所提出的經由改善Moravec方法的結果。

## Moravec 角檢測演演算法

這是最早使用來做角檢測的做法，他首先定義所謂的「角」就是那些自我相似程度低的點。這個演演算法檢查所有圖像中的像素，並考慮以該像素為中心點的一片範圍，該範圍與他周圍覆蓋最大的另一個範圍的相似度做為參考，而相似度通常是將兩個範圍的對應點計算誤差的平方和(SSD: Sum of Squared differences)，越小代表相似度越高。

舉例來說，如果一個像素是位於一片均勻強度的區塊時，這時每個像素周圍與其附近的像素的周圍都是相當類似的，因此相似度也高，不會被認為是角。而在圖形的邊界處的像素，若是與邊垂直方向處附近的像素，兩者的周圍圖像相似度就會比較低，然而若是與邊平行附近的像素，兩者的周圍圖像相似度就會比較高(看到的都是一條在同樣位置的邊)，因此也不會被認為是角。只有在那些與附近像素的周圍圖像都很不相似的像素，才會被認為是「角」。

如果我們將這種與周圍的相似程度量化(使用誤差平方和)，並且找出整個圖像的局部最大值(局部最不相似的點)，這些局部最大值就很有可能是我們想要偵測到的「角」。

## Harris & Stephens / Plessey / Shi–Tomasi 角檢測演演算法

Harris & Stephens改善了Moravec的方法，他們直接考慮每個像素沿著特定方向處的像素的SSD，而不是使用與周圍像素範圍的SSD。

不失一般性，我們假設現在是對一個灰階的2維圖像來做偵測。令這個圖為 $I$，考慮 $(u, v)$ 位置的像素與周圍的區塊，並選取一個固定向量 $(x, y)$ 作為該像素的參考像素，因此固定的這個向量 $(x, y)$ 所算出的所有差平方和的總和記做為 $S(x,y)$，而對於每個 $(u,v)$ 做不同加權，就可以得到

$$S(x,y) = \sum_u \sum_v w(u,v) \, \left( I(u+x,v+y) - I(u,v)\right)^2$$

如果使用泰勒展開式將 $I(u+x,v+y)$ 展開，假設$I_x$ 和 $I_y$是 $I$ 的偏微分，可以得到

$$I(u+x,v+y) \approx I(u,v) + I_x(u,v)x+I_y(u,v)y$$

並將原式改寫成

$$S(x,y) \approx \sum_u \sum_v w(u,v) \, \left( I_x(u,v)x + I_y(u,v)y \right)^2,$$

或是使用矩陣來簡化算式

$$S(x,y) \approx \begin{pmatrix} x & y \end{pmatrix} A \begin{pmatrix} x \\ y \end{pmatrix},$$

其中矩陣 $A$ 為

$$A = \sum_u \sum_v w(u,v) \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \end{bmatrix} = \begin{bmatrix} \langle I_x^2 \rangle & \langle I_x I_y \rangle\\ \langle I_x I_y \rangle & \langle I_y^2 \rangle \end{bmatrix}$$

式子中的矩陣是Harris矩陣，而括號代表對所有的 $(u,v)$ 取加權平均。

一個角(或者說是要點)可以被刻劃成該點 $(u,v)$ 使得 $S(u,v)$ 在不同的方向 $(x,y)$ 之間有很大的變異數的點。

根據分析 $A$ 的特徵值，以上的特性可以用下面的方式闡述：在角的點上，$A$應該要有兩個"大"的特徵值。

根據特徵值的大小，我們會有以下的推論：

- 如果 $\lambda_1 \approx 0$ 和 $\lambda_2 \approx 0$ 則這個像素 $(u,v)$ 只是一個普通的點。
- 如果 $\lambda_1 \approx 0$ 但 $\lambda_2$ 為一個大的正數，則這個像素 $(u,v)$ 位在邊上面。
- 如果 $\lambda_1$ 和 $\lambda_2$ 都是大的正數，則這個像素 $(u,v)$ 就是一個角點。

Harris 和 Stephens 注意到計算特徵值在運算量上面相當大，計算中需要使用到取平方根的操作，因此給出了另一個近似值 $M_c$，其中 $\kappa$ 是需要調整的重要參數。

$$M_c = \lambda_1 \lambda_2 - \kappa \, (\lambda_1 + \lambda_2)^2 = \operatorname{det}(A) - \kappa \, \operatorname{trace}^2(A)$$

因此該演演算法不需要真正的去做 $A$ 的特徵分解，而只需要去估計 $A$ 的行列式和跡。

Shi–Tomasi角檢測在假設一般圖像每個像素所給出的函數值通常是光滑(smooth)且穩定的(stable)，他直接去計算 $\min(\lambda_1, \lambda_2)$，有時又稱該方法為 Kanade-Tomasi 角檢測。

而 $\kappa$ 的值通常是由經驗決定，而通常在 0.04–0.15 的範圍裡是可解的。

另一種方法可以避開選擇 $\kappa$ 的困擾，這個方法是使用 Noble's角測度 $M_c'$，而角測度是計算所有特徵值的調和平均數

$$M_c' = 2 \frac{\operatorname{det}(A)}{\operatorname{trace}(A) + \epsilon},$$

其中 $\epsilon$ 是一個小的正常數。

而角的位置的共變異數矩陣即為$A^{-1}$，可以寫成

$$\frac{1}{\langle I_x^2 \rangle \langle I_y^2 \rangle - \langle I_x I_y \rangle^2} \begin{bmatrix} \langle I_y^2 \rangle & -\langle I_x I_y \rangle\\ -\langle I_x I_y \rangle & \langle I_x^2 \rangle \end{bmatrix}.$$

## Förstner 角檢測

在某些情況，會希望更精確地去計算角的位置，為了得到近似值，Förstner演演算法可以解出閉集上的角附近範圍中的所有切線與最接近這些切線的點，該演演算法依賴於在一個理想的角附近。

一個像素 $\mathbf{x'}$ 的切線 $T_{\mathbf{x'}}(\mathbf{x})$ 可以由數學式給出

$$T_\mathbf{x'}(\mathbf x)=\nabla I(\mathbf{x'})^{\top}(\mathbf{x}-\mathbf{x'})=0$$

其中 $\nabla I(\mathbf{x'})=[I_{\mathbf{x}},I_{\mathbf{y}}]^{\top}$ 是圖像 $I$ 在 $\mathbf{x'}$ 點的梯度向量。

而最靠近長方形範圍 $N$ 中所有切線的點 $\mathbf{x}_{0}$ 為

$$\mathbf{x}_{0}=\underset{\mathbf{x}\in \mathbb{R}^{2\times 1}}{\operatorname{argmin}}\int_{\mathbf{x'}\in N}T_{\mathbf{x'}}(\mathbf{x})^{2}d\mathbf{x'}$$

$\mathbf{x}_{0}$ 到切線 $T_{\mathbf{x'}}$ 的距離依照梯度向量的大小來加權，因此若經過該點的切線有較大的梯度的話，在加權上就會占較高的比重。

嘗試著解 $\mathbf{x}_{0}$：

$$\begin{align} \mathbf{x}_{0}&=\underset{\mathbf{x}\in \mathbb{R}^{2\times 1}}{\operatorname{argmin}} \int_{\mathbf{x'}\in N}(\nabla I(\mathbf{x'})^{\top}(\mathbf{x}-\mathbf{x'}))^{2}d\mathbf{x'}\\ &=\underset{\mathbf{x}\in \mathbb{R}^{2\times 1}}{\operatorname{argmin}}\int_{\mathbf{x'}\in N}(\mathbf{x}-\mathbf{x'})^{\top}\nabla I(\mathbf{x'})\nabla I(\mathbf{x'})^{\top}(\mathbf{x}-\mathbf{x'})d\mathbf{x'}\\ &=\underset{\mathbf{x}\in \mathbb{R}^{2\times 1}}{\operatorname{argmin}}\, (\mathbf{x}^{\top}A\mathbf{x}-2\mathbf{x}^{\top}\mathbf{b}+c) \end{align}$$

$A\in\mathbb{R}^{2\times 2},\textbf{b}\in\mathbb{R}^{2\times 1},c\in\mathbb{R}$ 分別為

$$\begin{align} A&=\int \nabla I(\mathbf{x'})\nabla I(\mathbf{x'})^{\top}d\mathbf{x'}\\ \mathbf{b}&=\int \nabla I(\mathbf{x'})\nabla I(\mathbf{x'})^{\top}\mathbf{x'}d\mathbf{x'}\\ c&=\int \mathbf{x'}^{\top}\nabla I(\mathbf{x'})\nabla I(\mathbf{x'})^{\top}\mathbf{x'}d\mathbf{x'}\\ \end{align}$$

要找出最小值可以經由計算式子對 $x$ 的微分，並讓微分後的式子等於0來解$x$

$$2A\mathbf{x}-2\mathbf{b}=0 \Rightarrow A\mathbf{x}=\mathbf{b}$$

注意到如果要能夠解上述等式， $A\in\mathbb{R}^{2\times 2}$ 必須要是可逆的，或是說 $A\in\mathbb{R}^{2\times 2}$ 必須要是滿秩的，而解可以寫成

$$x_{0}=A^{-1}\mathbf{b}$$

只有在範圍 $N$ 之中有角時才存在。

此外 Lindeberg給出一個方法去評估角的確定性，經由最小化誤差

$$\tilde{d}_{min} = \frac{c - b^T A^{-1} b}{\mbox{trace} A}$$

而這對所有尺度都適用，因此這方法可以自動去適應圖片中不同梯度大小所造成的影響與誤差。

此外

- $c$ 可以被視為最小平方誤差計算中的誤差，因此當 $c$ 為0時，將沒有誤差存在。
- 這個演演算法可以用來找出環狀特徵的中心，只需將上述演演算法中的梯度向量改成法線向量即可。

## 多尺度 Harris 算子

Harris 算子中的微分矩陣 $A$ 會牽涉到要計算圖像中的偏微分(Image derivatives) $I_x, I_y$，而這可以用附近點的線性組合來做逼近，而計算線性組合時通常也會需要根據附近圖像值的大小來做縮放，因此需要兩個額外的縮放參數：(i)一個局部縮放參數主要用來讓圖像平滑一點 (ii)一個積分參數用來累積在微分操作中的非線性算子，作為圖像的修正項。

若以 $I$ 代表原本圖形中的亮度，$L$ 為一個將圖像 $I$ 通過一個高斯濾波器讓圖像變模糊的結果。

$$g(x, y, t) = \frac {1}{2{\pi} t}e^{-(x^2+y^2)/2t}$$

而參數 $t$ 用來控制高斯濾波器的變異數

$$L(x, y, t)\ = g(x, y, t) * I(x, y)$$

令 $L_x = \partial_x L$ 和 $L_y = \partial_y L$ 為 $L$ 的偏導數，此外用一個高斯濾波器作為縮放大小，參數 $s$ 可供調整。

則multi-scale second-moment matrix 可以被定義為

$$\mu(x, y; t, s) = \int_{\xi = -\infty}^{\infty} \int_{\eta = -\infty}^{\infty} \begin{bmatrix} <pre> L_x^2(x-\xi, y-\eta; t) & L_x(x-\xi, y-\eta; t) \, L_y(x-\xi, y-\eta; t) \\ </pre> L_x(x-\xi, y-\eta; t) \, L_y(x-\xi, y-\eta; t) & L_y^2(x-\xi, y-\eta; t) \end{bmatrix} g(\xi, \eta; s) \, d\xi \, d\eta.$$

則我們也可以用前面類似的方法來計算 $\mu$ 的特徵值，並且定義多尺度 Harris corner 測度為

$$M_c(x, y; t, s) = \operatorname{det}(\mu(x, y; t, s)) - \kappa \, \operatorname{trace}^2(\mu(x, y; t, s))$$

注意到局部尺度參數 t 和積分尺度參數 s，這兩個參數通常會有相對關係 $s = \gamma^2 t$，其中 $\gamma$ 為一個參數通常介於 $[1, 2]$ 之間。

因此我們可以計算多尺度 Harris corner 測度 $M_c(x, y; t, \gamma^2 t)$，在給定的任何局部尺度 t 去做到多尺度角偵測。

一般實際應用上，多尺度角偵測通常會先做尺度選擇的步驟，可以參考正規化尺度拉普拉斯算子。

$$\nabla^2_{norm} L(x, y; t)\ = t \nabla^2 L(x, y, t) = t (L_{xx}(x, y, t) + L_{yy}(x, y, t))$$

常常在尺度選擇的部分用到，並且用來同時調整尺度以適應該角點。

- 空間範圍內多尺度角測度 $M_c(x, y; t, \gamma^2 t)$ 的最大值：

$$(\hat{x}, \hat{y}; t) = \operatorname{argmaxlocal}_{(x, y)} M_c(x, y; t, \gamma^2 t)$$

- 將圖像通過一個高斯濾波器 $\nabla^2_{norm}(x, y, t)$ (拉普拉斯算子)所得到的 $L$ 的局部最大值或最小值：

$$\hat{t} = \operatorname{argmaxminlocal}_{t} \nabla^2_{norm}L(\hat{x}, \hat{y}; t)$$