Esercizi parte 6 - DEI - Università degli Studi di Padova

Università degli Studi di Padova - Facoltà di Ingegneria
Corso di Laurea in Ingegneria Biomedica
A.A. 2006-2007
Laboratorio di
Elaborazione di Dati, Segnali e
Immagini Biomediche
(Parte 6)
Prof. Giovanni Sparacino
Dipartimento di Ingegneria dell’Informazione
Università di Padova
e-mail: [email protected]
web: http: www.dei.unipd.it/~gianni
1
IMAGE PROCESSING TOOLBOX: ALCUNI COMANDI UTILI
imshow (im)
mostra l’immagine in im, assumendo che im contenga già valori quantizzati
(interi tra 0 e 255)
uint8(A)
trasforma A (in doppia precisione) in una matrice di interi a 8 bit (0-255), saturando su 0 e 255
double(im)
passa da interi a valori in doppia precisione (utilizzabili per elaborazioni quali somme, radici,…)
imhist(im)
mostra l’istogramma dei valori di grigio dell’immagine contenuta nella matrice im
filter2(B,im)
filtra la matrice in im secondo la maschera contenuta in B
2
IMAGE PROCESSING TOOLBOX: ALCUNI COMANDI UTILI
im=mat2gray(A,[Amin,Amax])
trasforma la matrice A in una immagine in B/N im graficabile. I valori di im saranno compresi tra
0 e 1. Valori di A minori di Amin saturano sul nero, valori in A maggiori di Amax saturano sul
bianco
colormap(gray)
fa sì che la prossima figura venga riportata in livelli di grigio
imagesc(im)
prende il contenuto della matrice im, la scala, e ne fa un’immagine grafica
panoramica sull’ IPT
Dalla command window di Matlab lanciare: demo. Selezionare "Toolboxes" e quindi "Image
Processing". Provare le demo delle tecniche studiate a lezione (Histogram Equalization,
Intensity Adjustment, Filtering, ...)
3
ESERCIZIO 1
Premessa: I comandi in blu richiedono l’Image Processing Toolbox (IPT), per il quale il DEI ha solo 50
licenze. Se l’IPT non fosse disponibile sulla propria stazione, lavorare a coppie (oppure completare fuori
orario o a casa)
Caricare con load imdemos le immagini in B/N e controllare con whos il contenuto.
Scegliere un’immagine (ad esempio saturno) e assegnarla ad una matrice di servizio
(es. originale=saturn). Se l’IPT è disponibile, usare imshow per mostrare nel riquadro 1
l’immagine originale.
Sfruttando il comando double, ottenere una versione
dell’immagine in doppia precisione e quindi "lavorabile" algebricamente.
a) Fare il negativo dell’immagine scelta. Sfruttando l’istruzione uint8, creare una matrice
a 256 livelli di grigio graficabile. Se l’IPT è disponibile, mostrarla nel riquadro 2 usando il
comando imshow
b) Fare una versione schiarita. Sfruttando l’istruzione uint8, creare una matrice a 256
livelli di grigio graficabile. Se l’IPT è disponibile, mostrarla nel riquadro 3
c) Fare una versione scurita. Sfruttando l’istruzione uint8, creare una matrice a 256
4 livelli
di grigio graficabile. Se l’IPT è disponibile, mostrarla nel riquadro 4
d) Applicare una procedura di slicing, ad esempio portando a luminosità 200 tutti i pixel
aventi luminosità compresa tra 80 e 100. Sfruttando l’istruzione uint8, creare una
matrice a 256 livelli di grigio graficabile. Se l’IPT è disponibile, mostrare il risultato nel
riquadro 5
e) Se l’IPT è disponibile, usare il comando histeq per equalizzare l’immagine, e mostrare
il risultato nel riquadro 6
f) Costruire una maschera a supporto 3X3 per fare un filtro passa-basso e applicarla
mediante il comando filter2. Sfruttando l’istruzione uint8, creare una matrice a 256 livelli
di grigio graficabile. Se l’IPT è disponibile, mostrare il risultato nel riquadro 7
g) Costruire una maschera a supporto 6X6 per fare un filtro passa-basso e applicarla
mediante il comando filter2. Sfruttando l’istruzione uint8, creare una matrice a 256 livelli
di grigio graficabile. Se l’IPT è disponibile, mostrare il risultato nel riquadro 8
h) Costruire una maschera a supporto 3X3 per fare un filtro passa-alto e applicarla
mediante il comando filter2. Sfruttando l’istruzione uint8, creare una matrice a 256 livelli
di grigio graficabile. Se l’IPT è disponibile, mostrare il risultato nel riquadro 9
5
% esercitazione del 14 marzo 2007
% operatori puntuali e filtri 2D
close all
load imdemos
imoriginale=saturn; % scelgo l'immagine su cui lavorare (vd whos)
% slicing
soglia1=80;
soglia2=120;
livello=200;
for k=1:nr
for j=1:nc
[nr,nc]=size(imoriginale);
if im(k,j)>soglia1 & im(k,j)<soglia2
subplot(331)
im2(k,j)=livello;
imshow(imoriginale)
else im2(k,j)=im(k,j);
title('originale')
end
end
% procuro una versione "lavorabile" in doppia precisione dell'immagine
end
im=double(imoriginale);
imquantiz=uint8(im2);
subplot(335)
% negativo
imshow(imquantiz)
im2=255-im;
title('slicing')
imquantiz=uint8(im2);
subplot(332)
imshow(imquantiz)
title('negativo')
% schiarito
im2=sqrt(255*im);
imquantiz=uint8(im2);
subplot(333)
imshow(imquantiz)
title('schiarito')
% scurito
im2=(im.^2)/255;
imquantiz=uint8(im2);
subplot(334)
imshow(imquantiz)
title('scurito')
% filtraggio passa-basso selettivo
F=1/64*ones(8,8);
im2=filter2(F, im);
imquantiz=uint8(im2);
subplot(338)
imshow(imquantiz)
title('filtrato passa-basso selettivo')
% filtro sharpening
F=-ones(3,3);
F(2,2)=8;
im2=filter2(F, double(im));
imquantiz=uint8(im2);
subplot(339)
imshow(imquantiz)
title('filtrato passa-alto')
% equalizzazione
im2=histeq(uint8(im),64);
imquantiz=uint8(im2);
subplot(336)
imshow(imquantiz)
title('equalizzato')
% filtraggio passa-basso
F=1/9*ones(3,3);
im2=filter2(F, im);
imquantiz=uint8(im2);
subplot(337)
imshow(imquantiz)
title('filtrato passa-basso')
6
originale
negativo
s chiarito
s curito
s licing
equalizzato
filtrato pas s a-bas s o
filtrato pas s a-bas s o s elettivo
filtrato pas s a-alto
7
originale
negativo
s chiarito
s curito
s licing
equalizzato
filtrato pas s a-bas s o
filtrato pas s a-bas s o s elettivo
filtrato pas s a-alto
8
ESERCIZIO 2
Per l’esercizio 1, se IPT è disponible ottenere per ogni elaborazione l’istogramma dei
livelli di grigio usando il comando imhist (mettere gli istogrammi in figura 2 sui subplot
corrispondenti a quelli delle immagini in figura 1) . Osservare come i diversi operatori
puntuali si riflettano sull’istogramma
9
Soluzione
….
figure(1)
subplot(331)
imshow(imoriginale)
title('originale')
figure(2)
subplot(331)
imhist(imoriginale)
title('hist originale')
figure(1)
….
(saturno)
his t originale
his t negativo
his t s chiarito
1000
1000
1000
500
500
500
0
0
0
100
200
his t s curito
0
0
100
200
his t s licing
1000
1000
1000
500
500
500
0
0
0
100
200
his t pas s a-bas s o
0
100
200
his t equalizzato
0
100
200
his t pas s a-alto
0
0
100
200
his t pas s a-bas s o s el
1000
1500
1000
1000
500
500
500
0
0
0
100
200
0
0
100
200
0
100
200
10
ESERCIZIO 3
Caricare nel workspace il contenuto del file imdemos.mat e scegliere una
matrice immagine su cui lavorare.
Usare il comando double per trasformarla in una matrice su cui è possibile
fare operazioni algebriche.
a) Generare una matrice con elementi estratti da rumore bianco (usare media
nulla e SD “piccola”, es. 10) e costruire, mediante addizione, una versione
“rumorosa” dell’immagine scelta. Sfruttando l’istruzione uint8, creare una matrice a
256 livelli di grigio graficabile. Se l’IPT è disponibile, mostrare il risultato nel riquadro 2
b) Sfruttando l’istruzione randn seguita dalla round, generare una matrice che
simula rumore sale e pepe e costruire, mediante addizione, una versione
“rumorosa” dell’immagine scelta. Sfruttando l’istruzione uint8, creare una matrice a
256 livelli di grigio graficabile. Se l’IPT è disponibile, mostrare il risultato nel riquadro 3
11
c) Fare un filtraggio passa-basso (supporto 3X3) dell’immagine sporcata con
rumore gaussiano Se l’IPT è disponibile, mostrare il risultato nel riquadro 5
d) Fare un filtraggio passa-basso (supporto 3X3) dell’immagine sporcata con
rumore sale e pepe Se l’IPT è disponibile, mostrare il risultato nel riquadro 6
e) Usando il comando medfilt2, fare un filtraggio mediano (supporto 3X3)
dell’immagine sporcata con rumore gaussiano. Se l’IPT è disponibile, mostrare il
risultato nel riquadro 5 di figura 2 (nella prima riga di 3 riquadri mostrare l’immagine
originale e le due versioni rumorose)
f) Usando il comando medfilt2, fare un filtraggio mediano (supporto 3X3)
dell’immagine sporcata con rumore sale e pepe. Se l’IPT è disponibile, mostrare
il risultato nel riquadro 6 di figura 2
12
Passa-basso
originale
con rumore gaus s iano
con rumore s ale & pepe
pas s a-bas s o
pas s a-bas s o
13
Mediano
originale
con rumore gaus s iano
con rumore s ale & pepe
mediano
mediano
14
% esercitazione del 14 marzo 2007
% filtraggio 2D: passabasso vs mediano
clear all
close all
% carico immagini e scelgo quella di lavoro
load imdemos
originale=saturn;
% plotto originale
[nr,nc]=size(originale);
subplot(231)
imshow(originale)
% filtraggio passa-basso su immagine sporcata gaussiano
F=1/9*ones(3,3);
imtmp=filter2(F, imrumwn);
impb=uint8(imtmp);
subplot(235)
imshow(impb)
title('passa-basso')
% filtraggio passa-basso su immagine sporcata sale & pepe
F=1/9*ones(3,3);
im5=filter2(F, imrumsp);
im5=uint8(im5);
subplot(236)
imshow(im5)
title('passa-basso')
figure(2)
subplot(231)
imshow(originale)
title('originale')
subplot(232)
imshow(imqtzwn)
title('con rumore gaussiano')
subplot(233)
imshow(imqtzsp)
title('con rumore sale & pepe')
% procuro una versione "lavorabile"
im=double(originale);
title('originale')
% aggiungo rumore bianco
v=10*randn(nr,nc);
imrumwn=im+v;
imqtzwn=uint8(imrumwn);
subplot(232)
imshow(imqtzwn)
title('con rumore gaussiano')
% filtraggio mediano su immagine sporcata gaussiano
im5=medfilt2(double(imqtzwn),[3 3]);
im5=uint8(im5);
subplot(235)
imshow(im5)
title('mediano')
% aggiungo rumore sale e pepe
v=0.25*randn(nr,nc);
v=255*round(v);
imrumsp=im+v;
imqtzsp=uint8(imrumsp);
subplot(233)
imshow(imqtzsp)
title('con rumore sale & pepe')
% filtraggio mediano su immagine sporcata sale & pepe
im5=medfilt2(double(imqtzsp),[3 3]);
im5=uint8(im5);
subplot(236)
imshow(im5)
title('mediano')
15
Esercizio 4
a) Scrivere un codice Matlab che prepari le matrici 256X256 corrispondenti ai quattro
grating sinusoidali sotto. Se l’IPT è disponibile, mostrarle in quattro riquadri
grating verticale
grating orizzontale
grating s fas ato
grating s fas ato, frequenza lungo X maggiore
16
Soluzione a)
clear al close all
nc=256;
nr=256;
% grating con sfasamento
for j=1:nr
for k=1:nc
A(j,k)=floor(128+128*sin(2*pi*(1/X0*k+1/Y0*j)));
end
end
im=uint8(A);
subplot(223)
imshow(im)
title('grating sfasato')
% -grating sinusoidale a linee verticali
A=zeros(nr,nc);
X0=100;
x=(1:1:nr);
riga=floor(128+128*sin(2*pi/X0*x));
for k=1:nr
A(k,:)=riga;
end
im=uint8(A);
subplot(221)
imshow(im)
title('grating verticale')
% grating con sfasamento - frequenza maggiore lungo X
for j=1:nr
for k=1:nc
A(j,k)=floor(128+128*sin(2*pi*(4*1/X0*k+1/Y0*j)));
end
end
im=uint8(A);
subplot(224)
imshow(im)
title('grating sfasato, frequenza lungo X maggiore')
% -grating sinusoidale a linee orizzontali
A=zeros(nr,nc);
Y0=100;
y=(1:1:nc)';
colonna=floor(128+128*sin(2*pi/Y0*y));
for k=1:nc
A(:,k)=colonna;
end
im=uint8(A);
subplot(222)
imshow(im)
title('grating orizzontale')
17
b) Se IPT è disponible, ottenere per ogni grating l’istogramma dei livelli di grigio usando
il comando imhist (mettere gli istogrammi in figura 2 sui subplot corrispondenti a quelli
delle immagini in figura 1) . Come dovrebbe venire l’istogramma ? Perché nei casi
orizzontale e verticale non è conforme alle attese ?
his t grating verticale
his t grating orizzontale
1000
1000
500
500
0
0
0
50
100
150
200
250
his t grating s fas ato
0
1000
500
500
0
0
50
100
150
100
150
200
250
his t grating s fas ato, frequenza lungo X maggiore
1000
0
50
200
250
18
0
50
100
150
200
250