9.1 角點偵測與Harris演算法#

概念說明#

什麼樣的影像位置是好的特徵點？把一個小視窗在影像上微微平移 (u, v)，比較視窗內容的變化量：

E(u, v) = Σ_{(x,y)∈W} w(x,y) · [I(x+u, y+v) − I(x, y)]²

三種地形的行為截然不同：

平坦區：怎麼移 E 都 ≈ 0——沒有資訊，不是特徵。
邊緣：沿邊緣滑動 E ≈ 0，垂直邊緣滑動 E 大——只有一個方向有變化，定位有歧義（沿著邊走到哪都像）。
角點：任何方向移動 E 都大幅上升——兩個方向都「不安分」，這才是能唯一定位的特徵。

Harris 演算法 (1988) 把 E(u,v) 對小位移做一階泰勒展開，整理成二次型：

E(u, v) ≈ [u v] · M · [u v]ᵀ，　M = Σ w · [[Iₓ², IₓI_y], [IₓI_y, I_y²]]

M 稱為結構張量 (structure tensor) 或二階矩陣。它的兩個特徵值 λ₁、λ₂ 正好編碼了三種地形：

λ₁、λ₂ 都小 → 平坦
一大一小 → 邊緣
兩個都大 → 角點

直接算特徵值慢，Harris 用一個巧妙的單值代理（det M = λ₁λ₂，trace M = λ₁+λ₂）：

R = det(M) − κ·trace(M)²，κ ≈ 0.04 ~ 0.06

R 大 → 角點；R ≈ 0 或負 → 平坦/邊緣。只用加減乘就能分類，這是 Harris 流行四十年的原因。

性質：(1) 對旋轉不變——M 的特徵值不隨座標旋轉而變；(2) 對亮度平移不變——微分把常數消掉；(3) 不對尺度不變——放大影像後同一角點的 R 會變，需要多尺度（這正是 SIFT 的動機，見下一節）。(4) Harris 對模糊很敏感，宜在適度平滑的影像上計算梯度。

實務流程尾端通常做非極大值抑制（R 圖上取局部極大）與門檻，只保留最顯著的 N 個角點。

演算法步驟#

影像轉 float 灰階，可先輕度高斯平滑。
用 Sobel 求 Iₓ、I_y。
對 Iₓ²、I_y²、IₓI_y 各做高斯（或盒狀）加權求和，得 M 的三個分量場。
逐像素算 R = det − κ·trace²。
門檻 R > T，再做非極大值抑制（3×3 內取最大）。
輸出角點座標清單（可依 R 排序取前 N）。

Python 程式碼範例#

import numpy as np

def gauss(img, k=5, sigma=1.0):
    x = np.arange(k)-k//2
    g = np.exp(-x**2/(2*sigma**2)); g /= g.sum()
    t = np.apply_along_axis(lambda r: np.convolve(r, g, "same"), 1, img)
    return np.apply_along_axis(lambda c: np.convolve(c, g, "same"), 0, t)

def conv3(img, K):
    P = np.pad(img, 1, mode="edge")
    out = np.zeros(img.shape)
    for y in range(img.shape[0]):
        for x in range(img.shape[1]):
            out[y,x] = (P[y:y+3, x:x+3]*K).sum()
    return out

def harris(img, k=0.05, T=1e5, top=10):
    I = img.astype(np.float64)
    Ix = conv3(I, np.array([[-1,0,1],[-2,0,2],[-1,0,1]]))
    Iy = conv3(I, np.array([[-1,-2,-1],[0,0,0],[1,2,1]]))
    Ixx, Iyy, Ixy = gauss(Ix*Ix), gauss(Iy*Iy), gauss(Ix*Iy)
    det = Ixx*Iyy - Ixy**2
    tr  = Ixx + Iyy
    R = det - k*tr**2
    # 非極大值抑制
    H, W = R.shape
    pts = []
    for y in range(2, H-2):
        for x in range(2, W-2):
            if R[y,x] >= T and R[y,x] == R[y-2:y+3, x-2:x+3].max():
                pts.append((R[y,x], x, y))
    pts.sort(reverse=True)
    return pts[:top], R

# 測試影像：一個矩形（4 個角點）+ 一條直線邊
img = np.zeros((40, 40))
img[10:30, 8:26] = 200
img[:, 36:] = 150          # 垂直邊（非角點）

pts, R = harris(img, T=5e4)
print("偵測到的角點（R, x, y）:")
for r, x, y in pts:
    print(f"  R={r:9.0f}  ({x}, {y})")
print("矩形理論角點: (8,10) (25,10) (8,29) (25,29)")
ys, xs = np.where((img > 100))
print("垂直邊緣中點 R =", R[20, 36].round(0), "（應遠小於角點 R）")

執行後：Harris 找出的點會聚集在矩形的四個角附近（含 R 值排序），而右側的長直邊——梯度很大但只有單方向變化——的 R 值接近 0 或為負，被結構張量的兩個特徵值一大一小的特徵自然排除。這正是「角點 vs 邊緣」判別力的直接數值展示。

複雜度分析#

三個 Sobel 卷積 O(N)、三次高斯平滑 O(N)（可分離）、R 計算 O(N)、NMS O(N)（樸素 25 鄰居最大值）。總計 O(N)，常數約 10~20 次卷積。空間 O(N)（Ix、Iy、三個矩、R）。工程上：κ 常取 0.04；門檻 T 對不同影像需調整，常用「R 最大值的 1%」這類相對門檻；需要次像素精度時可對 R 做二次曲面擬合。想清楚一件事：Harris 的輸出是「點清單」，丟掉了尺度與方向資訊——補上這兩件事，就是 SIFT 的故事。

練習#

把測試影像旋轉 30° 後再跑 Harris，比較角點 R 值是否大致不變（旋轉不變性），但座標已跟著轉。
將影像放大 2 倍（雙線性）再跑，觀察 R 門檻需要調整——這驗證「Harris 不具尺度不變性」。
實作 Shi–Tomasi 準則（min(λ₁, λ₂) = (tr − √(tr²−4det))/2 取角點），與 Harris 比較對「弱角點」的檢出率。