反卷积复原算法

一、Richardson-Lucy 算法R-L算法是目前世界上应用最广泛的函数恢复技术之一，它是一 种迭代方法。MATLAB提供的deconvlucy ()函数还能够用于实现复 杂图像重建的多种算法中，这些算法都基于Lucy-Richardson最大化 可能性算法。R-L算法是一种迭代非线性复原算法，它是从最大似然公式推导 出来的，图像用泊松分布加以模型化的。当下面这个迭代收敛时模型 的最大似然函数就可以得到一个令人满意的方程：f (x, y) = f (x, y)y h(x, y)用*h( X, y) * fk (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 Algorithmclc clearclose allI=imread(E:); %彩色图像的像素为512*512I1=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_cell2);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次，加WTJ7 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首先提出的，在图像复原领域，由于维纳滤波计算量小，复原效果好，从而得到了广泛的应用和发展。维纳滤波最开始主要应用在一维信号处理里，取得了比较不错的效果。之 后，维纳滤波法也用于二维信号处理中，也取得了比较好的效果。维纳滤波器寻找一个统计误差函数：e 2 = E( f-f )2最小的估计f。E是期望值操作符，f是未退化的图像。该表达 式在频域可表示为1H (u, v)2H(u, v) H(u, v)2 + S (u, v) / S (u, v)F (u, v)二G (u, v)其中，H(u,v) 表示退化函数H (u, v)2 = H * (u, v)H(u, v)H * (u,v)表示H(u,v)的复共轭SJu,v)二N(u,v)|2表示噪声的功率谱S (u,v)二|F(u,v)|2表示未退化图像的功率谱比率S (u, v)/ S (u, v)称为信噪功率比。在IPT中维纳滤波使用函数 f耳deconvwnr 来实现的。维纳滤波能最佳复原的条件是要求已知模糊的系统函数，噪声功 率谱密(或其自相关函数)，原图像功率谱密度(或其自相关函数)。但实际上，原图像功率谱密度(或其自相关函数)一般难以获知，再 加上维纳滤波是将图像假设为平稳随机场的前提下的最佳滤波，而实 际的图像通常不能满足此前提。因此维纳滤波复原算法在实际中只能 获得次最佳实施，它更多的是具有理论价值，被用作度量其他算法性 能优劣的标杆。维纳滤波复原函数deconvwnr()的调用格式：J=deconvwnr(I,PSF，NCORR，ICORR)其中，I表示输入图像，PSF表示点扩散函数，NSR(默认值为0)、 NCORR和ICORR都是可选参数，分别表示信噪比、噪声的自相关函数、 原始图像的自相关函数。输出参数J表示复原后的图像。 维纳滤波复原源代码：%维纳滤波在图像复原中的应用clcclearclose allI=imread();%原始图像noise=5*randn(size(I)；% randn(l,lx)表示生成 1 *lx 的矩阵，矩阵的每个元素都是随机数noise二noisemin(min(noise)； % randn(size(I)是返回一个和 A有同样维数大小的随机数组J=double(I)+noise；R1=wiener2(J,10 10)；%未知噪声R2=wiener2(J,10 10,noise)； %已知噪声分布figuresubplot(2,2,l),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(非盲复原)； x文件（E）编辑 m（v）插入 工具桌面働窗口陛）帮助凹原贻图像理化罔豫宵复原也宵复原三、正则滤波另一个线性复原的方法称为约束的最小二乘方滤波，在IPT中称 为正则滤波，并且通过函数deconvreg来实现。在最小二乘复原处理中，常常需要附加某种约束条件。例如令Q为f的线性算子，那么最小二乘方复原的问题可以看成使形式为的函数，服从约束条件l_Hf引的最小化问题，这种有附加条件 的极值问题可以用拉格朗日乘数法来处理。寻找一个，使下述准则函数为最小：W(f)二 Qf 2 +x|g -Hf 2 - inil2式中九叫拉格朗日系数。通过指定不同的Q，可以得到不同的复 原目标。实验结果如下：些ALX文件解(E)查看M插入工具宾面口窗口 (W)帮眈(H)Original ImageGray ImageA=Blurred and NoisyJ LAGRAdecongg(AFSFNP)deconvreg(AlPSF1 t0 .VLAGRA)deconvreg(A PSF,J： 10*LAGRA)正则滤波所用的源代码：clcclear close allI=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来执行盲区卷积。实验如下：% 盲反卷积图像复原clcclearclose allI=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=7x7V=.0001;BlurredNoisy=imnoise(imfilter(I,PSF),gaussian,0,V);BlurredNoisy=double(BlurredNoisy);WT=zeros(size(I);WT(5:end-4,5:end-4)=1;% 从第五行到倒数第五行，第五列到倒数第五列全部置为1INITPSF=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);H Figure 1文件旧编轻 M(V)插入工具CD 桌面9)奮口(W)帮勛凹 “I叮傀巴哽摇日尽| 口餌I 口Original ImageGray imageA=Blurred and NoisyTrue PSFDeblured ImageRecovered PSr