利用LMedS初解計算雜訊干擾之二張校正影像之間的精準相機外部參數

國

立

交

通

大

學

資訊科學與工程研究所

碩

士

論

文

利用

LMedS 初解計算雜訊干擾之二張校正

影 像 之 間 的 精 準 相 機 外 部 參 數

Accurate Camera Pose Estimation from the Noisy

Calibrated Stereo Images Based on LMedS Initial Results

研 究 生：鄭喬文

指導教授：陳 稔 教授

利用

LMedS 初解計算雜訊干擾之二張校正影像之間的精準相機外部

參

數

Accurate Camera Pose Estimation from the Noisy Calibrated Stereo

Images Based on LMedS Initial Results

研 究 生：鄭喬文 Student：Chiao-Wen Cheng

指導教授：陳 稔 Advisor：Zen Chen

國 立 交 通 大 學

資 訊 科 學 與 工 程 研 究 所

碩 士 論 文

A Thesis

Submitted to Institute of Computer Science and Engineering College of Computer Science

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master

in

Computer Science July 2008

利用

LMedS 初解計算雜訊干擾之二張校正影像之間的精準相機外部

參數

研 究 生：鄭喬文

指導教授：陳 稔

國 立 交 通 大 學

資 訊 科 學 與 工 程 研 究 所

摘要 本論文的目的在於利用同一物體不同角度的2 張影像，以其中對應點的資訊 來估測精確的相機外部參數，首先以LMedS 方法隨機計算出多組外部參數解， 並取誤差值較小的前幾名解。用MST(Minimum Spanning Tree)為基礎的分群計算 各小區域的橢圓體，對橢圓體細切再分群，直到橢圓體退化，以求得改進解及附 近的error surface 變化(=橢圓體的 eigen values 及 eigenvectors 分布)，然後依橢圓 體的eigen values 及 eigenvectors 來決定各維度執行搜尋的順序，最後執行以黃金 切割搜尋法(GSS)為基礎的 5D 空間最佳解搜尋法，以找到 5D 全域空間的最佳解。

Accurate Camera Pose Estimation from the Noisy Calibrated Stereo

Images Based on LMedS Initial Results

Student：

Chiao-Wen

Cheng Advisor：Zen Chen

Institute of Computer Science and Engineering College of Computer

Science

National Chiao Tung University

Abstract

Estimating accurate camera pose with extrinsic parameters from this thesis uses information of correspondent points in two different views from the same object. First, using LMedS method to obtain a lot of solutions of extrinsic parameters, and take the front solutions with smaller median error. Use clustering based on MST(Minimum Spanning Tree)to compute the information of ellipsoids of the local regions. Fine cut the ellipsoids and clustering again until the ellipsoids degenerate.Then we will get the solution with refinement and the variation of errorsurface( = eigen values and eigenvectors of ellipsoid). Decide the executing sequence of each dimesion by eigen values and eigenvectors of ellipsoid. Finally, executing the 5D optimization search method base on Golden section search(GSS) to find the optimized solution in 5D solution space.

Key words：camera estimation、extrinsic parameters、 optimized searching in multi-dimesion

致謝

首先由衷感謝我的指導教授陳 稔博士這兩年來對我不辭辛勞的指導與教 誨，讓我的邏輯思考能力了有很大的啟發，同時也讓我學習到研究學問必須持有 的態度與執著；當在緊要關頭時，老師也會陪著我一起面對眼前的難題、細心地 分析問題始末，並指引出一條正確的解決之道。在此，向老師致上最深的敬意！ 同時也感謝PAIS 實驗室的所有成員在研究上給我許多幫助，並時常給予最 大的鼓勵與打氣，感謝和我一起渡過700 多個日子的文昭學長，學長的幫助，讓 我在碩士之路上克服許多難關。感謝親愛的父母給予我充足的精神和物質生活， 讓我能心無旁鶩地將學業完成；也感謝我的女朋友詩雲，在這段時間陪我一同分 享內心喜悅、傾聽失意的牢騷，並在旁支援我繼續走下去，謝謝妳。 最後，在感謝這段時間陪伴我的所有人，僅以此篇論文獻給你們！ 鄭喬文2008/8 于 新竹交大 PAIS 實驗室

目錄

中文摘要 ……… i 英文摘要 ……… ii 致謝 ……… iii 目錄 ……… iv 表目錄 ……… v 圖目錄 ……… vi 第1 章 緒論……….. 1 1.1 研究動機與目標……….. 1 1.2 相關研究 ……… 2 1.3 研究方法介紹 ……… 3 1.4 論文組織 ……… 4 第2 章 外部參數的求法及其表示法……….. 6 2.1 從 Essential matrix 分解出外部參數……….. 6 2.2 評估方法的設定……….. 8 2.3 外部參數的 5D 角度表示法...……… 10 第3 章 LMedS 方法求解與初解之改進………. 14 3.1 用 LMedS 得到初解及分群…… ……….. 14 3.2 初解之改進………. 17 第4 章 搜尋全域性正解的離散值 ……… ………... 21 4.1 全域搜尋法………. 21

4.2 Golden section search(GSS) ……… 24

第5 章 實驗與結果 ……… 29 5.1 半雜訊實驗結果 ………. 32 5.2 全雜訊實驗結果 ………. 53 5.3 真實影像實驗結果………. 65 第6章 結論與未來發展方向 ……… 69 6.1 結論 ……… 69 6.2 未來發展方向 ……… 69 參考文獻 ……… 70

表目錄

表2.1 相機位移矩陣 T 轉換成球座標的正負判別……….. 13

表3.1 LMedS 結果的前 50 名 5D 解及其誤差值………. 16

表3.2(a) cluster 的 R mean & eigenvalues & eigenvector……….. 18

表3.2(b) cluster 的 T mean & eigenvalues & eigenvector……….. 18

表3.3 細切之後前 50 名解的 weighted code….………. 20 表5.1 實驗所使用的 10 張模擬影像……….. 31 表5.2 半雜訊實驗(加入 3σ = 3%的雜訊)初始解及最佳解的誤差值紀錄.…. 47 表5.3 半雜訊實驗(加入 3σ = 3%的雜訊)初始解及最佳解重建的 3D 模型 和ground truth 的比較結果………. 50 表5.4 全雜訊實驗(加入 3σ = 0.5%的雜訊)初始解及最佳解的誤差值紀錄…. 54 表5.5 全雜訊實驗(加入 3σ = 0.5%的雜訊)初始解及最佳解重建的 3D 模型 和ground truth 的比較結果………. 56 表5.6 全雜訊實驗(加入 3σ = 1.5%的雜訊)初始解及最佳解的誤差值紀錄….. 60 表5.7 全雜訊實驗(加入 3σ = 1.5%的雜訊)初始解及最佳解重建的 3D 模型 和ground truth 的比較結果……….………... 62

圖目錄

圖1.1 相機定位估測……… 1 圖1.2 重建精確的 3D 模型需要精確的外部參數……… 2 圖1.3 大略流程圖……… 4 圖2.1 求外部參數的流程圖……… 6 圖2.2 Essential matrix 分解外部參數的四種情形……….. 8 圖2.3 評估外部參數解好壞的 3 種方法……… 8 圖2.4 相機座標相對球座標的轉換關係 ……….. 11 圖3.1、最小平方中值法示意圖……….. 14

圖3.1 Minimum Spanning Tree 示意圖……… 17

圖4.1 Golen section search(GSS)示意圖……….. 21

圖4.2(a) GSS 改進步驟(1)……… 23

圖4.2(b) GSS 改進步驟(2)……… 23

圖4.3 error surface(水平軸為 ωR，垂直軸為 θR)………..………. 25

圖4.4 error surface(水平軸為 θR，垂直軸為 ψR)………... 25

圖4.5 三個參數空間的參數設定為Rlast、(Rhor,R )、(ver Thor,Tver)………… 26

圖4.6 全域搜尋法執行順序流程圖……….…. 27

圖5.1 在立方塊(cube)可見的 3 個面上的對應點編號………. 29

圖5.2 立方體(cube)在 2D 上的影像示意圖……… 30

圖5.3 為表 5.2 的示意圖(x 軸為實驗編號，y 軸為 median error)……... 47

圖5.4 半雜訊實驗(加入 3σ = 3%的雜訊)重建 3D 後的 9 個邊長長度 和原3D 對應邊長相差絕對值的平均(cm)………... 50 圖5.5 半雜訊實驗(加入 3σ = 3%的雜訊)重建 3D 後的 9 個邊長長度 和原3D 對應邊長相差絕對值的標準差……… 51 圖5.6 半雜訊實驗(加入 3σ = 3%的雜訊)重建 3D 後的 12 個角度 和原3D 對應角度相差絕對值的平均(度)……….. 51 圖5.7 半雜訊實驗(加入 3σ = 3%的雜訊)重建 3D 後的 12 個角度 和原3D 對應角度相差絕對值的標準差………. 52 圖5.8 為表 5.4 的示意圖(x 軸為實驗編號 , y 軸為誤差值)……….... 55 圖5.9 全雜訊實驗(加入 3σ = 0.5%的雜訊)重建 3D 後的 9 個邊長長度 和原3D 對應邊長相差絕對值的平均(cm)………. 57 圖5.10 全雜訊實驗(加入 3σ = 0.5%的雜訊)重建 3D 後的 9 個邊長長度 和原3D 對應邊長相差絕對值的標準差………. 57 圖5.11 全雜訊實驗(加入 3σ = 0.5%的雜訊)重建 3D 後的 12 個角度 和原3D 對應角度相差絕對值的平均(度)……….. 58 圖5.12 全雜訊實驗(加入 3σ = 0.5%的雜訊)重建 3D 後的 12 個角度

