反卷积复原算法

一、Richardson-Lucy 算法
R-L 算法是目前世界上应用最广泛的函数恢复技术之一，它是一种迭代方法。

MATLAB 提供的deconvlucy （）函数还能够用于实现复杂图像重建的多种算法中，这些算法都基于Lucy-Richardson 最大化可能性算法。

R-L 算法是一种迭代非线性复原算法，它是从最大似然公式推导出来的，图像用泊松分布加以模型化的。

当下面这个迭代收敛时模型的最大似然函数就可以得到一个令人满意的方程：
1(,)(,)(,)[(,)](,)(,)
k k k g x y f x y f x y h x y h x y f x y ∧∧+∧=⊕* 其中，*代表卷积，⊕代表相关，∧f 代表未退化图像的估计，g
和h 和以前定义一样。

在IPT 中，L-R 算法由名为deconvlucy 的函数完成的。

deconvlucy()函数的调用格式：J=deconvlucy(I ，PSF ，NUMIT ，DAMPAR ，WEIGHT)。

其中，I 表示输入图像，PSF 表示点扩散函数。

其他参数都是可选参数：NUMIT 表示算法的迭代次数，默认为10次；DAMPAR 是一个标量，它指定了结果图像与原图像I 之间的偏离阈值表，默认值为0(无衰减）；WEIGHT 是一个与I 同样大小的数组，它为每一个像素分配一个权重来反映其重量，表示像素加权值，默认值为原始图像的数值。

图像复原源代码：
%% Deblurring Gray Images Using the Lucy-Richardson Algorithm
clc
clear
close all
I=imread('E:\'); % 彩色图像的像素为512*512
I1=rgb2gray(I); % 灰度图像的像素为512*512 % figure,imshow(I),title('Original color image');
% figure,imshow(I1),title('Original gray image');
I2=I1(1:2:end,1:2:end); % 图像的像素设置为256*256
figure,imshow(I2),title('Gray Image 256*256');
PSF = fspecial('gaussian',5,5); % 点扩散函数
Blurred = imfilter(I2,PSF,'symmetric','conv');
figure;
imshow(Blurred);
title('Gaussian Blurred');
V = ;
BlurredNoisy = imnoise(Blurred,'gaussian',0,V);
figure;
imshow(BlurredNoisy);
title('Blurred & Noisy');
K=size(I2);
WT=zeros(K);
WT(5:end-4,5:end-4)=1;
J1 = deconvlucy(BlurredNoisy,PSF);
% H1 = deconvlucy(BlurredNoisy,PSF,5); % 迭代5次
% H1_cell=deconvlucy({BlurredNoisy},PSF,5);
% H2_cell=deconvlucy(H1_cell,PSF);
% H2=im2uint8(H2_cell{2});
J2 = deconvlucy(BlurredNoisy,PSF,5,im2uint8(3*sqrt(V))); % 迭代5次
J3 =deconvlucy(BlurredNoisy,PSF,15,im2uint8(3*sqrt(V)));% 迭代15次
J4 =deconvlucy(BlurredNoisy,PSF,25,im2uint8(3*sqrt(V)));% 迭代25次
J5 =deconvlucy(BlurredNoisy,PSF,40,im2uint8(3*sqrt(V)));% 迭代40次
J6 =deconvlucy(BlurredNoisy,PSF,20,im2uint8(3*sqrt(V)),WT);% 迭代20次,加WT
J7 =
deconvlucy(BlurredNoisy,PSF,40,im2uint8(3*sqrt(V)),WT); % 迭代40次,加WT
%
figure, imshow(J1);
title('J1:deconvlucy(A,PSF)');
% figure, imshow(H1); title('H1:Restored Image NUMIT=5');
% figure,imshow(H2),title('H2:Restored Image NUMIT=15'); figure, imshow(J2);
title('J2:deconvlucy(A,PSF,NUMIT=5,DAMPAR)');
figure, imshow(J3);
title('J3:deconvlucy(A,PSF,NUMIT=15,DAMPAR)');
figure, imshow(J4);
title('J4:deconvlucy(A,PSF,NUMIT=25,DAMPAR)');
figure, imshow(J5);
title('J5:deconvlucy(A,PSF,NUMIT=40,DAMPAR)');
figure, imshow(J6),
title('J6:deconvlucy(A,PSF,NUMIT=20,DAMPAR,WEIGHT)'); figure, imshow(J7),
title('J7:deconvlucy(A,PSF,NUMIT=40,DAMPAR,WEIGHT)');
二、维纳滤波
维纳滤波法是由Wiener首先提出的，在图像复原领域，由于维纳滤波计算量小，复原效果好，从而得到了广泛的应用和发展。

维纳
滤波最开始主要应用在一维信号处理里，取得了比较不错的效果。

之后，维纳滤波法也用于二维信号处理中，也取得了比较好的效果。

维纳滤波器寻找一个统计误差函数：
}){(22∧
-=f f E e 最小的估计∧
f 。

E 是期望值操作符，f 是未退化的图像。

该表达式在频域可表示为
22(,)1(,)[](,)(,)(,)(,)/(,)f H u v F u v G u v H u v H u v S u v S u v η∧=+ 其中，
),(v u H 表示退化函数
),(),(),(2v u H v u H v u H *=
),(v u H *表示),(v u H 的复共轭
2
),(),(v u N v u S =η表示噪声的功率谱
2),(),(v u F v u S f =表示未退化图像的功率谱 比率(,)/(,)f S u v S u v η称为信噪功率比。

在IPT 中维纳滤波使用函数deconvwnr 来实现的。

维纳滤波能最佳复原的条件是要求已知模糊的系统函数，噪声功率谱密（或其自相关函数），原图像功率谱密度（或其自相关函数）。

但实际上，原图像功率谱密度（或其自相关函数）一般难以获知，再加上维纳滤波是将图像假设为平稳随机场的前提下的最佳滤波，而实际的图像通常不能满足此前提。

因此维纳滤波复原算法在实际中只能获得次最佳实施，它更多的是具有理论价值，被用作度量其他算法性能优劣的标杆。

维纳滤波复原函数deconvwnr()的调用格式：J=deconvwnr(I，PSF，NCORR，ICORR)
其中，I表示输入图像，PSF表示点扩散函数，NSR(默认值为0)、NCORR和ICORR都是可选参数，分别表示信噪比、噪声的自相关函数、原始图像的自相关函数。

输出参数J表示复原后的图像。

维纳滤波复原源代码：
% 维纳滤波在图像复原中的应用
clc
clear
close all
I=imread(''); % 原始图像
noise=5*randn(size(I)); % randn(1,lx)表示生成1*lx的矩阵，矩阵的每个元素都是随机数
noise=noise-min(min(noise)); % randn(size(I))是返回一个和A有同样维数大小的随机数组
J=double(I)+noise;
R1=wiener2(J,[10 10]); % 未知噪声
R2=wiener2(J,[10 10],noise); % 已知噪声分布
figure
subplot(2,2,1),imshow(uint8(I));title('原始图像');
subplot(2,2,2),imshow(uint8(J));title('退化图像');
subplot(2,2,3),imshow(uint8(R1));title('盲复原');
subplot(2,2,4),imshow(uint8(R2));title('非盲复原');
三、正则滤波
另一个线性复原的方法称为约束的最小二乘方滤波，在IPT中称为正则滤波，并且通过函数deconvreg来实现。

在最小二乘复原处理中，常常需要附加某种约束条件。

例如令Q
为f的线性算子，那么最小二乘方复原的问题可以看成使形式为
2
f
Q
的函数，服从约束条件
2
2
n
f
H
g=
-
∧
的最小化问题，这种有附加条件
的极值问题可以用拉格朗日乘数法来处理。

寻找一个∧
f，使下述准则函数为最小：
2
2
2
)
(n
f
H
g
f
Q
f
W-
-
+
=
∧
∧
∧
λ
式中λ叫拉格朗日系数。

通过指定不同的Q，可以得到不同的复原目标。

实验结果如下：
正则滤波所用的源代码：
clc
clear
close all
I=imread('E:\');
subplot(321), imshow(I),title('Original Image');
I1=rgb2gray(I);
subplot(322),imshow(I1),title('Gray Image');
PSF=fspecial('gaussian',7,10);
V=.01;
H=imfilter(I1,PSF);
BlurredNoisy=imnoise(H,'gaussian',0,V); NOISEPOWER=V*prod(size(I1));
[J, LAGRA] = deconvreg( BlurredNoisy, PSF,NOISEPOWER); K = deconvreg(BlurredNoisy, PSF,[],LAGRA/10);
K0=deconvreg(BlurredNoisy, PSF,[],LAGRA*10);
subplot(323),imshow(BlurredNoisy);
title('A=Blurred and Noisy');
subplot(324),imshow(J);
title('[J LAGRA]=deconvreg(A,PSF,NP)');
subplot(325);imshow(K);
title('deconvreg(A,PSF,[],*LAGRA)');
subplot(326);imshow(K0);
title('deconvreg(A,PSF,[],10*LAGRA)');
四、盲反卷积
在图像复原过程中，最困难的问题之一是，如何获得PSF的恰当估计。

那些不以PSF为基础的图像复原方法统称为盲去卷积。

它以MLE为基础的，即一种用被随机噪声所干扰的量进行估计的最优化策略。

工具箱通过函数deconvblind来执行盲区卷积。

实验如下：
%% 盲反卷积图像复原
clc
clear
close all
I=imread('E:\');
subplot(321),imshow(I),title('Original Image');
I=rgb2gray(I);
subplot(322),imshow(I),title('Gray Image');
PSF=fspecial('gaussian',7,10); % PSF=7x7
V=.0001;
BlurredNoisy=imnoise(imfilter(I,PSF),'gaussian',0,V);
BlurredNoisy=double(BlurredNoisy);
WT=zeros(size(I));
WT(5:end-4,5:end-4)=1; % 从第五行到倒数第五行，第五列到倒数第五列全部置为1
INITPSF=ones(size(PSF)); % INITPSF=ones(7x7)
FUN=inline('PSF+P1','PSF','P1');
[J,P]=deconvblind(BlurredNoisy,INITPSF,5,10*sqrt(V),WT,FUN, 0); % 迭代5次
[K,P]=deconvblind(BlurredNoisy,INITPSF,10,10*sqrt(V),WT,FUN ,0); % 迭代10次
[L,P]=deconvblind(BlurredNoisy,INITPSF,20,10*sqrt(V),WT,FUN ,0); % 迭代20次
subplot(323);imshow(mat2gray(BlurredNoisy));
title('A=Blurred and Noisy');
subplot(324);imshow(mat2gray(J));
title('True PSF');
subplot(325);imshow(mat2gray(K));
title('Deblured Image');
subplot(326);imshow(mat2gray(L));
title('Recovered PSF');。