数字图像处理上机

数字图像处理上机

班级：

学号：

姓名：

1. 产生右图所示亮块图像 f 1(x,y)（128×128大小，暗处=0，亮处=255），对其进行

FFT ：

（1） 同屏显示原图f 1和FFT(f 1)的幅度谱图； （2） 若令f 2(x,y)=（-1）

x+y

f 1(x,y)，重复以上过程，比较二者幅度谱的异同，简

述理由；

（3） 若将f 2(x,y)顺时针旋转45度得到f 3(x,y)，试显示FFT(f 3)的幅度谱，并与

FFT(f 2)的幅度谱进行比较。 解答： （1）图像：

（2）图像：

相同点：频谱的实质没有改变，幅度等都没有发生变化。

不同点：f2的频谱是对f1频谱的移位，它时f1的频谱从原点（0，0）移到了中心（64，64），而得到了一个完整的频谱。

（3）图像：

程序：f1=zeros(128,128);

for i=38:1:90

for j=58:1:70

f1(i,j)=255;

end

end

figure(1)

subplot(1,2,1);

imshow (f1);

subplot(1,2,2);

imshow (fft2(f1));

% f2(x,y)=（-1）^(x+y)* f1(x,y)

for i=1:1:128

for j=1:1:128

f2(i,j)=(-1)^(i+j)*f1(i,j);

end

end

figure(2);

subplot(1,3,1);

imshow (f2);

f3=imrotate(f2,-45,'bilinear');%将f2顺时针旋转45度

subplot(1,3,2);

imshow(fft2(f2));%显示f2的频谱

subplot(1,3,3);

imshow(fft2(f3));%显示f3的频谱

2. 对256 256大小、256级灰度的数字图像lena进行频域的理想低通、高通滤波，同屏显示原图、幅度谱图和低通、高通滤波的结果图。

程序：

figure(1);

fid=fopen('d:\img\lena.img','r');

data=(fread(fid,[256,256],'uint8'))';

subplot(2,2,1)

imagesc(data);

colormap(gray);

title('LENA','Color','r');

subplot(2,2,2);

imshow(fft2(data));

s=fftshift(fft2(data));

[M,N]=size(s); %分别返回s的行数到M中，列数到N中

n=2; %对n赋初值

%GLPF滤波，d0=5，15，30（程序中以d0=30为例）

d0=30; %初始化d0

n1=floor(M/2); %对M/2进行取整

n2=floor(N/2); %对N/2进行取整

for i=1:M

for j=1:N

d=sqrt((i-n1)^2+(j-n2)^2); %点（i,j）到傅立叶变换中心的距离

h=1*exp(-1/2*(d^2/d0^2)); %GLPF滤波函数

s(i,j)=h*s(i,j); %GLPF滤波后的频域表示

end

end

s=ifftshift(s); %对s进行反FFT移动

%对s进行二维反离散的Fourier变换后，取复数的实部转化为无符号8位整数

s=uint8(real(ifft2(s)));

subplot(2,2,3); %创建图形图像对象

imshow(s);

p=fftshift(fft2(data));

[M,N]=size(p); %分别返回p的行数到M中，列数到N中

n=2; %对n赋初值

%GLPF滤波d1=30

d1=30; %初始化d1

n3=floor(M/2); %对M/2进行取整

n4=floor(N/2); %对N/2进行取整

for i=1:M

for j=1:N

dd=sqrt((i-n3)^2+(j-n4)^2); %点（i,j）到傅立叶变换中心的距离

h1=1-exp(-1/2*(dd^2/d1^2)); %GHPF滤波函数

p(i,j)=h1*p(i,j); %GHPF滤波后的频域表示

end

end

p=ifftshift(p); %对p进行反FFT移动

%对s进行二维反离散的Fourier变换后，取复数的实部转化为无符号8位整数

p=uint8(real(ifft2(p)));

subplot(2,2,4); %创建图形图像对象

imshow(p);

