Algoritmi di ricostruzione tomografica

Gli algoritmi di ricostruzione tomografica sono un insieme di metodi matematici utilizzati per ottenere un'immagine della sezione o una ricostruzione tridimensionale di un soggetto a cui è stata effettuata una tomografia.

Questi algoritmi rielaborano i dati di proiezione contenuti nel sinogramma del soggetto sottoposto a scansione, che è anche il risultato della trasformata di Radon del soggetto stesso. Tali dati sono costituiti da una serie di proiezioni effettuate da diverse angolazioni attorno al soggetto da analizzare e vengono acquisiti misurando l'attenuazione dei raggi X (nel caso della TAC) mentre questi attraversano il soggetto in diverse direzioni. L'obiettivo di questi algoritmi è quindi quello di determinare l'attenuazione dei raggi X in ogni punto dello spazio di ricostruzione sotto analisi al fine di estrarne un'immagine.

Classificazione Gli algoritmi (o metodi) di ricostruzione tomografica possono essere divisi principalmente in due categorie:

Metodi analitici → algoritmi basati su considerazioni matematiche riguardo al problema, in particolare sull'antitrasformata di Radon. I metodi analitici sono poco onerosi in termini di calcolo, ma sono poco robusti in scenari caratterizzati da una bassa emissione di radiazioni o in cui il rapporto segnale rumore è basso. Metodi algebrici → considerano la ricostruzione come il risultato di un problema di ottimizzazione ottenuto in maniera iterativa. Questi ultimi metodi permettono lo sfruttamento di eventuali informazioni conosciute a priori riguardo l'oggetto da ricostruire ma, in virtù della loro natura iterativa e d'ottimizzazione, sono generalmente più impegnativi da un punto di vista computazionale rispetto ai metodi analitici, sebbene il calcolo possa essere parallelizzato e accelerato con l'ausilio di una GPU.

Metodi analitici Gli approcci analitici sono basati su considerazioni matematiche effettuate sul principio di funzionamento della tomografia, ovvero l’operazione di proiezione. Nel caso bidimensionale, l'obiettivo è estrarre una fetta (o slice) del soggetto sotto analisi corrispondente al piano su cui si trovano i raggi, sottesi tra la sorgente e il rivelatore (o detector). Per farlo, il soggetto è posto solitamente sull’asse di rotazione di sorgente e detector e questi catturano diverse proiezioni del soggetto da diverse angolazioni intorno ad esso.

Alla fine di questa procedura, i dati raccolti possono essere riorganizzati sovrapponendo ciascuna proiezione all'altra per creare il sinogramma, che si presenta come un'immagine in cui le righe corrispondono alle proiezioni ottenute a diversi angoli e le colonne corrispondono ai sensori che compongono il detector. Il nome sinogramma deriva dal fatto che il suo contenuto ricorda l'unione di diverse onde sinusoidali: supponendo che il soggetto della scansione sia un punto (quindi una funzione delta di Dirac in 2D), il risultato nel sinogramma sarebbe una singola onda sinusoidale, con un'ampiezza crescente all’allontanarsi dall'asse di rotazione. Il sinogramma può essere quindi considerato una combinazione di molteplici onde sinusoidali poiché il soggetto della scansione, se discretizzato, può essere visto come una combinazione di funzioni delta di Dirac, e dal momento che il sinogramma è anche il risultato della trasformata di Radon di un oggetto e dal momento che tale trasformazione è lineare, la combinazione di sinusoidi che compone il sinogramma è lineare.

Una singola proiezione contenuta nel sinogramma può essere formalizzata come:

Dove è l’oggetto sotto scansione, mentre rappresenta la direzione di un raggio del fascio di scansione all’angolo che incide sul pixel t del detector, e può essere formalizzato come:

L’intero sinogramma può essere invece ottenuto applicando la trasformata di Radon:

