Laboratorio di Neuroingegneria Segmentazione A.A . 2012-2013

1. Laboratorio di Neuroingegneria Segmentazione A.A. 2012-2013

2. Cosa faremo oggi • Caricamento di volumi MRI cerebrale • T1 • FLAIR • Calcolo della soglia di Otsu • Stima della mistura di Gaussiane

3. Caricamento Sono 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 NIFTI xT1=T1.img; % Estraggo una fetta campione xim=xT1(:,:,70); % Normalizzo di dati I=xim/max(xim(:));

4. Es1: Sogliatura Scegliere un valore casuale T0compreso 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

5. Es 1: Sogliatura Provare l’algoritmo seguente: % 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;

6. 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

7. 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)

8. 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

9. Es 4: mistura di Gaussiane % Funzione per calcolare una gaussiana unidimensionale % m,s parametri scalari functiony=Gauss(x,m,s) . . . % Funzione per calcolare la mistura di gaussiane % w,m,s parametri vettoriali (Nx1) functionMG=GaussMix(x,w,m,s) . . . % Funzione per il calcolo dei residui functionres=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);

10. Es 4: mistura di Gaussiane 1 % Calcolo l’istogramma [yh,xh]=hist(I,100); % Opzioni per l’ottimizzazione options=optimset('Display’,Iter'); % Numero di Gaussiane N= … ; % Inizializzazione dei parametri m0=[1:N]/N+1; s0=ones(1,N)*0.1; w0=ones(1,N)/N; beta0=[w0,m0,s0]; % Ottimizzazione betaopt=lsqnonlin(@GMMres,beta0,[],[],options,xh,yh);

11. Es 4: Segmentazione % Calcolare la pdf di tutti i voxel dell’immagine per ogni gaussiana dati=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));

12. Es 5: EM

13. Es 5: EM • Utilizzare la funzione per calcolare la pdf • Creare una funzione per calcolare medie, varianze e pesi campionari • Creare uno script che iteri il procedimento

14. Statistiche campionarie % Funzione per calcolare una gaussiana unidimensionale % m,s parametri scalari functiony=Gauss(x,m,s) . . . % Funzione per calcolare la mistura di gaussiane % w,m,s parametri vettoriali (Nx1) functionMG=GaussMix(x,w,m,s) . . . % Funzione per il calcolo dei residui function[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=…;

15. EM % Numero di Gaussiane N= … ; % Inizializzazione dei parametri m0=[1:N]/N+1; s0=ones(1,N)*0.1; w0=ones(1,N)/N; % Numero massimo di iterazioni MaxIter= … ; % Ciclo principale EM iter=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); … …