Modello della Telecamera, Calibrazione e Rettificazione

A.A. 2010/2011
Elaborazione dell'immagine M
Modello della Telecamera, Calibrazione e
Rettificazione
Samuele Salti, Luigi Di Stefano
Proiezione prospettica
•
Consideriamo
M=[x,y,z]T
– un punto dello spazio 3D, M=[x,y,z]T, le
cui coordinate sono espresse nel sistema
di riferimento (sdr) solidale con la
telecamera (sdr “standard”).
C
– La sua proiezione sul piano immagine I,
indicata con m=[u,v]T
•
Le equazioni, non lineari, che
esprimono le coordinate immagine in
funzione delle coordinate 3D sono date
da:
⎧
⎪⎪u =
⎨
⎪v =
⎪⎩
x
m=[u,v]T
u
c
f
z
v
I
y
f
x
z
f
y
z
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Spazio Proiettivo
•
Lo spazio fisico è uno spazio euclideo 3D (R3) in cui i punti, una volta fissato un
sistema di riferimento, si possono rappresentare con vettori tridimensionali.
– In questo spazio si dice che le rette parallele non si intersecano, o si intersecano “all’infinito”
– La rappresentazione di questi punti all’infinito non è possibile in questo spazio vettoriale.
•
•
Immaginiamo di aggiungere una coordinata alle nostre terne, per es.
( x y z ) diventa ( x y z 1)
e consideriamo entrambe le ennuple rappresentazioni valide dello stesso punto.
In aggiunta, non vincoliamo la quarta coordinata ad essere 1, ma definiamo
(x
•
•
•
y
z 1) ≡ ( 2 x 2 y 2 z 2 ) ≡ ( kx ky kz k ) ∀k ≠ 0
In questa rappresentazione, un punto dello spazio è quindi identificato da una classe
di equivalenza di quadruple, dove quadruple equivalenti differiscono solo per un
fattore moltiplicativo.
Questa è la cosiddetta rappresentazione in coordinate omogenee del punto (x, y, z).
Lo spazio ottenuto rappresentando i punti in coordinate omogenee è detto spazio
proiettivo, P3.
È immediata l’estensione a spazi euclidei di qualunque dimensione (Rn Pn)
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Punti all’infinito
• Gli unici punti che non hanno un corrispondente nello spazio euclideo, sono i
punti per cui la nuova coordinata vale 0, es (x, y, z, 0).
• L’ipotetico corrispondente (inesistente nello spazio euclideo) sarebbe
rappresentato dalle coordinate (x/0, y/0, z/0), ovvero “coordinate infinite”.
• Usando le coordinate omogenee è quindi possibile rappresentare
analiticamente i punti all’infinito come
tutti i punti la cui ultima coordinata vale 0.
• Il punto di coordinate (0, 0, 0, 0) non è definito.
– Questo punto NON è l’origine del sistema di riferimento scelto per lo spazio
euclideo (0, 0, 0), il quale ha coordinate omogenee (0, 0, 0, k), k≠0.
3
• Si può dimostrare che tutti i punti all’infinito di P giacciono su un piano, che
prende il nome di piano all’infinito.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Un esempio in R2
• Considerando il sistema che definisce l’intersezione di due rette
⎧ ax + by + c = 0
⎨
⎩a′x + b′y + c′ = 0
questo ha soluzioni (Teorema di Cramer)
−c b
− c ′ b′
x=−
a b
a ′ b′
a −c
a ′ − c′
y=−
a b
a′ b′
• Se le rette sono parallele, il denominatore si annulla, e quindi non c’è
soluzione (o impropriamente, la soluzione è un punto “all’infinito”).
• Se le rette sono distinte, almeno uno dei tre valori presenti nella formula è
comunque non nullo.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Un esempio in R2
⎛
⎜
⎜x =
⎜
⎜
⎝
−c b
−c′ b′
,
a b
a′ b′
a −c ⎞
⎟
a ′ − c′ ⎟
y=
a b ⎟
⎟
a′ b′ ⎠
≡
⎛ −c b
,
⎜
′
′
−
c
b
⎝
a −c
,
a ′ − c′
P2
a b⎞
⎟
a′ b′ ⎠
R2
• Mantenendo le tre coordinate separate, e usando la rappresentazione in
coordinate omogenee, si passa nello spazio proiettivo P2.
• In questo spazio proiettivo il sistema precedente ha sempre una soluzione,
quella indicata dalla formula in alto.
• Si può facilmente verificare che il punto all’infinito di un fascio di rette si
rappresenta in coordinate omogenee aggiungendo uno zero ad un qualsiasi
vettore euclideo diretto come il fascio:
⎡bc '− b ' c ⎤ ⎡ −b ⎤ ⎡b ⎤ ⎡ −b '⎤ ⎡b ' ⎤
⎢ a ' c − ac '⎥ ≈ ⎢ a ⎥ ≈ ⎢ − a ⎥ ≈ ⎢ a ' ⎥ ≈ ⎢ − a '⎥
⎢
⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎢⎣ 0
⎥⎦ ⎢⎣ 0 ⎥⎦ ⎢⎣ 0 ⎥⎦ ⎢⎣ 0 ⎥⎦ ⎢⎣ 0 ⎥⎦
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Riassumendo
n
• Ogni spazio euclideo R può essere esteso in una nuova struttura geometrica,
n
il cosiddetto spazio proiettivo P , rappresentandone i punti in coordinate
omogenee.
• La rappresentazione in coordinate omogenee prevede l’aggiunta di una nuova
coordinata k allo spazio vettoriale
– k ≠ 0 rappresenta punti effettivi di Rn le cui coordinate sono xi/k, i = 1…n
– k = 0 rappresenta punti ideali di Rn , i cosiddetti punti all’infinito.
• Questa nuova struttura permette quindi di rappresentare e trattare
analiticamente in maniera omogenea (senza introdurre eccezioni o casi speciali)
anche i punti all’infinito.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Proiezione prospettica in coord. omogenee (1)
• Abbiamo visto che sussiste una relazione non lineare tra le coordinate della
scena e quelle del punto immagine:
u=
f
f
x v= y
z
z
• Torniamo al nostro punto M e alla sua proiezione prospettica m, le cui
coordinate euclidee nel S.d.R. Standard e in quello immagine avevamo
definito rispettivamente come M = [x,y,z ]T e m = [u,v ]T .
• Espresse in coordinate omogenee (che indicheremo sempre con una tilde
sovrastante il vettore) le loro rappresentazioni diventano:
⎡u ⎤
= ⎢v ⎥
m
⎢ ⎥
⎢⎣ 1 ⎥⎦
Elaborazione dell'immagine M
⎡ x⎤
⎢ y⎥
=⎢ ⎥
M
⎢z⎥
⎢ ⎥
⎣1 ⎦
Samuele Salti, Luigi Di Stefano
Proiezione prospettica in coord. omogene (2)
• In coordinate omogenee (e quindi considerando la trasformazione come tra
spazi proiettivi), la proiezione prospettica diviene una trasformazione lineare.
⎡ x⎤
⎢f z⎥
⎡u ⎤ ⎢
⎥ ⎡ fx ⎤ ⎡ f
⎢ v ⎥ = ⎢ f y ⎥ = ⎢ fy ⎥ = ⎢ 0
⎢ ⎥ ⎢ z⎥ ⎢ ⎥ ⎢
⎢⎣ 1 ⎥⎦ ⎢
⎥ ⎢⎣ z ⎥⎦ ⎢⎣ 0
1
⎢
⎥
⎢⎣
⎥⎦
• Usando la notazione matriciale
= PM
km
0
f
0
⎡ x⎤
0 0⎤ ⎢ ⎥
y⎥
⎥
⎢
0 0⎥
⎢z⎥
1 0 ⎥⎦ ⎢ ⎥
⎣1 ⎦
Esempi
-VP fascio generico
- VP fascio // asse z
- VP fascio // piano immagine
o anche, dove ≈ indica “uguale a meno di un fattore di scala arbitrario”
≈ PM
m
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
La matrice di proiezione prospettica
rappresenta il modello geometrico della telecamera, e viene
• La matrice P
detta matrice di proiezione prospettica (PPM).
• Nel caso in cui le distanze sono misurate in unità di distanze focali ( f = 1) la
PPM diviene
⎡1 0 0 0⎤
P = ⎢0 1 0 0⎥ = [I | 0]
⎢
⎥
⎢⎣0 0 1 0⎥⎦
• Questa forma “estrema” della PPM rappresenta l’essenza della proiezione
prospettica. Ad essa si fa sovente riferimento come matrice di proiezione
prospettica “standard” o “canonica”.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Un modello più realistico
• Per rendere più realistico il nostro modello in un sistema di acquisizione reale
è necessario tenere conto:
– della digitalizzazione dell’immagine;
– della trasformazione rigida tra la telecamera e la scena (roto-traslazione).
Z
X
x
Y
O
u
C
v
R, T
f
c=[uo, vo]T
z
I
y
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Digitalizzazione
Δu
Δv
c (u0,,v0)
Δu=dimensione orizzontale del pixel
Δv=dimensione verticale del pixel
f
x
z
f
v= y
z
u=
1
Δu
1
→v=
Δv
→u =
f
x = ku
z
f
y = kv
z
f
x + u0
z
f
y + v0
z
La digitalizzazione viene considerata inserendo nelle formule di proiezione lo
scaling lungo i due assi dovuto alla quantizzazione del piano immagine
nonché la traslazione del piercing point (punto principale) dovuta alla scelta
del sistema di riferimento dell’immagine (angolo in alto a sinistra).
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
La matrice degli intrinseci
• Tenendo conto delle nuove relazioni, la PPM può essere riscritta come:
⎡ fku
⎢
P = ⎢ 0
⎢ 0
⎣
0
fkv
u0
v0
0
1
0 ⎤ ⎡ fku
⎥ ⎢
0⎥ = ⎢ 0
0 ⎥⎦ ⎢⎣ 0
0
fkv
0
u0 ⎤ ⎡ 1 0
⎥⎢
v0 ⎥ ⎢0 1
1 ⎥⎦ ⎢⎣0 0
0 0⎤
⎥
0 0 ⎥ = A [ I | 0]
1 0 ⎥⎦
• La matrice A, che modella le caratteristiche del sensore, è detta matrice dei
parametri intrinseci.
• I parametri intrinseci possono essere ridotti al minimo necessario ponendo
αu=fku , αv=fkv : si tratta della lunghezza focale espressa in pixel orizzontali e
verticali. I parametri intrinseci minimi sono quindi 4.
• Il modello più generale prevede un quinto parametro nella matrice degli
intrinseci, detto skew, che rappresenta l’angolo tra gli assi del sistema di
riferimento del sensore. La sua cotangente occuperebbe la posizione A[1,2],
ma in pratica è sempre 0 ( = ctg(π/2) ).
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Trasformazione rigida fra i S.d.R. (1)
• Quanto visto finora si basava sull’assunzione di un S.d.R. per la scena molto
particolare, coincidente con quello standard della telecamera.
• Più in generale il S.d.R. della scena sarà legato a quello standard da
– una rotazione intorno al centro ottico (modellabile da una matrice ortogonale R)
– uno spostamento (modellabile da un vettore di traslazione t)
• La relazione tra le coordinate di uno stesso punto nei due S.d.R. è dunque
⎡X ⎤
⎡x⎤
W = ⎢⎢Y ⎥⎥ , M = ⎢⎢ y ⎥⎥ ⇒ M = RW + T
⎢⎣ Z ⎥⎦
⎢⎣ z ⎥⎦
che, passando in coordinate omogenee, può essere riscritta come:
⎡X ⎤
⎡x ⎤
⎢Y ⎥
⎢ y⎥
= ⎢ ⎥,M
=⎢ ⎥⇒M
= ⎡R T⎤ W
= GW
W
⎢
⎥
⎢Z ⎥
⎢z ⎥
⎣ 0 1⎦
⎢ ⎥
⎢ ⎥
1
⎣ ⎦
⎣1 ⎦
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Trasformazione rigida fra i S.d.R. (2)
Fino ad adesso abbiamo visto che
= A [ I | 0] M
km
• Considerando anche la trasformazione fra i S.d.R. :
= GW
M
= A [ I | 0] GW
km
⎡R T⎤ = A [ I | 0] ⎢
km
W
⎥
⎣ 0 1⎦
• La forma generale di una PPM è quindi
P = A [ I | 0] G ovvero P = A [ R | T]
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Parametri estrinseci
• La matrice G, che modella la posizione della telecamera rispetto alla scena, è
detta matrice dei parametri estrinseci.
• Una generica matrice di rotazione (che ha 9 elementi) è definita da 3
parametri indipendenti, che corrispondono ai valori degli angoli di rotazione
rispetto agli assi del S.d.R. I parametri estrinseci sono quindi 3+3 = 6.
• La PPM più generale tiene quindi conto degli effetti del sensore, tramite A,
della posizione del S.d.R. mondo grazie a G e della proiezione prospettica
incarnata da [ I|0 ].
• Esistono anche effetti dovuti alle distorsioni introdotta dalle lenti che
possono essere modellati tramite un vettore di parametri, detti di distorsione
radiale e tangenziale. Questo vettore comunque non modifica la forma della
PPM, che si basa solo sul modello pin-hole puro.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Influenza della lenti
• La forma così ottenuta della PPM si basa sul modello pin-hole puro.
• In sistemi di acquisizione reali esistono anche effetti dovuti alle distorsioni
introdotte dalle lenti, specialmente per lenti a focale corta e basso costo. La
deviazione più importante è dovuta a fenomeni di distorsione radiale
(“curvatura” delle lenti). Effetti secondari sono introdotti da distorsioni
tangenziali(“decentramento” dei componenti di un sistema di lenti e difetti
di produzione).
• Questi fenomeni sono modellati tramite una relazione non lineare tra i punti
effettivamente “osservati” sul piano immagine e i punti ideali (non distorti)
⎛ x'⎞
⎛ x ⎞ ⎛ dx ⎞
L
r
=
( )⎜ ⎟ + ⎜ ⎟
⎜ ⎟
⎝ y '⎠
⎝ y ⎠ ⎝ dy ⎠
Barrel
dipendente dalla distanza r dal centro
Pincushion
di distorsione ( xc , y c )
r=
( x − xc ) + ( y − yc )
2
Elaborazione dell'immagine M
2
Samuele Salti, Luigi Di Stefano
Coefficienti di distorsione
• La funzione di distorsione radiale L(r) è definita solo per r positivi e L(0) = 1.
Questa funzione non lineare viene tipicamente approssimata tramite il suo
sviluppo di Taylor ( fino all’n-simo ordine, a seconda della precisione voluta )
L ( r ) = 1 + k1r 2 + k2 r 4 + k3 r 6 + ...
• Il vettore di distorsione tangenziale viene invece approssimato con
+ p2 ( r 2 + 2 x 2 ) ⎞
⎛ dx ⎞ ⎛⎜ 2 p1 xy
⎟
⎜ ⎟=
2
2
⎟
⎝ dy ⎠ ⎜⎝ p1 ( r + 2 y ) + 2 p2 xy
⎠
• I coefficienti per la correzione della distorsione radiale k1, k2,…, kn insieme al
centro di distorsione radiale ( xc , y c ) e i due coefficienti di distorsione
tangenziale p1e p2 ampliano e completano l’insieme dei parametri intrinseci del
modello standard di una telecamera. Tipicamente si assume, per semplicità,
che il centro di distorsione radiale coincida con il centro dell’immagine.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
La distorsione ottica nel processo di
formazione dell’immagine
• La distorsione ottica è di solito modellata come una trasformazione che
avviene dopo la proiezione delle coordinate 3D sul piano immagine. Dopo di
che, la matrice degli intrinseci applica una trasformazione affine all’immagine,
traslando coordinate fisiche sul piano immagine in coordinate pixel.
• Il modello di formazione dell’immagine prevede quindi i seguenti passi:
1.
Passare dal sistema di riferimento 3D mondo a quello 3D standard della
telecamera, come indicato dagli estrinseci:
M = Rw + T
2.
Proiettare prospetticamente, ovvero dividere per la terza coordinata:
3.
Applicare le deformazioni dovute al modello delle lenti:
⎛ x'⎞
⎛ x ⎞ ⎛ dx ⎞
L
r
=
(
)
⎜ ⎟
⎜ ⎟+⎜ ⎟
⎝ y '⎠
⎝ y ⎠ ⎝ dy ⎠
4.
Passare dal sistema di riferimento 2D standard a quello in coordinate pixel,
secondo la matrice dagli intrinseci:
T
x = x / z ,
y = y / z
m = A(x'
Elaborazione dell'immagine M
y ')
Samuele Salti, Luigi Di Stefano
Calibrazione (1)
• Abbiamo definito un modello analitico del processo di formazione
dell’immagine.
• Questo modello è rappresentato dalla PPM, che a sua volta può essere
decomposta in tre elementi indipendenti, i parametri intrinseci (la matrice A),
la rotazione R e la traslazione T.
• La calibrazione di un sistema di visione consiste nella stima il più possibile
accurata dei parametri che definiscono questo modello per ogni telecamera
che compone il sistema. A seconda della applicazione è sufficiente stimare
solo la PPM oppure è necessario stimare separatamente le matrici A, R e T.
• Idea alla base di ogni algoritmo di calibrazione: conoscendo la corrispondenza
tra proiezioni 2D e punti 3D di coordinate note, è possibile riscrivere
l’equazione della proiezione prospettica come un sistema lineare con i
parametri come incognite e quindi risolverlo.
• Per ottenere queste corrispondenze, si usano oggetti di forma nota
(scacchiere, pattern ripetitivi, …)
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Calibrazione (2)
• I metodi per la calibrazione possono
essere classificati in due categorie:
– Quelli che usano un’immagine di molti
(almeno 2) piani contenenti un pattern
noto.
– Quelli che usano molte (almeno 3)
immagini diverse di uno stesso
pattern piano.
• Nella pratica è molto più difficile
procurarsi oggetti 3D adatti, con piani
perfettamente ortogonali piuttosto che
un pattern planare che può invece essere
costruito con buona precisione.
• Sono disponibili tools standard per la
calibrazione (OpenCV, Matlab Camera
Calibration Toolbox..)
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Metodo di Zhang (1)
Acquisire n immagini di un
pattern con m corner interni
Per ogni immagine:
stima di una guess iniziale di
omografia Hi
Stima di una guess iniziale per
A a partire dalle Hi
Dato A e Hi stima di una
guess iniziale per Ri e Ti
Stima di una guess iniziale per
i parametri di distorsione k
Raffinamento di ogni
omografia Hi minimizzando
l’errore di riproiezione
Elaborazione dell'immagine M
Raffinamento di A, Ri, Ti, k
minimizzando l’errore di
riproiezione
Samuele Salti, Luigi Di Stefano
Metodo di Zhang (2)
• Sono noti
– il numero di corner interni della scacchiera piana, diversi lungo le due dimensioni.
– la lunghezza del lato dei quadrati che la compongono
• Gli angoli interni della scacchiera possono essere facilmente rintracciati
nell’immagine con algoritmi standard (es. Harris corner detector, eventuale
raffinamento sub-pixel per maggior precisione).
• In ogni frame si fissa il sistema di riferimento 3D nell’angolo in alto a sinistra
della scacchiera con il piano z=0 coincidente con il pattern e gli altri due assi
orientati come la scacchiera, in modo che sia mantenuta sempre la stessa
associazione tra assi e dimensioni della scacchiera (es x = righe, y = colonne).
– La terza coordinata sarà sempre 0.
– La x e la y si ricavano grazie alla lunghezza del lato dei quadrati, che è noto.
• Due passi per derivare i parametri:
– Guess iniziale ottenuta con ottimizzazione lineare (minimizza errore algebrico)
– Raffinamento ottenuto con ottimizzazione non lineare (min errore geometrico)
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Estrinseci ottenuti
• Notare che in ogni frame avremo una diversa stima dei parametri estrinseci.
– Se non interessa un S.d.R. 3D particolare, si può prendere quello associato ad un
frame a caso, p. es. il primo.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
P come omografia
• Con il S.d.R. scelto, in ogni frame si considerano solo punti 3D per cui z=0.
L’equazione di proiezione prospettica diventa quindi:
⎡ p1,1
=⎢ p
= Pw
km
⎢ 2,1
⎢⎣ p3,1
p1,2
p1,3
p2,2
p3,2
p2,3
p3,3
⎡ x⎤
p1,4 ⎤ ⎢ ⎥ ⎡ p1,1
⎥ y
⎢
p2,4 ⎥ ⎢ ⎥ = ⎢ p2,1
⎢0⎥
⎥
p3,4 ⎦ ⎢ ⎥ ⎢⎣ p3,1
⎣1⎦
p1,2
p2,2
p3,2
p1,4 ⎤ ⎡ x ⎤
⎥
′
p2,4 ⎥ ⎢ y ⎥ = Hw
⎢ ⎥
p3,4 ⎥⎦ ⎢⎣ 1 ⎥⎦
• H è nota come omografia, e rappresenta la forma più generale di
trasformazione lineare tra piani. Con w’ si è indicato il vettore (x; y; 1). Si può
pensare H come una semplificazione di P nel caso di oggetti planari.
• Se abbiamo una scacchiera con m corner, possiamo impostare m sistemi di 3
equazioni lineari come quello appena illustrato, in cui le coordinate 3D e 2D
dei punti sono note grazie alle corrispondenze individuate nell’immagine
i-sima e le 9 incognite sono gli elementi di Hi. In realtà, Hi, così come Pi, è
nota a meno di un fattore di scala, e quindi i parametro indipendenti sono
soltanto 8.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Stima di Hi (1)
• A partire dall’equazione dell’omografia, si ottiene
⎡hh1T1 ⎤
′ − h T2 w
′⎤
⎡ vh3T w
⎢ ⎥
⎢
⎥
= Hw
× Hw
′ ⇒ m
′ = 0 ⇒ ⎢ h1T w
′ − uh3T w
′ ⎥ = 0, H = ⎢hhT22 ⎥
km
T
⎢hh3T3 ⎥⎦
⎢uh T2 w
⎥
′
′
ovvero
h
w
v
−
1
⎣
⎣
⎦
⎡ 0T
⎢
T
′
w
⎢
⎢ ′T
⎣ −vw
′T
−w
0T
′T
uw
′T ⎤ ⎡ h1 ⎤
vw
⎥⎢ ⎥
′T ⎥ ⎢h 2 ⎥ = Ah = 0
−uw
0T ⎦⎥ ⎢⎣ h3 ⎥⎦
• Di queste tre equazioni in 9 incognite solo due sono linearmente indipendenti
e tipicamente si mantengono le prime due. Sfruttando tutte le corrispondenze
trovate, si costruisce per ogni immagine un sistema di 2m equazioni in 9
incognite e si ottiene una prima stima di Hi minimizzando un “errore
algebrico” dato dalla norma del vettore Ah ed imponendo l’ulteriore vincolo
che
h = 1 (Algoritmo DLT)
• La soluzione standard del problema precedente si ottiene mediante la
decomposizione ai valori singolari (SVD) della matrice A.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Stima di Hi (2)
• Usando questa stima iniziale, si raffina poi Hi risolvendo ai minimi quadrati
mediante l’algoritmo di Levenberg-Marquardt il problema non-lineare
j − Hi w
′j ,
min ∑ j m
2
Hi
j = 1.....m
• Questo passo di ottimizzazione corrisponde alla minimizzazione dell’errore
ottenuto confrontando, per ogni vertice della scacchiera, le coordinate pixel
effettivamente misurate sull’immagine con la proiezione mediante l’omografia
stimata delle coordinate 3D (con z=0) note. Questo tipo di errore è detto
“errore geometrico”.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Stima dei parametri intrinseci (1)
• Hi è nota a meno di un fattore di scala. La sua relazione con i parametri della
PPM è dunque
H = [ h1 h 2
h 3 ] = [ p1 p 2
p 4 ] = λ A [r1 r2
T]
• R è una matrice ortogonale, i suoi vettori sono ortonormali. Questo vincola i
parametri intrinseci ad obbedire alle seguenti relazioni:
r1Tr2 = 0
⇒ h1T A − T A −1h 2 = 0
r1Tr1 = r2Tr2
⇒ h1T A − T A −1h1 = h 2T A − T A −1h 2
in cui le incognite sono i parametri di B=A-TA-1 . Poiché A è triangolare
superiore, B è simmetrica e le incognite sono solo 6.
• Impilando queste due equazioni calcolate per ogni immagine, si ottiene un
sistema lineare 2n x 6, risolvibile se sono disponibile almeno 3 immagini (vedi
slide successiva per maggiori dettagli)
Elaborazione dell’immagine M
Samuele Salti, Luigi Di Stefano
Stima dei parametri intrinseci (2)
• Fatte le seguenti posizioni
⎛ B11
⎜
B = A − T A −1 = ⎜ B12
⎜B
⎝ 13
B12
B22
B23
hi = [ hi1 , hi 2 , hi 3 ] ,
B13 ⎞
T
⎟
B23 ⎟ , b = [ B11 , B12 , B22 , B13 , B23 , B33 ]
B33 ⎟⎠
v ij T = ⎡⎣ hi1h j1 , hi1h j 2 + hi 2 h j1 , hi 2 h j 2 , hi 3 h j1 + hi1h j 3 , hi 3h j 2 + hi 2 h j 3 , hi 3h j 3 ⎤⎦
si può osservare che
hiT B h j = v ij Tb
h1T Bh 2 = 0
⇒ v12Tb = 0
h1T Bh1 = h T2 Bh 2 ⇒ v11Tb = v 22 Tb ⇒ ( v11 − v 22 ) b = 0
T
• Quindi ogni immagine fornisce 2 equazioni nei 6 parametri indipendenti di B,
e con n immagini di calibrazione si ottiene un sistema lineare omogeneo
risolvibile ai minimi quadrati: Vb = 0
• Dai valori di b è possibile poi ricavare in forma chiusa gli intrinseci, cioè A.
Elaborazione dell’immagine M
Samuele Salti, Luigi Di Stefano
Stima dei parametri estrinseci
• Nota A, e avendo stimato Hi in ogni frame, è possibile derivare prima Ri e poi
Ti per ogni frame.
H i = ⎡⎣h1,i
h 2,i
h k ,i = λ Ark ,i
h3,i ⎤⎦ = λ A ⎡⎣r1,i
r2,i
Ti ⎤⎦
⇒ λrk ,i = A −1h k ,i , k = 1, 2
• Poiché rk,i ha norma unitaria (è un versore), risulta
λ = λrk ,i
• r3 si può ricavare dal vincolo di ortonormalità con r1 ed r2 come
• Infine Ti
r3,i = r1,i × r2,i
Ti =
Elaborazione dell’immagine M
1
λ
A −1h3,i
Samuele Salti, Luigi Di Stefano
Stima dei coefficienti di distorsione (1)
• Data la stima delle omografie, abbiamo sia le coordinate reali (distorte)
provenienti dall’immagine dei punti del pattern sia le coordinate ideali (nondistorte) date dalla proiezione dei punti 3D attraverso le omografie. Il metodo
di Zhang prevede di sfruttare queste informazioni per stimare i coefficienti
k1,k2 della distorsione radiale.
• Data la matrice degli intrinseci (già stimata), A, ed il modello della distorsione
ottica:
⎛ x'⎞
⎛ x ⎞
2
4 ⎛ x⎞
⎜ ⎟ = L ( r ) ⎜ ⎟ = (1 + k1r + k2 r ) ⎜ ⎟
⎝ y '⎠
⎝ y⎠
⎝ y⎠
va innanzitutto determinato il legame fra le coordinate pixel distorte (u’,v’) e
quelle ideali (u, v ) :
⎡ u '⎤
⎡ x ' ⎤ ⎛ αu
⎢ v ' ⎥ = A ⎢ y '⎥ = ⎜ 0
⎢ ⎥
⎢ ⎥ ⎜
⎢⎣1 ⎥⎦
⎢⎣1 ⎥⎦ ⎝⎜ 0
0
αv
0
Elaborazione dell'immagine M
u '− u0
⎧
u0 ⎞ ⎡ x ' ⎤ ⎪ x ' =
αu
⎟⎢ ⎥ ⎪
v0 ⎟ ⎢ y '⎥ → ⎨
v '− v0
⎪
⎟
=
y
'
⎢
⎥
1 ⎠ ⎣1 ⎦
⎪⎩
αv
u − u0
⎧
⎡u ⎤
⎡ x ⎤ ⎪ x =
αu
⎢v ⎥ = A ⎢ y ⎥ → ⎪
⎢ ⎥
⎢ ⎥ ⎨
v − v0
⎢⎣1 ⎥⎦
⎢⎣1 ⎥⎦ ⎪ y =
⎪⎩
αv
Samuele Salti, Luigi Di Stefano
Stima dei coefficienti di distorsione (2)
⎧ u '− u0
− u0 ⎞
2
4 ⎛u
=
1
+
k
r
+
k
r
( 1 2 )⎜ α ⎟
⎪
α
x
'
x
⎡ ⎤
⎪ u
⎝ u ⎠
2
4 ⎡ ⎤
=
1
k
r
k
r
+
+
→
(
)
1
2
⎢ y '⎥
⎢ y ⎥ ⎨
− v0 ⎞
⎣ ⎦
⎣ ⎦ ⎪ v '− v0
2
4 ⎛v
1
k
r
k
r
=
+
+
( 1 2 )⎜ α ⎟
⎪ α
⎝ v ⎠
⎩ v
•
⎧u ' = u + ( k1r 2 + k2 r 4 ) ( u − u0 )
⎪
E’ possibile quindi impostare un sistema in cui
⎨
2
4
v
'
=
v
+
k
r
+
k
r
(
) ( v − v0 )
⎪
1
2
le incognite sono i coefficienti di distorsione radiale: ⎩
⎡( u − u0 ) r 2
⎢
2
−
v
v
r
(
)
0
⎣
( u − u0 ) r 4 ⎤ ⎡ k1 ⎤ ⎡u '− u ⎤
=⎢
⎥
4⎥⎢
( v − v0 ) r ⎦ ⎣ k2 ⎦ ⎣ v '− v ⎥⎦
2
⎛ u − u0 ⎞ ⎛ v − v0 ⎞
r2 = ⎜
⎟ +⎜
⎟
α
α
⎝ u ⎠ ⎝ v ⎠
2
• Disponendo di m punti in n immagini, si può impostare un sistema lineare di
2mn equazioni in 2 incognite, risolvibile ai minimi quadrati.
Dk = d → k = D+ d = ( DT D ) DTd
-1
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Ottimizzazione non-lineare
• La procedura seguita fino a qui minimizza un errore algebrico, senza alcun
significato fisico.
• In presenza di rumore, risulta più robusto raffinare il risultato tramite una
Maximum Likelihood Estimate (MLE) basata su un errore geometrico.
• Facendo l’ipotesi di rumore indipendente e identicamente distribuito, la MLE
per il nostro modello risulta
n
m
∑∑ m
i =1 j =1
ij
ˆ ( A, k , R i , Ti , w j )
−m
2
• Questo problema è risolvibile ai minimi quadrati in modo non lineare grazie
all’algoritmo di Levenberg-Marquardt.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Calibrazione camera stereo
• Fino a qui abbiamo visto come si possono calibrare 2 telecamere
singolarmente, ovvero come ottenere
– Parametri intrinseci e coefficienti di distorsione
– Parametri estrinseci, ovvero la trasformazione rigida dal S.d.R. standard della
camera ad un S.d.R. mondo arbitrario (p.es. relativo al pattern di calibrazione)
• Calibrare un sistema di 2 telecamere contemporaneamente, implica conoscere
in aggiunta a queste informazioni la trasformazione rigida tra una
telecamera e l’altra, R e T.
• In teoria, la conoscenza degli estrinseci rispetto ad uno stesso sistema di
riferimento per tutte le telecamere (p.es. quelli relativi al primo frame per ogni
camera), consentirebbe di ottenere questa informazione, componendo le
trasformazioni Gi(Gj)-1
– N.B. Questo implica camere sincronizzate in HW o pattern fermo!
• Questo metodo è però poco robusto nei confronti del rumore e può portare a
risultati molto lontani dalla struttura reale del sistema di acquisizione.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Calibrazione camera stereo
• Come prima si effettua la stima di una guess iniziale per R e T e poi si
minimizza l’errore di riproiezione.
• Guess iniziale: mediana dei valori di Ri e Ti per ogni coppia di immagini,
ricavati da Gi(Gj)-1.
– Ri convertito in un vettori di 3 elementi con la trasformata di Rodriguez.
• Quindi ottimizzazione non lineare con l’algoritmo di Levenberg-Marquardt
minimizzando l’errore di riproiezione in entrambe le immagini
n
m
∑ ∑∑ m
k = L , R i =1 j =1
Elaborazione dell'immagine M
k
ij
ˆ ( A, k , R, T, w j )
−m
2
Samuele Salti, Luigi Di Stefano
SdR di un sistema stereo
• Ottenuta la calibrazione di un sistema stereo, si opera tipicamente una scelta
standard per la definizione del sistema di riferimento mondo della coppia di
camere.
• In particolare si sceglie il sistema di riferimento standard di una delle due
camere, per esempio la sinistra. Gli estrinseci della camera destra sono così
rappresentati dalle matrici di rotazione e traslazione appena calcolate. Le PPM
corrispondenti sono quindi
P L = A L ⎡⎣I 0 ⎤⎦ P R = A R ⎡⎣ R T ⎤⎦
z = zL
fL
cL
zR
R, T
x = xL
fR
xR
cR
Samuele Salti, Luigi Di Stefano
Rettificazione
•
•
•
•
•
•
Quando si lavora con sistemi stereo, è
normale fare l’assunzione di lavorare
con un sistema stereo laterale.
È impossibile ottenere questa configurazione tramite allineamento meccanico.
Se si è calibrata la coppia stereo è però
possibile rettificare via software le
immagini acquisite.
L’idea che sta dietro la rettificazione `e quella di definire due nuove PPM, P1e P2,
ottenute ruotando le matrici originali attorno ai loro centri ottici finché i piani
focali non diventano co-planari (e quindi anche i piani immagine).
Inoltre, per garantire che i punti coniugati abbiano le stesse coordinate verticali, si
impone che le nuove PPM abbiano gli stessi parametri intrinseci.
Tutto il processo si basa sul modello pin-hole puro: le trasformazioni calcolate
devono essere applicate su immagini da cui siano state rimosse le
distorsioni.
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Nuove PPM
• Noti Ai ,Ri e il centro ottico ci per entrambe le camere si procede come segue:
– Si assume una nuova matrice degli intrinseci Anew arbitraria: p. es. la media
delle Ai
– Si ricava la nuova matrice Rnew, che determina l’orientazione della telecamera,
anch’essa identica per entrambe le PPM. Tenendo presente che i vettori riga di
Rnew rappresentano gli assi X, Y, e Z, del sistema di riferimento della telecamera,
espressi in coordinate mondo, si può ottenere un sistema laterale imponendo
• nuovo asse X parallelo alla linea di base
• nuovo asse Y perpendicolare a nuovo X e
versore arbitrario k
• nuovo asse Z perpendicolare a X e Y
r1 = ( c 2 − c1 ) / c 2 − c1
r2 = k ∧ r1
r3 = r1 ∧ r2
k è un versore arbitrario, tipicamente il vecchio asse Z di una delle due telecamere.
• Le nuove PPM sono dunque A [ R
new
new
Elaborazione dell'immagine M
− R new c1 ] , A new [ R new
− R new c 2 ]
Samuele Salti, Luigi Di Stefano
Il centro ottico
• Si può dimostrare che ogni PPM è definita da una matrice di rango massimo,
pari a 3, e viceversa ogni matrice 3x4 di rango massimo definisce una
proiezione prospettica.
– Intuitivamente : il rango massimo è necessario, altrimenti la proiezione di un
punto non sarebbe un punto ma un sottospazio di dimensione maggiore (retta,
piano, …)
• Ricordando che (0,0,0) non rappresenta alcun punto in coordinate omogenee,
segue che un punto che appartiene al kernel (soluzioni del sistema omogeneo
associato) della matrice è un punto per cui la proiezione non è definita.
• Dal rango massimo segue che il kernel di questa matrice ha dimensione uno
nello spazio omogeneo, ovvero contiene un unico punto effettivo. Il centro
ottico è infatti l’unico punto dello spazio per cui la proiezione non è definita,
ed ha quindi equazione:
P ⎡c ⎤ = 0 ovvero c = − P −1p
4
⎢ ⎥
⎣1⎦
Elaborazione dell'immagine M
con P = ⎡⎣ P p 4 ⎤⎦
Samuele Salti, Luigi Di Stefano
Scomposizione della PPM
• I metodi di calibrazione visti forniscono già A, R e T separatamente.
• È comunque sempre possibile ricavare gli elementi che compongono la PPM.
Questo è possibile perché ogni matrice quadrata non singolare si può
decomporre nel prodotto tra una matrice ortogonale e una triangolare
superiore tramite la fattorizzazione QR.
−1
• Sia dunque P = UL la fattorizzazione QR dell’inversa di P, con U
ortogonale e L triangolare superiore.
• Poiché
P = ⎡⎣ P p 4 ⎤⎦ = ⎡⎣ AR AT ⎤⎦ ⇒ P −1 = R −1A −1 ⇒
R = U −1 , A = L−1
• Infine, per ricavare T
T = A −1p 4 = Lp 4
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Trasformazioni di rettificazione
• Per rettificare un’immagine
dobbiamo calcolare la
trasformazione che mappa il piano
immagine della telecamera reale in
quello della telecamera rettificata.
• Poiché la trasformazione che
abbiamo imposto attraverso il
calcolo delle nuove PPM
rappresenta una rotazione intorno
al centro ottico, possiamo scrivere
le equazioni dei raggi ottici come
o1 λo ∈ R, P o1 = [ Po1 p o1 ]
⎧⎪ w = c1 + λo Po−11m
⎨
−1
n1 λn ∈ R, P n1 = [ Pn1 p n1 ]
⎪⎩w = c1 + λn Pn1 m
da cui si ricava poi l’omografia di
rettificazione
n1 = λ Pn1Po−11m
o1 λ ∈ R
m
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Esempio
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano
Principali Riferimenti Bibliografici
1) R. Hartley and A. Zisserman, “Multiple View Geometry in Computer
Vision”, 2° Edition, Cambridge University Press, 2003.
2) Z. Zhang, “A flexible new technique for camera calibration”, IEEE Trans.
On Pattern Analysis and Machine Intelligence, 22(11):1330-1334, November
2000.
3) A. Fusiello, “Visione Computazionale – Appunti delle lezioni”,
http://profs.sci.univr.it/~fusiello
Elaborazione dell'immagine M
Samuele Salti, Luigi Di Stefano