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^[1].

Ricostruzione tomografica - Proiezione, Back projection e Filtered back projection.

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^[2] 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 il 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^[3].
• Metodi algebrici → considerano la ricostruzione come il risultato di un problema di ottimizzazione ottenuto in maniera iterativa^[3]. 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^[1].

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^[4][5].

${\displaystyle f(x,y)}$ è l'oggetto sotto scansione, ${\displaystyle p_{\theta }(t)}$ una proiezione o funzione detector, mentre ${\displaystyle L(\theta ,t)}$ è rappresentata dalla linea passante tra ${\displaystyle A}$ e ${\displaystyle B}$

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

${\displaystyle p_{\theta }(t)=\int _{L(\theta ,t)}f(x,y)ds}$

Dove ${\displaystyle f(x,y)}$ è l’oggetto sotto scansione, mentre ${\displaystyle L(\theta ,t)}$ rappresenta la direzione di un raggio del fascio di scansione all’angolo ${\displaystyle \theta }$ che incide sul pixel t del detector, e può essere formalizzato come^[5]:

${\displaystyle L(\theta ,t)={(x,y)\in \mathbb {R} \times \mathbb {R} :x\cos \theta +y\sin \theta =t}}$

L’intero sinogramma può essere invece ottenuto applicando la trasformata di Radon^[5]:

${\displaystyle {\mathcal {R}}f(\theta ,t)=\int {L(\theta ,t)}f(x,y)ds}$

È possibile applicare al sinogramma la trasformata di Radon inversa per riottenere la funzione ${\displaystyle {\textit {f}}}$ 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 ${\displaystyle \theta }$ e punti del rivelatore ${\displaystyle t}$^[4][5][3].

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 ${\displaystyle {\textit {f}}}$, il teorema della fetta di Fourier afferma che:

• effettuare la proiezione di tale funzione su una linea monodimensionale con un angolo ${\displaystyle \theta }$ (come nella trasformata di Radon) e poi effettuare la trasformata di Fourier monodimensionale di tale proiezione;
• effettuare la trasformata di Fourier bidimensionale di ${\displaystyle {\textit {f}}}$ e poi sezionare il risultato lungo una linea con un angolo ${\displaystyle \theta }$ che attraversa l'origine;

sono operazioni equivalenti, che portano allo stesso risultato^[4][3][6].

Ricostruzione tramite trasformata di Fourier.

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 ${\displaystyle \theta }$ 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 ${\displaystyle \theta }$. A questo punto, la funzione ${\displaystyle {\textit {f}}}$ può essere recuperata applicando la trasformata di Fourier bidimensionale inversa^[4][3][6].

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^[6].

Algoritmo di retroproiezione

Ricostruzioni con diverso numero di proiezioni del fantoccio di Shepp-Logan. Algoritmo di retroproiezione.

Fantoccio di Shepp-Logan

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^[7].

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^[8]. Questa operazione può essere espressa dalla seguente equazione:

${\displaystyle f_{bp}(x,y)=\int p_{\theta }(\rho )d\theta =\int p_{\theta }(x\cos \theta +y\sin \theta )d\theta }$

dove ${\displaystyle f_{bp}(x,y)}$ rappresenta l'immagine di uscita, ${\displaystyle p_{\theta }}$ una proiezione a un certo angolo e ciascun punto ${\displaystyle (x,y)}$ rispetta l'equazione ${\displaystyle \rho =x\cos \theta +y\sin \theta }$^[9][10].

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^[11].

Algoritmo di retroproiezione filtrata

Ricostruzioni con diverso numero di proiezioni del fantoccio di Shepp-Logan. 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^[3].

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^[11].

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^[3][11].

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:

${\displaystyle {\begin{cases}f_{fbp}(x,y)=\int q_{\theta }(x\cos \theta +y\sin \theta )d\theta \\\\q_{\theta }(\rho )=\int P_{\theta }(\omega )H(\omega )e^{j2\pi \omega t}d\omega \end{cases}}}$

