1、1第八章第八章圖形識別圖形識別、匹配與三維影像重匹配與三維影像重建建28.1 前言n8.2 統計圖形識別n8.3 影像間的匹配對應 n8.4匹配演算法原理n8.5三維影像重建n8.6 二維影像的深度計算n8.4.1 動態規劃式的BSSC解法n8.4.2 KMP演算法n8.6.1 稠密式視差估測n8.6.2 相機校正n8.3.1 Harris角點偵測法n8.3.2 SIFT關鍵點偵測法n8.3.3 點集合匹配法38.2 統計圖形識別n貝氏決策理論 假設有二類木頭，A和B，佔P(A)的比例而B佔P(B)的比例，。已知 。利用木頭的紋理 X 來評估該木頭的種類。()()1P AP B()()P AP
2、 B圖8.2.1 P(X|A)和P(X|B)的分佈圖 當 時，我們說該木頭為A類，仍會犯誤判的風險，畢竟 仍有機率值。=XX(|)(|)P X AP X B(|)P X B4)()()|()()()|(XPAPAXPXPXAPXAP我們有興趣的是給一個X值，該木頭屬於A或B的機率為何？依據貝氏法則，透過貝氏法則，這個事後機率就可由事前機率P(A)、P(B)、和 求得。)()|()()|()(BPBXPAPAXPXP此處。)|(XAP)|(AXP)|(BXP)|(XAP)|(XBP機率值XX圖8.2.2 P(A|X)和P(B|X)的分佈圖 當 ，這時可判斷該木頭為 A，畢竟冒的風險較低，也就是
3、。去掉 和 的 項，當 時，我們判斷該木頭為A。XX=)|()|(XBPXAP)|()|(XBPXerrorP)|(XAP)|(XBP)(XP)()|()()|(BPBXPAPAXP5將紋理由一維擴充到 d 維而將木頭的種類由2種擴充到 t 種。令第 i 類的識別器為 ，此處紋理向量 而 代表第 i 類木頭，。作用單調遞增函數於上，得ti 1()log(|)log()iiig XP X TP Tn識別器()(|)()iiig XP X T P T12(,)dXx xxiT圖8.2.3 識別器示意圖()tgX如果 為最大值，該木頭分類為 Tj ()jgX6視窗在影像上的移動，可得到強度變化情形：
4、n(1)平面：往任何方向移動僅造成小變化n(2)含一條邊：與邊平行的變化量小；反之則大n(3)含角點或獨立點：往任何方向變化皆大8.3 影像間的匹配對應8.3.1 Harris 角點偵測法(1)(2)(3)7令 利用高斯函數來平滑雜訊的影響，令綜合灰階梯度變化之影響可表示為 xx(1,0,1)ff yy(1,0,1)tff 2xAG 2yBG xyCG 222(y)ACxAxByCxyxCBy()()txy M xyE8n函數E是一種局部自我關聯的函數，矩陣M就是函數E的代表。矩陣M的兩個特徵值代表：和 皆很小：代表視窗內為平滑區 和 中，一大一小：代表含一邊的區域 和 皆很大：代表含角點的區
5、域111222圖8.3.1.1 矩陣M的特徵值所代表的意義9n影響值 R=det(M)-k(trace(M)2決定是否小區域內有角點，判斷的準則為 且 ：代表平滑區 ：代表含一邊的區域 ：代表有角點的區域利用以上方法可以將所有角點找出來。0R 0R 0R(a)原始影像(b)找出的角點集0R 10n將影像進行縮減取樣 n 次，例如n=2可分成3個影像尺度，進一步利用 m 個不同標準差的高斯函數，例如m=4，高斯函數標準差分別是 ，將不同高斯函數與3個影像尺度作迴積運算：接著對兩相鄰的尺度做DOG:23,kkk(,)(,)*(,)ttL x y kG x y kI x y1(,)(,)(,)*(,
6、)tttD x y kG x y kG x y kI x y1(,)(,)ttL x y kL x y k 8.3.2 SIFT關鍵點偵測法11圖8.3.2.1 DOG金字塔12考慮每個像素與周圍8個像素及上下層同位置周圍9個像素做比較，如果該像素為極值，則設為候選關鍵點。圖8.3.2.2 DOG金字塔找極值13以每個候選關鍵 X 點當原點，算出 其中 的三維座標，令座標時有極值點D()。代表一個極值點，如果內的三個座標值都小於0.5;若大於0.5則依照該方向來內插出新的候選關鍵整數點，繼續用新的候選關鍵點計算二階泰勒展開式求值。若 ，代表此區域的對比較低，應移除此關鍵點。22(0)1(0)(
7、)(0)2TTDDD XDXXXXX(,)TXx y22(0)(0)TDDXX0122(0)(0)TDDXX|()|0.03D14 除掉邊上的候選關鍵點，先求出赫斯矩陣 H的特徵值 和 可計算如下:可得xxxyxxxyyxyyxyyyDDDDHDDDD12det()detxxxyxyyyDDHIDD22()()0 xxyyxxyyxyDDD DD12trace()xxyyDDH21 2det()xxyyyyD DDH15令 且 ，得到得判斷式令 ，若判斷式成立，移除邊上的候選關鍵點。1212222212222122()()(trace()(1)det()HH 22(trace()(1)det(
8、)HH1016以剩餘關鍵點為中心的區塊，計算區塊內每個像素的梯度值 與梯度方向關鍵點的旋轉不變性(Rotation Invariant)。(,)m x y(,)x y22(,)(1,)(1,)(,1)(,1)m x yL xyL xyL x yL x y1(,)tan(,1)(,1)/(1,)(1,)x yL x yL x yL xyL xy17圖8.3.2.3 決定關鍵點的特徵向量18圖8.3.2.4 SIFT流程圖19(a)輸入的腦血管影像(b)SIFT演算法選取到的關鍵點圖8.3.2.520n另一種作法是尋找轉移矩陣 T，使 ToF 與 有最小誤差。令 ，代表 中的 與 中的 之距離n引
9、入最大可能(Maximum-likelihood)的觀念，得其中 p(d)為一常態分佈(Normal Distribution)。n匹配的精神就是找一個 T 使得上式有最大值。F()id T1im Fifif121(),(),()|)()mmiip d Td TdTTp d T121()ln(),(),()|)ln()mmiiL Tp d Td TdTTp d T8.3.3 點集合匹配法F218.4 匹配演算法原理 8.4.1 動態規劃式的BSSC解法nBSSC(Banded String-to-String Correction)之目的為在樣本和正本之間找出最大的匹配點與匹配點的位置。圖8.
10、4.1.1 匹配結果123458910PP PP PPPP1 23458 9 10TT TT TTT T將和想像成特徵值。因為視差的關係，只能與 範圍內的特徵匹配到。jPiTjPiT樣本=正本=給一個例子如下：圖8.4.1.1的搜尋範圍就是黑色粗邊框住的範圍而得到的最佳匹配就是鋸齒形的路徑上之黑圓點集。22n匹配算子 取代算子R 刪除算子D D(a)=的花費定為1插入算子I I()=a的花費定為1 a b，R(a)=b的花費定為1a=b，R(a)=b的花費定為0 圖8.4.1.1一共含有下列十四個運算(1)I()=P1 (2)R(T1)=P2 (3)R(T2)=P3 (4)D(T3)=(5)D
11、(T4)=(6)I()=P4 (7)R(T5)=P5 (8)R(T6)=P6 (9)D(T7)=(10)D(T8)=(11)I()=P7 (12)I()=P8 (13)R(T9)=P9 (14)R(T10)=P10 23n動態規劃核心式子 iTTTiT.1 21jPPPjP.1 21Edit0,k=k Editk,0=kEdit,min(Edit1,edit(,),ii jijTEdit,1edit(,),Edit1,1edit(,)jiji jPijT P表示將T1i轉換成P1j的花費。edit(Ti,Pj)表示R(Ti)=Pj的花費edit(Ti,)表示D(Ti)=的花費edit(,Pj)
12、表示I()=Pj的花費 nBSSC的時間複雜度為n為正本長度，為視差()n248.4.2 KMP演算法n樣本先進行事前處理 利用樣本中子字串(Substring)與前置字串(Prefix String)的吻合度，並記錄其吻合的長度於陣列 J 中。i1 2 3 4 5 6 7 8 9 10P i acaca a acacJ i 0 0 1 2 3 1 1 2 3 4P3=a=P1 J3=1P35=aca=P13 J5=3P710=acac=P14 J10=4nKMP字串匹配演算法T =cccacacaaacaccaa Ti Pi，。T413=P110。T5 P1，T6=P1 T68=P13 T9
13、 P4，J 4=2 P1J 4-1=P1=T9-J4+18=T8 陣列J 提供了一個跳躍的機制，讓匹配動作可一直往右前進。31 inKMP的時間複雜度為O(m+n)258.5 三維影像重建8.5.1 稠密式視差估測 n問題定義 稠密式視差估測是在L上的所有像素中和R上的所有像素中找到一個對應。n二階段式的分割與克服方法 PXYZOLXYZOR右影像(R)左影像(L)圖 8.5.1.1 稠密視差估測示意圖(Dense Disparity Estimation)26n第一階段的分割與克服方法1.利用Marr-Hildreth測邊法得到主要特徵像素集：22222212 1NNNNLLLLNNSSSS
14、NNN，L上第 N/2 列中的主要特徵集可表示為：R上第 N/2 列中的主要特徵集可表示為：22222212 1NNNNRRRRNNSSSSNNN，2.和 的匹配工作 2NLS2NRS N 222222LS,()min,NNNNNRmatchLRrightoccLRCvSwCSvSwCSvSw 22,NNleftoccLRCSvSw其中 22222222222222,1,1,1,1NNLRNNNNleftoccLRLROccNNNNrightoccLRLROccNNNNmatchLRLRSv SwCSvSwC SvSwCCSvSwC SvSwCCSvSwC SvSwC272/NLS2/NRS
15、12/NLS 22/NLS wSNL2/2/2/NNLmS 12/NRS 22/NRS vSNR2/2/2/NNRnS圖 8.5.1.2 和 的匹配2NLS2NRS其中 22212,*NNLRSv SwCll 2211NNLLlSvSv 2221NNRRlSwSw122211121 11mean112llNNLLRRpqISvpISwqll122222211121 111mean1mean2llNNLLRRpqISvpISwqll取代算子：28ABXCDX(a)右遮蔽花費(b)左遮蔽花費 圖 8.5.1.3 右遮蔽花費和左遮蔽花費示意圖 2212/2;CCOCkfthr0thr2*thr2/2
16、221thrf3.加入消去法則 若兩兩匹配位置的差，形成之序列為 ，則-2 的視差偏移不太正常，可以予以去除。假設在 列時的平均視差偏移序列為 。若新進來的視差偏移序列為 ，則 16 的視差偏移可予以去除。5,2,5,6,52Nk5,5,6,45,7,3,16,429n第二階段的分割與克服方法 1.在二段匹配的子區間中，再進一步找出個別的像素配對 1NN12/NLS12/NLS12/NLS12/NLS12/NLS12/NLS1ABC三維物體三維物體圖8.5.1.4 二子區間匹配 2.利用平均視差偏移量d，可將第二階段的工作轉成BSSC問題的解決上 22222,2,2,2,NNLRWNNLLRL
17、x yWCSvv SwwINx SvvyINx Swwy R算子：D算子：類似於前面定義的左遮蔽花費I 算子：類似於前面定義的右遮蔽花費W為事前定義好的小視窗12(1)NRS2(2)NRS2(3)NRSN30圖8.5.1.5 對應的BBSC搜尋空間31圖8.5.1.6 輸入的二張影像(a)影像L(b)影像R 圖8.5.1.7 得到的視差圖 圖8.5.1.8 重建後的三維五角大廈圖 328.5.2 相機校正n我們有三個座標系統，世界座標系統(World Coordinate System)、相機座標系統(Camera Coordinates System)和二維的影像系統(Image Syste
18、m)。根據三角比例關係，可得 其中 f 為焦距長。123456789 wxwywzxxrrrTyR yTRrrrTTzzrrrT txXfztyYfz33圖 8.5.2.1 相機座標系統和影像座標系統影像平面鏡心 理論上 為 投射到影像座標系統的理想位置。因為透鏡輻射效應，實際投射到的影像座標系統的位置應為 。滿足下式(,)ttX Y(,)x y(,)x y(,)aaXYaxtXEXaytYEY34上兩式中的誤差項 和 可表示成在光學透鏡的輻射失真(Radial Distortion)的影響下，內部參數滿足 ，我們令 ，只關心 的求解。前面所提的 為實數，但在數位影像的座標系統一般皆為整數座標
19、 。假設感應器(Sensor)的 方向共有 個感應器。是影像在x方向的解析度。令xEyE2412()xaEXrr2412()yaEYrr22aarXY21201(,)aaXY(,)IIX Y()X Y()xySSxNxxxxSddN35而 為影像的中心座標。和 的關係可表示如下其中 s 為待解的放大係數。已知 投影到影像座標系統的理想位置為 。進而得到同理可得且 (,)xyC C(,)aaXY(,)IIXYaIxxXXCsdaIyyYYCsd(,)x y(,),ttxyX Yffzz2121()()()IxxaxxIxxIxxxxXC dxfXEEzsXC dXC drssXdXdrss21y
20、yyd Yd Yrfz2222()()xxyIyyXdXdrd YCd Yss36 再利用三維世界座標系統和二維影像座標系統的關係可得旋轉矩陣 R 分別對 x 軸、y 軸和 z 軸達到任一角度的旋轉，可寫成下式 至此，已經六個外部參數(Extrinsic Parameters)分別為 、和 及五個內部參數(Intrinsic Parameters)分別為 、和 。接下來要利用數值的方法求解這些參數。21123789xxwwwxwwwzXdXdrrxr xr xTfssr xr xr xT45621789wwwyyywwwzr xr xr xTd Yd Yrfr xr xr xTcoscossi
21、ncossinsincoscossincoscoscoscossinsinsincossincoscossinsinsincoscoscosR xTyTyTfs1xCyC37先解出五個外部參數，在二維影像座標上取一經過原點的向量 ，且在相機座標上取一平行的向量 ，所以得令 則我們得到再利用五組世界座標 和五組底片上的真實座標 就可解出 、。(,)aaXY(,)x y456123(,)(,)00()()aaaaawwwyawwwxXYx yXyY xXr xr xr xTY rxr xr xT0wZ 111211415 yyawawaawawayxyyT rT rY xY yYX xX yXT
22、TT rT r(,)wwwxyz(,)aaXY11yT r12yT r1yxT T14yT r15yT r38 令則旋轉矩陣可改寫成利用 中每一行向量與列向量為單位長，可得下式其中 。111212114545yyyyrrT rT rCT rT rrr123456789yyyyrTrTrRrTrTrrrrR1 222212121 222245451 21 22222222 1 214251()1()1()1()(1)yyyyyyyyyrTrTTrrRrTrTTrrTrrTrrMT 22221245Mrrrr39 因為 R 為正交矩陣，任兩行的內積必為零。將 R 的第一行和第二行做內積，可以得到
23、此等式可解得 的兩個解，如下 前面曾推導過暫時將 和 設為零，得2421 52 4()10yyr rr rTMT 2yT1 2221 52 4221 52 44()2()yMMr rr rTr rr r45621789wwwyyywwwzr xr yr zTd Yd Yrfr xr yr zT1wZ4578wwyywwzr xr yTd Yfr xr yT40令 且 上式可改寫成再利用兩組以上的世界座標 和真實座標 解出 和 。然後以此 和 及令 為初始值來解線性方程式如此即可透過數值的解法求得 的逼近解。12 yyzftd Yt d YT145wwytr xr yT278wwtr xr y(
24、,)wwwxyz(,)aaXYfzTfzT102112yyztd Yd YrftT1(,)zf T418.6 二維影像的深度計算Epipolar線BC右邊相機焦點中心EFAEpipolar面物體上一點左影像右影像Stereo底線左邊相機焦點中心圖8.6.1 Epipolar面42範例1:已知兩部相機的焦距參見式(1.3.1)相同且均為 ，假設左邊的Epipolar線之長度為 而右邊的Epipolar線之長度為 。已知物件上的一點 ，試問 之值為何?解答:根據給定條件，圖8.6.2為對應的示意圖，圖中的B和C代表左相機和右相機。利用圖8.6.2中，和 ，透過兩相似三角形的邊比例關係，可得到flxrx(,)A x y zzBFGBDA CHICEA lxxzfrxxbzfrxCIHGFBlxfb(,)A x y zDE圖8.6.2 利用兩平行相機求z值43在上式中，我們假設兩不相機相距b，座標系統的原點設在B點。可進一步推得再利用 ，可得到深度為解答完畢lzxxf()rxzxbfflrzxzxbflrbfzxx
