Imaging in Biologia e Medicina

Istogrammi e "Point Operations"

L'istogramma di un'immagine rappresenta la distribuzione della luminosità dei pixel e riesce a dare una rappresentazione qualitativa e immediata dei problemi di esposizione di un'immagine. cameraman.tiff Il comando di Matlab che calcola l'istogramma è imhist e nella sua forma più semplice apre una finestra con l'istogramma e una barra di legenda della relazione intensità / frequenza. Il comando ammette anche 2 argomenti di output dove vengono memorizzati il numero di conteggi per ogni valore di pixel e l'array dei valori dei pixel nell'immagine. Per immagini con una bit-depth > 8 può essere utile usare la forma del comando con 2 argomenti, essendo il secondo argomento il numero di valori (bins ) da rappresentare nell'istogramma. istogramma dell'esempio L'istogramma non rappresenta un'immagine in modo univoco: è sempre possibile calcolare un istogramma della distribuzione di pixel, ma data una distribuzione non è possibile ricalcolare l'immagine ad essa associata, tuttavia l'istogramma permette di stimare due caratteristiche essenziali di un immagine il contrasto e la dinamica . Il contrasto può essere definito in vari modi, ma in ogni caso è essenzialmente proporzionale alla differenza tra la luminosità più alta e quella più bassa presenti nell'immagine. La dinamica è invece il numero dei livelli di intensità diversi presenti nell'immagine.

E' semplice determinare con Matlab la dinamica di un'immagine. Scaricate l'immagine xray-chest.png . L'esempio fa uso delle funzione find che cerca in un vettore o una matrice gli elementi che soddisfano una determinata condizione e la funzione numel che ritorna il numero di elementi (o oggetti) contenuti in un array o matrice.

Differenze nelle librerie grafiche e,in secondo ordine, nell'implementazione del codice portano a risultati differenti tra Matlab e Octave o tra diverse installazioni della stessa applicazione

%
% Determinazione della dinamica di un'immagine

% carichiamo quest'immagine derivata da un'altra 
% immagine generata da un sistema
% a raggy X per diagnostica medica

% prima di eseguire i comandi di questo esercizio
% leggete il significato della funzione 'find'
% digitando il comando
%
% help find
%

% questa linea deve essere eseguita solo con Octave
% e solo all'inizio della sessione di lavoro
    
% in uno script che deve preservare la compatibilità
% con matlab si potrebbe mettere questo test
%
%if (exist('OCTAVE_VERSION', 'builtin') ~= 0)
%    pkg load image
%end
%

imgrgb = imread('xray-chest.png');

class(imgrgb)
size(imgrgb)

% convertiamo l'immagine RGB in una immagine grayscale

img = rgb2gray(imgrgb);

imshow(img);
size(img)

% calcoliamo l'istogramma con la funzione imhist memorizzando 
% l'array dei conteggi dell'istogramma e l'array dei valori a
% cui i conteggi si riferiscono

[counts,values] = imhist(img);