圖5.13 為表 5.6 的示意圖(x 軸為實驗編號，y 軸為 median error)………… 60 圖5.14 全雜訊實驗(加入 3σ = 1.5%的雜訊)重建 3D 後的 9 個邊長長度 和原3D 對應邊長相差絕對值的平均(cm)………. 62 圖5.15 全雜訊實驗(加入 3σ = 1.5%的雜訊)重建 3D 後的 9 個邊長長度 和原3D 對應邊長相差絕對值的標準差………. 63 圖5.16 全雜訊實驗(加入 3σ = 1.5%的雜訊)重建 3D 後的 12 個角度 和原3D 對應角度相差絕對值的平均(度)……….. 63 圖5.17 全雜訊實驗(加入 3σ = 1.5%的雜訊)重建 3D 後的 12 個角度 和原3D 對應角度相差絕對值的標準差………. 64

第 1 章 緒論

1.1 研究動機與目標

相機定位估測(camera pose estimation) 是一個重要研究題目，主要是由影像 估測相機拍攝時對於某一座標系的的旋轉(rotation)及位移(translation)資訊，如 圖 1.1，可用在電腦視覺技術中的空間定位相關應用上，如：人臉辨識、門禁 系統、多張影像 3D 重建、機器人視覺(robot vision)，還有視覺式之同時定位 與製圖技術(visual simultaneous localization and mapping)等等。

估測準確性是一個很重要的指標，只有當準確度提高到某一個程度，才具有 實用價值，但對於原始資料有可能包含雜訊，則容易造成難以估測準確的結果。 對於此影像含有雜訊的問題，傳統的演算法僅能得到一個粗估的解，離正確解尚 有一段差距。在過去幾年的應用中，此問題尚未被重視，主要的原因在於影像解 析度約百萬相素的等級，就算影像對應點有雜訊的干擾，估測的解仍能得到不錯 的結果。但是隨著影像感測元件解析度的提昇，現在的相機畫素已達千萬相素的 等級，原本可忽略的影像誤差在如此高的解析度下開始會對結果產生影響，進而 造成後續資料處理(如影像式三維建模技術)的誤差。因此本論文的研究希望能提 昇影像辨認的學理與技術水準，增進其辨認的準確率及可用性。 本論文的目標將探討一個使用2 張影像及其對應特徵點的資訊來估測精確相 機外部參數的方法，即使在雜訊干擾下也能得到相機5D 參數離散空間(discrete solution space)的最佳解。預期研究成果可逐步應用於電玩娛樂、身份證明、人機 介面、安全監控、關卡控制、嫌犯追蹤等用途上。 圖1.1 相機定位估測

1.2 相關研究

本節將討論在相機定位估測的問題裡精準外部參數的必要性，還有目前在此 問題中所使用的解決方法。

在今年的CVPR 的 Y. Furukawa, and J. Ponce,”Accurate Camera Calibration from Multi-View Stereo and Bundle Adjustment”，討論到在建立 3D 模型時，越精 確的模型需要越精確的外部參數來達成，如下圖1.2。目前在使用機械手臂控制 下，實際外部參數最多精確到0.01°，而因此造成的誤差對於要建立精確模型還 是不足，所以有開發估測精準外部參數方法的必要性。

圖1.2 重建精確的 3D 模型需要精確的外部參數 一般而言，目前在估測精準外部參數的課題上可大致分成3 種方法： Linear methods：快速但較無法保證會收斂到正確解(最佳解)，主要應用於

real-time issue。例如：A. Ansar and K. Daniilidis, “Linear pose estimation from points or lines”, PAMI 2003。

Non-linear methods：包括 Gradient-descent, Gauss-Newton, Levenberg-Marquardt, conjugate gradient methods, etc。缺點是只能應用在有唯一解 的情況，當有區域極值(local minimum)時，結果有可能會陷 入區域極值中。

Stochastic methods：包括Genetic algorithm, Taguchi method, simulated annealing, etc。解決以上2種方法的缺點，但會需要以大量的計算來換 取準確的結果。例如：劉嘉哲, “使用直交退火演算法解相機 校正問題”, 逢甲大學資訊工程學系, 碩士論文,民92； “Global Optimization through Searching Rotation Space and Optimal Estimation of the Essential Matrix“,[ Richard I. Hartley and Fredrik Kahl , 2007]；Bundle Adjustment

[Yasutaka Furukawa and Jean Ponce , 2008]

此外，由於通常有許多的特徵點，所以會搭配robust estimation，除了可以剔 除outlier之外，還能稍微減少noise的影響(有機會挑出較小noise的點)。常用的如 M-estimator, Least Median Squares[Peter J. Rousseeuw ,1984], and RANSAC (RANdom Sample Consensus) [Fishler & Bolles, 1981]。本論文採用搭配LMedS的 方法先求出初解。

1.3 研究方法介紹

本論文要執行的工作包括建立模擬及實拍影像，以LMedS 求初始解，初解的 改進，以改進解為中心做GSS 為基礎的 5D 空間最佳解搜尋等工作，細項包括： (1) 用 2.3 節的方法，對 5 個角度的設定來取得模擬的影像，或是直接拍攝實物。 (2) 以 LMedS 方法隨機計算出多組外部參數解，並取誤差值較小的前幾名解。 (3) 用 MST(Minimum Spanning Tree)為基礎的分群計算各小區域的橢圓體。 (4) 對橢圓體細切再分群，直到橢圓體退化，求得改進解及附近的 error surface 變

化(=橢圓體的 eigen values 及 eigenvectors 分布)。

(5) 依橢圓體的 eigen values 及 eigenvectors 來決定各維度執行搜尋的順序。 (6) 執行以黃金切割搜尋法(GSS)為基礎的 5D 空間最佳解搜尋法，以找到 5D 全域

空間的最佳解。 執行概略流程如下：

Input：View1 & View2 的 對應特徵點

LMedS 求 t 次相機外部參數解當作初解

分析LMedS 前 p 個最小 median error 解，建立其 Minimum Spanning Tree(MST)

對MST 做分群(Clustering)，求各群 Covariance matrix 橢圓 體之eigen values 及 eigen vectors

對各群橢圓體細切，得到初解的改進 refinement，直至細 切橢圓體退化(至少一個 eigen value 接近 0)

啟用Golden section search(GSS)搜尋相機外部參數的最佳 解：

(1)轉換 5D Camera parameter space 為 3D Rotation space 與 2D translation space，(2)進而再用橢圓體 eigen values 及 eignen vectors 定義三個參數空間:1D rotation eigenspace +2D rotation space +2D translation space，(3)使用 Golden section line search 連環搜尋相機外部參數的最佳解

Output：最佳近似解& error value 圖1.3 大略流程

1.4 論文組織

本論文共分為六章，首章為本論文研究方向的概括性介紹，並介紹相關主題 中使用的各種方法；第二章介紹如何用2 張影像對應點的資訊求出外部參數解，

並將外部參數以5 個角度來表示，同時應用在建立模擬影像的方法中；第三章介 紹如何用LMedS 求出初始解，並以分群及細切的方法找到包含最佳解的小區域 橢圓體及初步改進解；第四章將說明如何以第3 章的改進解為中心，在附近範圍 以Golden Section Search(GSS)演算法來找到更加精確的外部參數解；第五章為本 論文實驗相數據及分析：第六章為本論文的結論探討及後續工作的規劃。

第 2 章 外部參數的求法及其表示法

本章將說明如何使用二張影像中的已知對應點(至少八點)求出基礎矩陣 (fundamental matrix) Forigin，經由 Forigin及兩張影像各自的內部參數 K 及 K’可計

算出essential matrix Eorigin。此外亦使用essential matrix 本身的特性來進行修正，

經修正後再分解出二張影像之間的rotation matrix R 及 translation T。利用修正後 的 Emodfied及內部參數可計算出較正確的fundamental matrix Fmodified，最後將外部

參數 R, T 以 5 個角度 ψR , θR , ωR , ψT , θT來表示，使解空間從6D 降為 5D，求外 部參數的完整流程如圖2.1 所示。 圖2.1 求外部參數的流程

2.1 外部參數的求法

令二張影像為 及 ，且有 組的對應點，此 組的對應點可利用現有 的特徵點擷取與對應之相關技術事先計算。令 為 上的對應點座 標， 為 上的對應點座標，其中 及 中的對應點與兩影 像之間的fundamental matrix 有以下關係： 1 View View₂ n n T i i i (u ,v ) P = View₁ T ' i ' i i (u ,v )

P′= View2 View1 View2

F ' ₀ i i P FP =

[

]

0 1 1 33 32 31 23 22 21 13 12 11 = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⇒ _i i ' i ' i v u F F F F F F F F F v u (2.1) 由上式可知，若將 令為1，僅需八組對應點即可利用 least-square method 來計 算 ，計算方法如下： 33 F F

