Università degli Studi di Roma “Sapienza”
Facoltà di Ingegneria
Laurea Magistrale in Ingegneria Informatica
Tesina del Corso di Visione e Percezione A.A. 2008-09
Studenti:
Luciano Blasi
Redjan Shabani
matr. 321716
matr. 1013173
Scopo del lavoro:
Stimare una nuvola di punti 3D di una scena e la posa della camera
che effettua le riprese, prendendo in esame corrispondenze di
punti tra immagini successive.
Sequenza di elaborazione:
1) Lettura della sequenza immagini ed estrazione delle features per ogni immagine
2) Calcolo delle corrispondenze tra features in coppie di immagini
3) Calcolo delle matrici fondamentali associate ad ogni coppia di immagini
4) Calcolo delle matrici delle camere in base alle stime delle matrici fondamentali
5) Triangolazione dei raggi dei punti corrispondenti per ottenere i punti dello spazio
6) Correzione degli errori sulle stime delle matrici delle camere e sulle stime dei punti
3D e ricostruzione proiettiva.
7) Passaggio dalla ricostruzione proiettiva alla ricostruzione affine
8) Passaggio dalla ricostruzione affine alla ricostruzione metrica
Riferimenti:
[1] Hartley R , Zisserman “A Multiple View Geometry In Computer Vision”, 2Ed.,2003
[2] A. Fusiello “Visione Computazionale”, Appunti delle lezioni ed. 2009
[3] Peter Kovesi Peter’s functions for Computer Vision
1) Lettura della sequenza immagini ed estrazione delle
features per ogni immagine
Le immagini sono state prese con pose camera “vicine” al fine di ottenere
features che fossero più facilmente coniugabili tra immagini successive.
Le features sono rilevate con il seguente algoritmo di Harris
(§ Alg. 9.2 in [2])
1) Calcolo delle derivate dell’immagine:
2) Calcolo di
e loro filtraggio con gaussiana
3) Per ogni punto calcolano la matrice delle misure di Harris:
4) Si calcolano le quantità :
5) Si effettua la soppressione dei non massimi locali facendo scorrere una
finestra centrata sul punto in esame.
6) Sogliatura sulla base di una frazione del valore massimo.
7) I punti risultanti sono i corner, punti di interesse.
L’implementazione (extractCorners.m) riceve un set di immagini
e restituisce gli insiemi di corner
Nota: in questo algoritmo l’operatore deve scegliere due parametri che hanno
impatto da valutare sulla qualità delle features estratte:
1) Grandezza del kernel gaussiano (o sua deviazione standard) che fissa la
scala spaziale alla quale si vogliono determinare i punti salienti
2) Soglia o numerosità dei punti salienti ottenuti.
2) Calcolo delle corrispondenze tra features in coppie di
immagini
La ricerca delle corrispondenze tra features su immagini diverse, è
basata sulla similarità.
L’implementazione (correlation.m) riceve in input le immagini
e i relativi set di corner
e restituisce i set di corrispondenze nella forma
dove
contiene le coppie corrispondenti nelle immagini
e
La funzione usa il codice matchbycorrelation.m descritto in [3],
Risultato del Harris-Corner-Detector applicato a due immagini consecutive
Risultato del calcolo delle corrispondenze dalle
precedenti immagini.
3) Calcolo delle matrici fondamentali associate ad ogni
coppia di immagini
La matrice fondamentale descrive tutta la geometria epipolare del sistema
definito dalla coppia di camere che osservano la scena. Vale:
La matrice fondamentale viene calcolata con l’algoritmo degli 8 punti
(§ Alg. 11.1 di [1]). Dati in input 8 o più corrispondenze di punti provenienti
da due immagini correlate, si eseguono i seguenti passi:
1) Normalizzazione: si trasformano le coordinate immagine x in base a
Dove la
è una matrice di traslazione e scalatura che
sposta il baricentro del set dei punti, nell’origine e scala le
coordinate degli
in modo che la somma delle
loro norme sia
.
2) Calcolo di
. Il calcolo di F passa attraverso la soluzione del
seguente sistema lineare:
Trovate le soluzioni del sistema con la decomposizione
SVD otteniamo:
In generale il
a causa del rumore. Superiamo il
problema decomponendo in SVD la
e ponendo il più piccolo
valore singolare a 0.
3) Denormalizzazione. La matrice
finale è:
L’algoritmo degli 8 punti è implementato in fundmat8p.m.
In letteratura non è consigliato andare avanti con la F così
calcolata. Si controlla per ogni match la relazione:
Se il risultato è “molto” differente dallo zero si segna come
outlier e si scarta dai match di interesse.
La funzione ranscfitfundmatrix.m di P.Kovesi [2], calcola la
matrice fondamentale degli inlier con l’ausilio di RANSAC.
Tale funzione richiede anche un parametro soglia che è
determinante per distinguere inliers da outliers.
Lo script multipleFundMat.m
è la versione plurale di
ranscfitfundmatrix.m. In
input richiede una serie di
set di corrispondenze.
Il calcolo restituisce una
serie di matrici fondamentali:
Visualizzazione dei feature correlati e raffinati
dopo il calcolo della matrice fondamentale
4) Calcolo delle matrici delle camere in base alle stime delle
matrici fondamentali
Le matrici di proiezione prospettica delle camere, definiscono la
relazione che sussiste tra punti 3D nel mondo e immagine di tali
punti nella camere considerate.
Date due matrici P,P’ associate ad una coppia di camere, la
matrice fondamentale che le lega è unica ed è data da:
dove
è la pseudo-inversa di P,
, con
Data la matrice F la coppia di camere ad essa associata è:
.
Nel caso di N viste successive è necessario calcolare le N matrici
delle camere a partire dalle matrici fondamentali. Per ogni
si calcolano le coppie di camere in forma canonica
.
Il sistema di riferimento dell’i-esima coppia sta nell’i-esima
camera.Bisogna quindi rototraslare per spostare il sistema di
riferimento nella prima camera. La funzione cameras.m riceve
in input un insieme di matrici fondamentali e restituisce l’insieme
delle coppie di matrici delle camere come descritto in precedenza.
5) Triangolazione dei punti corrispondenti per ottenere
i punti dello spazio
Con le matrici delle camere, si possono stimare i punti 3D dello
spazio dati i punti nelle immagini. Conoscendo l’immagine x di un
punto dello spazio X e la posa della camera, si costruisce il raggio
passante per il centro della camera ed il punto x. Analogamente si
costruisce il raggio associato a x’. L’intersezione dei due raggi è il
punto 3D. Dato che le immagini del punto non sono quelle del caso
ideale, i raggi difficilmente si intersecheranno tra loro. Si pone cosi
il problema di quale punto dello spazio scegliere come proiezione
delle due immagini. In [2] si descrive uno dei possibili metodi per
la triangolazione, chiamato metodo linear-eigen. Le matrici P e P’
vengono scritte nella seguente forma:
mettendo cosi in evidenza i vettori riga di
ciascuna delle due. Partendo dalle relazioni
di proiezione prospettica, si arriva alla relazione:
Nel caso ideale la soluzione e data dal nucleo
della matrice dei coefficenti. In realta i dati sono
affetti da rumore e si è costretti a determinare la
soluzione per mezzo della decomposizione SVD.
Ricostruzione proiettiva basata sulla coppia di foto viste prima
6) Correzione degli errori sulle stime delle matrici delle
camere e sulle stime dei punti 3D e ricostruzione proiettiva.
In presenza di molteplici viste l’errore in fase di triangolazione si cumula, sia per la
stima fatta alle matrici delle camere sia nel calcolo dei punti nello spazio proiettivo.
E’ necessario raffinare i risultati minimizzando l’errore nell’immagine. La minimizzazione
dell’errore viene fatta intervenendo sia sulle coordinate dei punti ottenuti per
triangolazione, sia sui parametri delle telecamere. La procedura e chiamata bundle
adjustment, descritta in [2]. Si cerca di spostare sia le N fotocamere che gli n punti 3D
affinchè la somma delle distanze tra il j-esimo punto riproiettato nell’i-esima fotocamera
ed il punto misurato sia il più piccolo possibile:
Nella Ri sarà necessario parametrizzare in modo che compaiano solo tre incognite
(es. gli angoli di Eulero), invece che tutti i nove elementi. Sempre in [2] si consiglia
applicare una strategia in due passi alternati:
Prima si tengono fermi i punti e si risolve rispetto alle telecamere e poi si fissano le
telecamere e si calcolano i punti 3D. Finita la procedura di minimizzazione degli errori,
le camere
e i punti
costituiscono la ricostruzione
proiettiva. Più precisamente, si dispone di una stima di struttura e moto a meno
di una trasformazione proiettiva.
Da ricostruzione proiettiva ad affine
Dati di ingresso, ricostruzione proiettiva in P3:
Dati di uscita, ricostruzione affine in A3:
(Pp,P’p,{Xpi})
(Pa,P’a,{Xai})
Le due ricostruzioni sono legate da un’opportuna trasformazione
proiettiva H, incognita, ma costante per tutti i punti.
Un piano che si trova all’infinito nella “vera” ricostruzione1,
viene rappresentato in P3 (nella nostra ricostruzione proiettiva)
da un 4-vettore omogeneo finito ∏. Trovare questo vettore ci
consente di determinare H, infatti:
1) In qualche modo determino le coordinate di questo ∏ nella mia
ricostruzione proiettiva P3
2) ∏ in A3 ha posizione canonica (0,0,0,1)T
3) Una H trasforma il piano ∏ in: H-T. ∏
4) Trovare H tale che H-T. ∏= (0,0,0,1)T, quindi
 I | 0
