Imaging in Biologia e Medicina

Calcolo della Soglia con Metodo Iterativo

L'algoritmo calcola una soglia di intensità luminosa per la segmentazione di un'immagine attraverso iterazioni successive. La soglia viene ricalcolata ad ogni iterazione fino a quando non diventa abbastanza stabile. In questo metodo viene fatta una stima iniziale a priori per una soglia e in base ad essa si fa una determinazione preliminare delle due due famiglie dei pixel di foreground e background , determinazione che cambierà ad ogni iterazione. Schematicamente:

1. Si fissa un valore iniziale di soglia T che si assume essere compreso tra le due distribuzioni di valori. Il valore della media globale può essere una scelta possibile
2. All'inizio di ogni iterazione viene salvato il valore calcolato della soglia per poterlo confrontare alla fine dell'iterazione. Se siamo alla prima iterazione la soglia da salvare è quella determinata a priori
3. Si determinano le due classi di pixel (foreground/background) e per ciascuna si calcola i valori medi m₁ e m₂
4. La nuova soglia viene ricalcolata come $$ T = \frac{m_{1} + m_{2}}{2} $$
5. La soglia appena calcolata viene confrontata con quella salvata all'inizio dell'iterazione. Se la differenza non è inferiore ad un valore dato allora si ritorna al secondo punto
6. Se la differenza è invece inferiore il ciclo termina e con l'ultima soglia calcolata si determina l'immagine binaria

La difficoltà sta tutta nella manipolazione effiente di queste due classi di elementi della matrice dell'immagine, in particolare il calcolo della media dei valori di ciascuna classe .

La funzione mean di Matlab/Octave calcola una media dei valori di una matrice matrice, in prima istanza producendo un array riga con le medie di ogni colonna. Usando la notazione (:) produce uno singolo scalare con la media di tutta la matrice. Quindi questa funzione così com'è quindi non può funzionare. Tuttavia Una soluzione è quella di partire dalla definizione di media

$$ m = \frac {\sum_{i=1}^{N} p_{i}}{N} $$

E combinarla con l'uso di matrici binarie per 'selezionare' di volta in volta solo i pixel di una classe rendendo ininfluenti ai fini del calcolo della media i pixel dell'altra classe.

Assumiamo di avere stabilito un stima della soglia: le due classi di pixel sono individuate dalle matrici binarie fore e back

img=imread('grani-riso.png');
riso = mat2gray(rgb2gray(img));
soglia = mean(riso(:));
fore = riso > soglia;
back = riso <= soglia;

Se moltiplichiamo la matrice riso con una della 2 matrici abbiamo un immagine dove mantengono la propria originaria intensità luminosa solo i pixel di una classe mentre gli altri sono posti a intensità zero (quindi neri). Per esempio

c1 = riso .* fore;
c2 = riso .* back;
subplot(1,2,1); imshow(c1)
subplot(1,2,2); imshow(c2)

Le matrici fore e back sono matrici binarie: i valori ammessi per i loro elementi sono solo 1 e 0. In particolare con 1 indichiamo se un determinato pixel appartiene ad una determinata classe. Di conseguenza la somma di tutti gli elementi di ciascuna è la numerosità della classe che rappresentano, mentre se sommiamo tutti gli elementi di c1 e c2 otteniamo un risultato a cui hanno contribuito solo i pixel di una determinata classe. Calcolare la media di ciascuna classe quindi è semplice e può essere condensato in una semplice funzione che può essere scritta in un m-file

function mediasel = mediasel(img,mask)
%
% Calcola la media solo dei valori in img
% che vengono selezionati da una matrice binaria
% 
%   Argomenti:
%
%       img:  immagine grayscale NxM
%       mask: maschera booleana (binaria) che
%             seleziona gli elementi che partecipano
%             al calcolo della media
%

% la funzione inizia controllando che le dimensioni
% dell'immagine e della maschera siano identiche. Se
% questa condizione non è verificata la funzione esce
% con un messaggio di errore.

