本論文要執行的工作包括建立模擬及實拍影像，以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，求外 由上式可知，若將 令為1，僅需八組對應點即可利用 least-square method 來計 算 ，計算方法如下：
F33
F
11
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：
以 Forigin計算 中對應點在 上的epipolar line 到 上相對應點的距離
加上 中對應點在 上的epipolar line 到 上相對應點的距離，計算出
View1 View2 View2
View2 View1 View1
n 個值並取中位數，即為 n 個值排序後的第⎡⎢n/ 2⎤⎥個值，公式如(2.3)~(2.5)；
影矩陣(projection matrix)，如下(2.6)：
Emodified View1 View2
[ ]
3 1 3 2
一開始外部參數為R(rotation matrix)及 T(translation matrix)，一般而言，R 可 以用三個角度：對x 軸旋轉角度 θx、對 y 軸旋轉角度 θy、對z 軸旋轉角度 θz來
而 T real=[Tx , Ty ,Tz]T表示2 台 camera 實際上在世界座標裡的相對位置，但以 2
θ θ ϕ
⎡⎣ c_final c_final c R= x y z ⎤⎦
用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）是將所有的資料作排序後，取中間的數值，並找出最小的那一組
因應而生。
我們先用LMedS 的方法在 2 張已校正影像的 n 組對應點中，重複 t 次隨機挑 選k 組(至少八組)對應點，經由計算 Fundamental matrix，加上已知的內部參數，
可得到Essential matrix，再求出 t 個外部參數解及其誤差值，如第 2 章所述，並 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 = eigenvaluemax cmax (3.2) 其它軸半長的長度依 eigenvalue的比例可得。
其 中 cmax =max{ ci = (xi−mean)TΣ−1(xi−mean) ,i= ~ }1 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組新解的誤差值。
初步細切後，同樣對誤差值較小的前50 名做分群，不過為盡量防止橢圓體退
化而且延橢圓體形狀來分群，所以之後在建立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 之間找到一極小值。步驟如下：
Fa Fc
Step 1：
c b2
b1
a α
F(α)
1 unimodal，我們需對 GSS 稍做修改來符合目前的情況，主要想法有二：(1)原本 在比較內部2 點的誤差值大小後，下一階段沒有選到的範圍會被捨棄不再細看，
但在非unimodal 的情況下無法保證沒有選到的範圍絕對不會有更小誤差值的 點，如圖4.2，所以我們在兼顧準確度及效率下從要捨棄的範圍中多選出一點來 觀察;(2)如果 2 端點的誤差值比內部選出 2 點的誤差值還小，則端點到選擇點之 間的範圍可能會有區域極小值，所以這段範圍要紀錄下來，之後執行GSS。
圖4.2(a) GSS 改進步驟(1)
Case 1：如果F a( )≥F b( )2 ，則照原來做法繼續對( , )b c1 執行GSS。
Case 2：如果 ，則認為 可能有區域極值，分別對 ，
執行GSS。
( ) ( )2
F a ≤F b ( , )a b1 ( , )a b1 ( , )b c1
改進前，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。因此接下來將針對五度空間該以如何的方式
由圖 4.3 可稍微了解整個解空間的複雜，其中會有許多的區域極值(local trap)，一般的搜尋法限於各種因素很難找到最佳解，所以我們採用 5 度空間的最 佳化搜尋，另外，由於每一維度的error surface 分佈並不一致，在某些維度上的 分佈較單純，甚至近乎unimodal function，較適合 GSS 的使用；某些維度的分佈 則過於複雜，無法順利地使用GSS。因此接下來將針對五度空間該以如何的方式