1. 实验目的
(1) 理解退化模型。
(2) 掌握常用的图像复原方法。
2. 实验内容
(1) 模拟噪声的行为和影响的能力是图像复原的核心。
(2) 空域滤波
实验一
1. 1 产生至少 2 种不同类型的噪声，并绘制原图像、加噪后图像及对应直方图于 一个图形窗口中[subplot(m,n,p)]。
1.2 对有椒盐噪声图像进行 5×5 方形窗口中值滤波。
附加内容： 自编程实现均值、中值、 自适应中值滤波器
1.1
>> f = imread('D:\ALLDOWNLOAD\实验三\实验三\lena.bmp');%读取图片到工作区
%MTALAB默认显示的图像f的原始图像
subplot(4,2,1),imshow(f),title('原始图像');
%MTALAB默认显示的图像f的直方图
subplot(4,2,2),imhist(f),title('原始图像直方图');
a=2;%均匀噪声
b=5;
noise= a + (b - a) * rand(100, 100) ;
blackIm=zeros( 100, 100);
noisedIm=noise + blackIm;
subplot(4,2,3);imshow(noisedIm,[]);title('加均匀噪声noisedIm后图像');
subplot(4,2,4 );histogram(noisedIm);title('均匀噪声noisedIm的直方图');
j = imnoise(f,'salt & pepper');%椒盐噪声
subplot(4,2,5);imshow(j,[]);title('加椒盐噪声后图像');
subplot(4,2,6 );histogram(j);title('加椒盐噪声后的直方图');
g=imnoise(f,'gaussian',0,1); %高斯噪声
subplot(4,2,7);imshow(g,[]);title('加高斯噪声后图像');
subplot(4,2,8 );histogram(g);title('加高斯噪声后的直方图');
1.2
f=imread('D:\ALLDOWNLOAD\实验三\实验三\house.bmp');
[m,n]=size(f);
subplot(4,2,1),imshow(f),title('原始图像');
j = imnoise(f,'salt & pepper',0.1);%椒盐噪声
subplot(4,2,2);imshow(j,[]);title('加椒盐噪声后图像');
w = ones(5)/(5*5); %建立一个5×5大小的均值滤波器
g1 = imfilter(j, w); %相关运算，默认边界充零
subplot(4,2,3),imshow(g1),title('加均值滤波器后图像');
g2=medfilt2(j,[5 5]);%建立大小为5x5的邻域大小的中值滤波器
subplot(4,2,4),imshow(g2),title('加中值滤波器后图像');
Smax=5;%确定最大滤波半径
u=zeros(m+2*Smax+1,n+2*Smax+1);
u(Smax+1:m+Smax,Smax+1:n+Smax)=j;
u(1:Smax,Smax+1:n+Smax)=u(1:Smax,1:n); %扩展上边界
u(1:m+Smax,n+Smax+1:n+2*Smax+1)=u(1:m+Smax,n:n+Smax); %扩展右边界
u(m+Smax+1:m+2*Smax+1,Smax+1:n+2*Smax+1)=u(m:m+Smax,Smax+1:n+2*Smax+1); %扩展下边界
u(1:m+2*Smax+1,1:Smax)=u(1:m+2*Smax+1,Smax+1:2*Smax); %扩展左边界
g3=u;
for i=Smax+1:m+Smax
for j=Smax+1:n+Smax
r=1; %初始滤波半径
while r~=Smax
W=u(i-r:i+r,j-r:j+r);
W=sort(W);
Zmin=min(W(:));
Zmax=max(W(:));
Zmed=W(uint8((2*r+1)^2/2));
if Zmin<Zmed && Zmed<Zmax %如果当前邻域中值不是噪声点，那么就用此次的邻域
break;
else
r=r+1; %否则扩大窗口，继续判断
end
end
if Zmin<u(i,j) && u(i,j)<Zmax %如果当前这个像素不是噪声，原值输出
g3(i,j)=u(i,j);
else %否则输出邻域中值
g3(i,j)=Zmed;
end
end
end
subplot(4,2,5);imshow(g3(Smax+1:m+Smax,Smax+1:n+Smax),[]);title('加入自适应中值滤波后图像');
实验二
2.2 选择一幅清晰的灰度图像，对该图像进行模糊化处理并加入不同强度的高斯 噪声，然后分别采用逆滤波、维纳滤波和约束最小二乘方滤波对退化图像进行复 原，比较各种图像复原方法的复原效果。
附加内容： 自编程实现不同截至频率的逆滤波并比较其复原效果。
- 逆滤波（高斯噪声）
%%读入图片
I = im2double(imread('D:\ALLDOWNLOAD\实验三\实验三\house.bmp'));% [0,1]
[M,~] = size(I);
figure;
subplot(2,3,1), imshow(I);
title('原始图像');
%模拟运动模糊H(u,v)
T=1;a=0.02;b=0.02;
v=[-M/2:M/2-1];u=v';
A=repmat(a.*u,1,M)+repmat(b.*v,M,1);
H=T/pi./A.*sin(pi.*A).*exp(-1i*pi.*A);
H(A==0)=T;% replace NAN
%得到模糊图像
F=fftshift(fft2(I));
FBlurred=F.*H;
%显示模糊图像
IBlurred =real(ifft2(ifftshift(FBlurred)));
subplot(2,3,2), imshow(uint8(255.*mat2gray(IBlurred)));
title('运动模糊图像');
%% 无噪声情况下的复原
FDeblurred=FBlurred./H;
IDeblurred=real(ifft2(ifftshift(FDeblurred)));
subplot(2,3,4), imshow(uint8(255.*mat2gray(IDeblurred)));
title('无噪情况下直接逆滤波');
%% Simulate Noise Model
noise_mean = 0;
noise_var = 1e-3;
noise=imnoise(zeros(M),'gaussian', noise_mean,noise_var);
FNoise=fftshift(fft2(noise));
%得到加噪的运动模糊图像
FBlurred_Noised=FNoise+FBlurred;
IBlurred_Noised=real(ifft2(ifftshift(FBlurred_Noised)));
subplot(2,3,3), imshow(uint8(255.*mat2gray(IBlurred_Noised)));
title('加噪运动模糊图像');
%直接逆滤波
FDeblurred2=FBlurred_Noised./H;
FH1=abs(FDeblurred2);
IDeblurred2=real(ifft2(ifftshift(FDeblurred2)));
subplot(2,3,5), imshow(uint8(255.*mat2gray(IDeblurred2)));
title ('有噪情况下直接逆滤波');
%设置截至半径为Radius的复原
Radius=33;
FDeblurred2=zeros(M);
for a=1:M
for b=1:M
if sqrt((a-M/2).^2+(b-M/2).^2)<Radius
FDeblurred2(a,b)= FBlurred_Noised(a,b)./H(a,b);
end
end
end
FH2=abs(FDeblurred2);
IDeblurred2=real(ifft2(ifftshift(FDeblurred2)));
subplot(2,3,6), imshow(uint8(255.*mat2gray(IDeblurred2)));
title(['加噪运动模糊图像半径为', num2str(Radius),'的圆内逆滤波复原']);
2.维纳滤波
%读取图片到工作区
I = im2double(imread('D:\ALLDOWNLOAD\实验三\实验三\house.bmp'));
subplot(3,2,1);imshow(I);title('Original Image (courtesy of MIT)');%原始图像
%模拟一个运动模糊，仿真相机移动信号
LEN = 21;
THETA = 11;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');
subplot(3,2,2);imshow(blurred);title('Blurred Image')%运动模糊图像
%仿真噪声，假定噪声为 0 ，恢复图像
noise_mean = 0;
noise_var = 0.0001;
blurred_noisy = imnoise(blurred, 'gaussian', noise_mean, noise_var);
subplot(3,2,3);imshow(blurred_noisy);title('Simulate Blur and Noise')%仿真噪声图像
%无噪声假设的条件下重构图像，使用最佳噪声/信号功率比恢复图像
estimated_nsr = 0;
wnr2 = deconvwnr(blurred_noisy, PSF, estimated_nsr);
subplot(3,2,4);imshow(wnr2);title('Restoration of Blurred, Noisy Image Using NSR = 0')%无噪声重构图像
%进行合理的噪声假设重构图像
estimated_nsr = noise_var / var(I(:));
wnr3 = deconvwnr(blurred_noisy, PSF, estimated_nsr);
subplot(3,2,5);imshow(wnr3);title('Restoration of Blurred, Noisy Image Using Estimated NSR');%合理噪声假设重构图像
3.约束最小二乘方滤波
% 读取图片到工作区
I = im2double(imread('D:\ALLDOWNLOAD\实验三\实验三\house.bmp'));
[hei,wid,~] = size(I);
subplot(2,3,1),imshow(I);
title('Original Image (courtesy of MIT)');
% 模拟一个运动模糊
LEN = 21;
THETA = 11;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');
subplot(2,3,2), imshow(blurred); title('Blurred Image');
% 反滤波器百分比
If = fft2(blurred);
Pf = psf2otf(PSF,[hei,wid]);
deblurred = ifft2(If./Pf);
subplot(2,3,3), imshow(deblurred); title('Restore Image')
% 模拟加性噪声的百分比
noise_mean = 0;
noise_var = 0.00001;
blurred_noisy = imnoise(blurred, 'gaussian', ...
noise_mean, noise_var);
subplot(2,3,4), imshow(blurred_noisy);title('Simulate Blur and Noise')
%尝试使用自制的约束最小二乘过滤器进行恢复
p = [0 -1 0;-1 4 -1;0 -1 0];
P = psf2otf(p,[hei,wid]);
gama = 0.001;
If = fft2(blurred_noisy);
numerator = conj(Pf);
denominator = Pf.^2 + gama*(P.^2);
deblurred2 = ifft2( numerator.*If./ denominator );
subplot(2,3,5), imshow(deblurred2);
title('Restoration of Blurred Using Constrained Least Squares Filtering');%使用约束最小二乘过滤恢复模糊
subplot(2,3,6); imshow(deconvreg(blurred_noisy, PSF,0)); title('Regul in Matlab');
