chap4 图像复原 part4

来源：互联网发布：gulp 合并js 编辑：程序博客网时间：2024/05/17 02:20

chap4

1. 约束最小二乘方滤波

用向量来表达卷积: gMN×1=HMN×MNfMN×1+ηMN×1 . 这个表达可以参考《反卷积和信号复原》

例子： x=[0.4,1,0.8]′,y=[0.3,0.6]′ 计算 x与 y 的卷积
x 0.4 1 0.8 这个过程可以用卷积不卷来解释。h(t)偏移 tao，乘以 x 在 t
y 0.3 0.6 tao 点的值，再求和。如此，取遍所有的 tao。“叠加原理”
0.12 0.3 0.24 左侧这个过程用矩阵来表达：
o.24 0.6 0.48 y=F(x,M)×h
g 0.12 0.54 0.84 0.48

⎡⎣⎢⎢⎢0.120.540.840.48⎤⎦⎥⎥⎥=⎡⎣⎢⎢⎢0.30.60.30.60.30.6⎤⎦⎥⎥⎥⎡⎣⎢0.410.8⎤⎦⎥ 但是，这是完全卷积的情况，即是：Ny=Nx+Nh−1 。对于图像处理等实际情况，通常是部分卷积，比如二维图像，卷积后的结果与原始图像是一样的大小。这时仅仅是取其中的一部分出来。原理如下：
⎡⎣⎢⎢⎢⎢⎢⎢y0y1y2...yL−1⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢h0h1h2...hM−1h0h1...hM−2hM−1h0...hM−3hM−2...............h0h1h2...hM−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢x0x1x2...xN−1⎤⎦⎥⎥⎥⎥⎥⎥ 这是完全卷积，从中取出一部分，变为：

⎡⎣⎢⎢⎢⎢⎢⎢y0y1y2...yL−1⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢hM−1hM−2hM−1hM−3hM−2............h1h2...hM−1h0h1...hM−2hM−1h0...hM−3hM−1h0......h0h1h0⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢x0x1x2...xN−1⎤⎦⎥⎥⎥⎥⎥⎥

故，y=FTh^×h . 若是将g=Hf+η 写成如上部分卷积的形式，其中g=MN×1,H=MN×MN,η=MN×MN

原理：

约束最小二乘方原理推导代码验证 这部分的推导和代码验证，公式等都参考左侧blog即可

2. 基于Lucy-Richardson算法的迭代非线性复原

逆滤波、维纳滤波、约束最小二乘法复原都是线性复原方法，L-R是非线性复原。
MATLAB自带函数的用法 J = deconvlucy(I,PSF,numit,dampar,weight,readout,subsmpl)
其中，numit最大迭代次数；
dampar 指定输出图像与输入图像的阈值偏差，低于此阈值出现阻尼。对于dampar的值低于原始值的像素，迭代被抑制。这抑制了这样的像素产生噪声，保留了其他地方的图像细节。
weight 权重反映了相机的录制质量。可以根据平场校正的数量来调整权重。
readout 对应于加性噪声和读出的相机噪声的方差。
subsmpl 子采样，当PSF在比图像更精细的subsmpl的网格上给出时使用。

J = DECONVLUCY(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT,SUBSMPL)；

J-复原图像；I-输入的模糊图像；PSF-点扩展函数；NUMIT-迭代次数，缺省值为10；

DAMPAR-偏差阈值，缺省值为0。图像通过迭代计算后，尤其在低信噪比的情况下，复原图像可能会出现一些斑点。这些斑点并不代表图像的真实信息，而是输出图像过于逼近噪声所产生的结果。为了控制这种现象的产生，采样DAMPAR作为收敛参数，参数指定收敛过程中复原图像与原始输入图像背离程度的阈值。对于那些超过阈值的数据，在反复迭代过程中将被禁止。

WEIGHT-像素记录能力，缺省值为原始图像的数值；可能包含坏像素的数据可能会随时间和位置的变化而变化。通过指定WEIGHT数组参数，可以忽略图像中某些指定的像素，需要忽略的像素对应的WEIGHT数组元素值为0.

