1

DEFORMAZIONI PROSPETTICHE

PARTE I: trasformazioni geometriche

1.1 Introduzione Abbiamo analizzato e cercato di correggere il problema delle distorsioni geometriche a cui sono

soggette le immagini ottenute tramite macchine fotografiche o fotocamere digitali quando il punto

di ripresa della scena non risulta frontale al soggetto ripreso.

Se si vuole, per esempio, fotografare un quadro appeso ad una parete o la facciata di un edificio è

opportuno che il piano della fotocamera su cui si forma l’ immagine sia esattamente parallelo alla

tela e alla parete, altrimenti si ottiene una foto in cui il soggetto, che nella realtà ha forma

rettangolare, risulta restringersi verso l’alto oppure avere la forma di un quadrilatero irregolare

come illustrano le seguenti foto.

Figura 1. 1 Foto affetta da deformazione del trapezio

2

Figura 1. 2Foto affetta da deformazione prospettica

In questi casi è utile poter modificare l’immagine in modo da simulare un punto di ripresa rispetto

al quale si può ottenere la corretta geometria della scena. Si vuole operare tale correzione mediante

una opportuna trasformazione geometrica della immagine che non ne alteri le caratteristiche

radiometriche.

Poichè si è lavorato con immagini bitmap a 24 bit si sono prese in esame trasformazioni che

operano sui singoli pixel dell’immagine. Affinchè tali trasformazioni siano accettabili è necessario

che esse preservino la continuità delle linee, i reciproci rapporti di posizione e proporzione degli

oggetti e non ne alterino le forme1.

1.2 Aspetti geometrici del processo di formazione d elle immagini

Per correggere le distorsioni geometriche a cui sono soggette le fotografie scattate da un punto di

ripresa non ottimale è utile analizzare la relazione geometrica tra la scena osservata e quella della

corrispondente immagine ottenuta con un apparecchio fotografico. Tale relazione si basa sulla

geometria proiettiva e sulla prospettiva.

Si può schematizzare la macchina fotografica come un sistema costituito da una camera oscura e

dotato di una apertura nella quale è posta una lente biconvessa. Questa focalizzando i raggi

luminosi provenienti dall’esterno fa formare sulla parete opposta della camera oscura l’immagine

degli oggetti inquadrati che viene poi fissata da una pellicola fotografica .

1 Cioè non devono modificare la forma degli elementi della scena in un modo che non rientra nella logica della variazione del punto di ripresa.

3

y,Y piano immagine

x,X punto oggetto

(X,Y,Z)

z,Z

F

(x,y) centro della lente

punto immagine

Figura 1. 3 Modello del sistema di formazione di un’immagine. Il sistema di coordinate x,y,z costituisce il sistema di coordinate fisso della macchina.

Il piano xy ove si forma l’immagine coincide con il fondo della camera oscura e il punto F

rappresenta il centro della lente biconvessa di distanza focale f. Il sistema di coordinate xyz

costituisce il sistema di coordinate fisso della macchina, l’asse Z coincide con l’asse ottico

passante per il centro del piano di formazione dell’immagine e il centro F della lente.

Le relazioni tra gli oggetti dello spazio reale e la loro immagine sono regolate dalle leggi

dell’ottica geometrica, secondo cui si suppone che i raggi luminosi si propaghino in linea retta .

Un generico punto dello spazio P (punto oggetto) di coordinate (X,Y,Z) viene proiettato attraverso

il punto F in un punto P’ del piano dell’immagine di coordinate (x,y) .

Dalle similitudini dei triangoli si ottiene che

fXx

f Z

fYy

f Z

=−

=−

(1.1)

in cui il punto proiettato non dipende in modo lineare dal punto oggetto. Per semplificare tali

relazioni si è preferito passare alle coordinate omogenee.

Sia v il vettore contenente le coordinate del punto oggetto

4

X

Y

Z

=

v

Il vettore omogeneo v’ corrispondente di v é

'

sX

sY

sZ

s

=

v

dove s reale è la costante di scala.

La matrice rappresentativa della proiettività centrata in F è

1 0 0 0

0 1 0 0

0 0 1 0

10 0 1

f

= −

A (1.2)

per cui si ha

' '=w Av (1.3)

con

'

sX

sY

sZ

sZs

f

= −

w

normalizzando w’ si ottiene

fX

f Z

fY

f Z

fZ

f Z

−

= − −

w

vettore in cui le prime due componenti risultano uguali alle relazioni (1.1) .

La proiezione inversa che da un punto immagine P’ di coordinate (x,y) torna indietro al punto

oggetto corrispondente P è data da

5

' '= -1v A w (1.4)

con

1 0 0 0

0 1 0 0

0 0 1 0

10 0 1

f

−

=

1A (1.5)

'

sx

sy

sz

s

=

w

in cui z é una variabile libera. Calcolando quindi la trasformazione prospettica inversa si ottiene in

coordinate omogenee

'

sx

sy

sz

szs

f

= +

v

Le corrispondenti coordinate cartesiane sono

fx

f z

fy

f z

fz

f z

−

= − −

v (1.6)

Oppure

fxX

f z

fyY

f z

fzZ

f z

= −

= −

= −

(1.7)

Risolvendo rispetto a z e sostituendo si ottiene

6

( )

( )

xX f Z

f

yY f Z

f

= − = −

(1.8)

Da queste relazioni si rileva immediatamente che, poichè z è arbitrario, il fibrato del punto

immagine (x,y) tramite la proiezione è costituito da tutti i punti della retta congiungente il punto

immagine con il fuoco F che hanno Z>f. Il processo di proiezione non è univoco e tutti i punti

della retta PF vengono proiettati nel punto P’. E’ quindi necessario specificare una delle coordinate

di un punto oggetto, per esempio Z, per determinare le altre due a partire da quelle del punto

immagine (x,y).

Inoltre l’oggetto osservato deve trovarsi ovviamente anteriormente alla macchina fotografica ossia

Z deve essere tale che Z>f.

Una posizione non ottimale del piano di formazione dell’immagine rispetto alla scena che si vuole

fotografare può dare origine nella foto a deformazioni geometriche che alterano le forme degli

oggetti fotografati.

1.3 Modello matematico discreto dell’immagine

Si è lavorato con immagini digitali bitmap a 24 bit adatte ad essere manipolate da un calcolatore.

Da un punto di vista matematico esse possono essere approssimate con una matrice a valori interi

del tipo:

(1,1) (1,2) (1,3) (1, )

(2,1) (2,2) (2,3) (2, )

( , )

( ,1) ( ,2) ( ,3) ( , )

G G G G K

G G G G K

j k

G J G J G J G J K

=

G

⋯

⋯

⋮ ⋮ ⋮ ⋮

⋮ ⋮ ⋮ ⋮

⋯

(1.9)

Con

1

1

j J

k K

≤ ≤≤ ≤

J indica il numero di righe di cui l’immagine è composta, K è il numero di pixel per ogni riga.

7

Nelle immagini monocromatiche l’intero G(j,k) è un valore compreso tra 0 e 255 e rappresenta

uno dei 256 livelli di grigio cioè di intensità luminosa

{ }( , ) : 0 255G j k I n n∈ = ∈ Ν ≤ ≤ (1.10)

in quelle a colori invece a ciascun elemento G(i,j) della matrice corrisponde una terna ordinata di

interi (r,g,b)

( , ) ( , , )G j k r g b=

ove r,g,b indicano rispettivamente le intensità di rosso, verde e blu la cui combinazione produce il

colore e assumono valori interi compresi tra 0 e 255 corrispondenti ai 256 diversi livelli di

intensità delle tre sorgenti di colore.

( , ) r g bG j k I I I∈ × × (1.11)

con

{ }{ }{ }

: 0 255

: 0 255

: 0 255

r

g

b

I n n

I n n

I n n

= ∈ Ν ≤ ≤

= ∈ Ν ≤ ≤

= ∈ Ν ≤ ≤

(1.12)

1.4 Trasformazioni geometriche

Sia G(j,k) con 1 j J≤ ≤ e 1 k K≤ ≤ l’immagine di output generata da una operazione geometrica

sull’immagine discreta di input F(p,q) con 1 1p P q Q≤ ≤ ≤ ≤ che operi sulla geometria della

scena mantenendone inalterate le caratteristiche radiometriche.

In generale F(p,q )e G(j,k) possono avere dimensioni diverse.

Quantitativamente tale trasformazione è definita con le formule

1

2

( , ) ( , )

( , )

( , )

G x y F p q

x t p q

y t p q

===

(1.13)

La prima formula indica l’uguaglianza dell’intensità luminosa nei punti corrispondenti delle due

immagini, le ultime due individuano le trasformazioni tra le coordinate dei punti che si

corrispondono.

8

E’ molto importante osservare che in generale la coppia (x,y), ottenuta applicando le

trasformazioni t1 e t2 alla generica coppia (p,q), non è una coppia di interi corrispondente ad un

elemento della matrice G(j,k).

Se immaginiamo l’immagine discreta come una griglia rettangolare2 i cui nodi corrispondono agli

elementi della matrice che la approssima, le trasformazioni t1 e t2 deformano la griglia

dell’immagine di input in una griglia irregolare come quella mostrata a tratto continuo nella figura

di destra.

Figura 1. 4 Il principio delle trasformazioni geometriche. La figura di destra è l’immagine di input con griglia quadrata regolare, quella di sinistra a tratto continuo è l’immagine trasformata con griglia irregolare, la griglia regolare tratteggiata individua i punti dell’immagine di output di cui si deve conoscere l’intensità.

La costruzione dell’immagine di output non si può ottenere direttamente con i nodi di quest’ultima

poichè gli elementi della matrice G corrispondono ai nodi della griglia tratteggiata ed è quindi

necessario che siano definiti i valori dell’intensità luminosa in questi punti.

Nella pratica per ovviare a questo problema si utilizzano le trasformazioni geometriche inverse di

t1 e t2 che indichiamo con t3 e t4

3

4

' ( , )

' ( , )

x t j k

y t j k

==

(1.14)

Ad ogni coppia ordinata (j,k) con 1 j J≤ ≤ e 1 k K≤ ≤ corrispondente ad un vertice della griglia

tratteggiata (o all’elemento G(j,k) della matrice) si associa la coppia di coordinate (x’,y’) . In

generale anche quest’ultima non corrisponde ad un nodo della griglia dell’immagine di input ed è

quindi necessario assegnargli un opportuno3 valore attraverso un processo di interpolazione.

2 Col termine griglia rettangolare indichiamo l’insieme compatto di R2 [1,J]×[1,K] in cui è evidenziato il sottoinsieme dei punti (x,y) con x o y intero. Chiamiamo nodo o vertice della griglia ogni punto della griglia (x,y ) con x e y interi. 3 Nel capitolo seguente si specificherà meglio questo punto.

9

Finora si è solo specificato che le trasformazioni geometriche coinvolte operano sui singoli pixel:

un punto dell’immagine di imput viene mappato in un punto dell’immagine di output. Affinchè tali

trasformazioni siano accettabili è anche necessario, come si è accennato nell’introduzione, che

esse preservino la continuità delle linee, i reciproci rapporti di posizione e proporzione degli

oggetti e non ne alterino le forme4.

E’ inoltre importante osservare che in generale tali trasformazioni sono definite basandosi su

sistemi di riferimento di coordinate cartesiane che non coincidono con il sistema di indici

dell’array associato all’immagine discreta ed è quindi necessario di volta in volta stabilire le

relazioni che legano questi due riferimenti.

Nella figura 1.5 è illustrato per esempio il caso in cui il riferimento cartesiano ha l’origine (0,0)

posizionata nell’angolo in basso a sinistra dell’immagine, mentre per l’array è l’angolo in alto a

sinistra di indici (1,1) che funge da riferimento.

k 1 2 K 1

j

yj

J

xk

Figura 1. 5 Esempio di immagine in cui il riferimento cartesiano (xk,yj) ha origine nell’angolo in basso a sinistra, mentre il riferimento (1,1) dell’array è l’angolo in alto a sinistra.

Le relazioni tra i due sistemi sono date da

1

21

2

k

j

x k

y J j

= −

= + − (1.15)

Più in generale si scrive

4 Cioè non devono modificare la forma degli elementi della scena in un modo che non rientra nella logica della variazione della prospettiva

10

( )

( )j

k

y u j

x v k

= =

(1.16)