11 21 33 _{9 1} 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 0 0 * * 0 ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' x n n n n n n where F F A f A F u u u v u v u v v v u v u u u v u v u v v v u v u u u v u v u v v v u v A u u u v u v u v v v u v u u u v u v ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ = 9 1 1 1 1 ' ' ' '1 n n n n n n Nx u v v v u v ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ (2.2) 接下來對A 做 SVD 分解，

其中f 的最佳解為對應 A 最小 eigen-value 的 unit eigen-vecter。

令 ，由已知的內部參數 K 及 K’及 ，可求出Essential matrix ⎡ ⎤ ⎢ = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ origin 11 12 13 21 22 23 31 32 33 F F F F F F F F F F ⎥ _F origin 'T Eorigin =K ForiginK，對 做SVD 分解得到 ，在有誤差的情 況下 ，對D 做修正，使 Eorigin origin = T E UDV ⎡ ⎤ ⎢ = ⎢ ⎢ ⎥ ⎣ ⎦ 0 0 0 0 0 a D b c ⎥ ⎥ 0 + ⎡ ⎤ ⎢ ⎥ = _⎢ + _⎥ ⎢ ⎥ ⎣ ⎦ modified ( ) / 2 0 0 0 ( ) / 2 0 0 0 0 T a b E U a b V K ， 同時計算出 _'T 1。 Fmodified =K − Emodified − 相同的Essential matrix 會有 4 種不同的外部參數(=R、T)分解可能[1]，計算如下： 對Emodified做SVD 分解得到 [ 1 2 3] T T Emodified =UDV = u u u DV ，令 ，R 有2 種可能為 − ⎡ ⎤ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ 0 1 0 1 0 0 0 0 1 W 1 T R =UWV 及 2 T T R =UW V ，T 有 2 種可能為T₁=u₃及T ₃，從幾 何意義上來看，如下圖2.2(a)~(d)，其中 A、B 代表 2 台相機，圖 2.2(a)表示 3D 點同時在2 台相機前方，圖 2.2(b)表示 3D 點在 2 台相機後方，圖 2.2(c)表示 3D 點在A 相機後方但在 B 相機前方，圖 2.2(d)表示 3D 點在 A 相機前方但在 B 相 機後方，而這4 種情形會產生相同位置的 2D 對應點，但僅有在圖 2.2(a)情形下 的外部參數(R、T)為我們所求。 = − 2 u

A B ● ● ● (a) B A ● ● (b) ● A B’ ● (c) ● ● (d) A B’ ● ● ● 圖2.2 Essential matrix 分解外部參數的四種情形

2.2 評估方法的設定

圖2.3 評估外部參數解好壞的 3 種方法 如圖2.3，我們要評估求出的外部參數解好壞有以下 3 種方法： 方法1：

以 F_origin計算 中對應點在 上的epipolar line 到 上相對應點的距離

加上 中對應點在 上的epipolar line 到 上相對應點的距離，計算出

1

View View₂ View₂ 2

n 個值並取中位數，即為 n 個值排序後的第⎡⎢n/ 2⎤⎥個值，公式如(2.3)~(2.5)； 令 為 上的對應點座標， 為 上的對應點座 標，其中， T i i i (u ,v ) P = View₁ P_i' ( ', ')= u v_{i i} T View₂ i

P 對應在View₂上的Epipolar Line 的係數為：

⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ = _{⎜ ⎟ ⎜}= _⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ origin 1 2 3 ' '( ) * ' 1 ' i i u Line i F v (2.3) ' 對應在 上的Epipolar Line 的係數為： i P View₁

(

)

(

)

= _origin = _{1 2 3} ( ) _i' _i' 1 * Line i u v F (2.4) 第i 個 epipolar line 到對應點的距離和的值如下：

(

)

(

)

= + + + 2 2 2 2 1 2 1 2 ( ) * 1 ' ' 1 * '( ) ( ) ( ) ' ' T i i i i Line i u v u v Line i Error i (2.5)

最後得到的Error value = median Error i i{ ( ), =1 ~n} 。

方法2：

計算方法同(1)但 F_origin改用F_modified。 方法3：

以 分解出外部參數，並以已知的內部參數 K 及 K’求出 及 的投 影矩陣(projection matrix)，如下(2.6)：

E_modified View₁ View₂

[ ]

= 1 | 0 PM K I , PM₂ =K R T'

[

|

]

(2.6) 接下來，由 、 、 、 重建出3D 點，再將 3D 點投影回 及 產生新的2D 點，計算新的 2D 點和 、 之間的距離，同樣在n 個值取中位數 當做誤差值。 i P P_i' PM₁ PM₂ View₁ View₂ i P P_i' 重建3D 點的方法如下： 1 View 上的點 和 之間有(2.7)式的關係，其中 Xi = (xi , yi ,zi ,1)T為原始的3D 點 i P PM₁ ⎛ ⎞ ⎛ ⎞ _{⎜ ⎟} ⎜ ⎟_{× ⎜ ⎟ =} ⎜ ⎟ _{⎜ ⎟} ⎜ ⎟ _{⎜ ⎟} ⎝ ⎠ ⎝ ⎠ i 1 ( ) 1 1 i i i i i x u y v PM z 0

3 1 3 2 1 1 1 1 iT 1 1 ( ) ( ) 0 , ( ) ( 1 1 1

where pm = ith row of PM

i i i i i i T T T T i i i i i x x x y y y u pm pm v pm pm z z z ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⇒ − = − ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ) 1 i i i x y z = T (2.7) 由(2.7)可以導出 3 1 1 1 3 2 1 1 3 1 2 2 3 2 2 2 4 1 0 0 * * , 0 ' 1 0 ' T T i i T T i i i T i i T T i x where x u pm pm y v pm pm A X A A z u pm pm v pm pm ⎡ − ⎤ ⎛ ⎞ ⎡ ⎤ ⎢ ⎥ ⎜ ⎟ ⎢ ⎥ ₋ ⎢ ⎥ ⎜ ⎟ ⎢ ⎥ = = =_{⎢ ⎥} ⎜ ⎟ ⎢ ⎥ − ⎢ ⎥ ⎜ ⎟ ⎢ ⎥ ₋ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ ⎣ ⎦ (2.8) 接下來對A 做 SVD 分解，

其中重建3D 點X_i的最佳解為對應A 最小 eigen-value 的 unit eigen-vecter。

重建出3D 點後，將其投影 及 產生新的2D 點 及

，計算新的2D 點和 、 之間的距離

1

View View₂ P_{i new_} =(u_{i new_} ,v_{i new_} )T

_ ' ( _ ', _ ')

T

i new i new i new

P = u v P_i P_i'

2 2 2

_ _ _ _

( ) ( i new i) ( i new i) ( i new' i') ( i new' i')

Error i = u −u + v −v + u −u + v −v 2 (2.9)

最後得到的Error value = median Error i i{ ( ), =1 ~n}

本論文所使用的評估方法為方法2，原因如下：經測試，以方法 1 所求出有 小誤差值的外部參數解，和正解有很大差距，和我們預期不同(有小誤差值的解 要接近正解)，所以不適用，而許多 5D 解以方法 2 和方法 3 所計算出的誤差值的 排序相當，但方法3 所需的步驟較多，故使用方法 2 來當做我們評估解好壞的方 法。

2.3 外部參數的 5D 角度表示法

一開始外部參數為R(rotation matrix)及 T(translation matrix)，一般而言，R 可 以用三個角度：對x 軸旋轉角度 θx、對 y 軸旋轉角度 θy、對z 軸旋轉角度 θz來

表示，如(2.6)式：

Rotation Matrix R = ( )Rθ_x ⋅ ( )Rθ_y ⋅ ( )Rθ_z

= (2.6)

1 0 0 cos( ) 0 sin( ) cos( ) sin( ) 0 0 cos( ) sin( ) 0 1 0 sin( ) cos( ) 0 0 sin( ) cos( ) sin( ) 0 cos( ) 0 0 1

y y z z x x z z x x y y θ θ θ θ θ θ θ θ θ θ θ θ ⎡ ⎤ − ⎡ ⎤ ⎡ ⎢ ⎥ ⎢ ₋ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎣ ⎦_{⎣ ⎦}⎣ ⎤ ⎥ ⎥ ⎥⎦

而 T real=[Tx , Ty ,Tz]T表示2 台 camera 實際上在世界座標裡的相對位置，但以 2 張平面的影像僅可獲得T 的方向資訊，T 的真實長度沒辦法求出，所以用 2 張平 面影像求出的T 只需用 2 個維度的值來表示即可。計算方法如 2-1 節，令求出的 T = [tx , ty ,tz]T，其長度為1，以球面座標系統的 2 個角度(ψT ,θT)表示，轉換方法 如後述。 R 若以(θx , θy, θz)表示，很難直觀的看出 camera 旋轉後的方向，而且由以下三 個原因，改用(ψR ,θR , ωR)表示之，轉換方法如後述。 „ 和 T 同樣以球面座標(經緯度)模式表示 „ 能簡單明確顯示相機旋轉的方向 „ 避免萬向鎖(Gimbal lock )問題 yc xc ystep3 ψR θR ψR xstep2 x0 zstep1 xstep1 ystep1 Step 2 Step 3 Step 1 r 1 z0 zc ystep2 zstep2 xstep3 Step4：旋轉ωR y0 圖2.4 相機座標相對球座標的關係 初始設定為物體的原點置於球座標的球心，其z 軸指向北極，₀ x 軸定在東經₀ 0°，y 軸定在東經 90°，且物體不動。 ₀ 相機一開始位置與物體一致，然後移動到球面上，半徑r1，移動順序如圖2.4 從Step1~Step3，使x 、_c y 的相機平面切於球面，相機光軸_c z 指向球心。此時相_c 機原點座標(x ,_c0 y ,_c0 z_c0)移到東經θ_R角，北緯(90-ϕ_R)角，其數學式為 0 1 0 1 0 1 sin cos sin sin cos c R c R c R x r y r z r R R ϕ θ ϕ θ ϕ = = = (2.7) 則相機對原始球座標的位移矩陣為T =_⎣⎡t_x t_y t_z⎦⎤T =

