Università degli Studi di Napoli – “Parthenope”
Dipartimento di Scienze e Tecnologie
Corso di Laurea in Informatica
TESI DI LAUREA TRIENNALE
Relatore: Dott. Livia Marcellino
Laureando: Ivan Osato (0124/250)
Anno Accademico 2013/2014







La struttura delle immagini DWI
Il metodo PCA locale overcomplete
Introduzione al GP-GPU
Ambiente di sviluppo CUDA
Il software parallelo sviluppato
Test e risultati
Conclusioni e Sviluppi futuri
Ivan Osato - Università degli Studi di Napoli "Parthenope"
2
Le immagini DWI (Diffusion Weighted Images) sono immagini di risonanza magnetica
(MRI) pesate in diffusione, capaci di misurare il moto microscopico delle molecole di
acqua in un tessuto biologico, permettendo quindi l’analisi della microstruttura bianca
cerebrale sia in condizioni normali che patologiche.

La diffusione rappresenta il moto caotico e disordinato
delle molecole di un mezzo, dovuto all’agitazione termica.

L’elemento preponderante è l’acqua che comprende il 6590% in volume dei tessuti biologici e svolge la funzione di
mezzo di trasporto dei composti biochimici.

Strutture con diffusione normale sono rappresentate
convenzionalmente più scure (segnale di risonanza più
attenuato) mentre laddove la velocità di diffusione è
minore (ischemia) sono rappresentate più chiare.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
3

Il più recente sviluppo della DWI è rappresentato dall’imaging in tensore di
diffusione (Diffusion Tensor Imaging, DTI), che è una tecnica che mette in
evidenza non soltanto l’entità della diffusione delle molecole di acqua nei tessuti
ma anche la direzione consentendo di superare i limiti di direzionalità di DWI.

Le immagini DTI si ottengono attraverso l’acquisizione di immagini pesate in
diffusione lungo diverse direzioni tramite l’applicazione di una sequenza EPI
(Echo Planar Imaging), cui viene fatta seguire in successione una serie di
gradienti con orientazioni variabili.
slice (3D)
direzioni (4D)
Ivan Osato - Università degli Studi di Napoli "Parthenope"
4

Le immagini DWI hanno per natura, un basso SNR (Signal to Noise Ratio) a
causa della scarsa ampiezza del segnale ed il pronunciato rumore termico.

Questo basso rapporto rende l’analisi estremamente complessa, pregiudicando
la stima dei parametri di diffusione quantitativi.
Esistono diverse tecniche di denoising
che possono migliorare la qualità dei
dati DWI come step di postelaborazione
Ivan Osato - Università degli Studi di Napoli "Parthenope"
5
Il metodo PCA locale per la rimozione del rumore Rician da immagini DWI, sfrutta la
ridondanza multi-direzionale delle immagini pesate in diffusione multicomponente,
per ridurre localmente le componenti principali meno significative attraverso un
approccio overcomplete.
L’elaborazione di ogni singolo voxel avviene localmente, all’interno di finestre locali
che possono sovrapporsi a quelle dei voxel adiacenti. In seguito alla sovrapposizione
delle patch (overlapping), elaborate separatamente, si ottengono diverse stime per
alcuni voxel da processare attraverso una media (overcomplete).
Ivan Osato - Università degli Studi di Napoli "Parthenope"
6
La Principal Component Analysis (PCA, Analisi per Componenti Principali) è una
trasformazione lineare ortogonale dei dati che vengono mappati in un nuovo sistema di
riferimento ortogonale in modo da massimizzare la varianza (rispetto ad ogni altra proiezione)
associata alla prima coordinata, poi quella associata alla seconda coordinata e così via. Le
proiezioni, considerate nell’ordine, prendono il nome di componenti principali.
.
La PCA consiste nel calcolo degli autovalori e degli autovettori della matrice di covarianza XTX
dei dati centrati:


.
Il primo autovettore di XTX è la prima componente principale, mentre il primo autovalore è
la varianza associata a tale componente principale, che è la massima possibile.
La varianza che rimane, è associata alle altre direzioni e quindi alle restanti componenti
principali.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
7
Per ogni punto xp con p = (i, j, z, k) nel dominio dell’immagine 4D, le patch 3D
circostanti xp in ogni direzione K, sono riordinate come vettore colonna di una
matrice X di (N x K) componenti (N elementi della patch x K direzioni).
PER OGNI MATRICE X
PCA LOCALE