H 




5) Applico H trovato a (Pp,P’p,{Xpi}) e trovo:
Pa=PpH-1, Pa’=Pp’H-1, Xai=HXpi.
Nei problemi di ricostruzione derivata da dati reali, esiste una “vera” ricostruzione che consiste nei
veri punti Xi e nelle vere camere P,P’ che hanno generato le osservazioni misurate. L’insieme dei punti
ricostruiti X’i differiscono dagli Xi da un’opportuna trasformazione (proiettiva nel presente caso).
1
1)
Determinazione delle coordinate di ∏ in P3
Il metodo per determinare le coordinate di ∏ in P3 cambia a seconda
di quali informazioni vogliamo (o possiamo) sfruttare (§ pag.268-271 di [1]).
Ci limitiamo a considerare qui il caso in cui l’informazione che
vogliamo sfruttare è il parallelismo di alcune linee della scena.
Selezioniamo “opportunamente” alcune coppie di linee che sappiamo
essere parallele nella realtà e determiniamo le coordinate in P3 del
loro punto di fuga V (prima le coordinate nelle immagini v e v’ e poi
triangolando con le P,P’ note, otteniamo le coordinate in P3 di V - §
pag.60 di [2]).
L’”opportunamente” evidenzia il fatto che i punti V non devono
essere allineati in P3 per poter determinare ∏ (§ eq. 3.3 di [1]).
E’ da notare che non è necessario trovare i punti di fuga in entrambe
le immagini. Infatti chiamando v un punto di fuga individuato nella
prima immagine e l’ la corrispondente linea nella seconda immagine,
dal momento che i punti di fuga soddisfano i vincoli epipolari posso
calcolare il corrispondente punto di fuga v’ nella seconda immagine
come l’intersezione di l’ e la linea epipolare Fv di v. La costruzione
del punto V dello spazio P3 può essere espresso algebricamente come
la soluzione delle due equazioni: ([v]x.P)V = 0 e (l’.P’)X = 0
queste equazioni esprimono il fatto che V mappa v nella prima
immagine e ad un punto su l’ nella seconda immagine.
2)
∏ in A3 ha posizione canonica (0,0,0,1)T
Dimostriamo ora che il vettore (0, 0, 0, 1)T rappresenta
il piano all’infinito, cioè quel piano contiene tutti e soli i
punti (vettori) ideali.
1)
Un punto all’infinito è rappresentato da un
vettore del tipo (x1, x2, x3, 0)T.
2)
La condizione che un punto appartenga ad un
piano è che sia nullo il prodotto scalare tra punto
e piano (che equivale a scrivere per esteso
l’equazione omogenea del piano passante per il
punto)
3)
(0,0,0,1)T.(x1, x2, x3, 0) = 0
c.d.d.
3)
Una H trasforma il piano ∏ in: H-T. ∏
Come dimostrato nella slide precedente, l’appartenenza del
punto x al piano ∏ si scrive:
∏.x=0
Sia H una matrice omogenea 4x4 non-singolare.
E’ immediato verificare che ∏.H-1.H.x= ∏.x=0
Quindi i punti H.x giacciono tutti sul piano ∏.H-1 (piano che
si può anche scrivere come H-T.∏)
Quindi se
H trasforma un punto x in punto Hx
H trasforma un piano ∏ in piano H-T.∏
4) Trovare H tale che H-T. ∏= (0,0,0,1)T
E’ immediato verificare che con
vale:
quindi:
 I | 0
