積分影像 (Integral Image)

由於模板內的影像要做積分後才能用於海森矩陣的計算，若沒有積分影像，則需進 行相當耗時的連加運算。積分影像[32]之目的在於加速矩形遮罩的摺積 (convolution) 運算，當矩形遮罩越大，重複運作在影像中位置相近的區域時，累加次數和範圍重疊的 部分就會越多，當累加計算的部分無法重複利用，就會造成浪費。

為了解決此問題，可透過積分影像的方式來加速，利用同一張影像，積分只需計算 一次的特性，矩形遮罩的摺積運算則可重複使用區塊像素累加成果，藉此大幅降低計算 量。如圖 3-8 所示，其中圖 3-8 (a) 代表從左上點往右下點所延伸之一組長方形深色面 積，深色面積內的灰階總和紀錄在右下紅點，其計算如下式：

𝐼_∑(𝑥, 𝑦) = ∑ ∑ 𝐼(𝑖, 𝑗)

𝑗≤𝑦

𝑗=0 𝑖≤𝑥

𝑖=0

(3.8)

當影像積分建立後，就可以利用查表的方式，快速拼湊出任意矩行區域的像素總和。如 圖 3-9 所示，然而其中的𝐼_∑1= ∑𝐵_𝑙1、𝐼_∑2 = ∑(𝐵_𝑙1+ 𝐵_𝑙2)、𝐼_∑3 = ∑(𝐵_𝑙1+ 𝐵_𝑙3)和𝐼_∑4 =

∑(𝐵_𝑙1+ 𝐵_𝑙2+ 𝐵_𝑙3+ 𝐵_𝑙4)若要計算區塊𝐵_𝑙4的像素總和，可藉由積分影像中位置的值𝐼_∑1、 𝐼_∑2、𝐼_∑3和𝐼_∑4得出：區塊𝐵_𝑙4像素總和 = 𝐼_∑4 − 𝐼_∑3− 𝐼_∑2+ 𝐼_∑1。

(a) (b) (c)

圖 3-8 積分影像示意圖 (a) 影像積分概念 (b) 影像灰階值 (c) 影像積分值

18

19

20

圖 3-10 盒濾波器替代高斯二階微分運算 (摘自[3])

圖 3-11 尺度空間構造示意圖(摘自[33])

SURF 採用9 × 9的盒濾波器作為初始尺度，產生之影像當成最底層的尺度，越往上 層濾波器大小也跟著增加，15 × 15、21 × 21、27 × 27 等等，如此可產生多組尺度空 間，以獲取更多特徵，過程中使用線性內插法，降低階和階之間 (octave-to-octave) 的尺 度變化，避免尺度取樣因變化大而過於粗糙。獲得尺度空間後，從這些尺度中搜尋可能 的特徵點，也就是搜尋空間中的極值(extrema)，每層影像上的每個點對其範圍內之對應 點，即鄰近3 × 3 × 3的空間中比較海森矩陣行列式的值，如圖 3-12，若該點的值大於鄰 近 26 個點，即判斷為特徵點。

21

圖 3-12 特徵點定位示意圖 3.3.3 SURF 敘述子

找到特徵點後，則需要對特徵點進行描述，由於尺度、旋轉不變特性，須找到特徵 點之主方向(major orientation)。SURF 採用哈爾小波濾波器 (圖 3-13) 來進行水平與垂 直方向的反應值計算，由計算出的反應值來訂定特徵點之方向。

圖 3-13 哈爾(Haar)小波濾波器

首先以特徵點為圓心，半徑大小 6𝜎(𝜎為特徵點所在尺度值)之圓形區域內，對每一 點進行哈爾小波遮罩運算(小波邊長為 4𝜎)，運算圓域內向訴之水平與垂直的一階導數 值，再以圓心角π/3的扇形區塊進行掃描一圈，統計同一掃描角度內之特徵點對於哈爾 小波的反應值，擁有最大的加總值的角度即為特徵點的主要方向。圖 3-14 所示，斑點 為哈爾小波遮罩運算一階導數值，灰色為扇形掃描視窗，箭頭長度為經過加總後的大小，

22

箭頭方向則表示特徵點方向。

圖 3-14 主方向示意圖