con u e v biettive.

Le trasformazioni geometriche elementari delle immagini sono le traslazioni, le rotazioni, i

cambiamenti di scala e le proiezioni.

1.4.1 Traslazione

La traslazione a corpo rigido di un’immagine implica che ogni suo punto venga spostato di una

certa distanza in una direzione fissata.

Si può quantificare la traslazione con un vettore t

x

y

t

t

Se si indicano con ( , )q pu v le coordinate cartesiane dell’immagine discreta di input F(p,q) e con

( , )k jx y quelle dell’immagine di output G(j,k) si ha

k q x

j p y

x u t

y v t

= +

(1.17)

Le relazioni tra gli indici degli array delle due immagini discrete si ottengono componendo la

traslazione con le relazioni che legano tali indici con le corrispondenti coordinate cartesiane.

Supponendo valide le equazioni (1.4) si ha

' ( )

'y

x

j p P J t

k q t

= − − −

= + (1.18)

Il simbolo di primo indica che, a meno che xt e yt non siano valori interi, j’ e k’ possono non

essere interi. E’ quindi necessario approssimarli all’intero più vicino anche se questo può

comportare la comparsa di spazi vuoti nell’immagine di output.

Per evitare ciò si preferisce operare con le trasformazioni inverse: per ogni indice (j,k) dell’array

di output G(j,k) si determina il corrispondente indice nell’immagine di input F(p,q) tramite

11

' ( )

'y

x

p j P J t

q k t

= + − +

= − (1.19)

ove il simbolo di primo indica ancora che p’ e q’ non sono necessariamente interi . E’ allora

opportuno interpolare i valori di F(p,q) in modo da ottenere una stima ragionevole di F’(p’,q’)5 da

assegnare poi a G(j,k) come si vedrà nel prossimo capitolo.

1.4.2 Trasformazioni di scala

La dimensione assoluta di un’immagine può essere cambiata moltiplicando ciascuna coordinata

per un fattore di scala. Applicando lo stesso fattore a tutte le coordinate si ottiene un’immagine

della stessa forma ma di dimensioni differenti, se invece vengono applicati differenti fattori alle

diverse componenti essa varierà sia nelle dimensioni che nella forma. La relazione vettoriale è data

da

0

0k x q

j y p

x s u

y s v

=

(1.20)

dove xs e ys , dette costanti di scalatura, devono essere reali e positive.

Se sono entrambe maggiori di uno si ha un ingrandimento, se sono minori un rimpicciolimento.

La relazione tra gli indici degli array delle immagini di output e input dell’esempio considerato per

la traslazione risulta ora essere

1 1 1'

2 2

1 1 1'

2 2

y

x

p j J Ps

q ks

= + − + +

= − +

(1.21)

Anche in questo caso non è detto che p’ e q’ siano interi ed è perciò necessario ricorrere

all’interpolazione per stimare il valore di F’(p’,q’) da assegnare a G(j,k).

5 L’applicazione F’(p’,q’) è una funzione definita su [0,P]×[0,Q] ⊂ 2

R a valori in I per le immagini monocromatiche e

in r g bI I I× × per quelle a colori tale che F’(p’,q’)=F(p’,q’) se p’,q’ sono interi altrimenti F’(p’,q’) è ottenuto

interpolando F(p,q).

12

1.4.3 Rotazione

La rotazione dell’immagine rispetto all’origine del sistema di assi cartesiano è data da

cos sin

sin cosk q

j p

x u

y v

θ θθ θ

− =

(1.22)

dove θ indica l’angolo di rotazione calcolato rispetto all’asse orizzontale. Si considera positiva la

rotazione antioraria.

La rotazione rispetto ad un generico punto può essere ottenuta traslando l’origine del sistema nel

punto prescelto, applicando la rotazione (1.22) e riportando poi con la traslazione inversa il centro

del sistema nella posizione iniziale.

1.4.4 Composizione di trasformazioni

Le operazioni di traslazione, rotazione e scalatura possono essere composte tra loro. In questi casi

è utile introdurre le coordinate omogenee perchè consentono la concatenazione di un numero

qualsiasi di trasformazioni elementari in un’unica matrice di trasformazione, operazione che non è

possibile compiere utilizzando le coordinate cartesiane perchè le traslazioni richiedono la somma

di matrici.

Si consideri in { }3 \ 0R la relazione di equivalenza data da:

( ) ( ) { }30 1 2 0 1 2 , , , , , \ 0x x x y y y∀ ∈R

i due vettori sono equivalenti

( ) ( )0 1 2 0 1 2 , , , ,x x x y y y∼

se e solo se ∃ λ ∈R tale che i ix yλ= per i= 0,1,2

e il relativo spazio quoziente i cui elementi sono indicati con { }( , , )x y z .

Con le coordinate omogenee si identifica la posizione del punto del piano P di coordinate

cartesiane (x,y) utilizzando uno qualunque dei vettori rappresentanti delle classe { }( , , )hx hy h con

h∈R .

Le coordinate cartesiane sono legate a quelle omogenee dalla relazione

13

hxx

hhy

yh

=

= (1.23)

Per semplicità, in genere, per rappresentare il punto proprio (x,y) si sceglie il vettore (x,y,1) come

rappresentante della classe.

Tutte le trasformazioni finora viste in coordinate omogenee sono ottenute con matrici 3×3 che

hanno il vantaggio di includere i coefficienti di traslazioni, rotazioni e scalature.

L’equazione della traslazione (1.17) viene convertita nel prodotto di matrici

1 0

0 1

1 0 0 1 1

x

y

x t u

y t v

=

(1.24)

quella della rotazione (1.22) in

cos sin 0

sin cos 0

1 0 0 1 1

x u

y v

θ θθ θ

− =

(1.25)

e quella della trasformazione di scala (1.20) in

0 0

0 0

1 0 0 1 1

x

y

x s u

y s v

=

(1.26)

Per combinare un insieme di trasformazioni omogenee individuali in un’unica equivalente

trasformazione si deve prendere il prodotto delle matrici individuali nell’ordine appropriato. Si

indica la generica trasformazione omogenea con la matrice non singolare

11 12 13

21 22 23

31 32 33

a a a

a a a

a a a

=

T (1.27)

Se si desidera applicare n trasformazioni 2 3,, , .... n1T T T T basta concatenare le operazioni in una

matrice come segue

2 1...= nv T T Tu (1.28)

dove v e u sono vettori omogenei.

14

1.4.5 Trasformazioni prospettiche

Lunghezze, angoli, parallelismi e forme spesso, quando vengono osservate o fotografate, appaiono

distorte. Per esempio una rotaia è formata da due binari tra loro paralleli che però appaiono ad un

osservatore come due linee che si avvicinano tra loro con l’aumentare della distanza fino a toccarsi

all’orizzonte.

All’inizio del capitolo si è definita la proiezione prospettica che approssima il modo con cui un

dispositivo fotografico costruisce l’immagine di un oggetto. Approfondiamo in questa sezione gli

aspetti geometrici di tale proiezione.

Un buon modello fisico per lo studio della prospettiva è il sistema formato da un oggetto bi o

tridimensionale (Ob) posto nello spazio tridimensionale, da un piano (π) dell’immagine (I) in cui

questa si forma e da un fuoco (O) posizionato dal lato opposto ad Ob rispetto a π. L’immagine di

ciascun punto x dell’oggetto è il punto x* intersezione tra il raggio congiungente x e il fuoco O,

detto raggio proiettante, e il piano π. Si ottiene così una corrispondenza uno a uno tra i punti

dell’oggetto e i punti dell’immagine.

oggetto O fuoco

immagine

piano dell’immagine

Figura 1. 6 Un modello fisico per la visione in prospettiva di oggetti tridimensionali

15

Il caso in cui O sia tra Ob e π può essere ricondotto al nostro modello poiché l’immagine che si va

a formare su π risulta per simmetria rovesciata ma identica a quella che si otterrebbe su un piano

π* parallelo e opposto, rispetto ad O, a π . 6

π∗ π

h

Ob I* O

h I

Figura 1. 7 L’immagine che si va a formare su ππππ risulta per simmetria rovesciata ma identica a quella che si otterrebbe su un piano ππππ* parallelo e opposto, rispetto ad O, a ππππ .

Questo lavoro si è limitato allo studio della prospettiva di oggetti che si possono assimilare, con

buona approssimazione, ad elementi piani. Le ragioni di questa limitazione vanno ricercate nel

fatto che il nostro scopo è definire trasformazioni dell’immagine che permettano in generale di

ottenere, a partire dalla stessa, un’immagine che simuli un punto di ripresa diverso rispetto a

quello da cui si è scattata la foto. Quando però si fotografa una scena tridimensionale tale

simulazione non è possibile in quanto la trasformazione geometrica che la opera non può essere

definita nello stesso modo per tutti i punti della foto. Un elemento 3D è interpretabile come

composizione di superfici 2D. L’immagine che si ha da un punto di ripresa A di tali superfici

dipende dalla loro posizione rispetto ad A. Quando si cambia punto di ripresa da A a B la nuova

visione che si ha da B di ciascuna di esse dipende dalla loro nuova posizione rispetto a B stesso.

La trasformazione geometrica che simula questo cambiamento deve quindi tenere conto che la

deformazione che l’immagine di ciascuna superficie deve subire è diversa da superficie a

superficie. Nella figura seguente per esempio le deformazioni che le superfici 1, 2 e 3 devono

subire per essere trasformate rispettivamente nella superficie 1’, 2’ e 3’ sono diverse tra loro.

6 Rientra in questo modello il sistema di formazione dell’immagine illustrato nel paragrafo 1.2

16

3 3’

2’

2 1 1’

Figura 1. 8 Il cubo elemento tridimensionale viene decomposto in superfici 2D Quando si vuole simulare un punto di ripresa diverso ciascuna di esse subisce una deformazione diversa per essere trasformata. Esempio la trasformazione che mappa 1 in 1’ preserva le linee verticali, mentre quella che mappa 3 in 3’ no.

Si dovrebbe quindi poter individuare nella foto ciascuna delle superfici che compongono

l’elemento 3D e definire per ognuna di esse una mappa di trasformazione diversa; è impossibile

ottenere un algoritmo, il cui costo computazionale sia accettabile, che realizzi tutto ciò.

Come illustra l’esempio nelle foto 1.9 e 1.10 l’applicazione a immagini di soggetti 3D delle

trasformazioni definite per scene piane che sono è studiate può generare deformazioni anomale

degli oggetti.

Figura 1. 9 Foto originale

17

Figura 1. 10 Correzione della foto 1.9 con proiezione : la facciata del palazzo viene trasformata in modo corretto, la torre, invece, essendo un elemento tridimensionale subisce una deformazione anomala

Infine si può osservare che nella riduzione prospettica dallo spazio 3D dell’oggetto a quello 2D

dell’immagine si perdono informazioni sulla scena: gli elementi della scena che sono “nascosti”

all’obiettivo della macchina non compaiono nell’immagine. In generale però, quando si cambia

punto di ripresa, elementi che prima risultavano “nascosti” potrebbero divenire “visibili”. Nella

figura sottostante si è evidenziato in verde gli elementi non visibili dal punto A e che sono invece

visibili dal punto B. Se si suppone di scattare la foto dal punto A per definire una trasformazione

dell’immagine che ne simuli la ripresa dal punto B è quindi necessario recuperare le informazioni

sugli elementi della scena visibili da B ma è impossibile poter risalire solo attraverso l’immagine

ottenuta da A alle informazioni relative alla zona nascosta rispetto ad A.

18

A

B

Figura 1. 11 Dal punto di ripresa B si possono vedere sia la circonferenza piccola che il lato del pentagono in verde che non sono visibili da A.

Il modello per lo studio della prospettiva di un oggetto piano permette di specificare anche il piano

π0 contenente l’oggetto (piano reale). In particolare i punti appartenenti all’intersezione tra π0 e π

(linea di terra lt) risultano punti fissi della mappa tra Ob ed I.

Piano dell’immagine

fuoco O

immagine

oggetto linea di terra

piano reale

Figura 1. 12 Modello piano per la prospettiva

L’analisi della figura di due rette di π0 , tra loro parallele e incidenti la linea di terra, introduce il

concetto di punto di fuga o punto all’orizzonte definito come il punto dell’immagine in cui le due

19

rette parallele appaiono incontrarsi; nel nostro modello è il limite delle immagini dei punti

