实验三 图像复原（源代码一站式复制粘贴）

一、实验目的

1.了解图像降质/复原处理的模型。
2.了解估计降质函数的基本原理。
3.掌握降质图像中常见噪声模型及参数估计方法、基本原理、实现步骤。
4.加深对几种常用的图像复原方法的理解。

二、实验原理

1. 图像复原的一般过程为：

2. 复原方法
   对于线性移不变系统，在空域中图像降质过程通常建模为如下卷积形式：

式中，表示观测图像，即模糊、有噪声的降质图像，表示图像的点扩散函数，表示加性噪声。若已知点扩散函数和加性噪声，就很容易反解出原图像函数。由于在空域中直接复原存在大规模计算问题，所以复原一般通过对图像进行傅里叶变换后，在频域中进行图像复原。
(1)无约束复原——维纳滤波(最小均方误差)
维纳滤波也称为最小二乘滤波，它是使原始图像与其恢复图像之间的均方误差最小的复原方法。对图像进行维纳滤波主要是为了消除图像中存在的噪声，通过此方法，维纳滤波的复原原理可表示为：

式中，为恢复图像的傅里叶变换，为降质图像的傅里叶变换。，表示噪声的功率谱，表示原图像的功率谱。功率谱函数与自相关函数是一对傅里叶变换。对于随机噪声而言，无法用确定的函数表示，也就不能用频谱来表示，在这种情况下，通常用功率谱来描述它的频域特性。
假设原图像的功率谱也为常数，上式可简化为

式中，为噪声与信号的功率比。若平稳随机过程的统计特征未知时，也可以用上式近似计算维纳滤波的传递函数，此时，近似为一个适当的常数。

(2) 约束复原——平滑约束复原(约束最小平方滤波)
约束最小平方滤波是一种比较容易实现的线性复原方法，约束复原除要求了解关于降质系统的传递函数之外，还需要知道某些噪声的统计特性或噪声与图像的某些相关情况。在最小平方约束复原中，要设法寻找一个最优估计，使得形式为的函数最小化。
在此准则下，可把图像的复原问题看作对求

目标泛化的最小值。通过最优化求解得到对应的复原图像的频域表示为：

式中，为降质函数，，表示的复共轭，为恢复图像的傅里叶变换，为降质图像的傅里叶变换，表示拉普拉斯模板的傅里叶变换—高通滤波器的传递函数：

平滑约束最小二乘复原将平滑性约束的正则项加入到复原过程中，因此，这是一种正则化滤波方法。
(3)自适应中值滤波
设表示输入图像在像素处的灰度值，表示输出图像在像素处的灰度值，表示中心在像素处的邻域像素集合，、和分别表示邻域像素集合中的最小灰度值、最大灰度值和中间灰度值，表示允许的最大模板尺寸。
自适应中值滤波包括两个阶段：阶段A和阶段B，具体的执行过程如下：
阶段A：；
；
若且，则转到阶段B
否则增大模板的尺寸
若模板的尺寸，则重复执行阶段A
否则输出。
阶段B：；
；
若且，则输出
否则输出

三、实验内容与要求

1.读入选择图像库中一幅灰度图像，对图像用“motion”及“disk”类型进行模糊降质，显示模糊前后的图像。
2.生成大小为的棋盘格式图像（自行查阅checkerboard函数的使用方法），对该图像进行模糊加噪，分别用点扩散函数、NSR、NCORR和ICORR为参数对降质图像进行恢复，显示并对比恢复结果。
3.对内容1中降质图像采用维纳滤波（deconvwnr）和最小二乘方滤波（deconvreg）的方法复原，显示复原后的图像。
4.编写自适应中值滤波adpmedfilt(g, Smax)，分析自适应中值滤波的优点。
5.从图像库中读取三幅灰度图像，对每幅依次添加椒盐噪声、高斯噪声、均匀分布噪声，观察图像的变化。然后依次用均值滤波，中值滤波、自适应中值滤波和基于局部区域统计特征的自适应滤波方法（wiener2）对噪声图像进行处理，并比较处理后的结果。通过比较总结出自适应滤波的优势和适用的滤波场合。

四、实验的具体实现