dove ${\displaystyle q_{\theta }}$ sono le proiezioni filtrate nel dominio delle frequenze ottenute dalla moltiplicazione del loro spettro ${\displaystyle P}$ con lo spettro del filtro ${\displaystyle H}$ ed ogni punto dell'immagine di uscita ${\displaystyle (x,y)}$ soddisfa l'equazione ${\displaystyle \rho =x\cos \theta +y\sin \theta }$.

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^[11].

Metodi algebrici

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

${\displaystyle Wv=p}$

dove ${\displaystyle v}$ è un vettore monodimensionale contenente tutti gli ${\displaystyle n}$ pixel del volume da ricostruire, ${\displaystyle p}$ è un vettore monodimensionale contenente tutti gli ${\displaystyle M}$ elementi del sinogramma (ciascun elemento rappresenta quindi il valore della misurazione effettuata da un pixel del detector in una determinata direzione di proiezione) e ${\displaystyle W}$ è una matrice bidimensionale di dimensione ${\displaystyle M\times n}$ chiamata matrice di proiezione^[3][12][13].

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^[13].

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

${\displaystyle p_{i}:=\sum _{j=0}^{n-1}w_{ij}v_{j}}$

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 ${\displaystyle W}$ 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^[13].

Questa formulazione matriciale è un modo per esprimere in forma di combinazione lineare ciò che accade quando un oggetto è sottoposto a tomografia ed è una formulazione equivalente a quella di una trasformata di Radon discreta: in entrambi i casi il risultato sarà il sinogramma del volume del soggetto sotto analisi. La forma inversa di questa formulazione algebrica è alla base degli approcci algebrici per ottenere una ricostruzione. Considerando tale formulazione inversa:

${\displaystyle u:=W^{T}p}$

e fissando come riferimento un singolo valore del vettore risultante ${\displaystyle u}$, l'operazione che si sta svolgendo è la seguente:

${\displaystyle u_{j}:=\sum _{i=0}^{m-1}w_{ij}p_{i}}$

che esprime come il valore associato a un certo pixel del volume sia una combinazione lineare di tutte le misurazioni correlate a quel pixel. Se la moltiplicazione della matrice ${\displaystyle W^{T}}$ con il vettore ${\displaystyle p}$ viene eseguita, viene svolta un'operazione equivalente a quella dell'algoritmo di retroproiezione, ottenendo una ricostruzione^[3][12][13].

Kernel di proiezione di Joseph (a), di striscia (b) e di linea (c).

I valori contenuti all'interno della matrice di proiezione dipendono dal setup di acquisizione e dalla sua geometria e, per la definizione di tali valori, è possibile ricorrere a diversi kernel di proiezione^[14]. Nel caso del kernel a linea, i pesi vengono definiti considerando la distanza percorsa da ogni raggio attraverso ciascun pixel del volume. Nel caso del kernel di Joseph, per ogni riga o colonna del volume del soggetto si considerano i due pixel più vicini dove il raggio interseca tale riga o colonna e i pesi vengono trovati mediante interpolazione lineare tra i pixel considerati. Un altro possibile kernel è il kernel a striscia, in cui i raggi non sono considerati come linee ma come strisce con una larghezza pari alla larghezza dei sensori del rivelatore e i pesi vengono trovati misurando l'area di ciascun pixel coperta dalla striscia^[13][14].

Ciascun kernel di proiezione presenta un compromesso tra precisione ed efficienza. Il kernel a striscia è quello che permette una rappresentazione più vicina e quindi più precisa a ciò che accade nella realtà, ma è anche il più complesso in termini di calcolo dei pesi. Per questo motivo, il kernel a linea viene di solito preferito^[13].

Tecnica di ricostruzione iterativa simultanea

La tecnica di ricostruzione iterativa simultanea (simultaneous iterative reconstruction technique o SIRT) è un approccio di ricostruzione algebrico e iterativo^[15][16] che, ad ogni passo, aggiorna la ricostruzione finale utilizzando contemporaneamente tutti i dati di proiezione. L'equazione che viene risolta ad ogni passo del ciclo di iterazione si può esprimere come^[17]:

${\displaystyle v^{(k+1)}=v^{(k)}+CW^{t}R(p-Wv^{(k)})}$