dell’oggetto disposti su due linee parallele che si allontanano sempre più dal fuoco.7

punto π

di fuga

O

fuoco

πo lt

Figura 1. 13 Il punto di fuga nell’immagine di due rette parallele ortogonali alla linea di terra è l’intersezione delle loro immagini. L’immagine di rette parallele alla linea di terra preserva tale parallelismo.

Geometricamente il punto all’orizzonte è determinato considerando i piani π1 e π2 passanti per O

e contenenti rispettivamente le rette r1 e r2, tra loro parallele.

punto π

di fuga l

O

π2

πo r2

π1

r1

Figura 1. 14 Individuazione del punto di fuga intersezione tra il piano immagine e la retta l a sua volta intersezione dei piani ππππ1, contenente r1 e O, e ππππ2 contenente r2 e O. 7 Per la definizione e la costruzione del punto all’infinito si è fatto riferimento a[13].

20

Tali piani, poichè hanno O in comune, si intersecano in una linea l contenente O e parallela ad r1 e

r2 e al piano π0. Il punto di fuga è l’intersezione di l con il piano dell’immagine . Infatti l

rappresenta la posizione limite dei raggi proiettanti dei punti di r1 e r2 che si allontanano sempre

più dal fuoco. Ciascun raggio proiettante appartiene a π1 e π2, il loro limite deve essere contenuto

in π1 ∩

π2.

Più in generale per individuare il punto all’orizzonte è sufficiente una sola linea r. Il piano

passante per O e r ed il piano per O, parallelo al piano dell’oggetto, si intersecano in l, che a sua

volta interseca il piano delle immagini nel punto all’orizzonte associato ad r e a tutte le linee di π0

ad essa parallele.

Finora abbiamo considerato linee dell’oggetto perpendicolari al piano dell’immagine e π

ortogonale a π0 : gli stessi risultati possono essere estesi al caso in cui le linee parallele non sono

ortogonali a π , ma comunque lo intersecano, e a piani con angolo di incidenza non retto.

Punto di fuga π

O

s* α

r2

s

r1

πo

Figura 1. 15 Configurazione generale per il punto di fuga: il piano reale e quello immagine non sono ortogonali e le rette r1 e r2 tra loro parallele non sono perpendicolari alla linea di terra

L’immagine di una retta s di π0 parallela alla linea di terra è una retta s* di π anch’essa parallela

alla linea di terra .

Se interpretiamo la prospettiva come una corrispondenza biettiva tra i punti del piano dell’oggetto

e i punti del piano dell’immagine, che associa ad x di π0 x* di π intersezione tra questi e il raggio

21

proiettante tra x e O, si pone il problema che tale applicazione sia definita in ogni punto di π0 e che

ogni punto di π sia immagine di un punto di π0 . Da qui la necessità di introdurre i punti all’infinito

(detti anche punti impropri o ideali) oggetto di analisi della geometria proiettiva e di considerare

coordinate omogenee.

E’importante evidenziare che una prospettiva o una composizione di prospettive hanno la

proprietà di mandare punti linearmente dipendenti in punti linearmente dipendenti (mappa cioè

rette in rette) e che la restrizione della prospettiva a una qualunque retta del piano π0 è ancora una

prospettiva. Da ciò si può concludere che tali trasformazioni sono particolari applicazioni

proiettive8 in cui il fuoco è il centro della trasformazione9.

Se si indica con

( , , )v v v vP x y z↔

il generico punto del piano reale e con

( , , )n n n nP x y z↔

il suo corrispondente nel piano immagine, la generica proiezione è definita dalla relazione

n vP P= T (1.29)

con

11 12 13

21 22 23

31 32 33

a a a

a a a

a a a

=

T (1.30)

matrice non singolare.

1.5 Correzione della distorsione trapezoidale

Questa procedura corregge la distorsione trapezoidale fenomeno di deformazione dell’immagine

che si verifica quando il piano dell’immagine non è parallelo al piano dell’oggetto ma forma con

questo un angolo α rispetto all’asse orizzontale dell’immagine. L’immagine risulta così più larga

nella zona inferiore e più stretta in quella superiore o viceversa (per ulteriori approfondimenti si

rimanda al paragrafo precedente sulla geometria proiettiva ).

8 Per la definizione di proiettività si veda appendice A 9 Per la dimostrazione che una qualunque applicazione biettiva tra due piani che goda delle proprietà appena enunciate sia una proiettività si veda il teorema 12 appendice A

22

Si è visto nel paragrafo 1.2 che dal punto di vista proiettivo il modello della formazione di

un’immagine con una fotocamera può essere schematizzato con una proiezione centrale (figura

1.3) in cui gli oggetti cambiano forma e dimensione in funzione della loro distanza dal centro della

proiezione corrispondente alla lente secondo le equazioni (1.3) e (1.4). Si suppone di lavorare con

foto di scene che possono essere approssimate come piane per cui è possibile definire il piano che

le contiene detto piano reale (pr).

pr

F

αααα

pf

Figura 1. 16 Il piano reale pr e quello dell’immagine pf non sono paralleli ma formano un angolo di incidenza αααα rispetto all’asse orizzontale dell’immagine

23

y pm

pf h

Ly f F

O

-Ly h’ pr

V

αααα

Figura 1. 17 Il piano reale pr e quello dell’immagine pf formano un angolo di intersezione αααα, il piano della mappa pm è parallelo a pr.

L’angolo tra questo e il piano dell’immagine (pf) è detto angolo di inclinazione (α), se è zero

l’immagine che si ottiene è verticale, altrimenti risulta distorta. Si suppone inoltre di essere nei

casi in cui la linea di terra lt, intersezione tra pf e pr è perpendicolare all’asse verticale sia della

scena che dell’immagine per cui le linee parallele ad lt non appaiono alterate.

Questo algoritmo procede alla correzione proiettiva della foto inclinata dal piano della foto ad un

opportuno piano parallelo al piano reale detto piano della mappa (pm). L’immagine ottenuta dopo

24

la correzione prospettica risulta a sua volta una proiezione prospettica rispetto ad F dell’oggetto

sul piano della mappa, poichè però quest’ultimo è scelto parallelo al piano reale, l’immagine

corretta mantiene le proporzioni dell’oggetto reale.

Si considerino sull’immagine di input e su quella di output due sistemi di coordinate cartesiane

ortonormali, rispettivamente v vOx y e ' n nO x y , aventi entrambi origine nel centro della foto, asse

delle ascisse orizzontale orientato da sinistra a destra e asse delle ordinate verticale orientato

dall’alto verso il basso.

(0,0)

x

y

Figura 1. 18 Sistema di coordinate cartesiane associate all’immagine

Poichè le dimensioni delle immagini sono finite si ha che

' '

' '

x v x

y v y

x n x

y n y

L x L

L y L

L x L

L y L

− ≤ ≤− ≤ ≤

− ≤ ≤− ≤ ≤

(1.31)

Con xL e yL dati e 'xL e 'yL da determinarsi.

La proiettività da definire mappa un punto vP dell’immagine di input di coordinate cartesiane

( , )v vx y in un punto nP dell’immagine di output di coordinate ( , )n nx y . Poichè si opera con

proiettività è opportuno passare alle coordinate omogenee: vP è individuato dal vettore ( , , )v v vx y z

legato a ( , )v vx y da

vv

v

vv

v

xx

z

yy

z

=

= (1.32)

25

nP da ( , , )n n nx y z legato a ( , )n nx y da

nn

n

nn

n

xx

z

yy

z

=

= (1.33)

Dato che i punti vP sono propri per semplicità si è posto

1vz =

si ha così

( , ,1)

( , , )v v v

n n n n

P x y

P x y z

↔↔

La proiezione in forma vettoriale assume la forma

n v=P MP (1.34)

1 1 1

2 2 2

3 3 3 1

n v

n v

n

x a b c x

y a b c y

z a b c

=

(1.35)

con , ,i i ia b c ∈R .

Oppure svolgendo il prodotto

1 1 1n v vx a x b y c= + + (1.36)

2 2 2n v vy a x b y c= + + (1.37)

3 3 3n v vz a x b y c= + + (1.38)

Per poter definire univocamente la mappa è necessario conoscere la focale con cui si è scattata la

foto e poter individuare all’interno di questa una linea retta l che nella realtà si sa essere verticale,

in questo modo, come si vedrà più avanti si può risalire all’angolo di inclinazione

26

V

l

y

Figura 1. 19 La retta evidenziata in rosso deve essere raddrizzata

Figura 1. 20 La linea rossa nella foto rappresenta un esempio di retta l che deve essere raddrizzata

Sia f la distanza focale e 1 1 2 2( , ), ( , )v v v vx y x y i punti che individuano l

• Si individua il punto di fuga6 indicato con V di coordinate (0, )v con

6 Per la definizione di punto di fuga si rimanda alla sezione 1.4.5.

27

2 11 1

2 1

v vv v

v v

y yv x y

x x

−= − +−

(1.39)

intersezione tra la retta l e l’asse verticale. Perchè il raddrizzamento possa

essere applicato deve valere

yv L> (1.40)

ossia la foto non deve ritrarre la linea di orizzonte della scena.

• La proiezione preserva la centralità della figura rispetto all’asse verticale y . Tutti i punti

con ascissa 0vx = sono mappati in punti con 0nx = indipendentemente dalla loro

ordinata.

Dall’equazione (1.36) sostituendo si ottiene

1 1 10 0 n v vx a b y c y= + + = ∀

se e solo se

1 0c = (1.41)

1 0b = (1.42)

• Secondo la teoria prospettica i punti con ordinata vy v= appartengono alla linea di

orizzonte e devono essere mappati in un punto improprio (con 0nz = ), indipendentemente

dall’ascissa.

Dall’equazione (1.38) sostituendo si ottiene

3 3 3 0n v vz a x b y c= + + =

vx∀ se e solo se

3 0a = (1.43)

3 3c b v= − (1.44)

• In base alla definizione di punto di fuga fissato il punto vP di coordinate ( , )v vx y la retta t

congiungente vP con V deve essere proiettata in una retta t’ parallela all’asse y. La scelta

della particolare parallela è legata alla scelta del particolare piano pm parallelo al piano

reale ed è quindi arbitraria

28

t’

t

V

•

Pv

y

Figura 1. 21 La retta t deve essere proiettata nella retta t’ parallela all’asse y

y

h pm

Ly F

o

-Ly h’

Figura 1. 22 Il piano pm è scelto in modo tale che h=h’

Si è optato di scegliere il piano pm per cui ogni punto vP del bordo inferiore dell’immagine

di imput di ordinata v yy L= è mappato in punto nP di pm avente distanza h’ dal fuoco F

pari a quella che F ha con vP (indicata con h)

'h h= (1.45)

Dalle similitudini dei triangoli e della simmetria della figura si può osservare che i punti di

coordinate ( , )v yx L vengono così mappati in punti con ascissa n vx x= . Questo equivale a

29

scegliere di mandare la retta t di figura 1.21 nella retta t’ intersecante t proprio nel punto di

ordinata v yy L= appartenente al bordo inferiore della foto.

E’ molto importante sottolineare l’arbitrarietà di questa scelta, una possibile alternativa

sarebbe stata per esempio optare per il piano pm che mantenesse fisse le ascisse dei punti

di ordinata 0vy = .

In conclusione si è scelta la proiezione tale che se v yy L= allora n vx x= .

Sostituendo nelle (1.36), (1.38) tali condizioni si ottiene

1 1

3( )v yn

n vn y

a x b Lxx x

z b L v

+= = =

−

vera se e solo se

1

3

1( )y

a

b L v=

− (1.46)

1 0b = (1.47)

da cui

( )

( )y

n vv

L vx x

y v

−=

− (1.48)

• I bordi orizzontali della foto di input devono essere mandati nei bordi orizzontali di quella

di output in modo tale che quest’ultima risulti centrata. Se v yy L= allora 'n yy L= e se

v yy L= − allora 'n yy L= − indipendentemente dalle ascisse.

Da

2 2 2

3( )n v v

nn v

y a x b y cy

z b y v

+ += =−

sostituendo si ottiene

2 2 2

3

'( )

v yy

y

a x b L cL

b L v

+ +=

− (1.49)

2 2 2

3

'( )

v yy

y

a x b L cL

b L v

− +− =

− − (1.50)

vx∀ da cui