[

x_c0 y_c0 z_c0

]

T。 再看Step3 結束後，相機座標系的三個軸方向如下：

θ θ ϕ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢₌ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ = − = × 其中 step2 c0 c0 c c0 c0 c0 c0 o 3 z0 R 2 z0 R y0 R step1 step1 0 x 3 c 3 x -x z = -v =0- y -y z -z R ( ) =R ( )R (180 )x , x = x z step step step step x x y x (2.8) 其次，再將此x 、_c y 對球面的切平面對此時的_c z 軸旋轉_c ω_R角，最後的x 、_c y 軸_c 方向算法如下： ω ω = c_final step3 c_final step3 x = ( )x y ( )y c c z R z R R R (2.9) 其中 (ω c z R R )為對z 方向軸旋轉_c ωR角度，計算如下： Given ⎡_⎣ ⎤_⎦ 1 2 3 c c c c z = z z z T * ( ) (cos )( ) (sin ) c z R R R R ω = +A ω I−A + ω A (2.10) where ∗ ⎡ ⎤ ⎡ − ⎤ ⎢ ⎥ ⎢ ⎥ =⎢ ⎥ =⎢ − ⎥ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ _{⎣ ⎦} ⎣ ⎦ 1 1 1 2 1 3 3 2 2 1 2 2 2 3 3 1 3 1 3 2 3 3 2 1 c c c c c c c c c c c c c c c c c c c c c c c c z z z z z z 0 z z z z z z z z , z 0 z z z z z z z z z 0 A A 最後相機座標系與原始參考球座標系之間的3D 旋轉矩陣 R 即可求得 (2.11) ⎡⎣ c_final c_final c R= x y _{z ⎤⎦} 以上為相機對球座標旋轉及位移的執行過程，假設相機#1 對球座標旋轉矩陣 為R，位移矩陣為₁ T₁，相機#2 對球座標旋轉矩陣為R ，位移矩陣為₂ T ，則可由₂ R 、₁ 2 R 、T₁、T₂求出相機#2 對相機#1 的相對旋轉矩陣R 及位移矩陣₁₂ ，即把相機 #1 的座標系設為 ， 12 T 1 0 0 0 1 0 0 0 1 R I ⎡ ⎤ ⎢ ⎥ = = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ 0 0 0 0 T ⎡ ⎤ ⎢ ⎥ = = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ，R 及₁₂ T₁₂算法如下： 12 2 1 12 2 2 1 1 T T R R R T T R R T = = − (2.12) 如果要以Essential matrix 分解出來的 T 計算出(ψT,θT)，轉換方法如下： 首先，因T 為 unit vector，即表示在(2.7)式中的 r1=1，令 T x y z T = ⎣⎡t t t ⎤_{⎦ ，其中} ，以反三角函數及tz算出ψT的值，再用ψT及tx、 ty算出θT，但反三角函數非1 對 1 函數，所以要同時用 tx, ty, tz三個值來判定ψT、 θT正負，計算角度大小及判定正負方法如下表： 0o ≤ ≤ϕ 180 , 180o − o ≤ ≤θ 180o

用tx，ty，tz正負判定ψT、θT範圍 tx ty tz ψT θT ≧0 ≧0 ≧0 0≦ψT≦90 0≦θT≦90 ≧0 ≦0 ≧0 0≦ψT≦90 -90≦θT≦0 ≦0 ≧0 ≧0 0≦ψT≦90 90≦θT≦180 ≦0 ≦0 ≧0 0≦ψT≦90 -180≦θT≦-90 ≧0 ≧0 ≦0 90≦ψT≦180 0≦θT≦90 ≧0 ≦0 ≦0 90≦ψT≦180 -90≦θT≦0 ≦0 ≧0 ≦0 90≦ψT≦180 90≦θT≦180 ≦0 ≦0 ≦0 90≦ψT≦180 -180≦θT≦-90 表2.1 相機位移矩陣 T 轉換成球座標的正負判別

第 3 章 LMedS 方法求解與初解之改進

本章將說明如何應用LMedS 方法來求得初解，首先在多組對應點中任意挑選 部份的對應點(至少八組)來計算 fundamental matrix，可以得到一組外部參數解(如 第2 章所述)，反覆進行此步驟直至產生出足夠多(如 4000 次)的外部參數解，每 一組解皆有相對應的誤差值(如 2.2 所述)，並在這些解當中，取出誤差值較小的 前幾個解。 在此假設全域最佳解與這些較佳的解會相當的接近，因此針對這些較佳的 解，利用Minimum Spanning Tree(MST)來進行分群(clustering)，得到若干個可能 包含全域最佳解的分群，並利用能包含各分群範圍的橢圓體來描述這些分群。接 著對這些橢圓體進行切割(refinement)，來探查橢圓體內是否有更佳的解，並以此 獲得改進的解。

3.1 用 LMedS 方法求初解及分群

在求外部參數時，如果沒有雜訊（noise）或 outlier 的干擾，最小平方法（least squares）便可以找到最好的答案。但在找到對應點，或者是要選取對應點的時候， 除了模擬資料（simulate data）之外，否則一定會有雜訊或 outlier。最小平方中 值法（LMedS）是將所有的資料作排序後，取中間的數值，並找出最小的那一組 資料，圖3.1 是任取兩個點來建構一條直線做說明：{（ ）| i =1, 2, 3, … m}， r 代表的是殘餘量（residual）：r = |ax+by+1|，在此為點到線的距離，共有 n 組 資料點，取了m 組點來畫線。 i i b a , Median ……… Median ……… 3 3,b a 2 2

,b

a

a ,m bm 1 1

,b

a

從 小排到大 2 11 r 2 12 r Min 2 1 n

r

(Residual) 2 21 r 2 22 r 2 2n

r

r

_mn2 圖3.1、最小平方中值法示意圖 由於最小平方中值法並不需要設定參數或是門檻值，像RANSAC 需要設定兩 個重要參數，而M-estimator需要給定足夠好的初始值，才有可能在幾個 iteration 後得到較好的結果；最小平方中值法主要是利用資料點群體力量，來找出結果， 所以在目前強健式估測法中，仍然常被拿來應用或比較。最小平方中值法有一個 相當重要的前提假設：outlier 不能佔超過全部資料的點一半，否則取到的中位數 會有問題（若超過一半則取到的點可能已經是錯誤資料了），因此又有別的方法

因應而生。 我們先用LMedS 的方法在 2 張已校正影像的 n 組對應點中，重複 t 次隨機挑 選k 組(至少八組)對應點，經由計算 Fundamental matrix，加上已知的內部參數， 可得到Essential matrix，再求出 t 個外部參數解及其誤差值，如第 2 章所述，並 從中挑出有較小誤差的解，以期望可降低noise 造成的影響。 由於此t 個外部參數解在 5D 解空間隨機分布，在此假設全域最佳解與這些較 佳的解會相當的接近。所以在此t 個外部參數解中，選出誤差較小的前 p 個解， 以期望在包含此p 個解的橢圓體範圍中能夠得到接近全域最佳解的初始解。 以下表3.1 為例，此例中 n=19, t=4000, p=50, k=8。 num ψR θR ωR ψT θT median error 與正解的 Euclidean distance 1 29.6682 141.4966 -121.405 74.1991 317.2275 6.6906 9.6442 2 32.7053 142.2839 -123.827 71.7361 317.8237 6.8135 6.0297 3 28.8298 143.0459 -124.37 69.4278 317.8257 7.8483 9.0224 4 27.6632 137.5118 -117.67 72.2759 313.0297 9.4849 12.6352 5 42.9453 141.0924 -123.253 69.5353 317.7244 9.4916 5.8443 6 36.0949 142.4432 -123.678 66.3573 317.6953 10.5248 3.1184 7 36.1698 141.7338 -122.765 70.313 318.3048 10.5706 2.605 8 41.2414 139.2393 -121.476 67.8386 315.594 13.5777 4.7286 9 39.5726 139.5241 -121.497 65.219 315.1695 14.2116 4.7397 10 40.441 139.1832 -121.745 72.2206 316.3319 16.9895 5.3676 11 26.8339 141.8992 -122.442 65.9585 316.4764 17.4494 10.8159 12 35.952 140.8474 -122.268 70.2366 316.5295 18.3644 2.3842 13 36.1866 139.7995 -121.099 68.0759 316.4651 19.4008 2.2254 14 40.1426 138.9176 -121.799 73.4008 315.8542 20.5001 6.2819 15 41.1813 138.0518 -120.737 72.0803 315.7656 21.2277 6.4749 16 38.5412 140.4795 -122.203 72.246 317.9028 23.5363 4.1221 17 42.0116 138.916 -121.46 59.9859 313.762 23.9649 10.4997 18 31.2562 140.4504 -123.067 82.2089 316.9913 26.2528 15.0929 19 27.7324 138.1323 -121.695 84.0006 314.6408 26.7347 18.7062 20 44.0845 138.0402 -121.22 69.1971 314.5459 27.0404 7.9877 21 39.2174 139.7523 -122.002 65.7418 314.9225 28.2056 4.2131