READOUT-噪声矩阵，缺省值为0。CCD摄像机获得数字图像的过程中噪声主要有两部分组成：一个是呈泊松分布的光子计算噪声；而另一个则是镜头呈高斯分布的读取噪声。Lucy-Richardson算法计算过程中，需要自己声明第二种噪声情况。READOUT即是用来指定这种读取噪声。其参数值通常是读取噪声变量和背景噪声变量的总和，数值的大小将指定一个能够确保所有数值为正数的偏移量。

SUBSMPL-子采样时间，缺省值为1。如果将采样不足的数据的复原过程，建立在一个好的网格操作基础上，则可以较大程度的提高重建效果。SUBSMPL来指定采样不足的比例。如果数据采样不足是由于图像获取过程中的镜头像素装仓？？问题产生的，那么每个像素观察到的PSF就是一个好的网格PSF。
Matlab中的L-R算法的步骤：
(1) PSF的FFT OTF尺寸与I一样大。
(2) 准备迭代参数

Reference：
[1] William Hadley Richardson, ‘Bayesian-Based Iterative Method of Image Restoration’, 1972, Journal of OSA.
[2]Lucy

算法原理及代码实现

References: pudn王教团CV作业，逆滤波、维纳滤波、R-L算法代码参考
原理： [陈波博士学位论文《自适应光学图像复原理论与算法研究》chap6.2.1]
RL算法是Poisson噪声模型下的极大似然估计图像复原方法，是从Bayes原理推导出来的。根据贝叶斯原理，p(O|I)=p(I|O)p(O)p(I) . 根据Helmholtz原理，最佳猜测原理通过最优化能量泛函实现。不同的最优化能量泛函的不同会导致使用的方法不同。MAP是指最大化p(I|O)p(O)时的目标O，ML是指最大化p(I|O)时的目标O。又根据假定图像噪声是Possion噪声或Gaussian噪声的不同分为：GML、GMAP、PML、PMAP。还有Myopic算法。其中，解决PML图像复原的方法就是：R-L算法。也称为EM算法（Expectation Maximization） .
f^k+1=f^k⋅F−1[H⋆F[gF−1(HF^k)]] 主要的公式是这一个。参照如下代码理解：
优缺点： 没有是否收敛的判断，好的结果可能被迭代掉了。同时，也是要已知PSF的情况下才可以使用。用同一个OTF，一个初始估计f^0然后迭代，直到达到最大迭代次数为止。

function resim = Lucy(ifbl, LEN, THETA, iterations, handle)%Performing Median Filter before restoring the blurred imageifbl = medfilt2(abs(ifbl));%Initialising the initial estimate to the blurred imageest = ifbl;%Create PSF of degradationPSF = fspecial('motion',LEN,THETA);%Convert psf to otf of desired size%OTF is Optical Transfer FunctionOTF = psf2otf(PSF,size(ifbl));i = 1;while i<=iterations    %Converting the estimate to frequency domain    fest = fft2(est);    %Multiplying OTF with the estimate in frequency domain    fblest = OTF.*fest;    %Converting the blurred image estimate to spatial domain    ifblest = ifft2(fblest);    %Calculating ratio of blurred image and estimate of the deblurred image    iratio = ifbl./ifblest;    %Converting the ratio to frequency domain    firatio = fft2(iratio);    %Calculating the correction vector    corrvec = OTF .* firatio;    %Converting the correction vector to spatial domain    icorrvec = ifft2(corrvec);    %Multiplying correction vector & estimate of deblurred image to find next estimate    aftercorr = icorrvec.*est;    est = aftercorr;    %Setting the waitbar indicating how much is completed    waitbar(i/iterations, handle);    i = i+1;end%Restored imageresim = abs(est);

自己做的RL算法，结果如下。添加了收敛情况观察。用的是下一次与上一次的2范数作为参考
这里写图片描述