2 0a = (1.51)

30

Sommando la (1.49) e la (1.50) si ottiene

2 2 2 2

3 3

0( ) ( )

y y

y y

b L c b L c

b L v b L v

+ − ++ =

− −

da cui

2

22

yb Lc

v= − (1.52)

allora

( )

2

2

3

yv

nv

Lb y

vy

b y v

−

=−

(1.53)

si ponga 2 3b bb= e si semplifichi

( )

2y

v

nv

Ly

vy b

y v

−

=−

(1.54)

• La deformazione prospettica ha alterato nell’immagine di input la proporzione tra altezza e

larghezza degli elementi della scena. La determinazione del coefficiente b è legata alla

ricostruzione nell’immagine di output delle proporzioni presenti nella scena reale.

Per risalire al rapporto “reale” tra altezze e larghezze degli oggetti si osserva che

nell’intorno di ciascun punto della foto di input, se si considerano grandezze infinitesime,

si può supporre con buona approssimazione che esso non sia stato alterato dalla

deformazione.

Quindi, affinchè nell’immagine di output siano ricostruite le proporzioni della scena reale,

è necessario che la proiezione che si sta definendo preservi il rapporto tra altezza e

larghezza infinitesime presenti nell’intorno del punto dell’immagine di input.

Per rendere meno pesanti i calcoli si consideri il centro dell’immagine di input di

coordinate (0,0) e nel suo intorno si prendano due incrementi infinitesimi ,v vdx dy

rispettivamente lungo l’asse delle ascisse e lungo l’asse delle ordinate. Se si indica con OP

il segmento di estremi (0,0) e (vdx ,0) e con OQ quello di estremi (0,0) e (0,vdy ) la

proiezione, per come è stata finora definita, deve mandare OP in un segmento infinitesimo

31

O’P’ parallelo all’asse delle ascisse dell’immagine di output e OQ in un segmento

infinitesimo O’Q’ appartenente all’asse delle ordinate della stessa; ove O’,P’,Q’ sono le

immagini di O,P,Q. Sia ndx la lunghezza di O’P’ e ndy quella di O’Q’, ndy dipende sia da

vdy che dall’angolo di inclinazione α, a parità di incremento vdy al crescere di α

ndy aumenta, si considera quindi il rapporto

cosn

v

dy

dy

α (1.55)

per esprimere il rapporto tra le due lunghezze a prescindere dall’angolo .

Affinchè la proiezione ricostruisca le proporzioni reali essa deve preservare la proporzione

tra gli incrementi verticale e orizzontale ossia deve essere tale da mantenere invariato il

rapporto

:v vdy dx (1.56)

cioè deve valere

: cos :v v n ndy dx dy dxα= (1.57)

e scambiando gli estremi si ottiene

cosn n

v v

dx dy

dx dy

α= (1.58)

Dalle equazioni (1.48) e (1.54) derivando rispetto a vx e vy si ottiene

yn

v v

v Ldx

dx y v

−=

− (1.59)

( )( )

2 2

2

yn

v v

b L vdy

dy v y v

−=

− (1.60)

che nell’intorno di (0,0) in cui abbiamo considerato gli incrementi valgono

yn

v

v Ldx

dx v

−= − (1.61)

( )2 2

3

yn

v

b L vdy

dy v

−= (1.62)

e sostituendo nella (1.58) si arriva a

32

( )2 2

3cos

y yb L v v L

v vα

− −= − (1.63)

Se si osserva nella figura 1.17 il triangolo OFV è facile vedere che

2 2

cosv

f vα =

+ (1.64)

da cui si ottiene

2 2

y

v f vb

v L

+=

+ (1.65)

Per concludere se si indica con h la proiezione che si è definita si ha

[ ] [ ]

( ) ( ): , , ', ' ', '

, ,

x x y y x x y y

v v n n

h L L L L L L L L

x y x y

− × − → − × −

֏

con

( )

2

2 2

( )

( )y

n vv

yv

ny v

L vx x

y v

Ly

vv f vy

v L y v

−= −

− + =

+ −

(1.66)

Restano da determinare le dimensioni dell’immagine di output ossia le costanti 'xL e 'yL . Poichè

si è scelto di mantenere fisse le ascisse dei punti del bordo inferiore dell’immagine di input si ha

'x xL L= (1.67)

Sostituendo nella (1.49) i valori definiti sopra si ricava

22 2

'( )

y yy

y y

L v Lf vL

v L L v

−+=

+ − (1.68)

Anche in questo caso come per le trasformazioni elementari è opportuno partire dai punti

dell’immagine discreta di output G(j,k) e determinare i punti corrispondenti F(p’,q’)

nell’immagine discreta di input applicando la trasformazione inversa della proiezione h composta

con i cambiamenti di coordinate che trasformano gli indici (j,k) in ( , )n nx y e ( , )v vx y in (p’,q’) e

interpolando nel caso in cui p’ e q’ non siano interi.

33

Nelle foto seguenti si possono osservare i risultati ottenuti implementando le relazioni trovate. In

tutti i casi l’immagine viene raddrizzata mantenendo inalterate le carratteristiche della scena: non

ci sono distorsioni o alterazioni degli elementi in essa contenuti. Si osserva però che parte

dell’immagine è stata tagliata, questa perdita di informazioni è diretta conseguenza del fatto che si

è scelto il bordo inferiore della foto di input come retta che mantiene invariate le ascisse. Si

potrebbe ridurre tale perdita scegliendo per esempio la linea passante per il centro dell’immagine.

Figura 1. 23 Immagine 1.20 corretta con la trasformazione del trapezio

Figura 1. 24 Immagine ottenuta correggendo la foto 1.1

34

1.6 Correzione della distorsione rotazione–trapezio

Questa procedura si interessa al caso in cui si fotografa una scena che può essere approssimata

come piana con una macchina fotografica ruotata lateralmente rispetto alla verticale della scena e

inclinata rispetto al piano della stessa. Se si considera il sistema di coordinate fisso della macchina

fotografica, introdotto nel paragrafo 1.2, il piano reale pr in questo caso risulta ruotato attorno agli

assi x (asse orizzontale) e z (asse ottico) di due angoli non nulli.

y x

z

Figura 1. 25 Il piano reale pr rispetto al sistema di riferimento fisso della macchina fotografica risulta ruotato di un angolo αααα attorno all’asse x e di un angolo ββββ attorno all’asse z. Se ββββ====0000 si rientra nel caso della deformazione del trapezio.

L’immagine che si ottiene risulta affetta dalla deformazione del trapezio e gli elementi che nella

scena reale erano verticali appaiono inclinati rispetto all’asse verticale dell’immagine.

Si indichi con α l’angolo di inclinazione tra il piano reale e quello dell’immagine (angolo di

rotazione attorno a x) e con β l’angolo tra l’asse verticale della scena reale e quello della foto

(angolo di rotazione attorno a z).

35

y’’ pf

y

y’

β pr

F

z

αααα

x

Figura 1. 26 Il piano reale pr e quello dell’immagine pf hanno un angolo di inclinazione αααα. In rosso sono tratteggiati i bordi della foto. L’asse y’ è l’asse verticale della scena reale la sua immagine in pf è l’asse y’’ che forma con l’asse verticale della foto y ( in verde ) un angolo ββββ diverso da zero se la macchina fotografica viene ruotata lateralmente ( cioè attorno all’asse z).

Figura 1. 27 Esempio di foto soggetta a deformazione del trapezio e rotazione rispetto alla verticale. Le rette in blu rappresentano le linee l1 ed l2 selezionate nell’immagine di input, l’asse y, in verde, è l’asse verticale della foto, in rosso sono indicate le linee parallele alla linea di terra. La retta y’’, in giallo, è l’immagine dell’asse verticale (y’) della scena reale ortogonale alle linee rosse. L’angolo tra y e y’’ è l’angolo di rotazione laterale ββββ, la loro intersezione il centro di tale rotazione.

36

Il piano reale pr e quello dell’immagine pf si intersecano lungo una linea di terra che risulta

ortogonale all’asse verticale della scena reale e forma un angolo pari a π+β con l’asse verticale

della foto (y) per cui nell’immagine le trasformate delle linee parallele alla linea di terra, che

preservano il parallelismo, costituiscono un fascio di rette parallele che forma un angolo π+β con

y ed è perpendicolare alla retta y’’ immagine dell’asse verticale della scena reale y’ . L’angolo tra

y’ e y è l’angolo β e il loro punto di intersezione è il centro C della rotazione laterale, in genere

distinto dal centro della foto O di coordinate (0,0).

La correzione prospettica che si deve applicare sull’immagine di input deve essere tale da

raddrizzare gli elementi inclinati, rendere orizzontali le linee parallele alla linea di terra e

correggere la deformazione del trapezio.

Si è scelto di operare questa trasformazione in due passi: uno che porta alla definizione di una

proiezione t che va a correggere la deformazione trapezoidale, l’altro che raddrizza la scena

applicando una opportuna rotazione r ai suoi elementi. Se si indica con I l’immagine di input e con

O quella di output, la trasformazione p che mappa I in O può essere vista come la composizione

delle due trasformazioni t e r

( ) ( ( ))O p I r t I= =

Il problema si riduce quindi alla definizione di una proiezione t del tipo definito nella sezione 1.5 e

di una rotazione r. Affinchè esse siano definite in modo univoco si deve essere in grado di risalire

agli angoli α e β.

Si considerino sull’immagine di input I e su quella di output O due sistemi di coordinate cartesiane

ortonormali, rispettivamente v vOx y e ' n nO x y , aventi entrambi origine nel centro della foto, asse

delle ascisse orizzontale orientato da sinistra a destra e asse delle ordinate verticale orientato

dall’alto verso il basso.

37

(0,0)

x

y

Figura 1. 28 Sistema di coordinate cartesiane associate all’immagine

Poichè le dimensioni delle immagini sono finite si ha che

' '

' '

x v x

y v y

x n x

y n y

L x L

L y L

L x L

L y L

− ≤ ≤− ≤ ≤

− ≤ ≤− ≤ ≤

Con xL e yL dati e 'xL e 'yL da determinarsi.

La trasformazione p da definire mappa un punto vP dell’immagine di input di coordinate

cartesiane ( , )v vx y in un punto nP dell’immagine di output di coordinate ( , )n nx y

( ) ( ( ))n v vP p P r t P= =

Per definire in modo univoco le mappe t e r è necessario :

o conoscere la focale con cui è stata scattata la foto,

o poter individuare all’interno di questa due rette che nella realtà si sa essere parallele tra

loro e verticali, ma che nell’immagine appaiono inclinate ed incidenti nel punto di fuga V,

o determinare la posizione del centro C della rotazione laterale.

La focale f e le rette l1 e l2 non sono sufficienti per definire in modo univoco le trasformazioni

infatti, come si può osservare nell’immagine sottostante, le rette l1 e l2 individuano il punto di

fuga V, ma è necessario conoscere anche il centro C della rotazione laterale per determinare

l’angolo β

38

l1

l2 V

β1 C1 β2 C2

y

Figura 1. 29 Le rette l1 e l2 individuano il punto di fuga coincidente con la loro intersezione ma è necessario determinare la posizione del centro di rotazione laterale C ( punto di intersezione tra y e y’’) per determinare ββββ. Nell’immagine è evidenziato come le stesse rette l1 e l2 possano corrispondere a rotazioni laterali di centro e angoli diversi (C1 ≠≠≠≠ C2 e ββββ1 ≠≠≠≠ ββββ).

Noto C infatti si può individuare la retta y’’ passante per C e per il punto V e determinare l’angolo

tra y e y” che coincide con β. La proiettività p deve essere tale da mappare y’’ nell’asse verticale

della foto di output ( 0nx = ) e ogni retta per V in una retta parallela a 0nx = .

Poichè si considerano immagini di cui si conosce solo la focale e le rette l1 e l2 i dati a nostra

disposizione sono insufficienti per individuare C all’interno della foto di input, si è scelto perciò di

far procedere l’utente per tentativi facendogli fornire un valore λ proporzionale alla distanza tra O

e C.

Dati la focale f, i punti (x1,y1),(x2,y2) che individuano l1, (x3,y3), (x4,y4) che individuano l2 e λ

si procede alla definizione di t e r:

• Si determina il punto di fuga V di coordinate (u,v) intersezione tra l1 ed l2 con

