1 / 15

Laboratorio di Neuroingegneria Segmentazione A.A . 2012-2013

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. Caricamento. Sono a disposizione due file NIFTI : T1.nii e FLAIR.nii

Download Presentation

Laboratorio di Neuroingegneria Segmentazione A.A . 2012-2013

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  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); … …

More Related