% gli indici al vettore 'counts' con conteggi > 0 (cioè valori
% di intensità usate nell'immagine) si possono ottenere con la
% funzione find.

find(counts > 0)
ans =

     1
     9
    20
    29
    31
    39
    48
    58
    59
    66
    76
    84
    86
    92
   102

% il numero di elementi che caratterizzano la dinamica è quindi

numel(find(counts > 0))
ans =  15

% e i valori dei conteggi presenti.

counts(find(counts > 0))
ans =

   40998
    1887
    2079
    2775
    3282
    4302
    5961
    7662
    9540
   15258
   19782
   25308
   36417
   18216
    1413

La dinamica può essere cambiata per esempio rinunciando a parte della risoluzione dell'immagine con la funzione imresize . Verificate la dinamica dell'immagine dopo che è stata ridotta di dimensione per un quarto con metodo di interpolazione bilineare

img2=imresize(img,0.75,'bilinear')

Point Operations

Le trasformazioni dette Point Operations eseguono manipolazioni dei valori di intensità dei pixel senza cambiare la geometria o la dimensione dell'immagine e il sono lo strumento basilare per modificare il contrasto. L'obiettivo della modifica del contrasto potrebbe essere semplicemente il bisogno di migliorare la leggibilità di un'immagine, oppure quello di separare più nettamente tra loro regioni diverse (background (sfondo) e foreground (primo piano)) in modo da costruire criteri per la loro selezione. Queste sono generalmente dette operazioni di thresholding (da threshold : soglia). Un'altra ragione di trasformazione dell'intensità è la cosiddetta correzione di gamma che cerca di compensare le non linearità del dispositivo di output. La loro forma matematica generale è

$$ a' \leftarrow g(a) $$ $$ a' \leftarrow g(I(u,v),u,v) $$

Dato un valore di intensità di un pixel queste trasformazioni lo ricalcolano univocamente in un altro valore, di fatto alterando l'istogramma originale e ridistribuendo l'istogramma in un modo determinato. Osserviamo l'effetto di alcune trasformazioni su un'immagine di prova. Carichiamo l'immagine cameraman.tiff e mostriamo la sua distribuzione dei valori di intensità.

Correzione logaritmica

L'espressione della correzione logaritmica è

$$ img_{out}(u,v) = C_{\alpha} * log(1 + \alpha * img_{in}(u,v)) $$

Dove la costante C va determinata in funzione di alpha in modo tale che la funzione costituisca una mappa del range di valori dell'immagine su sè stessa. Il ruolo del parametro alpha è quello di determinare la curvatura e di conseguenza la non linearità della trasformazione, ampliando il contrasto all'interno di un intervallo di luminosità a danno dei valori rimanenti

%
% Carichiamo l'immagine originale e verifichiamo
% la sua rappresentazione interna

% con octave potrebbe essere necessario caricare il package 'image'
%

if (exist('OCTAVE_VERSION', 'builtin') ~= 0)

    pkg load image

end

img = imread('cameraman.tiff');

class(img)
size(img)

% Una volta verificato che l'immagine sia già 
% rappresentata in scala di grigi (come un solo piano, 
% cioè una semplice matrice NxM) normalizziamola 
% nell'intervallo [0,1]

img = mat2gray(img);

% Valore massimo e valore minimo dell'immagine

max(img(:))
min(img(:))

% costruiamo la trasformazione logaritmica, che
% chiameremo logtr, usando una 'anonymous function' 
% di MATLAB/Octave. La linea di comando seguente definisce 
% la funzione logtr(v,alfa) dove alfa è uno scalare che
% determina la curvatura della funzione, mentre v può essere 
% sia uno scalare (un semplice numero) che una
% matrice (come quella che rappresenta un'immagine)

logtr=@(v,alfa) log(1 + alfa * v) / log(1 + alfa)

% Suddividiamo la finestra della grafica in 6 sotto-pannelli
% e fissiamo il primo pannello come pannello corrente

subplot(2,3,1)

% stampiamo l'immagine nella sua forma originale

imshow(img)

% Infine mostriamo le immagini trasformate per
% ognuno dei valori di alfa. Usiamo il comando 'for' per
% eseguire gli stessi comandi per ognuno dei valori di alfa.
% La variabile 'p' controlla quale pannello della finestra 
% grafica riceverà l'output grafico. Questo indice viene
% incrementato ad ogni ciclo con la linea p = p + 1 

% questa sequenza di comandi può essere scritta all'interno 
% di un file con estensione .m per essere eseguito 
% (e.g. histlog.m)
%
% L'uso dei file-m per salvare programmi (script)
% e il significato del ciclo 'for ... end' sono spiegati
% nella lezione 
%
% http://imaging.biol.unipr.it/download/L2021-7
%

p = 2;
for alfa=[1 50 100 200 500] 
    subplot(2,3,p) 
    imshow(logtr(img,alfa)); 
    p = p + 1; 
end

% -- histlog.m

Il comando for di MATLAB/Octave definisce un blocco di codice che viene ripetuto n volte, una volta per ognuno dei valori contenuti nel vettore definito all'inizio

Correzione esponenziale

La trasformazione esponenziale ha la forma

$$ img_{out} = C_{\alpha} * ((1+\alpha)^{img_{in}(u,v)} - 1) $$

Dove il carattere '^' rappresenta l'elevazione a potenza. La funzione anonima che costituisce il nucleo della trasformazione è

exptr = @(v,alpha) ((1 + alpha) .^ v - 1) / alpha

Notate l'operatore '.^' che ritorna una una matrice avente come elementi il risultato dell'elevazione a potenza calcolato per ogni elemento Si può verificare che exptr(v,alpha) ha la stessa proprietà di mappare l'intervallo [0,1] su se stesso, così che ogni immagine in forma normalizzata viene convertita in un'altra immagine avente la stessa proprietà.

Trasformazione sigmoidale

Un'altra trasformazione è quella della formula

$$ 1 \over {1 + (m / f)^E} $$

g = @(f,E,m) 1 ./ (1 + (m ./ (double(f) + eps)) .^ E)

che dipende da 3 argomenti:

• m: un valore di soglia (livello di transizione)
• E: un esponente che determina la rapidità della transizione della funzione
• f: un valore di intensità del pixel (può essere la matrice completa)

Con E abbastanza grande l'output è sostanzialmente un'immagine binaria e cambiando m nell'intervallo [0,1] si ottiene una selezione di regioni aventi luminosità maggiore del valore di m

eps è una variabile scalare o matrice built-in che contiene il valore più piccolo rappresentabile dal processore. La funzione di questa variabile è di evitare che il denominatore della frazione possa diventare 0.

Complemento di un immagine (negativo)

Nonostante la complementazione dei valori di intensità di un'immagine sia concettualmente banale si deve tenere conto della rappresentazione binaria usata. Per questo esiste in Octave la funzione imcomplement che esegue i controlli necessari ed inverte l'immagine nel modo opportuno

imshow(imcomplement(img))

Esercizio

Costruzione di una funzione per le point operation logaritmica , esponenziale sigmoidale