Ad ogni iterazione, il sinogramma della soluzione corrente viene generato e sottratto dal sinogramma di input dell'algoritmo: questo è espresso dalla porzione di espressione ${\displaystyle p-Wv^{(k)}}$, dove la proiezione diretta ${\displaystyle Wv^{(k)}}$ della soluzione corrente ${\displaystyle v^{(k)}}$ viene sottratta dai dati di proiezione disponibili ${\displaystyle p}$. In seguito, la differenza appena calcolata viene pesata tramite la matrice diagonale ${\displaystyle R}$, la quale contiene sulla diagonale principale un valore uguale all'inverso della somma di tutte le righe della matrice di proiezione ${\displaystyle W}$. Questa operazione assicura che i valori della differenza calcolata in precedenza siano ponderati da un valore che rappresenta l'inverso della lunghezza percorsa da ciascun raggio all'interno del volume, al fine di garantire la convergenza dell'algoritmo. Questa quantità viene poi retroproiettata mediante la pre-moltiplicazione con la matrice ${\displaystyle W^{T}}$ per riportare l'elemento calcolato nel dominio della ricostruzione e il risultato viene nuovamente pesato dalla matrice ${\displaystyle C}$, la quale contiene sulla diagonale un valore uguale all'inverso della somma di tutte le colonne della matrice di proiezione ${\displaystyle W}$. Infine il termine di aggiornamento appena calcolato viene aggiunto alla soluzione corrente, generando una nuova soluzione aggiornata^[17].

La soluzione iniziale utilizzata per inizializzare questo algoritmo è solitamente un'immagine nera, con valori pari a zero^[17][18].

Una ulteriore formulazione di questo algoritmo è la seguente:

${\displaystyle f_{j}^{(k+1)}=f_{j}^{(k)}+{\frac {1}{\sum _{i=1}^{M}\omega _{ij}}}\sum _{i=1}^{M}({\frac {\omega _{ij}(p_{i}-\sum _{h=1}^{N}\omega _{ih}x_{h}^{(k)})}{\sum _{h=1}^{N}\omega _{ih}}})\quad \forall j}$

In questa equazione, ${\displaystyle f_{j}^{(k)}}$ rappresenta la ricostruzione al passo ${\displaystyle k}$, ${\displaystyle \omega _{ij}}$ rappresenta il valore contenuto nella colonna ${\displaystyle i}$ e nella riga ${\displaystyle j}$ della matrice di proiezione ${\displaystyle W}$, ${\displaystyle p}$ sono i valori dei dati di proiezione e ${\displaystyle x}$ sono i valori contenuti nella ricostruzione^[18].

Tecnica di ricostruzione algebrica

Sequenza animata degli step di ricostruzione di ART, una iterazione.

La tecnica di ricostruzione algebrica (algebraic reconstruction technique o ART) è una tecnica di ricostruzione algebrica e iterativa il cui concetto è stato proposto inizialmente da Stefan Kaczmarz^[19] e successivamente adattato da Richard Gordon, Robert Bender e Gabor Herman^[20] per la ricostruzione di immagini^[21].

Questo algoritmo è in grado di calcolare una soluzione approssimata del sistema di equazioni lineari ${\displaystyle Wu=p}$ eseguendo iterativamente la formula:

${\displaystyle u^{k+1}=u^{k}+\lambda _{k}{\frac {p_{i}-\langle w_{i},u^{k}\rangle }{|w_{i}|^{2}}}w_{i}^{T}}$

dove ${\displaystyle a_{i}}$ è la riga ${\displaystyle i}$ della matrice ${\displaystyle W}$, ${\displaystyle u_{i}}$ è il componente numero ${\displaystyle i}$ del vettore ${\displaystyle u}$ e ${\displaystyle \lambda _{k}}$ è un parametro compreso tra ${\displaystyle 0}$ e ${\displaystyle 1}$ utilizzato per controllare la convergenza del sistema, dando la possibilità di controllare il bilanciamento tra qualità dell'output e tempo di calcolo^[20].

Tecnica di ricostruzione algebrica simultanea

La tecnica di ricostruzione algebrica simultanea (simultaneous algebraic reconstruction technique o SART) è un approccio di ricostruzione algebrico ed iterativo, evoluzione della tecnica di ricostruzione algebrica (ART), sviluppato da Anders Andersen e Avinash Kak nel 1984^[22].