( 1 3)( 2 1)( 4 3) 1( 2 1)( 4 3) 3( 4 3)( 2 1)

( 4 3)( 2 1) ( 2 1)( 4 3)

y y x x x x x y y x x x y y x xu

y y x x y y x x

− − − − − − + − −=− − − − −

(1.69)

39

( 3 1)( 4 3)( 2 1) 3( 2 1)( 4 3) 1( 4 3)( 2 1)

( 4 3)( 2 1) ( 2 1)( 4 3)

x x y y y y y y y x x y y y x xv

y y x x y y x x

− − − − − − + − −=− − − − −

(1.70)

• Si calcola la distanza CV secondo la formula

2 2(1 ) (1 )CV OV u vλ λ= + = + + (1.71)

Condizione per poter applicare t è che

x

y

u L

v L

> >

• La proiezione t che corregge la deformazione trapezoidale mappa le rette passanti per il

punto di fuga V in rette parallele alla retta per C e V. La sua equazione è del tipo (1.66)

definita nel paragrafo 1.5 per cui

[ ] [ ]

( ) ( ): , , ', ' ', '

, ', '

x x y y x x y y

v v

t L L L L L L L L

x y x y

− × − → − × −

֏

con

( )

2

2 2

( )'

( )

'

yv

v

yv

y v

L vx x

y v

Ly

vv f vy

v L y v

−= −

− + =

+ −

con

v CV= (1.72)

• Si determinano le dimensioni dell’immagine di output ossia le costanti 'xL e 'yL . Per la

(1.67) si ha

'x xL L=

e per la (1.68)

22 2

'( )

y yy

y y

L v Lf vL

v L L v

−+=

+ −

40

• Si definisce la rotazione r che muta il fascio di rette parallele alla retta per C e V e

contenente le trasformate tramite t di l1 ed l2 (indicate con t(l1)e t(l2)) nel fascio di rette

parallele all’asse verticale della foto di output. Data l’incertezza della posizione di C si

preferisce centrare tale rotazione in O . Per determinare l’angolo di rotazione γ che

trasformi t(l1) e t(l2) in rette verticali si calcolano gli angoli che t(l1) e t(l2) formano con la

verticale della foto, indicati rispettivamente con β1 e β2: la tangente di γ è uguale alla

media delle tangenti di β1 e β2. Siano (t(x1),t(y1)), (t(x2),t(y2)), (t(x3),t(y3)) e (t(x4),t(y4))

le immagini dei punti dati ottenute applicando t, β1 e β2 sono dati da

( 2) ( 1)

1( 2) ( 1)

t x t xtg

t y t yβ −=

− (1.73)

( 4) ( 3)

2( 4) ( 3)

t x t xtg

t y t yβ −=

− (1.74)

da cui

1 2

2

tg tgtg

β βγ += (1.75)

Posti

2

1cos

1 tgγ

γ=

+

2

sin1

tg

tg

γγγ

=+

r è così definita

[ ] [ ]: ', ' ', ' ', ' ', '

( ', ') (, )

x x y y x x y y

n n

r L L L L L L L L

x y x y

− × − → − × −

֏

con

cos ' sin '

sin ' cos 'n

n

x x y

y x y

γ γγ γ

= − = +

(1.76)

In conclusione p ha equazione

[ ] [ ]: , , ', ' ', '

( , ) ( , )

x x y y x x y y

v v n n

p L L L L L L L L

x y x y

− × − → − × −

֏

con

( ( ))

( ( ))n v

n v

x r t x

y r t y

= =

(1.77)

41

Anche in questo caso, come per la trasformazione del trapezio, è opportuno partire dai punti

dell’immagine discreta di output G(j,k) e determinare i punti corrispondenti F(p’,q’)

nell’immagine discreta di input applicando la trasformazione inversa della proiezione p, cioè

1 1 1p t r− − −= � 10, composta con i cambiamenti di coordinate che trasformano gli indici (j,k) in

( , )n nx y e ( , )v vx y in (p’,q’) e interpolando nel caso in cui p’ e q’ non siano interi.

Di seguito sono riportati esempi di foto corrette applicando la trasformazione (1.77)

Figura 1. 30 Correzione della foto 1.27 Con λλλλ =10.5%

10 Nella composizione di mappe si assume che l’ordine con cui si applicano va da destra a sinistra

42

Figura 1. 31 Immagine originale

Figura 1. 32 Correzione della foto 1.31 con λλλλ=7%

Si può osservare che nelle foto corrette compaiono zone nere, esse corrispondono a punti

dell’immagine di output che la trasformazione rotazione associa a punti che sono fuori

dall’immagine di input.

43

1.7 Correzione di tipo proiezione

Questa procedura viene applicata su immagini in cui si possono identificare due punti di fuga V1 e

V2 . Supponendo di lavorare con foto di scene che possono essere approssimate come piane per

cui è possibile definire il piano che le contiene detto piano reale pr; questo tipo di immagini si

ottiene quando il punto di ripresa è tale per cui pr risulta ruotato attorno agli assi x,y,z del sistema

di coordinate fisso della macchina fotografica di tre angoli non nulli.

y x

z

Figura 1. 33 Il piano reale pr rispetto al sistema di riferimento fisso della macchina fotografica risulta ruotato di un angolo αααα attorno all’asse x, di un angolo ββββ attorno all’asse y e di un angolo γγγγ attorno all’asse z. Se ββββ====γγγγ====0000 si rientra nel caso della deformazione del trapezio, se ββββ====0000 si è nel caso della deformazione rotazione –trapezio.

Anche in questo caso la formazione dell’immagine può essere schematizzata con una proiezione

centrale con centro nella lente della macchina: le rette verticali e orizzontali della scena reale

vengono mappate dalla proiezione rispettivamente in un fascio di rette passante per V1 e in uno

passante per V2. La trasformazione prospettica p che si definisce nella procedura deve essere tale

da deformare l’immagine in modo che sembri ottenuta da un punto di ripresa frontale alla scena, p

è quindi un’applicazione da pf a un piano pm parallelo al piano reale e trasforma tutte le rette per

V1 in rette verticali e quelle passanti per V2 in rette orizzontali.

Indicata con I l’immagine di input e con O quella di output si ha

( )O p I= (1.78)

44

Figura 1. 34 Esempio di immagine a cui applicare la correzione di tipo proiezione

La proiezione p è stata decomposta nel prodotto di due proiezioni h e g del tipo definito nel

paragrafo 1.5 :

( ( ))O g h I= (1.79)

h mappa le rette passanti per il punto di fuga V1 in rette verticali, g trasforma le rette per V2 in

rette orizzontali mantenendo invariate quelle verticali.

V1

V2 V2

Figura 1. 35 La correzione proiezione può essere eseguita applicando due trasformazioni tipo trapezio una che mappa le rette passanti per V1 in rette verticali, la seconda che manda le linee per V2 in rette orizzontali mantenendo invariate quelle verticali.

45

Per definire in modo univoco h e g è necessario conoscere la focale con cui si è scattata la foto e

poter individuare all’interno di questa una coppia di rette (l1 el2) che nella realtà si sa essere

parallele e verticali e che nell’immagine appaiono oblique e incidenti in quello che è il punto di

fuga V1 e una coppia di rette (l3 e l4) che nella realtà si sa essere parallele e orizzontali ma

appaiono nella foto oblique e incidenti nel secondo punto di fuga V2.

Si considerino sull’immagine di input e su quella di output due sistemi di coordinate cartesiane

ortonormali, rispettivamente v vOx y e ' n nO x y , aventi entrambi origine nel centro della foto, asse

delle ascisse orizzontale, orientato da sinistra a destra, e asse delle ordinate verticale, orientato

dall’alto verso il basso, come illustrato in figura 1.18.

Poichè le dimensioni delle immagini sono finite si ha che

" "

" "

x v x

y v y

x n x

y n y

L x L

L y L

L x L

L y L

− ≤ ≤− ≤ ≤

− ≤ ≤− ≤ ≤

Con xL e yL dati e "xL e " yL da determinarsi.

La proiettività da definire mappa un punto vP dell’immagine di input di coordinate cartesiane

( , )v vx y in un punto nP dell’immagine di output di coordinate ( , )n nx y .

Data la focale f e i punti (x1,y1), (x2,y2) che individuano l1,(x3,y3),(x4,y4) che danno l2,

(x5,y5),(x6,y6) per l3 e (x7,y7),(x8,y8) per l4

• Si individua il punto di fuga V1 intersezione di l1 e l2 di coordinate (u1,v1) con

( 1 3)( 2 1)( 4 3) 1( 2 1)( 4 3) 3( 4 3)( 2 1)

1( 4 3)( 2 1) ( 2 1)( 4 3)

y y x x x x x y y x x x y y x xu

y y x x y y x x

− − − − − − + − −=− − − − −

(1.80)

( 3 1)( 4 3)( 2 1) 3( 2 1)( 4 3) 1( 4 3)( 2 1)

1( 4 3)( 2 1) ( 2 1)( 4 3)

x x y y y y y y y x x y y y x xv

y y x x y y x x

− − − − − − + − −=− − − − −

(1.81)

Perchè il raddrizzamento h possa essere applicato deve valere

1

1x

y

u L

v L

> >

• Si considera la traslazione t1 che porta l’asse verticale 1vx u= in 0vx = di equazione

46

1v v

v v

x x u

y y

= − =

(1.82)

il punto di fuga V1 ha ora coordinate (0,v1)

• Si applica la proiezione h definita nella sezione 1.5 di equazione

[ ] [ ]

( ) ( ): , , ', ' ', '

, ', '

x x y y x x y y

v v

h L L L L L L L L

x y x y

− × − → − × −

֏

con

( )

2

2 2

( )'

( )

'

yv

v

yv

y v

L vx x

y v

Ly

vv f vy

v L y v

−= −

− + =

+ −

(1.83)

e con v=v1.

Le dimensioni dell’immagine dopo l’applicazione di h ossia le costanti 'xL e 'yL sono date

da

'x xL L=

secondo la (1.67) (si è scelto di mantenere fisse le ascisse dei punti del bordo inferiore

dell’immagine di input ) e da

22 2

'( )

y yy

y y

L v Lf vL

v L L v

−+=

+ −

secondo la (1.68)

• Si considerano le immagini dei punti (xi,yi) con i=5,6,7,8 ottenute applicando h e t1 cioè

(h(t1(xi)),h(t1(yi))) con i =5,6,7,8 . Per comodità si pone per i=5,6,7,8

' ( 1( ))

' ( 1( ))

xi h t xi

yi h t yi

==

Il punto di fuga V2 è l’intersezione tra l’immagine di l3 e quella di l4 . Le sue coordinate

sono (u2,y2) con

47

( 5 ' 7 ')( 6 ' 5 ')( 8' 7 ') 5 '( 6 ' 5')( 8 ' 7 ') 7 '( 8 ' 7 ')( 6 ' 5 ')

2( 8' 7 ')( 6 ' 5 ') ( 6 ' 5')( 8' 7 ')

y y x x x x x y y x x x y y x xu

y y x x y y x x

− − − − − − + − −=− − − − −

(1.84)

( 7 ' 5')( 8' 7 ')( 6 ' 5 ') 7 '( 6 ' 5')( 8' 7 ') 5'( 8' 7 ')( 6 ' 5')

2( 8' 7 ')( 6 ' 5') ( 6 ' 5 ')( 8' 7 ')

x x y y y y y y y x x y y y x xv

y y x x y y x x

− − − − − − + − −=− − − − −

(1.85)

Perchè il raddrizzamento g possa essere applicato deve valere

2 '

2 'x

y

u L

v L

> >

• Si considera la traslazione t2 che porta l’asse verticale ' 2y v= in ' 0y = di equazione

' '

' ' 2

x x

y y v

= = −

(1.86)

il punto di fuga V2 ha ora coordinate (u2,0).

• Si definisce la proiezione g che mappa il fascio di rette passanti per V2 nel fascio delle

rette orizzontali. L’applicazione g è ottenuta dalla applicazione definita nella sezione 1.5

invertendo il ruolo delle ascisse e delle ordinate, si ha cosi’

[ ] [ ]( ) ( )

: ', ' ', ' ", " ", "

', ' ,

x x y y x x y y

n n

g L L L L L L L L

x y x y

− × − → − × −

֏

con

( )

2

2 2