在獲得特徵點定位及主方向後，對特徵點建立特徵向量。以特徵點為中心，建立一 個邊長為 20σ之正方形區域，將正方形區域切割成4 × 4個方形子區域，在每個子方型區 域中以邊長5σ為取樣，並將哈爾小波濾波器以邊長為2σ進行運算，以左上角第一個起 點像素點起，每隔σ計算一次，在一個子區域內共計算了 25 次。定義𝐾_𝑥為哈爾小波濾波 器的水平反應值，𝐾_𝑦為哈爾小波濾波器的垂直反應值，會依照所選擇的特徵點而決定水 平與垂直的方向，再將𝐾_𝑥與𝐾_𝑦各別加總為∑ 𝐾_𝑥及∑ 𝐾_𝑦，為了得到向量長度取之絕對值其 分別加總為∑|𝐾_𝑥|、∑|𝐾_𝑦|共四個量值。有 16 個子區域，總計向量維度為4 × 4 × 4 = 64，

如圖 3-15 所示，最終得到 SURF 特徵點的特徵向量V_𝑖 = (∑ 𝐾_𝑥, ∑ 𝐾_𝑦, ∑|𝐾_𝑥| , ∑|𝐾_𝑦|)

23

圖 3-15 敘述子產生示意圖

3.4 特徵點匹配

所擷取到的影像序列經由 SURF 特徵點偵測演算法運算後，即可得到特徵點序列，

為能透過特徵點取得影像上的變化，首先要進行相鄰影像間的特徵點匹配。然而，直接 以特徵值進行匹配的結果含有錯誤的匹配與不可靠的對應關係，因此需要使用多種方法 進行篩選。

我們首先在特徵空間中進行初步的匹配，定義𝐹及𝐹′為兩張連續影像的特徵集合。

在兩張影像中特徵點χ ∈ F對應到χ′ ∈ F′之關係可表示為：

𝜒^′= argmin

𝜒 ̆𝜖𝐹

𝑑(𝜒, 𝜒^∗) (3.14)

其中𝑑：𝐹 × 𝐹^′→ ℝ 表示為一個選定的方法來衡量兩個特徵之間的相似性。對於不同型 態的特徵點敘述子，所採用的檢測相似度方法亦有所不同。對於二元型態的特徵點敘述

24

子 (binary feature descriptors) 使用漢明距離 (Hamming distance) 做為比對基礎；其它型 態的特徵點敘述子可能採用平方距離之和（sum of squared distances, SSD）或絕對距離 之和（sum of absolute distances, SAD）進行測量。在本研究中，在本研究中，在𝐹^′中採 用暴力搜尋 (brute force) 方法找出連續影像中的對應關係𝜒 → 𝜒^′。 改善特徵匹配後中，我們再以隨機抽樣一致（random consensus sampling, RANSAC）

技術[36]進行離群值排除。RANSAC 是一個非確定性的演算法，能從可能含有一定比例

25

定迭代次數結束後，RANSAC 以達成最大內群數量之模型作為最佳模型，並使用內群 數據再次最佳化該模型，以增進其準確性。其大致流程見圖 3-16。

在本研究的迭代運算中，是以本質矩陣 (essential matrix) 作為模型。匹配組合中的 兩個特徵點會有以下之關係：

𝑝̇𝐊^−⊤𝐄𝐊^−𝟏𝑞̇ < 𝜀_𝐹 (3.17)

點𝑝̇與𝑞̇為互相匹配的兩特徵點座標，𝐊為相機內部參數矩陣，𝐄為本質矩陣，𝜀_𝐹則為誤 差容許值，理想狀況下該值應該趨近於零。本質矩陣𝐄由選定的五個匹配組合基於五點 算法[22]進行評估而來，並使用該矩陣計算所有匹配組合之 Sampson 距離（詳述後章 4.2 節）作為篩選的閥值，此閥值於本研究設為 0.5 像素。

RANSAC 的優點是能夠對資料進行穩定的評估，即使在大量離群值時也可以精準 地估計其中的參數。然而該演算法是一個隨機過程，無法保證獲取的內群資料即為最佳 結果，且須謹慎定義參數、調整閥值，避免將過量或不足之資料剔除。

圖 3-16 RANSAC 架構圖 隨機抽樣設定內群

計算內群的模型

將陌生點帶入模型中 計算是否為內群

記下內群數量 重新計算模型

通過估計局內點與模型的錯誤率來評估模型

最好的模型即為結果

26

3.5 本體移動估計