3.对给定的两种128 128、256级灰度的数字图像（图像磁盘文件名分别为Fing_128.img（指纹图）和Cell_128.img（显微医学图像）进行如下处理：

（1）对原图像进行直方图均衡化处理，同屏显示处理前后图像及其直方图，比较异同，并回答为什么数字图像均衡化后其直方图并非完全均匀分布。

（2）对原图像加入高斯噪声，用4-邻域平均法平滑加噪声图像（图像四周边界不处理，下同），同屏显示原图像、加噪图像和处理后的图像。

①不加门限；

② 加门限 2(,)T f m n =，（其中∑∑

=i

j

j i f N

n m f ),(1),(2

）

解答： （1）

异同：由于原图像中目标物的灰度主要集中于低亮度部分，而且象素总数比较多，经过直方图均衡后，目标物的所占的灰度等级得到扩展，对比度加强，使整个图像得到增强。

数字图像均衡化后其直方图并非完全均匀分布的原因：由于原图像中目标物的灰度主要集中于低亮度部分，而且象素总数比较少，而所占的灰度等级比较多，因此图像的对比度比较好，亮度比较大，整体图像清晰。经过直方图均衡后，目标物的所占的灰度等级被压缩，对比度减弱，反而使目标物变的难以辨认。

数字图像均衡化后，其直方图并非完全均匀分布，这是因为图像的象素个数和灰度等级均为离散值；而且均衡化使灰度级并归，因此，均衡化后，其直方图并非完全均匀分布。

程序： figure(1);

fid=fopen('D:\img\cell_128.img','r'); %打开无格式文件 data1=(fread(fid,[128,128],'uint8'))'; %将打开的文件读入到data1

subplot(4,2,1); %将figure(1)分成4*2的8个子窗口, data11=uint8(data1); imshow(data11);%图象显示

title('CELL','Color','b'); %加标题

subplot(4,2,2);

title('原图像直方图');

imhist(data11);

subplot(4,2,3); %取第二个子窗口

data2=uint8(data1);%将灰度图象转换成uint8格式

b=histeq(data2); %直方图均衡化

imshow(b,256); %显示均衡化图象,256可缺省

title('均衡化','Color','b');

subplot(4,2,4);

imhist(b);

title('均衡化后图像直方图');

subplot(4,2,5)

fid=fopen('d:\img\fing_128.img','r'); %打开无格式文件

data3=(fread(fid,[128,128],'uint8'))';%将打开的文件读入到data3

data31=uint8(data3); %将灰度图象转换成uint8格式imshow(data31); %显示灰度图象

title('FING','Color','b');

subplot(4,2,6)

imhist(data31);

title('原图像直方图');

subplot(4,2,7);

data4=uint8(data3);%将灰度图象转换成uint8格式

d=histeq(data4); %直方图均衡化

imshow(d,256); %显示均衡化图象,256可缺省

title('均衡化','Color','b');

subplot(4,2,8);

imhist(d);

title('均衡化后原图像直方图');

（2）

程序：

% CELL

figure(2);

fid=fopen('D:\img\cell_128.img','r'); %打开无格式文件

data1=(fread(fid,[128,128],'uint8'))'; %将打开的文件读入到data1 I=uint8(data1);

I1=imnoise(I,'gaussian'); %加乘性噪声

H1=[0 1 0;1 0 1;0 1 0]/4; %4×4领域模板

J=imfilter(I,H1); %领域平均

subplot(2,4,1),imshow(I); %显示图像I

title('原图像');

subplot(2,4,2),imshow(I1);

title('加噪声后图像');

subplot(2,4,3),imshow(J); %不加门限平滑

title('不加门限平滑后图像');

%加门限后滤波

T= 2*sum(I1 (:))/128^2;

im_T = zeros(128,128);

for i = 1:128

for j = 1:128

if abs(I1(i,j) - J(i,j))>T

im_T(i,j) = J(i,j);

else

im_T(i,j) = I1(i,j);

end

end end

subplot(2,4,4); imshow(im_T); title('加门限后');

% FING

figure(2);

fid=fopen('D:\img\fing_128.img','r'); %打开无格式文件 data1=(fread(fid,[128,128],'uint8'))'; %将打开的文件读入到data1 I=uint8(data1);

I1=imnoise(I,'gaussian'); %加乘性噪声 H1=[0 1 0;1 0 1;0 1 0]/4; %4×4领域模板 J=imfilter(I,H1); %领域平均

subplot(2,4,5),imshow(I); %显示图像I title('原图像');

subplot(2,4,6),imshow(I1); title('加噪声后图像');

subplot(2,4,7),imshow(J); %不加门限平滑 title('不加门限平滑后图像'); %加门限后滤波 T= 2*sum(I1 (:))/128^2; im_T = zeros(128,128); for i = 1:128

for j = 1:128

if abs(I1(i,j) - J(i,j))>T im_T(i,j) = J(i,j); else

im_T(i,j) = I1(i,j); end end end

subplot(2,4,8); imshow(im_T); title('加门限后');

4. （1） 用Laplacian 锐化算子（分1=α和2=α两种情况）对256?256大小、256级灰度的数字图像lena 进行锐化处理，显示处理前、后图像。 （2） 若令

f

n m f n m g 2

1),(),(?-=α，

2(,)4(,)[(1,)(1,)(,1)(,1)

g m n f m n f m n f m n f m n f m n αα=--+++-++

)]

