图像处理第二次作业(2)

对一幅灰度图像：

(1)计算并画出此图像的中心化频率谱。

(2)分别用高斯低通和高斯高通滤波器对图像进行频域处理，并对结果进行分析。

(3)用频域拉普拉斯算子对此图像进行锐化处理，并对结果进行分析。

(1)程序

clear;

clc;

A=imread('lena.jpg','jpg');

B=rgb2gray(A);

figure;

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

B_d=double(B);

[m,n] = size(B);

%fft2函数实现

if ismatrix(B_d)%ismatrix:确定输入的是否为矩阵

if nargin==1%nargin用来判断输入变量的个数

f = fftn(B_d);%fftn:N维离散傅里叶变换

else

f = fftn(B_d,[m n]);

end

else

if nargin==1

f = fft(fft(B_d,[],2),[],1);

else

f = fft(fft(B_d,n,2),m,1);

end

end

C3=log(1+abs(f));

subplot(2,2,2);imshow(C3,[]);title('原图像频谱');

%计算并画出此图像的中心化频率谱

[m,n] = size(B);

for i=1:m

for j=1:n

C1(i,j)=((-1)^(i+j))*(B_d(i,j));

end

end

C2=log(1+abs(fft2(C1)));

C4=fftshift(log(1+abs(fft2(B_d))));

subplot(2,2,3);imshow(C2,[]);

title('编写函数实现频谱中心化')

subplot(2,2,4);imshow(C4,[])

title('matlab函数实现频谱中心化')

%分别用低通滤波和高通滤波对此图像进行频域处理

C5=fftshift(fft2(B_d));

P1=round(m/2);

Q1=round(n/2);

D0=30;%半径取30,随着D0的增加，低通滤波图像越来越模糊B_lf=C5;

B_hf=C5;

%低通滤波器

for i=1:m

for j=1:n

D=sqrt((i-P1)^2+(j-Q1)^2);

if D<=D0

flag=1;

else

flag=0;

end

B_lf(i,j)=flag*C5(i,j);

end

end

B_lf=uint8(real(ifft2(ifftshift(B_lf))));

figure;

%高通滤波器

for i=1:m

for j=1:n

D=sqrt((i-P1)^2+(j-Q1)^2);

if D>D0

flag=1;

else

flag=0;

end

B_hf(i,j)=flag*C5(i,j);

end

end

B_hf=uint8(real(ifft2(ifftshift(B_hf))));

subplot(1,3,1);imshow(B);title('原图像');

subplot(1,3,2);imshow(B_lf);

title('高斯低通滤波');

subplot(1,3,3);imshow(B_hf);

title('高斯高通滤波');

%用频域拉普拉斯算子对此图像进行锐化处理

for i=1:m

for j=1:n

D1=(i-P1)^2+(j-Q1)^2;

H(i,j)=(-4)*((pi)^2)*D1;

H1(i,j)=H(i,j)*C5(i,j);

end

end

H2=real(ifft2(ifftshift(H1)));

H2max=max(H2(:));

H3=H2(1:m,1:n);

k=1;

H4=B_d-k*H3/H2max;

figure;

subplot(1,2,1);imshow(B);title('原图像'); subplot(1,2,2); imshow(H4,[]);

title('拉普拉斯图像');

%使用moon图像

Task3=imread('moon.tif','tif');

Task3_1=double(Task3);

Task3_2=fftshift(fft2(Task3_1));

[m1,n1]=size(Task3);

P2=round(m1/2);

Q2=round(n1/2);

for i=1:m1

for j=1:n1

D1=(i-P2)^2+(j-Q2)^2;

H(i,j)=(-4)*((pi)^2)*D1;

H1(i,j)=H(i,j)*Task3_2(i,j);

end

end

H2=real(ifft2(ifftshift(H1)));

H2max=max(H2(:));

H3=H2(1:m1,1:n1);

k=1;

H4=Task3_1-k*H3/H2max;

figure;

subplot(1,2,1);imshow(Task3);title('原图像'); subplot(1,2,2); imshow(H4,[]);

title('拉普拉斯图像');

(2)实验结果