Keypoint e descrittori locali

1. Feature locali: l'idea e i requisiti

Per allineare due immagini (il problema lasciato aperto dal capitolo 6) la soluzione passa dalle feature locali. La pipeline è composta da tre fasi:

1. Localizzazione dei punti di interesse (keypoint);
2. Calcolo di un descrittore per ciascun punto;
3. Matching, cioè la corrispondenza tra descrittori di immagini diverse.

Le buone feature locali devono possedere quattro caratteristiche:

Per rendere l'approccio più robusto si usa solitamente un'analisi multi-scala, che permette di gestire meglio la variabilità delle immagini. Le scelte principali per la localizzazione sono Harris–Laplace, Laplacian of Gaussian, Difference of Gaussian e dense sampling; per i descrittori SIFT, SURF e BRIEF.

2. Quali punti scegliere: flat, edge e corner

I punti di interesse dovrebbero essere facilmente riconoscibili analizzando una piccola porzione di immagine nel loro intorno: uno spostamento anche piccolo della finestra di osservazione in qualsiasi direzione dovrebbe portare a un cambiamento notevole di intensità. Questo criterio divide i punti in tre classi:

• nelle regioni flat non si rilevano cambiamenti di intensità in nessuna direzione;
• in corrispondenza di edge si rilevano cambiamenti di intensità nella direzione ortogonale all'edge;
• nei corner si misurano cambiamenti significativi in tutte le direzioni.

Solo i corner sono candidati affidabili: un edge lascia ambigua la posizione lungo la sua direzione (il «problema dell'apertura»).

3. Il corner detector di Harris

Il corner detector di Harris [Harris88] individua i corner come i punti in cui si incontrano due edge, caratterizzati da un gradiente elevato in entrambe le direzioni (orizzontale e verticale). Il cambiamento di intensità all'interno di una finestra in seguito a una traslazione di u, v è:

E(u, v) = Σ_{x,y} w(x, y) · [ I(x+u, y+v) − I(x, y) ]²

dove w(x, y) è la window function (1 dentro la finestra, 0 fuori, oppure Gaussiana). Con l'espansione di Taylor del secondo ordine si arriva alla forma quadratica:

E(u, v) ≈ [u  v] · M · [u  v]ᵀ        M = Σ_{x,y} w(x, y) · ⎡ ∇x²    ∇x·∇y ⎤
                                                              ⎣ ∇x·∇y  ∇y²    ⎦

L'analisi degli autovalori λ₁, λ₂ di M permette di rilevare la presenza di corner (Tavola 7.2). Nella pratica si usa il coefficiente di cornerness, che non richiede il calcolo esplicito degli autovalori:

R = λ₁·λ₂ − α·(λ₁ + λ₂)² = det(M) − α·trace(M)²

con α costante normalmente in [0.04, 0.06]. Il valore di R è grande nei corner, piccolo nelle regioni flat e negativo in corrispondenza degli edge. L'algoritmo è quindi in tre passi:

1. calcolo del valore di cornerness R per ciascun pixel;
2. identificazione dei punti con score elevati (R > soglia);
3. selezione dei punti di massimo locale (soppressione dei non massimi).

Proprietà del detector di Harris

4. L'invarianza di scala: scale space, LoG e DoG

Date due immagini a scala diversa, regioni di dimensioni corrispondenti hanno aspetto simile; il problema è trovarle in modo indipendente nelle due immagini. L'idea: definire una funzione invariante per cambiamenti di scala (es. intensità media) sulla regione e prendere i punti di massimo di questa funzione: il rapporto tra le scale dei massimi corrisponde al fattore di scala tra le immagini.

Lo scale space è la rappresentazione multi-scala ottenuta applicando filtri Gaussiani con valori crescenti di σ — l'effetto visivo è una progressiva sfuocatura. La Gaussiana 2D:

G(σ) = 1/(2πσ²) · exp( −(x² + y²) / (2σ²) )

è un filtro separabile: il filtraggio 2D si sostituisce con due filtraggi 1D, uno per dimensione, con notevole risparmio di calcolo. Alcune strutture caratteristiche (edge, blob) si individuano con le derivate della Gaussiana. La funzione più interessante per la localizzazione di blob è il Laplacian of Gaussian, la derivata seconda:

LoG = ∂²G/∂x² + ∂²G/∂y²

Poiché la risposta del Laplaciano diminuisce all'aumentare della scala, per rendere la risposta invariante per scala si moltiplica per σ²:

LoG_norm(x, y) = σ² · ( ∂²G/∂x² + ∂²G/∂y² )

I localizzatori invarianti per scala visti a lezione:

I passaggi del LoG detector [Lind98]: 1) convoluzione dell'immagine con il Laplaciano normalizzato a diverse scale; 2) ricerca dei massimi nello scale-space della risposta del Laplaciano al quadrato.

