Kalibrace a rekonstrukce

Kalibrace a rekonstrukce 2. úloha z 33PVI

Vypracoval: Jan Doležel, [email protected]

a) Autokalibrace z rotace b) Rekonstrukce a kalibrace

Autokalibrace z rotace Zadání Nalezněte kalibrační matici pro scénu zadanou obrázky. Máme k dispozici šest obrázků z Karlova náměstí a pět korespondencí bodů mezi prvním snímkem a ostatními.

Teorie Vyjdeme ze základní rovnice kamery, kterou můžeme rozepsat:

X 1

X  1

X  1

αu = Ρ  = ( A | − AT )  = (KR | − KRT )  Kde K ∈ ℜ 3x 3 je kalibrační matice kamery, R ∈ ℜ 3x 3 je matice rotace a T ∈ ℜ 3 je střed kamery v souřadnicích světa. α ∈ ℜ je koeficient, vyjadřující, že body ležící na přímce procházející středem kamery, se zobrazí do jednoho bodu v obrázku. u ∈ ℜ 3 je bod obrazu odpovídající bodu ve světě X ∈ ℜ 3 . Kalibrační matice je pro všechny snímky stejná, mění se pouze rotace, která zachovává ortonormalitu. Díky tomu, nám stačí znát matici K, abychom mohli měřit úhel, který svírají dva body. K chceme ve tvaru horní trojúhelníkové matice:  k11  K = 0 0 

k12 k 22 0

k13   k 23  k 33 

Bez újmy na obecnost můžeme předpokládat, že T = 0. Potom můžeme pro dvě kamery psát:

α i u i = Ai X = KR i X α j u j = A j X = KR j X Již víme, že β j u j = H ij u i vyjadřuje převod mezi souřadnicemi bodů dvou kamer, −1

Hij je homografie. Matici Hij můžeme rozepsat H ij = H j H i , kde τ i H i = KR i = Ai , −1

stejně tak R ij = R j R i . Takže můžeme upravovat:

KR j

−1

H ij = H j H i = H ij

ij

Zavedeme H =

H ij

1

τ

= 3

j

−1

τ i R i K −1 =

τi K ( R j R i ) K −1 = δ ij KR ij K −1 j τ −1

δ ij KR ij K −1

(δ ) (K R ij

3

1

3

ij

K −1

)

1

= KR ij K −1 3

Protože R je ortonormální rotace: ij

ij T

ij

ij

I = R ij R ij = ( K −1 H K )T ( K −1 H K ) = K T H K −T K −1 H K T

ij T

ij

H K −T K −1 H = K −T K −1

Označíme si ω = K −T K −1 , ω ∈ ℜ 3 x 3 a ω T = (K −T K −1 ) = K −T K −1 = ω . Pro výpočet ω tedy potřebuji šest rovnic. T

 h11   h12 h  13

h21 h22 h23

h31  ω11 ω12 ω13  h11   h32  ω12 ω22 ω23  h21 h33  ω13 ω23 ω33  h31

=> =>

h12 h22 h32

h13   ω11 ω12 ω13     h23  =  ω12 ω22 ω23  = ω h33   ω13 ω23 ω33 

( h11h21h31 )ω ( h11h21h31 )T = ω11

=>

=>

ω kl = ( h1k h2 k h3k )ω ( h1l h2 l h3l )T , kde (k,l) = {(1,1), (1,2), (1,3), (2,2), (2,3), (3,3)} =>

ω kl = ω11 h1k h1l + ω12 ( h1k h2 l + h2 k h1l ) + ω13 ( h1k h3l + h3k h1l ) + ω22 h2 k h2 l + ω 23 ( h2 k h3l + h3k h2 l ) + ω33 h3k h3l Nyní můžeme sestavit matici A, aby Aω = 0 :

