数字图像处理图像复原

数字图像处理上机实习报告

（第四专题）

学生姓名：

班 级：

学 号：

指导老师：

专题四、图像复原

一．题目要求

设已知一图象p04-03-01，设计一高通滤波器将其退化，再对退化后的图象施加椒盐噪声。利用高通滤波器的系统函数作为退化函数（高通滤波器的截止频率自行设定），采用逆滤波方式和最小二乘方滤波将其重新复原，并比较复原结果。

二.算法设计

1.逆滤波方式复原图像

2.最小二乘方滤波方式复原图像

三.程序代码与分析

1.逆滤波方式复原图像

（1）测试代码

f = imread('D:\matlab2011\work\p04-03-01.bmp');

f = im2double(f);

%经过高通滤波后的图像

[g1,G] = imhpga(f,3);

%使用高斯逆滤波还原

[ga1,GA1] = imihpga(g1,2.5);

[ga2,GA2] = imihpga(g1,3);

[ga3,GA3] = imihpga(g1,3.5);

subplot(2,3,1);imshow(abs(f)); title('原图像')

subplot(2,3,3);imshow(abs(g1)); title('高通滤波')

subplot(2,3,4);imshow(abs(ga1)); title('D0=2.5')

subplot(2,3,5);imshow(abs(ga2)); title('D0=3')

subplot(2,3,6);imshow(abs(ga3)); title('D0=3.5')

（2）指数高通滤波器作为退化函数

function [g,G] = imhpga(f,D0)

%指数高通滤波器

F = fftshift(fft2(f));

[M,N] = size(F);

m = fix(M/2);

n = fix(N/2);

for u=1:M

for v=1:N

D = sqrt((u-m)^2+(v-n)^2);

% H = 1 - exp(-D^2/(2*D0^2));

H = exp(-D0^2/(D^2));

G(u,v)=F(u,v)*H;

end

end

G = ifftshift(G);

g = ifft2(G);

End

本程序中采用指数高通滤波器，利用该高通滤波器的系统函数作为退化函数，其传递函数为：

nvuDDevuH]),([0),(

（3）采用逆滤波方式将图像复原

function [f,F] = imihpga(g,D0)

%图像复原

G = fftshift(fft2(g));

[M,N] = size(G);

m = fix(M/2);

n = fix(N/2);

for u=1:M

for v=1:N

D = sqrt((u-m)^2+(v-n)^2);

%H = 1 - exp(-D^2/(2*D0^2));

H = exp(-D0^2/(D^2));

F(u,v)=G(u,v)/(H + 0.01);

end

end

F = ifftshift(F);

f = ifft2(F);

End

图像恢复处理的关键是建立图像的退化模型。逆滤波恢复法也称反向滤波法，是一种无约束恢复的图像恢复技术。),(vuH称为系统的传递函数。

通常在无噪声的理想情况下：),(),(),(vuHvuFvuG

则),(),(),(vuHvuGvuF，对该式进行傅里叶反变换可得到),(yxf

逆滤波复原基本原理可归纳如下：

A.对退化图像),(yxg做二维离散傅里叶变换，得到),(yxG；

B.计算系统点扩散函数),(yxh的二维傅里叶变换，得到),(vuH；

C.计算),(vuF；

D.计算),(vuF的傅里叶逆变换，求得复原图像),(yxf。

2.采用最小二乘方滤波将图像复原

（1）测试代码

clc;

A=imread('D:\matlab2011\work\p04-03-01.bmp');

subplot(121),imshow(A);title('原图像');

g=double(A);

G=fft2(g);

G=fftshift(G);

subplot(122),imshow(log(abs(G)),[-1, 10]);title('原图频谱');

%%butterworth高通滤波

fn=10; %fn为滤波器阶数

D0=0.5; %截止频率为1

D1=1; %截止频率为5

D2=5; %截止频率为10

[t1,T1] =Butter_hpbw(A,D1,fn);

%维纳滤波

[h0,H0] = hpbw_weinalvbo(J0,D0,fn,0.2);

figure(2);

subplot(221),imshow(t1);title('退化后图像');

subplot(222),imshow(T1);title('退化后频谱');

subplot(223),imshow(J1);title('添加椒盐噪声后图像');

subplot(224);imshow(uint8(real(h1)));title('维纳滤波复原图像');

(2)butterworth高通滤波器作为退化函数

function [g,G] =Butter_hpbw(f,D0,fn)

%Butterworth高通滤波

%f为原图，D0为高通截止频率，fn为滤波器阶数

%返回g为滤波后图像，G为其频谱

F = fftshift(fft2(double(f)));

[M,N] = size(F);

H=zeros(M,N);

m = fix(M/2);

n = fix(N/2);

for u=1:M

for v=1:N

D=sqrt((u-m)^2+(v-n)^2);

H(u,v)= 1/(1 + (D0/D)^(2*fn) );

end

end

G=F.*H;

g= ifft2(G);

g=uint8(abs(g)); %化为易显示的格式

end

本程序中采用butterworth高通滤波器，利用该高通滤波器的系统函数作为退化函数，n阶butterworth高通滤波器的传递函数为：

nvuDDvuH20)],(/[11),(

（3）采用最小二乘方滤波将图像复原

function [g,G] = hpbw_weinalvbo(f,D0,fn,K)

%Butterworth高通维纳滤波

%f为原图，D0为高通截止频率，fn为滤波器阶数,K为噪声对信号的功率谱度比

%返回g为高通维纳滤波后图像，G为其频谱

F0=fftshift(fft2(double(f)));

[M,N] = size(F0);

m = fix(M/2);

n = fix(N/2);

H=zeros(M,N);

H0=zeros(M,N);

H1=zeros(M,N);

for u=1:m

for v=1:n

H(u,v) = 1/(1 + (D0/D)^(2*fn) );

H0(u,v)=(abs(H(u,v)))^2;

H1(u,v)=H0(u,v)/(H(u,v)*(H0(u,v)+K));

end

end

index=find(H1==0); %找出H中为0的元素的位置,很关键的操作，否则很难看出效果

H1(index)=1; %重新赋值其大小影响复原后图像的亮度

G=H1.*F0;

g=ifft2(fftshift(F2));

end

维纳滤波也叫最小二乘方滤波，是一种有约束的恢复处理方法，也是频域恢复处理中的一种，它是使原始图像与恢复图像之间的均方差最小的恢复方法。

图像恢复准则为),(yxf和),(ˆyxf之间的均方误差2e达到最小，即

})],(ˆ),({[min22yxfyxfEe

当采用线性滤波来恢复图像时，问题变为寻找点扩散函数),(yxhw，使),(*),(),(ˆyxgyxhyxfw满足图像恢复准则。

满足这一要求的传递函数为

),(),(|),(|),(),(2*vuPvuPvuHvuHvuHfnw

则有

),(),(),(|),(|),(),(ˆ2*vuGvuPvuPvuHvuHvuFfn

采用维纳滤波复原基本原理可归纳如下：

A.对退化图像),(yxg做二维离散傅里叶变换，得到),(vuG；

B.计算系统点扩散函数),(yxh的二维傅里叶变换，得到),(vuH；

C.估算图像的功率谱密度),(vuPf和噪声的谱密度),(vuPn

D.计算),(ˆvuF；

E.计算),(ˆvuF的傅里叶逆变换，求得复原图像),(ˆyxf。

四.结果分析

1.逆滤波方法复原图像

（1）原图像及其频谱