chap4
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.48y=F(x,M)×h
g 0.12 0.54 0.84 0.48
故,
原理:
约束最小二乘方原理推导代码验证 这部分的推导和代码验证,公式等都参考左侧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原理推导出来的。根据贝叶斯原理,
优缺点: 没有是否收敛的判断,好的结果可能被迭代掉了。同时,也是要已知PSF的情况下才可以使用。用同一个OTF,一个初始估计
function resim = Lucy(ifbl, LEN, THETA, iterations, handle)
%Performing Median Filter before restoring the blurred image
ifbl = medfilt2(abs(ifbl));
%Initialising the initial estimate to the blurred image
est = ifbl;
%Create PSF of degradation
PSF = fspecial('motion',LEN,THETA);
%Convert psf to otf of desired size
%OTF is Optical Transfer Function
OTF = 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 image
resim = abs(est);
自己做的RL算法,结果如下。添加了收敛情况观察。用的是下一次与上一次的2范数作为参考
clc; clear;close all
f = 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中的算法:盲去卷积原理及在图像复原中的应用
这个代码编写失败了。不知道网络上是否有研究这个,又自己写的。我写的时候出现的问题有:1.psf尺寸比g小,在乘的时候矩阵大小不一,需要处理;2.写的代码不收敛,不清楚问题在哪。 先放着在!!!!
只能知道应用deconvblind的结果,和其中的原理。但具体实现没有做成功
下面是用的MATLAB自带函数deconvblind的代码和结果:
迭代50次的结果明显不好,所以并不是迭代次数越多越好!
%% 例5.11 盲去卷积 估计PSF deconvblind
clc
clear
PSF = 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次后的结果')