國立高雄大學資訊工程學系
碩士論文
不同影像特徵點演算法之單眼視覺測程比較研究
A Comparative Study of Monocular Visual Odometry
Using Different Image Feature Identification Algorithms
研究生：莊振錡 撰
指導教授：陳佳妍 博士
i
致謝
回首過去、展望未來。入學的景象仍記憶猶新，這一路走來學習到了許多的 知識、技術以及寶貴的經驗，也獲得了不少人的幫助。感謝伴隨在這段路程中的 師長、前輩、同儕、晚輩以及家人，將帶著這份感恩與回饋的心，朝向更美好的 未來前進。 能夠順利的完成研究，首先感謝陳佳姸老師的指導，很榮幸能夠成為佳姸老 師的學生，不僅能夠在學術上有不少的收穫，也得到許多各方面的機會能夠去表 現並學到無價的經驗。感謝簡祥任學長的幫助與付出，很感謝在理論方面的解惑 與技術層面的協助，能夠清楚理解觀念且對於實作的技術與技巧又更進一步。感 謝視覺與圖學實驗室的成員，共同經歷了大大小小的事情，也豐富了在實驗室的 生活。感謝系辦的陳淑真助理總是適時的給予許多方面的協助，使得多數的事情 都能有效率及圓滿地完成。感謝家人一路上的鼓勵及關心，能夠無後顧之憂的進 行研究工作。感謝陳銘志教授與殷堂凱教授能夠擔任學位口試委員，給予建議與 修正，使得研究的內容與成果更為完善。 在這一路上要感謝的人實在太多了，單純用這樣的內容篇幅是無法表達的， 也許有人說：「那就謝天吧！」，但我想不僅僅是謝天，給予這樣的機會與緣分能 夠認識與得到這麼多人的幫助，更應該是帶著這份感恩，在未來當能力所及時， 也可以伸出雙手提供適時的協助。 最後，僅將這本研究論文獻給所有要感謝的每一位。 莊振錡 謹致於 國立高雄大學 資訊工程學系 視覺與圖學實驗室 西元 2016 年 12 月ii
不同影像特徵點演算法之單眼視覺測程比較研究
指導教授：陳佳妍 博士 國立高雄大學資訊工程學系 學生：莊振錡 國立高雄大學資訊工程學系 摘要 要使移動中的裝置知道本身的移動狀態，需要由本體移動估測來得知。近年 來，基於影像特徵的本體移動估測技術，持續在 VO (Visual Odometry)，V-SLAM (Visual Simultaneously Localisation and Mapping)和 SfM (Structure-from-Motion)的 發展中扮演重要關鍵。在計算相機移動狀態時，影像中特徵點的檢測以及描述， 將對計算的精確度以及成本發揮關鍵作用。在本研究中將探討不同形態特徵點對於單眼視覺測程的結果，其中包含四種 常見的演算法 SIFT、SURF、ORB 和 BRISK，以及近期提出的 A-KAZE 特徵演 算法。在實驗中，使用 The KITTI Vision Benchmark Suit 提供之影像資料，產生 不同型態的特徵點，利用多項限制條件保留良好的特徵匹配後，再從計算之三維 座標、對應的二維座標與相機參數透過有效 N 點透視參數估測演算法(Efficient Perspective-n-Point, EPnP) 以計算相機的移動狀態，接者對於相機移動狀態中的 旋轉與位移進行非線性最佳化調整，以降低估策過程中的累積誤差並獲得最佳之 估測結果。本研究提出一個混合兩種對應關係(2D-2D 及 3D-2D)的最佳化調整方 法，藉由各自的對應數量做為最佳化調整中之參考依據。 完成單眼視覺測程後，首先比較特徵點偵測所需的時間及空間等數據，再對 所估測出之結果與真實路徑之準確率進行各方面的比較。在比較結果中，可發現 SURF 及 SIFT 演算法具有較好的表現，ORB 演算法犧牲準確率但提供快速的計 算，而 BRISK 以及 A-KAZE 則是表現出準確率與計算成本的折衷。
iii
A Comparative Study of Monocular Visual Odometry
Using Different Image Feature Identification
Algorithms
Advisor: Dr. Chia-Yen Chen
Department of Computer Science and Information Engineering National University of Kaohsiung
Student: Chen-Chi Chuang
Department of Computer Science and Information Engineering National University of Kaohsiung
ABSTRACT
The determination of the motion of a device as it is moving is commonly known as ego-motion estimation. In recent years, ego-motion estimation approaches based on image features play important roles in visual odometry (VO), visual simultaneous localization and mapping (V-SLAM), as well as structure from motion (SfM). In the determination of the camera motions, the detection and description of feature points within images are crucial to the accuracy and cost of calculation.
In this work, we investigate five different approaches to feature point detection and description and how they affect the outcome of monocular visual odometry. The selected approaches are SIFT, SURF, ORB and BRISK, as well the recently proposed A-KAZE method. In our work, the KITTI Vision Benchmark Suit is used to provide image data, with various types of feature points detected by the selected method and multiple constraints utilized to retain good feature matching. The matched features are used to calculate corresponding 3D and 2D coordinates. The coordinates, combined with camera parameters, are used in Efficient Perspective-n-Point to determine the camera pose estimation. Then we apply a proposed non-linear optimization adjustment approach, which integrates two types of correspondence (2D-2D and 3D-2D) to reduce the accumulative error within the estimation process. After monocular visual odometry, we analyze the time and storage costs of the different feature detection approaches and compare the ego-motion estimated by each approach
iv
with the reference GPS data to determine the performance of each approach.
Keywords: Visual Odometry, Ego-Motion Estimation, SIFT, SURF, ORB, BRISK, A-KAZE
v
目錄
致謝 ... i 摘要 ... ii ABSTRACT ... iii 目錄 ... v 圖目錄 ... viii 表目錄 ... x 符號定義 ... xi 第一章 緒論 ... 1 1.1 前言 ... 1 1.2 研究目的 ... 1 1.3 研究方法 ... 2 1.4 論文架構 ... 2 第二章 文獻探討 ... 4 第三章 影像特徵點 ... 7 3.1 SIFT 特徵點演算法 ... 7 3.1.1 尺度空間極值偵測 ... 8 3.1.2 特徵點定位 ... 10 3.1.3 特徵點定向 ... 12 3.1.4 特徵點描述子 ... 12 3.2 SURF 特徵點演算法 ... 13 3.2.1 積分影像... 13 3.2.2 快速海森檢測 ... 15 3.2.3 SURF 敘述子 ... 18 3.3 ORB 特徵點演算法 ... 19vi 3.3.1 FAST 特徵偵測 ... 19 3.3.2 Harris 角點偵測 ... 20 3.3.3 Oriented FAST ... 23 3.3.4 BRIEF 敘述子 ... 24 3.4 BRISK 特徵點演算法 ... 26 3.4.1 尺度空間測徵點偵測 ... 26 3.4.2 採樣方法與旋轉估計 ... 27 3.4.3 BRISK 特徵點描述 ... 29 3.5 A-KAZE 特徵點演算法 ... 30 3.5.1 非線性擴散濾波 ... 30 3.5.2 快速顯式擴散 ... 31 3.5.3 非線性尺度空間 ... 33 3.5.4 特徵偵測... 34 3.5.5 特徵描述... 34 第四章 單眼視覺移動估測 ... 36 4.1 前言 ... 36 4.1.1 相機模型... 36 4.1.2 影像校正與極幾何限制 ... 41 4.2 特徵匹配 ... 42 4.2.1 暴力描述子匹配 ... 42 4.2.2 匹配結果檢測 ... 42 4.2.3 強健特徵匹配 ... 43 4.3 單眼視覺測程 ... 45 4.3.1 視覺測程起始 ... 45 4.3.2 本體移動估測 ... 45
vii 4.3.3 非線性最佳化 ... 47 4.3.4 深度復原... 51 4.3.5 狀態更新... 51 第五章 實驗與結果分析... 53 5.1 特徵點偵測演算法比較 ... 53 5.2 單眼視覺測程結果 ... 61 5.3 單眼視覺測程結果分析 ... 65 第六章 結論與未來方向... 83 6.1 結論 ... 83 6.2 未來研究方向 ... 84 參考文獻 ... 85
viii
圖目錄
圖 2-1 基礎視覺測程流程概略圖 ... 5 圖 3-1 高斯尺度卷積與差分影像[4] ... 9 圖 3-2 特徵點定位示意圖[34] ... 10 圖 3-3 積分影像概念圖 ... 14 圖 3-4 透過積分影像計算區域像素總和 ... 14 圖 3-5 盒濾波器簡化圖[34]... 16 圖 3-6 SURF 尺度空間示意圖[34] ... 17 圖 3-7 特徵點定位示意圖[34] ... 18 圖 3-9 測量點分割測試 ... 20 圖 3-10 特徵點三種型態示意圖 ... 21 圖 3-11 依據λ1及λ2關係所產生之情形 ... 23 圖 3-12 隨機取樣 7 個點對示意圖 ... 25 圖 3-13 BRISK 尺度空間檢測示意圖[7] ... 27 圖 3-14 BRISK 特徵點採樣模型[7] ... 28 圖 4-1 本研究所實作之單眼視覺測程流程圖 ... 36 圖 4-2 相機針孔模型[29] ... 37 圖 4-3 形變示意圖 ... 39 圖 4-4 具形變之相機投影模型[34] ... 40 圖 4-5 極幾合限制示意圖[10] ... 41 圖 4-6 RANSAC 流程示意圖[34] ... 44 圖 4-7 N 點透視問題(PnP)示意圖... 46 圖 4-8 重投影誤差示意圖[34] ... 48 圖 4-9 極幾何誤差示意圖 ... 49 圖 5-1 2011_09_26_drive_0002_sync 影像序列示意 ... 54ix 圖 5-2 2011_09_26_drive_0091_sync 影像序列示意 ... 54 圖 5-3 2011_09_26_drive_0095_sync 影像序列示意 ... 55 圖 5-4 2011_09_26_drive_0002_sync 特徵點示意 ... 57 圖 5-5 2011_09_26_drive_0091_sync 特徵點示意 ... 59 圖 5-6 2011_09_26_drive_0095_sync 特徵點示意 ... 61 圖 5-7 以 2011_09_26_drive_0002_sync 影像序列執行結果盒狀圖 ... 62 圖 5-8 以 2011_09_26_drive_0091_sync 影像序列執行結果盒狀圖 ... 63 圖 5-9 以 2011_09_26_drive_0095_sync 影像序列執行結果盒狀圖 ... 64 圖 5-10 2011_09_26_drive_0002_sync 估測結果誤差圖 ... 66 圖 5-11 2011_09_26_drive_0091_sync 估測結果誤差圖 ... 66 圖 5-12 2011_09_26_drive_0095_sync 估測結果誤差圖 ... 67 圖 5-13 2011_09_26_drive_0002_sync 估測結果分段誤差 ... 70 圖 5-14 2011_09_26_drive_0091_sync 估測結果分段誤差 ... 73 圖 5-15 2011_09_26_drive_0095_sync 估測結果分段誤差 ... 76 圖 5-16 2011_09_26_drive_0002_sync 移動軌跡比較 ... 77 圖 5-17 2011_09_26_drive_0002_sync 路徑放大圖 ... 78 圖 5-18 2011_09_26_drive_00091_sync 移動軌跡比較 ... 79 圖 5-19 2011_09_26_drive_0091_sync 路徑放大圖 ... 80 圖 5-20 2011_09_26_drive_0095_sync 移動軌跡比較 ... 81 圖 5-21 2011_09_26_drive_0095_sync 路徑放大圖 ... 82
x
表目錄
表 3-1 特徵點演算法比較 ... 7 表 5-1 2011_09_26_drive_0002_sync 特徵點相關統計比較 ... 54 表 5-2 2011_09_26_drive_0091_sync 特徵點相關統計比較 ... 55 表 5-3 2011_09_26_drive_0095_sync 特徵點相關統計比較 ... 55 表 5-4 2011_09_26_drive_0002_sync 估測結果比較表 ... 62 表 5-5 2011_09_26_drive_0091_sync 估測結果比較表 ... 63 表 5-6 2011_09_26_drive_0095_sync 估測結果比較表 ... 64 表 5-7 2011_09_26_drive_0002_sync 分段位移誤差比較表 ... 68 表 5-8 2011_09_26_drive_0002_sync 分段旋轉誤差比較表 ... 69 表 5-9 2011_09_26_drive_0091_sync 分段位移誤差比較表 ... 71 表 5-10 2011_09_26_drive_0091_sync 分段旋轉誤差比較表 ... 72 表 5-11 2011_09_26_drive_0095_sync 分段位移誤差比較表 ... 74 表 5-12 2011_09_26_drive_0095_sync 分段旋轉誤差比較表 ... 75 表 5-13 2011_09_26_drive_0002_sync 估測與真實路徑比較 ... 78 表 5-14 2011_09_26_drive_0091_sync 估測與真實路徑比較 ... 80 表 5-15 2011_09_26_drive_0095_sync 估測與真實路徑比較 ... 82xi
符號定義
第三章 影像特徵點 3.1 SIFT 特徵點演算法 σ 高斯尺度 𝐺(x, y, σ) 可變尺度之高斯函數 𝐼(x, y) 影像𝐼在(x, y)位置的像素值 𝐿(x, y, σ) 高斯模糊影像 𝐷(x, y, σ) 高斯差分影像 x𝑑 高斯差分中紀錄的候選點 H 海森矩陣 α 矩陣的特徵值 β 矩陣的特徵值 𝑟𝑑 矩陣的特徵值比值 𝑚(𝑥, 𝑦) 梯度強度 𝜃(𝑥, 𝑦) 特徵方向 3.2 SURF 特徵點演算法 𝐼Σ(𝑥, 𝑦) 積分影像在位置(x, y)的值 G 高斯函數 𝑥̇ 特徵點偵測時影像上的點位置 Lxx, Lyy, Lxy 點𝑥̇在尺度𝜎之高斯二階偏微分於 X 軸與 Y 軸對應的摺積 Dxx, Dyy, Dxy X 軸與 Y 軸方向對應的盒子濾波器遮罩 𝐻𝑎𝑝𝑝𝑟𝑜𝑥 以盒子濾波器遮罩取代後的近似海森矩陣 𝜔𝐵 Lxx, Lyy, Lxy 與 Dxx, Dyy, Dxy 之轉換對應係數xii 𝐻𝑥, 𝐻𝑦 哈爾小波濾波器在 X 軸與 Y 軸的反應值 𝑉𝑖 特徵點 i 之特徵向量 3.3 ORB 特徵點演算法 𝑡 門閥值 Ip 點𝑝之灰階度 I𝑖 FAST 測量點分割測試圓周上的點 𝑢 表示水平的偏移量 𝑣 表示垂直的偏移量 E(u, v) 灰階度變化 W𝑥,y 視窗函數𝑤的位置 I𝑥,y 在(x, y)位置上像素的灰階 λ1 矩陣特徵值 λ2 矩陣特徵值 R Harris 判斷函數 m𝑝𝑞 二階矩 O FAST 特徵點中心 C FAST 特徵重心 (xi, yi) 二元測試點對 P(xi) xi的灰階值 P(yi) yi的灰階值 τ 二元測試函數 𝑓𝑛𝑑 N 個點對所產生的 N-bits 二元字串 3.4 BRISK 特徵點演算法 𝑐𝑖 第𝑖階層
xiii 𝑑𝑖 第𝑖內階層 𝒯 表示階層與內階層的縮放情形 (p𝑖, p𝑗) BRISK 採樣點之間的配對 𝐼(p𝑖, 𝜎𝑖) p𝑖經過高斯濾波後的結果 𝐼(p𝑗, 𝜎𝑗) p𝑗經過高斯濾波後的結果 g(p𝑖, p𝑗) 兩個採樣點之間的梯度變化 A 採樣點配對的集合 ℒ 長距離點對的集合 𝒮 短距離點對的集合 𝜎𝑚𝑖𝑛 長距離點對集合ℒ的閥值 𝜎𝑚𝑎𝑥 短距離點對集合𝒮的閥值 𝒢 特徵總體方向 ℒ𝑁 ℒ集合中的個數 𝒢𝑥 長距離點對集合在𝑥軸方向的梯度和 𝒢𝑦 長距離點對集合在𝑦軸方向的梯度和 ℬ BRISK 二元敘述子 p𝑖𝜃 像素點p𝑖環繞特徵點旋轉 𝜃後的結果 3.5 A-KAZE 特徵點演算法 𝐝𝐢𝐯 散度 𝛁 梯度 𝑳 影像亮度 𝐜(𝐱, 𝐲, 𝐭) 傳導函數 𝛁𝑳𝝈 𝑳𝝈經過高斯函數後的梯度 𝒕 尺度參數
xiv 𝐠𝟏, 𝐠𝟐 典型擴散係數 𝛕𝒋 可變的步階值 𝛕𝒎𝒂𝒙 最大的步階值 𝛉𝒏 循環反應時間 𝐴(𝐿𝑖) 擴散係數相關矩陣 𝐿𝑖 通過 FED 濾波的圖像 𝛰 A-KAZE 層 𝑆 A-KAZE 子層 第四章 單眼視覺移動估測 4.1 前言 𝑓 相機焦距 𝑝̇(𝑢̇ , 𝑣̇) 平面影像上一點 X 三維座標點 x 二維座標點 𝐊 相機矩陣包含相機內部參數 𝑓𝑢 X軸的等效焦距 𝑓𝑣 Y軸的等效焦距 (𝑢𝑐, 𝑣𝑐) 影像平面中心的座標位置 𝐑𝒄 旋轉矩陣 t𝑐 位移向量 𝑝̃(𝑢̃, 𝑣̃) 點𝑝̇(𝑢̇, 𝑣̇)投影於正規化影像平面後的點 𝑝̆(𝑢̆, 𝑣̆) 稜鏡形變後的點 𝑝̈(𝑢̈, 𝑣̈) 稜鏡形變修正後的點 p ~ 極幾何平面對應𝑝̃的點
xv 4.2 特徵匹配 𝐹, 𝐹′ 兩張連續影像的特徵集合 𝜒 𝐹集合中的特徵點 𝜒′ 𝐹′集合中的特徵點 4.3 單眼視覺測程 Γ EPnP 中的參考點[𝑐1𝑐2𝑐3𝑐4] ∅RPE 重投影誤差 𝜋𝐊 投影函數 (𝐅x)𝑖 表示為𝐅x的第𝑖個元素 [𝐭𝒎]× (t𝑥, t𝑦, t𝑧)T的斜對稱矩陣 𝛿 Sampson distance 距離 ∅EPI 極幾何上的對應誤差 α∅ 阻尼參數控制極限制的權重 𝛽∅ 自訂的權重值 𝑁RPE 3D-2D 關係間的對應數量 𝑁EPI 2D-2D 關係間的對應數量 𝑘, 𝑘′ 自由變動參數 a, a′ 反向投影射線方向向量 c′ 新的相機中心 𝐦𝑖,𝑘 第𝑘張影像中第𝑖個特徵點的觀測狀態向量 𝜔𝑖,𝑘 第𝑘張影像中第𝑖個特徵點的權重 𝑓𝑘−1,𝑘 第𝑘 − 1張影像到第𝑘張影像狀態轉換函數
1
第一章 緒論
1.1 前言
電腦視覺(Computer Vision)以及三維技術在近年來隨著硬體設備的進步以及 新技術的提出，朝向熱門及高度發展的趨勢。三維技術的運用不僅限於電腦科學 的領域，更擴展到工業、考古、建築、醫療及娛樂等多面向。而如何讓電腦視覺 系統了解到系統本身所在的位置及移動的過程，則需要由本體移動估測來得知。 在解決本體移動問題時，除了考量到本體移動估測結果的精準度之外，另外 在解決問題時所花費的計算成本，亦是值得考慮的重要議題。其中特徵點的種類 是影響本體移動估測結果的精確度以及計算成本的因素，因此如何選擇適當的特 徵點演算法，儼然成為值得研究的問題。1.2 研究目的
近年來，除了許多擁有三維技術的硬體裝置蓬勃發展，應用於各種層面之外， 軟體層面上亦有相當大的演進。本體移動估測方法即為許多學界及業界研究的主 題之一，在本體移動估測流程中，特徵點的偵測、描述及匹配皆為重要的一環。 偵測及描述特徵點的演算法自 Moravec 於 1980 年提出角點偵測方法[1]，再 經 1988 年 Harris[2]的改良更加穩定。而 David Lowe 於 1999 年提出尺度不便特 徵轉換演算法 (Scale Invariant Feature Transform, SIFT)[3]，並於 2004 年進行完善 的總結後[4]，Herbert Bay 等人於 2006 年基於 SIFT 演算法提出加速穩健特徵演 算法(Speed Up Robust Features, SURF)[5]。特徵點演算法在近年來有結合不同方法的展現，在 2011 年 Rublee 等學者提 出 ORB(Oriented FAST and Rotated BRIEF)演算法結合併改良兩種既有方法[6]， 同年 Leutenegger 等學者亦提出 BRISK(Binary Robust Invariant Scalable Keypoints) 演算法[7]。Alcantarilla 在 2012 年提出 KAZE 演算法[8]，翌年對該方法改良後提
2 出 A-KAZE 演算法[9]。 本研究之目的，在於以單眼視覺建立一個視覺測程系統，透過特徵偵測、視 覺測程估測及最佳化調整後，對於不同型態的特徵點所產生之結果，進行比較與 分析。首先，針對不同型態特徵點，比較偵測所需的時間及空間等數據，再對所 估測出之結果與真實路徑之準確率進行比較。
1.3 研究方法
為了比較使用不同演算法產生的特徵點，在視覺測程上的效果，首先要先行 建立一個視覺測程的系統，並透過設計的步驟降低系統影響結果的程度。相同情 況下，為了確保移動估測的準確性，故選用 KITTI 標準資料庫之影像資料進行 實驗的測試資料。 要建立視覺測程的系統並比較不同特徵點之結果，首先要先對影像進行特徵 點偵測，在此本研究選用了 SIFT、SURF、ORB、BRISK 及 A-KAZE 等演算法， 對影像序列進行偵測後即可得到特徵點序列。得到特徵點後，本研究對前後影像 進行配對，並透過相關的限制條件及方法增強配對效果。將匹配成功的特徵點坐 持續性的追蹤，再經由三角測量得到三維資訊，並記錄位置及角度狀態變化後， 可將每個位置資訊串連得到估測路徑。 在得到估測路徑後，即可與標準資料庫中所提供的實際路徑進行比較，推算 估測的準確程度。在實驗及分析中，將對路徑估測的結果進行多方面分析，並比 較不同特徵點偵測時所需要的成本。1.4 論文架構
本論文第一章為緒論，包含前言、研究目的及研究方法，其餘章節內容如下 敘述。3 行相關探討。 第三章為說明單眼視覺測程的實作流程。 第四章為影像特徵點，其中會針對 SIFT、SURF、ORB、BRISK 及 A-KAZE 等特徵點演算法，說明演算法相關的理論，包含特徵點的偵測、特徵點描述。 第五章為實驗及結果分析，對於第三章所提出的方法及第四章的相關演算法 理論，進行多組實驗後依據其結果進行分析。 第六章為結論與未來方向。
4
第二章 文獻探討
本研究之目的，在於以單眼視覺建立一個視覺測程系統，並對於不同型態的 特徵點所產生之結果，進行比較與分析。視覺測程(Visual Odometry, VO)的雛形 於 1980 年代由 Moravec[1]開始進行相關的研究，而 VO 的名稱是 Nister 在 2004 年所提出[10]，是以單相機或多相機作為唯一輸入來源，用來解決本體移動的問 題。解決系統自身在環境中的移動路徑問題，稱之為本體移動。 隨著視覺測程相關技術的逐年發展，逐漸可以分為以外觀為基礎或以特徵點 為基礎的兩種模式[11]。在以外觀為基礎的模式中，透過每一個像素的強度建立 密集的影像對應關係(例如[12], [13], [14])，相對的會造成大幅運算的需求。而在 以特徵點為基礎的模式中，採用影像抽像化的概念，透過影像中強度值建立稀疏 但是具有鑑別度的關鍵點並找出關鍵點的描述向量(例如[15], [16], [17], [18])。在 特徵空間中建立影像間像素點的對應關係，可依序用來解決本體移動的問題。 以特徵為基礎的視覺測程流程主要可以表示為圖 2-1。首先對於輸入的影像 序列進行特徵點的偵測與匹配，透過匹配成對之特徵點可推算出特徵點之深度資 訊，並將特徵點資訊記錄，用於後續影像的追蹤。利用所追蹤特徵點，在影像序 列中的移動變化及相機參數可估測出移動路徑，並最後進行最佳化處理使估測結 果更為完善。
5 圖 2-1 基礎視覺測程流程概略圖 特徵偵測(Feature Detection)步驟中，影像特徵辨識演算法由兩個階段組成， 先對影像做關鍵點偵測後，針對偵測到的關鍵點進行描述提取。在關鍵點偵測的 階段，透過使用單尺度或多尺度的變化對影像中以局部進行運算，若在運算的區 域中有某個像素點具有強烈的反應，一般認為這種情況的像素點是包含豐富的結 構資訊。對於在局部區域中心且具有最大程度的反應時，可視該點為關鍵點，為 了 進 一 步 增 加 影 像 特 徵 的 獨 特 性 ， 可 採 用 非 最 大 值 抑 制 (non-maximum suppression)方法來踢除過多的冗餘關鍵點。取得影像中的關鍵點後，下一步要對 關鍵點的局部特性或描述進行編碼，編碼後所得到的向量表示即為關鍵點描述提 取之結果。 「選擇何種影像特徵是最好的?」，在進行本體移動時會思考到這樣的問題 [19]。SIFT 演算法[3], [4]及 SURF 演算法[5]是熱門的選擇，多用於解決姿態估算 (pose-estimation )問題。ORB 演算法[6]近年被提出後，因其大幅度降低計算的需 求，所以多用於嵌入式機器人系統以及即時應用軟體(例如[20])。而 KAZE 演算 法[8]和 A-KAZE 演算法[9]皆為相當新穎的方法，並被認為具有更好的計算效 率。
6 SIFT、SURF 及 ORB 等特徵辨識演算法在一般的情形下，已經由[21]進行比 較與分析後，得到相關的評價與結論。在本研究中，是以基於單眼視覺測程架構 下，所產生的實驗結果，對應 SIFT、SURF、ORB、BRISK 及 A-KAZE 等演算 法，進行後續的比較與分析，根據實驗分析的結果來回應上述的問題。 移動估測(Motion Estimation)是視覺測程中重要的一環，使用相機影像進行 移動估測時，可依據特徵點資訊之類型分為以下三種方式： 2D-2D：二維影像特徵點在不同時間下的對應關係 3D-2D：三維結構與二維影像特徵點的對應關係 3D-3D：三維資料間的對齊 在本研究中，是建立於單眼視覺，因此以第一種方式(2D-2D)作為基礎估測方法。 透過單眼視覺所計算出之估測路徑缺乏實際長度單位，需要以影像序列中最初兩 張做為相對參考基準。
7
第三章 影像特徵點
本章將介紹在本研究中使用到的五種特徵點演算法：SIFT、SURF、ORB、 BRISK 以及 A-KAZE。相關比較如下表：
表 3-1 特徵點演算法比較
SIFT SURF ORB BRISK A-KAZE Origin:
Year and publication 1999 2006 2011 2011 2013 Scale invariant YES YES YES YES YES Rotation invariant YES YES YES YES YES Keypoint measure DoG Hessian FAST FAST Hessian Descriptor type Integer
vector Real vector
Binary string Binary string Binary string Descriptor length 128 bytes 256 bytes 256 bits 512 bits 64/256/486
bits
3.1 SIFT 特徵點演算法
尺度不變特徵轉換演算法(Scale Invariant Feature Transform, SIFT)為電腦視 覺中用來偵測影像中的特徵點以及描述特徵點的演算法，該演算法在 1999 年被 David Lowe[3]提出，並於 2004 年進行完善總結[4]。SIFT 演算法具有抗尺度， 旋轉以及仿射(Affine)等特性，並為後續許多演算法改良之基礎。
8 並透過方向定位描述特徵方向以及定義出特徵點之敘述。 3.1.1 尺度空間極值偵測 為了使影像在不同的縮放尺度下都能夠找出特徵點，SIFT 演算法盡可能在 不同尺度下的影像進行特徵點偵測，藉由該方式呈現抗尺度變化的特性。在許多 合理的假設情況下，Koenderink[22]及 Lindeberg[23]兩位學者在研究中分別提出， 將影像通過高斯模糊函數後最適合呈現影像尺度空間的變化結果。因此在 SIFT 演 算 法 中 建 立 影 像 金 字 塔 (Image Pyramid) 以 及 藉 由 高 斯 差 分 (Difference of Gaussian, DoG)來找出影像在每個尺度空間中之極值。 定義σ為高斯尺度，將可變尺度之高斯函數𝐺(x, y, σ)與影像𝐼(x, y)進行卷積 (Convolution)可得到高斯模糊影像𝐿(x, y, σ)，如下式： 𝐺(x, y, σ) = 1 2𝜋𝜎2𝑒−(𝑥 2+𝑦2)/2𝜎2 (3.1) 𝐿(x, y, σ) = 𝐺(x, y, σ) ∗ 𝐼(x, y) (3.2) 其中∗為卷積。 為了加快尋找極值的速度以及節省計算的成本，SIFT 演算法將鄰近兩個高 斯模糊影像相減，可得到高斯差分影像𝐷(x, y, σ)，如圖 3-1[4]。若定義相鄰影像 的高斯模糊尺度比值為𝑘，則高斯差分影像𝐷(x, y, σ)可由下式推得： 𝐷(x, y, σ) = ((x, y, 𝑘σ) − 𝐺(x, y, σ)) ∗ 𝐼(x, y) (3.3) = 𝐿(x, y, 𝑘σ)- 𝐿(x, y, σ) (3.4)
9 圖 3-1 高斯尺度卷積與差分影像[4] 在 SIFT 演算法中所建立的影像金字塔是依照圖像的大小來決定建立的階數 (Octave)，每一階是相同尺度但模糊程度不同的高斯模糊影像集合，相鄰影像模 糊關係為上述𝑘。而下一階則為兩倍的模糊程度，因此取樣之影像大小為原先之 一半。將同一階內之相鄰高斯模糊影像相減即可得到高斯差分影像𝐷(x, y, σ)，因 此透過每一階所產生之多組高斯影像差分即可進行影像極值之偵測。 在同一階中，取一點像素比較其鄰近 26 點，其中包含該像素所在影像周圍 8 點以及相鄰上下影像各 9 點，若該點像素值為區域內最大值或最小值，則可視 為極值如圖 3-2[34]。將偵測出的極值紀錄即為候選之特徵點，並進行特徵點定 位。
10 圖 3-2 特徵點定位示意圖[34] 3.1.2 特徵點定位 透過尺度空間極值偵測所得到的候選特徵點，需要在進一步進行篩選後才可 被認定為特徵點。使用高斯差分對於影像中的邊緣點以及雜訊較為敏感，因此需 要將這些不可靠的點刪除，保留高於給定閥值且非邊緣的點。 為了求得更準確的位置和尺度並達到次像素標準，同時去除低對比度，要透 過一個三維的二次方程式來進行，首先將高斯差分𝐷(x, y, σ)做泰勒展開式，結果 如(3.5)式： 𝐷(x𝑑) = 𝐷 +𝜕𝐷𝑇 𝜕x𝑑 x𝑑 + 1 2x𝑑𝑇 𝜕2𝐷 𝜕x𝑑2x𝑑 (3.5) 其中x𝑑為高斯差分中紀錄的候選點，而對x𝑑偏微分即為候選點之偏移量，如(3.6) 式： x̂ = −𝜕x𝜕2𝐷 𝑑2 −1 𝜕𝐷 𝜕x𝑑 (3.6)
11 將x̂帶入(3.5)式中可得到(3.7)式： 𝐷(x̂) = 𝐷 +12𝜕𝐷𝜕x𝑇 𝑑 x̂ (3.7) 若|𝐷(x̂)|小於所設定之閥值，則該點被視為對比度過低並將其剔除。 將對比度過低的點剔除後，對於疑似邊緣點的候選點，需要進行檢測並淘汰。 對每一個候選點建立一個2 × 2的海森矩陣(Hession Matrix) H，如(3.8)式： H = || 𝜕2𝐷 𝜕𝑥2 𝜕2𝐷 𝜕𝑦𝜕𝑥 𝜕2𝐷 𝜕𝑦𝜕𝑥 𝜕2𝐷 𝜕𝑦2 || (3.8) 可從建立的海森矩陣中解出矩陣的特徵值(Eigenvalue)α、β (α > β)，以及得到 (3.9)式與(3.10)式如下： 𝑇𝑟(H) = α + β = 𝜕2𝐷 𝜕𝑥2 + 𝜕2𝐷 𝜕𝑦2 (3.9) 𝐷𝑒𝑡(H) = αβ = 𝜕2𝐷 𝜕𝑥2 𝜕2𝐷 𝜕𝑦2 − 𝜕2𝐷 𝜕𝑥𝑦 𝜕2𝐷 𝜕𝑦𝑥 (3.10) 並定義一不等式來偵測疑似邊緣點，如(3.11)式： 𝑇𝑟(H) 𝐷𝑒𝑡(H)< (𝑟𝑑+ 1)2 𝑟𝑑 (3.11) 其中，𝑟𝑑為α 及 β的比值，作為不等式的閥值。若(3.11)式之不等式成立無法成立， 則將該點視為邊緣點並淘汰。 對於通過對比及邊緣測試的候選點，即可視為所求之特徵點。
12 3.1.3 特徵點定向 確認特徵點後，需要進行特徵點描述，如此特徵點才可以進行特徵匹配。SIFT 演算法如前述已具有抗尺度變化之特性，為了進一步能夠達到抗旋轉之效果，因 此需要得到特徵點的方向資訊以及計算出特徵點主方向。將每個特徵點的梯度強 度及方向予以計算，梯度強度𝑚(𝑥, 𝑦)及方向𝜃(𝑥, 𝑦)的計算方式如(3.12)式與(3.13) 式表示： 𝑚(𝑥, 𝑦) = √(𝐿(𝑥 + 1, 𝑦) − 𝐿(𝑥 − 1, 𝑦))2+ (𝐿(𝑥, 𝑦 + 1) − 𝐿(𝑥, 𝑦 − 1))2 (3.12) 𝜃(𝑥, 𝑦) = tan−1(𝐿(𝑥 + 1, 𝑦) − 𝐿(𝑥 − 1, 𝑦)/(𝐿(𝑥, 𝑦 + 1) − 𝐿(𝑥, 𝑦 − 1)) (3.13) 設計一個方向直方圖，將(3.12)式以及(3.13)式所得之梯度強度與方向進行統計。 對於二維的影像來說角度有 360 度，因此取每 10 度為一區間共分為 36 區間，而 對於每一個方向直方圖的區間會加入相對應的高斯函數作為區間內的權重值。最 後統計方向直方圖中每個區間的總和，找出具有最大強度值的區間，區間角度即 為該特徵點的主方向。 若在方向直方圖中，具有最大強度之 80%以上的角度會被視為一個尖峰 (Peak)。若尖峰的數量在同一個直方圖中多餘兩個，則表示該特徵點為多方向 (Multiple Orientation)之情況。針對多方向的特徵點，將會在相同的尺度及位置上 創造一個特徵點，但主方向不同。雖然多方向的特徵點情況只出現在 15%的特徵 點中，但這樣的情況會使比對更為穩定。 3.1.4 特徵點描述子 經由先前章節的描述方法後，特徵點已經具有位置、尺度及主方向等屬性， 在此將對於特徵點進行描述並產生描述子，作為後續進行特徵點比對的依據。
13 特徵點描述的第一步為取一特徵點做為中心，以中心開啟一個16 × 16大小 的視窗，再以4 × 4為單位切割成 16 個子區塊，每一子區塊加入高斯函數作為權 重。統計每一個子區塊內包含像素的梯度強度以及方向後，將結果歸納於具有 8 個方向區間的方向直方圖。
3.2 SURF 特徵點演算法
為了加速尺度不便特徵轉換演算法 (Scale-invariant feature transform, SIFT)的執 行速度及減少 SIFT 演算法產生特徵點敘述子的向量維度，Herbert Bay[5]等人於 2006 年提出加速穩健特徵演算法(Speed Up Robust Features, SURF)，該演算法更 加穩健及快速，為現今進行特徵點偵測已及特徵點描述的代表性演算法之一。
SURF 演算法透過積分影像(Integral Image)的方式進行加速運算，並利用可 變尺度的海森矩陣檢測(Hessian Detector)做快速偵測特徵後，再以哈爾(Haar)濾 波器找出特徵點主方向，最後將各子區域旋轉至同一方向對齊後便可產生出 SURF 特徵點敘述子。 3.2.1 積分影像 在一般情況下進行矩形遮罩的摺積(Convolution)運算時，若矩形遮罩越大需 運算的累加次數就會愈多，且若有相近的矩形遮罩重疊時，重疊得部分無法重複 利用，造成計算浪費的現象。 積分影像利用同一張影像積分只需計算一次的特性進行加速，解決重複計算 的問題。因此矩形遮罩的摺積運算即可直接使用區塊內的累加結果，達到大幅降 低計算量的效果。積分影像概念如圖 3-3 所示，計算公式如下(3.14)
14 圖 3-3 積分影像概念圖 𝐼∑(𝑥, 𝑦) = ∑ ∑ 𝐼(𝑖, 𝑗) 𝑗≤𝑦 𝑗=0 𝑖≤𝑥 𝑖=0 (3.14) 建置完成的積分影像，可以迅速透過查表的方式，找出任一區域內的像素總 和。如圖 3-4 所示，若要找出區域 D 的像素總和，可透過積分影像中該區域的鄰 近四個值𝐼∑1、𝐼∑2、𝐼∑3、𝐼∑4得出。因此 D 區域的像素總合為： 𝐼∑4− 𝐼∑3 − 𝐼∑2+ 𝐼∑1 (3.15) 圖 3-4 透過積分影像計算區域像素總和
15
3.2.2 快速海森檢測
在 SURF 演算法中使用海森矩陣(Hessian Matrix)來尋找極值。藉由矩陣行列 式值可以知道在連續函數中一階倒數所得的臨界點是否為區域間的極值。若給定 一個二階可微分的連續函數𝑓(𝑥, 𝑦)，則海森矩陣可表示為(3.16)： H = || 𝜕2𝑓 𝜕𝑥2 𝜕2𝑓 𝜕𝑦𝜕𝑥 𝜕2𝑓 𝜕𝑦𝜕𝑥 𝜕2𝑓 𝜕𝑦2 || (3.16) 行列式值𝑑𝑒𝑡(H)為(3.17)： 𝑑𝑒𝑡(H) =𝜕2𝑓 𝜕𝑥2 𝜕2𝐷𝑓 𝜕𝑦2 − 𝜕2𝑓 𝜕𝑥𝑦 𝜕2𝑓 𝜕𝑦𝑥 (3.17) 其中判斷是否為極值之分類如下(3.18)： { 𝑑𝑒𝑡(H) > 0 且𝜕2𝑓 𝜕𝑥2 > 0 則𝑓(𝑥, 𝑦)是區域極小點 𝑑𝑒𝑡(H) > 0 且𝜕2𝑓 𝜕𝑥2 < 0 則𝑓(𝑥, 𝑦)是區域極大點 𝑑𝑒𝑡(H) < 1 則𝑓(𝑥, 𝑦)是鞍點 (3.18) 依據上式(3.18)，若為鞍點表示函數在該點時，對 x 軸與 y 軸之一階導函數皆為 0，且兩方向之二階導數其正負號相反。而相對極小值以及相對極大值表示還數 在該點時一階導函數為 0 且兩方向之二階導數其正負號相同。有前述可知，凡是 區域內的極值，其海森矩陣行列式值必定大於 0。透過這樣的特性，進一步推展 至不同尺度的海森矩陣上：
16 𝐻(𝑥, 𝜎) = [𝐿𝐿𝑋𝑋(𝑥, 𝜎) 𝐿𝑋𝑌(𝑥, 𝜎) 𝑋𝑌(𝑥, 𝜎) 𝐿𝑌𝑌(𝑥, 𝜎)] (3.19) 在𝐻(𝑥, 𝜎)中的每一個元素是影像上點 x 在尺度𝜎下的高斯二階偏微分卷積： 𝐿𝑖𝑗(𝑥, 𝜎) = 𝜕2 𝜕𝑖𝜕𝑗𝐺(𝜎) ∗ 𝐼(𝑥, 𝑦), 𝑖, 𝑗 ∈ {𝑥, 𝑦}, (3.20) 其中G(𝜎)為高斯函數。
如同在 SIFT 演算法中為了加速運算，SURF 使用盒濾波器(Box Filter)將海森 矩陣進行簡化，如圖 3-5[34]所示： 圖 3-5 盒濾波器簡化圖[34] 將𝐷𝑥𝑥、𝐷𝑦𝑦及𝐷𝑥𝑦取代海森矩陣後，其行列式可近似表示為(3.21)： det(𝐻𝑎𝑝𝑝𝑟𝑜𝑥) = 𝐷𝑥𝑥𝐷𝑦𝑦− (ω𝐵𝐷𝑥𝑦)2 (3.21) 其中ω為平衡方程式的權重係數，若以大小為9 × 9的盒濾波來偵測，其結果將近 似於使用高斯二階導數取𝜎 = 1.2，此時ω = 0.9如(3.22)。
17 ω =|𝐿𝑥𝑦(1.2)|𝐹|𝐷𝑦𝑦(9)|𝐹 |𝐿𝑦𝑦(1.2)|𝐹|𝐷𝑥𝑦(9)|𝐹 = 0.912 ≈ 0.9 (3.22) SURF 演算法承自於 SIFT 演算法，因此同樣需要建立尺度空間。為了提高 影像取樣的計算效率，SURF 的影像尺度空間會隨著盒濾波器尺度而增加，有別 於 SIFT 演算法中原始影像透過影像金字塔式的縮放，較為不易失真，如圖 3-6[34]。 圖 3-6 SURF 尺度空間示意圖[34] 在尺度空間中同樣會被分割成不同的階數(Octave)，每一階會依尺度變化或是盒 濾波器尺度變化，細分為相同數目的層，並產生不同尺寸的濾波器。但隨著階數 的提高所偵測到的特徵點卻下降。為了降低曾與層之間的尺度變化，SURF 演算 法使用線性內插法，避免取樣因尺度變化過大而較為粗糙。 SURF 特徵點的影像尺度空間定位方法，將每一尺度空間下，取一點比較鄰 近 26 點的海森矩陣行列式值，若為最大即可判定為特徵點。相鄰 26 點為尺度空 間中之上下層各 9 點與所在層之周遭 8 點，如圖 3-7[34]所示。
18 圖 3-7 特徵點定位示意圖[34] 3.2.3 SURF 敘述子 SURF 演算法為了要具有旋轉不變的特性，因此需要將特徵點訂出方向。 SURF 採用哈爾(Haar)小波濾波器來進行水平與鉛直方向的反應值計算，將計算 的反應值來訂定特徵方向，如圖 3-8。 圖 3-8 哈爾(Haar)小波濾波器 以特徵點為圓心，取特徵點所在的尺度 6 倍長為半徑，對此區域內以圓心角π/3 的扇形區塊進行掃描一圈，統計同一掃描角度內之特徵點對於哈爾小波的反應值， 擁有最大的加總值的角度即為特徵點的主要方向。 在得到特徵點定位以及主方向後，要對特徵點建立特徵向量。以特徵點為中 心，設立一個邊常為所在尺度 20 倍的正方形區域，其方向為特徵點之方向。將 該正方形分割為大小4 × 4的子方型區域，每個子方型區域中，取 5× 5的取樣點， 並將哈爾小波濾波器的邊長為尺度的 2 倍。定義𝐻𝑥為哈爾小波濾波器的水平反
19 應值，而𝐻𝑦為哈爾小波濾波器的垂直反應值，水平與垂直方向會依照所選擇之 特徵點方向而定，再將𝐻𝑥及𝐻𝑦分別加總為 ∑𝐻𝑥 及 ∑𝐻𝑦 。而為了得到特徵強度變 化的資訊因此取𝐻𝑥及𝐻𝑦之絕對值並分別加總後為 ∑|𝐻𝑥| 及 ∑|𝐻𝑦| ，最終得到 SURF 特徵點的特徵向量𝑉𝑖，向量為度為4 × 4 × 4 = 64。 𝑉𝑖 = (∑𝐻𝑥, ∑𝐻𝑦, ∑|𝐻𝑥|, ∑|𝐻𝑦|) (3.23)
3.3 ORB 特徵點演算法
ORB (Oriented FAST and Rotated BRIEF)是 Rublee 等學者在 2011 發表於 ICCV (International Conference on Computer Vision)之演算法[6]，ORB 演算法是基 於 FAST (Features from Accelerated Segment Test)[24]及 BRIEF (Binary Robust Independent Elementary Features)[25]之結合與改良，其中改良 FAST 演算法進行 特徵點偵測與特徵方向，BRIEF 演算法提供特徵點描述子。
3.3.1 FAST 特徵偵測
ORB 演算法透過 FAST 演算法進行特徵點偵測，FAST 演算法將影像中每一 個像素皆做為待測量點 p ，並以 𝑝 為圓心圈選半徑為 3 像素之圓，取圓周上涵 蓋之 16 點像素進行分割測試(segment test)，如圖 3-9。
20
圖 3-9 測量點分割測試
將 𝑝 之灰階度(intensity)定義為 Ip ，取一 𝑡 做為門閥值(threshold value)，定義圓 周上每個像素 I𝑖 與待測點 𝑝 之關係，可表示為：
compared𝑖 = {
0, I𝑖 ≤ Ip− 𝑡 (darker)
1, I𝑖 > Ip+ 𝑡 (brighter) (3.24) FAST-9 定義為當圓周上有連續或累積多於 9 個相同關係時可將𝑝認定為候選特徵 點，而在 FAST-12 中可透過高速測試(high-speed test)來加快判定是否為特徵點。
在 FAST-12 中先判斷圖中 1, 9 像素是否為相同關係，若相同再判斷 5, 13 像 素，至少需要有三個相同關係才會判斷剩下的 12 個像素，否則直接判定為非候 選特徵點。 3.3.2 Harris 角點偵測 待測點 𝑝 經過 FAST 特徵偵測後雖然可以初步判斷是否為候選特徵點，所得 之候選特徵點可能有三種型態如圖 3-10，但無法確認該候選特徵點為何種，因此 需透過 Harris 角點偵測來進行判別。
21 圖 3-10 特徵點三種型態示意圖 Harris 角點偵測定義一個視窗𝑤對影像進行垂直與水平方鄉的移動偵測，透 過灰階變化的程度可找出欲求之角點，灰階變化(E(u, v))之定義如公式(3.25)。 E(u, v) = ∑ W(x, y)(𝐼(𝑥 + 𝑢, 𝑦 + 𝑣) − 𝐼(𝑥, 𝑦))2 x,y (3.25) 其中 W𝑥,y 表示視窗函數𝑤的位置 I𝑥,y 表示在(x, y)上像素的灰階度 𝑢 表示水平的偏移量 𝑣 表示垂直的偏移量 取Ix+u,y+v 透過泰勒展開式並捨去二階以上項，可得到： Ix+u,y+v ≈ I𝑥,y+ uIx+ vIy (3.26) 因此
E(u, v) ≈ ∑ W𝑥,y(I(x, y) + uIx(x, y) + vIy(𝑥, 𝑦) − I(x, y))2 x,y
22
E(u, v) ≈ ∑ W𝑥,y(uIx(x, y) + vIy(x, y)) 2 x,y (3.28) E(u, v) ≈ ∑ W𝑥,y(u2I x2(x, y) + v2Iy2(x, y) + uIx(x, y)vIy(x, y) x,y + vIy(x, y)uIx(x, y)) (3.29) E(u, v) ≈ ∑ 𝑊𝑥,𝑦([𝑢 𝑣] [ 𝐼𝑥 2(𝑥, 𝑦) 𝐼 𝑥(𝑥, 𝑦)𝐼𝑦(𝑥, 𝑦) 𝐼𝑦(𝑥, 𝑦)𝐼𝑥(𝑥, 𝑦) 𝐼𝑦2(𝑥, 𝑦) ] [𝑢𝑣]) 𝑥,𝑦 (3.30)
將 W𝑥,y 代入高斯模糊函數(Gaussian Blur)G
E𝑢,𝑣 ≈ [𝑢 𝑣] [Ix 2⊗ G I xIy⊗ G IyIx⊗ G Iy2 ⊗ G ] [𝑢𝑣] (3.31) 其中可定義一矩陣 M 𝑀 = [Ix 2⊗ G I xIy ⊗ G IyIx⊗ G Iy2 ⊗ G] (3.32) M 矩陣可得其特徵值(eigenvalue) λ1 及 λ2 ，表示水平及垂直方向之強度分布情形， 若值越大表示強度越強。 λ1 及 λ2 之關係所產生之對應結果如圖 3-11。
23 圖 3-11 依據λ1及λ2關係所產生之情形 如圖，可將 λ1 及 λ2 之關係歸類出三種情況： Flat： λ1 及 λ2 皆很小，表示其水平與垂直方像變化強度不明顯 Edge： λ1 及 λ2 其中一個值較大，表示在水平或垂直方像變化強度較顯著 Corner： λ1 及 λ2 兩者皆大，表示在水平及垂直方向變化強度明顯 為了加速判斷，Harris 定義一判斷函數 R R = det(M) − kTr(M)2 (3.33) det(M) = λ1λ2, Tr(M) = λ1+λ2, 0.04 ≤ k ≤ 0.06 (3.34) { Edge R < 0 Flat R = 0 Corner R > 0 (3.35) 3.3.3 Oriented FAST 為了解決 FAST 的演算法無法處理多尺度變化的影像，ORB 在進行特徵偵 測時先將影像建立一影像金字塔空間，該空間被分割為數階層，每一層代表不同 的尺度變化，並將每一層影像分別進行 FAST 特徵點偵測，提供尺度不變特徵。
24
然而 FAST 因缺乏特徵方向，因此在 ORB 中透過二階矩(moment)找出重心 C 並利用與特徵點中心 O 的向量來計算出特徵方向，二階矩 m𝑝𝑞計算可由公式 (3.36)： m𝑝𝑞 = ∑ 𝑥𝑝𝑦𝑞 𝑥,𝑦 𝐼(𝑥, 𝑦) (3.36) 重心 C 之計算為： C = ( m10 m00 m01 m00 ) (3.37) 依據上述可得一向量𝑂𝐶⃗⃗⃗⃗⃗ ，方向角度可表示為： θ = atan2(m10, m01) (3.38) 3.3.4 BRIEF 敘述子
在取得特徵點後，利用 BRIEF (Binary Robust Independent Elementary Features) 演算法進行特徵點描述，然而為了讓 BRIEF 演算法可以達到旋轉不變並降低雜 訊的影響，在進行 BRIEF 前先將所有特徵方向旋轉至統一方向確保方向一致： S = (xyi, … , xn i, … , yn) (3.41) Sθ = RθS (3.42) gn(P, θ) ∶= 𝑓𝑛(P)|(xi, yi) ∈ Sθ (3.43)
25 BRIEF 演算法選定一個大小長寬固定為 S 像素的區塊，以特徵點為區塊(P) 中心，對區塊內隨機抽取 N 個點對(xi, yi)，如圖 3-12 為對長寬為 5 像素的區塊 中隨機取樣 7 個點對。 圖 3-12 隨機取樣 7 個點對示意圖 將區塊中每個點對(xi, yi)中， xi 的灰階值 P(xi) 與 yi 的灰階值 P(yi) 進行比對二 元測試(τ)，其測試定義如下式： τ(P; xi, yi) ∶= {1, P(x0, P(xi) < P(yi) i) ≥ P(yi) (3.39) 對於 N 點對(xi, yi)的二元測試結果，可到 N 個二元結果，可形成長度為 N 的二 元串列，如圖中對 8 個點對依序進行二元測試，所得到的 0 或 1 的結果，可形成 長度為 8bits 的二元字串，以公式來表示 N 個點對所產生的 N-bits 二元字串。 𝑓𝑛𝑑(P) ∶= ∑ 2𝑖−1τ(P; x i, yi) 1≤𝑖≤𝑛𝑑 (3.40)
26
在進行特徵點匹配時，兩特徵點的二元字串差異即可推算出漢明距離，漢明距離 愈短表示兩特徵點愈相似。
3.4 BRISK 特徵點演算法
BRISK (Binary Robust Invariant Scalable Keypoints)同樣於 2011 年在 ICCV 中 被提出[7]，Leutenegger 等學者以 Mair 等學者所提出之 AGAST (Adaptive and Generic Accelerated Segment Test)演算法[27]為研究基礎，並將 FAST (Features From Accelerated Segment Test)演算法應用於尺度空間中的測量，使得 BRISK 能 在連續尺度空間中估測每個關鍵點的真實比例。 3.4.1 尺度空間測徵點偵測 在 BRISK 演算法的架構中，尺度空間金字塔包含n個階層(octaves) 𝑐𝑖以及n個 內階層(intra-octaves) 𝑑𝑖(𝑖 = {0, 1, … , 𝑛 − 1})，一般情況取𝑛 = 4。每一個階層𝑐𝑖是 透過對前一個階層𝑐𝑖−1取樣而來，取樣倍率為 0.5。每一個內階層𝑑𝑖位於階層𝑐𝑖以 及𝑐𝑖+1之間，對於第一個內階層𝑑0是對階層𝑐0以 1.5 倍率取樣而得到，其餘的內 階層𝑑𝑖是以前一個內階層𝑑𝑖−1以 0.5 倍率取樣而來，如圖 3-13[7]。因此可定義一 個尺度參數𝒯表示階層與內階層的縮放情形： 𝒯(𝑐𝑖) = 2𝑖 (3.44) 𝒯(𝑑𝑖) = 2𝑖⋅ 1.5 (3.45)
27 圖 3-13 BRISK 尺度空間檢測示意圖[7] BRISK 演算法在進行特徵點偵測時，對於階層與內階層都進行 FAST-16 測 試，並須滿足給定的閥值 T，用來找出可能的關鍵點。對於𝑐0需做特殊處理，使 用 FAST5-8 進行測試後得到虛擬層 d(-1)分數做為𝑐0下一層。接下來，對於這些 區 域 的 可 能 關 鍵 點 ， 需 要 在 尺 度 空 間 中 進 行 非 最 大 值 抑 制 (non-maxima suppression)：關鍵點需要在同一層中，滿足相鄰 8 個像素 FAST 分數最大的條件。 且所在層的上層和下層的分數也需要抑制。FAST 分數 s 被定義為考慮圖像關鍵 點仍為角點的最大閾值。 3.4.2 採樣方法與旋轉估計
相較於 SIFT 以及 SURF 演算法採用特徵向量來描述特徵點，BRISK 採用二 元串列，透過簡單的亮度比較測試後來進行描述。如同 BRIEF 演算法，這種二 元測試的方法已證實是非常有效的[25]。在 BISK 中確定了每個特徵點的特徵方 向，實現旋轉不變性是特徵點強健性的關鍵。
28 在 BRISK 中，以特徵點為中心取40 × 40像素大小的採樣區塊。在區塊中建 立環繞特徵點且不同半徑的同心圓，在每一個圓上以相同間格方式進行採樣包含 特徵點可得到 N 個點(N=60)，為了消除採樣點中的強度干擾，因此需要對採樣 點進行高斯濾波。BRISK 採樣模型[7]如下圖所示： 圖 3-14 BRISK 特徵點採樣模型[7] 圖 3-14 中，藍色圈表示採樣中心點進行高斯濾波，高斯濾波尺度σ𝑖與同心圓上 之像素間格正比；紅圈為經過高斯濾波後的採樣點。 對於採樣點之間的配對表示為(p𝑖, p𝑗)，p𝑖及p𝑗經過高斯濾波後的結果分別 表示為𝐼(p𝑖, 𝜎𝑖)及𝐼(p𝑗, 𝜎𝑗)，因此兩個採樣點之間的梯度變化g(p𝑖, p𝑗)可由下式進 行計算： g(p𝑖, p𝑗) = (p𝑖, −p𝑗) ⋅𝐼(p𝑖, 𝜎𝑖) − 𝐼(p𝑗, 𝜎𝑗) ‖p𝑖− p𝑗‖2 (3.46)
29
其中1 ≤ 𝑖, 𝑗 ≤ 𝑁。所有採樣點配對的集合可以表示為 A：
A = {(p𝑖, p𝑗) ∈ 𝐑𝟐× 𝐑𝟐|𝑖 < N ∧ 𝑗 < 𝑖 ∧ 𝑖, 𝑗 ∈ 𝑁} (3.47)
在 BRISK 中定義採樣點之間的距離，分為長距離點對(long-distance pairings) 的集合ℒ以及短距離點對(short-distance pairings)的集合𝒮： ℒ = {(p𝑖, p𝑗) ∈ A|p𝑖 − p𝑗 > 𝜎𝑚𝑖𝑛} ⊆ A (3.48) 𝒮 = {(p𝑖, p𝑗) ∈ A|p𝑖1 − p𝑗 < 𝜎𝑚𝑎𝑥} ⊆ A (3.49) 其中𝜎𝑚𝑖𝑛為長距離點對集合 ℒ 的閥值，在此定義為13.67𝑡；𝜎𝑚𝑎𝑥為短距離點對 集合𝒮的閥值，在此定義為9.75𝑡，𝑡為特徵點所在的尺度空間。 利用長距離點對集合 ℒ 可計算出特徵總體方向 𝒢： 𝒢 = [𝒢𝒢𝑥 𝑦] = 1 ℒ𝑁⋅ ∑ g(p(p 𝑖, p𝑗) 𝑖,p𝑗)∈ℒ (3.50) 其中 ℒ𝑁 為 ℒ 集合中的個數，𝒢𝑥為長距離點對集合在𝑥軸方向的梯度和，𝒢𝑦為長 距離點對集合在𝑦軸方向的梯度和。 3.4.3 BRISK 特徵點描述 為了確保 BRISK 建立的特徵點敘述子具備旋轉不變性，需要計算特徵點主 方向𝜃：
30 𝜃 = arctan2 (𝒢𝑦, 𝒢𝑥) (3.51) 特徵點主方向確定後，將短距離點對集合 𝒮 中每個點對進行灰階強度值比較， 產生二元敘述子 ℬ ： ℬ = {1, 𝐼(p𝑗𝜃, 𝜎𝑗) > 𝐼(p𝑖𝜃, 𝜎𝑖) 0, 𝑒𝑙𝑠𝑒 ∀(p𝑖𝜃, p𝑗𝜃) ∈ 𝒮 (3.52) 其中p𝑖𝜃為像素點p 𝑖環繞特徵點旋轉 𝜃後的結果。利用 3.4.2 節中所述，特過採樣 模式和距離閾值可獲得長度為 512bits 的二元序列。
3.5 A-KAZE 特徵點演算法
在 2012 年時 Alcantarilla 等學者提出 KAZE 特徵演算法[8]，KAZE 為日文單 字「風」的意思。大自然中的空氣流動通常是一種非線性的過程，因此 KAZE 演算法有別於一般 SIFT 演算法或 SURF 演算法等線性高斯金字塔尺度分析，而 是在圖像中進行非線性擴散處理。
A-KAZE 演算法是在 2013 年時由 Alcantarilla 等學者基於 KAZE 演算法進行 改良[9]，可使運算速度提升並減少運算量。A-KAZE 演算法如同 KAZE 演算法 使用非線性尺度空間來進行特徵點偵測以及特徵點描述。
3.5.1 非線性擴散濾波
非線性擴散濾波(Nonlinear Diffusion Filtering)描述在變動尺度中影像亮度的 變化做為某種程度的流動函數之散度，典型的擴散方程式可表示如下：
31 𝝏𝑳 𝝏𝒕 = 𝐝𝐢𝐯(𝐜(𝐱, 𝐲, 𝐭) ∙ 𝛁𝐋) (3.53) 其中𝐝𝐢𝐯代表散度、𝛁表是為梯度以及𝑳為影像亮度，𝒕在方程式中做為尺度參數， 當𝒕值愈大时會造成影像顯示愈簡單。 在同一尺度下，傳導函數定義為𝐜(𝐱, 𝐲, 𝐭)如下： 𝐜(𝐱, 𝐲, 𝐭) = 𝐠(|𝛁𝑳𝝈(𝒙, 𝒚, 𝒕)|) (3.54) 𝛁𝑳𝝈為𝑳𝝈經過高斯函數後的梯度，函數𝐠表示為擴散係數，Perona 及 Malik[26]提 出典型的擴散係數𝐠𝟏、𝐠𝟐分別為： 𝐠𝟏 = 𝐞𝐱𝐩 (− |𝛁𝑳𝝈|𝟐 𝝀𝟐 ) (3.55) 𝐠𝟐 = 𝟏 𝟏 +|𝛁𝑳𝝈|𝟐 𝝀𝟐 (3.56) 函數𝐠𝟏有助於保留對比度較高之邊緣，函數𝐠𝟐則易於保留較寬闊的區域，而參數 𝝀為限制因子來做為影像中邊緣保留的程度控制，𝝀值愈大保留的邊緣細節愈少。 在 A-KAZE 演算法中採用𝐠𝟐函數做為擴散係數。 3.5.2 快速顯式擴散 因為沒有用於非線性擴散方程式的解析解，因此需要透過非線性偏微分方程 (PDE)來進行離散化以求得近似解。而採用顯格式(explicit schemes)是簡單且容易 實現的方法，但因對於穩定性問題所以需要多次的迭代運算以達到期望的規模。 在 KAZE 演算法中以加法運算子分割(Additive Operator Splitting, AOS)來解決，
32
AOS 為半隱格式(semi-implicit schemes)方法，可用於任何步長並可創建非線性尺 度空間以作為特徵點偵測及特徵點描述。但 AOS 需要在每個時間步長(time step) 進行線性方程大型系統求解。
在 A-KAZE 演算法中使用快速顯擴散(Fast Explicit Diffusion, FED)來加快非 線性尺度空間中的特徵偵測。FED 結合顯格式與半隱格式的優點，是基於顯格 式的盒濾波器(box filters)分解。迭代盒濾波器(Iterated box filters)結果近似於高斯 函數，具有好的品質及容易實作。FED 將原本的𝐧次迭代擴散方程分為 𝐌 個新 的循環計算，在每一個循環中可變的步階值(step sizes) 𝛕𝒋 來自於盒濾波器的分 解： 𝛕𝒋= 𝛕𝒎𝒂𝒙 𝟐 𝐜𝐨𝐬𝟐(𝝅𝟐𝒋 + 𝟏 𝟒𝒏 + 𝟐) (3.57) 其中𝛕𝒎𝒂𝒙為最大的步階值，但仍然保有顯格式的穩定條件。對於對應一次 FED 的循環反應時間𝛉𝒏可由下式獲得： 𝛉𝒏 = ∑ 𝛕𝒋 𝒏−𝟏 𝒋=𝟎 = 𝛕𝒎𝒂𝒙 𝒏𝟐+ 𝒏 𝟑 (3.58) 典型的擴散方程式(3.53)離散化，可以使用顯格式透過向量矩陣來表示，如 下： 𝐿𝑖+1+ 𝐿𝑖 τ = 𝛢(𝐿𝑖)𝐿𝑖 (3.59) 其中𝐴(𝐿𝑖)為一矩陣與擴散係數有關，由 Weickert 所提出。而𝐿𝑖+1可直接經由通
33 過 FED 濾波的圖像 𝐿𝑖 及 𝐴(𝐿 𝑖)來計算，如下： 𝐿𝑖+1 = (𝐼 + 𝛕A(𝐿𝑖))𝐿𝑖 (3.60) 𝐼為單位矩陣。若給定𝐿𝑖+1,0 = 𝐿𝑖，則對於𝐧次且步階值為𝛕 𝒋的 FED 循環可得到下 式： 𝐿𝑖+1,𝑗+1 = (𝐼 + 𝛕 𝒋𝐴(𝐿𝑖)) 𝐿𝑖+1,𝑗, 𝑗 = 0, … , 𝑛 − 1 (3.61) 當每一次的 FED 執行完畢時，需要計算下一個新的𝐴(𝐿𝑖)矩陣。 3.5.3 非線性尺度空間
在建立非線性尺度空間(Nonlinear Scale Space)前，需要先定義一組演化時間。 尺度空間由層(octaves) 𝛰 及子層 (sub-levels) 𝑆 組成。對於每一個層及子層可 對應到相對的尺度 𝜎，對於尺度的轉換如下列公式所示： σ𝑖 = 2𝑜+𝑠/𝑆, 𝑜 𝜖 [0 … 𝛰 − 1], 𝑠 𝜖 [0 … 𝑆 − 1], 𝑖𝜖[0 … 𝑀] (3.62) 其中 M 是濾波影像的總數(𝛰 × 𝑆)。 對於非線性擴散濾波是以時間做為擴散的變化，因此需將尺度空間由像素單 位 σ𝑖 轉換至時間單位 𝑡𝑖 ，轉換關係可由下式得知。 𝑡𝑖 = 1 2σ𝑖2, 𝑖 = {0 … 𝑀} (3.63)
34 3.5.4 特徵偵測 在非線性尺度空間中，計算每一個通過 FED 濾波圖像 𝐿𝑖 的海森矩陣行列 式，而對於多尺度微分運算需要依據尺度進行正規化： 𝜎𝑖,𝑛𝑜𝑟𝑚= 𝜎𝑖/2𝑜𝑖 (3.64) 𝐿𝑖𝐻𝑒𝑠𝑠𝑖𝑜𝑛 = 𝜎 𝑖,𝑛𝑜𝑟𝑚2 (𝐿𝑖𝑥𝑥𝐿𝑖𝑦𝑦− 𝐿𝑖𝑥𝑦𝐿𝑖𝑥𝑦) (3.65) 為了計算二階導數，使用 Scharr 濾波器搭配步階值為𝜎𝑖,𝑛𝑜𝑟𝑚，Scharr 濾波器具有 近似旋轉不變性的優點。對於每一層的影像以3 × 3的視窗進行 Scharr 濾波器檢 測，並尋找大於給定閥值的極值點。通過的極值點再判斷是否為𝑖 + 1層及𝑖 − 1層 中𝜎𝑖 × 𝜎𝑖範圍的相對最大值。 為了加速尋找極值，以待測點為中心取鄰近 26 點進行海森矩陣運算，若其 行列式值皆大於鄰近點即可視為特徵點並以次像素來定位。 3.5.5 特徵描述
A-KAZE 在特徵點描述中提出 M-LDB (Modified-Local Difference Binary)做 為改良 LDB (Local Difference Binary)的方法，利用來自非線性尺度空間的梯度變 化和強度資訊。 LDB 描述方法在[28]中有詳細說明，並遵循與 BRIEF 相同的原 理，然而 LDB 是以區域間的平均值進行二元測試增加強健性，與 BRIEF 使用單 個像素不同。對於每一個比較測試，除了強度值之外，亦考量比較區域中的水平 和垂直導數的平均值，形成 3bits 的比較結果。 LDB 提出在更細的步驟中使用各種尺寸的網格，將區塊劃分為2 × 2、3 × 3 或4 × 4等網格。如果敘述子是直立的（非旋轉不變），這些分割區塊的平均值可
35 以使用積分圖像來快速計算，如在[28]中。然而，當考慮到特徵點的旋轉時，不 能使用積分圖像，並對旋轉後分割區塊中所有點在計算時間上可能相對昂貴。藉 由如 KAZE 中計算特徵點的主方向來獲得旋轉不變性後，LDB 依照特徵點主方 向進行網格旋轉。 M-LDB 以特徵尺度函數σ的步長對網格進行二次採樣，取代網格內每個分割 區塊中的所有像素取平均值，這種近似於平均值的結果在實驗中表現良好。尺度 依賴(scale-dependent)的採樣可使得敘述子對尺度變化具有強健性。M-LDB 使用 在特徵檢測步驟中計算的導數，減少建立描述子過程所需要的步驟數，並透過縮 減敘述子的大小，可望改善結果或至少減少計算負荷而不降低性能。
36
第四章 單眼視覺移動估測
本章將介紹單眼視覺測程系統中實作與基本原理，實作中的特徵點來自於 SIFT、SURF、ORB、BRISK 及 A-KAZE 等特徵演算法。單眼視覺測程系統實作 流程如圖 4-1 所示 圖 4-1 本研究所實作之單眼視覺測程流程圖4.1 前言
4.1.1 相機模型 透過相機對三維空間擷取影像時，會因為相機擺放的位置、方向以及內部的 元件屬性不同，而造成轉換平面影像時的差異。對於三維空間與二維平面影像的37
投影關係，可經由相機的幾何參數描述。在相機的幾何參數中，可分內部參數以 及外部參數，內部參數針對相機本身與內部元件的屬性，包含：鏡頭焦距(Focal Length)、稜鏡扭曲(Lens Distortion)、像素長寬比(Pixel Aspect Ratio)與影像主軸 點(Principle Point)；在外部參數的部分，因為相機內部座標系統與真實世界座標 系統不同，兩者兼具有一定的轉換關係，其關係可用一組旋轉矩陣(Rotation Matrix)與位移矩陣(Translation Matrix)來表達。 最常用的透視相機模型可以說是針孔投影系統，透過針孔(Pinhole)成像原理 將三維空間座標投影至平面影像上，如圖 4-2 可表示針孔相機模型[29]。在圖 4-2 中，假設𝑓為相機焦距，在空間中有一點𝑝，三維座標為(𝑥, 𝑦, 𝑧)。點𝑝經過透視投 影(Perspective Projection)後，所投射在平面影像上一點𝑝̇(𝑢̇ , 𝑣̇)，𝑢̇表示平面影像 中之X軸方向座標，而𝑣̇則是Y軸方向的座標。 圖 4-2 相機針孔模型[29] 給定空間中三維座標點X = (𝑥, 𝑦, 𝑧)T以及投影在影像平面之二維座標點x = (𝑢, 𝑣)T，則相對應之投影模型如下： 𝑓 𝑝(𝑥, 𝑦, 𝑧) Image Plane 𝑦𝐼𝑚𝑎𝑔𝑒 𝑥𝐼𝑚𝑎𝑔𝑒 𝑝̇(𝑢̇, 𝑣̇) 𝑧 𝑜 𝑦 𝑥 Optical Axis
38 ( 𝑢 𝑣 1) ∼ ( 𝑓𝑢 0 0 𝑓𝑣 0 0 𝑢𝑐 0 𝑣𝑐 0 1 0 ) ( 𝑥 𝑦 𝑧 1 ) = 𝐊(𝐈 𝟎) ( 𝑥 𝑦 𝑧 1 ) (4.1) 其中∼表示為等價，𝐊為相機矩陣包含相機內部參數，由一個3 × 3的上三角矩陣 表示如下： (𝑓0 𝑓𝑢 𝛾𝑣 0 0 𝑢𝑐 𝑣𝑐 1) (4.2)
在相機矩陣中， 𝑓𝑢 與 𝑓𝑣 表示為X軸與Y軸的等效焦距(Effective Focal Lengths)，𝑢𝑐 及𝑣𝑐代表影像平面中心的座標位置，而 𝛾 為影像平面中X軸與Y軸之夾角與相對直 角的偏移程度，通常可忽略為 0。 若假設X = (𝑥, 𝑦, 𝑧)T是一個靜態的點，即當相機移動到下一個影像位置時， 仍然保持在原本座標上。當相機移動時，可以得到新的投影關係： (𝑢 ′ 𝑣′ 1 ) ∼ 𝐊(𝐑𝒄 t𝑐) ( 𝑥 𝑦 𝑧 1 ) (4.3) 其中旋轉矩陣𝐑𝒄 ∈ SO(3)以及位移向量t𝑐 ∈ ℝ3共同組成歐幾里得轉換(Euclidean transformation)，用來表示相機的運動情形。 透過相機取得平面影像時，容易受到稜鏡的曲率影響而產生折射，造成影像 產生徑向形變(Radial Distortion)，由圖 4-3 可示意正常影像與徑向形變影像。
39
(a)正常影像
(b)
圖 4-3 形變示意圖 (a)正常影像 (b)徑向形變影像
解決影像形變，首先對點𝑝̇(𝑢̇, 𝑣̇)投影於正規化影像平面(Normalized Image Plane) 後，得到點𝑝̃(𝑢̃, 𝑣̃)如圖[34]所示：
40 圖 4-4 具形變之相機投影模型[34] 由(4.4)式中可求得點𝑝̃(𝑢̃, 𝑣̃)： (𝑢̃𝑣̃ 1 ) ∼ (1 00 1 0 0 0 0 0 0 1 0 ) (𝐑𝒄 t𝑐 0 1) ( 𝑥 𝑦 𝑧 1 ) (4.4) 形變後的點𝑝̆(𝑢̆, 𝑣̆)，由(4.5)式中計算： (𝑢̆𝑣̆ 1 ) ∼ (𝑢̃𝑣̃ 𝑟22𝑢̃𝑣̃+ 2𝑣̃2 0 0 𝑟2 + 2𝑢̃2 0 2𝑢̃𝑣̃ 0 0 1 ) ( 1 + 𝑘1𝑟2 + 𝑘 2𝑟4+ 𝑘3𝑟6+ ⋯ 𝑟1 𝑟2 1 ) (4.5) 其中𝑟2 = 𝑢̃2+ 𝑣̃2為徑向形變(Radial Distance)之平方，代表理想點與中心座標之 距離，(𝑘1, 𝑘2, 𝑘3… , 𝑟1, 𝑟2)為形變細數(Distortion Coefficients)。將𝑝̆(𝑢̆, 𝑣̆)投影回到 影像平面上，可得到修正後的點𝑝̈(𝑢̈, 𝑣̈)，轉換如(4.6)式：
Camera Image Plane 𝑥𝑤 𝑦𝑐 𝑧𝑐 𝑥𝑐 (0, 0) (0, 0)
Normalized Image Plane 𝑧𝑤 𝑦𝑤 World Camera 𝑧 = 1 𝑧 = 𝑓 𝑝 𝑝̇ 𝑝̃ 𝑝̆ 𝑝̈
41 (𝑢̈𝑣̈ 1 ) ∼ (𝑓0 𝑓𝑢 𝛾𝑣 0 0 𝑢0 𝑣0 1) ( 𝑢̆ 𝑣̆ 1 ) (4.6) 4.1.2 影像校正與極幾何限制 基於 4.1.1 之相機模形進行相機校正，其目的在於找出相機之內部參數及外 部參數，或得參數用來推測出影像資訊，在本研究中使用的影像資料來自於 KITTI 標準資料庫(The KITTI Vision Benchmark Suite)[37]，該網站所提供之影像 序列已經校正完成，故在本研究中不探討校正實作的過程。 在 2D-2D 的對應情況下， 移動估測的主要特色為極幾何限制 (epipolar constraint)，其概念為在連續圖像中，有一特徵點p̃可在下一張圖像中找到對應的 特徵點p~，可藉由下列公式計算出： 𝑝̃′⊤ 𝐸p̃ = 0 (4.7) p̃以及𝑝̃′為正規化座標，關於極幾何限制如圖 4-5 [10]表示： 圖 4-5 極幾合限制示意圖[10]
42
4.2 特徵匹配
影像序列經由特徵點偵測演算法運算後，可得到特徵點序列，而如何透過特 徵點推測影像的變化，首先需要進行相鄰影像間的特徵點匹配。然而，初步的特 徵點匹配會得到許多不明確或是錯誤的對應關係，因此需要使用多種的方法進行 過濾以及增加強健性。 4.2.1 暴力描述子匹配 在特徵空間中進行初步的匹配，定義𝐹及𝐹′為兩張連續影像的特徵集合。影 像中特徵點𝜒 𝜖 𝐹對應到𝜒′ 𝜖 𝐹′之關係可表示為： 𝜒′ = 𝑎𝑟𝑔𝑚𝑖𝑛 𝜒 ̆𝜖𝐹 𝑑(𝜒, 𝜒∗) (4.8) 其中𝑑：𝐹 × 𝐹′ → ℝ 表示為一個選定的方法來衡量兩個特徵之間的相似性。 對於不同型態的特徵點敘述子，所採用的檢測相似度方法亦有所不同。對 於 二 元 型 態 的 特 徵 敘 述 子 (binary feature descriptors) 使 用 漢 明 距 離 (Hamming distance)做為比對基礎；其它型態的特徵點敘述子可能採用平方距離 之和（Sum of Squared Distances, SSD）或絕對距離之和（Sum of Absolute Ddistances, SAD）進行測量。在本研究中，在𝐹′中採用暴力搜尋(brute-force)方法 找出連續影像中的對應關係𝜒 → 𝜒′。
4.2.2 匹配結果檢測
在初步的特徵點匹配完畢後，需要進一步使用限制條件剔除不明確的匹配組 合。給定𝜒̆′ 𝜖 𝐹′是與𝜒匹配的次佳結果，在𝜒 → 𝜒′的對應關係中如果滿足(4.9)式 時，可認定為不明確的匹配組合，並將之剔除。
43 𝑑(𝜒, 𝜒′) 𝑑(𝜒, 𝜒̆′)< 𝑟 (4.9) 其中𝑟 = (0, 1]為給定的的閥值，Lowe[4]提出定義𝑟 = 0.6。 除了對於比值的限制之外，下一步針對配對不一致的情形進行篩選，一致及 不一致的配對組合如圖 所示。將對應關係𝜒 → 𝜒′透過𝐹′到𝐹的向後配對，若滿足 (4.10)式時，表示不具備配對一致性。 argmin 𝜒 ̆𝜖𝐹 𝑑(𝜒̆, 𝜒′) ≠ 𝜒 (4.10) 經過(4.9)式及(4.10)式篩選後，可得到明確的特徵點匹配組合𝜒 ↔ 𝜒′。 4.2.3 強健特徵匹配 為了增強影相匹配，對於在𝐹中無法在 𝐹′內找到對應的特徵點，採用 KLT 追 蹤演算法[30]找尋在 𝐹′最相似的特徵點進行配對。KLT 演算法在欲尋找配對的特 徵點上開啟一個視窗，認定相鄰影像中亦有相對應的視窗，且視窗內的所有像素 位移皆相同，計算平方距離之和（Sum of Squared Distances, SSD）的最小值來找 到最有可能的對應點，得到增強特徵集合。增強特徵集合中，本研究使用隨機抽 樣一致（random consensus sampling, RANSAC）技術進行離群值排除。關於 RANSAC 大致流程如下圖 4-6[34]：
44 圖 4-6 RANSAC 流程示意圖[34] RANSAC 以隨機的方式選擇五個匹配組合，產生符合之模型。計算所有匹 配組合與模型的差異，重複上述步驟，再給定迭代次數結束或得到比預期更小之 差異值時，選定該模型做為標準。給定一閥值剔除錯誤匹配，留下穩健的特徵匹 配組合。在每次迭代運算中，產生之模型透過一本質矩陣(essential matrix)，匹配 組合中的兩個特徵點位置會有以下關係： 𝑝̇𝐸𝑞̇ < ε𝐹 (4.11) 點𝑝̇與𝑞̇為相互匹配的兩個特徵點座標，𝐸為本質矩陣，ε𝐹則為誤差容許值，理想 狀況下應趨近於 0。本質矩陣 E 是由選擇的五個匹配組合提供五點算法[31]進行 估計而來。並使用該矩陣來來評估所有匹配組合的 Sampson 距離，以 Sampson 距 離作為篩選的閥值，通常為 0.5。
45