1,()1,(++-+n m f n m f

则回答如下问题：

① ),(n m f 、),(1n m g 和),(2n m g 之间有何关系？ ② ),(2n m g 代表图像中的哪些信息？ ③ 由此得出图像锐化的实质是什么？

解答： （1）

源代码：

%laplacian 算子锐化

I=imread('D:\img\LENA256.bmp'); % 读入原图像 figure(1);

subplot(1,3,1);imshow(I);title('原始图像'); L=fspecial('laplac ian'); L1=[0 -1 0;-1 5 -1;0 -1 0]; L2=[0 -2 0;-2 9 -2;0 -2 0];

LP1=imfilter(I,L1,'replicate');% α=1时的拉普拉斯算子 LP2=imfilter(I,L2,'replicate');% α=2时的拉普拉斯算子

subplot(1,3,2);imshow(LP1);title('Laplacian 算子α=1锐化图像'); subplot(1,3,3);imshow(LP2);title('Laplacian 算子α=2锐化图像');

（2）

①因为),(),(22n m f n m g ?-=α，所以),(n m f 、),(1n m g 和),(2n m g 之间有以下关系：

)

,(),(),(21n m g n m f n m g +=

②),(2n m g 代表了原图像中的二阶梯度信息；),(1n m g 是边缘增强后的数字图像； ③由此可以得出：图像锐化的实质是将原图像与梯度信息叠加（梯度信息所占的比例由

α的大小决定，α值越大则梯度信息所占的比例越大），相当于对目标物的边缘进行了增强。

5. 分别利用Roberts 、Prewitt 和Sobel 边缘检测算子，对256?256大小、256级灰度的数字图像lena 进行边缘检测，显示处理前、后图像。 解答：

图像：

程序：

I=imread('D:\img\LENA256.bmp'); % 读入原图像

figure(1)

%Roberts梯度法锐化

subplot(2,2,1);imshow(I);title('原始图像');

J=double(I);

[IX,IY]=gradient(J); %计算梯度

A=sqrt(IX.*IX+IY.*IY);

subplot(2,2,2);imshow(A,[]);title('Roberts梯度法锐化图像'); %Prewitt算子锐化

S=imfilter(I,fspecial('Prewitt'));

subplot(2,2,3);imshow(S);title('Prewitt算子锐化图像');

%Sobel算子锐化

S=imfilter(I,fspecial('sobel'));

subplot(2,2,4);imshow(S);title('Sobel算子锐化图像');