Analisi comparativa di tecniche di registrazione applicate ad immagini di testa-collo e del seno in risonanza

Capitolo 1
Registrazione
Basi teoriche sul processo di registrazione
In questo capitolo vengono presentati i concetti base per la comprensione del
problema in modo da rendere chiaro l’uso dei tool elastix [2] e FSL[3] (che ver-
ranno presentati in seguito). La \Image Registration"  e un’importante strumen-
to nel campo delle immagini medicali. In molte situazioni cliniche le immagini
di un paziente sono realizzate allo scopo di analizzare la situazione temporale
del paziente. Queste immagini sono acquisite con, per esempio, scanner a raggi
X, scanner a risonanza magnetica , che si basa sull’interazione di un campo ma-
gnetico statico con le molecole d’acqua presenti all’interno di ciascun organo del
nostro corpo, Tomograa Computerizzata (o pi u semplicemente TAC), la quale
sfrutta l’interazione dei tessuti con radiazioni di opportuna lunghezza d’onda,
e scanner a ultrasuoni, per fornire la conoscienza dell’anatomia del soggetto[1].
La combinazione dei dati del paziente, mono-modali o multi-modali, spesso con-
tiene informazioni cliniche aggiuntive che non  e possibile avere con un’analisi
separata delle immagini. Come messo in evidenza da diversi specialisti, solo
attraverso la combinazione dei dati eterogenei acquisiti nei diversi studi  e pos-
sibile ottenere un reale guadagno a livello dell’informazione sulla quale i medici
possono basare le proprie valutazioni cliniche. A tal ne, deve essere trovata la
relazione spaziale fra le immagini. Il compito della registrazione  e quindi quello
di trovare una mappatura spaziale \uno ad uno" tra i voxel
1
di una immagine
e i voxel dell’altra immagine. Nelle prossime sezioni si introduce la descrizione
matematica del processo di registrazione.
1.1 Costruzione del processo di registrazione
Nel processo della registrazione sono convolte due immagini:
 xed image I
F
(x): immagine di partenza
 moving image I
M
(x): immagine deformata
1
Un voxel  e un elemento di volume che rappresenta un valore di intensit a di segnale o di
colore in uno spazio tridimensionale, analogamente al pixel che rappresenta un dato di un’im-
magine bidimensionale. I voxel vengono spesso usati come elemento base per la visualizzazione
e l’analisi di dati medici e scientici. Per semplicare, un voxel corrisponde ad un volume
corporeo che viene poi rappresentato da un pixel nell’immagine bidimensionale.
11
CAPITOLO 1. REGISTRAZIONE
Figura 1.1: Componenti basilari della registrazione
L’obiettivo della registrazione  e quello di trovare uno spostamento u(x) che
costruisca l’immagine I
M
(x + u(x)) spazialmente allineata a I
F
(x). Una for-
mulazione equivalente  e vedere la registrazione come un problema che consi-
ste nel trovare una trasformazione T(x) = x + u(x) che costruisca I
M
(T(x))
spazialmente allineata a I
F
(x). La qualit a dell’allineamento  e denita da una
misurazione della distanza o da una misura della somiglianzaS, come la somma
delle dierenze al quadrato (SSD), il rapporto di correlazione, o la misura di
informazione reciproca (mutual information measure: MI). Poich e il problema  e
malposto per la trasformazione non-rigida T, viene spesso introdotto un termi-
neP in modo da regolarizzare la mappatura. Comunemente il problema della
registrazione  e formulato come un problema di ottimizzazione dove la funzione
di costoC  e ridotta al minimo:
b
T = arg minC(T;I
F
;I
M
); con (1.1)
C(T;I
F
;I
M
) =S (T;I
F
;I
M
) + P(T) (1.2)
dove   e una costante per regolare l’andamento di trasformazioni nonrigide.
Per risolvere il problema della minimizzazione esistono due possibili appocci:
parametrico e non parametrico. In questa tesi sono presentate solo le tecniche
parametriche in quanto i tool utilizzati si basano su questo tipo di approccio.
Nel metodo parametrico il numero di possibili trasformazioni  e limitato dall’in-
troduzione di parametri (modelli) nella formula T. Il problema originale quindi
si ottimizza e diventa:
b
T
 = arg min
T C(T
 ;I
F
;I
M
) (1.3)
dove sta ad indicare che la trasformazione  e parametrizzata: il vettore  con-
terr a quindi tutti i parametri introdotti. Per esempio, quando la trasformazione
 e modellata in 2D (trasformazione rigida), il vettore dei parametri  conterr a
