Imaging in Biologia e Medicina

Thresholding

Le operazioni di thresholding hanno lo scopo di selezionare regioni diverse di un'immagine sulla base dell'intensità luminosa. Questa operazione spesso si avvale di un fase preliminare di manipolazione del contrasto (contrast stretching ) per marcare il più possibile differenze di luminosità all'interno dell'immagine. L'output di una operazione di thresholding è un'immagine binaria formata quindi da soli valori 0 o 1 (alcuni linguaggi rappresentano i valori logici con le espressione dell'algebra booleana 'true' e 'false')

Nell'immagine di prova sono raffigurati dei grani di riso e l'obiettivo dell'esercizio è di sperimentare le qualità (e i limiti) di alcuni metodi volti a separare gli oggetti (i grani) dallo sfondo scuro.

La funzione graythresh implementa alcuni algoritmi di calcolo automatico della soglia di un'immagine. Il metodo più popolare (che si ottiene per default chiamando graythresh con la matrice dell'immagine come unico argomento) è il metodo di Otsu. Nel metodo di Otsu si assume che sfondo e primo-piano appartengano a due popolazioni di valori statisticamente separate. Il valore di soglia è tale da minimizzare la varianza intra-classe (e di conseguenza massimizzare quella inter-classe)

Questa parte dell'esercitazione funziona con Octave perché l'omonima funzione di Matlab implementa solo l'opzione con il metodo di Otsu.

graythresh(f,'metodo di thresholding');
oppure
graythresh(f,'metodo di thresholding', opzione);

Dove f è la matrice dell'immagine e opzioni una variabile opzionale che può essere richiesta da uno specifico metodo.

%
% -- thresholding.m
%
% confronto del risultato di thresholding con vari metodi
%
%

pkg load image

img = imread('grani-riso.png');

% L'immagine dovrebbe essere NxMx3, quindi in RGB

size(img)

% suddividiamo la finestra in 8 diagrammi (2x4)

subplot(2,4,1)
imshow(img);
title('Originale');

% Per default graythresh calcola la soglia con il metodo di Otsu.
% Il metodo di Otsu definisce come soglia l'intensità luminosa per
% la quale la varianza intra popolazione dei pixel di background e 
% foreground e minima (e di conseguenza la varianza inter-popolazione
% è massima)

subplot(2,4,2)
otsu = graythresh (img);
imshow(im2bw(img,otsu));
title('Otsu');

% Il metodo di concavità cerca di determinare una soglia sulla
% base di un'analisi della struttura delle concavità dell'istogramma
% cercando il miglior punto di minimo secondo criteri di ottimizzazione
% volti a minimizzare l'influenza del rumore o di sorgenti spurie

subplot(2,4,3);
concavity=graythresh (img,'concavity');
imshow(im2bw(img,concavity));
title('concavity');

% 'intermodes' assume che l'istogramma abbia una forma bimodale (2 picchi)
% e determina la soglia come la media dei massimi

subplot(2,4,4)
intermodes=graythresh (img,'intermodes');
imshow(im2bw(img,intermodes))
title('intermodes');

% Il metodo di entropia massima calcola la soglia
% come il valore che separa l'istogramma in due regioni tali
% che l'entropia globale sia massima

subplot(2,4,5)
max_entropy=graythresh(img,'MaxEntropy');
imshow(im2bw(img,max_entropy))
title('MaxEntropy');

% Minimizzazione di una funzione dell'errore
% nella determinazione dell'istogramma come distribuzione
% normale bimodale

subplot(2,4,6)
min_error=graythresh(img,'MinError');
imshow(im2bw(img,min_error));
title('MinError');

% Calcolo della soglia con il vincolo
% che i momenti di ordine 0 (integrale), 1 (media), 2 (varianza)
% siano preservati nell'immagine binaria che si ottiene
% dalla selezione sulla base della soglia

subplot(2,4,7)
moment=graythresh(img,'moments');
imshow(im2bw(img,moment));
title('moments');

% 'MaxLikelihood': calcolo della soglia secondo massimizzazione
% della funzione di verosimiglianza statistica