clc; clear;close allf = checkerboard(8);%用棋盘图像来测试PSF = fspecial('motion',7,45);%用运动模糊来模拟gb = imfilter(f, PSF,'circular');noise = imnoise(zeros(size(f)),'gaussian',0,0.001);%加高斯噪声g = gb + noise;%图像+模糊+高斯噪声[J1, R1] = MyRL(g, PSF, 10);[J2, R2] = MyRL(g, PSF, 20);[J3, R3] = MyRL(g, PSF, 200);subplot(3,2,1);imshow(J1);title('迭代10次结果')subplot(3,2,2);plot(1:10,R1);title('收敛情况')subplot(3,2,3);imshow(J2);title('迭代20次结果')subplot(3,2,4);plot(1:20,R2);title('收敛情况')subplot(3,2,5);imshow(J3);title('迭代200次结果')subplot(3,2,6);plot(1:200,R3);title('收敛情况')%RL算法f_k+1 = f_k* F'(H* F(g/(F'(H*F_k)))function [J RstErr] = MyRL(g,PSF,NUMIT)%g—输入退化图像 PSF—输入模糊函数 NUMIT 迭代的最大次数est = g;%初始估计为f0本身OTF = psf2otf(PSF,size(g));%OTF 大小与图像一致，填充for i = 1: NUMIT    Fest = fft2(est);%估计的原始图像    Fblest = Fest .* OTF;%估计的模糊图像    ifblest = ifft2(Fblest);%空域的估计图像    iratio = g ./ ifblest;%比例    Fratio = fft2(iratio);    icorrec = ifft2(conj(OTF) .* Fratio);%    nextest = est .* icorrec;    RstErr(i) = norm(nextest - est);    est = nextest;    i = i+1;end    J = abs(est);end

3. 盲解卷积

[f,PSFe] = deconvblind(g, INITPSF, NUMIT, DAMPAR, WEIGHT) MATLAB自带盲解卷积函数的使用形式。其中，里面是使用R-L算法进行迭代复原。输入初始PSF的估计。其中，PSF估计受其初始推测尺寸的影响十分大，而其初始值的大小对复原结果影响不大。
参见blog中的算法：盲去卷积原理及在图像复原中的应用
f^k+1=f^k⋅F−1[H⋆kF[gF−1(H^kF^k)]] 这是在频域中的RL公式，当进行盲解卷积时，我们通过初始PSF和初始f去计算g_k+1,然后用g_k+1和f又去计算f_k+1.如此迭代即可。那么在上述公式中，只需要变换相应的参数即可。
h^k+1=h^k⋅F−1[F⋆k−1Fourier[gF−1(H^kF^k−1)]]
这个代码编写失败了。不知道网络上是否有研究这个，又自己写的。我写的时候出现的问题有：1.psf尺寸比g小，在乘的时候矩阵大小不一，需要处理；2.写的代码不收敛，不清楚问题在哪。先放着在！！！！

只能知道应用deconvblind的结果，和其中的原理。但具体实现没有做成功

下面是用的MATLAB自带函数deconvblind的代码和结果:
迭代50次的结果明显不好，所以并不是迭代次数越多越好！
这里写图片描述

%% 例5.11 盲去卷积 估计PSF deconvblindclcclearPSF = fspecial('gaussian',7,10);subplot(3,2,1);imshow(pixeldup(PSF,73),[]);title('原始PSF图像')f = checkerboard(8);SD = 0.01;g = imnoise(imfilter(f,PSF),'gaussian',0,SD^2);INITPSF = ones(size(PSF)); % 初始值(尺寸大小与原始PSF一样)DAMPAR = 10*SD;LIM = ceil(size(PSF,1)/2);WEIGHT = zeros(size(g));WEIGHT(LIM + 1:end - LIM, LIM + 1:end - LIM) = 1;NUMIT = 5;[fr PSFe] = deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);subplot(3,2,2);imshow(pixeldup(PSFe,73),[])title('使用盲去卷积估计PSF迭代5次后的结果')NUMIT = 10;[fr PSFe] = deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);subplot(3,2,3);imshow(pixeldup(PSFe,73),[])title('使用盲去卷积估计PSF迭代10次后的结果')NUMIT = 20;[fr PSFe] = deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);subplot(3,2,4);imshow(pixeldup(PSFe,73),[])title('使用盲去卷积估计PSF迭代20次后的结果')NUMIT = 50; % 并非迭代次数越多越好！！！[fr PSFe] = deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);subplot(3,2,5);imshow(pixeldup(PSFe,73),[])title('使用盲去卷积估计PSF迭代50次后的结果')

阅读全文

0 0