Questo algoritmo è strettamente legato anche all'algoritmo SIRT. La differenza principale tra SIRT e ART consiste nel fatto che in SIRT il passo di iterazione contiene sia la proiezione diretta sia la retroproiezione su tutti i valori del rivelatore, mentre nell'approccio ART si utilizza solo un singolo valore del rivelatore per effettuare l'aggiornamento. Ciò implica che l'approccio di ART genera ricostruzioni molto più rapide, tuttavia la convergenza di questo metodo non è garantita e ciò può rappresentare un problema in caso di misure rumorose^[23].

L'algoritmo SART può essere considerato come un compromesso tra questi due metodi: in questo approccio, ad ogni iterazione vengono aggiornati solo i valori corrispondenti ai punti del dominio della ricostruzione la cui proiezione appartiene a una determinata direzione di proiezione, combinando così la velocità del metodo ART con la stabilità di convergenza di SIRT^[22][23].

Nell'algoritmo SART, l'ordine in cui le direzioni di proiezione vengono considerate può essere qualunque, tuttavia è importante considerare che, se vengono scelte due direzioni di proiezione vicine in due iterazioni consecutive, questo comporterà poca aggiunta di informazione alla ricostruzione che viene costruita iterativamente. Per questo motivo, è sconsigliabile considerare le proiezioni in ordine sequenziale, ma sarebbero meglio usarle in ordine casuale^[23].

Il passo di aggiornamento di SART può essere formalizzato come segue e, se confrontato con la formula esplicita di SIRT, si può notare come siano molto simili e che la formula di SART ha maggiori vincoli legati al valore del rivelatore da utilizzare.

${\displaystyle f_{j}^{(k+1)}=f_{j}^{(k)}+{\frac {1}{\sum _{i\in I_{\theta }}\omega _{ij}}}\sum _{i\in I_{\theta }}({\frac {\omega _{ij}(p_{i}-\sum _{h=1}^{N}\omega _{ih}u_{h}^{(k)})}{\sum _{h=1}^{N}\omega _{ih}}})\quad \forall j}$

In questa equazione ${\displaystyle f_{j}^{(k)}}$ rappresenta la ricostruzione al passo ${\displaystyle k}$, ${\displaystyle \omega _{ij}}$ rappresenta il valore contenuto nella colonna ${\displaystyle i}$ e nella riga ${\displaystyle j}$ della matrice di proiezione ${\displaystyle W}$, ${\displaystyle p}$ sono i valori dei dati di proiezione e ${\displaystyle x}$ sono i valori contenuti nella ricostruzione^[22][23].

Numerose varianti di SART sono state proposte negli anni nella letteratura scientifica, tra le quali: SARTF, FA-SART^[24], OS-SART, VW-OS-SART^[25].

Metodo del gradiente coniugato

Il metodo del gradiente coniugato ai minimi quadrati (conjugate gradient least squared o CGLS) è un algoritmo algebrico utilizzato per risolvere un generico sistema di equazioni del tipo ${\displaystyle Ax=b}$ e può essere dunque utilizzato per la ricostruzione tomografica, dato che risolve anche il sistema di equazioni che descrivono la relazione tra i dati di proiezione e l'immagine da ricostruire^[26].

Durante il processo di ricostruzione, CGLS cerca di minimizzare la differenza quadratica tra i dati di proiezione misurati e quelli stimati iterativamente. Ad ogni iterazione, l'algoritmo calcola una direzione di ricerca coniugata che minimizza l'errore residuo. Successivamente, questa direzione di ricerca viene utilizzata per aggiornare l'immagine stimata e CGLS continua a iterare fino a quando l'immagine converge verso una soluzione approssimata. Questo algoritmo è noto per la sua velocità e capacità di gestire grandi quantità di dati, ma può essere sensibile al rumore nei dati di proiezione^[26][27].

CGLS è un algoritmo derivato dal più generale metodo del gradiente coniugato, ma altre varianti sono state proposte per risolvere il problema tomografico come ad esempio MCGT che applica il metodo del gradiente coniugato a un sistema di Tikhonov^[28].

Metodi basati sull'apprendimento profondo