clear 、、
clc
%% 1.读入选择图像库中一幅灰度图像，对图像用“motion”及“disk”类型进行模糊降质，显示模糊前后的图像。
img = imread('Fig0217(a).tif');
%  motion类型进行模糊降质
a = fspecial('motion',40,75);%对图像模糊建模
f1 = imfilter(img,a,'circular');%图像空域滤波
% noise = imnoise(zeros(size(img)),'gaussian',0,0.001);%生成噪声
% f1=imadd(MF,im2uint8(noise));%运动模糊图像
% “disk”类型进行模糊降质
b = fspecial('disk',8);%圆形区域均值滤波
f2 = imfilter(img,b,'circular');%图像空域滤波
% noise = imnoise(zeros(size(img)),'gaussian',0,0.001);%生成噪声
% f2=imadd(MF,im2uint8(noise));
figure
subplot(1,3,1),imshow(img),title('原图');
subplot(1,3,2),imshow(f1),title('运动模糊图像');
subplot(1,3,3),imshow(f2),title('圆形区域均值滤波图');
%% 2.生成大小为120*120的棋盘格式图像对该图像进行模糊加噪，分别用点扩散函数、NSR、NCORR和ICORR为参数对降质图像进行恢复，显示并对比恢复结果
I = checkerboard(120);%生成棋盘格式图像
%imshow(I)
%对图像进行模糊加噪
noise_mean = 0;
noise_var = 0.001;
PSF = fspecial('motion',20,15);
g=imfilter(I,PSF,'conv','circular');%得到PSF的退化图片
noise = imnoise(zeros(size(I)),'gaussian',noise_mean,noise_var)%添加噪声
I1 = g+noise;
%点扩散函数恢复, 生成逆滤波图像
J1=deconvwnr(I1,PSF);

%噪信功率比估计值还原 nsr
Sn = abs(fft2(noise)).^2;%噪声功率谱
nA = sum(Sn(:))/prod(size(noise));%噪声平均功率
Sf = abs(fft2(I)).^2;%图像功率谱
fA = sum(Sf(:))/prod(size(I));%图像平均功率
r = nA / fA;%计算噪声和信号比例
J2=deconvwnr(I1,PSF,r);%使用deconvwnr函数生成常数比率的维纳滤波图像

%NCORR ICORR
NCORR=fftshift(real(ifft2(Sn)));%噪声自相关函数
ICORR=fftshift(real(ifft2(Sf)));%原图像自相关函数
J3=deconvwnr(I1,PSF,NCORR,ICORR);%使用deconvwnr函数生成自相关函数的维纳滤波图像

figure
subplot(2,2,1),imshow(I1),title('模糊图像');
subplot(2,2,2),imshow(J1),title('逆滤波图像');
subplot(2,2,3),imshow(J2),title('常数比率的维纳滤波图');
subplot(2,2,4),imshow(J3),title('自相关函数的维纳滤波图');

%% 3.对内容1中降质图像采用维纳滤波（deconvwnr）和最小二乘方滤波（deconvreg）的方法复原，显示复原后的图像
img1 = deconvwnr(f1,a);
img2 = deconvwnr(f2,b);
img3 = deconvreg(f1,a,0.1);
img4 = deconvreg(f2,b,0.1);

figure
subplot(2,2,1),imshow(img1),title('维纳滤波f1');
subplot(2,2,2),imshow(img2),title('维纳滤波f2');
subplot(2,2,3),imshow(img3),title('最小二乘方滤波f1');
subplot(2,2,4),imshow(img4),title('最小二乘方滤波f2');


%% 5
%读取三幅图像
P1 = imread('Fig0217(a).tif');
P2 = imread('Fig0219(a).tif');
P3 = imread('Fig0228(a).tif');
%对每幅依次添加椒盐噪声、高斯噪声、均匀分布噪声，观察图像的变化。
%椒盐噪声
P11 = imnoise(P1,'salt & pepper',0.05);%噪声密度0.05
P21 = imnoise(P2,'salt & pepper',0.05);
P31 = imnoise(P3,'salt & pepper',0.05);
%高斯噪声
P12=imnoise(P1,'gaussian',0,0.001);
P22=imnoise(P2,'gaussian',0,0.001);
P32=imnoise(P3,'gaussian',0,0.001);
%均匀分布噪声
P13=imnoise(P1,'speckle');
P23=imnoise(P2,'speckle');
P33=imnoise(P3,'speckle');
figure
subplot(3,3,1),imshow(P11),title('椒盐噪声P1');
subplot(3,3,2),imshow(P21),title('椒盐噪声P2');
subplot(3,3,3),imshow(P31),title('椒盐噪声P3');
subplot(3,3,4),imshow(P12),title('高斯噪声P1');
subplot(3,3,5),imshow(P22),title('高斯噪声P2');
subplot(3,3,6),imshow(P32),title('高斯噪声P3');
subplot(3,3,7),imshow(P13),title('均匀分布噪声P1');
subplot(3,3,8),imshow(P23),title('均匀分布噪声P2');
subplot(3,3,9),imshow(P33),title('均匀分布噪声P3');
%然后依次用均值滤波，中值滤波、自适应中值滤波和基于局部区域统计特征的自适应滤波方法（wiener2）对噪声图像进行处理
h=fspecial('average',3);
%P1
c1 = filter2(h,P11)/255;%对椒盐噪声图进行均值滤波
c2 = medfilt2(P11,[3,3]);%对椒盐噪声图进行中值滤波
c3 = adpmedfilt(P11,100);%对椒盐噪声图进行自适应中值滤波
c4 = wiener2(P11,[5,5]);%对椒盐噪声图进行局部区域自适应
%高斯噪声图滤波处理
c5 = filter2(h,P12)/255;
c6 = medfilt2(P12,[3,3]);
c7 = adpmedfilt(P12,100);
c8 = wiener2(P12,[5,5]);
%均匀分布噪声图滤波处理
c9 = filter2(h,P13)/255;
c10 = medfilt2(P13,[3,3]);
c11 = adpmedfilt(P13,100);
c12 = wiener2(P13,[5,5]);