...    ij ij  ij ij ij ij ij ij ij ij ij ij 2 ⋅ h11 h21 2 ⋅ h11 h31 h21 h21 2 ⋅ h21 h31 h31 h31  ω11   h11 h11 − 1   ij ij ij ij ij ij ij ij ij ij ij ij ij ij  h ij h ij h11ij h22 + h21 h12 − 1 h11ij h32 + h31 h12 h21 h22 h21 h32 + h31 h22 h31 h32  ω12  11 12  ij ij  ij ij ij ij ij ij ij ij ij ij ij ij ij ij h11ij h23 + h21 h13 h11ij h33 + h31 h13 − 1 h21 h23 h21 h33 + h31 h23 h31 h33  ω13   h11 h13 ij ij ij ij ij ij ij ij ij ij  h ij h ij  ω  = 0 2 ⋅ h h 2 ⋅ h h h h − 1 2 ⋅ h h h h 12 22 12 32 22 22 22 32 32 32  12ij 12ij  22  ij ij ij ij ij ij ij ij ij ij ij ij ij ij ij ij h12 h23 + h22 h13 h12 h33 + h32 h13 h22 h23 h22 h33 + h32 h23 − 1 h32 h33  ω23   h12 h13   ij ij  ij ij ij ij ij ij ij ij ij ij  ω h h 2 ⋅ h h 2 ⋅ h h h h 2 ⋅ h h h h − 1 13 23 13 33 23 23 23 33 33 33  13 13  33    ...   K vypočítáme z ω pomocí Choleskyho faktorizace.

Vypracování Pro vypracování úlohy jsme použili systém Matlab. Navíc jsme pro lepší výsledky použili normalizaci obrázků do čtverce s hranou velikosti 1. Pokusil jsem se

vypočítat matici K i bez normalizace, ale v tomto případě nešla Choleskyho faktorizace provést, protože ω nebyla pozitivně definitní. S normalizací vyšlo K:  2844,1 − 149,31 427,05    K = 0 4964,4 517,18   0 0 1  

Pro matici K často platí (pro čtvercové pixely): a 0 b    K = 0 a c  0 0 d   

kde a odpovídá převodu jednotek světa do pixelů, d jak daleko je střed kamery od snímaného čipu a je-li d = 1, b a c by měli být hodnoty středu obrázku (v pixelech). V tomto tvaru můžeme ω spočítat z

ω kl = ω11 ( h1k h1l + h2 k h2l ) + ω13 ( h1k h3l + h3k h1l ) + ω 23 ( h2 k h3l + h3k h2l ) + ω33 h3k h3l a (k,l) = {(1,1), (1,3), (2,3), (3,3)}. Takto vypočítané K vychází: 0 612,41  3559   K = 0 3559 331,04   0 0 1  

Vzhledem k tomu, že korespondence nebyly úplně přesné a zavedli jsme pro K velké omezení, b a c vyšli hodně mimo střed obrázku (426, 568).

Rekonstrukce a kalibrace Zadání Proveďte rekonstrukci scény z poskytnutých obrazů, nalezněte transformační matice kamer.

Obr.: Dva z osmi poskytnutých obrazů scény

Teorie Na začátku známe pouze souřadnice některých bodů v obrázcích - uij , které jsme si museli sami zjistit (16 okótovaných bodů a vrchol jehlanu), a chceme znát souřadnice odpovídajících bodů Xi. Situaci mezi dvěma kamerami popisuje obrázek:

Vyjdeme z těchto rovnic: (Y − C ) − b = Y − C ' x − b = x'

∃α > 0...α u β = x β

x β = A(Y − C )

∃α ' > 0...α ' u' β ' = x' β '

x' β ' = A' (Y − C ')

A −1 x β − A −1 b β = A' −1 x' β ' x β − b β = AA' −1 x' β '

α u β − b β = α ' AA' −1 u' β ' kde Y jsou souřadnice X v soustavě světa; β a β’ jsou souřadné soustavy kamer se středem v C, C’; x a x’ jsou souřadnice X v soustavách β a β’; b a b’ jsou souřadnice kamery se středem v C’ v soustavě β (vzhledem k C); A, A’ jsou matice přechodu od soustav β a β’ k souřadné soustavě světa. Následující úpravy jsou již neekvivalentní v tom smyslu, že vedou jen vpřed:

α u β − b β = α ' AA' −1 u' β ' =>

=>

α b β × u β − b β × b β = α ' b β × ( AA' −1 u' β ' )

α u β (b β × u β ) = α ' u β (b β × ( AA' u' β ' )) −1

=>

=>

−1