L'influenza del rumore di Poisson nella ricostruzione fa sì che la rete U-net non riesca a ricostruire un oggetto simile a una lesione ad alto contrasto.

I metodi di apprendimento profondo sono ampiamente utilizzati per la ricostruzione di immagini e sono stati applicati in diversi ambiti, tra cui la riduzione del rumore nelle tomografie a bassa dose di radiazioni, la ricostruzione con numero sparso di viste, la tomografia ad angolo limitato e la riduzione degli artefatti dovuti alla presenza di metalli. Una panoramica eccellente può essere trovata nell'edizione speciale della rivista IEEE Transactions on Medical Imaging^[29]. Ad esempio, una delle maggiori tipologie di algoritmi di ricostruzione tramite deep learning applica reti neurali di post-processing per migliorare una ricostruzione già esistente, ottenuta in precedenza tramite metodi di ricostruzione convenzionali. Un esempio di applicazione è la riduzione degli artefatti utilizzando la rete U-Net nella tomografia ad angolo limitato^[30].

In un'immagine ricostruita mediante un metodo completamente basato sui dati usando come input il sinogramma, possono tuttavia verificarsi degli errori. Pertanto, l'integrazione di operatori noti nella progettazione dell'architettura delle reti neurali potrebbe risultare vantaggiosa, come descritto nel concetto di apprendimento di precisione^[31]. Ad esempio, la ricostruzione diretta di immagini dai dati di proiezione può essere appresa dalla struttura dell'algoritmo di retroproiezione filtrata^[32]. Un altro esempio è la costruzione di reti neurali mediante l'iterazione di algoritmi di ricostruzione^[33]. Oltre all'apprendimento di precisione, un approccio alternativo per migliorare la qualità delle ricostruzioni ottenute tramite deep learning è utilizzare metodi di ricostruzione convenzionali applicati a una precedente ricostruzione basata su deep learning^[34].

Software di ricostruzione tomografica

Per la ricostruzione tomografica sono disponibili toolbox open source come PYRO-NN^[35], TomoPy^[36], CONRAD^[37], ODL^[38], ASTRA toolbox^[39][40] e TIGRE^[41].

TomoPy è una toolbox open source in Python progettata presso l'Argonne National Laboratory appositamente per l'elaborazione di dati tomografici e le operazioni di ricostruzione di immagini. TomoPy include diversi algoritmi di ricostruzione, che possono essere eseguiti su workstation multicore e strutture di calcolo di grandi dimensioni^[42].

CONRAD è una toolbox open source specializzata nella tomografia a fascio conico. Offre una vasta gamma di funzionalità per l'acquisizione e l'elaborazione di immagini a raggi X, compresa la ricostruzione tomografica^[37].

ODL (Operator Discretization Library) è una toolbox open source per l'apprendimento automatico, la ricostruzione tomografica e altre applicazioni di elaborazione delle immagini. Supporta diverse geometrie di acquisizione, algoritmi di ricostruzione e metodi di regolarizzazione^[38].

TIGRE è una toolbox MATLAB sviluppata congiuntamente dall'Università di Bath e dal CERN per la ricostruzione di immagini CBCT. Offre una vasta gamma di algoritmi di ricostruzione e sfrutta l'accelerazione GPU per ridurre i tempi di calcolo^[41].

ASTRA toolbox è una toolbox MATLAB e Python di primitive GPU ad alte prestazioni per la tomografia 2D e 3D, sviluppato dal 2009 al 2014 dal laboratorio di visione iMinds-Vision Lab dell'Università di Anversa e successivamente sviluppato congiuntamente da imec-VisionLab (ex iMinds-VisionLab) e CWI di Amsterdam. Questa toolbox supporta la tomografia parallela, a ventaglio e a fascio conico oltre a permettere un posizionamento altamente flessibile della sorgente e del rivelatore. Nel 2016, ASTRA è stato integrato nel framework TomoPy^[43]. Questa integrazione ha permesso agli utenti di ASTRA di sfruttare le funzionalità di TomoPy per il filtraggio dei dati e la correzione degli artefatti, mentre ha permesso agli utenti di TomoPy di usufruire di un maggior numero di algoritmi di ricostruzione. Gli algoritmi di ricostruzione disponibili tramite l'integrazione tra TomoPy e ASTRA comprendono FBP, Gridrec, ART, SIRT, SART, BART, CGLS, PML, MLEM e OSEM^[44].

