图像复原 opencv 图像复原算法

图像复原的常用方法主要包括：逆滤波复原、维纳滤波复原、约束最小二乘法复原、Lucy-Richardson复原和盲解卷积复原。

逆滤波复原

f(x,y)表示输入图像，即为理想没有退化的图像，g(x,y)时退化后观察得到的图像，n(x,y)为加性噪声。通过傅里叶变换到频域后为：G(u,v) = F(u,v)H(u,v)+N(u,v)

图像复原的目的是给定G(u,v)和退化函数H(u,v)，以及关于加性噪声的相关知识，得到原图像F(u,v)的估计图像F’(u,v)，使得图像尽可能地逼近原图像F(u,v)。复原一幅图像的最简单的方法为：F’(u,v) = G(u,v)/H(u,v)，然后通过F’(u,v)的傅里叶反变换得到图像的估计值，称为逆滤波。逆滤波是一种非约束复原方法，非约束复原是指在以知退化图像G(u,v)的情况下，根据对退化模型H(u,v)和噪声N(u,v)的一些知识，做出对原图像的估计F’(u,v)，使得某种事先确定得误差准则为嘴小。在得到误差最小的过程中，没有任何约束条件。

直接逆滤波时：由于存在噪声，所以退化图像公式 F’(u,v) = G(u,v)+N(u,v)/H(u,v)，为了避免H(u,v)的值太小，可以在逆滤波时附加一些限制，只在原点附近的有跟邻域内进行复原，称为伪逆滤波

逆滤波器对图像进行复原：

%建立函数fftfilter()进行图像的频域滤波
function Z=fftifilter(X,H)
F=fft2(X,size(H,1),size(H,2));%傅里叶变换
Z = H.*F;%滤波
Z=ifftshift(Z);
Z=abs(ifft2(Z));%傅里叶反变换
Z=Z(1:size(X,1),1:size(X,2));
end

>> I = imread('E:\persional\matlab\images\cameraman.tif');
>> I = im2double(I);
>> [m,n] = size(I);
>> M=2*m;N=2*n;
>> u=-m/2:m/2-1;
>> v=-n/2:n/2-1;
>> [U,V] = meshgrid(u,v);
>> D=sqrt(U.^2+V.^2);
>> D0=130;%截止频率
>> H=exp(-(D.^2)./(2*(D0^2)));%高斯低通滤波器
>> N=0.01*ones(size(I,1),size(I,2));
>> N=imnoise(N,'gaussian',0,0.01);%添加噪声
>> J=fftifilter(I,H)+N;%频域滤波并加噪声
>> figure,
>> subplot(221),imshow(I);
>> subplot(222),imshow(J,[ ]);
>> HC=zeros(m,n);
>> M1=H>0.1;%频率范围
>> HC(M1)=1./H(M1);
>> K=fftifilter(J,HC);%逆滤波
>> HC=zeros(m,n);
>> M2=H>0.01;
>> HC(M2)=1./H(M2);
>> L=fftifilter(J,HC);
>> subplot(223),imshow(K,[ ]);
>> subplot(224),imshow(L,[ ]);

维纳滤波复原

J = deconvwnr(I,PSF,NSR)
 PSF：点扩展函数
 NSR：信噪比，信号均值与背景标准偏差的比值J = deconvwnr(I,PSF,NCORR,ICORR)
 NCORR：噪声的自相关函数
 ICORR：原始图像的自相关函数

通过维纳滤波对运动模糊图像进行复原：

>> I = imread('E:\persional\matlab\images\lena.bmp');
>> I = rgb2gray(I);
>> I = im2double(I);
>> LEN = 25;%运动位移像素
>> THETA = 20;%运动角度
>> PSF = fspecial('motion',LEN,THETA);%产生PSF
>> J = imfilter(I,PSF,'conv','circular');%运动模糊
>> NSR = 0;
>> K = deconvwnr(J,PSF,NSR);%维纳滤波复原
>> figure,
>> subplot(131),imshow(I);
>> subplot(132),imshow(J);
>> subplot(133),imshow(K);

通过维纳滤波对含有噪声的运动模糊图像进行复原：