24 40.3016 137.4462 -120.388 74.8545 315.675 30.5953 8.3298 25 37.5271 140.8688 -122.647 79.4313 318.0791 30.6358 11.0618 26 38.5522 136.7595 -119.898 77.9604 315.1862 30.8471 10.9956 27 36.5342 144.2826 -126.074 62.7031 318.0967 32.9034 7.6481 28 33.0913 137.455 -120.328 82.0094 314.8837 33.2035 14.9936 29 30.3879 143.2053 -124.89 75.1401 317.5832 33.9414 10.2232 30 37.9135 139.0839 -120.574 71.2325 316.901 34.4179 3.9101 31 36.6619 141.0545 -123.984 76.6479 316.8203 35.0664 8.4294 32 36.5823 139.7174 -121.091 68.635 317.0821 35.1031 1.9552 33 32.7626 141.1524 -125.432 57.9692 315.6235 35.1527 11.8952 34 36.1102 141.4645 -123.204 63.0327 315.4508 36.0942 5.8658 35 41.574 140.6714 -122.425 73.6503 318.8567 36.4867 6.9624 36 35.6502 139.8613 -122.237 61.0371 316.0397 36.9253 7.7284 37 41.6971 137.9025 -121.734 73.2427 314.1947 38.0765 7.8595 38 43.509 136.4972 -121.111 75.9487 313.9754 38.1183 11.2981 39 36.6372 137.454 -119.779 65.4183 314.4691 38.6548 6.0128 40 41.8516 135.1086 -119.111 73.3729 311.6951 39.1149 11.0038 41 28.8276 140.0757 -120.71 70.9901 316.8792 40.0534 9.0647 42 33.8341 138.2457 -122.813 86.4725 315.5252 40.202 18.6781 43 38.6968 132.1669 -118.152 79.2233 309.5387 40.6239 16.5381 44 39.2984 137.5257 -122.517 81.5958 314.5076 41.1943 14.0438 45 39.4081 141.5816 -123.476 78.4222 317.7742 42.6323 10.3245 46 39.925 133.0681 -119.809 79.712 310.2415 42.6791 15.894 47 32.8718 141.3361 -128.416 55.1156 315.3274 44.0129 15.3785 48 26.8601 140.9493 -122.368 61.3312 316.2955 44.4691 12.6475 49 38.3801 134.2182 -119.886 81.7184 313.1246 46.231 15.7173 50 21.7216 148.5722 -133.452 57.1188 315.7457 46.3699 23.5099 表3.1 LMedS 結果的前 50 名 5D 解及其誤差值 因為解空間很複雜，不止包含一個極小值，還可能分成許多的區域(local trap)，各含有一個區域極值，而我們取得的t個外部參數解就分布在這些區域中， 為了取得有可能包含全域最佳解的這些區域範圍(local regions)，接下來，以解的 5D值計算Euclidean distance，以此為依據來建立Minimum Spanning Tree(MST)， 將t個外部參數解之間距離的遠近表達出來，為之後分群做準備。如下圖3.1，x 軸為1~50點的編號，距離最近的2點其編號相鄰且有連接線，y軸即為Euclidean

distance，紀錄有連接在一起的點之間距離。

圖3.1 Minimum Spanning Tree 示意圖

以適當的gap 將 MST 分成數個群(clusters)，當只有分出 1 個群(cluster)時，則 直接對此群(cluster)做細切，細切方法如 3.2 節所述，以求得改進的解，但是當分 出數個群(clusters)時，必須同時對這些群(clusters)做細切再判斷較有可能包住正 解的cluster：

„ 若 cluster 是離正解較遠的區域最小值(local trap or local minimum)，該 cluster 的橢圓體依eigen values 及 eigen vectors 做細切，其誤差值改進的空間不大。 „ 若 cluster 較有可能包住正解，該 cluster 橢圓體依 eigen values 及 eigen vectors

做細切，可得到有較多改進的新解。

3.2 初解之改進

本節將說明如何分析解分群之Covariance matrix 求得包住各群(clusters)的橢 圓體，及如何對橢圓體細切來求得誤差值較小的解(=改進解)。 由於5D 的橢圓體結構太過複雜，難以分析，將互相之間較獨立的 R、T 分開 計算出covariance matrix ΣR、ΣT，如下(3.1)式 _ _ 1 1 ( )( n T R i R R i R i x mean x mean n = Σ =

∑

− − _R) (3.1a) _ _ 1 1 ( )( n T T i T T i T i ) T x mean x mean n = Σ =

∑

− − (3.1b)

並求出R、T 各自的 eigenvalues、eigenvectors 及 mean 如下表 3.2。 R ψR θR ωR mean 37.8616 139.3263 -121.8653 eigenvalue _{4.2152 2.4175 0.6472} 軸1 軸2 軸3 eigenvector -0.8384 -0.5411 0.066 0.4893 -0.6937 0.5286 -0.2403 0.4754 0.8463

表3.2(a) cluster 的 R mean & eigenvalues & eigenvector

T ψT θT mean 71.9678 315.8556 eigenvalue _{5.2282 1.9846} 軸1 軸2 eigenvector -0.9888 0.1491 0.1491 0.9888

表3.2(b) cluster 的 T mean & eigenvalues & eigenvector 接下來計算橢圓體各軸的長短，最長軸半長的長度如下：

the longest axis of ellipsoid = eigenvalue_max c_max (3.2)

其它軸半長的長度依 eigenvalue的比例可得。 其 中 _max max{ ( )T 1( ) , 1 i i i c = c = x −mean Σ− x −mean i= ~ }n ， 為 群 中 各 點 的 Mahalanobis distance 的最大值。 在計算出橢圓體各軸長度後，因為我們比較關心誤差值較小的部份，即是橢 圓體較中間的部份，所以只對橢圓體中間2/3 的部份做細切，細切方法如下： (a)對 5 軸各自等分成 4 等份，其中每一等份的長度設為(L1 , L2 , L3 , L4 , L5 )

(b)在每軸上取得 5 個等分點：(meani-2Li , meani-Li , meani, meani+Li, meani+2Li) ,

i=1~5

(c)5 軸上每一軸的 5 個值即可搭配得出在橢圓體內部 2/3 的 55組新解，計算出 55組新解的誤差值。

化而且延橢圓體形狀來分群，所以之後在建立MST 時，其依據改以 5D 的 weighted distance = [(n1i-n1j)2+( n2i-n2j)2+( n3i-n3j)2+( n4i-n4j)2+( n5i-n5j)2]，其中 nki = {±1, ±2, 0} , nkj = {±1, ±2, 0} , k=1~5，如下例表 3.3： num 軸1 軸2 軸3 軸4 軸5 1 2 0 1 -1 1 2 2 -1 1 -1 2 3 -1 -1 1 1 1 4 -1 -2 1 1 2 5 -2 -2 -2 2 1 6 1 0 -2 0 0 7 -2 -1 -1 2 0 8 -2 -1 -2 2 0 9 0 -1 -1 0 1 10 1 0 -1 0 0 11 2 1 2 -1 0 12 -1 0 2 1 0 13 -2 0 0 2 -1 14 0 -2 -2 0 2 15 1 -1 -2 0 1 16 1 1 -1 0 -1 17 -2 1 1 2 -2 18 0 0 0 0 0 19 2 0 2 -1 1 20 0 -2 0 1 2 21 0 0 1 1 0 22 2 -1 0 -1 2 23 0 -1 0 1 1 24 -1 -1 2 1 1 25 -1 -2 0 1 2 26 1 2 0 0 -2 27 -1 1 2 1 -1 28 0 -1 1 1 1 29 1 1 0 0 -1

31 0 -2 -1 1 2 32 0 1 -1 0 -1 33 1 -2 -2 0 2 34 -2 1 0 2 -2 35 0 0 -1 0 0 36 -1 -1 2 2 1 37 0 -2 -1 0 2 38 -2 0 -1 2 -1 39 0 1 2 1 -1 40 2 -1 2 -1 2 41 1 -1 2 -1 2 42 1 -2 0 1 2 43 0 0 -2 0 0 44 1 0 2 -1 1 45 2 2 -1 0 -2 46 1 2 1 0 -2 47 -1 -2 -2 2 1 48 1 -1 1 -1 2 49 0 0 2 1 0 50 -1 -1 -2 2 0 表3.3 細切之後前 50 名解的 weighted code 分群後，計算ΣR、ΣT，如果橢圓體退化(即最小的 eigenvalue 接近 0)，則不再 繼續細切，以目前最佳解為中心執行全域搜尋法，如果橢圓體沒退化，則算出橢 圓體各軸長度，並再次細切。

第 4 章

搜尋全域性正解的離散值

本章將說明如何以第 3 章的改進解為中心，在附近範圍以 Golden Section Search(GSS)演算法來找到更加精確的外部參數解，前提條件為正解位在以改進 解為中心的某特定範圍內。