Note

1. (EN) Gabor Herman, Fundamentals of computerized tomography: Image reconstruction from projection, Springer, 2009, pp. 27-36.
2. (EN) Paolo Russo, Handbook of X-ray Imaging Physics and Technology, CRC Press, 2018, p. 638.
3. (EN) Paolo Russo, Handbook of X-ray Imaging Physics and Technology, CRC Press, 2018, pp. 669-709.
4. (EN) Gabor Herman, Fundamentals of computerized tomography: Image reconstruction from projection, Springer, 2009, pp. 101-124.
5. (EN) ASTRA Toolbox, Analytical projection. URL consultato il 6 luglio 2023.
6. (EN) ASTRA Toolbox, Fourier Slice Theorem. URL consultato il 6 luglio 2023.
7. (EN) Stanley Deans, The Radon Transform and Some of Its Applications, New York, John Wiley & Sons, 1983, pp. 11-16.
8. (EN) D. E. Dudgeon e R. M. Mersereau, Multidimensional digital signal processing, Prentice-Hall, 1984, pp. 363-378.
9. (EN) Computed Tomography, su www.tristanvanleeuwen.github.io. URL consultato il 5 luglio 2023.
10. (EN) Computed Tomography and ASTRA Toolbox training course, su www.visielab.uantwerpen.be. URL consultato il 5 luglio 2023.
11. (EN) ASTRA Toolbox, Filtered Backprojection (FBP). URL consultato il 6 luglio 2023.
12. (EN) A. Kak e M. Slaney, Principles of Computerized Tomographic Imaging, New York, IEEE Press, 1999, pp. 276–277, 284.
13. (EN) ASTRA Toolbox, Algebraic Tomography. URL consultato il 6 luglio 2023.
14. (EN) Fu, Jian, Zhenzhong, Liu e Wang, Jingzheng, Multi-Mounted X-Ray Computed Tomography, in PloS one, vol. 11, aprile 2016, pp. 1-15, DOI:10.1371/journal.pone.0153406.
15. (EN) M. Beister, D. Kolditz e W. A. Kalender, Iterative reconstruction methods in x-ray ct, in Physica Medica, vol. 288, n. 2, 2012, pp. 94-108, DOI:10.1016/j.ejmp.2012.01.003.
16. (EN) L. L. Geyer, U. J. Schoepf e F. G. Meinel, State of the art: Iterative ct reconstruction techniques., in Radiology, vol. 276, n. 2, 2016, pp. 339-357, PMID 26203706.
17. (EN) Tessa Van, Sarah Wuyts e Maggie Goossens, ITERATIVE RECONSTRUCTION ALGORITHMS The implementation of iterative reconstruction algorithms in MATLAB., luglio 2007, pp. 1-9.
18. (EN) ASTRA Toolbox, Simultaneous Iterative Reconstruction Techniquse (SIRT). URL consultato il 6 luglio 2023.
19. (DE) S.Kaczmarz, Angenäherte Auflösung Von Systemen Linearer Gleichungen, in Bulletin del’Académie Polonaise des Sciences et Lettres, A35, 1937, pp. 355-357.
20. (EN) Richard Gordon, Robert Bender e Gabor T. Herman, Algebraic Reconstruction Techniques (ART) for three-dimensional electron microscopy and X-ray photography, in Journal of Theoretical Biology, vol. 29, n. 3, dicembre 1970, pp. 471-476, DOI:10.1016/0022-5193(70)90109-8.
21. (EN) Touraj Nikazad, Algebraic Reconstruction Methods (PDF), su diva-portal.org, 5-8.
22. (EN) A.H. Andersen e A.C. Kak, Simultaneous Algebraic Reconstruction Technique (SART): A superior implementation of the ART algorithm, in Ultrasonic Imaging, vol. 6, n. 1, gennaio 1984, pp. 81-94.
23. (EN) ASTRA Toolbox, Other Algebraic Reconstruction Methods. URL consultato il 6 luglio 2023.
24. Ziying Zhou, Yugang Li, Fa Zhang e Xiaohua Wan, FASART: An iterative reconstruction algorithm with inter-iteration adaptive NAD filter, in io-medical materials and engineering, vol. 26, Suppl 1, pp. S1409–S1415, DOI:10.3233/BME-151439.
25. (EN) Pan Jinxiao, Tie Zhou e Yan Han, Variable Weighted Ordered Subset Image Reconstruction Algorithm, in International Journal of Biomedical Imaging, vol. 2006, 2006, pp. 1–7, DOI:10.1155/IJBI/2006/10398.
26. (EN) M. Soleimani e T. Pengpen, Introduction: a brief overview of iterative algorithms in X-ray computed tomography, in Philosophical transactions. Series A, Mathematical, physical, and engineering sciences, vol. 373, n. 2043, 2015, pp. 1-6, DOI:10.1098/rsta.2014.0399.
27. (EN) Karl Henrik May, Andreas Keil e Georg Von Freymann, The Conjugate Gradient Least Square Algorithm in Terahertz Tomography, in IEEE Access, settembre 2021, pp. 14216 -142178, DOI:10.1109/ACCESS.2021.3116801.
28. (EN) Qi Wang e Huaxiang Wang, A modified conjugate gradient method based on the Tikhonov system for computerized tomography (CT), in ISA Transactions, vol. 50, n. 2, 2011, pp. 256-261, DOI:10.1016/j.isatra.2010.11.001.
29. (EN) Ge Wang, Jong Chu Ye, Klaus Mueller e Jeffrey A. Fessler, Image reconstruction is a new frontier of machine learning, in IEEE Transactions on Medical Imaging, vol. 37, n. 6, 2018, pp. 1289–1296, DOI:10.1109/TMI.2018.2833635.
30. (EN) Jawook Gu e Jong Chul Ye, Multi-scale wavelet domain residual learning for limited-angle CT reconstruction, Fully3D, 2017, pp. 443–447.
31. (EN) Andreas K. Maier, Christopher Syben, Bernhard Stimpel, Tobias Wuerfl, Mathis Hoffman, Frank Schebesch, Weilin Fu, Leonid Mill, Lasse Kling e Silke Christiansen, Learning with known operators reduces maximum error bounds, in Nature Machine Intelligence, vol. 1, n. 8, 2019, pp. 373–380, DOI:10.1038/s42256-019-0077-5, PMID 31406960.
32. (EN) Tobias Wuerfl, Mathis Hoffmann, Vincent Christlein, Katharina Breininger, Yixing Huang, Mathias Unberath e Andreas Maier, Deep Learning Computed Tomography: Learning Projection-Domain Weights from Image Domain in Limited Angle Problems, in IEEE Transactions on Medical Imaging, vol. 37, n. 6, 2018, pp. 1454–1463, DOI:10.1109/TMI.2018.2833499, PMID 29870373.
33. (EN) J. Adler e O. Öktem, Learned Primal-Dual Reconstruction, in IEEE Transactions on Medical Imaging, vol. 37, n. 6, 2018, pp. 1322–1332, DOI:10.1109/TMI.2018.2799231, PMID 29870362.
34. (EN) Y. Huang, A. Preuhs, G. Lauritsch, M. Manhart, X. Huang e A. Maier, Data Consistent Artifact Reduction for Limited Angle Tomography with Deep Learning Prior, Machine Learning for Medical Image Reconstruction, 2019, pp. 101–112, DOI:10.1007/978-3-030-33843-5_10.
35. (EN) C. Syben, M. Michen, B. Stimpel, S. Seitz, S. Ploner e A. Maier, PYRO-NN: Python Reconstruction Operators in Neural Networks, in Medical Physics, vol. 46, n. 11, 2019, pp. 5110–5115, DOI:10.1002/mp.13753, PMID 31389023.
36. (EN) D. Gursoy D, F. De Carlo, X. Xiao X e C. Jacobsen C, TomoPy: A framework for the analysis of synchrotron tomographic data, in Journal of Synchrotron Radiation, vol. 22, n. 5, 2014, pp. 1188–1193, DOI:10.1107/S1600577514013939, PMID 25178011.
37. (EN) A. Maier, H. G. Hofmann, M. Berger, P. Fischer, C. Schwemmer, H. Wu, K. Mueller, J. Hornegger, J. Choi, C. Riess, A. Keil e A. Farhig, CONRAD - A software framework for cone-beam imaging in radiology, in Medical Physics, vol. 40, n. 11, 2013, p. 111914, DOI:10.1118/1.4824926, PMID 24320447.
38. (EN) The ODL Contributors, Operator Discretization Library Documentation, su https://odlgroup.github.io. URL consultato il 7 luglio 2023.
39. (EN) Van Aarle, W., Palenstijn, W.J., De Beenhouwer, J., Altantzis T., Bals S., Batenburg K. J., and J. Sijbers, The ASTRA Toolbox: a platform for advanced algorithm development in electron tomography, in Ultramicroscopy, vol. 157, October 2015, pp. 35-47, DOI:10.1016/j.ultramic.2015.05.002, PMID 26057688.
40. (EN) W. Van Aarle, W J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. De Beenhouwer, K. J. Batenburg, and J. Sijbers, Fast and flexible X-ray tomography using the ASTRA toolbox, in Optics Express, vol. 24, n. 22, 2016, pp. 25129-25147, DOI:10.1364/OE.24.025129, PMID 27828452.
41. (EN) Ander Biguri, Manjit Dosanjh, Steven Hancock e Manuchehr Soleimani, TIGRE: a MATLAB-GPU toolbox for CBCT image reconstruction, in Biomedical Physics & Engineering Express, vol. 2, n. 5, 8 settembre 2016, p. 055010, DOI:10.1088/2057-1976/2/5/055010.
42. (EN) Bicer T., Gursoy D., Kettimuthu R., De Carlo F., e Foster I., Optimization of tomographic reconstruction workflows on geographically distributed resources, in Journal of Synchrotron Radiation, vol. 23, n. 4, 2016, pp. 997–1005, DOI:10.1107/S1600577516007980, PMID 27359149.
43. (EN) Pelt D.M., Gursoy D., Batenburg K.J., De Carlo F., Palenstijna W.J., and Sijbers J., Integration of TomoPy and the ASTRA toolbox for advanced processing and reconstruction of tomographic synchrotron data, in Journal of Synchrotron Radiation, vol. 23, n. 3, 2016, pp. 842–849, DOI:10.1107/S1600577516005658, PMID 27140167.
44. (EN) ASTRA Toolbox, ASTRA Toolbox Docs, su https://www.astra-toolbox.com/. URL consultato il 6 luglio 2023.