>> I = imread('E:\persional\matlab\images\lena.bmp');
>> I = im2double(I);
>> LEN = 21;%运动位移像素
>> THETA = 11;%运动角度
>> PSF = fspecial('motion',LEN,THETA);
>> J = imfilter(I,PSF,'conv','circular');%通过卷积对图像进行滤波，产生模糊图像
>> noise_mean=0;%均值
>> noise_var=0.0001;%方差
>> K = imnoise(J,'gaussian',noise_mean,noise_var);%添加高斯噪声
>> NSR1=0;
>> L1 = deconvwnr(K,PSF,NSR1);%维纳滤波复原
>> NSR2=noise_var/var(I(:));
>> L2=deconvwnr(K,PSF,NSR2);%图像复原
>> figure,
>> subplot(221),imshow(I);
>> subplot(222),imshow(K);
>> subplot(223),imshow(L1);%NSR为0的复原图像
>> subplot(224),imshow(L2);%采用真实的NSR时得到复原的图像

通过图像的自相关信息进行复原：

>> I = imread('E:\persional\matlab\images\lena.bmp');
>> I = im2double(I);
>> LEN = 20;%运动位移像素
>> THETA=10;%运动角度
>> PSF=fspecial('motion',LEN,THETA);%产生PSF
>> J = imfilter(I,PSF,'conv','circular');%运动模糊
>> noise=0.03*randn(size(I));
>> K = imadd(J,noise);%添加噪声
>> NP = abs(fft2(noise)).^2;
>> NPPower=sum(NP(:))/prod(size(noise));
>> NCORR=fftshift(real(ifft2(NP)));%噪声自相关函数
>> IP=abs(fft2(I)).^2;
>> IPower=sum(IP(:))/prod(size(I));
>> ICORR=fftshift(real(ifft2(IP)));%函数自相关函数
>> L=deconvwnr(K,PSF,NCORR,ICORR);%图像复原
>> figure,
>> subplot(221),imshow(I);%原图
>> subplot(222),imshow(J);%退化图像
>> subplot(223),imshow(K);%噪声图像
>> subplot(224),imshow(L);%图像复原

约束最小二乘法复原

J = deconvreg(I,PSF)
 J = deconvreg(I,PSF,NOISEPOWER)
 J = deconvregl(I,PSF,NOISEPOWER,LRANGE)
 J = deconvregl(I,PSF,NOISEPOWER,LRANGE,REGOP)
 [J,LAGRA] = deconvreg(IPSF,…)
 PSF：点扩展函数
 NOISEPOWER：噪声的强度，默认为0
 LRANGE：拉格朗日算子的搜索范围，默认为[10-9,109]
 REGOP：约束算子
 LAGRA：LAGRA为最终采用的拉格朗日算子

约束最小二乘法进行图像复原

>> I = imread('E:\persional\matlab\images\lena.bmp');
>> I = im2double(I);
>> PSF=fspecial('gaussian',8,4);
>> J = imfilter(I,PSF,'conv');%图像退化
>> v=0.04;
>> K = imnoise(J,'gaussian',0,v);%添加噪声
>> NP = v*prod(size(I));
>> L = deconvreg(K,PSF,NP);%图像复原
>> figure,
>> subplot(221),imshow(I);
>> subplot(222),imshow(J);%退化图像
>> subplot(223),imshow(K);%添加噪声
>> subplot(224),imshow(L);

通过拉格朗日算子进行图像复原：

>> I = imread('E:\persional\matlab\images\lena.bmp');
>> I = im2double(I);
>> PSF=fspecial('gaussian',12,6);
>> J = imfilter(I,PSF,'conv');
>> v=0.04;
>> K = imnoise(J,'gaussian',0,v);
>> NP = v*prod(size(I));
>> [L,LAGRA]=deconvreg(K,PSF,NP);%图像复原
>> edged=edgetaper(K,PSF);%提取边缘
>> M1=deconvreg(edged,PSF,[],LAGRA)%图像复原;
>> M2=deconvreg(edged,PSF,[],LAGRA*30);%增大拉格朗日算子
>> M3=deconvreg(edged,PSF,[],LAGRA/30);%减小拉格朗日算子
>> figure,
>> subplot(231),imshow(I);
>> subplot(232),imshow(K);
>> subplot(233),imshow(edged);
>> subplot(234),imshow(M1);
>> subplot(235),imshow(M2);
>> subplot(236),imshow(M3);

Lucy-Richardson 复原

