Transcript 第八章圖形識別、匹配與三維影像重建

第八章
圖形識別、匹配與三維影像重
建
1
8.1 前言


8.2 統計圖形識別
8.3 影像間的匹配對應

8.3.1 Harris角點偵測法

8.3.3 點集合匹配法





8.3.2 SIFT關鍵點偵測法

8.4匹配演算法原理
8.4.1 動態規劃式的BSSC解法

8.4.2 KMP演算法
8.5三維影像重建
8.6.1 稠密式視差估測

8.6.2 相機校正
8.6 二維影像的深度計算
2
8.2 統計圖形識別

貝氏決策理論
假設有二類木頭，A和B，Ａ佔P(A)的比例而B佔P(B)的比
例， P( A)  P( B)  1 。已知 P( A)  P( B) 。利用木頭的紋理 X
來評估該木頭的種類。
當 X = X 時，P( X | A)  P( X | B) ，
我們說該木頭為A類，仍會犯誤
判的風險，畢竟 P( X | B) 仍有機
率值。
P( X | A)
機
率
值
P( X | B)
X
圖8.2.1 P(X|A)和P(X|B)的分佈圖
X
3
我們有興趣的是給一個X值，該木頭屬於A或B的機率為何？
依據貝氏法則，
P( A  X ) P( X | A) P( A)
P( A | X ) 

P( X )
P( X )
此處 P( X )  P( X | A) P( A)  P( X | B) P( B) 。
透過貝氏法則，P( A | X ) 這個事後機率就可由事前機率P(A)、P(B)、
P( X | A) 和 P( X | B) 求得。
當 X = X，P( A | X )  P( B | X )，
這時可判斷該木頭為 A，畢竟
冒的風險較低，也就是
P(error | X )  P( B | X ) 。去掉
P( A | X ) 和 P( B | X ) 的 P(X ) 項，
當 P( X | A) P( A)  P( X | B) P( B)
時，我們判斷該木頭為A。
P( A | X )
機
率
值
P( B | X )
X
圖8.2.2 P(A|X)和P(B|X)的分佈圖
X
4

識別器
將紋理由一維擴充到 d 維而將木頭的種類由2種擴充到 t 種。
令第 i 類的識別器為 gi ( X )  P( X | Ti ) P(Ti )，此處紋理向量
X  ( x1, x2 , , xd ) 而 Ti 代表第 i 類木頭， 1  i  t 。
作用單調遞增函數於上，得
gi ( X )  log P( X | Ti )  log P(Ti )
g1 ( X )
如果 g j ( X ) 為最大值，
該木頭分類為 Tj
g2 (X )
MAX
X
Tj
gct ((XX))
圖8.2.3 識別器示意圖
5
8.3 影像間的匹配對應
8.3.1 Harris 角點偵測法
視窗在影像上的移動，可得到強度變化情形：



(1)平面：往任何方向移動僅造成小變化
(2)含一條邊：與邊平行的變化量小；反之則大
(3)含角點或獨立點：往任何方向變化皆大
(1)
(2)
(3)
6
令
x f  f  (1,0,1)  x
y f  f  (1,0,1)t  y
利用高斯函數來平滑雜訊的影響，令
C   x y  G
B  2y  G
A  2x  G
綜合灰階梯度變化之影響可表示為
 A C  x 
Ax  By  2Cxy  ( x y) 
 
C
B

 y 
 ( x y)M ( x y)t  E
2
2
7

函數E是一種局部自我關聯的函數，矩陣M就是函數E的代表。
矩陣M的兩個特徵值代表：
 1 和2皆很小：代表視窗內為平滑區
 1 和2中，一大一小：代表含一邊的區域
 1 和2皆很大：代表含角點的區域
圖8.3.1.1 矩陣M的特徵值所代表的意義
8

影響值 R=det(M)-k×(trace(M))2決定是否小區域內有角點，
判斷的準則為
 R  0 且 R  0 ：代表平滑區
 R  0：代表含一邊的區域
0：代表有角點的區域
 R
利用以上方法可以將所有角點找出來。
(a)原始影像
(b)找出的角點集
9
8.3.2 SIFT關鍵點偵測法

將影像進行縮減取樣 n 次，例如n=2可分成3個影像尺度，進一
步利用 m 個不同標準差的高斯函數，例如m=4，高斯函數標準
差分別是  , k , k 2 , k 3，將不同高斯函數與3個影像尺度作迴積
運算：
L( x, y, k t )  G( x, y, k t )* I ( x, y)
接著對兩相鄰的尺度做DOG:
D( x, y, k t )  (G( x, y, k t 1 )  G( x, y, k t ))* I ( x, y)
 L( x, y, k t 1 )  L( x, y, k t )