Bibliografia

• L. Faggioni, F. Paolicchi, E. Neri, Elementi di tomografia computerizzata, Springer, 2010, ISBN 8847016967.
• (EN) P. Russo, Handbook of X-ray Imaging Physics and Technology, CRC Press, 2018, ISBN 9781498741521.
• (EN) Gabor Herman, Fundamentals of computerized tomography: Image reconstruction from projection, Springer, 2009, ISBN 978-1852336172.
• (EN) Stanley Deans, The Radon Transform and Some of Its Applications, New York, John Wiley & Sons, 1983, ISBN 978-0486462417.

Voci correlate

• Raggi X
• Tomografia computerizzata
• Trasformata di Radon

Collegamenti esterni

• (EN) Vision Lab - University of Antwerp, Computed Tomography and ASTRA Toolbox training course, su visielab.uantwerpen.be. URL consultato il 6 luglio 2023.
• (EN) Features — TomoPy fa2d648c269267ff7ea8f5ee6f219722cde58bc3 documentation, su tomopy.readthedocs.io. URL consultato il 6 luglio 2023.
• (EN) The ASTRA Toolbox — ASTRA Toolbox 2.1.0 documentation, su www.astra-toolbox.com. URL consultato il 6 luglio 2023.
• (EN) TIGRE toolbox, su tomopedia.github.io. URL consultato il 6 luglio 2023.
• (EN) A. C. Kak and Malcolm Slaney, Principles of Computerized Tomographic Imaging, su slaney.org. URL consultato il 6 Luglio 2023.

Portale Matematica

Portale Medicina