以改進解為中心，將5D 的相機參數空間(camera parameter space)分成 3D 旋 轉空間(Rotation space)及 2D 位移空間(Translation space )，進而再用橢圓體的各 軸長及方向定義三個參數空間:1D rotation eigen-space +2D rotation eigen-space +2D translation eigen-space，並以橢圓體各軸所對應的 eigen values 大小來決定 5 個維度在最佳化搜尋時的搜尋順序。

接著使用以GSS 為基礎進行 5D 空間的最佳解搜尋，以期找出具全域最小誤 差值的最佳解。

4.1 黃金切割搜尋法(Golden Section Search)

Golden Section Search(GSS)是在一維空間中進行最佳解的搜尋，與 Binary Search 不同的是 GSS 每次以黃金比例數值 γ≒0.382 當做切割比例，可以幾何速 率找到最佳解，但GSS 僅適用於 unimodal function，表示在搜尋範圍中僅有一個 極大或極小值的情況，GSS 才能找到準確的值。若尋找的非一維空間或是非 unimodal function，GSS 就無法適用，因此在下一節，解決不能使用在多維空間 的方法是將5D 全域依 eigenvalue 大小依序在各單一維度上進行 GSS 以找到全域 最佳解，而且本節也將提出改良後的GSS 以適用在非 unimodal function 的情況。 傳統GSS 的演算法如下述：

圖4.1 Golen section search(GSS)示意圖

如圖4.1，點 a 及點 c 為已知 2 點，其誤差值為 及 ，我們要用GSS 在 ac 之間找到一極小值。步驟如下： a F F_c Step 1： c b2 b1 a α F(α)

1 2 1 1 2 2 ( ) (1 )( ) ( ) ( ) b a c a b a c a F F b F F b γ γ = + − = + − − = = Step 2：

While (c-a) >ε(a+c)

Case 1： If (F₁≤F₂) set { 2 2 1 2 1 1 2 1 1 ( ) ( ) c b b b F F b a b a F F b γ = = = = + − = } Case 2： else set { 1 1 2 1 2 2 1 2 2 (1 )( ) ( ) a b b b F F b a c b F F b γ = = = = + − − = } Step 3：Set α =(a c+ ) /2，為最後求得的解。 以上為傳統GSS 的做法，然而，目前複雜的值域空間中任何單一維度都是非 unimodal，我們需對 GSS 稍做修改來符合目前的情況，主要想法有二：(1)原本 在比較內部2 點的誤差值大小後，下一階段沒有選到的範圍會被捨棄不再細看， 但在非unimodal 的情況下無法保證沒有選到的範圍絕對不會有更小誤差值的 點，如圖4.2，所以我們在兼顧準確度及效率下從要捨棄的範圍中多選出一點來 觀察;(2)如果 2 端點的誤差值比內部選出 2 點的誤差值還小，則端點到選擇點之 間的範圍可能會有區域極小值，所以這段範圍要紀錄下來，之後執行GSS。

圖4.2(a) GSS 改進步驟(1) 改進想法(1)的做法如下： 在傳統做法的Step2 中，把目前搜尋的區域分成 2 個部份：( , )a b₂ 及( , )b c₁ ，由 1 b 或 決定接下來搜尋的區域要選擇 或是 ，若選擇 範圍，則 就不再細看並捨棄掉。但如圖4.2，( , 範圍的最佳解就在 這一段， 所以新的方法中，我們在 取一點 2 b ( , )a b₂ ( , )b c₁ ( , )a b₂ 2 ( , )b c a c) ( , )b c₂ 2 ( , )b c b₂'= +b₁ γ(c b )− ₂ ，計算 ，結果分 成2 個 case： 2 ( ' ) F b Case 1：如果 ，則認為( , 沒有細看的必要，照原來做法繼續對 執行GSS。 2 ( ' ) ( ) F b ≥F b₁ b c) 1 a b ) b c) 2 2 ( , )a b Case 2：如果 ，則認為 可能有區域極值，分別對( , ，( , 執行GSS。 2 ( ' ) ( ) F b ≤F b ( , )b c_{2 2 2} b2 b1 c F(α) a b2 ’ b2 b1 c F(α) a α α 圖4.2(b) GSS 的改進步驟(2) 改進想法(2)的做法如下： 在傳統做法中，因為已預先假設函數為unimodal 且極小值 ( , 之間，所以一開 始 及 不用考慮，但在非unimodal 的情況下，一開始 及 的值 也要和 及 一起比較。以圖4.3 為例，在傳統做法中 為選擇要細 ) a c ( ) F a F c( ) F a( ) F c( ) 1 ( ) F b F b( )₂ ( , )b c₁

Case 1：如果F a( )≥F b( )₂ ，則照原來做法繼續對( , )b c₁ 執行GSS。 Case 2：如果 ，則認為 可能有區域極值，分別對 ， 執行GSS。 2 ( ) ( ) F a ≤F b ( , )a b₁ ( , )a b₁ ( , )b c₁ 改進前，GSS 的複雜度計算如下： 假設為±4 度的範圍，精確度到 0.01 度，則800*0.618n ≤ ， 1 計算次數 = iteration 數 * 每個 iteration 的計算次數 = n * 1 = 14(每個維度) 改進後，GSS 的複雜度計算如下： 假設為±4 度的範圍，精確度到 0.01 度，則800*0.618n ≤ ， 1 計算次數 = iteration 數 * 每個 iteration 的計算次數 +後來增加的步驟 = n * 1 + 3 = 17(每個維度)

4.2 5 度空間最佳解搜尋法

由圖 4.3 可稍微了解整個解空間的複雜，其中會有許多的區域極值(local trap)，一般的搜尋法限於各種因素很難找到最佳解，所以我們採用 5 度空間的最 佳化搜尋，另外，由於每一維度的error surface 分佈並不一致，在某些維度上的 分佈較單純，甚至近乎unimodal function，較適合 GSS 的使用；某些維度的分佈 則過於複雜，無法順利地使用GSS。因此接下來將針對五度空間該以如何的方式 進行切割做相關的探討。 由於2D 空間的 T，和 3D 空間的 R 之間關係較獨立，所以可以分開計算。但 是在T 的 2D 空間中就有兩種可能的先後執行順序，先固定 ψT 再求 θT 的區域 最小值，或是先固定 θT 再求 ψT 的區域最小值。另一方面，在 3D 空間上的 R 則有六種不同的先後執行順序，與T 結合則有 12 種可能。

經過觀察，在eigenvalue 較小的維度上，由 error surface 可看到等高線較密 而error 的變化較大，比較可能有唯一的 local minimum，linear search 也較容易 找到最低點，所以先搜尋eigenvalue 較小的維度能得到較正確的 value 值，並可 以提供給之後的搜尋的維度，減少因error surface 過於複雜導致 GSS 誤判，進而 無法準確找到最佳解。

舉例說明，當ψR 的 eigenvalue 最小，而 θR 及 ωR 的 eigenvalue 差不多時， 在ωR 與 θR 所構成的 error surface(如圖 4.3 所示)，可以看出其 error surface 較為 複雜。而在由ψR 與 θR 所構成的 error surface(如圖 4.4 所示)，可以看出 error surface 較為單純，較容易使用GSS 來計算最佳解。

圖4.3 error surface(水平軸為 ωR，垂直軸為 θR)

求解較有利，三個參數空間的參數設定如圖4.5，1D rotation eigen-space 設為

last

R ，2D rotation eigen-space 設為(Rhor,R )，2D translation eigen-space 設為ver

(Thor,Tver)。

R-3D

圖4.5 三個參數空間的參數設定為Rlast、(Rhor,R )、(ver Thor,Tver)

實際執行時，如圖4.6，Rlast維度最後搜尋，因此在Rlast維度上以GSS 找出的

最佳解即為全域最佳解，而要執行GSS 則必須能計算出Rlast維度上任一點的誤

差值。

last

R 維度上任一點的誤差值即是固定R_last維度為某一值，在2D rotation eigen space 的Rhor維度上以GSS 找出的最小值，舉例說明，Rlast= r3這一點的誤差值為

在R_last= r3的情形下，用GSS 尋找Rhor維度上的最小值，同樣的道理，要執行GSS

則必須能計算出R_hor維度上任一點的誤差值，如圖4.6 上 r2這一點的誤差值即為

固定R_last= r3，Rhor= r2，並在R 維度上以 GSS 找出最小值。 ver

同理，如圖4.6R 維度上一點 rver 1的誤差值為固定R 的 3D 旋轉空間為(r1,r2 ,r3) 的情形下，在2D translation eigen-space 中找出的最小值。依序搜尋，在 維度 上t2點的誤差值為固定4D 的角度(r1 ,r2 , r3,t2)，並在 維度上以GSS 找出的最 小值。 hor T ver T

從前面定義可知， 的維度eigen value 較小，其 error surface 上的等高線較 密而誤差值變化較大，比較可能有唯一的極小值，找出的值比較準確可信，也因 此在 上做GSS 搜尋時所依據每一點的誤差值是可靠的，以此層遞上去，最後 才對 ver T hor T last R 維度執行GSS 搜尋。

T-2D

R-2D

T

T

v

R

R

_v veerr veerr

R

R

llaasstt

T

T

_hhoorr

R

R

h_hoorr

圖4.6 全域搜尋法執行順序流程 全域搜尋的詳細演算法如下說明： begin initialize

[

]

⎪ ⎩ ⎪ ⎨ ⎧ + − = − ← − ← − ← − ← ver ver hor ver hor last T ver T ver ver T hor hor R ver ver R hor hor R last last Range T , Range T T , Range T T , Range R R , Range R R , Range R R

for Rlast ←Rlast +1

for Rhor ←Rhor +1

for Rver ←Rver +1

for Thor ←Thor +1

do GSS alone Tver axis with error function

) ,