% La funzione all() verifica che la matrice passata come
% argomento sia fatta solo da valori 1. In questo caso 
% significa che il risultato del confronto tra size(img) 
% e size(mask) deve essere vero per entrambe le dimensioni
% delle matrici

    if (~all(size(img) == size(mask))) 
        error('La matrice immagine e la maschera devono avere identiche dimensioni')
    end

% inoltre la matrice 'mask' deve essere booleana

    if (~islogical(mask))
        error('Il secondo argomento deve essere una matrice booleana');
    end

% l'immagine deve essere in formato double normalizzata
% nell'intervallo [0,1] (solitamente generata con 'mat2gray')

    if (~isfloat(img))
        error('L immagine deve essere normalizzata e in formato double')
    end
    
% la matrice selezione preserva i pixel selezionati
% dalla maschera mentre i rimanenti sono posti a zero. 

    selezione = img .* mask;

% si calcola l'accumulatore della media. I pixel
% che non fanno parte della classe non contribuiscono,
% ma non contano nemmeno nella determinazione del
% denominatore della frazione della media

% per calcolare la somma usiamo la notazione 
% che permette di individuare un insieme sparso
% di elementi della matrice tramite una matrice
% booleana

    somma = sum(selezione(:));

% La somma degli elementi della maschera è il numero
% dei pixel che appartengono alla classe individuata
% dalla maschera stessa

    N = sum(mask(:));
    if (N > 0)
        mediasel = somma / N;
    else
        mediasel = 0;
    end
end

% -- mediasel.m

E' ora possible eseguire tutta l'iterazione, ricalcolare la nuova soglia e confrontarla con quella iniziale per determinare se sono sufficientemente vicine.

function soglia_calcolata = iterthresh (img,max_iterazioni,target)
%
% modello di script per la determinazione della soglia
% con metodo iterativo.
%
%  soglia_calcolata = iterthresh(img)
%
% viene usato un algoritmo iterativo per la determinazione della
% soglia. L'algoritmo usa 0.1 come differenza tra soglie successive
% che termina l'esecuzione. In ogni caso vengono eseguite al massimo
% 100 iterazioni
%
%  soglia_calcolata = iterthresh(img,max_iterazioni)
% 
% Stesse funzionalità ma il secondo argomento specifica un
% nuovo numero massimo di iterazioni
%
%  soglia_calcolata = iterthresh(img,max_iterazioni,target)
%
% Questo script è solo un modello: per essere funzionante
% deve essere completato
%

% passo 1 dell'algoritmo

% nella variabile soglia memorizziamo la media globale
% dell'intensità luminosa di img

    if (nargin == 1)
        max_iterazioni = 100;
        target = 0.1;
    elseif (nargin == 2)
        target = 0.1;
    end

   % valore iniziale della soglia
    
    soglia = mean(img(:));
    ultima_soglia = 0;
    
   % contatore del numero di iterazioni
   
    niter = 0;

    % eseguiamo il ciclo while fino a quando 
    % la differenza tra la ultime 2 soglie calcolate
    % non è minore di 'target' e fino a quando il
    % valore di niter > max_iterazioni.
    % L'espressione della condizione è resa più complicata
    % dal fatto che comunque vogliamo eseguire almeno
    % una iterazione

    while (((abs(soglia-ultima_soglia) > target) && ...
            (niter < max_iterazioni)) || (niter == 0))

        % salviamo il valore corrente della soglia
        % così che alla fine della iterazione
        % possiamo confrontarlo con il nuovo valore
        % della soglia ricalcolata
        
        ultima_soglia = soglia;
        
        % determiniamo le 'maschere' booleane
        % sulla base della soglia
        
        fore = img > soglia;
        back = img <= soglia;
               
        m1 = mediasel(img,fore);
        m2 = mediasel(img,back);

        % calcoliamo la nuova soglia. 
        % Togliendo il ; alla fine di questa linea
        % potete vedere la sequenza delle soglie 
        % calcolate
        
        soglia = (m1 + m2) / 2;
        
        % incrementiamo il numero delle iterazioni

        niter = niter + 1;

    end

    soglia_calcolata = soglia; 
end

% -- iterthresh.m