5. SIFT: la localizzazione dei keypoint

La Scale Invariant Feature Transform (SIFT) [Lowe99] estrae punti salienti basandosi su un'approssimazione del LoG denominata Difference of Gaussian:

LoG = σ² · ( G_xx(x, y, σ) + G_yy(x, y, σ) )        DoG = G(x, y, kσ) − G(x, y, σ)

Si costruisce una rappresentazione multi-risoluzione convolvendo l'immagine con Gaussiane di scala variabile, L(x, y, kσ) = G(x, y, kσ) ∗ I(x, y), e i punti si localizzano come estremi della funzione DoG:

D(x, y, σ) = [ G(x, y, kσ) − G(x, y, σ) ] ∗ I(x, y) = L(x, y, kσ) − L(x, y, σ)

Perché la DoG? Due motivazioni: le immagini smoothed servono comunque per il calcolo dei descrittori, quindi la localizzazione richiede in aggiunta semplici sottrazioni; inoltre la DoG è un'ottima approssimazione della LoG, necessaria per la vera invarianza per scala.

La piramide e le ottave

Il calcolo dello scale-space e delle DoG è computazionalmente oneroso: si usa una struttura piramidale che sfrutta L(x, y, 2σ) = L′(x, y, σ), dove L′ è un livello dello scale-space a metà risoluzione. Lo scale-space si divide in ottave, ciascuna con mappe ottenute convolvendo con filtri di scala nel range [σ, 2σ]. Il numero di livelli per ottava è regolato dal parametro s (intervalli), con moltiplicatore di scala k = 2^(1/s); servono s + 3 livelli di convoluzioni gaussiane per individuare i punti su un'intera ottava.

Una volta calcolate tutte le DoG, i punti di interesse sono gli estremi (massimi e minimi) delle differenze: un pixel è estremo se è maggiore (o minore) di tutti gli 8 vicini nella stessa immagine e dei 9 + 9 vicini delle scale adiacenti — 26 confronti totali.

Localizzazione accurata e filtri