un angolo di rotazione e una traslazione nelle direzioni x e y. Miglioriamo la
scrittura dell’equazione 1.3:
b  = arg min
 C( ;I
F
;I
M
) (1.4)
da questa equazione diventa chiaro che il problema originale 1.1  e semplicato.
In gura 1.1 si vedono le componenti generali della registrazione parametrica in
12
CAPITOLO 1. REGISTRAZIONE
uno schema a blocchi. Diverse componenti possono essere ritrovate dalle equa-
zioni 1.1 e 1.4; Alcune invece verranno introdotte nel seguito del capitolo. Come
elemento di base abbiamo le immagini. Il signicato di immagine in questo parti-
colare procedimento assume un’importanza fondamentale. Le immagini trattate
vengono generate dall’acquisizione sica attraverso vari metodi. Per questo mo-
tivo devono essere ricche di informazioni in modo da immagazzinare esattamente
il collegamento fra spazio sico (reale) e quello dei voxel digitalizzati. Poi abbia-
mo la funzione di costoC, o \metrica", che denisce la qualit a dell’allineamento.
Come accennato in precedenza, il costo della funzione consiste in una misura
di somiglianzaS e un termine regolatoreP. La denizione della misura di so-
miglianza introduce il componente \sampler". La procedura di ottimizzazione
nello schema a blocchi  e introdotto invece per risolvere il problema 1.4. Il mec-
canismo di ottimizzazione necessita della componente di interpolazione basata
sul valore dell’intensit a in quanto il processo di unione non avviene per voxel
(quindi non si basa sulla posizione dell’immagine). Un’altra cosa implicita nello
schema  e l’utilizzo di strategie multi-risoluzione per velocizzare la registrazione
che verr a spiegata in questo capitolo.
1.2 Metriche
Di seguito sono descritte in modo indicativo alcune metriche. Fra parentesi  e
indicato anche il relativo comando con il tool elastix:
 Somma delle Dierenze Quadrate (SSD) : (AdvancedMeanSquares)
SSD( ;I
F
;I
M
) =
1
j F
j
X
xi2 F
(I
F
(x
i
) I
M
(T
 (x
i
)))
2
(1.5)
con  F
dominio della xed image ej F
j numero di voxel.
 Coeciente di Correlazione Normalizzato (NCC) : (AdvancedNormalizedCorrelation)
