Pubblichiamo questo contributo scritto da Davide Poggiali su una tematica che crediamo sia di grande interesse per i nostri lettori. Ringraziamo Davide per questo articolo e inseriamo qui di seguito una sua breve presentazione.
Davide Poggiali è laureato in Matematica e ha fatto il Dottorato in Neuroscience. Nella vita ha fatto (in ordine sparso) il babysitter, il docente, il portiere diurno, il barista. Attualmente è assegnista di ricerca al Padova Neuroscience Center. Suona chitarra, charango e ukelele.
La TAC raccontata da un matematico
Hai 5 anni e finalmente è arrivata domenica. Il celebre ovetto con sorpresa è già sulla credenza e tu lo aspetti con impazienza da tutta la settimana. Ti arrampichi per cercare di prenderlo, ma i tuoi ti trattengono.
“Dopo pranzo, se mangi tutto”.
Allora resti lì a guardarlo, magari lo metti davanti a una lampadina accesa, nella speranza di trovare un indizio: quale sorpresa ci sarà dentro? Come faccio a scoprirlo senza aprire l’ovetto, disobbediendo così ai miei?
Ma anche: è possibile almeno in teoria avere la sorpresa senza aprire l’uovo?
Oggi, con l’aiuto di alcune leggi fisiche e con un po’ di matematica, cercheremo di rispondere a queste domande.
1. Come faccio a scoprire che sorpresa c’è dentro l’ovetto senza aprirlo?
La risposta a questa domanda viene dalla fisica delle onde elettromagnetiche: la luce che ci permette di vedere l’ovetto non è infatti che una dei tanti tipi di onda elettromagnetica presente in natura. Il fatto che non riusciamo a vedere la sorpresa nasce dal fatto che pellicola e cioccolata riflettono tutte le onde di luce visibile. Per “vedere” la sorpresa occorre quindi utilizzare delle onde più penetranti, ossìa con una frequenza maggiore di qualche ordine di grandezza rispetto alla luce visibile.
1.1 I raggi X
Dobbiamo ora fare un salto nel tempo puntando il nostro orologio all’anno 1895, ossia alla scoperta e utilizzo dei raggi X da parte del fisico tedesco Wilhelm Conrad Röntgen.
La mano della moglie di Rontgen raggi X. Fonte: wikipedia.
Le lastre a raggi X ci hanno permesso di vedere all’interno oggetti che non possiamo o vogliamo aprire, primo fra tutti il corpo di esseri viventi.
1.2 La legge di Beer-Lambert
La legge fisica che regola il passaggio di un raggio elettromagnetico attraverso un mezzo fisico è la legge di Beer-Lambert. Chiamiamo $$I(x)$$ l’intensità di un raggio, ossia la potenza per unità di superficie, e $$A(x)$$ l’attenuazione della materia attraversata, ossia la probabilità che un raggio subisca nel punto $$x$$ una deviazione dal suo percorso naturale. La legge di Beer-Lambert dice che
$$dI = – A(x) I(x) dx$$
ossia che per percorrere un percorso lungo $$dx$$ l’intensità decresce in modo direttamente proporzionale all’attenuazione del mezzo, del raggio stesso e della lunghezza percorso.
Questa formula ci dà già un’idea di cosa sta succedendo nel caso del nostro ovetto: la carta di alluminio, della cioccolata e del contenitore di plastica in cui è chiusa la sorpresa non sono trasparenti. Quindi l’ attenuazione è alta e quasi tutti i raggi di luce visibile deviano, non permettendo al nostro occhio di vedere cosa c’è dentro. Non basta quindi mettere una lampadina dietro l’ovetto, occorre aumentare l’intensità del raggio che attraversa l’oggetto.
Passiamo quindi il nostro ovetto ai raggi X, otteniamo una lastra di questo tipo.
Un ovetto scansionato con una TAC. Fonte: youtube.
E scopro che dentro c’è una macchinina. Peccato solo che da bambino non abbia avuto l’accesso ad una sorgente di raggi X.
2. Riesco ad ottenere la sorpresa senza aprire l’uovo?
Ovviamente, no! Ma posso barare, stampandomi in 3D la sorpresa in barba ai miei (i quali ormai non si occupano più della mia educazione).
Problema: dall’esterno riesco solo a ruotare attorno all’ovetto e non ottengo nulla che io possa stampare in 3D, ma solo una serie di proiezioni.
Un ovetto scansionato con una TAC al variare dell’angolo di scansione (sx) e dell’altezza di scansione (dx). Fonte: youtube.
2.1 La trasformata di Radon
Per poter stampare in 3D la sorpresa devo poter ricostruire la mia scansione come se avessi tagliato l’oggetto a fette e fotografato ogni fetta. E qui interviene la matematica!
Passando all’integrale la legge di Beer-Lambert
$$dI = – A(x) I(x) dx$$
risulta
$$\int_{x_0}^{x_1}A(x)dx =-\int_{x_0}^{x_1}\frac{dI}{I}=-log\left(\frac{I(x_1)}{I(x_0)}\right).$$
Conosco l’intensità del raggio emesso dalla sorgente e posso misurare l’intensità del raggio uscente; riesco a ricostruire la mappa di attenuazione A(x)?
La risposta è no, i dati a nostra disposizione non sono sufficienti. Posso però ruotare attorno all’oggetto ed effettuare diverse scansioni, come nell’immagine precedente. Ho quindi a disposizione delle proiezioni al variare di un angolo di rotazione, da cui cerco di ricostruire l’attenuazione.
Formalizziamo un po’ il nostro problema: sia
$$f : \Omega \subseteq \mathbb{R}^2 \longrightarrow \mathbb{R}$$
allora le sue proiezioni sono rappresentate dalla trasformata di Radon, ossia la funzione
$$\displaystyle \mathcal{R} f (t,\theta) =\int_{\ell(t,\theta)} f =\int_{\mathbb{R}} R_\theta \begin{pmatrix}t \\ s\end{pmatrix} ds$$
dove $$R_\theta = \begin{pmatrix}\cos\theta& -\sin\theta \\ \sin\theta & \cos\theta\end{pmatrix} $$
è la matrice di rotazione in verso orario.
L’immagine ottenuta campionando $$\mathcal{R}f$$ è detta sinogramma.
Schema per la trasformata di Radon: per ogni angolo fissato si ha una funzione a singola variabile, ottenuta per integrazione.
Il nostro scopo è quindi, data la trasformata di una funzione $$\mathcal{R}f$$, ricostruire la funzione $$f$$, ossia invertire la trasformata di Radon. Si tratta quindi di un cosiddetto problema inverso, un problema di cui conosciamo il meccanismo che genera le nostre misure e di cui vogliamo stimare i dati che hanno generato tali misure.
Per testare i nostri algoritmi utilizzeremo un’immagine simulata detta phantom di cui calcoreremo il sinogramma. Dal sinogramma ricostruiremo poi l’immagine.
import matplotlib.pyplot as plt
%matplotlib inline
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale
image = shepp_logan_phantom()
image = rescale(image, scale=0.4, mode='reflect')
imshow(image, cmap=plt.cm.Greys_r)
Questo è il phantom di Shepp-Logan, uno dei più usati in simulazione. Si tratta semplicemente di una somma di funzioni caratteristiche di diverse ellissi; il risultato somiglia vagamente ad un cervello.
theta = np.linspace(0., 359., max(image.shape))
sinogram = radon(image, theta=theta, circle=True)
imshow(sinogram, cmap=plt.cm.Greys_r)
Se applichiamo la trasformata di Radon all’immagine precedente ecco i dati di proiezione simulati. Il sinogramma che vedete in figura ha caratteristiche in comune con il sinogramma dell’ovetto: entrambi sono infatti una sovrapposizione di “onde” che si intersecano.
Come viene calcolata la trasformata di Radon dalla funzione radon che abbiamo visto in azione?
Ci sono diversi modi per approssimare $$\mathcal{R}f$$, il modo più comune è ruotare l’immagine in senso orario di un angolo $$\theta$$ per poi stimare l’integrale “per colonne”.
Costruzione del sinogramma al variare dell’angolo $$\theta$$
2.2 La retroproiezione
Come fare ad invertire la trasformata di Radon? Andiamo per tentativi. Chiediamoci per ogni punto $$(x,y)$$ del piano quali sono le rette $$\ell(t,\theta)$$ che sono passate per quel punto. La risposta è semplice: sono le rette di equazione $$t = x\cos\theta+ y\sin\theta$$ al variare dell’angolo $$\theta$$. Proviamo a fare una media del sinogramma valutato in queste rette. Trattandosi di funzioni continue, faremo una media integrale.
Sia quindi $$g: S \times \mathbb{R} \longrightarrow \mathbb{R}$$ definiamo la sua Backprojection (BP) come
$$\mathcal{R}^* g =\frac{1}{|S|} \int_{S} g(x\cos\theta+ y\sin\theta, \theta) \;d\theta$$
dove $$S$$ è l’insieme in cui varia $$\theta$$, di solito $$[0,2\pi]$$ oppure $$[0,\pi]$$.
Purtroppo non basta; infatti $$\mathcal{R}^*$$ non è la trasformata inversa di $$\mathcal{R}$$, ma ci somiglia; è infatti il suo operatore aggiunto, un po’ come la trasposta $$A^t$$ non è in generale l’inversa di una matrice $$A$$, ma può essere utile per calcolarla.
reconstruction_BP = iradon(sinogram, theta=theta,
filter = None, circle=True)
imshow(reconstruction_BP, cmap=plt.cm.Greys_r)
Backprojection del sinogramma precedente.
Come vedete otteniamo un’immagine simile al sinogramma di partenza, ma troppo sfocata; alcune delle ellissi più piccole non sono facilmente visibili. Occorre quindi lavorare un po’ sulla BP per ottenere qualcosa in più.
Come viene calcolata la BP? Di solito si ruota il sinogramma calcolato per colonne (per fissare le rette verticali) per poi sommare il risultato, e ottenere la media.
2.3 Un po’ di teoremi
Osservando il rapporto fra $$\mathcal{R}$$ e la trasformata di Fourier $$\mathcal{F}$$, viene fuori questo risultato interessante
Teor 1. Central Slice Theorem $$\forall f$$ Radon e Fourier-trasformabile e $$\forall(S,\theta) \in S\times\mathbb{R}$$
$$\mathcal{F}_2 f(S\cos\theta, S\sin\theta)=\mathcal{F}(\mathcal{R}f)_{(S,\theta)}$$
Dove $$\mathcal{F}_2$$ è la trasformata di Fourier 2-dimensionale, la trasformata $$\mathcal{F}$$ è applicata al solo parametro lineare $$t$$ per ogni angolo $$\theta$$.
Da questo Teorema possiamo ricavare il secondo, che finalmente definisce l’inversa della trasformata di Radon.
Teor 2. FBP formula
$$ f(x,y) =\frac{1}{2} \mathcal{R}^* \left[ \mathcal{F}^{-1} \left(|S| \mathcal{F}(\mathcal{R}f){(S,\theta)} \right) \right]{(x,y)}$$
2.4 La retroproiezione filtrata
Abbiamo quindi trovato la Filtered Back Projection (FBP); il nome deriva dal fatto che si tratta della BP di un filtro in frequenza applicato al sinogramma; se infatti chiamassimo $$\tilde{S} = \mathcal{F}^{-1}(|S|)$$ risulterebbe che
$$|S| \mathcal{F}(\mathcal{R}f) = \mathcal{F}(\tilde{S}) \mathcal{F}(\mathcal{R}f) = \mathcal{F} \left( \tilde{S} \ast \mathcal{R}f \right)$$
ossia, sostituendo nella formula FBP
$$ f = \tilde{S} \ast \mathcal{R}f .$$
Essendo però $$\tilde{S} = \mathcal{F}^{-1}(|S|)$$ impossibile da calcolare numericamente (l’integrale diverge!) lo si sostituisce con un filtro passa-basso, dato che il rumore si concentra sulle alte frequenze; si usa quindi $$\tilde{S} = \mathcal{F}^{-1}(|S| \chi_{[-L,L]} )$$ con $$L>0$$ un parametro reale da tarare.
Filtro del nostro sinogramma.
Backprojection del sinogramma filtrato.
reconstruction_FBP = iradon(sinogram, theta=theta,
filter = 'shepp-logan', circle=True)
imshow(reconstruction_FBP, cmap=plt.cm.Greys_r)
Filtered Backprojection del sinogramma.
Come vedete usando la FBP riesco ad ottenere un’immagine più nitida, in cui posso distingere anche i particolari più piccoli.
3. Ora di merenda
Torniamo al tanto agognato ovetto. Fetta per fetta abbiamo un sinogramma, da cui riusciamo a ricostruire la fetta.
Una volta ricostruite tutte le fette ho ottenuto un’immagine 3D, da cui posso isolare i punti della sorpresa (con un po’ di pazienza per le rifiniture)
per ottenere così una mappa di punti che posso passare ad una stampante 3D
fatto! Finalmente ho in mano la sorpresa (o meglio una sua copia) e tutto questo con soltanto: una laurea in matematica, competenze di programmazione, una TAC, un computer abbastanza potente. Alla faccia dei miei!
Qualche riferimento (in ordine di apparizione)
[1] https://it.wikipedia.org/wiki/Spettroelettromagnetico
[2] https://it.wikipedia.org/wiki/RaggiX
[3] https://it.wikipedia.org/wiki/LeggediLambert-Beer
[4] T. Sauer, Big Data Challenges in Imaging from Computer Tomography https://alice.rm.cnr.it/ocs/index.php/MASCOT2018/MASCOT2018/paper/view/338
[5] https://scikit-image.org/docs/dev/autoexamples/transform/plotradon_transform.html
[6] T. Feeman, the mathematics of medical imaging, Springer https://www.springer.com/gp/book/9783319226644
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
Ancora nessun commento