( ')'

'' ( ')

' ' ' '

( ' ')'

( ' ')

x

nx

xn

Lx

vv f vx

v L x v

L vy y

x v

− + = + −

− = −

(1.87)

e con ' 2v u= .

Le dimensioni dell’immagine di output ottenuta applicando la g sono

" 'y yL L= (1.88)

poichè g tiene fisse le ordinate dei punti del lato verticale destro dell’immagine

e

48

2 2 2( ') ' ' ( ')

"' ' ( ' ')

x xx

x x

f v L v LL

v L L v

+ −=+ −

(1.89)

ottenuta dalla (1.68) sostituendo l’ordinata con l’ascissa

In conclusione p ha equazione

[ ] [ ]

( ) ( ): , , ", " ", "

, ,

x x y y x x y y

v v n n

p L L L L L L L L

x y x y

− × − → − × −

֏

con

( 2( ( 1( ))))

( 2( ( 1( ))))n v

n v

x g t h t x

y g t h t y

= =

(1.90)

Anche in questo caso, come per le trasformazioni precedenti, è opportuno partire dai punti

dell’immagine discreta di output G(j,k) e determinare i punti corrispondenti F(p’,q’)

nell’immagine discreta di input applicando la trasformazione inversa della proiezione composta

con i cambiamenti di coordinate che trasformano gli indici (j,k) in ( , )n nx y e ( , )v vx y in (p’,q’) e

interpolando nel caso in cui p’ e q’ non siano interi.

L’immagine seguente è un esempio di foto corretta con questo algoritmo.

.

Figura 1. 36 Foto 1.34 corretta con proiezione

49

PARTE II : interpolazione dei livelli di grigio

2.1 Introduzione

Le tecniche di interpolazione sono di fondamentale importanza per la ricostruzione delle immagini

digitalizzate e quindi per la realizzazione delle trasformazioni geometriche descritte nella prima

parte. Queste, infatti, per essere operate necessitano, oltre che di una trasformazione spaziale,

anche di una interpolazione dei livelli delle intensità del colore. Indicata con t la proiezione che

trasforma l’immagine di input F nell’immagine di output G il processo di ricostruzione in genere

“opera a ritroso” partendo dai punti di G e applicando la trasformazione inversa t-1 per determinare

il valore corrispondente in F da assegnargli

1( , ) ( ( , ))G j k F t j k−=

Poichè di solito t-1 mappa la generica coppia di indici (j,k) di G in una coppia di coordinate

1( , ) ( ', ')t j k p q− = che non è una coppia di indici di F risulta necessario associare a questa un

opportuno valore 1( ', ') ( ( , ))F p q F t j k−= ottenuto interpolando i livelli delle intensità del colore di

F in un intorno di (p’,q’).

Il colore di ciascun punto è definito attraverso una terna di valori (r,g,b) punto dello spazio , ,R G BS .

Si è scelto di interpolare i livelli delle intensità dei colori fondamentali rosso, verde e blu in modo

indipendente, si parla così di interpolazione dei livelli di grigio riferendosi alle intensità di

ciascun colore. Naturalmente il valore approssimato ottenuto per interpolazione deve essere

anch’esso un punto di , ,R G BS . Se questa condizione non è verificata lo si sostituisce con il punto di

, ,R G BS le cui coordinate sono gli interi più vicini alle intensità calcolate. L’arrotondamento

dell’intensità di ciascun colore causa un errore che viene diffuso scaricando l’errore stesso sui

punti dell’immagine adiacenti al punto considerato. In questo caso si è scelto di scaricare l’errore

sul pixel successivo al punto considerato. Al fine di dare una formulazione più rigorosa del

problema dell’interpolazione dei livelli di grigio e di analizzare i metodi usati per risolverlo si

danno alcune nozioni generali sul problema dell’interpolazione.

50

2.2 Interpolazione e migliore approssimazione

Il problema dell’approssimazione di funzioni, di importanza fondamentale nella matematica

applicata, consiste nella “sostituzione” di una funzione f con una funzione scelta nell’ambito di

una fissata classe di funzioni in modo tale che il comportamento delle due soddisfi specificate

condizioni imposte dal problema, intuitivamente si dice che il loro andamento deve essere il più

simile possibile.

Una tipica situazione in cui si presenta un’operazione di questo tipo è la seguente: la funzione f da

approssimare non è nota, ma di essa si conoscono alcuni valori su un insieme di punti e si vogliono

avere indicazioni sul suo andamento in altri punti . In base al tipo e alla quantità dei dati e alle

caratteristiche del problema si può scegliere tra due differenti approcci: il primo cerca di costruire

la funzione che passi per certi punti assegnati (interpolazione), il secondo cerca di costruirne una

che si scosti di poco dai dati (migliore approssimazione).

In generale un problema interpolazione è posto nei seguenti termini: dati n+1 punti ix distinti con

i=0,1,...n si vuole determinare una funzione ϕ , appartenente ad una certa famiglia di funzioni Φ,

che in essi soddisfi m condizioni di interpolazione (vincoli che la ϕ e le sue derivate devono

soddisfare nei punti ix ) .

Nel caso in cui le condizioni siano del tipo

( ) 0,1....i ix y i nϕ = =

e la famiglia di funzioni Φ sia definita da

0 1( , , ,... )nx a a aϕΦ =

si cercano i valori 0 1, ,... na a a tali che

0 1 ( , , ,... ) 0,1,...i n ix a a a y i nϕ = =

Se ϕ dipende linearmente dai parametri ossia

0 1 0 0 1 1( , , ,... ) ( ) ( ) .... ( )n n nx a a a a x a x a xϕ ϕ ϕ ϕ= + + +

il problema è detto di interpolazione lineare.

Nel caso in cui la famiglia delle funzioni Φ sia uno spazio lineare su R di dimensione n+1, che

per comodità indichiamo con V, e se con V* si denota il corrispondente spazio duale, la

51

definizione può essere così generalizzata: assegnati n+1 funzionali lineari 0 1, ,.. nL L L in V* e dato

l’insieme di numeri reali { } 0,1..i i ny

= si cerca l’elemento v di V tale che

i=0,1,...ni iL v y= (2.1)

Se nV P= e ( )i iL v v x= ove { } 0,1..i i nx

=è un insieme di n+1 punti distinti si parla di interpolazione

polinomiale.

Scelta una base di V { }0 1, ,.... nϕ ϕ ϕ lo stesso problema può essere scritto in forma matriciale infatti

ogni elemento di V può essere scritto come combinazione lineare degli elementi della base

0

n

i ii

aϕ=∑

e la definizione di v che soddisfa le condizioni (2.1) equivale alla determinazione dei parametri ia

che soddisfano le

0

i=0,1,...nn

i j j ij

L a yϕ=

=

∑ (2.2)

da cui

=Ga p (2.3)

con a vettore dei parametri, p vettore dei valori iy e

0 0 0

0

( ) ( )

( ) ( )

n

n n n

L L

L L

ϕ ϕ

ϕ ϕ

=

G

…

⋮ ⋱ ⋮

⋯

(2.4)

detta matrice di Gram.

Il problema della migliore approssimazione può essere così enunciato: sia V uno spazio lineare

dotato di una norma • e W un suo sottospazio generato dagli elementi 1 2, ,.... nϕ ϕ ϕ di V. Il

sottospazio W ha come elementi le combinazioni lineari

1

n

i ii

aϕ=∑

e ha dimensione finita. Assegnato un elemento f di V si tratta di determinare un elemento w* di

W, da considerare migliore approssimazione di f nella metrica di V, in modo tale che la distanza

da f sia minima cioè

* f w f w w W− ≤ − ∀ ∈

52

In altre parole si tratta di determinare il minimo della funzione

0 11 0

( , ,.... )n

n i id a a a f a ϕ=

= −∑

se [ ]0 1* , * ,... * na a a minimizza d si dice che w* è l’elemento di migliore approssimazione di f.

Centrale nella risoluzione di questo tipo di problemi è la scelta della norma di V e del tipo di

funzione approssimante, cioè del sottospazio W in cui ricercare il w*, più adatti. Di solito si sceglie

la funzione approssimante tra i polinomi per la loro efficienza computazionale. Nel caso di

2( , )V L a b≡ con (a,b)⊂R si parla di migliore approssimazione nel senso dei minimi quadrati e si

ricerca il polinomio np P∈ che minimizza la norma del vettore residuo r di componenti

( ) ( )i if x p x−

1

22

220

( , ) ( ( ) ( ))m

i ii

r d f p f x p x=

= = − ∑ (2.5)

Se il problema è discreto, ossia si lavora su un insieme finito { }1 2, ,....... mS x x x= dell’intervallo

(a,b) e per ogni funzione f definita su S si considera la seminorma

1

22

2,1

( ( ))m

S ii

f f x=

= ∑�

la migliore approssimazione si ottiene cercando il polinomio np P∈ tale che

2, 2,* S S nf p f p p P− ≤ − ∀ ∈� �

2.2.1 Alcuni risultati sull’unisolvenza dei problem i di interpolazione in nR

Sia V uno spazio lineare di funzioni reali di dimensione finita n+1 e V* il suo duale. Dati

0 1, ,.. nL L L appartenenti a V* e l’insieme { } 0,1..i i ny

=di reali in base alla definizione date nel

paragrafo precedente risolvere il problema di interpolazione in V significa trovare la funzione

v∈V tale che i iL v y= per i=0,1,..n.

Si dice che 0 1, ,.. nL L L sono unisolventi in V se esiste una ed una sola soluzione al problema.

0 1, ,.. nL L L sono unisolventi se e solo se sono indipendenti in V* infatti vale

53

Lemma 2.1

Sia V spazio lineare di dimensione n . Se 1,.. nv v sono indipendenti in V e 1,.. nL L sono indipendenti

in V* allora

( ) 0i jL v ≠ (2.6)

Viceversa se 1,.. nv v o 1,.. nL L sono indipendenti e vale la (2.6) allora anche l’altro insieme è

indipendente.

Dimostrazione. Si supponga che ( ) 0i jL v = allora anche ( ) 0j iL v = .

Il sistema

1 1 1 2 2 1 1

1 1 2 2

( ) ( ) ...... ( ) 0

( ) ( ) ...... ( ) 0

n n

n n n n n

a L v a L v a L v

a L v a L v a L v

+ + + =⋅⋅

+ + + =

avrebbe una soluzione non banale 1,... na a . Allora per i =1,2,...n

1 1 2 2( ...... )( ) 0n n ia L a L a L v+ + + =

e poichè 1,.. nv v è una base di V

1 1 2 2( ...... )( ) 0n na L a L a L v+ + + =

per ogni v∈V, per cui 1 1 2 2 ...... 0n na L a L a L+ + + = e quindi 1,.. nL L sono linearmente dipendenti

contro le ipotesi.

Teorema 2.2

Sia V uno spazio lineare di dimensione n e 1,.. nL L n elementi di V*. Il problema di interpolazione

(2.1) ammette una e una sola soluzione per n valori arbitrari 1 2, ,... ny y y se e solo se gli iL sono

indipendenti in V*.

Dimostrazione. Per il lemma sopra se 1,.. nv v e 1,.. nL L sono indipendenti ( ) 0i jL v ≠ e il sistema

1 1 2 2( ...... ) i=1,2,..ni n n iL a v a v a v y+ + + =

o

54

1 1 2 2( ) ( ) ...... ( )i i n i n ia L v a L v a L v y+ + + = (2.7)

ha una soluzione 1,... na a e l’elemento

1

n

i ii

a v=∑ (2.8)

risolve l’interpolazione lineare .

Viceversa se il problema di interpolazione ha soluzione per ogni 1 2, ,... ny y y la (2.7) ha soluzione

per ogni 1 2, ,... ny y y e quindi ( ) 0i jL v ≠ . Allora per il lemma gli iL sono indipendenti.

Il determinante ( )i jL v è detto determinante di Gram generalizzato.

Si considerino ora i casi in cui il problema di interpolazione è polinomiale, cioè ( )i iL v v x= con v

ricercato tra i polinomi.

In una dimensione l’unisolvenza è garantita dal seguente risultato

Teorema 2.3

Si indichi con nP l’insieme dei polinomi p(x) con x ∈R a coefficienti reali di grado inferiore o

uguale ad n dati dall’espressione

10

( ) n

ii

i

p x a x a=

= ∈∑ R (2.9)

Dati n+1 punti di collocazione ( , )i ix y con i=0,1,2,..n e i jx x≠ se i≠j si può sempre determinare in

nP il polinomio p tale che

( ) 0,1,....i ip x y i n= = (2.10)

La dimostrazione dell’esistenza e unicità di tale polinomio si ottiene costruttivamente

introducendo in nP la base dei polinomi interpolanti di Lagrange che per ogni i soddisfano la

condizione

,

1 se i=k( )

0 se i ki k i kL x δ = = ≠

(2.11)

∀k=0,1,..n

la cui rappresentazione esplicita è

55

0,

( ) i=0,1,...nn

ji

j j i i j

x xL x

x x= ≠

−=

−∏ (2.12)

e rappresentando p così

( )0

( )n

i i xi

p x y L=

=∑ (2.13)

p è detto polinomio di Lagrange.

Questo teorema non può essere esteso al caso di interpolazioni polinomiali in più variabili in

quanto è stato dimostrato che in ogni insieme aperto di nR con n≥2 non è sempre possibile, dati n

punti arbitrari nix ∈R distinti, trovare una combinazione lineare di 0 1,... np p − (base dello spazio

dei polinomi di ordine minore o uguale ad n su nR ) che assumono assegnati valori su questi punti.

Teorema 2.4 (di Haar)

Sia S un insieme in nR con n≥2 che contenga un punto interno p e siano 1,... nf f (n>1) funzioni

definite su S e continue in un intorno di p .Questo insieme di funzioni non può essere unisolvente

su S.

Dimostrazione. Sia B una bolla centrata in p sufficientemente piccola affinchè le1,... nf f siano ivi

continue. Si scelgano n punti distinti in B 1 2, ,... np p p e si assuma che ( ) 0i jf p ≠ . Si fissino

3,... np p e si muovano 1p e 2p in modo continuo in b in maniera che si scambino di posizione tra

loro. In questo modo si è introdotto uno scambio di colonne in ( )i jf p il cui segno di conseguenza

cambia. Poichè le funzioni sono continue deve esserci qualche posizione intermedia di 1p e 2p

per cui il valore del determinante si annulla.

Da quanto detto si può quindi concludere che in generale per l’interpolazione di funzioni in più

variabili dare un numero di condizioni pari ai parametri del polinomio non garantisce la sua

costruibilità. In questi casi è opportuno risolvere il problema trovando il polinomio di migliore

approssimazione.

56

2.2.2 Interpolazione polinomiale in una e due dimens ioni in forma matriciale

Molto utile è la rappresentazione in forma matriciale del problema di interpolazione polinomiale in

una dimensione. Scelta come base per ( )nP x { }21, , ,.. nx x x si ha che il generico polinomio 0

ni

ii

a x=∑

deve soddisfare

0

( )n

ii j j

i

a x y=

=∑ (2.14)

con j=0,1,..n e { }0,... nx x n+1 punti distinti a cui corrispondono i valori { }0,... ny y da cui

=Va p (2.15)

con a vettore dei parametri, p vettore degli jy e

0 01

1

n

nn n

x x

x x

=

V

⋯

⋮ ⋱ ⋮

⋯

(2.16)

detta matrice di Vandermonde che per il teorema 2.2 ha determinante diverso da zero se e solo se i

punti sono distinti.

La stessa definizione può essere estesa ai problemi in due dimensioni: si scelga come base di

2 ( , )n

P x y l’insieme { }, 0,..

i j

i j nx y

=, le condizioni di interpolazione diventano

, 0

ni j

k ij k ki j

z a x y=

= ∑ (2.17)

con k=0,1,..m e { }( , )k kx y punti tra loro distinti . Indicato con a il vettore dei parametri incogniti

lungo s e con p quello dei valori assegnati si ottiene

Ma = p (2.18)

con

1 11

1

n n

n nm m

x y

x y

=

M

…

⋮ ⋱ ⋮

⋯

(2.19)

57

In base a quanto detto prima anche se i punti sono distinti e m=s M potrebbe avere determinante

nullo e il problema non ha soluzione. Se m>s il sistema è sovradeterminato e conviene cercare il

polinomio di migliore approssimazione che minimizzi lo scarto quadratico medio

2

0

( ( , ) )m

k k kk

p x y z=

−∑ (2.20)

Indicato con r il vettore residuo

r = p - Ma (2.21)

scrivendo la condizione necessaria di minimo (porre nulle le derivate di r rispetto ija ) si ottiene

T TM Ma = M p (2.22)

detto sistema delle equazioni normali da cui

-1T Ta = M M M p (2.23)

in cui la matrice -1T TM M M è detta pseudoinversa di M.

2.2.3 Interpolazione integrale

Come si è visto la scelta di punti distinti è condizione necessaria ma non sufficiente per

l’unisolvenza dell’interpolazione polinomiale in due dimensioni. Si può dimostrare che tale

condizione risulta sufficiente nel caso dell’interpolazione integrale in una e due dimensioni. Per

questo tipo di interpolazione in una dimensione si va a determinare il polinomio p(x) che soddisfa

le condizioni

( )i

i

x h

i

x h

p x dx z+

−

=∫ (2.24)

con 1=0,1,…n .

Questa interpolazione porta al medesimo risultato dell’interpolazione polinomiale in R .

Teorema 2.5

In una dimensione l’interpolazione integrale di un polinomio di grado n è equivalente a quella

polinomiale.

58

Dimostrazione. Si considerino n+1 punti distinti 0 1. ,...., nx x x a cui sono associati i valori

0 1. ,...., ny y y e si determini il polinomio interpolante di ordine n in una variabile

0 1( ) ... nnp x a a x a x= + + +

tale che ( )i ip x y= per i =0,1,…n, in base a quanto affermato prima tale polinomio esiste ed è

unico. Si consideri poi il polinomio p’(x ) di ordine n che soddisfa le seguenti condizioni

'( )i

i

x h

i

x h

p x dx y+

−

=∫

con i=0,1,…n. I coefficienti di p’ risultano in generale diversi da quelli di p. Si consideri ora il

polinomio q(x) ottenuto calcolando l’integrale

( ) ( )1 1

0

( ) '( )1

i

i

x h ni ii

ix h

aq x p s ds x h x h

i

++ +

=−

= = + − − +∑∫

che risulta ancora un polinomio di ordine n in quanto i termini di grado n+1 si semplificano tra

loro.

In base alle condizioni imposte su p’ e alla definizione di q si ha ( )i iq x y= per i=0,1,...n e per

l’unicità del polinomio interpolante si può concludere che q(x)=p(x) e quindi si ha l’equivalenza

delle due interpolazioni.

In due dimensioni l’interpolazione integrale va a determinare il polinomio di ordine 2n

00 10 0 01 11 1( , ) ... .... ....n n n nn n nnp x y a a x a x a y a xy a x y a x y= + + + + + + + + +

tale che

( , )ij

ij

Q

p x y dxdy z=∫ (2.25)

per i,j=0,1,…n ove con ijQ si indica il quadrato di lato l centrato nel punto ( , )i jx y a cui è

associato il valore ijz . Dalla (2.25) si ottiene

1 12 2

2 2

, 0 1 1

h k

h k

l lx yi j

l ln x y

ij hki j

x y

a zi j

+ ++ +

− −

=

=

+ +∑ (2.26)

con h,k=0,1,..n.

Indicato con a il vettore dei coefficienti di p(x,y) e con z quello dei valori ijz la (2.26) in forma

matriciale diventa

59

=Ma z (2.27)

Se i punti ( , )h kx y sono tra loro distinti la matrice M è non singolare e il sistema ammette una ed

una sola soluzione

−= 1a M z (2.28)

ossia esiste ed è unico il polinomio soluzione dell’interpolazione integrale.

Si nota che questo tipo di interpolazione non è invariante rispetto ad una rotazione degli assi di

riferimento.

2.2.4 Interpolazione bilineare e biquadratica

In questo lavoro si sono considerati metodi di interpolazione polinomiale in due dimensioni con

polinomi del primo e del secondo ordine poichè risultano essere i meno costosi da un punto di

vista computazionale e non danno origine a fenomeni di oscillazione.

Dal teorema 2.3 si ricava che per definire univocamente il polinomio interpolatore di ordine n in

una dimensione si necessita di n+1punti di collocazione distinti.

Se n=1 l’interpolazione è detta lineare e servono due punti 0x e 1x distinti a cui assegnare dei

valori. La base dei polinomi interpolanti di Lagrange è data da

10

0 1

01

1 0

x xL

x x

x xL

x x

−=−

−=−

(2.29)

da cui ∀x∈[ ]0 1,x x si ha

0 1( ) (1 )p x dx y dxy= − + (2.30)

con

0

1 0

x xdx

x x

−=−

(2.31)

La costruzione del polinomio interpolante una funzione f di più variabili viene di solito effettuata

sfruttando le proiezioni sui sottospazi.

Sia [ ] [ ], ,a b c d× il dominio rettangolare della funzione f da interpolare e sia

0 1

1 2

...

...

n

m

a x x x b

c y y y d

= < < < =

= < < < =

60

una partizione del dominio . Il polinomio interpolatore di Lagrange è dato da

0 0

( , ) ( , ) ( ) ( )n m

i j i ji j

L x y f x y x yϕ ψ= =

=∑∑ (2.32)

dove ( ), ( )i jx yϕ ψ sono polinomi di grado n e m soddisfacenti le condizioni

( ) , 0,1,...

( ) j, 0,1,...i k ik

j k ik

x i k n

y k m

ϕ δψ δ

= == =

(2.33)

( , )L x y può essere generato come risultato dell’applicazione della trasformazione

0

: ( , ) ( )n

x i ii

P g g x y xϕ=

→∑ (2.34)

sull’immagine ( )yP f della trasformazione

0

: ( , ) ( )m

y j jj

P d d x y yψ=

→∑ (2.35)

sulla funzione f cioè

( ( ))x yL P P f= (2.36)

xP e yP sono proiezioni e per la scelta delle ( ), ( )i jx yϕ ψ sono commutative,quindi variando

l’ordine con cui si interpola, prima per righe e poi per colonne o viceversa, si ottiene sempre lo

stesso risultato finale.

Se n=m=1 l’interpolazione è detta bilineare. I polinomi di Lagrange hanno la forma

0

1

( ) 1

( )

x dx

x dx

ϕϕ

= −=

(2.37)

con

0

1 0

x xdx

x x

−=−

(2.38)

e

0

1

( ) 1

( )

y dy

y dy

ψψ

= −=

(2.39)

con

0

1 0

y ydy

y y

−=−

(2.40)

da cui

0 0 1 1( ( , )) ( ) ( , ) ( ) ( , )xP f x y x f x y x f x yϕ ϕ= + (2.41)

61

0 0 1 1( ( , )) ( ) ( , ) ( ) ( , ) i=0,1y i i iP f x y y f x y y f x yψ ψ= + (2.42)

e unendo le espressioni

[ ] [ ]0 0 0 0 1 0 1 1 0 1 0 1 1 1( , ) ( ) ( , ) ( ) ( , ) ( ) ( ) ( , ) ( ) ( , )L x y y f x y y f x y x y f x y y f x yϕ ψ ψ ϕ ψ ψ= + + + (2.43)

[ ] [ ]0 0 0 1 1 0 1 1( , ) (1 ) (1 ) ( , ) ( , ) (1 ) ( , ) ( , )L x y dx dy f x y dyf x y dx dy f x y dyf x y= − − + + − + (2.44)

L’esistenza e unicità del polinomio è garantita dall’esistenza e unicità dei polinomi interpolatori

lungo le righe e le colonne.

(0,0) (x,0) (1,1) (x,y) • (1,0) (x,1) (1,1)

Figura 2. 1 L’interpolazione bilineare operata interpolando linearmente per righe e poi per colonne

Se n=m=2 l’interpolazione è detta biquadratica. In questo caso servono nove punti di

interpolazione quindi una partizione del dominio [ ] [ ], ,a b c d× data da

0 1 2

1 2 2

a x x x b

c y y y d

= < < =

= < < = (2.45)

62

(0,0) a (1,0) (2,0) (0,1) b (1,1) (2,1) • (x,y) (0,2) c (1,2) (2,2)

Figura 2. 2 L’interpolazione biquadratica coinvolge nove punti di interpolazione per semplicità in figura sono indicati solo i loro indici. Si interpola prima lungo le righe i punti a,b,c poi questi lungo la colonna in blu con parabole

I polinomi interpolatori di Lagrange hanno equazioni

( )( )( )( )( )( )

( )( )( )( )

( )( )

1 20

0 1 0 2

0 21

1 0 1 2

1 02

2 1 2 0

( )

( )

( )

x x x xx

x x x x

x x x xx

x x x x

x x x xx

x x x x

ϕ

ϕ

ϕ

− −=

− −

− −=

− −

− −=

− −

(2.46)

e

( )( )( )( )( )( )

( )( )( )( )

( )( )

1 20

0 1 0 2

0 21

1 0 1 2

1 02

2 1 2 0

( )

( )

( )

y y y yy

y y y y

y y y yy

y y y y

y y y yy

y y y y

ψ

ψ

ψ

− −=

− −

− −=

− −

− −=

− −

(2.47)

e il polinomio interpolatore diventa

, 0

( , ) ( ) ( )n

ij i ji j

p x y z x yϕ ψ=

= ∑ (2.48)

con ijz valore associato al punto ( , )i jx y . Anche in questo caso l’esistenza e unicità del polinomio

è data dall’esistenza e unicità dei polinomi lungo le righe e le colonne e la commutatività delle

63

proiezioni garantisce che scambiando l’ordine di interpolazione per righe e colonne il risultato non

cambia.

2.3 Interpolazione dei livelli di grigio

Come si è già accennato nell’introduzione, si è scelto di interpolare le intensità di ciascuno dei

colori fondamentali (rosso, verde, blu) in modo indipendente rispetto alle intensità degli altri due

colori . In questo modo ci si è ridotti a tre problemi indipendenti di interpolazione dei livelli di

intensità luminosa (rispettivamente del rosso,del verde e del blu). Per ciascuno di essi dato un

punto di coordinate (p’,q’) dell’immagine di input si va a calcolare un valore approssimato

dell’intensità di grigio ottenuto interpolando le intensità dei punti dell’immagine appartenenti ad

un certo intorno di (p’,q’). Ovviamente, la scelta del numero e del tipo di punti che vengono

coinvolti dipende dal metodo di interpolazione usato. Molto importante è sottolineare che, poichè i

livelli di grigio sono numeri interi compresi tra 0 e 255, anche il valore approssimato deve

appartenere all’insieme

{ }: 0 255I n n= ∈ ≤ ≤ℕ

I metodi che verranno analizzati non sempre garantiscono che questa condizione sia soddisfatta ed

è quindi necessario imporre controlli . Se l’intensità luminosa è maggiore di 255 o minore di 0 la si

pone rispettivamente uguale a 255 o 0, se è un valore reale tra 0 e 255 la si arrotonda all’intero più

vicino. . L’arrotondamento dell’intensità di ciascun colore causa un errore che viene diffuso

scaricando l’errore stesso sui punti dell’immagine adiacenti al punto considerato. In questo caso si

è scelto per semplicità di disperdere l’errore solo in avanti cioè sul pixel successivo al punto

considerato appartenente alla stessa riga.

La scelta più semplice di interpolazione è chiamata interpolazione di ordine zero e consiste

nell’associare al punto (p’,q’) l’intensità del punto dell’immagine che ha distanza euclidea minore.

Questo algoritmo ha il vantaggio di essere computazionalmente semplice ma risultano evidenti

effetti di scalatura nelle immagini in cui il livello di grigio cambia significativamente da un pixel

all’altro e può portare ad un errore spaziale dell’ordine di 2

2 pixel per unità .

L’approssimazione può essere migliorata utilizzando l’interpolazione polinomiale .

Per mantenere costi computazionali accettabili ed evitare fenomeni di oscillazione

dell’interpolante si sono usati polinomi interpolanti del primo e del secondo ordine.

64

2.3.1 Interpolazione bilineare

Si consideri l’immagine di input F(p,q) e il punto (p’,q’) in cui si vuole calcolare il valore del

polinomio interpolante p(x,y) che, in base a quanto detto precedentemente sull’interpolazione

bilineare, ha equazione

[ ] [ ]0 0 0 1 1 0 1 1( , ) (1 ) (1 ) ( , ) ( , ) (1 ) ( , ) ( , )P x y dx dy f x y dyf x y dx dy f x y dyf x y= − − + + − + (2.49)

Per definire in modo univoco p devono essere assegnati quattro punti dell’immagine non allineati

tre a tre si considerano allora i punti di F(p,q) che costituiscono i quattro vertici del quadrato

contenente (p’,q’). Per comodità si considera un riferimento locale centrato nel centro del quadrato

C rispetto al quale i vertici hanno coordinate { }, 0.5,0.5

( , )i j

i j=−

poichè i punti dell’immagine formano

un reticolo di lato 1.

Il quadrato è così determinato : se p’-[p’]<0.5 si considera come vertice di coordinate (-0.5,0.5) il

punto con ascissa [p’] altrimenti con ascissa [p’]+1 e con ordinata [q’] se q’-[q’]<0.5 altrimenti

con ordinata [q’]+1.

Nell’immagine seguente è illustrata la scelta dei punti interpolati e il punto (p’,q’) che rispetto al

nuovo sistema di riferimento ha coordinate (x,y)

[ ] [ ][ ] [ ]

' ' 0.5 se p'- p' 0.5

' ' 1 se p'- p' 0.5

x p p

x p p

= − − <

= − − > (2.50)

[ ] [ ][ ] [ ]

' ' 0.5 se q'- q' 0.5

' ' 1 se q'- q' >0.5

y q q

y q q

= − − < = − −

(2.51)

65

(-0.5,-0.5) a (0.5,-0.5) • (x,y) (0,0) (-0.5,0.5) b (0.5,0.5)

Figura 2. 3 Le intersezioni delle linee rosse sono i punti di F(p,q), si interpolano prima per riga i punti intersezione della linea blu di ascissa x con le linee rosse orizzontali indicati con a e b e poi tali valori tra loro

Indicata con dx e dy le quantità

0.5

0.5

dx x

dy y

= += +

(2.52)

la (2.49) diventa

[ ]

[ ]( , ) (1 ) (1 ) ( 0.5, 0.5) ( 0.5,0.5)

+ (1 ) (0.5, 0.5) (0.5,0.5)

P x y dx dy I dyI

dx dy I dyI

= − − − − + − +

− − + (2.53)

L’interpolazione in orizzontale e verticale sono lineari ma la loro composizione risulta una

superficie non lineare che si adatta bene ai quattro pixel vicini .

2.3.2 Interpolazione biquadratica

Si consideri l’immagine di input F(p,q) e il punto (p’,q’) in cui si vuole calcolare il valore del

polinomio interpolante p(x,y) che, in base a quanto detto precedentemente sull’interpolazione

biquadratica, ha equazione

, 0

( , ) ( ) ( )n

ij i ji j

p x y I x yϕ ψ=

= ∑

dove i polinomi interpolanti di Lagrange ,i jϕ ψ sono definiti nelle (2.46) e (2.47) e ,i jI è

l’intensità luminosa nel punto ( , )i jx y .

66

Affinchè p sia definito univocamente come specificato prima si devono assegnare nove punti

nell’intorno di (p’,q’), si sceglie la griglia 3 3× dei nove punti più vicini e per comodità li si indica

con { }, 1,0,1

( , )i j i jx y

=−. Si considera un sistema di coordinate locali centrato in 0 0( , )x y per cui i

punti da interpolare diventano { }, 1,0,1

( , )i j

i j=−

e (p’,q’) ha coordinate (x,y) date da

[ ][ ]

' ' 0.5

' ' 0.5

x p p

y q q

= − −

= − − (2.54)

sostituendo si ottiene

1

, 1

( , ) ( ) ( )ij i ji j

p x y I x yϕ ψ=−

= ∑ (2.55)

con

1

0

1

( 1)( )

2( ) ( 1)( 1)

( 1)( )

2

z zz

z z z

z zz

ζ

ζ

ζ

−−=

= + − ++=

(2.56)

per z=x,y e ,ζ ϕ ψ= .

2.3.3 Interpolazione integrale

Si consideri l’immagine di input F(p,q) e il punto (p’,q’) in cui si vuole calcolare il valore che

l’interpolazione integrale gli associa. In base a quanto detto precedentemente poichè siamo nel

caso bidimensionale l’interpolazione integrale va a determinare il polinomio di equazione

2 2 2 2 2 200 10 20 01 11 21 02 12 22( , )p x y a a x a x a y a xy a x y a y a xy a x y= + + + + + + + +

tale che

( , )ij

ij

Q

p x y dxdy z=∫ (2.57)

per i,j=-1,0,1, ove con ijQ si indicano i nove quadrati di lato 1 centrati nei nove punti di F(p,q) più

vicini a (p’,q’). Poichè si fa riferimento ad un sistema di coordinate centrato nel punto di F(p,q) il

cui quadrato contiene (p’,q’), i punti presi in esame hanno coordinate (i,j) con i,j=-1,0,1.

67

Figura 2. 4 In blu è segnato il quadrato ( ', ')p qQ di lato 1 centrato nel punto (p’,q’), in rosa i quadrati ijQ con

i,j=-1,0,1 centrati nei 9 punti di F(p,q) più vicini a (p’,q’).

Il valore ijz rappresenta l’intensità del grigio del punto (i,j).

Dalla (2.25) si ottiene

0.50.5 112

0.50.5

, 0 1 1

ji kh

jihk ij

h k

yxa z

h k

++ ++−−

=

=+ +∑ (2.58)

con i,j=-1,0,1.

Indicato con a il vettore dei coefficienti di p(x,y) e con z quello dei valori ijz la (2.58) in forma

matriciale diventa

=Ma z (2.59)

Poichè i nodi sono tra loro distinti la matrice M è non singolare e invertibile, il sistema ammette

una ed una sola soluzione

−= 1a M z (2.60)

con -1M matrice i cui elementi sono indicati in figura 2.5.

68

0.0017362 -0.0451390 0.0017361 -0.0451390 1.1736113 -0.0451390 0.0017362 -0.0451390 0.0017361

0.0208333 0.0000000 -0.0208333 -0.5416667 0.0000000 0.5416667 0.0208333 0.0000000 -0.0208333

-0.0208333 0.0416667 -0.0208333 0.5416667 -1.0833334 0.5416667 -0.0208333 0.0416667 -0.0208333

0.0208333 -0.5416667 0.0208333 0.0000000 0.0000000 0.0000000 -0.0208333 0.5416667 -0.0208333

0.2500000 0.0000000 -0.2500000 0.0000000 0.0000000 0.0000000 -0.2500000 0.0000000 0.2500000

-0.2500000 0.5000000 -0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 -0.5000000 0.2500000

-0.0208333 0.5416667 -0.0208333 0.0416667 -1.0833334 0.0416667 -0.0208333 0.5416667 -0.0208333

-0.2500000 0.0000000 0.2500000 0.5000000 0.0000000 -0.5000000 -0.2500000 0.0000000 0.2500000

0.2500000 -0.5000000 0.2500000 -0.5000000 1.0000000 -0.5000000 0.2500000 -0.5000000 0.2500000

Figura 2. 5 Matrice -1M

Determinati i coefficienti hka di p(x,y) il valore da associare alle coordinate (p’,q’) si ottiene calcolando l’integrale

( ', ')

0.50.5 112 2

0.50.5

, 0 , 0 1 1p q

dydx kh

dyh k dxhk hk

h k h kQ

yxa x y dxdy a

h k

++ ++−−

= =

=+ +∑ ∑∫ (2.61)

ove con ( ', ')p qQ si indica il quadrato di lato 1 centrato in (p’,q’) di coordinate dx,dy rispetto al

sistema cartesiano considerato.

69

2.4 Analisi dei risultati

Dalle prove eseguite interpolando con i metodi di interpolazione sopra descritti si sono ottenuti

risultati soddisfacenti utilizzando l’interpolazione biquadratica e quella integrale. Come si può

osservare dalle immagini sottostanti di questi due metodi quello integrale risulta dare una

definizione dei contorni leggermente migliore, ma ha lo svantaggio di essere più lento.

Figura 2. 6 Immagine originale Figura 2. 7 Immagine ottenuta usando l’interpolazione biquadratica

Figura 2. 8 Immagine ottenuta usando Figura 2. 9 Immagine interpolata con

l’interpolazione integrale Photoshop