給予 𝑡 − 1 影像之三維空間點雲(𝑥_𝑖, 𝑦_𝑖, 𝑧_𝑖)，與該點雲於 𝑡 時二維影像中的投影位 置 (𝑢_𝑖, 𝑣_𝑖)，在相機參數已知的情況下，估計兩個時間點的空間轉換關係，稱為 N 點透 視問題 (perspective-n-point problem, PnP)，如圖 3-17 所示。本研究以高效 N 點透視演 算法 (efficient perspective-n-point algorithm, EPnP)[35] 解出此問題。

EPnP 建構出超定性線性方程 (over-determined linear system) 以解出相機姿態。第 一步，在空間中選擇任意4 個能包含所有(𝑥_𝑖, 𝑦_𝑖, 𝑧_𝑖)的參考點{𝑐₁, 𝑐₂, 𝑐₃, 𝑐₄, }，三維空間點

27

觀察𝑁 ≥ 6個點與關係後，根據以上兩式可建構出齊次線性方程：

[

𝑓_𝑢𝛼₁^⊤ 0 (𝑢_𝑐− 𝑢₀)𝛼₁^⊤ 0 𝑓_𝑣𝛼₁^⊤ (𝑣_𝑐 − 𝑣₀)𝛼₁^⊤

⋮ ⋮ ⋮

𝑓_𝑢𝛼_𝑁^⊤ 0 (𝑢_𝑐− 𝑢₀)𝛼_𝑁^⊤ 0 𝑓_𝑣𝛼_𝑁^⊤ (𝑣 − 𝑣₀)𝛼_𝑁^⊤]

( Γ₁^′⊤

Γ₂^′⊤

Γ₃^′⊤

) = 0 (3.22)

該方程可表示為𝐀𝐱 = 0，其中未知數𝐱 = [Γ₁^′ Γ₂^′ Γ₃^′]^𝑇而𝐀為式 (3.22) 左側矩陣。該線性 方程可透過𝐴的奇異值分解 (singular value decomposition, SVD)，或者以更有效的方式計 算使𝐀^⊤𝐀𝜉 = 0的向量𝜉~[Γ₁^′ Γ₂^′ Γ₃^′]^𝑇。由於𝐀是一𝑁 × 12之矩陣，當𝑁非常大時，計算成 本也非常高，故改以計算𝐀^⊤𝐀使大小降為12 × 12進行有效率的運算。

要求出轉換矩陣𝐌，須先還原尺度係數。該係數可透過解開𝜉中的[𝑐₁^′ 𝑐₂^′ 𝑐₃^′ 𝑐₄^′]並與 參考點{𝑐₁, 𝑐₂, 𝑐₃, 𝑐₄, }比較彼此距離獲得。𝐌 可由 Horn 所提之封閉式[39]以關係𝑐_𝑗 → 𝑐_𝑗^′, ∀𝑗 ∈ {1,2,3,4}求出。

圖 3-17 N 點透視問題 PnP 示意圖（摘自[41]）

28

3.6 最佳化

以 EPnP 線性的方式進行計算相機移動參數𝐌 = (𝐑 𝐭)，產生的結果不夠穩健，且 誤差容易隨著移動的距離而累積，產生路徑偏移，因此需要進一步做最佳化調整，使重 投影誤差 (reprojection error) 最小化，如圖 3-18 所示，其中𝑃為影像平面上特徵點根據 相機內外部參數投影至三維空間中的座標，𝑝̇是特徵點在影像中的二維座標，當𝑃根據相 機參數投影至影像平面時，與𝑝̇之間的幾何距離差距，即為重投影誤差。依據重投影誤 差，定義出的最小平方目標如下式：

Φ(𝐑, 𝐭) = ∑‖(𝑢_𝑖, 𝑣_𝑖) − 𝜋_𝐊(𝐑 ∙ (𝑥_𝑖, 𝑦_𝑖, 𝑧_𝑖)^⊤+ 𝐭)‖²

𝑖

(3.23)

式子中𝜋_𝐊：ℝ³ → ℝ²是投影函數，使用相機矩陣𝐊將三維座標映射到二維平面座標的投 影函式，將所得到的齊次座標轉換為笛卡爾平面。

