Laboratorio di Neuroingegneria

Segmentazione

A.A. 2012-2013

Cosa faremo oggi

• Caricamento di volumi MRI cerebrale– T1– FLAIR

• Calcolo della soglia di Otsu

• Stima della mistura di Gaussiane

CaricamentoSono a disposizione due file NIFTI: T1.nii e FLAIR.nii

Dopo averli caricati, dal momento che le unità di misura dei dati MRI sono arbitrari, conviene:

- estrarre una fetta campione su cui lavorare

- normalizzare i dati:

% Estraggo i dati dell’immagine dalla struttura NIFTIxT1=T1.img;

% Estraggo una fetta campionexim=xT1(:,:,70);

% Normalizzo di datiI=xim/max(xim(:));

Es1: Sogliatura

Scegliere un valore casuale T0 compreso fra il minimo ed il massimo valore dell’immagine

Calcolare l’immagine binaria I>T

Calcolare la media m0 dei valori contenuti nella regione per cui I<T ed m1 della regione I>T

Es 1: SogliaturaProvare l’algoritmo seguente:

tt

t

t

t

TTfinché

mmT

TIImeanm

TIImeanm

1

101

1

0

2

)|(

)|( % Inizializzo T. . .

while (ec)T=Tnew;m1=mean(I(I<T));m2=mean(I(I>=T));Tnew=0.5*(m1+m2)ec=(T~=Tnew))

end;

Es 2: Soglia di Otsu

Calcolare la soglia di Otsu e confrontarla la segmentazione con quella ottenuta con l’esercizio precedente.

Suggerimenti:

- Creare un insieme di soglie da provare campionando le intensità dell’immagine

- Per ogni Ti calcolare la varianza interclasse

- Trovare la soglia corrispondente alla massima varianza interclasse

Es 3: Istogramma

Calcolare l’istogramma dell’immagine con la funzione hist() di Matlab e disegnarci sopra le soglie trovate con i due algoritmi precedenti

Suggerimento:

La funzione hist di Matlab può:

- visualizzare l’istogramma con un numero definito di campioni

- restituire i valori dell’istogramma ed i campioni nei quali è stato calcolato

% Possibili sintassi di hist

hist(I)hist(I,N)[y,x]=hist(I,N)y=hist(I,x)

Es 4: mistura di Gaussiane 1

1) Creare una funzione per calcolare la pdf

2) Creare una funzione per calcolare la somma di N pdf

3) Fit non lineare

4) Plottare l’istogramma e la mistura di gaussiane calcolata

5) Segmentare l’immagine

Es 4: mistura di Gaussiane % Funzione per calcolare una gaussiana unidimensionale% m,s parametri scalarifunction y=Gauss(x,m,s). . .

% Funzione per calcolare la mistura di gaussiane% w,m,s parametri vettoriali (Nx1)function MG=GaussMix(x,w,m,s). . .

% Funzione per il calcolo dei residuifunction res=MGres(par,x,y)

N=length(par)/3;w=par(1:N);m=par(N+1:2*N);s=par(2*N+1,3*N);

MG=GaussMix(x,w,m,s);

res=sum((MG-y).^2);

Es 4: mistura di Gaussiane 1

% Calcolo l’istogramma[yh,xh]=hist(I,100);

% Opzioni per l’ottimizzazione options=optimset('Display’,Iter');

% Numero di GaussianeN= … ;

% Inizializzazione dei parametrim0=[1:N]/N+1;s0=ones(1,N)*0.1;w0=ones(1,N)/N;beta0=[w0,m0,s0];

% Ottimizzazionebetaopt=lsqnonlin(@GMMres,beta0,[],[],options,xh,yh);

Es 4: Segmentazione

% Calcolare la pdf di tutti i voxel dell’immagine per ogni gaussianadati=I(:);

…y(mode,:)=w(mode)*Gauss(dati,m(mode),s(mode));…

% Calcolare per ogni voxel qual è la gaussiana col valore più elevato% Sugg: utilizzare opportunamente la funzione max()mgau=…

% Riportare il vettore alla dimensione dell’immagine:Iseg=reshape(mgau,size(I));

N

i

g

ii

i

ewgGM1

5.02

)(

N

j

g

j

j

g

ii

j

j

i

i

ew

e

gp

1

5.0

5.0

2

2

1

)(

M

jji

M

jtjji

ti

M

jji

M

jjji

ti

M

jji

ti

gp

ggp

gp

ggp

M

gp

w

1

1

21

1

1

11

11

)(

)ˆ)((ˆ

)(

)(

ˆ

)(

ˆ

Es 5: EM

Es 5: EM

1. Utilizzare la funzione per calcolare la pdf

2. Creare una funzione per calcolare medie, varianze e pesi campionari

3. Creare uno script che iteri il procedimento

Statistiche campionarie% Funzione per calcolare una gaussiana unidimensionale% m,s parametri scalarifunction y=Gauss(x,m,s). . .

% Funzione per calcolare la mistura di gaussiane% w,m,s parametri vettoriali (Nx1)function MG=GaussMix(x,w,m,s). . .

% Funzione per il calcolo dei residuifunction [wnew,mnew,snew]=MGupdate(dati,w,m,s,mode)

y=Gauss(dati,m(mode),s(mode));MG=GaussMix(dati,w,m,s);

p=y/MG;

wnew=…;mnew=…;snew=…;

EM% Numero di GaussianeN= … ;

% Inizializzazione dei parametrim0=[1:N]/N+1;s0=ones(1,N)*0.1;w0=ones(1,N)/N;

% Numero massimo di iterazioniMaxIter= … ;

% Ciclo principale EMiter=1;dati=I(:);

% fino a che il numero di iterazioni non supera MaxIter% aggiorno i valori dei parametri di ogni gaussiana……

[wnew,mnew,snew]=MGupdate(dati,w,m,s,mode);……