数字图像处理作业

一、基础部分

1. 产生右图所示图像 f1(m,n)，其中图像大小为256×256，中间亮

条为128×32，暗处=0，亮处=100。对其进行FFT： ① 同屏显示原图f1(m,n)和FFT(f1)的幅度谱图；

② 若令f2(m,n)=（-1）m+n f1(m,n)，重复以上过程，比较二者幅度谱的异同，

简述理由；

③ 若将f2(m,n)顺时针旋转90度得到f3(m,n)，试显示FFT(f3)的幅度谱，并

与FFT(f2)的幅度谱进行比较；

④ 若将f1(m,n) 顺时针旋转90度得到f4(m,n)，令f5(m,n)=f1(m,n)+f4(m,n)，

试显示FFT(f5)的幅度谱，并指出其与 FFT(f1)和FFT(f4)的关系； ⑤ 若令f6(m,n)=f2(m,n)+f3(m,n)，试显示FFT(f6)的幅度谱，并指出其与

FFT(f2)和FFT(f3)的关系，比较FFT(f6)和FFT(f5)的幅度谱。

程序如下： clc clear all m=256; n=256; f1=zeros(m,n);

f1(65:192,113:144)=100; FFTf1=fft2(f1); ff1=log(1+abs(FFTf1));

figure(1),subplot(131),imshow(f1),title('原图像f1');

subplot(132),imshow(abs(FFTf1),[]),title('原图像的傅里叶变换'); subplot(133),imshow(abs(ff1),[]),title('以对数坐标显示'); for x=1:m for y=1:n

f2(x,y)=(-1)^(x+y)*f1(x,y);

end end

FFTf2=fft2(f2); ff2=log(1+abs(FFTf2));

figure(2),subplot(131),imshow(f2),title('原图像乘-1的m+n次方得到f2'); subplot(132),imshow(abs(FFTf2),[]),title('f2的傅里叶变换'); subplot(133),imshow(abs(ff2),[]),title('以对数坐标显示'); f3=imrotate(f2,-90); FFTf3=fft2(f3); ff3=log(1+abs(FFTf3));

figure(3),subplot(131),imshow(f3),title('f2顺时针旋转90度得到f3'); subplot(132),imshow(abs(FFTf3),[]),title('f3的傅里叶变换'); subplot(133),imshow(abs(ff3),[]),title('以对数坐标显示'); f4=imrotate(f1,-90); f5=f1+f4; FFTf4=fft2(f4); FFTf5=fft2(f5); ff4=log(1+abs(FFTf4)); ff5=log(1+abs(FFTf5));

figure(4),subplot(131),imshow(f4),title('f1顺时针旋转90度得到f4'); subplot(132),imshow(abs(FFTf3),[]),title('f4的傅里叶变换'); subplot(133),imshow(abs(ff3),[]),title('以对数坐标显示'); figure(5),subplot(131),imshow(f5),title('f5=f1+f4');

subplot(132),imshow(abs(FFTf5),[]),title('f5的傅里叶变换'); subplot(133),imshow(abs(ff5),[]),title('以对数坐标显示'); f6=f2+f3; FFTf6=fft2(f6); ff6=log(1+abs(FFTf6));

figure(6),subplot(131),imshow(f6),title('f6=f2+f3');

subplot(132),imshow(abs(FFTf6),[]),title('f6的傅里叶变换');

subplot(133),imshow(abs(ff6),[]),title('以对数坐标显示'); 结果分析：

① 同屏显示原图f1(m,n)和FFT(f1)的幅度谱图；

f1的幅度谱

② 若令f2(m,n)=(-1)^(m+n)*f1(m,n)，重复以上过程，比较二者幅度谱的异同，简述

理由；

f2的幅度谱

f1(m,n)的幅度谱中，能量集中在四个角上，f2(m,n)的幅度谱能量集中在中心位置附近。