I punti individuati nello scale-space sono localizzati sui pixel, che per loro natura sono discreti e non rappresentano esattamente l'estremo dell'intorno: si effettua un'interpolazione tra i valori dell'intorno. Seguono una soglia di contrasto (si scartano gli estremi poco pronunciati) e il filtraggio degli edge (i punti su edge sono mal localizzati lungo la direzione dell'edge).

6. Orientazione e descrittore SIFT

Orientazione dei keypoint

Per costruire descrittori invarianti rispetto alla rotazione (planare) della scena, si assegna a ogni punto di interesse un'orientazione: si individua la direzione con gradiente maggiore. Si prende come riferimento l'immagine smoothed L alla scala più vicina a quella del punto e si calcolano magnitudo e orientazione del gradiente:

m(x, y) = √( [L(x+1, y) − L(x−1, y)]² + [L(x, y+1) − L(x, y−1)]² )
θ(x, y) = tan⁻¹( [L(x, y+1) − L(x, y−1)] / [L(x+1, y) − L(x−1, y)] )

Si costruisce un istogramma delle orientazioni a 36 bin, ciascuno che copre 10°. Ogni pixel contribuisce al bin della propria orientazione con un valore pari alla magnitudo pesata con una Gaussiana (σ = 1.5 volte la scala del punto) centrata nel punto stesso. Il massimo assoluto dell'istogramma è l'orientazione assegnata; se esistono massimi locali che superano l'80% del massimo assoluto, vengono generati nuovi keypoint con quelle orientazioni (stessa posizione e scala).

Il descrittore SIFT

I descrittori SIFT si calcolano sulle mappe di orientazione e magnitudo del gradiente in una finestra attorno al punto di interesse: la scala del keypoint seleziona le mappe, l'orientazione orienta la finestra (invarianza a rotazione). I valori di magnitudo vengono pesati con una Gaussiana con σ pari a metà della larghezza della finestra. I punti della regione sono partizionati in sottoregioni, e a ciascuna è associato un istogramma di 8 bin delle distribuzioni pesate della magnitudo per orientazione. Per robustezza si usa un'interpolazione trilineare per distribuire i valori del gradiente anche nei bin adiacenti.

Segue la normalizzazione per l'invarianza all'illuminazione: dapprima tutto il vettore viene normalizzato a 1 (variazioni lineari); poi, per eliminare l'effetto delle variazioni non lineari (es. saturazione della camera), si applica una soglia a 0.2 — si mantengono gli elementi con valore non superiore a 0.2 — e si normalizza di nuovo a 1. I parametri suggeriti: 4×4 istogrammi da 8 bin ciascuno → descrittore di dimensione 128.

7. SURF: box filter e immagine integrale

La tecnica SURF (Speeded-Up Robust Features) [Bay08] localizza e descrive i punti di interesse in modo estremamente efficiente mantenendo le proprietà di invarianza di SIFT: rotazioni, cambiamenti di scala e di illuminazione e prospettive limitate. L'idea è approssimare i filtri tradizionali (LoG) con box filter applicabili molto velocemente grazie all'immagine integrale.

L'immagine integrale

L'immagine integrale in posizione (r, c) è la somma dei valori dei pixel sopra e a sinistra di (r, c):

II(r, c) = Σ_{r′≤r, c′≤c} I(r′, c′)

Si calcola con una sola scansione dell'immagine originale usando le formule ricorsive S(r,c) = S(r,c−1) + I(r,c) e II(r,c) = II(r−1,c) + S(r,c), con S(r,−1) = 0 e II(−1,c) = 0. Con l'immagine integrale la somma in qualsiasi rettangolo richiede solo 4 letture: per il rettangolo D, II(4) + II(1) − II(2) − II(3). Le feature a 2, 3 o 4 rettangoli usano rispettivamente 6, 8 e 9 valori di riferimento.

Localizzazione e descrittore

La localizzazione si basa sulla matrice Hessiana H = [L_xx L_xy; L_xy L_yy], dove L_xx(x,y,σ) è il LoG dell'immagine; l'algoritmo ne usa un'approssimazione con box filter D_xx, D_yy, D_xy (scala più piccola: Gaussiana con σ = 1.2). Il determinante approssimato:

det(H_approx) = D_xx · D_yy − w · D_xy²       (w = 0.9)

La risposta ai filtri va normalizzata rispetto alla dimensione del filtro; i keypoint sono i massimi locali della mappa del determinante, con soppressione dei non massimi in un intorno 3×3×3 e interpolazione. A differenza di SIFT non serve una piramide di immagini: lo scale-space si ottiene ridimensionando i box filter e applicandoli all'immagine originale (stessa complessità a ogni scala), diviso in ottave con filtri di dimensione crescente.

L'orientazione si assegna calcolando la risposta pesata a due filtri Haar-like (in x e y) in un intorno circolare di dimensione 6s: ogni punto dell'intorno diventa un punto nello spazio 2D (risposta x, risposta y), e l'orientazione dominante si stima sommando le risposte su una serie di finestre di dimensione π/3.

Il descrittore parte da una regione quadrata centrata nel punto e orientata secondo l'orientazione assegnata, di dimensione 20s; la regione è divisa in 4×4 sottoregioni, in ciascuna si calcola la risposta ai filtri Haar-like su una griglia di punti (d_x orizzontale, d_y verticale), pesata con una Gaussiana. Per ogni sottoregione:

v = ( Σ d_x , Σ d_y , Σ |d_x| , Σ |d_y| )

La concatenazione su tutte le sottoregioni produce un descrittore di dimensione 64, normalizzato a vettore unitario per l'invarianza al contrasto. La versione estesa (dimensione 128) separa le somme di d_x e |d_x| per d_y < 0 e d_y > 0 (e viceversa per d_y).

8. FAST e BRIEF: l'efficienza estrema

FAST keypoints

Il keypoint FAST è un corner identificato sulla base delle variazioni locali di luminosità: un punto rappresenta un corner se ha un colore molto diverso (più chiaro o più scuro) da quello dei pixel circostanti. L'algoritmo:

1. consideriamo un pixel p e indichiamo con I_p la sua luminosità;
2. fissiamo una soglia T (es. T = 0.2 × I_p);
3. selezioniamo i 16 pixel in un intorno circolare di raggio 3;
4. p è un corner se ci sono N punti consecutivi con luminosità maggiore di I_p + T oppure minore di I_p − T (N tipicamente 12, 11 o 9 → FAST-12, FAST-11, FAST-9);
5. il procedimento si itera per tutti i pixel.

Per velocizzare la ricerca si controllano prima i pixel 1, 5, 9 e 13: il pixel p è candidato corner solo se almeno 3 di questi 4 soddisfano una delle due condizioni. I corner FAST sono spesso clusterizzati: dopo la detection serve una non-maxima suppression.

Varianti per raggiungere invarianza: per la rotazione si assegna un'orientazione al keypoint individuando un centroide di intensità: calcolati i momenti della patch m_pq = Σ_{x,y∈B} xᵖ yᵍ I(x, y) con p, q ∈ {0, 1}, il centroide è C = (m_10/m_00, m_01/m_00); la direzione del keypoint è quella che congiunge il centro geometrico O con C: θ = arctan(m_01 / m_10). Per la scala si cerca i corner dentro una image pyramid (immagini riscalate a diverse risoluzioni).

BRIEF: il descrittore binario

BRIEF (Binary Robust Independent Elementary Features) [Cal12] crea descrittori di dimensioni limitate, efficaci ed efficienti: una stringa binaria ottenuta valutando le differenze di intensità tra coppie di pixel nell'intorno del keypoint. Il test τ su una patch p di dimensione S×S:

τ(p; x, y) = 1 se p(x) < p(y)      altrimenti 0

dove p(x) è l'intensità in una versione smoothed di p in posizione x = (u, v)ᵀ. Si effettuano n_d test binari su coppie di pixel, ottenendo:

f_{n_d}(p) = Σ_{1 ≤ i ≤ n_d} 2^(i−1) · τ(p; x_i, y_i)

I valori di n_d suggeriti sono 128, 256 e 512. Il test τ considera solo le informazioni dei pixel in esame ed è molto sensibile al rumore: è preferibile uno smoothing preliminare (per efficienza, un box filter applicato al punto di interesse invece della Gaussiana). Le distribuzioni per il campionamento delle coppie: uniforme in [−S/2, S/2], Gaussiana isotropica, campionamento in due fasi (prima x_i, poi y_i con Gaussiana centrata in x_i), griglia polare, e il caso x_i = (0,0)ᵀ con y_i su una griglia polare approssimata.

9. Feature matching: dai descrittori ai match presunti

Il matching tra due immagini avviene in tre passi:

1. rilevare i keypoint e calcolarne i descrittori (SIFT, SURF, BRIEF, ...);
2. determinare un set di presunti match: ricerca esaustiva o strutture dati;
3. sulla base dei migliori match determinare la trasformazione ottimale (RANSAC).

Come confrontare due set di punti con descrittori locali? La ricerca esaustiva confronta ogni feature di un'immagine con tutte le feature dell'altra; per limitare la complessità si restringe il numero di confronti con hashing (una funzione hash organizza i descrittori in bucket di elementi «simili») o con index-based nearest neighbor search (es. KD-tree). Il problema di fondo: la presenza di outlier.

Come selezionare i match presunti

Per ciascun keypoint, il suo nearest neighbor è il keypoint a distanza euclidea minima. Molte feature non avranno un match corretto (background, rumore, occlusioni). Una soglia globale sulla distanza (SSD(patch₁, patch₂) < thr) non dà buoni risultati; è possibile scartare alcuni match errati osservando la distanza dal nearest neighbor (1-NN) in rapporto alla distanza dal secondo keypoint più simile (2-NN) — il cosiddetto ratio test.

10. RANSAC: la trasformazione ottimale in presenza di outlier

Come calcolare la trasformazione ottimale a partire dai match presunti, che contengono outlier? RANSAC (RAndom SAmple Consensus) è un metodo iterativo che stima un modello (es. i parametri di trasformazione) da un insieme di osservazioni in presenza di outlier. L'idea: i dati si dividono in inliers (dati la cui distribuzione è compatibile con i parametri del modello) e outliers (non rappresentati correttamente dal modello). Dati N punti (osservazioni) generati da un modello con parametri Θ, l'algoritmo cerca di stimare Θ:

Ripeti K volte:
  1. scegli n campioni casuali;
  2. adatta i parametri Θ agli n campioni;
  3. per ciascuno dei restanti N−n punti, calcola la distanza dal modello
     e conta il numero c di inliers.
Output: il modello Θ con c più elevato.

Ogni tipo di modello ha un minimal set: il numero minimo di esempi per calcolare il modello (es. linea → 2 punti). Per le trasformazioni di immagine il minimal set è di s coppie di punti: per l'omografia servono 4 coppie — la matrice H = [a b c; d e f; g h 1] si ricalcola a ogni iterazione e si contano gli inlier con SSD(p′ᵢ, H·pᵢ) < ε.

I parametri

• s: numero iniziale di punti nel minimal set (4 coppie per l'omografia);
• t: soglia di distanza che definisce un punto come inlier;
• K: numero di iterazioni necessarie per un modello affidabile.

Per la soglia: se il rumore nella posizione dei pattern è Gaussiano con deviazione standard σ, la distanza segue una distribuzione con g gradi di libertà; fissando t² = 3.84·σ², con probabilità del 95% tutti gli inlier avranno distanza d < t. Per il numero di iterazioni, con e = percentuale ipotetica di outlier e probabilità di successo p:

K > log(1 − p) / log(1 − (1 − e)^s)

Il numero K non dipende dal numero totale di punti. Con s = 4 e p = 0.99, la tabella vista a lezione dà: e = 0.05 → 3; e = 0.25 → 13; e = 0.40 → 34; e = 0.50 → 72 iterazioni.

Vantaggi e svantaggi

Verifica le tue conoscenze