deconvlucy()采用加速收敛的Lucy-Richardson 算法对图像进行复原
J = deconvlucy(I,PSF)
 J = deconvlucy(I,PSF,NUMIT)
 J = deconvlucy(I,PSF,NUMIT,DAMPAR)
 J = deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT)
 J = deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT)
 J = deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT,SUBSMPL)
 PSF：扩展函数
 NUMIT：算法重复的次数，默认值为10
 DAMPAR：偏差阈值，默认值为0
 WEIGHT：像素的加权值，默认为原始图像的初始值
 READOUT：噪声矩阵，默认值为0
 SUBSMPL：采样时间，默认值为1

运动模糊图像：

>> I = imread('E:\persional\matlab\images\lena.bmp');
>> I = im2double(I);
>> LEN = 30;
>> THETA = 20;
>> PSF = fspecial('motion',LEN,THETA);
>> J = imfilter(I,PSF,'circular','conv');%图像退化
>> K = deconvlucy(J,PSF,5);%复原，5次迭代
>> L = deconvlucy(J,PSF,15);%复原，15次迭代
>> figure,
>> subplot(221),imshow(I);
>> subplot(222),imshow(J);%图形退化
>> subplot(223),imshow(K);
>> subplot(224),imshow(L);

高斯噪声：

>> I = imread('E:\persional\matlab\images\lena.bmp');
>> I = im2double(I);
>> PSF = fspecial('gaussian',7,10);
>> v = 0.0001；
>> J = imnoise(imfilter(I,PSF),'gaussian',0,v);%高斯噪声
>> WT = zeros(size(I));
>> WT(5:end-4,5,5:end-4);%像素权重
>> K = deconvlucy(J,PSF,20,sqrt(v));%图像复原
>> L = deconvlucy(J,PSF,20,sqrt(v),WT);
>> figure,
>> subplot(221),imshow(I);
>> subplot(222),imshow(J);
>> subplot(223),imshow(K);
>> subplot(224),imshow(L);

盲解卷积复原

前面的图像复原方法，需要预先知道退化图像的PSF，盲解卷积复原方法，不需要知道PSF，而对PSF进行估计。盲解卷积复原算法的优点是在对退化图像无先验知识的情况下，仍然能够进行复原

[J,PSF] = deconvblind(I,INITPSF)
 [J,PSF] = deconvblind(I,INITPSF,NUMIT)
 [J,PSF] = deconvblind(I,INITPSF,NUMIT,DAMPAR)
 [J,PSF] = deconvblind(I,INITPSF,NUMIT,DAMPAR,WEIGHT)
 [J,PSF] = deconvblind(I,INITPSF,NUMIT,DAMPAR,WEIGHT,READOUT)
 INITPSF：PSF输入参数INITPSF的估计值
 NUMIT：重复的次数
 DAMPAR：偏移阈值，默认值为0
 WEIGHT：像素加权值，默认为原始图像的数值
 READOUT：噪声矩阵

运动模糊图像：

>> I = imread('E:\persional\matlab\images\lena.bmp');
>> I = im2double(I);
>> LEN=20;
>> THETA=20;
>> PSF=fspecial('motion',LEN,THETA);
>> J = imfilter(I,PSF,'circular','conv');%运动模糊
>> INITPSF=ones(size(PSF));%输入估计值
>> [K,PSF2] = deconvblind(J,INITPSF,30);%图像复原
>> figure,
>> subplot(221),imshow(PSF,[]);%显示PSF
>> subplot(222),imshow(PSF2,[]);%显示估计的PSF
>> subplot(223),imshow(J);%退化图像
>> subplot(224),imshow(K);%复原图像

对退化图像进行复原：

>> I = checkerboard(8);
>> PSF = fspecial('gaussian',7,10);
>> v = 0.001;
>> J = imnoise(imfilter(I,PSF),'gaussian',0,v);%添加高斯噪声
>> INITPPSF = ones(size(PSF));%输入估计值
>> WT = zeros(size(I));
>> WT(5:end-4,5:end-4)=1;:
>> [K,PSF2] = deconvblind(J,INITPPSF,20,10*sqrt(v),WT);%图像复原
>> figure,
>> subplot(131),imshow(I);
>> subplot(132),imshow(J);
>> subplot(133),imshow(K);