Laboratorio di Strumentazione per Bioimmagini
A.A. 2010-2011
Esercizio 1: soluzione
Creiamo un’immagine fantoccio con la funzione “phantom” di Matlab:
% creo immagine
>> I = phantom();
Definiamo tre diversi campionamenti angolari: il primo campionando con 360 angoli tra 0 e 360
gradi, il secondo campionando ogni 5 gradi (73 angoli tra 0 e 360 gradi) e il terzo campionando
ogni 10 gradi (37 angoli tra 0 e 360 gradi):
% definisco diversi passi di campionamento angolare
>> angoli_1 = 0: 360/359 : 360; % campiono con 360 angoli
>> angoli_2 = 0: 5 : 360; % campiono ogni 5 gradi
>> angoli_3 = 0: 10 : 360; % campiono ogni 10 gradi
Calcoliamo ora i sinogrammi relativi ai campionamenti appena definiti, usando la funzione Matlab
“radon”:
% calcolo sinogrammi
>> s1 = radon(I, angoli_1);
>> s2 = radon(I, angoli_2);
>> s3 = radon(I, angoli_3);
Visualizziamo i tre diversi sinogrammi ottenuti:
>> figure(1)
>> subplot(131)
>> imagesc(s1)
>> title('Sinogramma con 360 angoli tra 0 e 360 gradi')
>> subplot(132)
>> imagesc(s2)
>> title('Sinogramma con 73 angoli tra 0 e 360 gradi')
>> subplot(133)
>> imagesc(s3)
>> title('Sinogramma con 37 angoli tra 0 e 360 gradi')
>> colormap(hot)
Come si può vedere al diminuire della frequenza angolare di campionamento diminuisce il numero
di proiezioni nel sinogramma.
Calcoliamo ora le immagini ricostruite, mediante retroproiezione filtrata, a partire dai tre diversi
sino grammi precedentemente ottenuti:
% calcolo immagine ricostruita
>> Irec1 = iradon(s1, angoli_1);
>> Irec2 = iradon(s2, angoli_2);
>> Irec3 = iradon(s3, angoli_3);
E le visualizziamo:
>> figure(2)
>> subplot(131)
>> imagesc(Irec1)
>> title('Immagine ricostruita dal sinogramma 1')
>> subplot(132)
>> imagesc(Irec2)
>> title('Immagine ricostruita dal sinogramma 2')
>> subplot(133)
>> imagesc(Irec3)
>> title('Immagine ricostruita dal sinogramma 3')
>> colormap(gray)
Appare evidente come al diminuire della frequenza di campionamento angolare peggiori la qualità
dell’immagine ricostruita (artefatti e scarsa resa di alcune regioni).
Esercizio 2: soluzione
La funzione che vogliamo creare dovrà ricevere in ingresso un singolo profilo proiettato e l’angolo
corrispondente di proiezione. Restituirà in uscita l’immagine ottenuta retroproiettando il profilo:
function img = retroproiezione_singola(s, theta)
%
%
%
%
%
%
%
%
%
%
Funzione che calcola la retroproiezione (non filtrata) di un singolo
profilo.
INPUT:
- s : singolo profilo
- theta: angolo di proiezione
OUTPUT:
- img : immagine ricostruita (retroproiezione)
La lunghezza del profilo rappresenta la massima lunghezza attraverso l’immagine, ossia la
lunghezza della sua diagonale, quindi a partire da questa possiamo ricavare la lunghezza del lato
dell’immagine che vogliamo ricostruire:
N = 2*floor( length(s)/(2*sqrt(2)) ); % lunghezza lato dell'immagine
Ora che conosciamo il lato dell’immagine allochiamo lo spazio in memoria per contenere la
retroproiezione e definiamo gli assi x e y dell’immagine ricostruita:
img = zeros(N); % allochiamo spazio in memoria per l'immagine
% definiamo gli assi x e y per l'immagine ricostruita
[x y] = meshgrid(1:N, N:-1:1);
x = x - floor(N/2); % spostiamo l’origine dell’asse x al centro dell’immagine
y = y - floor(N/2); % spostiamo l’origine dell’asse y al centro dell’immagine
Ora definiamo alcune variabili utili che ci permetteranno di calcolare direttamente la
retroproiezione:
costheta = cos(theta);
sintheta = sin(theta);
len = length(s); % lunghezza del profilo (proiezione)
ctrIdx = ceil(len/2); % indice del centro della proiezione
Infine retroproiettiamo il profilo:
% retroproiezione
proj = s;
% questa in iradon viene filtrata prima di essere
% retroproiettata("spalmata")
t = round(x*costheta + y*sintheta);
img = img + proj(t + ctrIdx); % retroproiezione del profilo
img = img / max(img(:)); % normalizza valori immagine tra zero e 1
La variabile proj contiene il profilo p(k) non filtrato. Nella retroproiezione filtrata, s viene appunto
filtrato prima di essere retroproiettato (“spalmato”) nell’immagine.
Nella matrice t, tutti i pixel (x, y) che stanno su una stessa retta x cos( )  y sin( )  r hanno valore
uguale a r. Ne risulta che nella matrice t c’è una mappa di tutte le rette che ricoprono il dominio
dell’immagine da ricostruire. Tutte le rette hanno pendenza tale da risultare perpendicolari alla
direzione di acquisizione della proiezione p(k). Vediamo sotto l’esempio di una matrice t calcolata
per un’immagine 200 x 200 ad un angolo di 45 gradi:
L’assegnazione img = img + proj(t + ctrIdx); replica il profilo proj fino a costruire
un’immagine 200 x 200. In ogni pixel (x, y) di tale immagine mette il valore assunto da proj nel suo
indice pari al valore assunto da t(x, y) + ctrIdx.
Infine l’immagine img viene normalizzata affinché assuma valori compresi tra 0 e 1.
Verifichiamo ora la correttezza della funzione, retroproiettando tre profili corrispondenti agli angoli
0, pi/4; pi/2 calcolati su un’immagine di lato 200 avente valore zero ovunque tranne in un quadrato
centrale di lato 100 con valore 1:
>> I = zeros(200);
>> I(51:150,51:150) = 1;
>> theta1 = 0;
>> theta2 = pi/4;
>> theta3 = pi/2;
>> p1 = radon(I, theta1*180/pi);
>> p2 = radon (I, theta2*180/pi);
>> p3 = radon (I, theta3*180/pi);
(NB: Volendo potremmo calcolare le singole proiezioni con la funzione creata nell’esercizio 2 dello
scorso laboratorio…)
>> img1 = retroproiezione_singola(p1, theta1);
>> img2 = retroproiezione_singola(p2, theta2);
>> img3 = retroproiezione_singola(p3, theta3);
E ora visualizziamo le immagini retroproiettate:
>>
>>
>>
>>
figure(1)
imagesc(img1)
colormap(gray)
title('Retroproiezione per il profilo con theta = 0')
>>
>>
>>
>>
figure(2)
imagesc(img2)
colormap(gray)
title('Retroproiezione per il profilo con theta = pi/4')
>>
>>
>>
>>
figure(3)
imagesc(img3)
colormap(gray)
title('Retroproiezione per il profilo con theta = pi/2')
Verifichiamo infine la correttezza di tali retroproiezioni confrontandole con le corrispondenti
ottenute usando la funzione iradon (che però esegue una retroproiezione filtrata). Nella
documentazione di Matlab si vede che la retroproiezione di un singolo profilo richiede l’utilizzo
della funzione iradon nel seguente modo:
>> img1_irad = iradon([p1 p1], [theta1*180/pi theta1*180/pi])/2;
>> img2_irad = iradon([p2 p2], [theta2*180/pi theta2*180/pi])/2;
>> img3_irad = iradon([p3 p3], [theta3*180/pi theta3*180/pi])/2;
>>
>>
>>
>>
figure(4)
imagesc(img1_irad);
colormap(gray)
title('Retroproiezione filtrata per il profilo con theta = 0')
>>
>>
>>
>>
figure(5)
imagesc(img2_irad);
colormap(gray)
title('Retroproiezione filtrata per il profilo con theta = pi/4')
>>
>>
>>
>>
figure(6)
imagesc(img3_irad);
colormap(gray)
title('Retroproiezione filtrata per il profilo con theta = pi/2')
Esercizio 3: soluzione
La funzione richiesta riceverà in ingresso il sinogramma s e il corrispondente passo di
campionamento angolare thetas usato per la sua acquisizione:
function img = retroproiezione(s, thetas)
%
%
%
%
%
%
%
%
%
%
%
Funzione che calcola la retroproiezione (non filtrata) di un'immagine
dato il sinogramma e il passo di campionamento angolare usato per
acquisirlo.
INPUT:
- s : sinogramma
- thetas: passo di campionamento angolare
OUTPUT:
- img : immagine ricostruita (retroproiezione)
D’ora in poi considereremo il sinogramma s come una matrice che ha come colonne i diversi profili
(proiezioni) ottenuti per diversi angoli theta.
Di fatto l’esercizio si risolve facilmente a partire dalla funzione scritta per l’esercizio 2, con minime
variazioni nella fase finale di retroproiezione (che ora coinvolge più profili anziché uno solo).
La prima parte riportata di seguito è del tutto analoga alla prima parte dell’esercizio 2 (si veda sopra
per un commento più dettagliato):
N = 2*floor( size(s,1)/(2*sqrt(2)) ); % lunghezza lato dell'immagine
img = zeros(N); % allochiamo spazio in memoria per l'immagine
% definiamo gli assi x e y per l'immagine ricostruita
[x y] = meshgrid(1:N, N:-1:1);
y = y - floor(N/2);
x = x - floor(N/2);
costheta = cos(thetas);
sintheta = sin(thetas);
len = size(s,1); % lunghezza di ogni proiezione
ctrIdx = ceil(len/2); % indice del centro delle proiezioni
Ora inseriremo la retroproiezione realizzata precedentemente per un singolo profilo, all’interno di
un ciclo for. In tal modo calcoleremo iterativamente un’immagine retroproiettata a partire da ogni
profilo contenuto nel sinogramma s e sommeremo (integrazione lungo theta) le diverse immagini:
% retroproiezione
for i=1:length(thetas)
proj = s(:,i); %ogni colonna del sinogramma è un profilo da retroproiettare
t = round(x*costheta(i) + y*sintheta(i));
img = img + proj(t + ctrIdx); % calcoliamo e sommiamo le immagini
end
img = img / max(img(:)); % normalizza valori immagine tra zero e 1
Testiamo ora la funzione creata. Calcoliamo la retroproiezione di un’immagine di lato 200 avente
valore zero ovunque tranne in un quadrato centrale di lato 100 con valore 1:
>> I = zeros(200);
>> I(51:150,51:150) = 1;
Calcoliamo il sinogramma s corrispondente ad un passo di campionamento angolare thetas usando
la funzione radon:
>> thetas = 0:pi/179:pi; % definiamo passo di campionamento angolare
>> s = radon(I, thetas*180/pi); % calcoliamo sinogramma
(NB: Volendo potremmo calcolare il sinogramma con la funzione creata nell’esercizio 3 dello
scorso laboratorio…)
Calcoliamo la retroproiezione usando la funzione retroproiezione appena sviluppata:
>> img = retroproiezione(s, thetas);
Calcoliamo la retroproiezione filtrata usando la funzione Matlab iradon:
>> img_IRAD = iradon(s, thetas*180/pi);
Confrontiamo le due immagini ottenute:
>>
>>
>>
>>
figure(1)
imagesc(img)
colormap(gray)
title('Immagine retroproiettata senza filtraggio')
>>
>>
>>
>>
figure(2)
imagesc(img_IRAD)
colormap(gray)
title('Retroproiezione filtrata ottenuta con iradon')