È possibile applicare al sinogramma la trasformata di Radon inversa per riottenere la funzione e quindi la ricostruzione, tuttavia bisogna tenere presente che ciò funzionerebbe perfettamente solo se la trasformata di Radon, o la scansione CT, fossero effettuate in maniera continua e non discreta, con un numero infinito di angoli e punti del rivelatore .

Un altro elemento matematico è costituito dal teorema della fetta di Fourier, che collega tra loro la trasformata di Radon e la trasformata di Fourier bidimensionale. Dato una funzione bidimensionale , il teorema della fetta di Fourier afferma che: effettuare la proiezione di tale funzione su una linea monodimensionale con un angolo (come nella trasformata di Radon) e poi effettuare la trasformata di Fourier monodimensionale di tale proiezione; effettuare la trasformata di Fourier bidimensionale di e poi sezionare il risultato lungo una linea con un angolo che attraversa l'origine; sono operazioni equivalenti, che portano allo stesso risultato.

Questo teorema evidenzia la relazione che esiste tra la trasformata di Fourier bidimensionale e la trasformata di Radon: conoscendo le proiezioni per tutti gli angoli intorno a un oggetto, è possibile infatti ottenere la sua trasformata di Fourier bidimensionale semplicemente effettuando la trasformata di Fourier unidimensionale di tutte le proiezioni e poi organizzandole in uno spazio bidimensionale come linee che passano attraverso l'origine inclinate dell’angolo . A questo punto, la funzione può essere recuperata applicando la trasformata di Fourier bidimensionale inversa.

Dal momento che i dati di proiezione a diversi angoli sono ciò che produce una tomografia, questo algoritmo può essere usato per ottenere una ricostruzione. Tuttavia, sussistono dei problemi con l'utilizzo di questo approccio. Innanzitutto, la quantità di misurazioni ottenute da una tomografia è finita, il che porta ad ottenere un dominio di Fourier della ricostruzione che è campionato in maniera non uniforme: i campioni della trasformata di Fourier prodotti da questo approccio si trovano infatti distribuiti lungo cerchi concentrici invece che seguire una griglia, una caratteristica richiesta per algoritmi come la trasformata di Fourier veloce per ottenere l'antitrasformata di Fourier. Il secondo problema risiede nel fatto che la distribuzione campionaria del dominio di Fourier è molto più densa vicino all'origine rispetto alle regioni esterne: ciò significa che le frequenze più basse sono sovra rappresentate rispetto alle frequenze più alte e quando verrà applicata la trasformata di Fourier inversa l’immagine finale sarà sfocata.

Algoritmo di retroproiezione

L'algoritmo di retroproiezione (o algoritmo di back-projection) è un metodo di ricostruzione analitico equivalente all'operazione inversa di proiezione (o antitrasformata di Radon), che nella tomografia computerizzata avviene quando i raggi X proiettano il soggetto sotto scansione sui diversi pixel del detector.

In questo algoritmo, ciascun punto dell'immagine di uscita riceve un valore proporzionale alla proiezione sul detector del corrispondente punto del dominio dell'oggetto scansionato, per tutti i diversi angoli di proiezione. Per ottenere questo risultato, ciascuna proiezione viene propagata sull'immagine d'uscita lungo una direzione in accordo con il proprio angolo di acquisizione, così che diversi valori provenienti da diversi angoli di acquisizione si sommino nei punti dell'immagine producendo il risultato desiderato. Questa operazione può essere espressa dalla seguente equazione:

dove rappresenta l'immagine di uscita, una proiezione a un certo angolo e ciascun punto rispetta l'equazione .

Sebbene l'algoritmo di retroproiezione richieda una capacità di calcolo limitata per produrre un risultato, il suo svantaggio risiede nella produzione di immagini affette da sfocatura. Se si prova ad applicare l'algoritmo ad uno stesso soggetto di partenza ma si effettuano ricostruzioni con un numero crescente di proiezioni, si potrà notare come l'immagine di uscita con più proiezioni risulti più delineata rispetto a quelle con meno proiezioni, tuttavia ancora sfocata rispetto al soggetto scansionato. La ragione di questo effetto è dovuta ad una sovra rappresentazione delle basse frequenze rispetto alle alte frequenze nei dati di proiezione, e può essere risolto filtrando i dati come avviene nell'algoritmo di retroproiezione filtrata.