10
圖8.3.2.1 DOG金字塔
11
考慮每個像素與周圍8個像素及上下層同位置周圍9個像素
做比較，如果該像素為極值，則設為候選關鍵點。
圖8.3.2.2 DOG金字塔找極值
12
以每個候選關鍵 X 點當原點，算出
DT (0)
1 T  2 D(0)
D( X )  D(0) 
X X
X
2
X
2
X
其中 X  ( x, y,  )T 的三維座標，令座標
DT (0)  2 D(0)

X
X 2
時有極值點D( )。
0
1
  D (0)  D(0)
 

2

X

 X
代表一個極值點，如果 內的三個座標值都小於0.5; 若大於0.5
則依照該方向來內插出新的候選關鍵整數點，繼續用新的候選關
鍵點計算二階泰勒展開式求值。
若 | D( ) | 0.03，代表此區域的對比較低，應移除此關鍵點。
2
T
13
除掉邊上的候選關鍵點，先求出赫斯矩陣
 Dxx Dxy   Dxx Dxy 
H 



D
D
D
D
yy   xy
yy 
 yx
H的特徵值 1 和 2 可計算如下:
  Dxx  
Dxy  
det( H   I )  det  
 
 Dxy
D


yy


  2  (Dxx  Dyy )  (Dxx Dyy  Dxy 2 )  0
可得
1  2  Dxx  Dyy  trace(H )
12  Dxx Dyy  Dyy 2  det(H )
14
令 1  2 且1  2 ，得到
(trace( H ))2 (1  2 )2 (2  2 )2 (  1)2



2
det( H )
12
2