C = XTX
Calcolo autovettori di C in W
Y = XW
Thresholding di Y in Ŷ
con τ = γ * σ
XD = ŶW-1 + μ
OVERLAPPING DI XD
Le colonne della matrice sono riposizionate
nell’immagine di output (sommate ai valori già
presenti), nella posizione dove erano state estratte
come patch 3D dall’immagine rumorosa.
MEDIA
Ivan Osato - Università degli Studi di Napoli "Parthenope"
8
Ivan Osato - Università degli Studi di Napoli "Parthenope"
9
L’unità di elaborazione grafica (Graphics Processing Unit, GPU) è una tipologia
particolare di coprocessore composto da più multiprocessori paralleli, organizzati
secondo uno schema architetturale finalizzato al paradigma SIMD (Single Instruction
Multiple Data).
Queste unità nascono nell’ambito della Computer
Graphics per l’esecuzione di operazioni grafiche
(rendering) che lavorano su una grande quantità di
dati.
La GP-GPU (General Purpose Computing on GPU),
si propone di impiegare le GPU per applicazioni
parallele High Performance General Purpose,
sfruttando le loro capacità di elaborazione in parallelo
per compiti che vanno oltre la Computer Graphics.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
10
Come usare un ambiente ideato per la grafica in
applicazioni General Purpose?
NVIDIA ha modificato la GPU rendendola interamente
programmabile
per
le
applicazioni
scentifiche
e
aggiungendovi il supporto per linguaggi ad alto livello come C
e C++ nell’ambiente di sviluppo CUDA.
Il modello architetturale CUDA considera
la CPU e la GPU come due macchine
distinte dette host
e device, che
eseguono rispettivamente la parte di
codice sequenziale e parallela (kernel).
Ivan Osato - Università degli Studi di Napoli "Parthenope"
11
L’algoritmo parallelo sviluppato, prevede una griglia computazionale 3D in
cui ad ogni thread di un blocco è associato un voxel di competenza (i, j, z, k)
in tutte le K direzioni dell’immagine DWI.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
12
Versione 1 - full shared memory: uso totale di memoria condivisa in tutte le principali
operazioni di PCA locale.
VANTAGGI
SVANTAGGI
1. Bassa latenza
2. Risparmio maggiore di tempo
Capienza ridotta
KERNEL
• computeCovariance<<<-;-;->>>
• computePCA<<<-;-;->>>
• weightedAverage<<<-;-;->>>
Versione 2 - shared and global memory combination: gestione oculata della shared
memory in combinazione all’utilizzo della global memory.
VANTAGGI
SVANTAGGI
Capienza maggiore
Tempi d’accesso più alti
KERNEL
• computeCovariance<<<-;-;->>>
• computeProdYEXW<<<-;-;->>>
• computeProdXDEYWT<<<-;-;->>>
• weightedAverage <<<-;-;->>>
Ivan Osato - Università degli Studi di Napoli "Parthenope"
13
Questo kernel ha il compito di calcolare la matrice di covarianza C per ogni voxel (i, j, z).
La configurazione di thread utilizzata è la seguente:
BlockDim =
GridDim
=
Ogni thread del blocco, in parallelo, copia i dati di propria competenza dall’immagine
rumorosa per costruire la matrice X in shared memory, poi calcola la media di ciascuna
colonna della propria X (per normalizzare X) e per finire esegue il prodotto righe per
colonne per il calcolo della matrice di covarianza C.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
14

Calcola la Principal Component
Analysis in maniera overcomplete
usando la shared memory.

La configurazione di griglia è la
medesima del kernel precedente.

Agisce
in
global
memory
sull’immagine di output per
eseguire l’overlapping.
Nella fase di inizializzazione, i thread del blocco in parallelo costruiscono la matrice
X, il vettore di media Mean e gli autovettori W in shared memory. Successivamente,
eseguono tutte le operazioni di PCA locale e memorizzano i dati restaurati
nell’immagine di output (overlapping).
Ivan Osato - Università degli Studi di Napoli "Parthenope"
15
Nella versione 2, il kernel computePCA<<<-;-;->>> è splittato in due kernel che
eseguono separatamente i prodotti per occupare meno spazio in shared memory.

computeProdYEXW<<<-;-;->>> Ciascun thread di un blocco, in parallelo,
prima effettua il prodotto tra matrici per generare Y, poi esegue il thresholding
per la rimozione del rumore. Il prodotto fa un uso combinato di shared memory
(X, W) e global memory (Y).

computeProdXDEYWT<<<-;-;->>> Ciascun thread di un blocco, in parallelo,
prima esegue il prodotto tra matrici per la determinazione della matrice XD
restaurata, poi effettua overlapping per collocare i dati restaurati nell’immagine
di output. Il prodotto fa un uso combinato di shared memory (MEAN, W) e
global memory (Y, XD).
Ivan Osato - Università degli Studi di Napoli "Parthenope"
16