NCC( ;I
F
;I
M
) =
P
xi2 F
(I
F
(x
i
) I
F
)(I
M
(T
 (x
i
)) I
M
)
q
P
xi2 F
(I
F
(x
i
) I
F
)
2
P
xi2 F
(I
M
(T
 (x
i
)) I
M
)
2
(1.6)
con la media dei valori di grigio I
F
=
1
j F
j
P
xi2 F
(I
F
(x
i
)) e I
M
=
1
j F
j
P
xi2 F
(I
M
(T
 (x
i
)).
 Mutua Informazione (MI): (AdvancedMattesMutualInformation)
MI( ;I
F
;I
M
) =
X
m2L
M
X
f2L
F
p(f;m; ) log
2
 p(f;m; )
p
F
(f)p
M
(m; )
 (1.7)
dove L
F
e L
M
sono insiemi contenenti intensit a regolarizzate, p  e la pro-
babilit a congiunta discreta, e p
F
e p
M
sono le probabilit a marginali di
xed image e moving image .
 Mutua Informazione Normalizzata (NMI): (NormalizedMutualInformation)
La NMI  e denita da NMI = (H(I
F
) +H(I
M
))=H(I
F
;I
M
), con H che
denota l’entropia. Questa denizione pu o essere comparata con la de-
nizione di MI. Riscrivendo cos  MI in termini di entropia H: MI =
H(I
F
) +H(I
M
) H(I
F
;I
M
).
13
CAPITOLO 1. REGISTRAZIONE
 Kappa Static (KS): (AdvancedKappaStatic)
KS( ;I
F
;I
M
) =
2
P
xi2 F
1
I
F
(xi)>0;I
M
(T (xi))>0
P
xi2 F
1
I
F
(xi)>0
+ 1
I
M
(T (xi))>0
(1.8)
dove 1  e la funzione che vale 1 se le condizioni in apice sono vericate.
La misura SSD  e adatta per immagini con una distribuzione uguale delle
intensit a (es. per immagini catturate con la stessa modalit a di scanner). La
misura NNC  e meno rigorosa e consiste in una relazione lineare fra i valori di
intensit a di xed image con moving image. Pu o quindi essere usata spesso.
Anche la misura MI  e molto generale:  e solo una relazione fra le distribuzioni
probabilistiche di intensit a fra xed image e moving image. Non  e adatta per
immagini monomodali ma solo per quelle multimodali. La misura NMI funziona
bene anche su immagini monomodali al contrario di MI. Mentre per concludere
KS  e particolarmente adatta per la registrazione di immagini binarie (molto
utile per la registrazione di immagini segmentate).
1.3 Campionamento delle immagini
In formule come la 1.5 si trovano espressioni del tipo
P
xi2 F
questo signica che
 e necessario un ciclo in ogni punto dell’immagine xed. La fase di campionamen-
to serve per limitare computazionalmente questo processo. Il campionamento
pu o essere eettuato in diversi modi, tra i quali presentiamo i pi u utilizzati:
 Full: il campionamento full semplicemente seleziona tutti i voxel con
coordinate x
i
nella xed image .
 Grid: il campionamento grid denisce una griglia regolare sull’imma-
gine xed e seleziona le coordinate x
i
sopra la griglia. In questo mo-
do il campionamento Grid eettua un sottocampionamento dell’immagi-
ne xed. La grandezza della griglia (o equivalentemente del fattore di
sottocampionamento) viene denita dall’utente.
 Random: un campionamento Random seleziona casualmente un numero
denito dall’utente di voxel dell’immagine xed, le cui coordinate forma-
no le x
i
. Tutti i voxel hanno la stessa chance di essere selezionati. Un
campione pu o essere selezionato pi u volte.
 Random Coordinate: uguale al Random tuttavia le coordinate casuali
non sono limitate dalla posizione dei voxel. Coordinate fra voxel possono
anche essere selezionate. I valori I
F
(x
i
) di quelle locazioni fra pi u voxel
saranno calcolati per interpolazione.
Mentre a prima vista il full sampler sembra la scelta pi u ovvia, in pratica
non  e sempre vero in quanto richiede un costo computazionale pi u elevato per
immagini molto grandi.
1.4 Interpolazione
Come detto in precedenza, durante la fase di ottimizzazione del valoreI
M
(T
 (x))
la valutazione non viene eettuata sulla base della posizione dei voxel. Questo
14
CAPITOLO 1. REGISTRAZIONE
signica che  e necessaria una fase di interpolazione sull’intensit a nelle posizioni
che risultano essere fra pi u voxel. Esistono diversi metodi per l’interpolazione,
che variano in velocit a e qualit a. Di seguito sono riportati alcuni esempi (come
in precedenza fra parentesi  e presente il comando elastix):
 Nearest neighbour: (NearestNeighborInterpolator) Questa  e la tecnica
pi u semplice. Richiede poche risorse ma  e di bassa qualit a. Come lette-
ralmente esprime il metodo, il valore di ritorno sar a il valore dell’intensit a
del voxel pi u vicino.
 Linear: (LinearInterpolator) Il valore di ritorno  e una media pesata di
tutti i voxel vicini, con la distanza di ogni voxel presa come peso.
 N -th order B-spline: (BSplineInterpolator o BSplineInterpolatorFloat)
Pi u alto  e l’ordine maggiore sar a la qualit a del risultato, ma anche mag-
giore sar a la quantit a di tempo di computazione. Di fatto l’interpolazione
nearest neighbour (grado del polinomio N=0) e l’interpolazione lineare
(grado del polinomio N=1) vengono ancora molto utilizzate per la loro
velocit a di esecuzione.
Durante la fase di registrazione l’interpolazione migliore  e quella lineare (gra-
do del polinomio N=1) in quanto  e l’interpolazione che ore il miglior trade-o
fr a qualit a e velocit a. Per generare il risultato nale invece, i.e. la deformazione
risultato dalla registrazione,  e richiesto tipicamente l’utilizzo di un alto grado
di interpolazione (es. grado polinomiale N=3).
1.5 Trasformazioni
I modelli di trasformazione usati in T
 determinano che tipo di deformazione
possiamo trattare fra l’immagine xed e l’immagine moving. Di seguito sono
riportate le trasformazioni in ordine di \essibilit a":
 Translation: (TranslationTransform)
T
 (x) = x + t (1.9)
con t vettore di traslazione. Il vettore dei parametri sar a semplicemente
denito da  = t.
 Rigid: (EulerTransform)
T
 (x) = R(x c) + t + c (1.10)
con R matrice di rotazione e c centro di rotazione. L’immagine  e trattata
come un corpo rigido: pu o essere traslata e ruotata ma non pu o essere
scalata.
 Similarity: (SimilarityTransform)
T
 (x) =sR(x c) + t + c (1.11)
con s parametro di scalatura. Questo signica che l’immagine  e trattata
come un oggetto che pu o essere traslato, ruotato e scalato.
15
CAPITOLO 1. REGISTRAZIONE
 Ane : (AneTransform)
T
 (x) = A(x c) + t + c (1.12)
dove A  e una matrice senza restrizioni.
 B-splines: (BSplineTransform)
Per la categoria delle trasformazioni \non rigide" le B-Spline [4] sono
spesso utilizzate come una parametrizzazione:
T
 (x) = x +
X
x
k
2Nx
p
k
 3
(x x
k
) (1.13)
con x
k
i punti di controllo,  3
il cubo della B-Spline polinomiale mul-
tidimensionale, p
k
il vettore dei coecienti B-Spline (in senso lato, gli
spostamenti dei punti di controllo), eN
x
l’insieme di tutti i punti di con-
trollo nel supporto compatto della B-Spline in x. I punti di controllo
x
k
sono deniti su una griglia che viene applicata alla xed image . In
questi punti avremo le direzioni e i moduli di scalatura per eettuare la
modellazione della moving image .
 Thin-plate splines: (SplineKernelTransform) Thin-plate splines  e un’al-
tra trasformazione non-rigida. La trasformazione si basa su un’insieme di
punti di interesse nella xed image e nella moving image . Lo spostamento
dei punti di interesse d
k
=x
m
 x
f
forma il vettore dei parametri  . La
posizione dei punti di interesse nell’immagine xed  e data dall’utente. La
trasformazione  e espressa come somma di una componente ane e di una
componente \non-rigida":
T
 (x) = x + Ax + t +
X
x
k
c
k
G(x c
x
k
) (1.14)
dove gli x
k
sono le posizioni dei punti di interesse nell’immagine xed,
G(r) la funzione di base e c
k
sono i coecienti corrispondenti ad ogni
punto di interesse. I coecienti c
k
e gli elementi A e t sono costruiti per
lo spostamento di ogni punto di interesse. La scelta specica per ogni
funzione base G(r) determina il \comportamento sico".
1.6 Ottimizzazioni
Per risolvere il problema dell’ottimizzazione 1.4, i.e. per ottimizzare il vetto-
re dei parametri di trasformazione, comunemente si impiega una strategia di
ottimizzazione iterativa:
 k+1
= k
+a
k
d
k
(1.15)
con d
k
la \direzione cercata" al passo k e a
k
fattore scalare di guadagno che
controlla la grandezza del passo nella direzione cercata. Di seguito sono illustrati
due importanti metodi:
 Gradient descent (GD): (StandardGradientDescent o RegularStepGradientDescent)
16
CAPITOLO 1. REGISTRAZIONE
questo metodo prende la direzione cercata come il negativo del gradiente
della funzione di costo:
 k+1
= k
 a
k
g( k
) (1.16)
con g( k
) =
@C
@ valutata nella posizione corrispondente k
. Diverse scelte
esistono per il fattore di guadagno a
k
.
 Robbins-Monro (RM): (StandardGradientDescent o FiniteDierenceGradientDescent)
La ottimizzazione RM calcola la derivata della funzione di costo g( k
)
con una approssimazione ~ g
k
 k+1
= k
 a
k
~ g
k
(1.17)
tale approssimazione  e potenzialmente pi u veloce nella computazione. Na-
turalmente per funzionare l’errore assoluto fra g( k
) e ~ g
k
deve essere
minimo.
Di seguito riportiamo i comandi per altre strategie di ottimizzazione che pos-
sono essere utilizzate in elastix: FullSearch, ConjugateGradient, ConjugateGradientFRPR,
QuasiNewtonLBFGS, RSGDEachParameterApart, SimultaneousPerturbation,
CMAEvolutionStrategy.
1.7 Multi-Risoluzione
Si distinguono due strategie gerarchiche per la tecnica multirisoluzionale delle
immagini: riduzione della complessit a dei dati e riduzione della complessit a delle
trasformazioni.
1.7.1 Riduzione della complessit a dei dati
 E comune iniziare il processo di registrazione utilizzando immagini con una bassa
complessit a, es., immagini che hanno subito un ltraggio di smoothing. Questo
incrementa le chance di successo della registrazione. Consideriamo ora una serie
di immagini ottenute con un uso incrementale del ltro di smoothing. Se le im-
magini non sono solo ltrate dallo smoothing ma sono anche sotto-campionate, i
dati non hanno solo una complessit a minore ma sono eettivamente in quantit a
minore. Da ora in poi serie di immagini smooth e sotto-campionate le chiame-
remo piramidi. Esistono molti tipi di piramidi: piramidi di Gauss, piramidi di
Laplace, piramidi spline, piramidi wavelet ecc. Sicuramente la pi u comune  e la
piramide di Gauss:
1. Gaussian pyramid: (FixedRecursiveImagePyramid e MovingRecursiveImagePyramid)
Applica smoothing e sotto-campionamento
2. Gaussian scale space: (FixedSmoothingImagePyramid e MovingSmoothingImagePyramid)
Applica smoothing e non applica il sotto-campionamento
17