為求得式(3.23)的最小平方解，本研究使用萊文貝格-馬夸特演算法 (Levenberg-Marquardt algorithm, LM)[44]。LM 演算法為一種迭代演算法，能夠解決非線性最小化的 數值解，通常用於處理非線性規劃 (nonlinear programming) 以及最小平方曲線擬合 (least squares curve fitting) 問題。LM 演算法結合梯度下降法 (gradient descent algorithm) 與高斯-牛頓演算法 (Gauss-Newton algorithm) 優點，結合上述兩種方法並尋找最佳解，

梯度下降法收斂速度較緩慢，靠近極小值時速度減慢，也容易陷入區域解而失去找到最

佳解的機會；高斯-牛頓演算法則有較佳的收斂速度，但也會因收斂太快，錯過最佳解。

29

圖 3-18 重投影誤差示意圖（摘自[41]）

3.7 狀態更新

本研究在追蹤特徵點的觀察中，為了提高視覺測程與本體移動的穩健性，在相機移 動的整個影像序列中保持追蹤特徵的狀態。追蹤特徵的過程中，對每個追蹤中的特徵點 採用有效的整合技術，以維持特徵狀態的移動權重平均。對於真實狀態的計算，經由下 計算式：

𝐦̅_𝑖,𝑡 =𝜔̅_{𝑖,𝑡−1}⋅ 𝑓_{𝑡−1,𝑡}(𝐦̅_{𝑖,𝑡−1}) + 𝜔_𝑖,𝑡⋅ 𝐦_𝑖,𝑡

𝜔̅_𝑖,𝑡 (3.24)

令𝐦_𝑖,𝑡是第t張影像中第𝑖個特徵點的觀測狀態向量，本研究僅考慮特徵點之 3D 座標，而 𝜔_𝑖,𝑡𝜖[0, 1]為權重，表示測量狀態的可靠度。其中

𝜔̅_𝑖,𝑡 = 𝜔̅_{𝑖,𝑡−1}+ 𝜔_𝑖,𝑡 (3.25)

為移動權重 (running weigh)，𝑓_{𝑡−1,𝑡}為狀態的轉換函數，代表從前一個時間𝑡 − 1轉移到 目前時間點 𝑡。

30

第四章 多重模型本體移動估測

傳統基於特徵點匹配的本體移動估計，是採取最小化重投影誤差的策略，如前章 3.6 節所述。本章提出一新式模型，同時採用單眼視覺測程中常用之光度誤差 (photo-consistency)模型與對極幾何限制，以多重目標函數作為最佳化相機姿態之調整依據，強 化雙眼視覺測程對離群值與雜訊之容錯性，增強移動估測的穩健性。

4.1 投影誤差模型

如前章所述，給予 𝑡 − 1 影像之三維空間點雲(𝑥_𝑖, 𝑦_𝑖, 𝑧_𝑖)，與該點雲於 𝑡 時二維影 像中的投影位置 (𝑢_𝑖, 𝑣_𝑖)，重投影誤差之平方合目標函數定義為：

∅_RPE(𝐑, 𝐭) = ∑‖(𝑢_𝑖, 𝑣_𝑖) − 𝜋_𝐊(𝐑 ∙ (𝑥_𝑖, 𝑦_𝑖, 𝑧_𝑖)^⊤+ 𝐭)‖²

𝑖

(4.1)

式 (4.1) 可作為相機姿態估計 (𝐑, 𝐭) 精準與否之依據，並在點雲測量與投影觀測函數 誤差呈現高斯分布時，被證實能提供最大機率估計[45] (maximum likelihood estimation, MLE)。然而，由視差圖雜訊產生之三維點座標重新投影後可能造成極大誤差，往往高 於數十甚至數百像素，使得重投影誤差模型在數值運算上相對失準，因此常輔以而外限 制 (如極幾何[34]) 抑制此類雜訊。

4.2 極幾何模型

本研究除了對重投影誤差進行最佳化調整之外，使最佳化結果更為穩健，考慮了特 徵點二維影像的對應關係 (𝑢, 𝑣) ↔ (𝑢^′, 𝑣^′) 之極線誤差。如圖 4-1 所示，在影像𝐼_𝑡−1上

本研究除了對重投影誤差進行最佳化調整之外，使最佳化結果更為穩健，考慮了特 徵點二維影像的對應關係 (𝑢, 𝑣) ↔ (𝑢^′, 𝑣^′) 之極線誤差。如圖 4-1 所示，在影像𝐼_𝑡−1上