Calcola la media pesata per ogni
voxel dell’immagine di output
restaurata.

La configurazione di griglia è la
medesima dei kernel precedenti.

Non fa uso di shared memory.
Ogni thread di un blocco accede in memoria globale ai suoi K voxel di
competenza ed esegue la media dividendo il valore del voxel (calcolato
come somma delle stime di voxel dal processo di overlapping) per il numero
di sovrapposizioni.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
17
In questo primo test, si riportano i risultati numerici dell’errore di approssimazione
commesso sulle immagine restaurate rispetto alle immagini originali (ground-truth).
L’errore viene calcolato attraverso l’indice PSNR (Peak Signal-Noise Ratio) che è
in grado di valutare la qualità delle immagini.
ground-truth
(z = 0, k = 0)
ground-truth
(z = 0, k = 5)
ground-truth
(z = 0, k = 9)
Simulated DW-MRI Brain Data Sets: http://www.nitrc.org/projects/sim_dwi_brain/
Ivan Osato - Università degli Studi di Napoli "Parthenope"
18
Si riportano i valori di PSNR per le tre immagini campione utilizzate, al variare del
livello di rumore.
Immagine K = 5
Immagine K = 0
Livello
Rumore
Fattore
Threshold
Truth/Noise
Truth/Denoise
Livello
Rumore
Fattore
Threshold
Truth/Noise
Truth/Denoise
3%
3,5
39,3627
40,4992
3%
3,5
40,0504
46,4654
5%
4
35,1357
36,0129
5%
4
36,5721
43,7716
7%
5,5
32,5834
33,3899
7%
5,5
33,606
42,1731
9%
5,5
30,676
32,8739
9%
5,5
31,471
40,4865
Immagine K = 9
I valori tabellari dimostrano
che l’algoritmo è in grado di
produrre risultati di output
con un discreto grado di
accuratezza.
Livello
Rumore
Fattore
Threshold
Truth/Noise
Truth/Denoise
3%
3,5
40,2151
46,8808
5%
4
36,6531
44,0942
7%
5,5
33,7361
42,1145
9%
5,5
31,5565
40,5904
Aumentando il livello di rumore,
il PSNR tende a diminuire
poichè
l’immagine
è
più
degradata.
I valori (Truth/Denoise) sono più
alti rispetto a (Truth/Noise) dato
che l’immagine è restaurata ed
ha una qualità maggiore.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
19
Risultati di output che si ottengono dall’esecuzione dell’algoritmo per un z fissato (z = 0),
esaminando le direzioni k = 0, 5, 9 con i livelli di rumore 3% e 5%.
LIVELLO RUMORE 3%
LIVELLO RUMORE 5%
Ivan Osato - Università degli Studi di Napoli "Parthenope"
20
nvprof --print-gpu-trace ./myApp
Lo strumento di profiling nvprof consente di raccogliere e visualizzare i dati relativi alle
caratteristiche di esecuzione di un software parallelo eseguito su GPU in ambiente CUDA.
Duration
Grid size
Block size
DSMem
Kernel Name
1.10s
(88, 88, 7)
(2, 2, 3)
8.48 KB
computeCovariance
7.21s
(88, 88, 7)
(2, 2, 3)
34.40 KB
computePCA
7.32ms
(88, 88, 7)
(2, 2, 3)
0 KB
weightedAverage
Output del tool, eseguendo la
seconda versione dell’algoritmo
parallelo su 10 direzioni, con la
configurazione di (2, 2, 3) thread
per blocco e (88, 88, 7) blocchi per
un confronto diretto con la prima
versione.
Output del tool, eseguendo la prima
versione dell’algoritmo parallelo su 10
direzioni, con l’unica configurazione
possibile di (2, 2, 3) thread per blocco
e (88, 88, 7) blocchi.
Duration
Grid size
Block size
DSMem
Kernel Name
1.10s
(88, 88, 7)
(2, 2, 3)
8.48 KB
computeCovariance
1.04s
(88, 88, 7)
(2, 2, 3)
8.48 KB
computeProdYEXW
5.42s
(88, 88, 7)
(2, 2, 3)
8.48 KB
computeProdXDEYWT
7.27ms
(88, 88, 7)
(2, 2, 3)
0 KB
weightedAverage
La differenza sostanziale tra le due versioni dell’algoritmo parallelo, sta nel numero di kernel
utilizzati per ottenere il risultato, nella quantità di shared memory utilizzata e nel numero di
copie dati da host a device e viceversa.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
21
Sono stati effettuati vari test, variando per
ogni esecuzione il numero di direzioni
dell’immagine DWI (partendo da 2 direzioni
ed arrivando fino a un massimo di 10
direzioni entro i limiti di memoria della
macchina utilizzata).
Dataset size
(voxel)
GFlops CPU
GFlops GPU
1300992
0,03
0,33
1951488
0,05
0,63
2601984
0,13
1,6
3252480
0,31
3,8
3902976
0,7
8,16
4553472
1,41
17
5203968
2,6
31
5854464
4,3
51,51
6504960
7
83,14
Per 10 direzioni, le performance della CPU sono
pari a 7 GFlops mentre per la GPU sono 83,14
GFlops, ovvero circa 12 volte maggiori rispetto a
quelle della CPU ottenendo un guadagno
prestazionale pari a 91,6%.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
22
Nei diversi test, si è variato il numero di direzioni dell’immagine DWI (numero di voxel
processati). Le due versioni parallele sono state eseguite con la stessa configurazione
confrontabile (per motivi di spazio) di (2, 2, 3) thread per blocco e (88, 88, 7) blocchi di griglia.
Dataset
size
(voxel)
Time
CPU
1300992
17
1,04
1,35
1951488
25
1,72
2,03
2601984
31
2,82
3,1
3252480
39
4,23
4,61
3902976
49
6,6
6,7
4553472
57
9,1
9,74
5203968
68
15,01
12,8
5854464
82
18,4
17,4
6504960
95
23,55
24
Time
Time
GPU v1 GPU v2
Con 10 direzioni, il tempo impiegato dalla CPU è di 95 secondi
contro i 24 secondi della GPU (i tempi delle due versioni sono
approssimativamente uguali). Pertanto, il tempo impiegato dalla
GPU è circa 4 volte minore rispetto al tempo impiegato dalla CPU,
con un guadagno pari al 75%.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
23
Thread
per
blocco
Direzioni
Time
GPU
v2
Time
GPU
v1
2x2x3
2
1,35
1,04
3x3x3
2
1,43
1,1
4x4x3
2
1,45
1,04
5x5x3
2
1,5
1,27
2x2x3
10
24
23,55
3x3x3
10
31,8
-
4x4x3
10
35,6
-
5x5x3
10
29,30
-
L’algoritmo parallelo è stato eseguito al variare
delle direzioni, processando i voxel dell’immagine
con diverse configurazioni di griglia, nei limiti dello
spazio a disposizione sulla macchina utilizzata.
Aumentando le direzioni, il carico di lavoro
aumenta a prescindere dalla configurazione
utilizzata. Per ogni direzione, la configurazione
ottimale è sempre quella che utilizza 2 x 2 x 3
thread per blocco.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
24
CONCLUSIONI
Con il lavoro svolto in questa tesi, si è realizzata una implementazione parallela del
metodo “PCA Locale Overcomplete (OLPCA)” per la rimozione del rumore Rician
da immagini pesate in diffusione, in ambiente GPU-CUDA. Sono state sviluppate
due versioni parallele, per sopperire alle limitazioni dell’architettura utilizzata. I
risultati sperimentali, di confronto con l’algoritmo sequenziale, hanno evidenziato un
guadagno in termini di GFlops addirittura pari al 91,6% ed una riduzione del tempo di
esecuzione uguale al 75%.
SVILUPPI FUTURI