0 = u β (b β × ( AA' u' β ' ))

Jelikož vektorový součin lze vyjádřit jako matici lineárního zobrazení (pro T kterou navíc platí: [ y ]× = −[ y ]× ), můžeme poslední rovnici psát ve tvaru: 0 = u βT [bβ ]× AA' −1 u β, '

Můžeme zavést matici F = [bβ ]× AA' −1 , kterou nazýváme fundamentální maticí. Body e, e’ z obrázku nazýváme epipóly (obraz druhé kamery v prvním obrázku, respektive naopak), přímky l, l’ epipoláry (spojují obrazy bodů s epipóly) a rovinu π epipolární rovinou. Hodnost matice F je 2, protože matice A, A’ mají hodnost 3 a hodnost matice vektorového součinu je 2 (právě tehdy, když b ≠ 0). Pro epipóly platí: eβT F = λbβT F = λbβT [bβ ]× AA' −1 = 0

Feβ, ' = λ ' Fbβ ' = λ ' FA' A −1bβ = λ ' [bβ ]× AA' −1 A' A −1bβ = λ ' [bβ ]× bβ = 0

Epipóly představují levý a pravý nulový prostor matice F. Pro epipoláry platí následující rovnice: 0 = l ' T u' β '

=>

l ' = u βT F

0 = u βT F

=>

l = Fu' β '

Ze souřadnic bodů u β = (u, v,1)T a u' β = (u' , v' ,1)T a 0 = u βT Fu' β ' můžeme spočítat matici F:

 u1u'1   u2 u ' 2    un u ' n

u1v '1 u2 v ' 2 ... un v ' n

u1 u2

v1u'1 v2 u ' 2

v1v '1 v2 v ' 2

v1 v2

u'1 u'2

v '1 v'2

un

vn u ' n

vn v ' n

vn

u' n

v'n

    1   1    1     

f11   f12  f13   f 21  f 22  = 0  f 23   f 31  f 32   f 33 

Takto vypočítaná matice F má však hodnost 3. My chceme, aby měla hodnost dva. Docílíme toho tak, že provedeme její SVD rozklad, vynulujeme prvek na pozici σ33 a znovu mezi sebou vynásobíme. Nyní se dostáváme k rekonstrukci scény. Mějme dvě kamery s projekčními maticemi P a P’:

α i xi = PX i α i ' xi ' = P ' X i

=>

 PX i   αi

T

  P' X i  F    αi '

  = 0 

=>

X iT P T FP ' X i = 0

Když označíme Q = P T FP' , Q ∈ ℜ 4 x 4 . Snadno se přesvědčíme, že Q T = −Q (a hlavně P T FP ' = − P' F T P ) a Qii = 0 . Nyní můžu mezi P a Xi vložit H −1 H tak, aby PH −1 = ( I | 0 ) : ~~

~

αx = PX = PH −1 HX = P X = ( I | 0) X ~ ~

~

α ' x' = P' X = P' H −1 HX = P ' X = ( A | a ) X

Zvolili jsme tedy soustavu světa shodnou se soustavou první kamery. Máme tedy: ( I | 0 ) T F ( A | a ) = −( A | a ) F T ( I | 0 )

I   A   F ( A | a ) = −  F ( I | a ) a 0  AT F T 0   FA Fa     = − T T 0 0 a F 0     Vidíme, že a T F F = 0 a Fa = 0 , a protože i Fe' = 0 , je a = αe' . Můžeme tedy psát ~ ~ ( A | a ) X = ( A | αe' ) X . Dále vidíme, že FA = − AT F T . Pokud zvolíme za A = SF T , S ∈ ℜ 3 x 3 zjistíme, že: AT F T = − FA => FS T F T = − FSF T => S T = −S . Ale takové S je například maticový součin a pro nás nejvhodnější je [e']× . Dostáváme tedy:

~

αx = ( I | 0) X

α ' x' = ([e']× F T | αe' ) X

~

Vidíme, že druhá kamera má střed v nevlastním bodě ve směru epipólu. Z těchto dvou kamer můžeme získat homogenní souřadnice bodů ve světě:  [xi ]× Pi    X i = 0  [x' i ]× P' i 

Z naměřených bodů a homogenních souřadnic bodů ve světě můžeme spočítat všechny projekční matice ( xij = ( x1ji , x2ji , x3ji )T ):

α i j xi j = P j X i  ...  j T  x2 i X i xj XT  3i i  0000  ... 

−x X 0000 x3ji X iT j 1i

T i

  0000  p1jT    − x1ji X iT  p2jT  = 0  − x2ji X iT  p3jT   

Využijeme projekční rovnici ještě jednou a vypočítáme alfy. Nyní můžeme znovu přepočítat všechny projekční matice a homogenní souřadnice bodů ve světě – budeme tak mít jejich lepší aproximaci. Na část projekční rovnice vlevo od rovnítka provedeme SVD rozklad a P a X vypočítáme třeba jako P = U D a X = DV T . V matici P ∈ ℜ 4 x 3 M , kde M je počet kamer, máme pod sebou projekční matice všech kamer a v matici X ∈ ℜ Nx 4 , kde N je počet bodů, jsou ve sloupcích uloženy homogenní souřadnice bodů ve světě. Nyní již jen z rovnice αYi = HX i a několika vybraných bodů zjistíme H, která nám nalezené body transformuje do nějaké „pěkné“ souřadné soustavy. Taková je například ortonormální s jednotkou délky o velikosti hrany jehlanu, s počátkem v jednom z bodů podstavy. V této soustavě známe souřadnice pěti bodů: 4 bodů podstavy a vrcholu jehlanu. K výpočtu H potřebujeme taky pět bodů, ale žádné čtyři nesmí ležet v rovině. Přidáme tedy ještě jeden, u kterého známe souřadnice – nevlastní bod ve směru osy z (v obrázku nahoru). Jeho souřadnice v obrázku nalezneme jako průsečík dvou rovnoběžek (rovnoběžných s osou z). Matici H pak vypočteme:

 ...  T  Y2i X i Y X T  3i iT  Y4i X i  0000   0000   0000  ... 

− Y1i X

T i

0000 0000 Y3i X iT Y4i X iT 0000

0000 − Y1i X iT 0000 − Y2i X iT 0000 Y4i X iT

  0000  0000  h1T    − Y1i X iT  h2T  =0 0000  h3T    − Y2i X iT  h4T   − Y3i X iT   

Pomocí této matice můžeme všechny body zobrazit v námi zvolené soustavě.

Vypracování K vypracování byl i tentokrát použit program Matlab. Nejdříve jsme určili souřadnice 17ti zvolených bodů ve všech obrázcích pomocí funkce ginput, znormalizovali do čtverce o hraně délky 1 a pak postupovali podle výše uvedeného. Vyskytnul se problém při vytváření matice H, protože pouze z pěti bodů nešlo matici vytvořit. Musel jsem přidat ještě jeden bod podstavy. Ačkoli lineárně závislý na ostatních, matici již šlo vypočítat správně. U výsledné rekonstrukce velice závisí na přesnosti naměření bodů, ze kterých scénu vytváříme, proto se výsledky velmi liší dle použitých obrazů v jednotlivých fázích a bodů, ze kterých vytváříme přímky a počítáme nevlastní bod. Uvádím dvě scény, které vypadají nejlépe:

Scéna z kombinace 1-7

Scéna z kombinace 1-8

Při použití přiloženého skriptu, můžete vyzkoušet i jiné kombinace. Kombinace využívající obrazy 4, 5 nebo 6 se ukazují jako nepoužitelné. Ostatní kombinace dávají různě přijatelné výsledky.

Závěr Nalezl jsem kalibrační matici pro obrázky z Karlova náměstí a rekonstruoval scénu s jehlanem. U obou obrázků jsem zaznamenal velký přínos normalizace vstupních bodů na čtverec o hraně s délkou jedna. Bez použití normalizace nebylo možné získat kalibrační matici vůbec a při rekonstrukci byla výsledná scéna více deformovaná (některé, původně, kolmice se staly téměř rovnoběžkami).

Ukázka reprojekce bodů ve 3. obraze (modrá první vypočítané Pj, červená druhé přepočítání Pj). Vpravo dole obraz nevlastního bodu – průsečík přímek, tvořících strany horního čtverce