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 matr. 321716Redjan Shabani matr. 1013173

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

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

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.

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 delleloro norme sia .

2) Calcolo di . Il calcolo di F passa attraverso la soluzione delseguente sistema lineare:

Trovate le soluzioni del sistema con la decomposizioneSVD otteniamo:

In generale il a causa del rumore. Superiamo ilproblema decomponendo in SVD la e ponendo il più piccolovalore 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 comeoutlier e si scarta dai match di interesse. La funzione ranscfitfundmatrix.m di P.Kovesi [2], calcola lamatrice 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 diranscfitfundmatrix.m. In input richiede una serie di set di corrispondenze. Il calcolo restituisce una serie di matrici fondamentali:

Visualizzazione dei feature correlati e raffinatidopo 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 ognisi 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 diriferimento nella prima camera. La funzione cameras.m ricevein input un insieme di matrici fondamentali e restituisce l’insiemedelle 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 dellospazio 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 perla 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 nucleodella matrice dei coefficenti. In realta i dati sonoaffetti da rumore e si è costretti a determinare lasoluzione 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 fotocameraed 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 ricostruzioneproiettiva. Più precisamente, si dispone di una stima di struttura e moto a menodi una trasformazione proiettiva.

Da ricostruzione proiettiva ad affine

1 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).

Dati di ingresso, ricostruzione proiettiva in P3: (Pp,P’p,{Xpi}) Dati di uscita, ricostruzione affine in A3: (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

5) Applico H trovato a (Pp,P’p,{Xpi}) e trovo: Pa=PpH-1, Pa’=Pp’H-1, Xai=HXpi.

0|I

H

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 ilpunto)

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=0Quindi 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 HxH trasforma un piano ∏ in piano H-T.∏

4) Trovare H tale che H-T. ∏= (0,0,0,1)T

E’ immediato verificare che con

0|I

H

vale: HT . (0,0,0,1) = ∏

quindi: 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

1

1AH

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

3) Applico H trovato a (Pa,P’a,{Xai}) e trovo: Pe=PaH-1, Pe’=Pa’H-1, Xei=HXai.

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.