H 

  
HT . (0,0,0,1) = ∏
H-T . ∏= (0,0,0,1)T
ATTENZIONE:
questa formula non funziona se la coordinata finale di ∏ è zero (§ A4.2-pag.580 di [1]).
Da ricostruzione affine a metrica
Dati in ingresso, ricostruzione affine in A3:
(Pa=[M|m],Pa'=[M'|m'],{Xai})
Dati di uscita, ricostruzione metrica in E3: (Pe,P’e,{Xei})
I passi dell’algoritmo, sono:
1) In qualche modo individuo su qualche immagine la conica
assoluta w
2) Calcolo A, partendo da w e M con la fattorizzazione di
Cholesky dall'equazione: AAT = (MTwM)-1 e trovo H
 A1 
H 

1

3) Applico H trovato a (Pa,P’a,{Xai}) e trovo:
Pe’=Pa’H-1, Xei=HXai.
Pe=PaH-1,
1) Individuazione della conica assoluta w
Per determinare i coefficienti della matrice simmetrica w che
rappresenta la conica (assoluta) si usano combinazione delle
possibili tre fonti di vincoli seguenti:
1) Vincoli che derivano da ortogonalità sulle scene:
Un esempio è il punto di fuga per la direzione verticale (una
persona in piedi) e la linea di fuga del piano terra orizzontale (tra
loro subortogonali). Oppure un piano della scena fotografato
contenente informazioni metriche come una griglia quadrata.
2) Vincoli che provengono dalla conoscenza dei parametri interni.
Se la matrice di calibrazione di una camera è uguale a K,
allora l’immagine della conica assoluta è w = K-T.K-1
3) Vincoli che si presentano dalle stesse camere in tutte le
immagini (viste multiple).
Una delle proprietà della conica assoluta è che la sua proiezione
nelle immagini dipende solo dalla matrice di calibrazione della
camera e non dalla posa della camera.