Analisi comparativa di tecniche di registrazione applicate ad

Università degli Studi di Verona
FACOLTÀ DI SCIENZE MATEMATICHE, FISICHE E NATURALI
Corso di Laurea Triennale in Informatica Multimediale
Analisi comparativa di tecniche di registrazione applicate ad
immagini di testa-collo e del seno in risonanza
Autori:
Relatore:
Alessio Montagnini vr076776
Omar Zandonà vr077840
Ch.ma Prof. Gloria Menegaz
Correlatore:
Ing. Francesca Pizzorni Ferrarese
Anno Accademico 2009–2010
Indice
1 Registrazione
1.1 Costruzione del processo di registrazione . . . . . . . . .
1.2 Metriche . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Campionamento delle immagini . . . . . . . . . . . . . .
1.4 Interpolazione . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Trasformazioni . . . . . . . . . . . . . . . . . . . . . . .
1.6 Ottimizzazioni . . . . . . . . . . . . . . . . . . . . . . .
1.7 Multi-Risoluzione . . . . . . . . . . . . . . . . . . . . . .
1.7.1 Riduzione della complessità dei dati . . . . . . .
1.7.2 Riduzione della complessità delle trasformazioni
.
.
.
.
.
.
.
.
.
2 La Risonanza magnetica e metodi di acquisizione delle
gini
2.1 Risonanza magnetica nucleare . . . . . . . . . . . . . . . .
2.1.1 Tempo di rilassamento longitudinale-T1 . . . . . .
2.1.2 Tempo di rilassamento trasversale o rilassamento
spin-T2 . . . . . . . . . . . . . . . . . . . . . . . .
2.2 La tecnica di diffusione . . . . . . . . . . . . . . . . . . . .
2.3 DWI: Diffusion-Wieghted-Imaging . . . . . . . . . . . . .
2.4 DTI: Diffusion-Tensor-Imaging . . . . . . . . . . . . . . .
2.5 DCE: Diffusion-Contrast-Enhanched . . . . . . . . . . . .
2.6 Metodi a singola modalità verso metodi multimodali . . .
2.7 Misure di similitudine per la registrazione delle immagini
3 Strumenti ed algoritmi utilizzati
3.1 Standard e formati di immagini mediche .
3.2 Il programma 3D Slicer . . . . . . . . . .
3.3 Elastix . . . . . . . . . . . . . . . . . . . .
3.4 Fsl . . . . . . . . . . . . . . . . . . . . . .
3.5 Filtro Anisotropico . . . . . . . . . . . . .
3.5.1 Algoritmi di Perona e Malik . . . .
3.5.2 Codice . . . . . . . . . . . . . . . .
3.6 Region Growing . . . . . . . . . . . . . . .
3.7 Strumento per l’analisi oggettiva dei dati
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
11
11
13
14
14
15
16
17
17
18
imma19
. . . . 19
. . . . 19
spin. . . . 20
. . . . 20
. . . . 21
. . . . 22
. . . . 22
. . . . 23
. . . . 23
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
25
25
26
28
32
33
33
36
37
38
INDICE
4 Realizzazione e sviluppo
4.1 Immagini del seno . . . . .
4.1.1 Traccia di Sviluppo .
4.1.2 Risultati . . . . . . .
4.2 Immagini di testa-collo . . .
4.2.1 Traccia di Sviluppo .
4.2.2 Risultati . . . . . . .
5
.
.
.
.
.
.
.
.
.
.
.
.
Conclusioni
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
41
41
43
50
62
63
67
79
A Codice del filtro anisotropico 3D
83
Bibliografia
87
4
Elenco delle tabelle
3.1
3.2
Metriche disponibili in elastix . . . . . . . . . . . . . . . . . . .
Interpolazioni disponibili in elastix . . . . . . . . . . . . . . . . .
4.1
4.2
4.3
4.4
Tabella generale dei parametri . . . . . . . . . . . . . . . . . . .
Raffinazioni successive in elastix . . . . . . . . . . . . . . . . . .
Parametri utilizzati in FSL . . . . . . . . . . . . . . . . . . . . .
Risultati ottenuti nella registrazione Rigida in elastix, corrispondenti alle prime 12 prove riportate in tabella 4.1 . . . . . . . . .
4.5 Risultati ottenuti nella registrazione Affine in elastix, corrispondenti alle prove da 13 a 24 riportate in tabella 4.1 . . . . . . . .
4.6 Risultati ottenuti nella trasformazione BSpline in elastix . . . .
4.7 Risultati ottenuti nella fase di raffinamento dei risultati in elastix
4.8 Risultati ottenuti in FSL . . . . . . . . . . . . . . . . . . . . . .
4.9 Tabella risultati dopo eliminazione del rumore . . . . . . . . . . .
4.10 Risultati ottenuti in elastix . . . . . . . . . . . . . . . . . . . . .
4.11 Risultati ottenuti in elastix . . . . . . . . . . . . . . . . . . . . .
4.12 Risultati ottenuti in elastix . . . . . . . . . . . . . . . . . . . . .
5
29
29
44
45
46
51
53
55
57
58
59
72
74
76
Elenco delle figure
1.1
Componenti basilari della registrazione . . . . . . . . . . . . . . .
12
3.1
Tipiche immagini di risonanza magnetica del cervello. Da sinistra
verso destra, è raffigurata rispettivamente la ricostruzione assiale,
sagittale e coronale. . . . . . . . . . . . . . . . . . . . . . . . . .
Esempio di applicazione di un filtro anisotropico . . . . . . . . .
Funzioni g(∗) proposte da Perona e Malik . . . . . . . . . . . . .
Funzioni di flusso . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
33
34
35
3.2
3.3
3.4
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13
4.14
4.15
4.16
4.17
4.18
Immagine moving di tipo dwi, a destra si evidenzia la lesione
oggetto di idagine . . . . . . . . . . . . . . . . . . . . . . . . . . .
Immagine di riferimento di tipo dce-21 . . . . . . . . . . . . . . .
Nella parte evidenziata troviamo l’errore propagato dalla DWI. .
Registrazione dell’immagine simmetrica su DCE e successivo riutilizzo della trasformazione generata su immagine DWI. . . . . .
Registrazione dell’immagine simmetrica su DCE (utilizzando come passo intermedio di registrazione l’immagine DWI) e successivo riutilizzo della trasformazione generata su immagine DWI. .
Eliminazione dell’errore con applicazione di maschera generata utilizzando l’algoritmo di region growing sull’immagine DCE
filtrata con il filtro Anisotropico 3D . . . . . . . . . . . . . . . . .
Applicazione del filtro anisotropico sull immagine dce-21 . . . . .
Maschera dell’immagine dce-21 . . . . . . . . . . . . . . . . . . .
Grafico in norma 2: parametrizzazioni rigide utilizzando elastix
nelle ascisse ed errore nelle ordinate . . . . . . . . . . . . . . . .
Registrazione rigida . . . . . . . . . . . . . . . . . . . . . . . . .
Confronto tra registrazione rigida e volume di riferimento . . . .
Grafico in norma 2: parametrizzazioni affine utilizzando elastix
nelle ascisse ed errore nelle ordinate . . . . . . . . . . . . . . . .
Registrazione affine . . . . . . . . . . . . . . . . . . . . . . . . . .
Confronto tra registrazione affine e volume di riferimento . . . .
Grafico in norma 2: tutte le parametrizzazioni con elastix nelle
ascisse ed errore nelle ordinate . . . . . . . . . . . . . . . . . . .
Registrazione BSpline . . . . . . . . . . . . . . . . . . . . . . . .
Grafico degli errori in norma 2 dei test con “flirt”: numero delle
parametrizzazioni nelle ascisse ed errore nelle ordinate . . . . . .
Registrazione affine con applicazione della maschera per l’eliminazione del rumore . . . . . . . . . . . . . . . . . . . . . . . . . .
7
42
43
46
47
47
48
49
50
51
52
52
53
54
54
56
56
59
60
ELENCO DELLE FIGURE
4.19 Confronto tra registrazione con eliminazione del rumore e volume
di riferimento . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.20 Grafico degli errori: 8 immagini dce dello stesso paziente nell’asse
delle ascisse e valore fra 0 e 1 di errore nell’asse delle ordinate. .
4.21 Risultato utilizzando la parametrizzazione Rigida6 del seno della
2◦ paziente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.22 Confronto fra dce e risultato utilizzando la parametrizzazione
Rigida6 del seno della 2◦ paziente. . . . . . . . . . . . . . . . . .
4.23 Acquizione testa e collo DWI . . . . . . . . . . . . . . . . . . . .
4.24 Acquizione testa e collo T2 . . . . . . . . . . . . . . . . . . . . .
4.25 Acquizione testa e collo DCE . . . . . . . . . . . . . . . . . . . .
4.26 Primo approccio di registrazione utilizzando volume intermedio
T1 (T2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.27 Secondo approccio di registrazione utilizzando volume intermedio
T1 (T2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.28 Registrazione Affine DWI to T2 . . . . . . . . . . . . . . . . . . .
4.29 Confronto fra registrazione DWI-T2 e T2 . . . . . . . . . . . . .
4.30 Registrazione Affine T2 to DCE . . . . . . . . . . . . . . . . . . .
4.31 Confronto fra registrazione T2 e DCE . . . . . . . . . . . . . . .
4.32 Registrazione Affine DCE to T2 . . . . . . . . . . . . . . . . . . .
4.33 Grafico degli errori in norma 2 dei test con elastix : numero delle
parametrizzazioni nelle ascisse ed errore nelle ordinate . . . . . .
4.34 Trasformazione Rigida . . . . . . . . . . . . . . . . . . . . . . . .
4.35 Confronto fra trasformazione rigida-rigida DWI-DCE e DCE . .
4.36 Grafico degli errori in norma 2 dei test con elastix : numero delle
parametrizzazioni nelle ascisse ed errore nelle ordinate . . . . . .
4.37 Trasformazione Affine-Affine . . . . . . . . . . . . . . . . . . . . .
4.38 Confronto fra trasformazione affine-affine DWI-DCE e DCE . . .
4.39 Grafico degli errori in norma 2 dei test con elastix : numero delle
parametrizzazioni nelle ascisse ed errore nelle ordinate . . . . . .
4.40 Trasformazione Affine-BSpline . . . . . . . . . . . . . . . . . . . .
4.41 Confronto fra trasformazione Affine-BSpline DWI-DCE e DCE .
8
61
62
62
63
64
64
65
66
66
68
68
69
69
70
72
73
73
74
75
75
76
77
77
Introduzione
L’elaborazione digitale delle immagini in medicina è un campo in continuo sviluppo. In particolare le tecniche di registrazione di tali immagini sta diventando
sempre più uno strumento di importanza rilevante sia dal punto di vista diagnostico che dal punto di vista scientifico. L’intento di questa tesi è in primo luogo
descrivere e studiare vie di realizzazione di tale procedimento sulle immagini di
risonanza magnetica del seno e del cervello umano. In secondo luogo confrontare
i diversi approcci utilizzati in relazione ai diversi tools adottati. Le acquisizioni
delle immagini tramite la risonanza magnetica o MRN stanno avendo interessanti sviluppi grazie all’introduzione delle tecniche di diffusione che permettono
di avere informazioni qualitative diverse di una stessa parte anatomica. Nel corso del nostro elaborato abbiamo utilizzato due diverse tecniche di acquisizione:
DWI e DCE, le cui peculiarità verranno in seguito analizzate.
Dopo aver descritto nel primo capitolo i principi basilari della registrazione
delle immagini da un punto di vista matematico e da un punto di vista implementativo passeremo alla descrizione dei fondamenti della risonanza magnetica,
delle diverse acquisizioni e dell’imaging medico attuale. Successivamente elencheremo gli strumenti e i tool utilizzati nella realizzazione della nostra esperienza
andando a considerare per ognuno di essi efficienza e complessità. Concluderemo
nel capitolo “Realizzazione e Sviluppo” con la presentazione di alcuni risultati (i
più significativi) ottenuti. I risultati saranno accompagnati da un’analisi completa, da una descrizione del codice utilizzato e da una valutazione matematica
sviluppata per ottenere un giudizio formale che sarà ripreso e ampliato nelle
conclusioni finali.
9
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 verranno presentati in seguito). La “Image Registration” è un’importante strumento 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 magnetico statico con le molecole d’acqua presenti all’interno di ciascun organo del
nostro corpo, Tomografia Computerizzata (o più 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 contiene informazioni cliniche aggiuntive che non è 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 è possibile ottenere un reale guadagno a livello dell’informazione sulla quale i medici
possono basare le proprie valutazioni cliniche. A tal fine, deve essere trovata la
relazione spaziale fra le immagini. Il compito della registrazione è quindi quello
di trovare una mappatura spaziale “uno ad uno” tra i voxel1 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:
• fixed image IF (x): immagine di partenza
• moving image IM (x): immagine deformata
1 Un
voxel è un elemento di volume che rappresenta un valore di intensità di segnale o di
colore in uno spazio tridimensionale, analogamente al pixel che rappresenta un dato di un’immagine bidimensionale. I voxel vengono spesso usati come elemento base per la visualizzazione
e l’analisi di dati medici e scientifici. Per semplificare, 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 è quello di trovare uno spostamento u(x) che
costruisca l’immagine IM (x + u(x)) spazialmente allineata a IF (x). Una formulazione equivalente è vedere la registrazione come un problema che consiste nel trovare una trasformazione T(x) = x + u(x) che costruisca IM (T(x))
spazialmente allineata a IF (x). La qualità dell’allineamento è definita da una
misurazione della distanza o da una misura della somiglianza S, come la somma
delle differenze al quadrato (SSD), il rapporto di correlazione, o la misura di
informazione reciproca (mutual information measure: MI). Poichè il problema è
malposto per la trasformazione non-rigida T, viene spesso introdotto un termine P in modo da regolarizzare la mappatura. Comunemente il problema della
registrazione è formulato come un problema di ottimizzazione dove la funzione
di costo C è ridotta al minimo:
b = arg min C(T; IF , IM ),
T
con
C(T; IF , IM ) = −S(T; IF , IM ) + γP(T)
(1.1)
(1.2)
dove γ è 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 è limitato dall’introduzione di parametri (modelli) nella formula T. Il problema originale quindi
si ottimizza e diventa:
b µ = arg min C(Tµ ; IF , IM )
T
Tµ
(1.3)
dove µ sta ad indicare che la trasformazione è parametrizzata: il vettore µ conterrà quindi tutti i parametri introdotti. Per esempio, quando la trasformazione
è modellata in 2D (trasformazione rigida), il vettore dei parametri µ conterrà
un angolo di rotazione e una traslazione nelle direzioni x e y. Miglioriamo la
scrittura dell’equazione 1.3:
b = arg min C(µ; IF , IM )
µ
µ
(1.4)
da questa equazione diventa chiaro che il problema originale 1.1 è semplificato.
In figura 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 equazioni 1.1 e 1.4; Alcune invece verranno introdotte nel seguito del capitolo. Come
elemento di base abbiamo le immagini. Il significato di immagine in questo particolare procedimento assume un’importanza fondamentale. Le immagini trattate
vengono generate dall’acquisizione fisica attraverso vari metodi. Per questo motivo devono essere ricche di informazioni in modo da immagazzinare esattamente
il collegamento fra spazio fisico (reale) e quello dei voxel digitalizzati. Poi abbiamo la funzione di costo C, o “metrica”, che definisce la qualità dell’allineamento.
Come accennato in precedenza, il costo della funzione consiste in una misura
di somiglianza S e un termine regolatore P. La definizione della misura di somiglianza introduce il componente “sampler”. La procedura di ottimizzazione
nello schema a blocchi è introdotto invece per risolvere il problema 1.4. Il meccanismo di ottimizzazione necessita della componente di interpolazione basata
sul valore dell’intensità in quanto il processo di unione non avviene per voxel
(quindi non si basa sulla posizione dell’immagine). Un’altra cosa implicita nello
schema è l’utilizzo di strategie multi-risoluzione per velocizzare la registrazione
che verrà spiegata in questo capitolo.
1.2
Metriche
Di seguito sono descritte in modo indicativo alcune metriche. Fra parentesi è
indicato anche il relativo comando con il tool elastix:
• Somma delle Differenze Quadrate (SSD): (AdvancedMeanSquares)
X
1
(IF (xi ) − IM (Tµ (xi )))2
(1.5)
SSD(µ; IF , IM ) =
|ΩF |
xi ∈ΩF
con ΩF dominio della fixed image e |ΩF | numero di voxel.
• Coefficiente di Correlazione Normalizzato (NCC): (AdvancedNormalizedCorrelation)
P
xi ∈ΩF (IF (xi ) − IF )(IM (Tµ (xi )) − IM )
N CC(µ; IF , IM ) = qP
P
2
2
xi ∈ΩF (IF (xi ) − IF )
xi ∈ΩF (IM (Tµ (xi )) − IM )
(1.6)
P
con la media dei valori di grigio IF = |Ω1F | xi ∈ΩF (IF (xi )) e IM =
P
1
xi ∈ΩF (IM (Tµ (xi )).
|ΩF |
• Mutua Informazione (MI): (AdvancedMattesMutualInformation)
X X
p(f, m; µ)
M I(µ; IF , IM ) =
p(f, m; µ) log2
(1.7)
pF (f )pM (m; µ)
m∈LM f ∈LF
dove LF e LM sono insiemi contenenti intensità regolarizzate, p è la probabilità congiunta discreta, e pF e pM sono le probabilità marginali di
fixed image e moving image .
• Mutua Informazione Normalizzata (NMI): (NormalizedMutualInformation)
La NMI è definita da N M I = (H(IF ) + H(IM ))/H(IF , IM ), con H che
denota l’entropia. Questa definizione può essere comparata con la definizione di MI. Riscrivendo cosı̀ MI in termini di entropia H: M I =
H(IF ) + H(IM ) − H(IF , IM ).
13
CAPITOLO 1. REGISTRAZIONE
• Kappa Static (KS): (AdvancedKappaStatic)
P
2 xi ∈ΩF 1IF (xi )>0,IM (Tµ (xi ))>0
KS(µ; IF , IM ) = P
xi ∈ΩF 1IF (xi )>0 + 1IM (Tµ (xi ))>0
(1.8)
dove 1 è la funzione che vale 1 se le condizioni in apice sono verificate.
La misura SSD è adatta per immagini con una distribuzione uguale delle
intensità (es. per immagini catturate con la stessa modalità di scanner). La
misura NNC è meno rigorosa e consiste in una relazione lineare fra i valori di
intensità di fixed image con moving image. Può quindi essere usata spesso.
Anche la misura MI è molto generale: è solo una relazione fra le distribuzioni
probabilistiche di intensità fra fixed image e moving image. Non è 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 è particolarmente adatta per la registrazione di immagini binarie (molto
utile per la registrazione di immagini segmentate).
1.3
Campionamento delle immagini
P
In formule come la 1.5 si trovano espressioni del tipo xi ∈ΩF questo significa che
è necessario un ciclo in ogni punto dell’immagine fixed. La fase di campionamento serve per limitare computazionalmente questo processo. Il campionamento
può essere effettuato in diversi modi, tra i quali presentiamo i più utilizzati:
• Full: il campionamento full semplicemente seleziona tutti i voxel con
coordinate xi nella fixed image .
• Grid: il campionamento grid definisce una griglia regolare sull’immagine fixed e seleziona le coordinate xi sopra la griglia. In questo modo il campionamento Grid effettua un sottocampionamento dell’immagine fixed. La grandezza della griglia (o equivalentemente del fattore di
sottocampionamento) viene definita dall’utente.
• Random: un campionamento Random seleziona casualmente un numero
definito dall’utente di voxel dell’immagine fixed, le cui coordinate formano le xi . Tutti i voxel hanno la stessa chance di essere selezionati. Un
campione può essere selezionato più 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 IF (xi ) di quelle locazioni fra più voxel
saranno calcolati per interpolazione.
Mentre a prima vista il full sampler sembra la scelta più ovvia, in pratica
non è sempre vero in quanto richiede un costo computazionale più elevato per
immagini molto grandi.
1.4
Interpolazione
Come detto in precedenza, durante la fase di ottimizzazione del valore IM (Tµ (x))
la valutazione non viene effettuata sulla base della posizione dei voxel. Questo
14
CAPITOLO 1. REGISTRAZIONE
significa che è necessaria una fase di interpolazione sull’intensità nelle posizioni
che risultano essere fra più voxel. Esistono diversi metodi per l’interpolazione,
che variano in velocità e qualità. Di seguito sono riportati alcuni esempi (come
in precedenza fra parentesi è presente il comando elastix):
• Nearest neighbour: (NearestNeighborInterpolator) Questa è la tecnica
più semplice. Richiede poche risorse ma è di bassa qualità. Come letteralmente esprime il metodo, il valore di ritorno sarà il valore dell’intensità
del voxel più vicino.
• Linear: (LinearInterpolator) Il valore di ritorno è 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ù alto è l’ordine maggiore sarà la qualità del risultato, ma anche maggiore sarà la quantità 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à di esecuzione.
Durante la fase di registrazione l’interpolazione migliore è quella lineare (grado del polinomio N=1) in quanto è l’interpolazione che offre il miglior trade-off
frà qualità e velocità. Per generare il risultato finale invece, i.e. la deformazione
risultato dalla registrazione, è 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 fixed e l’immagine moving. Di seguito sono
riportate le trasformazioni in ordine di “flessibilità”:
• Translation: (TranslationTransform)
Tµ (x) = x + t
(1.9)
con t vettore di traslazione. Il vettore dei parametri sarà semplicemente
definito 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 è trattata
come un corpo rigido: può essere traslata e ruotata ma non può essere
scalata.
• Similarity: (SimilarityTransform)
Tµ (x) = sR(x − c) + t + c
(1.11)
con s parametro di scalatura. Questo significa che l’immagine è trattata
come un oggetto che può essere traslato, ruotato e scalato.
15
CAPITOLO 1. REGISTRAZIONE
• Affine: (AffineTransform)
Tµ (x) = A(x − c) + t + c
(1.12)
dove A è una matrice senza restrizioni.
• B-splines: (BSplineTransform)
Per la categoria delle trasformazioni “non rigide” le B-Spline [4] sono
spesso utilizzate come una parametrizzazione:
X
Tµ (x) = x +
pk β 3 (x − xk )
(1.13)
xk ∈Nx
con xk i punti di controllo, β 3 il cubo della B-Spline polinomiale multidimensionale, pk il vettore dei coefficienti B-Spline (in senso lato, gli
spostamenti dei punti di controllo), e Nx l’insieme di tutti i punti di controllo nel supporto compatto della B-Spline in x. I punti di controllo
xk sono definiti su una griglia che viene applicata alla fixed image . In
questi punti avremo le direzioni e i moduli di scalatura per effettuare la
modellazione della moving image .
• Thin-plate splines: (SplineKernelTransform) Thin-plate splines è un’altra trasformazione non-rigida. La trasformazione si basa su un’insieme di
punti di interesse nella fixed image e nella moving image . Lo spostamento
dei punti di interesse dk = xm − xf forma il vettore dei parametri µ. La
posizione dei punti di interesse nell’immagine fixed è data dall’utente. La
trasformazione è espressa come somma di una componente affine e di una
componente “non-rigida”:
X
Tµ (x) = x + Ax + t +
ck G(x − cfix
(1.14)
k )
xk
dove gli xk sono le posizioni dei punti di interesse nell’immagine fixed,
G(r) la funzione di base e ck sono i coefficienti corrispondenti ad ogni
punto di interesse. I coefficienti ck e gli elementi A e t sono costruiti per
lo spostamento di ogni punto di interesse. La scelta specifica per ogni
funzione base G(r) determina il “comportamento fisico”.
1.6
Ottimizzazioni
Per risolvere il problema dell’ottimizzazione 1.4, i.e. per ottimizzare il vettore dei parametri di trasformazione, comunemente si impiega una strategia di
ottimizzazione iterativa:
µk+1 = µk + ak dk
(1.15)
con dk la “direzione cercata” al passo k e ak 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 − ak g(µk )
(1.16)
∂C
con g(µk ) = ∂µ
valutata nella posizione corrispondente µk . Diverse scelte
esistono per il fattore di guadagno ak .
• Robbins-Monro (RM): (StandardGradientDescent o FiniteDifferenceGradientDescent)
La ottimizzazione RM calcola la derivata della funzione di costo g(µk )
con una approssimazione g̃k
µk+1 = µk − ak g̃k
(1.17)
tale approssimazione è potenzialmente più veloce nella computazione. Naturalmente 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 possono 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à dei dati e riduzione della complessità delle
trasformazioni.
1.7.1
Riduzione della complessità dei dati
È comune iniziare il processo di registrazione utilizzando immagini con una bassa
complessità, es., immagini che hanno subito un filtraggio di smoothing. Questo
incrementa le chance di successo della registrazione. Consideriamo ora una serie
di immagini ottenute con un uso incrementale del filtro di smoothing. Se le immagini non sono solo filtrate dallo smoothing ma sono anche sotto-campionate, i
dati non hanno solo una complessità minore ma sono effettivamente in quantità
minore. Da ora in poi serie di immagini smooth e sotto-campionate le chiameremo piramidi. Esistono molti tipi di piramidi: piramidi di Gauss, piramidi di
Laplace, piramidi spline, piramidi wavelet ecc. Sicuramente la più comune è 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
CAPITOLO 1. REGISTRAZIONE
1.7.2
Riduzione della complessità delle trasformazioni
La seconda strategia multirisoluzionale è iniziare la registrazione con meno gradi di libertà per il modello della trasformazione. Il grado di libertà di una
trasformazione equivale alla lunghezza (numero di elementi) del vettore µ dei
parametri.
Un esempio di questo si può ottenere applicando trasformazioni di crescente
difficoltà computazionale, se consideriamo una strategia a tre livelli possiamo
ottenere: al primo livello una trasformazione rigida (che produce solo 3 parametri di traslazione in µ), al secondo una trasformazione affine e al terzo una
trasformazione B-Spline.
18
Capitolo 2
La Risonanza magnetica e
metodi di acquisizione delle
immagini
In questo capitolo descriviamo alcuni metodi di acquisizione delle immagini relative alla risonanza magnetica. Citeremo brevemente i fondamenti su cui si basa
la risonanza, i cui principi fisici esulano dallo scopo della nostra trattazione, ed
elencheremo le principali tecniche di acquisizione su cui si basano le immagini
che abbiamo utilizzato e che fanno riferimento alla tecnica di diffusione. Andremo poi a definire le peculiaritá di metodi a singola modalità e dei metodi
multimodali per poi descrivere le misure di similitudine per la registrazione delle
immagini in medicina
2.1
Risonanza magnetica nucleare
La Risonanza Magnetica Nucleare (RMN o, raramente, RNM), in inglese Nuclear Magnetic Resonance (NMR), è una tecnica di indagine sulla materia basata sulla misura della precessione dello spin di protoni o di altri nuclei dotati
di momento magnetico quando sono sottoposti ad un campo magnetico[5]. Le
indagini mediche che sfruttano la RMN son dette anche tomografia a risonanza magnetica e danno informazioni diverse rispetto alle immagini radiologiche
convenzionali: il segnale di densità in RMN è dato infatti dal nucleo atomico dell’elemento esaminato, mentre la densità radiografica è determinata dalle
caratteristiche degli orbitali elettronici degli atomi colpiti dai raggi X.
2.1.1
Tempo di rilassamento longitudinale-T1
Il tempo di rilassamento longitudinale o rilassamento spin-reticolo o spin-lattice
o T1, è una costante di tempo della risonanza magnetica nucleare, che è un
fenomeno della fisica nucleare sfruttato per tecniche d’indagine della materia,
anche in campo biomedico a scopo diagnostico.
Il segnale raccolto al termine dell’irraggiamento dei nuclei con onde elettromagnetiche è costituito da onde aventi la medesima frequenza caratteristica della
19
CAPITOLO 2. LA RISONANZA MAGNETICA E METODI DI
ACQUISIZIONE DELLE IMMAGINI
precessione nucleare. L’andamento del tempo di questo segnale è determinato
da due costanti di tempo chiamate T1 e T2.
Se si perturba un campo magnetico con una magnetizzazione longitudinale,
a 180◦ , i nuclei atomici immersi nel campo invertiranno il proprio vettore di
magnetizzazione, fenomeno che avviene in aumento nel corso del tempo.
In ogni momento l’apparecchio RMN misura i valori della magnetizzazione
residua a un determinato tempo di eco (TE) dal quale, mediante una formula
esponenziale, risulta facile risalire al valore di T1. La scelta del tempo TE è
molto importante, poiché ne viene il contrasto che l’immagine può ottenere. In
genere T1 ha un tempo più lento rispetto a T2 poiché, essendo incrementale,
occorre un tempo sufficientemente lungo per permetterne l’apprezzamento (ma
non troppo lungo da eliminare il contrasto).
2.1.2
Tempo di rilassamento trasversale o rilassamento
spin-spin-T2
Il tempo di rilassamento trasversale o rilassamento spin-spin è una costante di
tempo della risonanza magnetica nucleare, è un fenomeno della fisica nucleare
sfruttato per tecniche d’indagine della materia e in campo biomedico a scopo
diagnostico.
Se si perturba una popolazione di atomi immersi in un campo magnetico,
quindi con spin allineati allo stesso, tramite un impulso a radiofrequenza in modo
da deflettere la direzione dell’orientamento totale degli spin di 90◦ rispetto al
campo magnetico principale, e quindi si azzera l’impulso di eccitazione, i nuclei
atomici immersi nel campo magnetico perderanno gradualmente il sincronismo
di precessione in misura diversa a seconda della disomogeneità della materia
stessa, durante il loro riallineamento al campo magnetico. Questo è correlabile
alla magnetizzazione iniziale e al tempo t, detto tempo “di risonanza” o “di eco”
(TE).
In ogni momento l’apparecchio RMN misura i valori della magnetizzazione
residua a un determinato tempo TE dal quale, mediante una formula logaritmica, risulta facile risalire al valore di T2. La scelta del tempo TE è molto
importante, poiché ne viene il contrasto che l’immagine può ottenere.
Essendo la formula del decadimento logaritmica, si dice comunemente che il
T2 di una sostanza è pari al tempo impiegato a ridurre la sua magnetizzazione
trasversale al 36,79% del suo valore originario (essendo 1/e=0,3679...), ma da
un punto di vista chimico il T2 di una sostanza è determinato dalla libertà di
movimento delle molecole in essa contenute: in acqua pura il T2, determinato
sul rilassamento degli atomi di idrogeno, è massimo (alcuni secondi), mentre
raggiunge valori brevi per i solidi cristallini (nell’ordine dei microsecondi).
2.2
La tecnica di diffusione
La RMN di diffusione [6] è una tecnica sensibile ai movimenti Browniani di
translazione delle molecole d’acqua su piccole distanze (diffusione). Il crescente interesse per tale tecnica è legato al fatto che le immagini di DWI (Diffusion Wieghted Imaging) evidenziano ad esempio le variazioni della mobilità dei protoni dell’acqua indotte dall’ischemia entro pochi minuti dall’attacco
ischemico.
20
CAPITOLO 2. LA RISONANZA MAGNETICA E METODI DI
ACQUISIZIONE DELLE IMMAGINI
Poiché le immagini pesate in diffusione sono molto sensibili al movimento,
sono essenziali sequenze ultraveloci per generare immagini diagnostiche. L’uso di
tecniche ecoplanari ha permesso di acquisire una singola slice in meno di 150 ms,
con la possibilità di eseguire uno studio completo del cervello in circa due minuti.
Per ottenere una sequenza in diffusione è necessario poter applicare dei potenti
gradienti di diffusione. Il primo gradiente pone fuori fase i protoni e viene detto
“dephasing gradient”, il secondo li ripone in fase se non vi è stato movimento
delle molecole dell’acqua. Il principio basilare della misurazione sta nel fatto
che ogni movimento delle molecole dell’acqua nel tempo di osservazione risulta
in una perdita di segnale e, quindi, in una diminuita intensità delle immagini.
In una immagine di diffusione, le strutture con veloce diffusione appariranno
più scure in quanto soggette ad una più alta attenuazione del segnale, mentre
le strutture con una minore velocità di diffusione appariranno più chiare.
2.3
DWI: Diffusion-Wieghted-Imaging
Un’altra tecnica di misurazione della diffusione è l’imaging pesato in diffusione
(Diffusion-weighted imaging, o DWI). Questa tecnica permette la misurazione
della distanza di diffusione delle molecole d’acqua. Più breve è questa distanza,
più chiara appare la regione considerata. In seguito ad una ischemia cerebrale
ad esempio le immagini DWI sono molto sensibili ai cambiamenti patofisiologici
che avvengono nella lesione. Si pensa che l’aumento delle barriere alla diffusione
delle molecole d’acqua come risultato dell’edema citotossico (rigonfiamento delle cellule), sia responsabile dell’incremento del segnale in una scansione DWI.
Altre teorie propongono che l’effetto sia dovuto a cambiamenti nella permeabilità cellulare o al venir meno della flusso citoplasmatico dipendente dall’ATP.
L’aumento del segnale DWI appare entro 5-10 minuti dall’insorgenza dei sintomi dell’attacco ischemico (in contrasto con la tomografia computerizzata, che
normalmente identifica i cambiamenti nei tessuti con un ritardo di 4-6 ore) e
rimane per un periodo fino a due settimane. La TC, per la sua scarsa sensibilità
all’ischemia acuta, è normalmente usata per verificare che non ci sia emorragia,
che impedirebbe l’uso dell’attivatore tissutale plasminogeno (t-PA).
È stato anche proposto che le misure di diffusione tramite MRI potrebbero
essere in grado di identificare cambiamenti molto piccoli nella diffusione dell’acqua extracellulare, il che potrebbe avere applicazioni nel campo della risonanza
magnetica funzionale: il corpo cellulare di un neurone si ingrandisce quando
conduce un potenziale d’azione, impedendo di conseguenza la naturale diffusione delle molecole d’acqua. Dai risultati su modelli animali e dalle informazioni
ottenute da pazienti con ictus (e parzialmente anche nei tumori) tale tecnica ha
le premesse per rivoluzionare la attuale diagnostica per immagini. Questo tipo
di imaging misura la diffusione delle molecole d’acqua nei tessuti biologici. In
un mezzo isotropico (ad esempio in un bicchiere d’acqua), le molecole di liquido
si muovono di moto browniano casuale. Invece nei tessuti biologici la diffusione
può essere anisotropica. DWI è quindi una modifica delle normali le tecniche di
risonanza magnetica, ed è un approccio che utilizza la misura del moto browniano delle molecole. L’ acquisizione RMN utilizza il comportamento dei protoni
in acqua per generare contrasto tra le caratteristiche clinicamente rilevanti di
una parte anatomica particolare. La natura versatile della RMN è dovuta a
questa capacità di produrre contrasto. Pesata in T1, le molecole di acqua in
21
CAPITOLO 2. LA RISONANZA MAGNETICA E METODI DI
ACQUISIZIONE DELLE IMMAGINI
un campione vengono eccitati con l’imposizione di un campo magnetico forte.
Pesata in T2, il contrasto è prodotto misurando la perdita di coerenza o di sincronia tra i protoni dell’acqua. Quando l’acqua è in un ambiente dove si può
liberamente diffondere, il rilassamento tende a richiedere più tempo. In alcune
situazioni cliniche, questo può generare contrasto tra una superficie patologica
ed il tessuto sano circostante.
2.4
DTI: Diffusion-Tensor-Imaging
I tessuti umani riescono a sopravvivere solo entro uno stretto intervallo di temperature in cui la maggior parte dei componenti dei tessuti è allo stato liquido.
Inoltre, i tessuti mostrano a livello microscopico una struttura altamente disomogenea; infatti le membrane cellulari e i vari organelli ostacolano il movimento
libero dell’acqua e di altre molecole. Pertanto la misura della mobilità dell’acqua
può risultare un valido strumento per descrivere la struttura dei tessuti su scala
microscopica, ben oltre la capacità di risoluzione delle usuali metodiche di imaging. L’imaging pesato in diffusione (DWI) e l’imaging del tensore di diffusione
(DTI) sono tecniche di Risonanza Magnetica, sensibili alle proprietà diffusive
delle molecole d’acqua e si presentano come importanti strumenti per la diagnosi
anche nella pratica clinica. Queste metodiche permettono di ottenere immagini
in cui l’intensità del segnale è legata al movimento casuale delle molecole d’acqua, grazie all’inclusione nella sequenza di intensi impulsi di gradiente di campo
magnetico, applicati prima e dopo un impulso a radio frequenza di 180◦ . Il
risultato è una diminuzione del segnale, che può essere ricondotta al coefficiente
di diffusione D. La mappa della variazione dell’intensità del segnale (DWI) può
fornire utili informazioni solo sulla diffusione lungo la direzione lungo la quale è
stato applicato il gradiente di campo magnetico. In molti tessuti, come nella sostanza bianca cerebrale, la diffusione è anisotropa, ovvero la mobilità dell’acqua
dipende dalla direzione, essendo questo un processo tridimensionale, e pertanto
si descrive tramite un tensore, detto tensore di diffusione (D). È possibile determinare i sei elementi indipendenti del tensore di diffusione con l’acquisizione di
almeno sei immagini pesate in diffusione acquisite lungo direzioni non collineari
e di un’immagine di riferimento non pesata in diffusione.
L’imaging di diffusione è un MRI metodo che produce in vivo immagini di
risonanza magnetica di tessuti biologici ponderata con le caratteristiche locali
di diffusione dell’acqua.
2.5
DCE: Diffusion-Contrast-Enhanched
Gli studi DCE (Diffusion Contrast Enhanched) permettono di tracciare nel tempo la distribuzione del mezzo di contrasto paramagnetico nell’organo, al fine di
poter individuare valori di vascolarizzazione e permeabilità caratteristici di tessuti tumorali soggetti ad angiogenesi. Solitamente l’impiego di tale tecnica diagnostica avviene valutando in modo semiquantitativo l’andamento nel tempo del
segnale MR. L’utilizzo clinico della DCE-MRI si concentra sulla localizzazione
di tumori cerebrali e mammellari, infiammazioni, ischemie e demielinizzazioni.
Per quanto riguarda l’utilizzo di questa tecnica per la diagnosi dei tumori
essa consente di:
22
CAPITOLO 2. LA RISONANZA MAGNETICA E METODI DI
ACQUISIZIONE DELLE IMMAGINI
• delineare le zone tumorali (radioterapia)
• stimare il grado del tumore (malignità)
• modulare la radioterapia in base alle mappe di Ktrans
• monitorare il trattamento alla chemio e/o radioterapia
2.6
Metodi a singola modalità verso metodi multimodali
Un’altra classificazione può essere quella tra i metodi a singola modalità e quelli multi modali. I metodi a singola modalità tendono a registrare immagini
nella stessa modalità acquisita dallo stesso tipo di scanner/sensore, mentre i
metodi multimodali di registrazione d’immagini tendono a fare la registrazione
d’immagini acquisite dai diversi tipi di scanner o sensori.
I metodi di registrazione multi-modale usati nell’imaging medico come immagini di un paziente sono spesso ottenuti da diversi scanner. Gli esempi includono
la registrazione delle immagini della TC/MRI cerebrali o dell’intero corpo, come
la PET/CT, a scopi di localizzazione dei tumori, la registrazione d’immagini col
contrasto (TC e RMN) rispetto ad immagini senza contrasto (TC e RMN) per
la segmentazione di parti specifiche dell’anatomia, e la registrazione d’immagini dell’ecografia (ultrasuoni) rispetto a immagini della TC per la localizzazione
della prostata nella radioterapia.
2.7
Misure di similitudine per la registrazione
delle immagini
Il confronto tra immagini simili viene ampiamente usato nell’imaging medico.
Una misura della similitudine delle immagini quantifica il grado di somiglianza
tra i pattern d’intensità in due immagini. La scelta di una particolare modalità
di misurazione delle similtudini dipende dal tipo e modalità d’acquisizione delle
immagini che devono essere sottoposte a registrazione. Alcuni esempi comuni
di somiglianza delle immagini includono la correlazione incrociata, la Mutual
Information vista precedentemente, la somma del quadrato delle differenze di
intensità, e la quota di uniformità dell’immagine.
L’informazione mutua e l’informazione mutua normalizzata sono i metodi
di misura della similarità delle immagini più usati per la registrazione delle
immagini multimodali. La correlazione incrociata, la somma delle differenze
d’intensità elevate al quadrato e le uniformità del radio dell’immagine vengono
comunemente usate per la registrazione d’immagini in mono-modalità.
23
Capitolo 3
Strumenti ed algoritmi
utilizzati
In questo capitolo verranno descritti tool utilizzati e alcuni algoritmi sviluppati.
Il nostro intento è quello di dare una visione chiara dei processi e delle strategie
elaborate nel corso della nostra ricerca.
3.1
Standard e formati di immagini mediche
Iniziamo quindi con il riportare i due formati utilizzati per la rappresentazione
delle immagini biomedicali:
• DICOM: lo standard DICOM[7](Digital Imaging and Communications in
Medicine, immagini e comunicazione digitali in medicina) definisce i criteri per la comunicazione, la visualizzazione, l’archiviazione e la stampa
di informazioni di tipo biomedico quali ad esempio immagini radiologiche. Lo standard DICOM è pubblico, nel senso che la sua definizione è
accessibile a tutti. La sua diffusione si rivela estremamente vantaggiosa
perché consente di avere una solida base di interscambio di informazioni tra apparecchiature di diversi produttori, server e PC, specifica per
l’ambito biomedico. I dati radiologici rappresentabili come immagini o
le immagini vere e proprie che vengono archiviate secondo lo standard
DICOM sotto forma di file vengono comunemente chiamate immagini DICOM. L’errore più comune che viene fatto nell’interpretazione del termine
è che queste siano assimilabili ad altri formati di compressione dell’immagine (es. JPEG, GIF, etc.). In verità lo standard DICOM applicato alla
codifica dei file non è nient’altro che un metodo per incapsulare i dati e
per definire come questi debbano essere codificati o interpretati, ma non
definisce alcun nuovo algoritmo di compressione. La maggior parte delle volte, l’immagine viene archiviata in forma non compressa, secondo la
codifica con la quale viene prodotta, ma esistono molti software che sono
in grado di produrre o interpretare file DICOM contenenti dati compressi
secondo vari algoritmi (JPEG, JPEG Lossless, JPEG Lossy, vari algoritmi
dello standard JPEG2000, ecc.).
25
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
Un file DICOM oltre all’immagine vera e propria,include anche un “Header”. Le informazione contenute nell’Header DICOM sono molteplici,per
esempio: nome e cognome del paziente, il tipo di scansione, posizione e
dimensione dell’immagine ecc. Le informazioni dell’Header vengono scritte in XXX byte del file DICOM,la dimensione ovviamente varia a seconda
della quantità di informazioni memorizzate. Tutte le informazioni memorizzati nell’Header vengono catalogate in gruppi di elementi, detti anche
“Tag DICOM” .
Esistono numerosi programmi che consentono la visualizzazione di immagini DICOM su normali PC, spesso anche liberamente scaricabili via rete
internet. Molti di questi permettono anche di eseguire sulle immagini
elementari operazioni di post-processing, quali misure lineari e di area.
• NII: I file in formato Nifti mantengono le stesse proprietà dei file DICOM
con il vantaggio-svantaggio di avere un unico file per la rappresentazione
del volume; la conversione da un formato all’altro è effettuabile con il
programma 3DSlicer. Noi abbiamo utilizzato due programmi a linea di
comando: DiNifti e DiffUnpack. Esempio di conversione:
diff-unpack API/00010001 dce-1.nii
Querying dicom files... 100% Found 30 single images in the
series.
Volume dimension: 88 128 30 Voxel size: 2.344 2.344 5.000
Writing output files... 100%
Figura 3.1: Tipiche immagini di risonanza magnetica del cervello. Da sinistra
verso destra, è raffigurata rispettivamente la ricostruzione assiale, sagittale e
coronale.
3.2
Il programma 3D Slicer
3D Slicer (Slicer)[8] è un pacchetto di software gratuito, open source per la visualizzazione scientifica e l’analisi delle immagini. Slicer permette di creare “stacks
26
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
di immagini”, partendo da immagini DICOM (Come fa SPM), trasformandole in blocchi di dati tri, quadri e pentadimensionali (ad esempio mostrando la
contemporanea variazione di due parametri, in un reticolo tridimensionale nel
tempo), ad esempio in filmati 3D riguardanti la fMRI, la PET, ecc.
3D Slicer permette anche di confrontare due o più risonanze magnetiche
di uno stesso soggetto nel tempo, per determinare se è apparsa una lesione
(nel individuo sano), se sono apparse nuove lesioni, oppure se queste si sono
ingrandite o ridotte.
3D Slicer viene utilizzato in molte applicazioni medicali, tra queste lo studio
dell’autismo, della sclerosi multipla, del lupus eritematoso sistemico, del carcinoma prostatico, della schizofrenia, nella biomeccanica ortopedica, nella BPCO,
nella malattia cardiovascolare, nella neurochirurgia e nella radioterapia.
Le capacità di Slicer 3D includono:
• Capacità di leggere e di scrivere (per anonimizzare) le immagini DICOM
e un buon numero di altri formati
• Visualizzazione interattiva delle immagini, triangolazione di modelli di
superficie 3D, e volume rendering.
• Editing manuale.
• Fusione e co-“registering” (fusione dei dati di due immagini in tempi
diversi) usando algoritmi di trasformazione rigida e non-rigida
• Segmentazione automatica
• Analisi e visualizzazione dei dati di diffusione del tensore di imaging.
• Tracking di dispositivi per procedure mediche guidate dall’imaging.
Slicer è compilato per l’utilizzo in molteplici piattaforme, includendo Windows, Linux, e Mac OS X.
Viene distribuito sotto una licenza BSD, gratuita e open source. La licenza
di 3D-Slicer non pone restrizioni all’uso del software. Comunque, gli autori
non dichiarano formalmente che il software possa essere utile per alcun compito
particolare. La responsabilità di adempiere alle leggi e alle regolamentazioni
nazionali ricade interamente sull’utente. Slicer non è mai stato approvato per
l’uso clinico in USA o in qualsiasi altro paese.
La piattaforma di Slicer fornisce funzionalità per la segmentazione di immagini e algoritmi di analisi delle immagini. Sono supportati molti formati di
file d’immagine standard, e l’applicazione integra capacità di interfaccia per il
software della ricerca biomedica.
Slicer è stato usato in molteplici ricerche cliniche. Nella ricerca sulla terapia
guidata da immagini, Slicer si utilizza frequentemente per costruire e visualizzare raccolte di dati provenienti dalla MRI che in seguito possono essere disponibili sia per pianificare l’intervento che in ambito intra-operatorio, permettendo
l’acquisizione di coordinate spaziali per il monitoraggio della posizione e del
lavoro degli strumenti operatori. In effetti, Slicer è stato già impiegato nella
“image-guided therapy”, ed attualmente le sue funzioni vengono aumentate e
perfezionate da molti team di ricerca in collaborazione con ingegneri informatici.
In aggiunta alla ricostruzione di modelli 3D dalle immagini RMN convenzionali, Slicer è stato usato anche per presentare informazioni derivate dalla fMRI
27
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
(Imaging di Risonanza Magnetica funzionale) che consiste nell’usare la MRI per
calcolare il flusso sanguigno nelle varie aree del cervello, in relazione più o meno
diretta all’attività neurale encefalica oppure del midollo spinale.
Il campo più avanzato in MRN è quello dello studio del tensore di diffusione “DTI” e “DWI” (usando la RMN per misurare le restrizioni o facilitazioni
date al libero movimento dell’acqua dai fasci di tessuto, permettendo cosı̀ la
ricostruzione di complesse strutture tridimensionali di fasci nervosi, nervi, vasi encefalici, il movimento del liquor nei ventricoli e acquedotti, ecc.) anche
in associazione con programmi di “trattografia” come TrackVis. Le procedure
di diffusione vengono impiegate dall’elettrocardiografia, anche sotto sforzo, per
studiare spasmi e stenosi transienti delle coronarie. Ad esempio, il pacchetto
DTI integrato in 3D-Slicer permette la conversione e l’analisi delle immagini
DTI. I risultati di queste analisi possono essere integrati con i risultati forniti
dalle analisi di MRI morfologica, l’angiografia in RMN e la fMRI.
3.3
Elastix
Elastix[9]è un software open source, basato su ITK[7]. Il software consiste
in una collezione di algoritmi che sono comunemente utilizzati per risolvere i
problemi della registrazione. elastix utilizza un approccio parametrico nella
registrazione, questo significa che impostando le configurazioni di input posso
velocemente ottenere risultati più o meno buoni a seconda delle immagini utilizzate. L’interfacciamento è a linea di comando, questo permette l’utilizzo di
script che semplificano il processamento di elevate quantità di dati. Sempre
nel pacchetto elastix scaricabile all’indirizzo http://elastix.isi.uu.nl/ è presente il programma Transformix che accompagna il processo di registrazione permettendo di effettuare trasformazioni riutilizzando i parametri generati da una
registrazione.
Alcuni esempi di utilizzo:
elastix -f ... -m ... -out out1 -p param1.txt
elastix -f ... -m ... -out out2 -p param2.txt -t0 out1/TransformParameters.0.txt
elastix -f ... -m ... -out out3 -p param3.txt -t0 out2/TransformParameters.0.txt
La chiamata al programma richiede l’impiego di input che sono determinati
dai comandi:
• “-f”: immagine di riferimento;
• “-m”: immagine da mappare;
• “-out”: path della directory di destinazione (i risultati saranno contenuti
in questa cartella);
• “-p”: file di testo contenente i parametri da utilizzare;
• “-t0”: (non obbligatorio) file contenente i parametri generati dalla registrazione, il file viene utile se si utilizza lo strumento trasformix.
I possibili parametri con il comando -p sono già stati riportati durante l’esposizione delle tecniche di registrazione.
Rivediamo questi comandi in comode tabelle:
28
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
Nome
Somma delle Differenze Quadrate (SSD)
Coefficiente di Correlazione Normalizzato (NCC)
Mutua Informazione (MI)
Mutua Informazione Normalizzata (NMI)
Kappa Static (KS)
Comando
AdvancedMeanSquares
AdvancedNormalizedCorrelation
AdvancedMattesMutualInformation
NormalizedMutualInformation
AdvancedKappaStatic
Tabella 3.1: Metriche disponibili in elastix
Nome
Nearest neighbour
Linear
N -th order B-spline
Comando
NearestNeighborInterpolator
LinearInterpolator
BSplineInterpolator o BSplineInterpolatorFloat
Tabella 3.2: Interpolazioni disponibili in elastix
Nome
Translation
Rigid
Similarity
Affine
B-splines
Thin-plate splines
Comando
TranslationTransform
EulerTransform
SimilarityTransform
AffineTransform
BSplinesTransform
SplineKernelTransform
Introduciamo brevemente il programma Transformix presente nel pacchetto
elastix. Il programma permette di utilizzare i parametri di uscita da una
trasformazione. Di seguito riportiamo la sintassi:
transformix -in inputImage -out outputDirectory -tp TransformParameters.txt
dove:
• -in l’immagine di input cioè quella su cui si vogliono applicare i parametri;
• -out cartella dove si salveranno i risultati (diario di computazione, calcoli
effettuati, immagine risultato);
• -tp file dei parametri creati da una registrazione che si vogliono utilizzare
(naturalmente questa procedura può dare anche risultati inutili: si pensi
ad esempio ad una applicazione di trasformazione generata utilizzando A
come immagine fixed e B come immagine moving, nel caso in cui A=B il
risultato sarà una immagine A mappata in modo errato su A stessa. Si
genera quindi un risultato poco interessante e quindi non utilizzabile per
nessuno scopo).
Di seguito è riportato un esempio di file dei parametri utilizzati nelle varie
registrazioni con i relativi commenti:
5
// Vi en e d e f i n i t o i l t i p o d i p i x e l i n t e r n o ,
// u s a t o p e r l a c o m p u t a z i o n e i n t e r n a
// I n g e n e r a l e e ’ s e t t a t o a f l o a t
// NB: q u e s t o non e ’ i l t i p o d e l l ’ immagine !
( FixedInternalImagePixelType ” f l o a t ”)
29
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
( MovingInternalImagePixelType ” f l o a t ”)
10
// La d i m e n s i o n e d e l l ’ immagine f i x e d e d e l l ’ immagine moving
( FixedImageDimension 3 )
( MovingImageDimension 3 )
15
// S p e c i f i c a r e s e s i v u o l e t e n e r c o n t o d e l l a
coseni
// d e l l ’ immagine .
// Raccomandato : t r u e .
( UseDirectionCosines ” true ”)
cosi detta direzione
dei
// ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ Componenti p r i n c i p a l i ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
20
25
30
35
// I s e g u e n t i componenti u s u a l m e n t e s o n o s e t t a t i i n q u e s t o modo :
( R e g i s t r a t i o n ” M u l t i R e s o l u t i o n R e g i s t r a t i o n ”)
( I n t e r p o l a t o r ” L i n e a r I n t e r p o l a t o r ”)
( ResampleInterpolator ” F i n a l L i n e a r I n t e r p o l a t o r ”)
( Resampler ” D e f a u l t R e s a m p l e r ” )
// S e t t a g g i d e i metodi p e r l a r i s o l u z i o n e :
( FixedImagePyramid ” F i x e d R e c u r s i v e I m a g e P y r a m i d ” )
( MovingImagePyramid ” MovingRecursiveImagePyramid ” )
// Componenti s e g u e n t i s o n o molto i m p o r t a n t i :
// L ’ o t t i m i z z a t o r e A d a p t i v e S t o c h a s t i c G r a d i e n t D e s c e n t (ASGD) l a v o r a i n
// g e n e r a l e molto bene . La t r a s f o r m a z i o n e e l a m e t r i c a s o n o i m p o r t a n t i
// e n e c e s s a r i d e f i n i r l i i n m an ie r a a d a t t a p e r l e d i v e r s e a p p l i c a z i o n i
( Optimizer ” AdaptiveStochasticGradientDescent ”)
( Transform ” A f f i n e T r a n s f o r m ” )
( Metric ” NormalizedMutualInformation ”)
// ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ T r a s f o r m a z i o n i ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
40
// S c a l a l a r o t a z i o n e comparandola con l a t r a s l a z i o n e ( devono
// e s s e r e n e l l o s t e s s o i n s i e m e d i v a l o r i )
( AutomaticScalesEstimation ” true ”)
45
// Automaticamente s i suppone c h e l a t r a s l a z i o n e i n i z i a l e e ’ a l l i n e a t a
// con i l c e n t r o g e o m e t r i c o d e l l ’ immagine f i x e d e d e l l ’ immagine moving
( AutomaticTransformInitialization ” true ”)
50
// Le t r a s f o r m a z i o n i s o n o c o m b i n a t e p e r a d d i z i o n e o p e r c o m b i n a z i o n e .
// I n g e n e r a l e , l a c o m p o s i z i o n e e ’ l a m i g l i o r o p z i o n e
// Questo comando non i n f l u e n z a d i molto i l r i s u l t a t o
( HowToCombineTransforms ”Compose ” )
// ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ Misura d i S i m i l a r i t a ’ ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
55
60
// Numero d i l i v e l l i d i g r i g i o i n o g n i l i v e l l o d i r i s o l u z i o n e ,
// p e r l a mutual i n f o r m a t i o n . Lavora bene s e s e t t a t o a 16 o 3 2 .
( NumberOfHistogramBins 3 2 )
// Se usiamo una maschera , q u e s t o e ’ molto i m p o r t a n t e .
// Se l a maschera v i e n e u t i l i z z a t a p e r i m p o s t a r e una r e g i o n e
// d i i n t e r e s s e , s e t t a r e a f a l s e .
( ErodeMask ” f a l s e ” )
// ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ M u l t i r i s o l u z i o n e ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
65
// Numero d i r i s o l u z i o n i . 1 . e ’ s u f f i c i e n t e p e r l i v e l l o d i d e f o r m a z i o n e
// p i c c o l o . 3 o 4 e f f e t t u a una buona r i s o l u z i o n e . Per immagini g r a n d i
// e ad e l e v a t a d e f o r m a z i o n e 5 o 6 e ’ una buona s c e l t a .
( NumberOfResolutions 4)
70
75
// I l f a t t o r e d i s o t t o c a m p i o n a m e n t o / b l u r r i n g p e r l e immagini p i r a m i d a l i .
// Per d e f a u l t , l e immagini s o n o s o t t o c a m p i o n a t e d i un f a t t o r e 2
comparato
// comparato con l a r i s o l u z i o n e s u c e s s i v a .
// Un e s e m p i o d i r i c a m p i o n a m e n t o p i r a m i d a l e d i una immagine 2D e ’ i l
seguente :
/ / ( ImagePyramidSchedule 8 8 4 4 2 2 1 1 )
30
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
( ImagePyramidSchedule 8 8 8
4 4 4
2 2 2
1 1 1 )
// ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ O t t i m i z z a t o r e ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
80
// Massimo numero d i i t e r a z i o n i p e r o g n i l i v e l l o d i r i s o l u z i o n e :
// p e r o t t e n e r e i n una t r a s f o r m a z i o n e r i g i d a un buon l a v o r o s i s e t t a
200 −500.
( MaximumNumberOfIterations 3 0 0 )
85
// P as s o d i o t t i m i z z a z i o n e i n mm. Per d e f a u l t e ’ u s a t a l a d i m e n s i o n e d e l
voxel .
/ / ( MaximumStepLength 1 . 0 )
// ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ Campionamento d e l l ’ immagine ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
90
// Numero d i campionamenti u s a t i n e l l a c o m p u t a z i o n e d e l l a
// mutual i n f o r m a z i o n . Con A d a p t i v e S t o c h a s t i c G r a d i e n t D e s c e n t 2000
// s o n o s u f f i c i e n t i .
( NumberOfSpatialSamples 2 0 4 8 )
95
// R i p u l i s c e q u e s t i campionamenti i n o g n i i t e r a z i o n e , e
// i n modo random .
( NewSamplesEveryIteration ” true ”)
( ImageSampler ”Random ” )
100
105
110
li
seleziona
// ∗∗∗∗∗∗∗∗∗∗∗∗∗ I n t e r p o l a z i o n e e r i c a m p i o n a m e n t o ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
// Ordine d e l l ’ i n t e r p o l a z i o n e B−S p l i n e u s a t a d u r a n t e l a r e g i s t r a z i o n e /
ottimizzazione .
// Puo ’ aumentare l a a c c u r a t e z z a d e l r i s u l t a t o s e s e t t a t a a 3 .
// E ’ s c o n s i g l i a t o s e t t a r l o a z e r o . Mentre s e t t a t o ad uno e f f e t t u a una
interpolazione
// l i n e a r e . I n m o l t i c a s i a p p l i c a t i q u e s t a e ’ una buona s c e l t a .
( BSplineInterpolationOrder 3)
// Ordine d i i n t e r p o l a z i o n e B−S p l i n e u s a t a d u r a n t e l a d e f o r m a z i o n e
finale .
// 3 da ’ una buona a c c u r a t e z z a ; e ’ raccomandata i n m o l t i c a s i .
// 1 da ’ una a c c u r a t e z z a p e g g i o r e d e l l a p r e c e d e n t e ( e f f e t t u a una
interpolazione lineare )
// 0 da ’ una a c c u r a t e z z a p e g g i o r e , ma e ’ a p p r o p r i a t a p e r immagini
binarie .
( FinalBSpline Interpolat ionOrder 3)
115
// V a l o r e d e i p i x e l e s t e r n i
( DefaultPixelValue 0)
d e l l ’ immagine s e t t a t o p e r d e f a u l t a 0 .
120
// S e t t a r e q u e s t o comando quando s i v u o l e g e n e r a r e l ’ immagine moving
deformata .
// S i puo ’ s e t t a r e ad e s e m p i o a f a l s o s e s i e ’ i n t e r e s s a t i
// s o l a m e n t e a l l a d e f o r m a z i o n e f i n a l e non r i g i d a .
( WriteResultImage ” true ”)
// Tipo d i p i x e l d e l f o r m a t o d e l l ’ immagine d e l l a
effettuata .
( ResultImagePixelType ” short ”)
( ResultImageFormat ” n i i ” )
registrazione
I parametri utilizzati sono di chiaro significato. Il programma prende in input
l’immagine e i parametri di trasformazione e produce in output il risultato all’interno di una directory. Naturalmente sta all’utente accertarsi di effettuare una
operazione di trasformazione “sensata”: questo significa che i parametri generati da una registrazione potranno essere utilizzati solo in un contensto coerente.
Transformix verrà utilizzato in più di una occasione nella nostra trattazione,
per affrontare alcune problematiche che vedremo in seguito.
31
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
3.4
Fsl
La FMRIB Software Library (con acronimo: FSL), che si può tradurre in italiano “libreria software per la risonanza magnetica funzionale”, consiste in una serie
di programmi che contengono strumenti per l’analisi delle immagini e programmi
di statistica per l’analisi dei dati provenienti da immagini di risonanza magnetica
funzionale, strutturale e di risonanza magnetica a tensore di diffusione.
I programmi FSL sono disponibili sia come file binario pre-compilato che
come codice sorgente disponibile sia per Mac OS X che per i PC. Vengono
concessi gratuitamente per l’utilizzo non commerciale.
Di seguito i programmi utilizzati ed una breve descrizione:
• BET -brain extraction-: (strumento per l’estrazione cerebrale) segmenta
(separa) il cervello dal non-cervello generando dati in forma di maschere
binarie, fornisce una modellazione delle superfici del cranio e dello scalpo.
• FLIRT: è uno strumento di registrazione delle immagini lineare.
• FNIRT: è uno strumento di registrazione delle immagini non lineare.
• FSLMATHS: è uno strumento per le operazioni matematiche da effettuare
sulle immagini, ammette diversi comandi come -add (addizione) o -mul
(moltiplicazione) che rappresentano gli operatori. Nel progetto è stato
utilizzato per la moltiplicazione di una registrazione con una maschera
allo scopo di eliminare rumore.
• FSLView: è uno strumento per la visualizzazione delle immagini nei formati Nifti Dicom ecc. È stato utilizzato anche per la costruzione di modelli
tridimensionali dei volumi usati.
Gli strumenti più utilizzati sono “Flirt” e “Fnirt”. “Flirt” è il programma
principale per effettuare la trasformazione affine (o in generale quella lineare).
Le opzioni principali sono: un volume di input (-in) e uno di riferimento (-ref );
Il calcolo della trasformazione affine che registra l’immagine di input sull’immagine di riferimento produce in output con il comando (-omat) la matrice affine e
con il comando(-out) l’immagine di input trasformata. Altri comandi che abbiamo utilizzato sono stati (-cost) che specifica la metrica da utilizzare (-interp)
per specificare l’interpolazione (-forcescaling) per ricampionare l’immagine di
input su quella di riferimento. Il programma “Fnirt” si utilizza invece per la registrazione non lineare. A seguire un esempio dell’utilizzo di questo programma
(limitatamente al nostro utilizzo):
fnirt -v –in=“dwi.nii” –ref=“dce.nii” –splineorder=2 –iout=“nonlineare.nii” –cout=“non-lineare-warp.nii”
I comandi -in e -ref sono rispettivamente per l’immagine moved e per l’immagine fixed; il comando -splineorder serve per impostare il grado della trasformazione non lineare. Mentre i comandi di output sono -iout per la produzione
dell’immagine risultato e -cout per la produzione del warp (analogamente alla matrice affine prodotta da flirt). Per concludere il comando -v produce in
output durante la fase di registrazione tutte le informazioni sviluppate e tutti i
passi eseguiti in tempo reale.
32
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
3.5
Filtro Anisotropico
Il filtro anisotropico (in inglese Anisotropic filtering, abbreviato quindi in AF) è
un metodo utilizzato nella grafica computerizzata per accrescere la qualità delle
immagini che ritraggono texture poste su superfici inclinate rispetto al punto di
osservazione (figura 3.2). Come il filtro bilineare e trilineare esso riduce l’effetto
di aliasing, inoltre causa meno sfuocamento in corrispondenza degli angoli di
visuale più estremi conservando maggiore dettaglio.
Figura 3.2: Esempio di applicazione di un filtro anisotropico
L’uso del filtro anisotropico è abbastanza costoso in termini di risorse ed
è diventato una caratteristica standard delle schede video a disposizione dei
normali consumatori negli anni novanta.
Questo filtro è ormai diventato comune nelle schede grafiche moderne, attivabile sia dall’utente attraverso le impostazioni dei driver, sia dalle applicazioni
grafiche o dai videogiochi.
3.5.1
Algoritmi di Perona e Malik
L’idea base di questo algoritmo è semplice[11]. Se si effettua la convoluzione tra
l’immagine e un kernel gaussiano si ottiene l’immagine di partenza filtrata con un
passa basso. Se si itera questa operazione sulle immagini che man mano si vanno
generando si ottiene una serie di immagini ognuna delle quali ha tre parametri:
x e y che sono le dimensioni, e t che indica il numero di iterazioni per ottenere
l’immagine in questione. Quindi l’immagine I(x, y, 0) corrisponde all’originale
mentre le immagini I(x, y, t) sono quelle ottenute mediante la convoluzione;
all’aumentare di t le immagini sono sempre più sfocate.
Questo procedimento è quello che viene chiamato filtraggio lineare, perchè è
indipendente dalle caratteristiche locali dell’immagine.
Inoltre questa convoluzione corrisponde, come dice Koenderink, alla soluzione della funzione di diffusione:
It = ∆I = (Ix x + Iy y)
Inoltre Koenderink [12] giustifica l’utilizzo della funzione di diffusione definendo
due concetti che sono concordanti con la sua scelta. Il primo è il criterio di
causalità (che viene rispettato anche da Perona e Malik). Il secondo è il criterio
33
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
di omogeneità, per il quale il “blurring” deve essere spazio-invariante. Questo
criterio viene contestato da Perona e Malik, che lo considerano frutto di una
semplificazione. È anzi su questo punto che sta la differenza tra filtraggio lineare
e quello non lineare.
Anzichè eseguire un “blurring” costante, questo viene modulato in funzione
delle caratteristiche locali dell’immagine, che vengono rappresentate mediante
il coefficiente di diffusione c(x, y, t). La funzione di diffusione diventa quindi:
It = ∇(c(x, y, t)∇I)
Per fare in modo che il filtro agisca maggiormente all’interno delle regioni e non
nei contorni, la quantità c deve essere prossima a zero nei contorni e prossima ad
uno nelle zone da filtrare. Una quantità che può rappresentare bene un contorno
è il gradiente: maggiore è il modulo del gradiente nel punto in questione e
maggiore sarà la probabilità che si tratti di un contorno. La quantità c può
essere scelta quindi in funzione del gradiente:
c = g(k ∇I(x, y, t) k)
le funzioni proposte da Perona (3.1) e Malik (3.2), per la corrispondenza tra il
modulo del gradiente e la quantità c sono le seguenti:
g(∇I) = e−(
k∇Ik
K
2
)
(3.1)
1
g(∇I) =
1+
k∇Ik
K
2
(3.2)
Il parametro k deve essere calcolato in relazione alle caratteristiche dell’imma-
Figura 3.3: Funzioni g(∗) proposte da Perona e Malik
gine considerata, e regola la pendenza delle due funzioni.
34
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
Come si può vedere in figura 3.3, queste funzioni hanno un valore massimo pari ad uno, in corrispondenza del gradiente nullo. All’aumentare del
gradiente le funzioni tendono a zero, in modo che per valori alti di gradiente,
corrispondenti ad un contorno, il coefficiente di diffusione sia circa nullo.
Per capire meglio il funzionamento dell’algoritmo è conveniente definire una
funzione di flusso:
φ(∇I) = g(∇I) · ∇I
La funzione di diffusione può essere quindi vista come:
It = ∇ · (φ(∇I))
Il flusso è massimo per valori di gradiente prossimi a k (figura 3.4). Quindi
se si vuole eliminare il rumore dell’immagine bisogna scegliere un valore di k
prossimo al gradiente del rumore. La certezza che applicando questo filtro viene
Figura 3.4: Funzioni di flusso
rispettato il criterio di causalità viene data dal teorema del massimo principio
citato da Perona e Malik.
La funzione di diffusione può essere discretizzata, e resa molto semplice e
facilmente implementabile, nella seguente maniera:
∂
∂
∂
∂
∂
I(x, y, t) =
c(x, y, t) I(x, y, t) +
c(x, y, t) I(x, y, t) ≈
∂x
∂x
∂x
∂y
∂y
1
∆x
≈
, y, t (I(x + ∆x, y, t) − I(x, y, t)) +
c x+
(∆x)2
2
1
∆x
c
x
−
,
y,
t
(I(x,
y,
t)
−
I(x
−
∆x,
y,
t))
+
−
(∆x)2
2
35
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
1
∆y
+
c x, y +
, t (I(x, y + ∆y, t) − I(x, y, t)) +
(∆y)2
2
∆y
1
c
x,
y
−
,
t
(I(x,
y,
t)
−
I(x,
y
−
∆y,
t))
=
−
(∆y)2
2
= φeast + φwest + φnorth + φsouth
I(x, y, t + ∆t) ≈ I(x, y, t) + ∆t · (φeast + φwest + φnorth + φsouth )
In cui i valori di φ possono essere approssimati, per ogni punto con:
φeast = g(|∇east I|) · (∇east I)
φwest = g(|∇west I|) · (∇west I)
φnorth = g(|∇north I|) · (∇north I)
φsouth = g(|∇south I|) · (∇south I)
In cui:
∇east I(i,j) = I(i,j+1) − I(i,j)
∇west I(i,j) = I(i,j−1) − I(i,j)
∇north I(i,j) = I(i−1,j) − I(i,j)
∇south I(i,j) = I(i+1,j) − I(i,j)
Mentre il termine ∆t viene solitamente posto a 1/4.
3.5.2
Codice
Lo script “anisodiff3d.m” applica il filtro anisotropico, secondo le formule di
Perona (3.1) e Malik (3.2), prendendo in ingresso:
• VOL: immagine 3D in scala di grigi;
• NUM-ITER: numero di iterazioni;
• DELTA-T: integrazione costante (0 <= ∆t <= 3/44) (di solito impostato
sempre al valore massimo di 3/44);
• KAPPA: valore di soglia che modifica la rigidità del metodo;
• OPTION: se impostato ad 1 applica l’algoritmo di Perona (3.1) altrimenti
se impostato a 2 applica l’algoritmo di Malik (3.2).
• VOXEL-SPACING: vettore 3x1 che rappresenta la grandezza di un voxel
(in millimetri) nell’immagine 3d (x, y, z), se l’immagine è di tipo DICOM
l’informazione può essere ricavata dal file dicominfo.
Per l’analisi completa dello script sviluppato si rimanda all’appendice. Di
seguito riportiamo invece un esempio di utilizzo:
36
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
function
5
10
5
10
f i l t r a t a = f i l t r a ( immagine )
img=immagine . img ;
d e l t a t = 3/44;
num iter = 2;
kappa = 8 0 ;
option = 2;
v o x e l s p a c i n g = ones ( 3 , 1 ) ;
p r i m o s t e p = a n i s o d i f f 3 D ( img , n u m i t e r ,
voxel spacing ) ;
immagine . img=p r i m o s t e p ;
f i l t r a t a =immagine ;
return
d e l t a t , kappa , o p t i o n ,
%Lo s c r i p t c a r i c a l ’ immagine , l o f i l t r a , c r e a l a maschera con r e g i o n
growing e s a l v a
clear all
d c e=l o a d n i i ( ’ c o n v e r t i t i / d c e _ 1 3 . n i i ’ ) ;
r i s u l= f i l t r a ( dce ) ;
r i s u l t a t o=c r e a m a s c h e r a ( r i s u l . img , 6 ) ;
d c e . img=r i s u l t a t o ;
clear risultato
r i s u l t a t o=d c e ;
s a v e n i i ( r i s u l t a t o , ’ maschera_13aniso . nii ’ ) ;
%Questo v i e n e r i p e t u t o p e r o g n i immagine da f i l t r a r e
...
3.6
Region Growing
La tecnica di “region growing” si basa sul raggruppamento automatico di regioni
di pixel che hanno una proprietà in comune. Ad esempio gruppi di pixel che
hanno la medesima luminosità. Suddividiamo l’immagine in cluster ossia in
zone con pixel aventi luminosità uguale. Si definisce poi una “distanza” ossia
una misura fra due cluster diversi. Se tale misura diventa inferiore ad un valore
prefissato allora i due cluster presi in esame si uniscono in un nuovo cluster
altrimenti si continua la ricerca fino a che non si finiscono tutte le coppie di
cluster. La ricerca parte da un punto preso nella regione di interesse.
Nel nostro progetto questa tecnica è stata utilizzata per lo sviluppo di alcune
maschere 3D.
Di seguito riportiamo il codice per l’applicazione bidimensionale:
5
10
f u n c t i o n [ g , NR, SI , TI ] = r e g i o n g r o w ( f , S , T)
%REGIONGROW Perform s e g m e n t a t i o n by r e g i o n g r o w i n g .
%
[ G, NR, SI , TI ] = REGIONGROW( F , SR , T) .
S can be an a r r a y ( t h e
%
same s i z e a s F) w i t h a 1 a t t h e c o o r d i n a t e s o f e v e r y s e e d p o i n t
%
and 0 s e l s e w h e r e .
S can a l s o be a s i n g l e s e e d v a l u e . S i m i l a r l y ,
%
T can be an a r r a y ( t h e same s i z e a s F) c o n t a i n i n g a t h r e s h o l d
%
v a l u e f o r e a c h p i x e l i n F . T can a l s o be a s c a l a r , i n which
%
c a s e i t becomes a g l o b a l t h r e s h o l d .
%
%
On t h e output , G i s t h e r e s u l t o f r e g i o n growing , w i t h e a c h
%
r e g i o n l a b e l e d by a d i f f e r e n t i n t e g e r , NR i s t h e number o f
%
r e g i o n s , S I i s t h e f i n a l s e e d image u s e d by t h e a l g o r i t h m , and TI
%
i s t h e image c o n s i s t i n g o f t h e p i x e l s i n F t h a t s a t i s f i e d t h e
%
threshold test .
15
20
f = double ( f ) ;
% I f S i s a s c a l a r , o b t a i n t h e s e e d image .
i f numel ( S ) == 1
S I = f == S ;
S1 = S ;
else
% S i s an a r r a y . E l i m i n a t e d u p l i c a t e , c o n n e c t e d s e e d l o c a t i o n s
% t o r e d u c e t h e number o f l o o p e x e c u t i o n s i n t h e f o l l o w i n g
% s e c t i o n s o f code .
37
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
S I = bwmorph ( S , ’ s h r i n k ’ , I n f ) ;
J = f i n d ( SI ) ;
S1 = f ( J ) ; % Array o f s e e d v a l u e s .
25
end
35
TI = f a l s e ( s i z e ( f ) ) ;
f o r K = 1 : l e n g t h ( S1 )
s e e d v a l u e = S1 (K) ;
S = a bs ( f − s e e d v a l u e ) <= T ;
TI = TI | S ;
end
40
% Use f u n c t i o n i m r e c o n s t r u c t w i t h S I a s t h e marker image t o
% obtain the r e g i o n s c o r r e s p o n d i n g to each seed i n S . Function
% bwlabel a s s i g n s a d i f f e r e n t i n t e g e r to each connected r e g i o n .
[ g , NR] = b w l a b e l ( i m r e c o n s t r u c t ( SI , TI ) ) ;
30
A seguire lo script ottimizzato per il nostro utilizzo (con estensione al 3D):
5
10
15
f u n c t i o n r i s u l t a t o=c r e a m a s c h e r a ( immagine , t )
%Crea Maschera
%i n p u t :
%immagine=immagine 3d
%t l i v e l l o d i t h r e s h o l d
%
%o u t p u t :
%r i s u l t a t o=maschera 3d
clear risultato ;
[ s i z e 1 s i z e 2 s i z e 3 ] = s i z e ( immagine ) ;
r i s u l t a t o=z e r o s ( s i z e 1 , s i z e 2 , s i z e 3 ) ;
maschera=z e r o s ( s i z e 1 , s i z e 2 ) ;
maschera ( c e i l ( s i z e 1 / 2 ) , c e i l ( 3 ∗ s i z e 2 / 4 ) ) =1;
f o r k =1:1: s i z e 3
r i s u l t a t o ( : , : , k ) = r e g i o n g r o w ( immagine ( : , : , k ) , maschera , t ) ;
end
return
3.7
Strumento per l’analisi oggettiva dei dati
Per la valutazione dei risultati ottenuti abbiamo utilizzato una funzione chiamata Tester scritta in codice C++ e basata sulle librerie di ITK (Insight Segmentation and Registration Toolkit) che permette di analizzare la registrazione
effettuata in base a tre parametri significativi i cui range vanno da 0 a 1;
• Errore calcolo NMI: Errore che indica il valore della Mutual Information.
• Errore calcolo RMSDifferences: Errore che indica la differenza di luminosità relative delle immagini registrate.
• Errore calcolo EdgeOverlap : Errore che indica il valore di sovrapposizione
degli edge delle due immagini registrate.
L’esecuzione del programma è computazionalmente costosa e il risultato va interpretato a seconda delle immagini a registrate. Ad esempio nella nostra elaborazione teniamo conto principalmente del terzo valore restituito, quello che
si riferisce alla sovrapposizione degli edge che, nel nostro caso rappresenta il
parametro più significativo in quanto la natura stessa delle immagini utilizzate è sensibilimente diversa (i.e. DWI/DCE) Nel caso del valore relativo alla
MutualInformation avremo un risultato positivo in cui esso si avvicinerà ad 1
mentre nel caso degl’altri due parametri calcolati essi risulteranno positivi se
tenderanno a 0. Per normalizzare i risultati e avere anche una visualizzazione
38
CAPITOLO 3. STRUMENTI ED ALGORITMI UTILIZZATI
grafica dell’analisi dei dati (questo argomento sarà trattato approfonditamente
nel capitolo 4) abbiamo creato una funzione (Matlab) che di seguito riportiamo:
%E r r o r e Mutual I n f o r m a t i o n
MI = [ ] ;
5
10
15
20
%E r r o r e D i f f e r e n z a
DI = [ ] ;
Intensita
%E r r o r e EdgeOverlap
EO= [ ] ;
%N o r m a l i z z a z i o n e
matrix n or m=z e r o s ( s i z e (MI) , 3 ) ;
m a t r i x=c a t ( 3 , MI , DI ,EO) ;
matrix n or m ( : , 1 )=normc(1− m a t r i x ( : , 1 ) ) ;
matrix n or m ( : , 2 )=normc ( m a t r i x ( : , 2 ) ) ;
matrix n or m ( : , 3 )=normc ( m a t r i x ( : , 3 ) ) ;
%G r a f i c i
x=[1:36];
figure
p l o t ( x , m a tr i x n or m ( : , 1 ) , ’ - g ’ , ’ L i n e W i d t h ’ , 2 )
h o l d on
p l o t ( x , m a tr i x n or m ( : , 2 ) , ’ - r ’ , ’ L i n e W i d t h ’ , 2 )
25
p l o t ( x , m a tr i x n or m ( : , 3 ) , ’ - b ’ , ’ L i n e W i d t h ’ , 2 )
30
g r i d on
%l e g e n d ( ’ MI ’ , ’NRMS’ , ’EDGE’ )
l e g e n d ( ’ MI ’ , ’ NRMS ’ , ’ E d g e O v e r l a p ’ )
yl a be l ( ’ Error ’ )
x l a b e l ( ’ Registration Case ’ )
39
Capitolo 4
Realizzazione e sviluppo
In questo capitolo tratteremo nel dettaglio tutte le analisi e le operazioni svolte
nella nostra ricerca. Descriveremo come a partire dalle immagini di Risonanza
Magnetica di diversa acquisizione e natura arriveremo a creare una buona registrazione attraverso gli algoritmi e i parametri messi a disposizione dai tool visti
in precedenza e che tratteremo separatamente solo dopo aver effettuato un approfondimento sui metodi di “raffinamento” che abbiamo seguito (sia per FSL
e elastix) ai fini di evitare la creazione di una mole di dati troppo elevata e di
difficile analisi. Le operazioni di svolgimento, i dati ottenuti e le le registrazioni
effettuate sia relative alle immagini del seno e a quelle di testa-collo, verranno
qui di seguito elaborati e rappresentati.
4.1
Immagini del seno
Le immagini del seno dei 5 soggetti analizzati sono organizzate in un dataset,
forniteci dalla Fondazione IEO [13], contenenti le acquisizioni di tipo DWI all’interno della relativa cartella (breast-diff-weigthed) e 8 serie di immagini DCE
acquisite a distanza di 60 secondi l’une dalle altre, dopo l’iniezione del liquido
di contrasto.
• dce-13 : Acquisizione contemporanea all’inserimento del liquido di contrasto
• dce-14 : Acquisizione dopo 1 minuto dall’inserimento del liquido
• dce-19 : Acquisizione dopo 2 minuti dall’inserimento del liquido
• dce-21 : Acquisizione dopo 3 minuti dall’inserimento del liquido (volume
di riferimento)
• dce-23 : Acquisizione dopo 4 minuti dall’inserimento del liquido
• dce-25 : Acquisizione dopo 5 minuti dall’inserimento del liquido
• dce-27 : Acquisizione dopo 6 minuti dall’inserimento del liquido
• dce-29 : Acquisizione dopo 7 minuti dall’inserimento del liquido
41
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Il contrasto maggiore che permette una migliore visualizzazione delle patologie (in questo caso un nodulo al seno), avviene a distanza di 2-3 minuti
dall’inserimento del liquido (dce-13), quindi abbiamo utilizzato come immagini
di riferimento (o campione rappresentativo) quelle ottenute dopo tale intervallo di tempo (nel nostro dataset sono indicate con l’identificativo dce-21). La
visualizzazione della patologia in modo chiaro e ben definito ci ha permesso
di utilizzare quest’ultima come parametro di confronto con le registrazioni che
verranno effettuate, in quanto rappresenta una variante anatomica sostanziale.
Le immagini forniteci sono nel formato DICOM descritto nel capitolo 3 e per
essere utilizzate dai tool sono state convertite nel formato nii attraverso, come
abbiamo visto, il comando dinifti oppure attraverso diff-unpack.
Figura 4.1: Immagine moving di tipo dwi, a destra si evidenzia la lesione oggetto
di idagine
42
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.2: Immagine di riferimento di tipo dce-21
4.1.1
Traccia di Sviluppo
Lo sviluppo è stato strutturato con l’obiettivo di ricercare la registrazione migliore per i volumi in nostro possesso utilizzando come immagine di riferimento
la DCE e come immagine di input la DWI. Tale procedimento può essere inteso
e pensato come registrazione multimodale per l’aumento della capacità informativa dell’immagine DWI[15]. Per quanto riguarda elastix possiamo scomporre
la nostra ricerca in tre fasi. In primis ci siamo occupati di effettuare un set di
registrazioni utilizzando i parametri di maggior rilevanza (“brute force”). Una
volta ottenuto il risultato migliore abbiamo continuato con il settaggio dei parametri meno rilevanti (che hanno quindi un impatto minore sulla registrazione)
in modo da migliorare ancor più la qualità della registrazione (procedimento
“top down”). Come passo finale al risultato ottenuto abbiamo applicato diversi
procedimenti di alto livello per aumentare al massimo il livello di qualità della
registrazione. Per la valutazione ci siamo concentrati in un primo tempo sul
risultato visivo della registrazione confrontando direttamente l’output ottenuto
con la DCE e poi con l’ausilio del tester già presentato all’inizio (valutazione
oggettiva). Per quanto riguarda FSL il procedimento è stato effettuato in modo
simile. Sono state sviluppate diverse registrazioni con diversi parametri in modo
da ottenere un criterio di registrazione per questo tipo di immagini volumetriche. Procedendo con ordine di seguito in tabella troviamo i parametri utilizzati
nella prima fase di registrazione con elastix:
43
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
N◦
1
2
2
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
Trasformazione
Rigida
Rigida
Rigida
Rigida
Rigida
Rigida
Rigida
Rigida
Rigida
Rigida
Rigida
Rigida
Affine
Affine
Affine
Affine
Affine
Affine
Affine
Affine
Affine
Affine
Affine
Affine
BSpline
BSpline
BSpline
BSpline
BSpline
BSpline
BSpline
BSpline
BSpline
BSpline
BSpline
BSpline
Interpolazione
Nearest neighbour
Nearest neighbour
Nearest neighbour
Nearest neighbour
Linear
Linear
Linear
Linear
BSpline
BSpline
BSpline
BSpline
Nearest neighbour
Nearest neighbour
Nearest neighbour
Nearest neighbour
Linear
Linear
Linear
Linear
BSpline
BSpline
BSpline
BSpline
Nearest neighbour
Nearest neighbour
Nearest neighbour
Nearest neighbour
Linear
Linear
Linear
Linear
BSpline
BSpline
BSpline
BSpline
Campionamento
Random
Random
RandomCoordinate
RandomCoordinate
Random
Random
RandomCoordinate
RandomCoordinate
Random
Random
RandomCoordinate
RandomCoordinate
Random
Random
RandomCoordinate
RandomCoordinate
Random
Random
RandomCoordinate
RandomCoordinate
Random
Random
RandomCoordinate
RandomCoordinate
Random
Random
RandomCoordinate
RandomCoordinate
Random
Random
RandomCoordinate
RandomCoordinate
Random
Random
RandomCoordinate
RandomCoordinate
Tabella 4.1: Tabella generale dei parametri
44
Metrica
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
MI
NMI
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
N◦
1
2
3
4
5
6
7
8
9
Metrica
Somma delle Differenze Quadrate (SSD)
Somma delle Differenze Quadrate (SSD)
Somma delle Differenze Quadrate (SSD)
Mutua Informazione Normalizzata (NMI)
Mutua Informazione Normalizzata (NMI)
Mutua Informazione Normalizzata (NMI)
Coefficiente di Correlazione Normalizzato (NCC)
Coefficiente di Correlazione Normalizzato (NCC)
Coefficiente di Correlazione Normalizzato (NCC)
Numero iterazioni
2
4
3
2
4
3
2
4
3
Tabella 4.2: Raffinazioni successive in elastix
Dove MI sta per “Mutual Information” mentre NMI per “Normalized Mutual
Information”. È facile notare che l’idea di fondo è la prova di alcuni ingressi
(parametri di input al tool elastix) provando tutte le combinazioni possibili.
Come vedremo nella sezione Risultati gli output migliori prodotti da queste
registrazioni sono stati ottenuti con la trasformazione Affine, in particolare la
registrazione che ha ottenuto la miglior valutazione è la numero 18. L’influenza
maggiore nella scelta dei parametri è data dalla scelta del tipo di trasformazione.
Come passo ulteriore abbiamo quindi riutilizzato il set di parametri utilizzati
nel tentativo migliore modificando di volta in volta alcuni parametri.
Le modifiche riportate partendo dalla configurazione numero 18 sono visibili
in tabella 4.2.
Questo ultimo tentativo non ci ha portato ad un miglioramento del risultato ma ha confermato la configurazione 18 (naturalmente alzando il numero di
iterazioni per l’ottimizzazione otteniamo un miglioramento del risultato). Aumentando il numero di iterazioni per l’ottimizzazione come abbiamo effettuato
in questi test ci siamo accorti che le performance peggiorano di un fattore non
significativo. Un fattore invece importante per ottenere un rapporto risultatoprestazioni efficiente è quello di utilizzare un campionamento Random. Come
già detto durante le premesse impostare un campionamento Full durante il processo di registrazione produce un risultato migliore ma peggiora drasticamente i
tempi di calcolo. In tabella 4.3 riportiamo la combinazione dei parametri utilizzata per la prova effettuata con FSL. Come passo successivo abbiamo applicato
la trasformazione 18 a tutte le fasi di propagazione del liquido di contrasto
(a tutte le immagini DCE da noi in possesso) verificando cosı̀ se veramente i
risultati ottenuti possono essere ritenuti accettabili per questo tipo di immagini.
Una prima considerazione finale sulla registrazione ottenuta ci ha fatto notare che in alcune parti del volume risultato è presente un artefatto dell’immagine
su un lato (figura 4.3).
L’errore si è propagato dalla immagine DWI. Questa distorsione è probabilmente dovuta ad un movimento del soggetto durante il processo di acquisizione
della immagine DWI. L’obiettivo che ci siamo prefissati in questa fase è stato
quello di vedere quanto ha influito sul processo di registrazione questa distorsione. Per fare ciò abbiamo ideato un meccanismo che ci ha permesso di ottenere la
trasformazione da applicare alla DWI senza registrare direttamente l’immagine
DWI con rumore. La nostra idea è stata quella di creare un “atlante” sviluppato
sfruttando la simmetria del corpo umano. Abbiamo utilizzato la parte sinistra
45
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
N◦
1
2
3
4
5
6
7
8
9
10
11
12
Interpolazione
Nearest Neighbour
Nearest Neighbour
Nearest Neighbour
Nearest Neighbour
Lineare
Lineare
Lineare
Lineare
Trilineare
Trilineare
Trilineare
Trilineare
Metrica
Mutulal Information
Correlation Ratio
Normalised Correlation
Normalised Mutual Information
Mutulal Information
Correlation Ratio
Normalised Correlation
Normalised Mutual Information
Mutulal Information
Correlation Ratio
Normalised Correlation
Normalised Mutual Information
Tabella 4.3: Parametri utilizzati in FSL
Figura 4.3: Nella parte evidenziata troviamo l’errore propagato dalla DWI.
del volume DWI, che a differenza della parte destra non è stata rovinata nel
processo di acquisizione, per ricostruire simmetricamente la parte danneggiata.
In altre parole ci siamo creati una nuova immagine DWI replicando simmetricamente la parte sinistra non danneggiata. Questa nuova immagine, che da ora in
poi chiameremo immagine “simmetrica”, ci è servita per avere una riproduzione
ideale di una immagine DWI “perfetta” (per perfetta intendiamo naturalmente
una immagine senza deformazioni o rumore). L’immagine simmetrica è stata
poi registrata sull’immagine DCE. La trasformazione ottenuta da quest’ultima
registrazione è stata applicata poi all’immagine DWI. In figura 4.4 possiamo
vedere questo primo processo di registrazione.
In questa registrazione la trasformazione generata serve per mappare l’immagine simmetrica sull’immagine DCE. Quindi tale trasformazione non può ancora
essere riutilizzata per mappare effettivamente l’immagine DWI sulla DCE. Ecco
perchè, come è possibile vedere in figura 4.5, abbiamo introdotto un’ ulteriore
46
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.4: Registrazione dell’immagine simmetrica su DCE e successivo
riutilizzo della trasformazione generata su immagine DWI.
registrazione. Ora viene effettuata una prima registrazione utilizzando l’immagine simmetrica come immagine moved e l’immagine DWI come immagine fixed
riportando cosı̀ l’immagine simmetrica in un sistema di riferimento uguale alla
DWI. Il risultato di questa prima registrazione viene riutilizzato come immagine moved nella seconda registrazione con immagine fixed DCE. In altre parole
si effettua una nuova registrazione fra il risultato della prima registrazione e
l’immagine DCE generando cosı̀ la trasformazione necessaria per passare dalla
DWI alla DCE.
Figura 4.5: Registrazione dell’immagine simmetrica su DCE (utilizzando come
passo intermedio di registrazione l’immagine DWI) e successivo riutilizzo della
trasformazione generata su immagine DWI.
Ricordiamo che lo scopo di utilizzare questo meccanismo è quello di verificare
se l’errore (difetto) presente nell’acquisizione DWI influenza la registrazione. Il
metodo per verificare ciò è quello di eliminare il rumore. Abbiamo cercato quindi
di sviluppare un modo che fosse il più possibile automatizzato non introducendo
nessuna operazione manuale per l’eliminazione dell’errore nella immagine DWI.
Ricapitolando utilizziamo una immagine senza rumore prodotta solamente utilizzando la parte sinistra della DWI. L’immagine simmetrica sfrutta appunto la
simmetria del corpo umano per ottenere cosı̀ una immagine molto simile alla
DWI. L’immagine simmetrica viene registrata sulla immagine DWI. Il risultato
viene poi registrato sulla DCE. I parametri prodotti saranno la trasformazione
da applicare alla DWI per ottenere il risultato finale. Si nota che durante la
47
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
prima registrazione “simmetrica to dwi” l’immagine di riferimento è proprio la
DWI che contiene la distorsione. Si nota anche però dai risultati che la registrazione produce una buona immagine che non contiene la distorsione DWI e che
quindi sicuramente non introduce errore nel processo di registrazione.
Otteniamo cosı̀ una nuova registrazione “dwi to dce” che contiene ancora la
distorsione ma che è stata generata in modo libero da essa. Prima di confrontare il risultato con la registrazione diretta “dwi-dce” decidiamo di eliminare
definitivamente l’errore che ora non ci serve più. Per l’eliminazione dell’errore abbiamo utilizzato una maschera generata con l’algoritmo di region growing
sull’immagine DCE (filtrata con il filtro anisotropico per un risultato più pulito). La moltiplicazione diretta (effettuata con il tool fslmaths di FSL) della
maschera ottenuta con i risultati delle registrazioni “dwi to dce” ci dà il risultato. In figura 4.6 vediamo una estensione di 4.5 con l’aggiunta di eliminazione
dell’errore con maschera.
Figura 4.6: Eliminazione dell’errore con applicazione di maschera generata utilizzando l’algoritmo di region growing sull’immagine DCE filtrata con il filtro
Anisotropico 3D
I risultati ottenuti verranno prodotti in dettaglio nella sezione Risultati ma
possiamo anticipare che la differenza fra le due registrazioni (registrazione diretta e registrazione utilizzando simmetrica) è minima, questo sta a significare che
anche se l’immagine DWI contiene una distorsione questa non produce errori
durante la fase di registrazione e quindi si ripercuote solo minimamente.
Per l’applicazione del filtro anisotropico ed in seguito per la creazione della
maschera 3D (applicata come parametro di registrazione) per i motivi descritti
precedentemente (eliminazione della distorsione) sono stati creati degli script
Matlab per automatizzare sia i processi di caricamento delle immagini che le
operazioni di filtraggio e di segmentazione (ricordiamo che tutte le immagini
sono salvate nel formato .nii che in Matlab viene descritto con una struttura con
4 campi, da cui da uno di essi si estrae la matrice che rappresenta l’immagine). Il
48
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
filtraggio anisotropico viene applicato al volume dell’acquisizione dce-21 tramite
l’algoritmo visto nel capitolo precedente.
Figura 4.7: Applicazione del filtro anisotropico sull immagine dce-21
I seed per il Region Growing sono stati calcolati inizialmente (sempre sulla
l’acquisizione dce-21) in modo manuale attraverso individuazione delle zone di
interesse. Questo procedimento però non permetteva di creare una maschera
accettabile per tutte le fette del volume in quanto alcune di esse (come naturale
nelle acquisizioni di risonanza) non contenevano più le informazioni utilizzate
per fissare i nostri seed. Considerato questo incoveniente si è proceduto in modo
complementare a quanto appena detto: è stata quindi costruita una maschera
sullo sfondo (background) per ogni fetta del volume, al fine di conservare ogni
informazione utile per ognuna di essa.
49
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.8: Maschera dell’immagine dce-21
4.1.2
Risultati
I risultati qui mostrati sono rappresentati nelle tabelle e nelle immagini seguenti divisi secondo le diverse trasformazioni utilizzate nelle registrazioni. Dopo
un’analisi soggettiva della registrazione e, attraverso uno strumento fornito dal
programma 3DSlicer, da un confronto diretto con l’immagine di riferimento 4.2
abbiamo convertito il risultato ottenuto nel formato DICOM. La riconversione a
questo formato ci permette di eseguire il programma Tester che, come descritto
nel capitolo precedente, permette un’analisi oggettiva e quantitativa in termini
di Mutual Information, Differenza di Intensità ed Edge Overlap. Per quanto
riguarda le immagini relative al seno i primi due indici sono poco significativi
perchè la registrazione è multimodale (es. DWI su DCE). Interessante invece è
il confronto fra i gradienti che quantifica l’errore della sovrapposizione (si effettuano calcoli sulla distanza dei contorni delle due immagini). Successivamente,
abbiamo estrapolato i dati ottenuti dall’esecuzione del programma e attraverso
uno script Matlab (visto sempre nel capitolo 3), normalizzati e visualizzati in
un grafico per avere un’informazione generale del lavoro effettuato.
Trasformazione Rigida
Con elastix i risultati ottenuti dal confronto della registrazione con il volume
di riferimento utilizzato (dce-21) sono riportati nella tabella seguente:
50
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Numero
1
2
3
4
5
6
7
8
9
10
11
12
Mutual Information
0.107703
0.109381
0.108304
0.108314
0.112827
0.112956
0.112985
0.113186
0.111759
0.111026
0.111065
0.111634
Differenza di intensità
0.293771
0.293170
0.293523
0.293346
0.288628
0.288534
0.288789
0.288471
0.291824
0.291507
0.292640
0.291492
Edge Overlap
0.273984
0.273946
0.273648
0.273647
0.260448
0.260400
0.260655
0.260681
0.265628
0.265875
0.266425
0.266284
Tabella 4.4:
Risultati ottenuti nella registrazione Rigida in elastix,
corrispondenti alle prime 12 prove riportate in tabella 4.1
Figura 4.9: Grafico in norma 2: parametrizzazioni rigide utilizzando elastix
nelle ascisse ed errore nelle ordinate
Confrontando i risultati riportati in figura 4.9, anche alla luce di quelli ottenuti con le altre modalità di trasformazione, possiamo vedere come questo
tipo di trasformazione produca buone registrazioni (numeri: 5,6,7,8) quando
abbiamo utilizzato un’interpolazione lineare che permette al contrario dell’interpolazione BSpline di non alterare i valori d’intensità dell’immagine ottenuta
e inoltre conservare una buona risoluzione, al contrario dell’interpolazione NearestNeighbour. Anche ad un confronto visivo è possibile notare la bontà della
registrazione ottenuta.
51
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.10: Registrazione rigida
Figura 4.11: Confronto tra registrazione rigida e volume di riferimento
Trasformazione Affine
La trasformazione affine ha fornito i risultati migliori da un punto di vista soggettivo e anche dai risultati ottenuti grazie al Tester (grafico 4.12) possiamo
notare che si ripresenta la situazione vista nel caso della trasformazione rigida.
Infatti le registrazioni migliori sono state effettuate utilizzando l’interpolazione
lineare tra i punti dell’immagine fissa e quella mobile (Numero 5,6,7,8).
52
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Numero
1
2
3
4
5
6
7
8
9
10
11
12
Mutual Information
0.112423
0.109469
0.095731
0.109826
0.112104
0.113391
0.112048
0.113216
0.110861
0.112423
0.109930
0.109905
Differenza di intensità
0.286743
0.287128
0.274874
0.287638
0.281345
0.284664
0.281941
0.283576
0.284156
0.286743
0.282332
0.286566
Edge Overlap
0.267715
0.276869
0.281985
0.276208
0.262907
0.261173
0.262403
0.261561
0.268853
0.267715
0.26924
0.267502
Tabella 4.5:
Risultati ottenuti nella registrazione Affine in elastix,
corrispondenti alle prove da 13 a 24 riportate in tabella 4.1
Figura 4.12: Grafico in norma 2: parametrizzazioni affine utilizzando elastix
nelle ascisse ed errore nelle ordinate
In particolare la registrazione numero 6 è risultata particolarmente buona
sia da un punto di vista quantitativo che qualitativo e rappresenta il punto di
partenza per un’ulteriore raffinamento del processo di registrazione. Possiamo
vedere dalle immagini seguenti come la registrazione rispecchia i contorni e
le informazioni del volume usato come riferimento nonostante la presenza del
fenomeno di ghosting che tratteremo, in termini di risultati ottenuti, in seguito.
Come vedremo analizzando le immagini della registrazione deformabile (o
53
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.13: Registrazione affine
Figura 4.14: Confronto tra registrazione affine e volume di riferimento
BSpline), le trasformazioni rigide e affine, che comportano una registrazione
lineare, sono ottimali per effettuare una registrazione multimodale. La trasformazione B-Spline rispetto ad una trasformazione lineare possiede molti più
gradi di libertà. Ricordando alcuni dei fondamenti teorici la registrazione può
essere vista come il problema di trovare una minimizzazione per una data metrica. La trasformazione BSpline permette di trovare tale minimizzazione con
54
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
un elevato grado di libertà. Questo causa alcune volte una deformazione eccessiva dell’immagine moving. Per capire questo processo è necessario considerare
lo spostamento di un punto (nell’immagine moving) in una posizione dell’immagine fixed in modo da minimizzare il valore dell’errore secondo la metrica
selezionata. Questo processo idealmente buono perde di “stabilità” se considerato globalmente. Infatti se ogni punto si muove indipendentemente dagli
altri otteniamo un risultato errato. In realtà la trasformazione BSpline è uno
strumento “potente” (che attraverso uno sviluppo più approfondito di quello
affrontato da noi) può essere settato e limitato per arginare questo effetto di
deformazione. Considerando le prove effettuate possiamo concludere che trasformazioni BSpline funzionano solo se effettuate con registrazioni monomodali.
Cioè quando le informazioni (i valori di intensità) contenute nelle immagini di
input sono molto simili e sono coerenti.
Trasformazione BSpline
Come premesso sopra, la trasformazione BSpline (figura 4.16) da un punto di
vista soggettivo, crea una registrazione errata. Ma naturalmente da un punto
di vista oggettivo e numerico, la registrazione ha avuto successo in alcuni casi
anche più delle precedenti Affine e Rigida. Il valore delle metriche si è minimizzato al massimo: per esempio se consideriamo la metrica “Mutual Information”
utilizzata nelle prove 1,5 e 9 otteniamo un evidente miglioramento nel campo
“Mutual Information” del programma Tester (grafico 4.15). Tutto questo non
ha però prodotto il risultato migliore sperato. La minimizzazione dell’errore
dettato dalla Metrica ha causato una eccessiva traslazione dei voxel coinvolti
causando una eccessiva deformazione.
Numero
1
2
3
4
5
6
7
8
9
10
11
12
Mutual Information
0.155688
0.141059
0.151939
0.141717
0.156958
0.148296
0.158923
0.150727
0.167372
0.155688
0.169669
0.16408
Differenza di intensità
0.238078
0.213782
0.206935
0.203801
0.212619
0.204777
0.203963
0.214584
0.236604
0.238078
0.233286
0.238342
Edge Overlap
0.268917
0.265236
0.269406
0.265269
0.263418
0.260844
0.263209
0.262315
0.273168
0.268917
0.272744
0.270281
Tabella 4.6: Risultati ottenuti nella trasformazione BSpline in elastix
55
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.15: Grafico in norma 2: tutte le parametrizzazioni con elastix nelle
ascisse ed errore nelle ordinate
Figura 4.16: Registrazione BSpline
56
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Raffinamento
Come abbiamo visto, la registrazione migliore è stata ottenuta con la trasformazione affine e nella nostra ricerca il risultato migliore si è avuto nella registrazione
affine numero 6. A partire da questa configurazione abbiamo iterato un processo di raffinamento attraverso la modifica dei parametri meno importanti ma
comunque non trascurabili per la registrazione. In parametri settati in questo
nuovo set di registrazioni, riguardano, come abbiamo visto, il tipo di metrica e
il numero di iterazioni dell’ottimizzazione. Anche in questo caso abbiamo utilizzato sia l’approccio soggettivo che quello oggettivo analizzando e raggruppando
i dati ricavati dalla funzione Tester come segue:
Numero
0
1
2
3
4
5
6
7
8
Mutual Information
0.0743407
0.0921767
0.0708796
0.0982557
0.11347
0.108603
0.0918939
0.0999608
0.0923162
Differenza di intensità
0.263661
0.25877
0.261671
0.28485
0.284559
0.292234
0.276577
0.262783
0.268545
Edge Overlap
0.285718
0.278022
0.285021
0.264495
0.261021
0.26032
0.273187
0.273123
0.276109
Tabella 4.7: Risultati ottenuti nella fase di raffinamento dei risultati in elastix
Dai dati ottenuti possiamo affermare che la migliore registrazione (con trasformazione affine) ottenuta avviene utilizzando la metrica della Mutua Informazione Normalizzata con 3 iterazioni di ottimizzazione (quindi quella di partenza).
In figura 4.13 abbiamo l’immagine finale ottenuta e in 4.14 un confronto con la
dce-21.
Fsl
I risultati ottenuti con FSL sono molto simili ai risultati ottenuti con elastix.
Tutte le considerazioni fatte nel caso elastix si rispecchiano su FSL. Anche
qui la registrazione avviene in modo corretto solo utilizzando una registrazione
lineare. Di seguito riportiamo una tabella contenente i risultati ottenuti con il
programma Tester utilizzando i parametri precedentemente introdotti:
57
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Numero
1
2
3
4
5
6
7
8
9
10
11
12
Mutual Information
0.103874
0.092814
0.0944517
0.0964465
0.103093
0.0822152
0.0916931
0.0951526
0.0967916
0.0919036
0.094093
0.0966201
Differenza di intensità
0.295117
0.299786
0.295776
0.295663
0.293827
0.298358
0.294432
0.294354
0.290351
0.294853
0.291023
0.290855
Edge Overlap
0.274743
0.273641
0.274376
0.274356
0.266619
0.266058
0.266805
0.26648
0.261482
0.261059
0.261673
0.261463
Tabella 4.8: Risultati ottenuti in FSL
Considerando questi valori e rapportandoli ai valori ottenuti con elastix a
parità di parametri notiamo una elevata similarità. Anche ad una visione soggettiva dei risultati rapportati a quelli di elastix otteniamo la conferma che i
due tool funzionano utilizzando algoritmi simili. Nel grafico 4.17 possiamo osservare come l’errore tende a migliorare verso le ultime configurazioni. Tuttavia
l’errore visualizzato è in norma 2 questo evidenzia notevolmente una differenza
di valutazione fra le varie parametrizzazioni: in realtà non c’è molta differenza
nei risultati ottenuti. I tempi di computazione sono accettabili se confrontati con i tempi ottenuti in elastix. Questo ci porta a concludere che a parità
di complessità computazione i due tool si comportano in modo analogo sia per
quanto riguarda il tempo di esecuzione sia da un punto di vista qualitativo. Non
è possibile invece stilare una descrizione dei risultati della trasformazione non
lineare “fnirt” in quanto le conclusioni ottenute non mutano dalle conclusioni
di elastix.
58
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.17: Grafico degli errori in norma 2 dei test con “flirt”: numero delle
parametrizzazioni nelle ascisse ed errore nelle ordinate
Operazioni di eliminazione del rumore
Nella tabella seguente sono riportati alcuni valori degli errori ottenuti con diverse configurazioni già citate in precedenza:
Descrizione
1◦ prova
2◦ prova
3◦ prova
4◦ prova
Mutual Information
0.115566
0.106294
0.120321
0.107307
Differenza di intensità
0.0.278307
0.278288
0.27905
0.279403
Edge Overlap
0.23167
0.235993
0.238735
0.235962
Tabella 4.9: Tabella risultati dopo eliminazione del rumore
Dove:
• 1◦ prova - Applicazione della maschera ottenuta dalla dce senza filtraggio
anisotropico direttamente al risultato migliore
• 2◦ prova - Applicazione della maschera ottenuta dalla dce con filtraggio
anisotropico direttamente al risultato migliore
• 3◦ prova - Applicazione della maschera ottenuta dalla dce senza filtraggio
anisotropico al risultato ottenuto con atlante
59
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
• 4◦ prova - Applicazione della maschera ottenuta dalla dce con filtraggio
anisotropico al risultato ottenuto con atlante
I risultati confermano chiaramente il fatto che il rumore presente nella parte
destra dell’immagine dwi non influenza direttamente il risultato della registrazione. Questo è dovuto al fatto che l’applicazione della maschera non migliora i
risultati in modo significativo. La trasformazione che ci ha permesso di ottenere
il risultato migliore è stata la BSpline (numero 35 della tabella 4.1).
Figura 4.18: Registrazione affine con applicazione della maschera per
l’eliminazione del rumore
60
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.19: Confronto tra registrazione con eliminazione del rumore e volume
di riferimento
Applicazione della tecnica ad altri volumi
Per verificare l’effettiva correttezza dei dati ottenuti abbiamo ottenuto i parametri che hanno prodotto i risultati migliori su tutte le dce (possediamo 8 immagini
dce rappresentanti le varie fasi di acquisizione dopo l’immissione del liquido di
contrasto). In questo modo possiamo verificare che le conclusioni redatte nella
dce-21 siano veritiere anche nelle altre immagini. Come mostra il grafico 4.20 (la
valutazione dell’errore è positiva se il valore tende a 0) i risultati ottenuti sono
completamente uniformi su tutte le immagini dce. Questo significa che la parametrizzazione selezionata con la dce-21 funziona anche con le altre 7 immagini
dce. Dopo aver fatto questo primo passo di riutilizzo della parametrizzazione
ottenuta abbiamo provato ad espandere il campo di applicazione anche ad immagini dello stesso tipo di altri pazienti. Il risultato è stato molto positivo: nelle
figure 4.21 e 4.22 (rispettivamente il risultato migliore del 2◦ seno utilizzato e
un confronto con la rispettiva dce) è possibile prendere atto di ciò. La configurazione migliore ottenuta con questi nuovi pazienti è risultata essere molto
simile alla configurazione migliore del primo seno esaminato. Infatti otteniamo
una buona registrazione solo con una trasformazione di tipo lineare. In figura
4.21 è stata utilizzata la configurazione Rigida6. Anche con questo tipo di volumi la configurazione BSpline non dà buoni risultati deformando eccessivamente
l’immagine dwi.
61
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.20: Grafico degli errori: 8 immagini dce dello stesso paziente nell’asse
delle ascisse e valore fra 0 e 1 di errore nell’asse delle ordinate.
Figura 4.21: Risultato utilizzando la parametrizzazione Rigida6 del seno della
2◦ paziente.
4.2
Immagini di testa-collo
Le immagini di testa-collo, forniteci dalla Fondazione IRCCS [14], sono organizzate in un dataset molto ampio che contiene le varie acquisizione organizzate in
cartelle ordinate cronologicamente (giorno dell’acquisizione). Le immagini contengono lo sviluppo temporale della situazione clinica di un paziente che delinea
la presenza di elementi tumorali nei linfonodi.
Solo alcune cartelle contengono le informazioni necessarie alla nostra ricerca
e quindi dopo una comunque non rapida analisi dell’immagini abbiamo ridotto il
nostro dataset a quattro cartelle, le uniche che contenevano tutte le acquisizioni
62
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.22: Confronto fra dce e risultato utilizzando la parametrizzazione
Rigida6 del seno della 2◦ paziente.
utili (DCE,DWI,T1 e T2) :
• 08apr08 (cartella di riferimento)
• 26aug08
• 12may09
• 19aug09
Ognuna di queste cartelle si suddivide a sua volta in cartelle contenenti, come
sopra accennato, le diverse modalità di acquisizione. Le immagini DWI come
visto anche nell’altro scenario risultano a bassa risoluzione, le immagini DCE
hanno una risoluzione migliore ma risultano comunque inferiori da un punto di
vista di informazioni che si possono ricavare, alle immagini morfologiche date
dall’aquisizione T1 e quella T2. Le immagini DCE sono circa 2100 per ogni cartella principale, e rappresentano 70 volumi di acquisizione diversi (ogni volume
è costituito da 30 immagini ciascuno). Ogni volume rappresenta una sequenza
di immagini che descrivono la regione anatomica che parte dalle spalle e risale
fino alla base del cervello. Nel nostro nuovo dataset descriviamo ogni volume
con un nome identificativo dce+10.nii per esempio. La cartella di riferimento,
oovero l’aquisizione alla diagnosi della malattia, contiene le informazioni più
utili da un punto di vista medico (tumore ai linfonodi e quindi relativo ingrandimento) e consentono anche a noi di fissare delle caratteristiche e dei parametri
importanti per valutare la buona riuscita della registrazione. Le altre cartelle
contengono un miglioramento graduale delle condizioni fisiche del paziente che
ne dimostrano la guarigione.
4.2.1
Traccia di Sviluppo
Per le immagini di testa-collo abbiamo seguito una traccia di sviluppo simile a
quella utilizzata per il seno. Per questo dominio applicativo ci siamo concentrati
come per il seno su un insieme di parametri. Le registrazioni sono state effettuate
utilizzando tutte le combinazioni di questo insieme di parametri. Tale insieme
è lo stesso considerato per la registrazione del seno (tabella 4.1). In questo
caso i primi risultati hanno evidenziato che nessuna registrazione diretta fra
l’immagine DWI e DCE dava risultati qualitativamente accettabili. Un ipotesi
possibile per questo insuccesso è la grossa differenza di intensità e risoluzione
fra i due volumi usati.
63
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.23: Acquizione testa e collo DWI
Figura 4.24: Acquizione testa e collo T2
64
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.25: Acquizione testa e collo DCE
Abbiamo tentato di aggirare l’ostacolo sfruttando un supporto intermedio
con una qualità informativa maggiore. Il supporto ci è stato dato dalle immagini
di acquisizione T1 e/o T2 le cui caratteristiche fisiche sono elecante nel relativo
paragrafo. Queste acquisizioni danno un importante accrescimento dal punto di
vista della risoluzione.
Abbiamo sviluppato principalmente due linee di sviluppo:
• La prima (figura 4.26) consiste in:
1. Registrazione dell’immagine DWI sull’immagine di riferimento T2
(T1)
2. Registrazione dell’immagine T2 sull’immagine di riferimento DCE.
3. Trasformazione attraverso il tool trasformix utilizzando la (1) con i
parametri della (2).
Il confronto in questo caso verrà effettuato tra il risultato ottenuto nel
punto (3) e l’immagine di riferimento DCE.
• La seconda consiste in:
1. Registrazione dell’immagine DWI sull’immagine di riferimento T2
(T1)
2. Registrazione dell’immagine DCE sull’immagine di riferimento T2.
65
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.26: Primo approccio di registrazione utilizzando volume intermedio T1
(T2)
In questo caso invece il confronto avverrà direttamente fra i due risultati
ottenuti (figura 4.27)
Figura 4.27: Secondo approccio di registrazione utilizzando volume intermedio
T1 (T2)
Il fattore comune fra i due processi è quello di applicare delle registrazioni
intermedie fra le diverse modalità di acquisizione come si evince dalle figure 4.26
e 4.27. Le registrazioni intermedie sono state effettuate con la stessa tecnica
utilizzata per le registrazioni delle immagini del seno. In particolare abbiamo
applicato 36 registrazioni della DWI sulla T2, 36 registrazioni della T2 sulla
DCE e 36 registrazioni della DCE sulla T2.
Nel caso del primo approccio per il processo di trasformazione (trasformix)
abbiamo eseguito una combinazione dei risultati ottenuti nelle registrazioni intermedie. Questa combinazione sarà descritta in dettaglio nella relativa sezione
“Risultati”. Nel secondo approccio invece abbiamo eseguito semplicemente un
confronto fra i due insiemi di risultati intermedi. In tutti e due gli approcci
abbiamo eseguito poi un’analisi oggettiva attraverso il programma “Tester”.
Come detto all’inizio di questa sezione il volume usato per le elaborazioni
mostra la patologia al massimo della sua espansione. I vari tipi di acquisizione riescono ad enfatizzare in maniera diversa queste informazioni in base alla
66
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
loro natura e quindi abbiamo riscontrato una difficoltà maggiore nella fase di
registrazione rispetto a volumi con presenza minore del fattore patologico. Per
affrontare queste difficoltà abbiamo tentato di utilizzare diverse tecniche di filtraggio e di mascheramento che purtroppo non hanno dato gli effetti sperati.
In particolare abbiamo provato l’utilizzo di particolari maschera create con la
tecnica del region growing sui vasi: l’idea è stata quella di trovare nel modo più
automatico possibile punti in comune fra le due immagini in input nel processo
di registrazione. Inoltre, abbiamo tentato di creare una maschera analogamente
a quanto fatto con il seno (utilizzando quindi la stessa tecnica di region growing applicata al background) in modo tale da contenere l’allargamento laterale
dell’immagine T2.
4.2.2
Risultati
I risultati qui mostrati sono rappresentati nelle tabelle e nelle immagini seguenti
divisi secondo le diverse trasformazioni utilizzate nelle registrazioni. Abbiamo
anche in questo caso estrapolato i dati ottenuti dall’esecuzione del programma
Tester e attraverso uno script Matlab li abbiamo normalizzati e visualizzati in
un grafico per avere una descrizione generale del lavoro effettuato. A differenza
delle immagini del seno qui abbiamo una diversa modalità di descrizione dei
risultati all’interno del processo di registrazione totale. In particolare descriveremo uno alla volta i passi intermedi effettuati e poi descriveremo i tentativi di
miglioramento che hanno avuto successo o meno.
Registrazione della DWI sulla T2
Come si può notare dalla figura 4.23 e dalla figura 4.25 l’acquisizione DWI
presenta una risoluzione minore dal punto di vista informativo della acquisizione DCE. Come già detto l’acquisizione T2 è densa di informazioni ed è ad
alta risoluzione quindi questa peculiarità ci consente di utilizzarla come supporto nel processo di registrazione, come passo intermedio nella registrazione
totale DWI to DCE. Come ben descritto nella traccia di sviluppo costruiamo
una registrazione usando come immagini fissa la T2 e come immagine mobile
la DWI, questa prima registrazione viene effettuata in entrambi i processi di
realizzazione (refconf1 e 4.27). Per fare ciò, abbiamo utilizzato la stessa tecnica
utilizzata nell’altro esperimento cioè abbiamo effettuato 36 registrazioni frutto
della combinazione lineare dell’insieme dei parametri descritti in tabella 4.1.
La miglior registrazione in questo primo passo è stata ricavata con la trasformazione Affine (parametrizzazione numero 18 della tabella 4.1) ed è visibile
in figura 4.28. Considerata la notevole differenza di risoluzione e in generale
di contenuto informativo possiamo dire che il primo passo intermedio ha dato
buoni risultati e soprattutto ci fornisce una prima trasformazione che ci avvicina all’immagine DCE e ci da quindi una ottima base per effettuare la seconda
registrazione.
67
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.28: Registrazione Affine DWI to T2
(a) Confronto tratto cortico spinale
(b) Confronto visione assiale
Figura 4.29: Confronto fra registrazione DWI-T2 e T2
Registrazione della T2 sulla DCE
Nel caso del primo approccio (figura 4.26) costruiamo poi una seconda registrazione usando come immagini fissa la DCE e come immagine mobile la T2.
Per fare ciò, abbiamo utilizzato la stessa tecnica utilizzata nell’altro esperimento cioè abbiamo effettuato 36 registrazioni frutto della combinazione lineare
dell’insieme dei parametri descritti in tabella 4.1.
Anche in questo caso la miglior registrazione in questo primo passo è stata
ricavata con la trasformazione Affine (parametrizzazione numero 18 della tabella 4.1) ed è visibile in figura 4.30. Questa seconda trasformazione ci ha fornito
i parametri per il passaggio finale da applicare alla registrazione eseguita precedentemente (DWI to T2). In figura 4.31 si nota che la registrazione non è
68
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
ottima in quanto i bordi dei linfonodi (soprattutto nella parte sinistra dell’immagine) non sono perfettamente “sovrapponibili”. Questo errore si ripercuoterà
nel passaggio finale come vedremo nei prossimi paragrafi.
Figura 4.30: Registrazione Affine T2 to DCE
(a) Confronto tratto cortico spinale
(b) Confronto visione assiale
Figura 4.31: Confronto fra registrazione T2 e DCE
69
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Registrazione della DCE sulla T2
Nel caso del secondo approccio (figura 4.27) costruiamo invece una seconda registrazione usando come immagini fissa la T2 e come immagine mobile la DCE.
Per fare ciò, abbiamo utilizzato la stessa tecnica utilizzata nell’altro esperimento cioè abbiamo effettuato 36 registrazioni frutto della combinazione lineare
dell’insieme dei parametri descritti in tabella 4.1.
La registrazione ottenuta non ha soddisfatto le nostre aspettative in quanto
da un punto di vista soggettivo le immagine risultante non dava delle informazioni essenziali per un possibile confronto con la prima registrazione (intensità
troppo forte che nascondeva i casi patologici).
Figura 4.32: Registrazione Affine DCE to T2
L’insuccesso di questa possibile via di sviluppo, cioè quella basata sul confronto diretto fra le registrazioni DWI to T2 e DCE to T2 ci ha indirizzato a
concentrarci sul primo approccio descritto (figura 4.26). Nei prossimi paragrafi saranno descritti i risultati ottenuti dall’ultimo passo di questo processo: la
trasformazione con il tool trasformix. Verranno quindi elencate le combinazioni
effettuate e verranno esposti i motivi delle scelte adottate. Per limitare i tempi
di computazione e di valutazione dei risultati (una registrazione con relativa
fase di testing può durare anche alcune ore) abbiamo cercato di generare delle trasformazioni in qualche modo collegate fra loro. Dopo una attenta analisi
abbiamo quindi scelto di applicare trasformazioni dello stesso tipo delle registrazioni del primo passo in modo da non stravolgere lo sviluppo della registrazione.
Di seguito elenchiamo le trasformazioni effettuate:
• Trasformazione Rigida-Rigida: Trasformazione di tipo Rigida (con
i parametri generati dalla registrazione con trasformazione Rigida T2
to DCE) sul risultato della registrazione con trasformazione Rigida nel
passaggio DWI to T2.
70
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
• Trasformazione Affine-Affine: Trasformazione di tipo Affine (con i parametri generati dalla registrazione con trasformazione Affine T2 to DCE)
sul risultato della registrazione con trasformazione Affine nel passaggio
DWI to T2.
• Trasformazione BSpline-BSpline: Trasformazione di tipo BSpline (con
i parametri generati dalla registrazione con trasformazione BSpline T2
to DCE) sul risultato della registrazione con trasformazione BSpline nel
passaggio DWI to T2.
• Trasformazione Affine-BSpline: Trasformazione di tipo BSpline (con
i parametri generati dalla registrazione con trasformazione Rigida T2
to DCE) sul risultato della registrazione con trasformazione Affine nel
passaggio DWI to T2.
71
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Trasformazione Rigida-Rigida
Con elastix i risultati prodotti da questa trasformazione calcolati dal programma Tester sono riportati nella tabella seguente, dove abbiamo rappresentato
solo i dati relativi alle trasformazioni che hanno generato da un punto di vista soggettivo un buon risultato (l’esperienza maturata con il lavoro effettuato
sulle immagini del seno ci ha insegnato che effettuare dei confronti oggettivi
su immagini dal contenuto informativo diverso è controproducente sia a livello
computazionale che a livello di espressività dei dati ottenuti).
Come si evince dalla tabella, la trasformazione che ci ha fornito i risultati
Numero
1
2
3
4
5
6
7
8
9
10
11
12
Mutual Information
—–
0.148491
—–
—–
—–
0.164397
—–
0.164039
0.149607
—–
—–
—–
Differenza di intensità
—–
0.139755
—–
—–
—–
0.137748
—–
0.137423
0.13887
—–
—–
—–
Edge Overlap
—–
0.489498
—–
—–
—–
0.478968
—–
0.481055
0.485994
—–
—–
—–
Tabella 4.10: Risultati ottenuti in elastix
Figura 4.33: Grafico degli errori in norma 2 dei test con elastix : numero delle
parametrizzazioni nelle ascisse ed errore nelle ordinate
72
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
migliori è stata la numero 6 (Fig. 4.34). Un primo approccio per effettuare
il confronto, per ottenere una valutazione soggettiva del risultato, è quello di
confrontare i tratti cortico spinali delle immagini coinvolte. In particolare i due
tratti cortico spinali devono essere coincidenti confrontandoli sul piano sagittale. Quando il risultato è “accettabile” si passa ad una seconda valutazione
più approfondita guardando le due visioni assiali. Nel nostro caso una buona
sovrapposizione dei linfonodi rappresenta un risultato adeguato alle aspettative,
in quanto strutture oggetto di indagine per queste acquisizioni. I risultati riportati mostrano come preannunciato che l’errore introdotto nella registrazione
T2 to DCE si ripercuote sul risultato finale peggiorandolo.
Figura 4.34: Trasformazione Rigida
(a) Confronto tratto cortico spinale
(b) Confronto visione assiale
Figura 4.35: Confronto fra trasformazione rigida-rigida DWI-DCE e DCE
73
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Trasformazione Affine-Affine
Con elastix i risultati prodotti da questa trasformazione calcolati dal programma Tester sono riportati nella tabella seguente:
Numero
1
2
3
4
5
6
7
8
9
10
11
12
Mutual Information
—–
—–
—–
—–
—–
0.183706
—–
0.178569
0.163765
0.15804
0.163352
0.161221
Differenza di intensità
—–
—–
—–
—–
—–
0.126964
—–
0.12808
0.13178
0.134077
0.130708
0.13207
Edge Overlap
—–
—–
—–
—–
—–
0.479107
—–
0.477231
0.484573
0.485488
0.4844887
0.485207
Tabella 4.11: Risultati ottenuti in elastix
Come si evince dalla tabella, la trasformazione che ci ha fornito i risultati mi-
Figura 4.36: Grafico degli errori in norma 2 dei test con elastix : numero delle
parametrizzazioni nelle ascisse ed errore nelle ordinate
gliori è stata la numero 6 (Fig. 4.37). Per il confronto diretto con gli strumenti
messi a disposizione dal programma 3DSlicer utilizziamo lo stesso approccio
descritto nella sezione della trasformazione Rigida-Rigida
74
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.37: Trasformazione Affine-Affine
(a) Confronto tratto cortico spinale
(b) Confronto visione assiale
Figura 4.38: Confronto fra trasformazione affine-affine DWI-DCE e DCE
Trasformazione BSpline-BSpline
Abbiamo già visto nella sezione dedicata al seno come la trasfromazione non
lineare per immagini di diversa acquisizione, non fornisce risultati apprezzabili.
Sia con FSL che con elastix le immagini prodotte non soddisfano gli obiettivi
da noi cercati. Questa tecnica che ricordiamo viene riproposta due volte nell’arco del processo da noi seguito amplifica l’errore “esponenzialmente”. Non
abbiamo quindi effettuato verifiche con il programma Tester per non eseguire
una computazione lenta e inutile.
Trasformazione Affine-BSpline
Con elastix i risultati prodotti da questa trasformazione calcolati dal programma Tester sono riportati in tabella 4.12.
75
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Numero
1
2
3
4
5
6
7
8
9
10
11
12
Mutual Information
—–
—–
—–
—–
—–
0.72514
—–
0.159616
0.16349
—–
—–
—–
Differenza di intensità
—–
—–
—–
—–
—–
0.130623
—–
0.13255
0.131528
—–
—–
—–
Edge Overlap
—–
—–
—–
—–
—–
0.480812
—–
0.488556
0.489551
—–
—–
—–
Tabella 4.12: Risultati ottenuti in elastix
Figura 4.39: Grafico degli errori in norma 2 dei test con elastix : numero delle
parametrizzazioni nelle ascisse ed errore nelle ordinate
Come si evince dalla tabella, la trasformazione che ci ha fornito i risultati migliori è stata la numero 6 (Fig. 4.40). I risultati riportati mostrano come anche
in questo caso l’errore introdotto nella registrazione T2 to DCE si ripercuote sul
risultato finale peggiorandolo. Abbiamo già discusso di come la trasformazione
BSpline non produce un risultato accettabile, anche qui, pur usando questo tipo di trasformazione solo per una parte del processo complessivo, otteniamo i
medesimi risultati.
76
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
Figura 4.40: Trasformazione Affine-BSpline
(a) Confronto tratto cortico spinale
(b) Confronto visione assiale
Figura 4.41: Confronto fra trasformazione Affine-BSpline DWI-DCE e DCE
77
Conclusioni
Gli obiettivi principali perseguiti durante questa tesi sono stati i seguenti:
1. Analisi e comprensione dello stato dell’arte nell’analisi delle immagini biomedicali ed, in particolare, del ruolo svolto dalla registrazione all’interno
di essa.
2. Studio di due tool parametrici “Elastix” ed “Fsl” e loro confronto.
3. Sviluppo dell’approccio più adatto per la registrazione multimodale di
immagini del seno e del distretto (testa e collo) di alcuni casi patologici, per
ottenere nei limiti del possibile una metodologia “standard” riapplicabile
a situazioni cliniche dello stesso tipo.
Lo studio dello stato dell’arte relativo alla registrazione ci ha dato le conoscenze adeguate, da un lato, per comprendere come utilizzare i tool parametrici presentati. Dall’altro lato, invece, ci ha dato la possibilità di individuare
problematiche ancora aperte su tale tecnica.
In primo luogo ci siamo concentrati sulle acquisizioni DWI e DCE del seno di un particolare paziente. Il nostro obiettivo è stato quello di effettuare
la registrazione dell’immagine DWI sull’immagine DCE. Questo ci ha portato
allo sviluppo di diverse vie risolutive che in molti casi non hanno portato al
raggiungimento dei risultati attesi: in particolar modo la nostra conclusione è
che in un approccio di questo genere di tipo intramodale la soluzione migliore
è quella della trasformazione di tipo lineare. Per trasformazione lineare intendiamo naturalmente tutti gli approcci rigidi e affini presentati. Abbiamo visto
come pur parlando di parti anatomiche identiche la registrazione multimodale
non può mai essere effettuata nello stesso modo. Possiamo quindi tracciare un
insieme ristretto di parametri da utilizzare ma non possiamo in alcun modo fissarli, questo perchè la diversità umana influenza notevolmente l’intero processo
(soprattutto in soggetti malati, dove la fenomenologia può essere molto diversa). Per esempio possiamo notare quanto detto nelle figure 4.2 e 4.21. In nostro
lavoro di valutazione dei risultati ottenuti è avvenuto complementariamente in
due modi: una prima valutazione di tipo soggettivo con il confrotto diretto del
risultato sull’immagine di riferimento utilizzando particolari strumenti offerti
dal programma 3DSlicer e una seconda valutazione più oggettiva attraverso la
valutazione di alcune metriche più o meno significative.
La conclusione più importante e considerevole a cui siamo arrivati è come la
trasformazione non lineare, pur essendo “sulla carta” la miglior soluzione, non
si sia rilevata utile ai nostri scopi (figura 4.16). Il motivo, per niente intuitivo, è
da ricercare nel fatto che la trasformazione BSpline utilizzata è talmente efficace che modifica l’immagine moving in maniera “drastica”, deformandola. Tale
79
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
trasformazione localmente ha la capacità di trovare la corrispondenza dei voxel
fra le due immagini considerate in modo da minimizzare al massimo la metrica
utilizzata. Questo processo, idealmente ottimale, se esteso sul risultato globale
fà comprendere come il risultato possa diventare errato. In una registrazione
che coinvolge immagini con differenze molto evidenti di intensità (l’esempio più
evidente di questo è quando si pensa al tratto cortico spinale visto nella sezione
di test-collo: assume valore di intensità opposti nelle acquisizioni DWI e DCE)
questo effetto è devastante ecco perchè nelle nostre prove di tipo multimodale la
registrazione non lineare non ha dato i risultati sperati. Questa conclusione ha
un’altro effetto essenziale sulla valutazione dei risultati: i risultati migliori sono
stati ritenuti tali perchè non solo avevano il più basso errore di Edge Overlap
ma anche perchè soddisfacevano una valutazione soggettiva diretta. Nelle trasformazioni BSpline l’errore si abbassava notevolmente rispetto a quella lineare
(appunto perchè effettivamente la metrica era stata minimizzata al massimo) ma
le immagini risultanti non permettevono un riscontro adeguato con le immagini
di rifermento. La tecnica BSpline si è rilevata utile invece quando abbiamo affrontato il problema della registrazione monomodale fra l’immagine simmetrica
e la DWI. In questo caso l’informazione presente nelle due immagini era molto
simile e un approccio non lineare ci ha permesso di ottenere un risultato ottimo.
Questa registrazione monomodale si è rilevata utile anche per evidenziare come
non solo il tipo di trasformazione ma anche altri parametri come l’interpolazione
e la metrica cambiano drasticamente a seconda del tipo di immagini coinvolte.
Nella seconda parte della tesi è stato trattato il problema della registrazione
della stessa modalità di acquisizione ma di un’altra parte anatomica: la regione
tasta-collo. Questa parte di sviluppo è stata molto più difficoltosa della precedente in quanto le immagini in nostro possesso non erano di elevatà qualità e
risoluzione. Per generare una registrazione corretta siamo dovuti ricorrere ad
un passaggio intermedio utilizzando una immagine ad alta risoluzione (acquisizione T1 e/o T2) per migliorare la qualità e per aumentare le informazioni
necessarie nel processo. In ogni modo abbiamo potuto consolidare le conclusioni già ottenute nel seno, per cui la registrazione che anche in questo caso è
stata di tipo multimodale è avvenuta utilizzando lo stesso insieme di parametri.
Per questo tipo di registrazione ci siamo concentrati solo su un paziente questo
perchè avevamo la possibilità di utilizzare volumi acquisiti in diversi periodi.
Uno dei punti deboli dello stato dell’arte è sicuramente il fatto che non
tutte le parti anatomiche del corpo umano sono trattate in modo completo ed
esaustivo in egual misura, nel nostro caso possiamo dire che il problema della
registrazione del cervello è molto più discusso rispetto a quello del seno.
Altro punto fondamentale della nostra tesi è l’utilizzo simultaneo dei due tool
FSL ed Elastix. A parità di parametri utilizzati (quando naturalmente questo è
stato possibile) abbiamo ottenuto sostanzialmente i medesimi risultati a meno
di qualche piccolo caso particolare. Questo ci indica come i due tool si basano su
algoritmi per l’effettuazione della trasformazione molto simili. Tuttavia mentre
l’utilizzo di Elastix porta quasi sempre a risultati concreti (talvolta anche errati),
l’utilizzo di parametri non corretti in FSL non dà sempre un output visibile sul
quale dedurre possibili soluzioni. Notoriamente FSL ha maturato nella sua
creazione una particolare tendenza per le immagini del cervello, questo però
non ha migliorato le nostre registrazioni in quanto ricordiamo che le immagini
in nostro possesso erano quelle di testa e collo (all’interno di FSL troviamo
particolari tool che permettono l’estrazione diretta del cervello dai volumi in
80
CAPITOLO 4. REALIZZAZIONE E SVILUPPO
ingresso per effettuare eventuali maschere, tale tecnica non ha funzionato nelle
nostre prove). Anche da un punto di vista dell’utilizzo, che ricordiamo è per
tutti e due a linea di comando, risulta essere nettamente migliore Elastix che
permette il richiamo di file “.txt” contenenti le parametrizzazioni volute. Inoltre
con Elastix è possibile definire una quantità maggiore di parametri rispetto a
FSL, che vanno dalla configurazione del tipo di campionamento alla modalità
di ottimizzazione. In questo modo Elastix acquisisce una certa flessibilità non
garantita da FSL che limita (nella versione con l’interfaccia grafica in modo
notevole) la possibilità di “movimento” all’interno delle configurazioni possibili.
Concludendo possiamo quindi considerare Elastix migliore rispetto a FSLma
questo giudizio è derivato dalla nostra esperienza e quindi di natura personale
ma comunque supportato dai risultati ottenuti in questa ricerca.
Ovviamente tutte le soluzioni proposte sono probabilmente migliorabili a
condizione di aumentare il costo computazionale dell’intero processo: si pensi ad
esempio ad un campionamento “Full” al posto di quello “Random” utilizzato. Il
nostro studio si è soffermato solo sull’utilizzo di due tool ma esiste un panorama
molto ampio di programmi che risolvono questo tipo di problematiche.
Le tecniche di registrazione, di validazione e di analisi dei risultati possono
essere il punto di partenza per un ulteriore processo di raffinamento e quindi
per un possibile sviluppo futuro. In particolare il processo sviluppato per la
registrazione delle immagini del seno, l’utilizzo dell’atlante come possibile strumento di miglioramento per alcuni problemi delle immagini è già ampiamente
adottato nelle attuali ricerche avanzate e ben descritto in letteratura [16] [17]
[18]. Mentre nella nostra ricerca abbiamo cercato di sfruttare la simmetria del
corpo come strumento di ottimizzazione per la definizione dei contorni e per
l’eliminazione del rumore, l’idea di utilizzare questa tecnica è utile per delineare
caratteristiche e proprietà fisiologiche che possono essere tralasciate usando un
processo “standard”. Le varie tecniche di registrazione sono in continuo sviluppo e stanno diventando una componente costante nell’imaging medico. La
registrazione dell’immagini biomediche costituisce quindi un’importante strumento diagnostico è il suo sviluppo e seguito parallelamente dall’analisi medica
dei risultati.
Da un punto di vista prettamente informatico questa esperienza che ha coinciso anche con il periodo tirocinio svolto all’interno del laboratorio VIPS della
facoltà di Scienze MM.FF.NN dell’Università di Verona, ci ha permesso di sviluppare algoritmi e tecniche avanzate di elaborazione digitale dell’immagini. In
particolare abbiamo affinato strumenti e conoscenze sulla manipolazione delle
immagini biomedicali e sull’utilizzo dei tool e programmi usati attualmente anche in ambiente di ricerca e industriale. Per un buon conseguimento dei risultati
abbiamo considerato varie aspetti dei strumenti utilizzati, come la complessità
degli algoritmi utilizzati, il modello di validazione costruito e il procedimento
di pianificazione del progetto. A conclusione di questo elaborato possiamo affermare di ritenerci soddisfatti per le conoscenze acquisite durante il periodo
di lavoro, dove abbiamo anche sperimentato concretamente il lavoro in team,
inoltre abbiamo avuto la possibiltà di utilizzare software e concetti teorici visti
durante il ciclo di studi in maniera approfondita su un campo di ricerca in pieno
sviluppo.
81
Appendice A
Codice del filtro
anisotropico 3D
5
10
15
20
25
30
f u n c t i o n d i f f v o l = a n i s o d i f f 3 D ( v o l , n u m i t e r , d e l t a t , kappa , o p t i o n ,
voxel spacing )
%ANISODIFF3D C o n v e n t i o n a l a n i s o t r o p i c d i f f u s i o n
%
DIFF VOL = ANISODIFF3D(VOL, NUM ITER, DELTA T, KAPPA, OPTION,
VOXEL SPACING) p e r f o m s
%
c o n v e n t i o n a l a n i s o t r o p i c d i f f u s i o n ( Perona & Malik ) upon a s t a c k o f
gray s c a l e images .
%
A 3D network s t r u c t u r e o f 26 n e i g h b o r i n g n o d e s i s c o n s i d e r e d f o r
d i f f u s i o n conduction .
%
%
ARGUMENT DESCRIPTION :
%
VOL
− g r a y s c a l e volume d a t a (MxNxP) .
%
NUM ITER − number o f i t e r a t i o n s .
%
DELTA T − i n t e g r a t i o n c o n s t a n t ( 0 <= d e l t a t <= 3 / 4 4 ) .
%
U s u a l l y , due t o n u m e r i c a l s t a b i l i t y t h i s
%
p a r a m e t e r i s s e t t o i t s maximum v a l u e .
%
KAPPA
− g r a d i e n t modulus t h r e s h o l d t h a t c o n t r o l s t h e
conduction .
%
OPTION
− c o n d u c t i o n c o e f f i c i e n t f u n c t i o n s p r o p o s e d by
Perona & Malik :
%
1 − c ( x , y , z , t ) = exp ( −( n a b l a I / kappa ) . ˆ 2 ) ,
%
p r i v i l e g e s high−c o n t r a s t e d g e s o v e r low−
c o n t r a s t ones .
%
2 − c ( x , y , z , t ) = 1 . / ( 1 + ( n a b l a I / kappa ) . ˆ 2 ) ,
%
p r i v i l e g e s wide r e g i o n s o v e r s m a l l e r o n e s
.
%
VOXEL SPACING − 3 x1 v e c t o r column w i t h t h e x , y and z
dimensions of
%
the voxel ( m i l i m e t e r s ) . In p a r t i c u l a r , only
c u b i c and
%
a n i s o t r o p i c v o x e l s i n t h e z−d i r e c t i o n a r e
considered .
%
When d e a l i n g w i t h DICOM images , t h e v o x e l
spacing
%
d i m e n s i o n s can be e x t r a c t e d u s i n g MATLAB’ s
dicominfo ( . ) .
%
%
OUTPUT DESCRIPTION :
%
DIFF VOL − ( d i f f u s e d ) volume w i t h t h e l a r g e s t s c a l e −
space parameter .
% Convert i n p u t volume t o d o u b l e .
vol = double ( vol ) ;
% Useful variables .
[ rows c o l s p a g s ] = s i z e ( v o l ) ;
83
APPENDICE A. CODICE DEL FILTRO ANISOTROPICO 3D
35
% PDE ( p a r t i a l d i f f e r e n t i a l
d i f f v o l = vol ;
clear vol
40
45
equation )
initial
condition .
% Center voxel d i s t a n c e s .
x = voxel spacing (1) ;
y = voxel spacing (2) ;
z = voxel spacing (3) ;
dx = 1 ;
dy = 1 ;
dz = z /x ;
dd = s q r t ( dxˆ2+dy ˆ 2 ) ;
dh = s q r t ( dxˆ2+dz ˆ 2 ) ;
dc = s q r t ( ddˆ2+dz ˆ 2 ) ;
50
55
60
65
70
75
80
% 3D
h1 =
h2 =
h3 =
h4 =
h5 =
h6 =
c o n v o l u t i o n masks − f i n i t e
z e r o s ( 3 , 3 , 3 ) ; h1 ( 2 , 2 , 2 ) =
z e r o s ( 3 , 3 , 3 ) ; h2 ( 2 , 2 , 2 ) =
z e r o s ( 3 , 3 , 3 ) ; h3 ( 2 , 2 , 2 ) =
z e r o s ( 3 , 3 , 3 ) ; h4 ( 2 , 2 , 2 ) =
z e r o s ( 3 , 3 , 3 ) ; h5 ( 2 , 2 , 2 ) =
z e r o s ( 3 , 3 , 3 ) ; h6 ( 2 , 2 , 2 ) =
differences .
−1; h1 ( 2 , 2 , 1 )
−1; h2 ( 2 , 2 , 3 )
−1; h3 ( 2 , 1 , 2 )
−1; h4 ( 2 , 3 , 2 )
−1; h5 ( 3 , 2 , 2 )
−1; h6 ( 1 , 2 , 2 )
=
=
=
=
=
=
1;
1;
1;
1;
1;
1;
=
=
=
,1)
,1)
,1)
,1)
,1)
1;
1;
1;
=
=
=
=
=
1;
1;
1;
1;
1;
h7 = z e r o s ( 3 , 3 , 3 ) ;
h8 = z e r o s ( 3 , 3 , 3 ) ;
h9 = z e r o s ( 3 , 3 , 3 ) ;
h10 = z e r o s ( 3 , 3 , 3 )
h11 = z e r o s ( 3 , 3 , 3 )
h12 = z e r o s ( 3 , 3 , 3 )
h13 = z e r o s ( 3 , 3 , 3 )
h14 = z e r o s ( 3 , 3 , 3 )
h7 ( 2 , 2 , 2 ) =
h8 ( 2 , 2 , 2 ) =
h9 ( 2 , 2 , 2 ) =
; h10 ( 2 , 2 , 2 )
; h11 ( 2 , 2 , 2 )
; h12 ( 2 , 2 , 2 )
; h13 ( 2 , 2 , 2 )
; h14 ( 2 , 2 , 2 )
−1; h7 ( 3 , 1 , 1 )
−1; h8 ( 2 , 1 , 1 )
−1; h9 ( 1 , 1 , 1 )
= −1; h10 ( 3 , 2
= −1; h11 ( 1 , 2
= −1; h12 ( 3 , 3
= −1; h13 ( 2 , 3
= −1; h14 ( 1 , 3
h15
h16
h17
h18
=
=
=
=
zeros
zeros
zeros
zeros
(3
(3
(3
(3
,3
,3
,3
,3
,3)
,3)
,3)
,3)
;
;
;
;
h15 ( 2
h16 ( 2
h17 ( 2
h18 ( 2
,2
,2
,2
,2
,2)
,2)
,2)
,2)
=
=
=
=
−1;
−1;
−1;
−1;
h15 ( 3
h16 ( 1
h17 ( 3
h18 ( 1
,1
,1
,3
,3
,2)
,2)
,2)
,2)
=
=
=
=
1;
1;
1;
1;
h19
h20
h21
h22
h23
h24
h25
h26
=
=
=
=
=
=
=
=
zeros
zeros
zeros
zeros
zeros
zeros
zeros
zeros
(3
(3
(3
(3
(3
(3
(3
(3
,3
,3
,3
,3
,3
,3
,3
,3
,3)
,3)
,3)
,3)
,3)
,3)
,3)
,3)
;
;
;
;
;
;
;
;
h19 ( 2
h20 ( 2
h21 ( 2
h22 ( 2
h23 ( 2
h24 ( 2
h25 ( 2
h26 ( 2
,2
,2
,2
,2
,2
,2
,2
,2
,2)
,2)
,2)
,2)
,2)
,2)
,2)
,2)
=
=
=
=
=
=
=
=
−1;
−1;
−1;
−1;
−1;
−1;
−1;
−1;
h19 ( 3
h20 ( 2
h21 ( 1
h22 ( 3
h23 ( 1
h24 ( 3
h25 ( 2
h26 ( 1
,1
,1
,1
,2
,2
,3
,3
,3
,3)
,3)
,3)
,3)
,3)
,3)
,3)
,3)
=
=
=
=
=
=
=
=
1;
1;
1;
1;
1;
1;
1;
1;
% Anisotropic diffusion .
f o r t = 1: num iter
85
90
95
100
105
% F i n i t e d i f f e r e n c e s . [ i m f i l t e r ( . , . , ’ conv ’ ) can be r e p l a c e d by convn
( . , . , ’ same ’ ) ]
% Due t o p o s s i b l e memory l i m i t a t i o n s , t h e d i f f u s i o n
% w i l l be c a l c u l a t e d a t e a c h page / s l i c e o f t h e volume .
f o r p = 1 : pags −2
d i f f 3 p p = d i f f v o l ( : , : , p : p+2) ;
aux = i m f i l t e r ( d i f f 3 p p , h1 , ’ c o n v ’ ) ; n a b l a 1 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h2 , ’ c o n v ’ ) ; n a b l a 2 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h3 , ’ c o n v ’ ) ; n a b l a 3 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h4 , ’ c o n v ’ ) ; n a b l a 4 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h5 , ’ c o n v ’ ) ; n a b l a 5 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h6 , ’ c o n v ’ ) ; n a b l a 6 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h7 , ’ c o n v ’ ) ; n a b l a 7 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h8 , ’ c o n v ’ ) ; n a b l a 8 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h9 , ’ c o n v ’ ) ; n a b l a 9 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h10 , ’ c o n v ’ ) ; n a b l a 1 0 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h11 , ’ c o n v ’ ) ; n a b l a 1 1 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h12 , ’ c o n v ’ ) ; n a b l a 1 2 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h13 , ’ c o n v ’ ) ; n a b l a 1 3 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h14 , ’ c o n v ’ ) ; n a b l a 1 4 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h15 , ’ c o n v ’ ) ; n a b l a 1 5 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h16 , ’ c o n v ’ ) ; n a b l a 1 6 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h17 , ’ c o n v ’ ) ; n a b l a 1 7 = aux ( : , : , 2 ) ;
aux = i m f i l t e r ( d i f f 3 p p , h18 , ’ c o n v ’ ) ; n a b l a 1 8 = aux ( : , : , 2 ) ;
84
APPENDICE A. CODICE DEL FILTRO ANISOTROPICO 3D
110
115
120
125
130
135
140
145
150
155
160
165
170
175
aux
aux
aux
aux
aux
aux
aux
aux
=
=
=
=
=
=
=
=
imfilter
imfilter
imfilter
imfilter
imfilter
imfilter
imfilter
imfilter
( diff3pp
( diff3pp
( diff3pp
( diff3pp
( diff3pp
( diff3pp
( diff3pp
( diff3pp
, h19 , ’ c o n v ’ ) ;
, h20 , ’ c o n v ’ ) ;
, h21 , ’ c o n v ’ ) ;
, h22 , ’ c o n v ’ ) ;
, h23 , ’ c o n v ’ ) ;
, h24 , ’ c o n v ’ ) ;
, h25 , ’ c o n v ’ ) ;
, h26 , ’ c o n v ’ ) ;
nabla19
nabla20
nabla21
nabla22
nabla23
nabla24
nabla25
nabla26
% Diffusion function .
i f o p t i o n == 1
c1 = exp ( −( n a b l a 1 / kappa ) . ˆ 2 ) ;
c2 = exp ( −( n a b l a 2 / kappa ) . ˆ 2 ) ;
c3 = exp ( −( n a b l a 3 / kappa ) . ˆ 2 ) ;
c4 = exp ( −( n a b l a 4 / kappa ) . ˆ 2 ) ;
c5 = exp ( −( n a b l a 5 / kappa ) . ˆ 2 ) ;
c6 = exp ( −( n a b l a 6 / kappa ) . ˆ 2 ) ;
c7 = exp ( −( n a b l a 7 / kappa ) . ˆ 2 ) ;
c8 = exp ( −( n a b l a 8 / kappa ) . ˆ 2 ) ;
c9 = exp ( −( n a b l a 9 / kappa ) . ˆ 2 ) ;
c 1 0 = exp ( −( n a b l a 1 0 / kappa ) . ˆ 2 ) ;
c 1 1 = exp ( −( n a b l a 1 1 / kappa ) . ˆ 2 ) ;
c 1 2 = exp ( −( n a b l a 1 2 / kappa ) . ˆ 2 ) ;
c 1 3 = exp ( −( n a b l a 1 3 / kappa ) . ˆ 2 ) ;
c 1 4 = exp ( −( n a b l a 1 4 / kappa ) . ˆ 2 ) ;
c 1 5 = exp ( −( n a b l a 1 5 / kappa ) . ˆ 2 ) ;
c 1 6 = exp ( −( n a b l a 1 6 / kappa ) . ˆ 2 ) ;
c 1 7 = exp ( −( n a b l a 1 7 / kappa ) . ˆ 2 ) ;
c 1 8 = exp ( −( n a b l a 1 8 / kappa ) . ˆ 2 ) ;
c 1 9 = exp ( −( n a b l a 1 9 / kappa ) . ˆ 2 ) ;
c 2 0 = exp ( −( n a b l a 2 0 / kappa ) . ˆ 2 ) ;
c 2 1 = exp ( −( n a b l a 2 1 / kappa ) . ˆ 2 ) ;
c 2 2 = exp ( −( n a b l a 2 2 / kappa ) . ˆ 2 ) ;
c 2 3 = exp ( −( n a b l a 2 3 / kappa ) . ˆ 2 ) ;
c 2 4 = exp ( −( n a b l a 2 4 / kappa ) . ˆ 2 ) ;
c 2 5 = exp ( −( n a b l a 2 5 / kappa ) . ˆ 2 ) ;
c 2 6 = exp ( −( n a b l a 2 6 / kappa ) . ˆ 2 ) ;
e l s e i f o p t i o n == 2
c1 = 1 . / ( 1 + ( n a b l a 1 / kappa ) . ˆ 2 ) ;
c2 = 1 . / ( 1 + ( n a b l a 2 / kappa ) . ˆ 2 ) ;
c3 = 1 . / ( 1 + ( n a b l a 3 / kappa ) . ˆ 2 ) ;
c4 = 1 . / ( 1 + ( n a b l a 4 / kappa ) . ˆ 2 ) ;
c5 = 1 . / ( 1 + ( n a b l a 5 / kappa ) . ˆ 2 ) ;
c6 = 1 . / ( 1 + ( n a b l a 6 / kappa ) . ˆ 2 ) ;
c7 = 1 . / ( 1 + ( n a b l a 7 / kappa ) . ˆ 2 ) ;
c8 = 1 . / ( 1 + ( n a b l a 8 / kappa ) . ˆ 2 ) ;
c9 = 1 . / ( 1 + ( n a b l a 9 / kappa ) . ˆ 2 ) ;
c 1 0 = 1 . / ( 1 + ( n a b l a 1 0 / kappa ) . ˆ 2 )
c 1 1 = 1 . / ( 1 + ( n a b l a 1 1 / kappa ) . ˆ 2 )
c 1 2 = 1 . / ( 1 + ( n a b l a 1 2 / kappa ) . ˆ 2 )
c 1 3 = 1 . / ( 1 + ( n a b l a 1 3 / kappa ) . ˆ 2 )
c 1 4 = 1 . / ( 1 + ( n a b l a 1 4 / kappa ) . ˆ 2 )
c 1 5 = 1 . / ( 1 + ( n a b l a 1 5 / kappa ) . ˆ 2 )
c 1 6 = 1 . / ( 1 + ( n a b l a 1 6 / kappa ) . ˆ 2 )
c 1 7 = 1 . / ( 1 + ( n a b l a 1 7 / kappa ) . ˆ 2 )
c 1 8 = 1 . / ( 1 + ( n a b l a 1 8 / kappa ) . ˆ 2 )
c 1 9 = 1 . / ( 1 + ( n a b l a 1 9 / kappa ) . ˆ 2 )
c 2 0 = 1 . / ( 1 + ( n a b l a 2 0 / kappa ) . ˆ 2 )
c 2 1 = 1 . / ( 1 + ( n a b l a 2 1 / kappa ) . ˆ 2 )
c 2 2 = 1 . / ( 1 + ( n a b l a 2 2 / kappa ) . ˆ 2 )
c 2 3 = 1 . / ( 1 + ( n a b l a 2 3 / kappa ) . ˆ 2 )
c 2 4 = 1 . / ( 1 + ( n a b l a 2 4 / kappa ) . ˆ 2 )
c 2 5 = 1 . / ( 1 + ( n a b l a 2 5 / kappa ) . ˆ 2 )
c 2 6 = 1 . / ( 1 + ( n a b l a 2 6 / kappa ) . ˆ 2 )
end
=
=
=
=
=
=
=
=
aux ( :
aux ( :
aux ( :
aux ( :
aux ( :
aux ( :
aux ( :
aux ( :
,:
,:
,:
,:
,:
,:
,:
,:
,2)
,2)
,2)
,2)
,2)
,2)
,2)
,2)
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
% D i s c r e t e PDE s o l u t i o n .
d i f f v o l ( : , : , p+1) = d i f f v o l ( : , : , p+1) + . . .
delta t ∗ ( . ..
( 1 / ( dz ˆ 2 ) ) ∗ c1 . ∗ n a b l a 1 + ( 1 / ( dz ˆ 2 ) ) ∗ c2 . ∗ n a b l a 2 +
...
( 1 / ( dx ˆ 2 ) ) ∗ c3 . ∗ n a b l a 3 + ( 1 / ( dx ˆ 2 ) ) ∗ c4 . ∗ n a b l a 4 +
...
85
APPENDICE A. CODICE DEL FILTRO ANISOTROPICO 3D
( 1 / ( dy ˆ 2 ) ) ∗ c5 . ∗ n a b l a 5 +
...
...
( 1 / ( dc ˆ 2 ) ) ∗ c7 . ∗ n a b l a 7 +
...
( 1 / ( dc ˆ 2 ) ) ∗ c9 . ∗ n a b l a 9 +
+ ...
( 1 / ( dh ˆ 2 ) ) ∗ c 1 1 . ∗ n a b l a 1 1
nabla12 + . . .
( 1 / ( dh ˆ 2 ) ) ∗ c 1 3 . ∗ n a b l a 1 3
nabla14 + . . .
...
( 1 / ( dd ˆ 2 ) ) ∗ c 1 5 . ∗ n a b l a 1 5
nabla16 + . . .
( 1 / ( dd ˆ 2 ) ) ∗ c 1 7 . ∗ n a b l a 1 7
nabla18 + . . .
...
( 1 / ( dc ˆ 2 ) ) ∗ c 1 9 . ∗ n a b l a 1 9
nabla20 + . . .
( 1 / ( dc ˆ 2 ) ) ∗ c 2 1 . ∗ n a b l a 2 1
nabla22 + . . .
( 1 / ( dh ˆ 2 ) ) ∗ c 2 3 . ∗ n a b l a 2 3
nabla24 + . . .
( 1 / ( dh ˆ 2 ) ) ∗ c 2 5 . ∗ n a b l a 2 5
nabla26 ) ;
180
185
190
( 1 / ( dy ˆ 2 ) ) ∗ c6 . ∗ n a b l a 6 +
( 1 / ( dh ˆ 2 ) ) ∗ c8 . ∗ n a b l a 8 +
( 1 / ( dh ˆ 2 ) ) ∗ c 1 0 . ∗ n a b l a 1 0
+ ( 1 / ( dc ˆ 2 ) ) ∗ c 1 2 . ∗
+ ( 1 / ( dc ˆ 2 ) ) ∗ c 1 4 . ∗
+ ( 1 / ( dd ˆ 2 ) ) ∗ c 1 6 . ∗
+ ( 1 / ( dd ˆ 2 ) ) ∗ c 1 8 . ∗
+ ( 1 / ( dh ˆ 2 ) ) ∗ c 2 0 . ∗
+ ( 1 / ( dh ˆ 2 ) ) ∗ c 2 2 . ∗
+ ( 1 / ( dc ˆ 2 ) ) ∗ c 2 4 . ∗
+ ( 1 / ( dc ˆ 2 ) ) ∗ c 2 6 . ∗
end
% I t e r a t i o n wa r ni n g .
f p r i n t f ( ’\ rIteration %d\n ’ , t ) ;
beep on
beep , p a u s e ( 0 . 2 ) , beep
beep o f f
195
200
end
205
% ”End o f Program ” w ar n in g .
beep on ;
beep , p a u s e ( 0 . 5 ) , beep , p a u s e ( 0 . 5 ) , beep , p a u s e ( 0 . 5 ) ,
beep , p a u s e ( 0 . 5 ) , beep , p a u s e ( 0 . 5 ) , beep , p a u s e ( 0 . 5 ) ;
86
Bibliografia
[1] Roger P.Woods: concetti basilari sulla registratura di immagini-UCLA
School of Medicine.
[2] Elastix: a toolbox for rigid and nonrigid registration of images,
http://elastix.isi.uu.nl/index.php
[3] Fsl: a toolbox for rigid and nonrigid registration of brain images,
http://www.fmrib.ox.ac.uk/fsl/index.html
[4] D. Rueckert, L. I. Sonoda, C. Hayes, D. L. G. Hill, M. O. Leach, and D.
J. Hawkes. Nonrigid registration using free-form deformations: Application
to breast MR images. IEEE Trans. Med. Imag., 18(8):712 – 721, 1999.
[5] Wikipedia: fondamenti della Risonanza Magnetica Nucleare
http://www.wikipedia.org/wiki/Risonanza_magnetica
[6] Wikipedia: Imaging a risonanza magnetica,
http://it.wikipedia.org/wiki/Imaging_a_risonanza_magnetica
[7] Digital Imaging in Communications and Medicine (DICOM) – Part1 : Introduction and Overview – National Electrical Manufacturers Association
Rosslyn, Virginia USA.
[8] 3DSlicer: a multi-platform, free and open source software package for
visualization and medical image computing, http://www.slicer.org/
[9] Manuale di Elastix, Stefan Klein and Marius Staring, January 6-2010.
[10] ITK: Insight Segmentation and Registration Toolkit, www.itk.org
[11] Decodifica di files DICOM e restauro di immagini di Risonanza Magnetica,
Salvatore La Bua, Calogero Crapanzano, Pietro Amato, Università degli
Studi di Palermo.
[12] J. Koenderink, The structure of images, Biol. Cybern., (1984) vol. 50, 363370.
[13] Fondazione IEO istituto europeo di oncologia,
http://www.ieo.it/italiano/index.asp
[14] Fondazione IRCCS istituto nazionale dei tumori,
http://www.istitutotumori.mi.it/default.asp
87
BIBLIOGRAFIA
[15] Materiale del corso di Elaborazione di immagini in medicina tenuto dal
prof.Andrea Giacchetti, Università degli Studi di Verona.
[16] Non-rigid Atlas-to-Image Registration by Minimization of ClassConditional Image Entropy, Emiliano D’Agostino, Frederik Maes,
Dirk Vandermeulen, Paul Suetens.
[17] Tissue-Based Affine Registration of Brain Images to form a Vascular Density Atlas, Derek Cool, Dini Chillet, Jisung Kim, Jean-Phillipe Guyon,
Mark Foskey, Stephen Aylward.
[18] Whole-brain functional magnetic resonance imaging mapping of acute nociceptive responses induced by formalin in rats using atlas registration-based
event-related analysis, Yen-Yu I. Shih, You-Yin Chen, Chiao-Chi V. Chen,
Jyh-Cheng Chen, Chen Chang, Fu-Shan Jaw.
88