Implementazione multi-GPU per incrementare ulteriormente il guadagno.
Testing del software parallelo su schede video con memorie più capienti.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
25



Thread per
blocco
Time GPU
v1
Time GPU
v2
2x2x3
1,04
1,35
3x3x3
1,1
1,43
4x4x3
1,04
1,45
5x5x3
1,27
1,5
La differenza nei tempi d’esecuzione tra le due versioni è decisamente bassa.
I tempi della seconda versione parallela, sono influenzati dalla quantità maggiore di
trasferimenti dati in memoria e dall’accesso alla global memory.
I risultati confermano ancora una volta che la configurazione ottimale è certamente quella
che prevede l’uso di blocchi da 2 x 2 x 3 thread.
Ivan Osato - Università degli Studi di Napoli "Parthenope"
28
L’algoritmo è stato eseguito al variare delle
direzioni, rilevando per ogni test il tempo di
esecuzione necessario alla computazione
degli autovettori su host.
Dataset size
(voxel)
Direzioni
Time W
1300992
2
0,9
1951488
3
1,31
2601984
4
2
3252480
5
3,1
3902976
6
4,45
4553472
7
6,4
5203968
8
8,25
5854464
9
11,30
6504960
10
15
Il tempo d’esecuzione aumenta all’aumentare delle direzioni dell’immagine DWI con un
andamento quadratico!
Ivan Osato - Università degli Studi di Napoli "Parthenope"
29