图像复原

1图像复原的而理论模型

定义：在成像过程中，由于成像系统各种因素的影响，可能使获得的图像不是真实景物的完善影像。图像在形成、传播和保存过程中使图像质量下降的过程，称为图像退化。图像复原就是重建退化的图像，使其最大限度恢复景物原貌的处理。图像复原只能尽量使图像接近其原始图像，但由于噪声干扰等因素，很难精确还原。

图像增强与图像复原的区别：图像增强的目的是消除噪声，显现那些被模糊了的细节或简单地突出一幅图像中读者感兴趣的特征，，不考虑图像质量下降的原因。图像复原是利用退化现象的某种先验知识，建立退化现象的数学模型，再根据模型进行反向的推演运算，以恢复原来的景物图像。因而，图像复原可以理解为图像降质过程的反向过程。建立图像复原的反向过程的数学模型，就是图像复原的主要任务。经过反向过程的数学模型的运算，要想恢复全真的景物图像比较困难。所以，图像复原本身往往需要有一个质量标准，即衡量接近全真景物图像的程度，或者说，对原图像的估计是否到达最佳的程度。

2噪声模型：噪声主要来源于图像的获取和传输过程。

（1）图像传感器的工作情况受各种因素的影响，如图像获取中的环境条件和传感器元器件自身的质量。

（2）图像在传输过程中主要由于所用传输信道被干扰而受到噪声污染。

噪声种类：高斯噪声、瑞利噪声、伽马噪声、指数噪声、均匀分布噪声、脉冲噪声（椒盐噪声）

在matlab命令窗口输入命令绘制概率密度图：

>> x=-4:.1:4;
>> subplot(321)
>> Y1=show_noise_pdf('gaussian',x,0,1);
>> plot(x,Y1);
>> title('高斯');
>> subplot(322)
>> Y2=show_noise_pdf('uniform',x,-3,3);
>> plot(x,Y2);
>> title('均匀');
>> subplot(323)
>> Y3=show_noise_pdf('salt & pepper',x);
>> plot(x,Y3);
>> title('椒盐');
>> subplot(324)
>> Y4=show_noise_pdf('rayleigh',x,1);
>> plot(x,Y4);
>> title('瑞利');
>> subplot(325)
>> Y5=show_noise_pdf('exp',x,1);
>> plot(x,Y5);
>> title('指数');
>> subplot(326)
>> Y6=show_noise_pdf('gamma',x,2,5);
>> plot(x,Y6);
>> title('伽马');

在matlab命令窗口中调用add_noise函数为图像添加噪声，并调用matlab自带函数hist绘制灰度直方图

>>I=imread('square.bmp');
>> J1=add_noise(I,'gaussian',0,10);
>> subplot(321)
>> hist(double(J1(:)),100)
>> title('高斯');
>> subplot(321)
>> J2=add_noise(I,'uniform',-20,20);
>> hist(double(J2(:)),100)
>> title('均匀');
>> subplot(323)
>> J3=add_noise(I,'salt & pepper',0.02);

>> hist(double(J3(:)),100)
>> title('椒盐');
>> subplot(324)
>> J4=add_noise(I,'rayleigh',30);
>> hist(double(J4(:)),100)
>> title('瑞利');
>> subplot(325)
>> J5=add_noise(I,'exp',15);
>> hist(double(J5(:)),100)
>> title('指数');
>> subplot(326)
>> J6=add_noise(I,'gamma',2,10);
>> hist(double(J6(:)),100)
>> title('伽马');

3空间滤波

空间滤波器种类：均值滤波器（适于高斯噪声）、统计排序滤波器（适于椒盐噪声）、自适应局部噪声消除滤波器（适于所有噪声）、自适应中值滤波器（适于椒盐噪声）。

在图像中加入0.05的高斯噪声，用中值滤波和自适应滤波进行复原：

>> I=imread('lena.bmp');
>> I0=imnoise(I,'salt & pepper',0.01);
>> I1=medfilt2(I0,[3,3]);
>> I2=adp_median(I0,7);
>> subplot(221);
>> imshow(I);
>> title('原图');
>> subplot(222);
>> imshow(I0);
>> title('椒盐噪声')
>> subplot(223);
>> imshow(I1);
>> title('中值滤波')
>> subplot(224);
>> imshow(I2);
>> title('自适应滤波');

4逆滤波复原

（1）退化

%读取原始图像
I=imread('lena.bmp');
figure(1);
subplot(121)
imshow(I)
title('原始图像');
f=im2double(I);

%傅里叶变换
F=fft2(f);
F=fftshift(F);

%执行退化
[M,N]=size(F);
[u,v]=meshgrid(1:M,1:N);
H=exp(-0.0025*((u-M/2).^2+(v-N/2).^2).^(5/6));
F=F.*H;

%傅里叶变换
X=ifftshift(F);
x=ifft2(X);
subplot(122)
x=uint8(abs(x)*256);
imshow(x)
imwrite(x,'lena_t.bmp');
title('退化图像');

（2）复原

%逆滤波
I0=imread('lena_t.bmp');

%阈值为128
I_new1=rev_filter(I0,H,128);

%阈值为108
I_new2=rev_filter(I0,H,108);

%阈值为78
I_new3=rev_filter(I0,H,78);