subplot(2,4,8)
max_likelihood=graythresh(img,'MaxLikelihood');
imshow(im2bw(img,max_likelihood));
title('MaxLikelihood');

% --- thresholding.m

Nonostante l'apparente semplicità dell'istogramma di quest'immagine non tutti i metodi appaiono soddisfacenti. In particolare ci sono alcuni grani nella regioni inferiore che non vengono efficacemente identificati senza sacrificare l'efficacia della soglia nella parte superiore dell'immagine. Eppure anche ad una sommaria valutazione ciascun grano di riso appare ben definito e separato dagl'altri. La ragione della difficoltà sta nel fatto che lo sfondo, seppure omogeneo, ha un gradiente di luminosità dal basso verso l'alto e quindi, anche se la nostra percezione vede chiaramente i grani in ciascuna regione dell'immagine, in realtà i grani di riso della parte inferiore hanno una luminosità molto vicina a quella dello sfondo nella zona superiore e quindi non possono essere chiaramente separati sulla base di un confronto con un parametro globale (la soglia)

Un altro metodo di calcolo della soglia si basa sulla determinazione dei percentili dell'istogramma che stimiamo siano dominio dello sfondo rispetto agli oggetti in primo piano. Questo è un metodo che richiede tentativi diversi (e per questo dipende dalla valutazione soggettiva dello sperimentatore), ma permette di poter esplorare l'effetto della soglia in diversi intervalli.

Questo esempio deve essere salvato in un file (e.g. percthresh.m) nella vostra cartella di lavoro, eventualmente modificato e quindi da Octave lanciato digitando il suo nome (senza la finale .m)

%
% -- percthresh.m
%
% confronto del risultato di thresholding con vari metodi
%
%

pkg load image
img = imread('grani-riso.png');

% L'immagine dovrebbe essere NxMx3, quindi in RGB

size(img)

% suddividiamo la finestra in 8 diagrammi (2x4)

imshow(img);
title('Originale');

% usiamo la forma graythresh(img,'percentile',perc)
% per determinare la soglia. Il valore di soglia calcolato
% corrisponde ad un certo valore di percentile della distribuzione
% cumulativa dell'istogramma (integrale dell'istogramma 
% da 0 a un certo valore normalizzato al numero totale di
% pixel presenti). Il valore della soglia 'perc' è 

% calcoliamo le soglie nell'intevallo 0.2-0.9 per 8 valori

figure

pc1 = 0.2
pc2 = 0.9

for p = [1:1:8]

    perc = pc1 + (pc2 - pc1) * p / 8;
    soglia = graythresh(img,'percentile',perc);
    subplot(2,4,p)
    caption = sprintf('Soglia: %f',perc);

    imshow(im2bw(img,soglia));
    title(caption)

end

% -- percthresh.m

Checker Shadow Illusion

Illusione di Adelson Un esempio interessante che spiega bene la nostra difficoltà percettiva a valutare oggettivamente certe caratteristiche di luminosità è la checker shadow illusion inventata da Edward H. Adelson. L'illusione si può creare attraverso una scena virtuale dove è mostrato un oggetto cilindrico che proietta un'ombra verso una scacchiera, una superficie costituita quindi da caselle chiare intervallate a caselle scure. Chiedetevi se la casella B appare alla vostra percezione più luminosa della casella A. In una immagine senza l'ombra del cilindro e in illuminazione omogenea la casella A apparirebbe più scura di B, ma tuttavia anche nell'ombra B continua ad sembrare più chiara. La seconda immagine rivela che questa è un illusione se si segue con l'occhio il ponte a luminosità costante aggiunto all'immagine che collega le caselle A e B. Dimostrazione (Se non riuscite a vedere con chiarezza questo fatto usate 2 fogli di carta per isolare il ponte dal resto dell'immagine) Potete trovare una rappresentazione interattiva dell'illusione di Adelson sul già citato sito di Michael Bach

Problemi di thresholding come quello dell'immagine dei grani di riso hanno bisogno di un preprocessing dell'immagine in modo da ridurre il gradiente di luminosità dello sfondo, riducendo così la variabilità locale e potendo quindi usare una soglia globale per la selezione, oppure dei metodi di selezione della soglia adattivi