Algoritmo di retroproiezione filtrata

L’algoritmo di retroproiezione filtrata (o di filtered back-projection) è un metodo analitico ed un’evoluzione dell’algoritmo di retroproiezione che mitiga gli effetti di sfocatura.

Operativamente, è diviso in due parti: fase di filtraggio e fase di retroproiezione. La prima è una fase di pre-elaborazione dei dati, che serve a preparare i dati prima che vengano rielaborati dalla seconda, ed è equivalente all’algoritmo di retroproiezione. La prima fase infatti consiste nel filtrare tutte le proiezioni presenti nel sinogramma al fine di attenuare le basse frequenze. In questo modo, quando i dati saranno ricostruiti nella seconda fase, l’immagine finale non apparirà sfocata.

Per effettuare questa operazione possono essere impiegati diversi tipi di filtri passa alto, ma spesso è utilizzato il filtro di Ram-Lak che è in grado di attenuare le basse frequenze, amplificare le alte e tagliare tutte le frequenze oltre una certa soglia definita come frequenza di taglio (o di cut-off). Il filtro viene applicato nel dominio delle frequenze. Per prima cosa, viene applicata la trasformata di Fourier alle proiezioni, poi il risultato di questa operazione viene moltiplicato per il filtro effettuando quindi il filtraggio. Infine, il risultato viene riportato nel dominio iniziale tramite l’antitrasformata di Fourier.

Una volta effettuato il filtraggio, i dati vengono retroproiettati seguendo i passi dell’algoritmo di retroproiezione. Questo metodo può essere formalizzato con le seguenti equazioni:

dove sono le proiezioni filtrate nel dominio delle frequenze ottenute dalla moltiplicazione del loro spettro con lo spettro del filtro ed ogni punto dell'immagine di uscita soddisfa l'equazione .

Provando ad applicare questo algoritmo ad uno stesso soggetto, ma effettuando ricostruzioni con un numero crescente di proiezioni e comparando i risultati dell'esperimento con quelli ottenuti utilizzando l’algoritmo di retroproiezione, si potrà notare come, a parità di numero di proiezioni utilizzate, il risultato dell’algoritmo di retroproiezione filtrata appare più definito e meno sfocato. Ciò è dovuto al fatto che le basse frequenze non sono più sovra-rappresentate rispetto alle alte poiché soppresse dal filtro.

Metodi algebrici Gli approcci algebrici si basano su un sistema di equazioni in forma matriciale che rappresenta l'operazione di proiezione. Questo sistema può essere scritto come:

dove è un vettore monodimensionale contenente tutti gli pixel del volume da ricostruire, è un vettore monodimensionale contenente tutti gli elementi del sinogramma (ciascun elemento rappresenta quindi il valore della misurazione effettuata da un pixel del detector in una determinata direzione di proiezione) e è una matrice bidimensionale di dimensione chiamata matrice di proiezione.

Il ruolo della matrice di proiezione è quello di mettere in relazione le misurazioni o proiezioni effettuate con il volume da ricostruire. Il contenuto della matrice può quindi essere interpretato come la quantità con cui ciascun pixel del volume scansionato contribuisce a ciascuna singola misurazione contenuta in ogni proiezione.

Quando l'operazione espressa dall'equazione riportata viene eseguita è come se, considerando un singolo valore del vettore risultante , l'operazione che si sta svolgendo sia:

che è una combinazione lineare di tutti i pixel del volume del soggetto, ciascuno pesato da un valore che esprime quanto contribuisce alla misurazione. Poiché una singola riga della matrice corrisponde fisicamente a un raggio che attraversa il volume sotto analisi e quindi non tutti i pixel del volume entreranno in contatto con esso, la matrice di proiezione risulta essere estremamente sparsa.