%阈值为48
I_new4=rev_filter(I0,H,48);
si=zeros(M,N,1,4,'uint8');
si(:,:,1)=I_new1;
si(:,:,2)=I_new2;
si(:,:,3)=I_new3;
si(:,:,4)=I_new4;

%绘图
figure
montage(si)
title('阈值分别为128,108,78,48')

5维纳滤波复原

定义：维纳滤波只能解决退化函数，没有加性噪声的问题，维纳滤波又称最小均方误差滤波，综合考虑了退化函数和噪声。

对256*256的lena灰度图像进行退化处理，再添加高斯噪声，对得到的退化图像进行复原步骤如下：

（1）用imnoise（）函数添加均值为0，方差为0.001的高斯噪声

（2）将逆滤波和维纳滤波复原效果进行对比。

%维纳滤波与逆滤波
clear;
clc;

%%
I=imread('lena.bmp');

%傅里叶变换
[m,n]=size(I);
FI=fft2(I);
FI=fftshift(FI);

%退化
k=0.0025;
u=1:m;
v=1:n;
[u,v]=meshgrid(u,v);
H=exp((-k).*(((u-m/2).^2+(v-n/2).^2).^(5/6)));
G=FI.*H;

%添加噪声
I0=real(ifft2(fftshift(G)));
I1=imnoise(uint8(I0),'gaussian',0,0.001);
figure(1)
imshow(I1);
imwrite(I1,'lena_wn.bmp')

%%
I1=imread('lena_wn.bmp');

%%逆滤波
I_new=rev_filter(I1,H,48);
figure(2);
imshow(I_new)
title('逆滤波结果');

%维纳滤波
K=0.05;
I_new1=wn_filter(I1,H,48,K);
figure(3);
imshow(I_new1);
title('维纳滤波结果');

从图中可以看出，取同样的半径时，维纳滤波所得图像比逆滤波消除噪声的效果更高

6有约束最小二乘复原

%约束最小二乘复原
%%
clear,clc
close all

%%产生退化图像
I=checkerboard(8);

%运动模糊的点扩散函数
PSF=fspecial('motion',7,45);

fprintf('点扩散函数:\n');
disp(PSF)

%对图像进行运动模糊滤波
Im1=imfilter(I,PSF,'circular');

%添加高斯噪声
noise=imnoise(zeros(size(I)),'gaussian',0,0.001);
Im=Im1+noise;
%%维纳滤波
Iw=deconvwnr(Im,PSF,0.02);

%约束最小二乘滤波
I=edgetaper(I,PSF);
Iz=deconvreg(Im,PSF,0.2,[1e-7,1e7]);

%%绘图
subplot(221);
imshow(I,[])
title('原始图像');

subplot(222);
imshow(Im,[])
title('退化图像');

subplot(223)
imshow(Iw,[])
title('维纳滤波');

subplot(224)
imshow(Iz,[])
title('最小二乘滤波');

7Lucky-Richardson复原

%run lucy
%L-R算法图像复原

%%
clear,clc
close all

%%
%棋盘格图像
I=checkerboard(8);

%点扩散函数
PSF=fspecial('gaussian',7,10);

%方差为0.0001
SD=0.01;
In=imnoise(imfilter(I,PSF),'gaussian',0,SD^2);

%%
%使用Lucy-Richardson算法对图像复原
Dampar=10*SD;
LIM=ceil(size(PSF,1)/2);
Weight=zeros(size(In));

%权值weight数组的大小是64*64
%并且有值为0的4像素宽的边界，其余像素都是1
Weight(LIM+1:end-LIM,LIM+1:end-LIM)=1;

%迭代次数为5
NumIt=5;
%利用deconvlucy来复原
J1=deconvlucy(In,PSF,NumIt,Dampar,Weight);

%迭代次数为10
NumIt=10;
J2=deconvlucy(In,PSF,NumIt,Dampar,Weight);

%迭代次数为20
NumIt=20;
J3=deconvlucy(In,PSF,NumIt,Dampar,Weight);

%迭代次数为100
NumIt=100;
J4=deconvlucy(In,PSF,NumIt,Dampar,Weight);

%%绘图
subplot(231);
imshow(I);
title('原图');

subplot(232);
imshow(In);
title('退化图像');

subplot(233);
imshow(J1);
title('迭代5次');

subplot(234);
imshow(J2);
title('迭代10次');

subplot(235);
imshow(J3);
title('迭代20次');

subplot(236);
imshow(J4);
title('迭代100次');

8盲去卷积图像复原

%run_blind
%盲去卷积复原
%%
clear,clc
close all

%产生棋盘格图像并进行退化
I=checkerboard(8);
PSF=fspecial('gaussian',7,10);
V=.0001;

%退化
BlurredNoisy=imnoise(imfilter(I,PSF),'gaussian',0,V);

%%复原
%权值
Weight=zeros(size(I));
Weight(5:end-4,5:end-4)=1;

%点扩散函数的估计值
InitPSF=ones(size(PSF));
[J, P]=deconvblind(BlurredNoisy,InitPSF,20,10*sqrt(V),Weight);
subplot(221);
imshow(BlurredNoisy);
title('退化图像');
subplot(222);
imshow(PSF,[]);
title('点扩展图像');
subplot(223);
imshow(J);
title('盲去卷积复原结果');
subplot(224);
imshow(P,[]);
title('输出的估计点扩展函数');