根据傅里叶变换的频域移位特性，f1(m,n)乘(-1）^(m+n),在频域里面对应其幅度谱平移M/2和N/2，即进行了频谱中心化。

③ 若将f2(m,n)顺时针旋转90 度得到f3(m,n)，试显示FFT(f3)的幅度谱，并与FFT(f2)

的幅度谱进行比较；

f3的幅度谱

从结果可以看出f3的幅度谱相对于f2的幅度谱也顺时针旋转了90度，原因是傅里叶的旋转不变性。

④ 若将f1(m,n) 顺时针旋转90 度得到f4(m,n)，令f5(m,n)=f1(m,n)+f4(m,n)，试显示

FFT(f5)的幅度谱，并指出其与 FFT(f1)和FFT(f4)的关系；

f1的幅度谱

f4的幅度谱

f5的幅度谱

关系：FFT(f5)=FFT(f1)+FFT(f4)

⑤ 若令f6(m,n)=f2(m,n)+f3(m,n)，试显示FFT(f6)的幅度谱，并指出其与FFT(f2)和

FFT(f3)的关系，比较FFT(f6)和FFT(f5)的幅度谱。

f6的幅度谱

关系：FFT(f6)=FFT(f2)+FFT(f3)，FFT(f6)是FFT(f5)中心化后的结果。

2. 产生教材104页题图4.18（右图）所示的二值图像（白为1，黑为0），编程实现习题4.18所要求的处理（33的平均滤波和中值滤波）功能（图像四周边界不考虑，处理结果按四舍五入仍取0或1），显示处理前后的图像，比较其异同。 程序如下： clc clear all

M=zeros(256,256);%构造一个所需图片大小的全零矩阵 for i=1:2:8 for j=1:2:8

M(1+32*(i-1):32*i,1+32*(j-1):32*j)=1; M(1+32*i:32*(i+1),1+32*j:32*(j+1))=1; end end N=M;H=M; l=(1/9)*ones(3,3); for i=2:255 for j=2:255

k=M(i-1:i+1,j-1:j+1);

N(i,j)=round(sum(sum(l.*k)));

H(i,j)=median(reshape(k,1,9)); end end

figure(1);subplot(1,2,1),imshow(M),title('原始图像'); subplot(1,2,2),imshow(N),title('平滑滤波结果'); figure(2);subplot(1,2,1),imshow(M),title('原始图像'); subplot(1,2,2),imshow(H),title('中值滤波结果'); 运行结果如图所示：

3*3平滑滤波结果

中值滤波结果

从程序运行结果可以看出处理前后图像不变。

3. 产生教材104页题图4.16所示的灰度图像（白为255，黑为0），分别加入高斯白噪声和椒盐噪声，再分别进行33的平均滤波和中值滤波，显示原图像、加噪图像和滤波结果图像，并比较四种滤波结果。 程序如下： clc clear all

im=zeros(256,256); for i=0:8 a=29+24*i;

im(24:233,a:a+7)=255; end

im1=imnoise(im,'gaussian'); im2=imnoise(im,'salt & pepper'); avgim1=im1; medim1=im1; avgim2=im2; medim2=im2;

avg=1/8*[1,1,1;1,0,1;1,1,1]; for x=2:255 for y=2:255

M1=im1(x-1:x+1,y-1:y+1); M2=im2(x-1:x+1,y-1:y+1); avgim1(x,y)=sum(sum(avg.*M1)); avgim2(x,y)=sum(sum(avg.*M2)); medim1(x,y)=median(reshape(M1,1,9)); medim2(x,y)=median(reshape(M2,1,9)); end end

figure(1),subplot(131),imshow(im),title('原图像');

subplot(132),imshow(im1),title('加入高斯白噪声后的图像'); subplot(133),imshow(im2),title('加入椒盐噪声后的图像');

figure(2),subplot(221),imshow(avgim1),title('高斯噪声图像均值滤波结果'); subplot(222),imshow(medim1),title('高斯噪声图像中值滤波结果'); subplot(223),imshow(avgim2),title('椒盐噪声图像均值滤波结果'); subplot(224),imshow(medim2),title('椒盐噪声图像中值滤波结果'); 程序运行结果如下：

对原图像加入两种噪声后结果

分别采用均值滤波和中值滤波对两幅噪声图像进行滤波处理

从结果来看，均值滤波器和中值滤波器都能对噪声进行滤除，但是均值滤波会使图像边缘变得模糊，中值滤波能够很好地保持图像的边缘。均值滤波器对高斯噪声图像滤波效果较好，中值滤波对椒盐噪声图像的滤波结果好于均值滤波器对椒盐噪声图像的滤波结果。

5. 编程实现教材214页所给图像门限化分割的迭代阈值算法，实现对某一灰度图像的二值化。

程序如下： clc clear all

im=imread('花朵.jpg'); if ndims(im)==3

im=rgb2gray(im); %读取图像，如果是彩色的转化为灰度图 end

[rows cols]=size(im);

[counts,x] = imhist(im);%计算图像的灰度直方图 L=uint8(size(counts,1)); immax=max(max(im)); immin=min(min(im));

T1=(immax+immin)/2; %以图像最大灰度和最小灰度的平均值初始化T1 Imean1=uint8(sum(counts(1:T1).*(x(1:T1)))/sum(counts(1:T1))); Imean2=uint8(sum(counts(T1+1:L).*x(T1+1:L))/sum(counts(T1+1:L))); T2=(Imean1+Imean2)/2;

%T1将图像直方图分为两部分，以各部分的灰度加权平 %均值初始化T2

threshold=0; %初始化阈值，令其为0 while(1) if T1==T2

threshold=T1; %设置循环结束条件，得出最终阈值 break; else T1=T2;

Imean1=uint8(sum(counts(1:T1).*(x(1:T1)))/sum(counts(1:T1))); Imean2=uint8(sum(counts(T1+1:L).*x(T1+1:L))/sum(counts(T1+1:L))); T2=(Imean1+Imean2)/2; end end

newim=(im>=threshold); %对图像进行阈值分割，得到二值图像 figure(1),subplot(1,2,1),imshow(im),title('原灰度图像');%图像显示 subplot(1,2,2),imshow(newim),title('阈值分割图像'); 程序运行结果如下：

阈值分割二值化图像结果