(R_last,R_hor,R_ver T_hor,T_ver Err

find minimum Tver_min

until GSS stopping criterion until hor T hor hor T Range T = +

do GSS alone Thor axis with error function

) (R_last,R_hor,R_ver,T_hor,T_{ver_min} Err

find minimum Thor_min

until GSS stopping criterion

+ =

) (R_last,R_hor,R_ver,T_{hor_min},T_{ver_min} Err

find minimum Rver_min

until GSS stopping criterion until hor R hor hor R Range R = +

do GSS alone Rhor axis with error function Err(Rlast,Rhor,Rver_min,Thor_min,Tver_min)

find minimum R_{hor_min}

until GSS stopping criterion until last R last last R Range R = +

do GSS alone R_last axis with error function Err(R_last,R_{hor_min},R_{ver_min},T_{hor_min},T_{ver_min})

find minimum Rlast_min

until GSS stopping criterion

return R_{last_min}

第 5 章實驗結果

在本論文的研究中，對最後求出來的全域解結果影響最大的因素為error surface 的變化，究其根源即是一開始對應點的雜訊大小，所以實驗分成 3 個方 面，前2 項為模擬實驗，第 3 項為實拍影像的結果，而模擬實驗又分成半雜訊實 驗及全雜訊實驗。 半雜訊實驗中，第2 張影像中至少有一半的對應點完全沒有雜訊，我們可以 確實掌握正解的位置(ground truth)，以此來驗證最佳解搜尋法所求出的 5D 解和 正解在解空間中是否相近，而全雜訊實驗則是模擬真實的情況，在第2 張影像中 所有的對應點都會加入雜訊，但因為正解的確實位置難以得知，最後以全域最佳 解重建3D 物體，觀察重建物體的邊長及角度來和原物體比較，以得知解的好壞。 實驗以一個立方體(cube)當做 3D 物體，對應點編號定義如下圖 5.1： 圖5.1 在立方塊(cube)可見的 3 個面上的對應點編號 模擬實驗的雜訊假設為高斯分布(Gaussian distribution)，雜訊大小參考立方體 (cube)在 2D 的投影邊長，將立方體(cube)在 2D 投影的最大邊長的 n%設定為 3σ， 以下圖5.2 為例，cube 在 2D 上最大邊長為 p p_{1 13}(= 289 pixels)，取其 3%(=8.659 pixels)為雜訊高斯分布的 3σ(mean = 0, 標準差= σ)。

圖5.2 立方體(cube)在 2D 上的影像 驗證結果是否可靠的步驟如下： 1. 以 19 組對應特徵點及外部參數解重建 3D 點(見 2.2 節) 2. 如圖 5.1，設定立方體可見的 9 個邊長為p p 、_{1 3} p p 、_{1 13} p p 、_{1 7} p p 、_{3 9} p p 、_{3 15} 7 9 p p 、p p 、_{13 15} p p 、_{15 19} p p ，12 個角度為_{9 19} ∠(JJJJK JJJJKP P_{1 3},P P_{1 7}) ∠ JJJJK JJJJKP P P P ∠ JJJJK JJJJK P PJJJJK JJJJJK P P 、 P P )、 ) P P 、∠ P P )、∠ P P ) 1 3 1 13 ( , 7 9 7 1 ( , ( _{15 3}, _{15 19} (JJJJK JJJJJK_{19 9}, _{19 15} 、∠(JJJJK JJJJKP P_{3 1},P P_{3 15}) P PJJJJK JJJJK P P ∠ JJJJK JJJJJK P PJJJJK JJJJK P P 、∠ P P)、 ) P P 、∠ _{3 9} P P_{3 1})、∠ P P ) 3 15 3 9 ( , 15 3 15 13 ( , ( , (JJJJK JJJJK_{9 3}, _{9 19} 、∠(JJJJK JJJJKP P_{9 7},P P_{9 3})、∠(JJJJJK JJJJP P_{13 15},P P_{13 1}K)。 3. 將重建 cube 的邊長放大，使 9 個邊長的平均值和 ground truth cube 邊長平均

相同(=40cm) 4. 計算重建 3D 後的 9 個邊長長度和 ground truth 對應邊長差的絕對值的平均 及標準差。 5. 計算重建 3D 後的 12 個角度和 ground truth 對應角度(=90o)差的絕對值的平 均及標準差。 實驗使用10 張模擬影像如下：

image 1 (cube 最大邊長=)

image 2 image3

image 4 image 5 image 6

image 7 image 8 image 9

image 10

5.1 半雜訊實驗

實驗的主要步驟如下： (1) 選取要實驗的 2 張影像 (2) 參考第 2 張影像中 cube 的最大邊長(取 3%為 3σ 的雜訊)，分別產生 9 個 x、y 座標的高斯雜訊(mean = 0, 標準差= σ)，並隨機加入第 2 張影像中的 9 個對應 點。 (3) 在 19 組對應點中隨機選取 10 組對應點來計算出外部參數解，以 LMedS 的方 法產生4000 組相機外部參數解當作初解。

(4) 分析有前 50 個最小誤差值的解，建立其 Minimum Spanning Tree(MST) (5) 對 MST 做分群(Clustering)，求各群 Covariance matrix 橢圓體之 eigen values 及

eigen vectors

(6) 對各群橢圓體細切，得到初解的改進(refinement)，直至橢圓體退化(至少一個 eigen value 為 0)。

(7) 啟用 Golden section search(GSS)搜尋相機外部參數的最佳解：(a)轉換 5D Camera parameter space 為 3D Rotation space 與 2D translation space，(b)進而 再用橢圓體eigen values 及 eignen vectors 定義三個參數空間:1D rotation eigenspace +2D rotation space +2D translation space，(c)使用 Golden section line search 連環搜尋相機外部參數的最佳解 (8) 以找到的最佳外部參數解重新建立 3D 點，計算重建 3D 後 cube 的邊長長度 和ground truth 對應邊長差的絕對值的平均及標準。 以下為實驗詳細過程及結果：

實驗一：使用合成影像資料，並於部份特徵對應點加入雜訊。

輸入兩張解析度為1024(水平)x768(垂直)的影像 Imgae1 及 Image2，影像中各有 19 點特徵對應點，如下圖所示。其中對 View2 的特徵點中的 9 點加入雜訊。

View 1 (Imgae1) View 2(Image2) View 1 及 View 2 之間正確的外部參數設定如下：

R T

0.6597 0.416 0.6258 -0.8894 -0.4356 -0.4669 0.7696 0.1203 將正確的外部參數轉成以5D 的表示法：

ψR θR ω ψT θT

39.6843 101.4559 -46.9876 83.0882 296.3790

在此實驗中，設定雜訊為normal distribution (mean = 0, 標準差= σ)，雜訊的 大小是以cube 在影像中投影的最大邊長(約 300 pixels)的 3%為三個標準差 (3σ 約等於 9.8324 個 pixels)。View 2 的 19 個特徵點當中加入雜訊後的結果如下表 所示。 特徵點編號在X 軸方向的雜訊大小在 X 軸方向的雜訊大小 1 -3.72704 -4.4568 2 0 0 3 0 0 4 0 0 5 0.424926 2.168381 6 -0.86687 -4.00717 7 0 0 8 0.037143 -2.13098 9 -3.26928 4.424792 10 3.759026 -2.25978 11 0 0 12 0 0 13 0 0 14 0 0 15 -1.092 -2.78659 16 -1.80507 -2.79693 17 0 0 18 -2.8215 -3.96807 19 0 0 利用LMEDS 計算後的最佳 50 個的解資料如下：

n distance 1 37.5531 99.8379 -45.9876 84.8507 295.4805 3.7591 3.4748 2 36.0189 98.1725 -44.4523 84.5594 294.7667 3.9509 5.9504 3 34.2524 98.9948 -44.8329 88.6701 295.9604 4.1066 8.4581 4 34.2524 98.9948 -44.8329 88.6701 295.9604 4.1066 8.4581 5 41.2622 101.9975 -46.919 82.9204 296.2237 4.2725 1.6853 6 34.297 99.0591 -44.8852 88.7017 295.9782 4.4058 8.4178 7 38.4302 100.2002 -46.0319 82.9454 295.7285 4.7959 2.1229 8 42.6694 102.7989 -48.2871 81.3871 296.8592 4.9876 3.9405 9 37.9832 100.8055 -46.1263 85.7353 296.3444 5.1658 3.3267 10 37.9856 101.4448 -46.9422 85.262 296.6763 5.2509 2.7752 11 39.0285 102.4825 -47.6271 85.1247 297.2173 6.3983 2.5968 12 37.1346 98.3514 -44.7518 83.4899 294.4052 6.4641 5.0195 13 36.7468 100.8538 -46.1187 87.6357 296.6693 6.8248 5.5237 14 40.9361 101.3863 -46.3703 83.1226 295.6708 6.83 1.5671 15 37.9234 100.0645 -45.4888 86.5125 295.7134 7.3285 4.4105 16 41.836 99.0254 -45.587 81.7029 293.4182 7.704 4.815 17 33.7608 101.5377 -45.7873 84.3183 298.3971 7.7908 6.4901 18 37.7916 100.6108 -46.0085 81.8255 296.3045 7.9153 2.6182 19 40.2755 101.4804 -46.5675 79.43 296.0941 7.9743 3.7403 20 39.034 102.1562 -47.1281 85.59 296.9954 8.0272 2.7517 21 43.0691 105.5251 -49.7908 81.4445 299.102 8.0621 6.7815 22 43.7877 101.4997 -47.5225 79.625 295.1069 8.1503 5.5442 23 35.3524 101.5015 -45.9431 82.5251 297.7128 8.26 4.6856 24 38.7429 98.4797 -45.5841 83.59 294.1984 8.3243 4.0891 25 38.0458 103.2024 -47.4385 84.4089 298.4583 8.4548 3.465 26 40.3055 101.4212 -46.8773 79.5141 296.0962 8.461 3.6405 27 40.4335 102.1987 -47.1764 84.4869 296.8915 8.8006 1.8352 28 39.1469 98.6854 -45.8275 84.2113 294.2206 9.0006 3.9026 29 39.2479 102.6772 -47.6295 85.8929 297.5577 9.1633 3.369 30 38.4597 100.6116 -46.405 82.8523 296.2495 9.2038 1.62 31 41.1537 100.7348 -46.0954 79.4781 295.3784 9.3482 4.1844 32 40.4995 100.0858 -45.2736 79.661 294.933 9.3832 4.395

33 37.8155 99.4717 -45.4176 81.8303 295.3497 9.4327 3.5406 34 40.3305 101.4924 -46.6713 81.1709 296.5068 9.617 2.0521 35 39.7327 99.2624 -46.3051 84.7055 294.5661 9.799 3.344 36 35.7458 101.6899 -46.2794 88.629 297.6977 9.9243 6.9648 37 43.52 103.6721 -48.8987 80.9707 297.532 10.1068 5.3934 38 38.1212 100.6823 -45.5551 87.3407 295.9162 10.2978 4.8365 39 36.1849 100.4142 -46.2364 86.0084 296.2496 10.3609 4.7371 40 41.4011 101.8832 -47.0462 80.0084 296.5861 10.498 3.5582 41 38.611 102.1176 -46.8123 86.9478 297.2641 10.591 4.1594 42 40.6352 101.7184 -46.7316 79.3856 296.346 10.6395 3.8404 43 39.9681 102.0977 -47.3653 85.8105 296.649 10.7301 2.8494 44 42.2912 102.8843 -48.4694 84.4834 297.0423 10.7593 3.6631 45 40.2495 99.0221 -46.7412 84.3908 293.9635 10.804 3.7196 46 39.6722 98.4493 -46.0531 83.5563 293.5848 10.8583 4.2356 47 39.4335 102.2773 -47.2754 83.3435 297.4144 10.9246 1.3991 48 37.9033 103.3648 -47.8614 86.2076 298.3001 11.0188 4.5827 49 39.7422 102.4623 -47.7064 86.5921 296.9656 11.0297 3.7622 50 40.2511 101.8242 -47.812 81.3859 296.2698 11.1722 2.0115 其中的median error 定義如下：

median error = median of 19 對應特徵點的 error

LMEDS 初始解的最小值為 3.7591，對這 50 點進行細切割，先分別對 R, T 進行 計算橢圓體的長短軸相關資訊。如下所示： R ψR θR ωR mean 39.7412 101.3659 -46.6089 eigenvalue _{1.1996 0.9086 0.2333} 軸1 軸2 軸3 eigenvector -0.6952 -0.6875 0.21 0.5798 -0.3635 0.7292 -0.425 0.6287 0.6513 T ψT θT mean 83.7217 296.6401 eigenvalue _{2.951 0.5204} 軸1 軸2 eigenvector -0.9922 -0.1244 -0.1244 0.9922 R 橢圓體最長軸的半長 = Cmax eigenvaluemax _{= 2.8900}

取最大切割單位 = ( 2 / 3 * half of 最長軸 * 1/2 ) = 0.9633 T 橢圓最長軸的半長 = Cmax eigenvaluemax _{= 6.2600}

取最大切割單位 = ( 2 / 3 * half of 最長軸 * 1/2 ) = 2.0867 其它軸切割單位依cluster 的 eigenvalue 的比例設定 Æ對 Cluster 1 R 細切的切割間隔設為 = (0.9633 0.7296 0.1873) T 細切的切割間隔設為 = (2.0867 0.3680) 對橢圓體細切後，可得到一個初步的refinement 結果 ψR θR ω ψT θT medi err 38.9598 101.2675 -46.6782 83.5729 296.3703 2.1176 以此解為中心，在每一個dimension ± 4 度的範圍內用新方法找出一個最佳近似 解最佳近似解. 搜尋過程中所觀察之 2D error surface contour (水平軸為 ωR , 垂 直軸為θR)如下圖所示：

用GSS 方法找到的最佳解的離散值如下： ψR θR ωR ψT θT medi err 39.6869 101.4677 -46.9934 83.0828 296.3891 0.0676 對照以下所附正解而言，各相機角度參數誤差約在0.01 度， media error = 0.0676 pixels. ψR θR ω ψT θT 39.6843 101.4559 -46.9876 83.0882 296.3790 以找到的最佳解的離散值所重建後的3D 結果如下： Error Function 的函數值 (pixels) 重建3D 後 的邊長長度 和ground truth 對應邊 長相差絕對 值的平均 (cm) 重建3D 後 的邊長長度 和ground truth 對應邊 長相差絕對 值的標準差 重建3D 後 的角度和 ground truth 對應角度相 差絕對值的 平均(度) 重建3D 後的 角度和ground truth 對應角 度相差絕對 值的標準差

初步改進後

的解 1.9021 0.3281 0.2656 0.6307 0.3877

最佳解

實驗二：使用合成影像資料，並於部份特徵對應點加入雜訊。

輸入兩張解析度為1024x768 的影像 Image1 及 Image3，影像中各有 19 點特徵對 應點，如下圖所示。其中對View2 的特徵點中的 9 點加入雜訊。

View 1(Image1) View 2(Image3) View 1 及 View 2 之間正確的外部參數設定如下： R T 0.8627 -0.1846 -0.4708 0.6826 0.3878 0.839 0.3816 -0.6314 0.3245 -0.5118 0.7955 0.368 將正確的外部參數轉成以5D 的表示法： ψR θR ω ψT θT 37.3003 140.9711 -122.3799 68.4080 317.2328 在此實驗中，設定雜訊為normal distribution(mean = 0, 標準差= σ)，雜訊的 大小是以cube 在影像中投影的最大邊長的 3%為三個標準差 (3σ 約等於 9.452 個pixels)。View 2 的 19 個特徵點當中加入雜訊後的結果如下表所示。 特徵點編號在X 軸方向的雜訊大小在 X 軸方向的雜訊大小 1 0 0 2 -3.72704 -4.4568 3 0 0 4 0.424926 2.168381 5 0 0 6 -0.86687 -4.00717 7 0 0

10 -3.26928 4.424792 11 0 0 12 3.759026 -2.25978 13 0 0 14 -1.092 -2.78659 15 0 0 16 -1.80507 -2.79693 17 0 0 18 -2.8215 -3.96807 19 0 0

利用LMEDS 計算後的 median error 最小的前 50 個解資料如下：

num ψR θR ω ψT θT median error 與正解的 Euclidea n distance 1 29.6682 141.4966 -121.405 74.1991 317.2275 6.6906 9.6442 2 32.7053 142.2839 -123.827 71.7361 317.8237 6.8135 6.0297 3 28.8298 143.0459 -124.37 69.4278 317.8257 7.8483 9.0224 4 27.6632 137.5118 -117.67 72.2759 313.0297 9.4849 12.6352 5 42.9453 141.0924 -123.253 69.5353 317.7244 9.4916 5.8443 6 36.0949 142.4432 -123.678 66.3573 317.6953 10.5248 3.1184 7 36.1698 141.7338 -122.765 70.313 318.3048 10.5706 2.605 8 41.2414 139.2393 -121.476 67.8386 315.594 13.5777 4.7286 9 39.5726 139.5241 -121.497 65.219 315.1695 14.2116 4.7397 10 40.441 139.1832 -121.745 72.2206 316.3319 16.9895 5.3676 11 26.8339 141.8992 -122.442 65.9585 316.4764 17.4494 10.8159 12 35.952 140.8474 -122.268 70.2366 316.5295 18.3644 2.3842 13 36.1866 139.7995 -121.099 68.0759 316.4651 19.4008 2.2254 14 40.1426 138.9176 -121.799 73.4008 315.8542 20.5001 6.2819 15 41.1813 138.0518 -120.737 72.0803 315.7656 21.2277 6.4749 16 38.5412 140.4795 -122.203 72.246 317.9028 23.5363 4.1221 17 42.0116 138.916 -121.46 59.9859 313.762 23.9649 10.4997 18 31.2562 140.4504 -123.067 82.2089 316.9913 26.2528 15.0929 19 27.7324 138.1323 -121.695 84.0006 314.6408 26.7347 18.7062