figure
subplot(3,4,1),imshow(c1),title('椒盐-均值');
subplot(3,4,2),imshow(c2),title('椒盐-中值');
subplot(3,4,3),imshow(c3),title('椒盐-自适应中值');
subplot(3,4,4),imshow(c4),title('椒盐-局部自适应');
subplot(3,4,5),imshow(c5),title('高斯-均值');
subplot(3,4,6),imshow(c6),title('高斯-中值');
subplot(3,4,7),imshow(c7),title('高斯-自适应中值');
subplot(3,4,8),imshow(c8),title('高斯-局部自适应');
subplot(3,4,9),imshow(c9),title('均匀-均值');
subplot(3,4,10),imshow(c10),title('均匀-中值');
subplot(3,4,11),imshow(c11),title('均匀-自适应中值');
subplot(3,4,12),imshow(c12),title('均匀-局部自适应');

%P2
d1 = filter2(h,P21)/255;%对椒盐噪声图进行均值滤波
d2 = medfilt2(P21,[3,3]);%对椒盐噪声图进行中值滤波
d3 = adpmedfilt(P21,100);%对椒盐噪声图进行自适应中值滤波
d4 = wiener2(P21,[5,5]);%对椒盐噪声图进行局部区域自适应
%高斯噪声图滤波处理
d5 = filter2(h,P22)/255;
d6 = medfilt2(P22,[3,3]);
d7 = adpmedfilt(P22,100);
d8 = wiener2(P22,[5,5]);
%均匀分布噪声图滤波处理
d9 = filter2(h,P23)/255;
d10 = medfilt2(P23,[3,3]);
d11 = adpmedfilt(P23,100);
d12 = wiener2(P23,[5,5]);

figure
subplot(3,4,1),imshow(d1),title('椒盐-均值');
subplot(3,4,2),imshow(d2),title('椒盐-中值');
subplot(3,4,3),imshow(d3),title('椒盐-自适应中值');
subplot(3,4,4),imshow(d4),title('椒盐-局部自适应');
subplot(3,4,5),imshow(d5),title('高斯-均值');
subplot(3,4,6),imshow(d6),title('高斯-中值');
subplot(3,4,7),imshow(d7),title('高斯-自适应中值');
subplot(3,4,8),imshow(d8),title('高斯-局部自适应');
subplot(3,4,9),imshow(d9),title('均匀-均值');
subplot(3,4,10),imshow(d10),title('均匀-中值');
subplot(3,4,11),imshow(d11),title('均匀-自适应中值');
subplot(3,4,12),imshow(d12),title('均匀-局部自适应');

%P3
e1 = filter2(h,P31)/255;%对椒盐噪声图进行均值滤波
e2 = medfilt2(P31,[3,3]);%对椒盐噪声图进行中值滤波
e3 = adpmedfilt(P31,100);%对椒盐噪声图进行自适应中值滤波
e4 = wiener2(P31,[5,5]);%对椒盐噪声图进行局部区域自适应
%高斯噪声图滤波处理
e5 = filter2(h,P32)/255;
e6 = medfilt2(P32,[3,3]);
e7 = adpmedfilt(P32,100);
e8 = wiener2(P32,[5,5]);
%均匀分布噪声图滤波处理
e9 = filter2(h,P33)/255;
e10 = medfilt2(P33,[3,3]);
e11 = adpmedfilt(P33,100);
e12 = wiener2(P33,[5,5]);

figure
subplot(3,4,1),imshow(e1),title('椒盐-均值');
subplot(3,4,2),imshow(e2),title('椒盐-中值');
subplot(3,4,3),imshow(e3),title('椒盐-自适应中值');
subplot(3,4,4),imshow(e4),title('椒盐-局部自适应');
subplot(3,4,5),imshow(e5),title('高斯-均值');
subplot(3,4,6),imshow(e6),title('高斯-中值');
subplot(3,4,7),imshow(e7),title('高斯-自适应中值');
subplot(3,4,8),imshow(e8),title('高斯-局部自适应');
subplot(3,4,9),imshow(e9),title('均匀-均值');
subplot(3,4,10),imshow(e10),title('均匀-中值');
subplot(3,4,11),imshow(e11),title('均匀-自适应中值');
subplot(3,4,12),imshow(e12),title('均匀-局部自适应');