得判斷式
(trace( H ))2 ( ' 1)2
<
det( H )
'
令  '  10，若判斷式成立，移除邊上的候選關鍵點。
15
以剩餘關鍵點為中心的區塊，計算區塊內每個像素的梯度值
m( x, y) 與梯度方向  ( x, y )
m( x, y) 
 L( x  1, y)  L( x 1, y)   L( x, y  1)  L( x, y 1)
2
2
 ( x, y)  tan1 (L( x, y 1)  L( x, y 1)) / (L( x 1, y)  L( x 1, y))
關鍵點的旋轉不變性(Rotation Invariant)。
16
圖8.3.2.3 決定關鍵點的特徵向量
17
找出DOG金字塔中的極值
移除在對比較低區域裡的候選關鍵點
移除位於邊上的候選關鍵點
配置關鍵點的主要方向
產生特徵向量
圖8.3.2.4 SIFT流程圖
18
(a) 輸入的腦血管影像
(b) SIFT演算法選取到的關鍵點
圖8.3.2.5
19
8.3.3 點集合匹配法

另一種作法是尋找轉移矩陣 T ，使 ToF 與 F 有最小誤差。
令 di (T )，1  i  m ，代表 F 中的 f i 與 F 中的 fi 之距離
m
p(d1 (T ), d 2 (T ),

, d m (T ) | T )   p(di (T ))
i 1
引入最大可能(Maximum-likelihood)的觀念，得
L(T )  ln p(d1 (T ), d 2 (T ),
m
, d m (T ) | T )   ln p(di (T ))
i 1
其中 p(d) 為一常態分佈(Normal Distribution)。

匹配的精神就是找一個 T 使得上式有最大值。
20
8.4 匹配演算法原理
8.4.1 動態規劃式的BSSC解法

BSSC(Banded String-to-String Correction)之目的為在樣本和正本
之間找出最大的匹配點與匹配點的位置。
給一個例子如下：
樣本= PP
1 2 P3 P4 P5    P8 P9 P10
正本= TT
1 2T3T4T5    T8T9T10
將 Pj 和 Ti 想像成特徵值。因為視差 
的關係，Pj 只能與 Ti  範圍內的特
徵匹配到。
圖8.4.1.1的搜尋範圍就是黑色粗邊
框住的範圍而得到的最佳匹配就是
鋸齒形的路徑上之黑圓點集。
P1 P2 P3 P4 P5 P6 P7 P8 P9 P10
T1
T2
T3
T4
T5
T6
T7
T8
T9




 



  
T10


|   3 |
圖8.4.1.1 匹配結果
21

匹配算子
a  b，R(a)=b的花費定為1
取代算子R  a = b，R(a)=b的花費定為0
刪除算子D  D(a)=  的花費定為1
插入算子I  I( )=a的花費定為1
圖8.4.1.1一共含有下列十四個運算
(1) I( )=P1
(6) I()=P4
(2)R(T1)=P2
(3)R(T2)=P3
(7)R(T5)=P5
(8)R(T6)=P6
(4)D(T3)=  (5)D(T4)= 
(9)D(T7)=  (10)D(T8)= 
(11)I( )=P7 (12) I( )=P8 (13)R(T9)=P9 (14)R(T10)=P10
22

動態規劃核心式子
Edit[i, j]  min(Edit[i 1, j]  edit(Ti , ),
Edit[i, j 1]  edit(, Pj ), Edit[i 1, j 1]  edit(Ti , Pj ))
表示將T[1…i]轉換成P[1…j]的花費。
T [1...i]  T1T2 ...Ti
P[1... j]  P1 P2 ...Pj
Edit[0, k]=k
Edit[k, 0]=k
edit(Ti, Pj)表示R(Ti)=Pj的花費
)表示D(Ti)=  的花費
edit(  , Pj)表示I()=Pj的花費
edit(Ti,

BSSC的時間複雜度為 (n)
n為正本長度，  為視差
23
8.4.2 KMP演算法

樣本先進行事前處理
利用樣本中子字串(Substring)與前置字串(Prefix String)的吻合度，
並記錄其吻合的長度於陣列 J [ ]中。
i 1 2 3 4 5 6 7 8 9 10
P[3]=a=P[1] J[3]=1
P [i] a c a c a a a c a c
P[3…5]=aca=P[1…3] J[5]=3
P[7…10]=acac=P[1…4] J[10]=4 J [i] 0 0 1 2 3 1 1 2 3 4

KMP字串匹配演算法
T [ ]=cccacacaaacaccaa
T[i]  P[i]，1  i  3。T[4…13] = P[1…10]。
T[5]  P[1]，T[6] = P[1]  T[6…8] = P[1…3]
T[9]  P[4]，J [4]=2 P[1…J [4]-1]=P[1]=T[9- J[4]+1…8]=T[8]
陣列J [ ]提供了一個跳躍的機制，讓匹配動作可一直往右前進。

KMP的時間複雜度為O(m + n)
24
8.5 三維影像重建
8.5.1 稠密式視差估測
(Dense Disparity Estimation)

問題定義
稠密式視差估測是在L上的所有
像素中和R上的所有像素中找到
一個對應。

P
二階段式的分割與克服方法
左影像 (L)
Y
OL
Z
右影像 (R)
Y
X
'
OR
X'
Z'
圖 8.5.1.1 稠密視差估測示意圖
25

第一階段的分割與克服方法
1. 利用Marr-Hildreth測邊法得到主要特徵像素集：
L上第 N/2 列中的主要特徵集可表示為：

S LN 2  S LN 2 1 S LN 2  2 
S LN 2 N N 2

，1  NN 2  N
R上第 N/2 列中的主要特徵集可表示為：

S RN 2  S RN 2 1 S RN 2  2 
N 2
N 2
2. SL 和 SR 的匹配工作



S

S RN 2 N N 2 ， 1  N N 2  N




C SLN 2  v  , S RN 2 ( w)  min Cmatch S LN 2  v  , S RN 2  w  , Crightocc S LN 2  v  , S RN 2  w  ,
Cleftocc
N 2
L
 v  , S RN 2  w  
其中 Cleftocc  S LN 2  v  , S RN 2  w   C  S LN 2  v  , S RN 2  w  1   COcc




Crightocc S LN 2  v  , S RN 2  w  C S LN 2  v  1 , S RN 2  w   COcc




Cmatch S LN 2  v  , S RN 2  w  C S LN 2  v  1 , S RN 2  w  1  CS N 2  v ,S N 2  w
L
R
26
S LN / 2
S LN / 2 1 S LN / 2 2
~
~
S RN / 2 1
S RN / 2 2 
CS N 2  v , S N 2  w   2 * l1  l2
L
~
~
取代算子：
S LN / 2 m N / 2 
S LN / 2 w
R
~
~
其中
S RN / 2
S RN / 2 v 
l1  SLN 2  v   SLN 2  v 1
~
~
l2  SRN 2  w  SRN 2  w 1
S RN / 2 n N / 2 
N 2
圖 8.5.1.2 SLN 2 和 SR
的匹配

11 l
1 l
mean    p1 1 I L S LN 2  v  1  p   q21 I R S RN 2  w  1  q 
2  l1
l2





2
2
1  1 l1
1 l2
N 2
N 2
I
S
v

1

p

mean

I
S
w

1

q

mean

 L L  
 R R  
2  l1 p 1
l2 q 1

2  
 


 


27
f
A
X
D
C
X
2*thr
B
(a) 右遮蔽花費
(b) 左遮蔽花費
thr
圖 8.5.1.3 右遮蔽花費和左遮蔽花費示意圖
COCC  k  f

2
1

  22 / 2; thr

0
thr

2
1

  22 / 2
3. 加入消去法則
若兩兩匹配位置的差，形成之序列為  5, 2,5,6,5  ，則 -2
的視差偏移不太正常，可以予以去除。假設在 N 2  k 列時
的平均視差偏移序列為  5,5,6, 4  。若新進來的視差偏移序
列為  5,7,3,16, 4  ，則 16 的視差偏移可予以去除。
28
三維物體

第二階段的分割與克服方法
A
B
C
1. 在二段匹配的子區間中，
再進一步找出個別的像素配對
2. 利用平均視差偏移量d，
可將第二階段的工作轉
1
成BSSC問題的解決上
SLN/2
1
SLN/2
1
SLN/2
1 N
圖8.5.1.4 二子區間匹配
R算子： C '  S LN 2  v   vˆ, S RN 2  w  wˆ  
W 2

x , y W 2

N/ 2
2
1(1) SSRNLN/2
1
1(3)NN
(2)
SSLNRN/2
11 SSRNL 2


I L N 2  x, S LN 2  v   vˆ  y  I R N 2  x, S LN 2  w   wˆ  y

W為事前定義好的小視窗
D算子：類似於前面定義的左遮蔽花費
I 算子：類似於前面定義的右遮蔽花費
29
圖8.5.1.5 對應的BBSC搜尋空間
30
圖8.5.1.6 輸入的二張影像
(a) 影像L
(b) 影像R
圖8.5.1.7
得到的視差圖
圖8.5.1.8
重建後的三維
五角大廈圖
31
8.5.2 相機校正

我們有三個座標系統，世界座標系統(World Coordinate
System) 、相機座標系統(Camera Coordinates System)和二維的
影像系統(Image System) 。
 x
 xw 
 y  R  y   T
 
 w
 z 
 zw 
 r1
R   r4
 r7
r2
r5
r8
r3 
r6 
r9 
Tx 
T  Ty 
Tz 
根據三角比例關係，可得
Xt  f
x
z
Yt  f
y
z
其中 f 為焦距長。
32
影像平面
鏡心
圖 8.5.2.1 相機座標系統和影像座標系統
理論上 ( X t , Yt ) 為( x, y ) 投射到影像座標系統的理想位置。
因為透鏡輻射效應，( x, y ) 實際投射到的影像座標系統的位置應
為 ( X a , Ya )。 滿足下式
X a  Ex  X t
Ya  Ey  Yt
33
上兩式中的誤差項 Ex 和 Ey 可表示成
Ey  Ya (1r 2  2r 4 )
Ex  X a (1r 2  2r 4 )
r
X a2  Ya2
在光學透鏡的輻射失真(Radial Distortion)的影響下，內部參數滿足
2
1，我們令  2  0，只關心 1 的求解。
前面所提的 ( X a , Ya ) 為實數，但在數位影像的座標系統一般皆為
整數座標 ( X I , YI )。假設感應器(Sensor)的 X (Y ) 方向共有 S x (S y ) 個感應
器。N x 是影像在x方向的解析度。令
d x'  d x
Sx
Nx
34
而 (Cx , Cy ) 為影像的中心座標。( X a , Ya ) 和 ( X I , YI ) 的關係可表示如下
X I  Cx  s
Xa
d 'x
YI  C y  s
其中 s 為待解的放大係數。
Ya
dy

已知 ( x, y ) 投影到影像座標系統的理想位置為 ( X t , Yt )   f , f
進而得到
同理可得
且

( X  Cx )d 'x
x
f  X a  Ex  I
 Ex
z
s
( X I  C x )d 'x ( X I  C x )d 'x 1r 2


s
s
X d 'x X d 'x 1r 2


s
s
d y Y  d y Y  1r 2  f
2
x
z
y
。
z
y
z
2
 Xd 'x 
 Xd 'x 
2
2
r 
  (d y (YI  Cy ))  
  (d y Y )
 s 
 s 
35
再利用三維世界座標系統和二維影像座標系統的關係可得
Xd 'x Xd 'x 1r 2
r x  r x  r x T

f 1 w 2 w 3 w x
s
s
r7 xw  r8 xw  r9 xw  Tz
d y Y  d y Y1r  f
2
r4 xw  r5 xw  r6 xw  Ty
r7 xw  r8 xw  r9 xw  Tz
旋轉矩陣 R 分別對 x 軸、 y 軸和 z 軸達到任一角度的旋轉，可寫成
下式
cos cos 

R    sin cos   cos sin  cos 
 sin sin   cos sin  cos 
sin cos 
cos cos 
 cos sin   sin sin  cos 
 sin  
cos  sin  
cos  cos  
至此，已經六個外部參數(Extrinsic Parameters)分別為 、、、Tx、Ty
和 Ty 及五個內部參數(Intrinsic Parameters)分別為 f、 s、1、Cx 和 C y。
接下來要利用數值的方法求解這些參數。
36
先解出五個外部參數，在二維影像座標上取一經過原點的向量
( X a , Ya ) ，且在相機座標上取一平行的向量 ( x, y ) ，所以得
( X a , Ya )  ( x, y)  0
X a y  Ya x  0
X a (r4 xw  r5 xw  r6 xw  Ty )  Ya (r1 xw  r2 xw  r3 xw  Tx )
令 Z w  0 則我們得到
Ya xw
Ya yw Ya
 X a xw
 Ty1r1 
 1 
 Ty r2 
 X a yw  Ty1Tx   X a
 1 
 Ty r4 
 T 1r 
 y 5
再利用五組世界座標 ( xw , yw , zw ) 和五組底片上的真實座標 ( X a , Ya ) 就
可解出 Ty1r1、Ty1r2、Ty1Tx、Ty1r5 、Ty1r4 。
37
令
 r1
C
 r4
則旋轉矩陣可改寫成
r2  Ty1r1 Ty1r2 
   1
1 
T
r
T
r5   y 4
y r5 

 rT
1 y

R   r4Ty
 r
 7
r2Ty
r5Ty
r8
r3 

r6 
r9 

利用 R 中每一行向量與列向量為單位長，可得下式

rT
1 y


R
r4Ty


12
2
2
2


  1  Ty (r1  r4 )

 
2
2
2

 1  T (r1  r2 )  



12
2
2
 1  Ty2 (r4  r5 )  

 

(1  MTy2 )1 2 

2
y
r2Ty
r5Ty
12
2
2
2

 1  Ty (r2  r5 ) 


2
2
12
2
其中 M  r1  r2  r4  r5 。
38
因為 R 為正交矩陣，任兩行的內積必為零。將 R 的第一行和第二行
做內積，可以得到
(r1 r5  r2 r4 )2 Ty4  MTy2  1  0
此等式可解得 Ty2 的兩個解，如下
12
M   M  4(r1 r5  r2 r4 ) 
T 
2(r1 r5  r2 r4 )2
2
2
2
y
前面曾推導過
d y Y  d y Y1r  f
2
r4 xw  r5 yw  r6 zw  Ty
r7 xw  r8 yw  r9 zw  Tz
暫時將 1 和 Z w 設為零，得
dyY  f
r4 xw  r5 yw  Ty
r7 xw  r8 yw  Tz
39
令 t1  r4 xw  r5 yw  Ty 且 t2  r7 xw  r8 yw 上式可改寫成
t1

f

 d y Y     t2 d y Y
Tz 
再利用兩組以上的世界座標 ( xw , yw , zw ) 和真實座標 ( X a , Ya )
解出 f 和 Tz 。然後以此 f 和 Tz 及令 1  0 為初始值來解線性方程式
t1
d y Y  d y Y1r  f
t2  Tz
2
如此即可透過數值的解法求得 ( f , Tz , 1 ) 的逼近解。
40
8.6 二維影像的深度計算
A
物體上一點
Epipolar面
Epipolar線

左影像
左邊相機
焦點中心
E

B
Stereo底線

右影像
F
C
圖8.6.1 Epipolar面
右邊相機
焦點中心
41
範例1:已知兩部相機的焦距[參見式(1.3.1)]相同且均為 f，假設左邊
的Epipolar線之長度為 xl 而右邊的Epipolar線之長度為 xr。已知物件
上的一點 A( x, y, z )，試問 z 之值為何?
解答:
根據給定條件，圖8.6.2為對應的示意圖，圖中的B和C代表左相機和
右相機。利用圖8.6.2中， BFG  BDA 和 CHI  CEA ，透過兩
相似三角形的邊比例關係，可得到
D
E


 A( x, y, z)
x xl

z f
x  b xr

z
f
xr
F
xl

f

H
G
B
b


I
C
圖8.6.2 利用兩平行相機求z值
42
在上式中，我們假設兩不相機相距b，座標系統的原點設在B點。可
進一步推得
zxl
x
f
x  ( zxr  bf ) f
再利用 zxl  zxr  bf ，可得到深度為
bf
z
xl  xr
解答完畢
43