python 盲去卷积算法 盲解卷积图像复原

一、图像复原与图像增强的区别

图像增强的目的是消除噪声，显现那些被模糊了的细节或简单地突出一幅图像中读者感兴趣的特征，不考虑图像质量下降的原因。图像复原是利用退化现象的某种先验知识，建立退化现象的数学模型，再根据模型进行反向的推演运算，以恢复原来的景物图像。因此图像复原可以理解为图像降质过程的反向过程。建立图像复原的反向过程的数学模型是图像复原的主要任务。

二、逆滤波复原

1、基本原理

f(x,y)表示输入图像，即理想的、没有退化的图像，g(x,y) 是退化后观察得到的图像，n(x,y)为加性噪声。通过傅立叶变换到频域后为：

图像复原的目的是给定G(u,v)和退化函数H(u,v),以及关于加性噪声的相关知识，得到原图像F(u,v)的估计图像F’(u,v),使该图像尽可能地逼近原图像F(u,v)。用于复原一幅图像的最简单的方法是构造如下的公式：

然后通过F’(u,v)的傅立叶反变换得到图像的估计值，称为逆滤波。逆滤波是一种非约束复原方法。非约束复原是指在已知退化图像G(u,v)的情况下，根据对退化模型H(u,v)和噪声N(u,v)的一些知识，做出对原图像的估计F’(u,v),使得某种事先确定的误差准则为最小。在得到误差最小的解的过程中，没有任何约束条件。对于直接逆滤波，由于存在噪声的影响，退化图像的估计公式为：

在进行逆滤波时，如果某个区域H(u,v)为0或非常小，而N(u,v)不为0且不是很小，则上式中的第2项往往比第1项大得多，从而使噪声放大，产生较大的误差。为了避免H(u,v)的值太小，可以在逆滤波时加一些限制，只在原点附近的有限邻域内进行复原，称为伪逆滤波。

2、matlab实现

% 通过逆滤波器对图像进行复原
close all;
clear all;
clc;
I=imread('cameraman.tif');
I=im2double(I);

% 构造高斯低通滤波器
s=fftshift(fft2(im2double(I)));
[a,b]=size(s);
D0=100; % 将高斯低通滤波器的截止频率D0设置为100
a0=round(a/2);
b0=round(b/2);
H=zeros(a,b);
for i=1:a
    for j=1:b
        distance=sqrt((i-a0)^2+(j-b0)^2);    % 根据高斯低通滤波器公式H(u,v)=e^-[D^2(u,v)/2*D0^2] 
        H(i,j)=exp(-(distance^2)/(2*(D0^2))); % exp表示以e为底的指数函数
    end
end

N=0.01*ones(a,b);   % 本例中a=b=256,此处亦有a=size(I,1),b=size(I,2)
N=imnoise(N,'gaussian',0,0.001); % 添加高斯噪声
J=fftfilter(I,H)+N; % 调用频域滤波函数并加入噪声

figure;
subplot(121),imshow(I);
title('原始图像');
subplot(122),imshow(J,[]); % 显示退化图像(采用高斯低通滤波器对图像进行退化,并添加均值为0,方差为0.001的高斯噪声进一步对图像进行退化)
title('退化图像');

HC=zeros(a,b);
M1=H>0.1;
HC(M1)=1./H(M1);
K=fftfilter(J,HC);% 逆滤波操作

HC=zeros(a,b);
M2=H>0.01;
HC(M2)=1./H(M2);
L=fftfilter(J,HC);% 逆滤波操作

figure;
subplot(121),imshow(K,[]);% 逆滤波复原,频率大
title('逆滤波时频率范围较大得到的图像');
subplot(122),imshow(L,[]);% 逆滤波复原,频率小
title('逆滤波时频率范围较小得到的图像');


function Z = fftfilter(X,H)
% 图像的频域滤波处理
% X为输入图像
% H为滤波器
% Z为输出图像
F=fft2(X,size(H,1),size(H,2));% 傅里叶变换
Z=H.*F;                       % 频域滤波
Z=abs(ifft2(ifftshift(Z)));   % 傅里叶反变换
Z=Z(1:size(X,1),1:size(X,2));
end

实现效果：

三、维纳滤波复原

1、基本原理

维纳滤波最早由Wiener首先提出，并应用于一维信号，取得了很好的效果。后来该算法又被引入二维信号处理，也取得相当满意的效果，尤其是在图像复原领域。由于维纳滤波器的复原效果好，计算量较低，并且抗噪性能优良，因而在图像复原领域得到了广泛的应用。

在MATLAB软件中，采用函数deconvwnr()进行图像的维纳滤波复原。该函数的调用格式如下：

2、matlab实现

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

close all;
clear all;
clc;
I=imread('onion.png');
I=rgb2gray(I);
I=im2double(I); % 注意一定要将图像类型转换为double类型,方便函数deconvwnr使用
LEN=25;
THETA=20;
% h = fspecial('motion',len,theta)返回一个过滤器,使其在与图像卷积后近似相机的线性运动
% len指定运动的长度,theta指定逆时针方向的运动角度（以度为单位）
% 过滤器将成为水平和垂直运动的矢量。默认的len是9,默认的theta是0,其对应于9个像素的水平运动
PSF=fspecial('motion',LEN,THETA);
J=imfilter(I,PSF,'conv','circular');% 表示加入运动模糊
NSR=0;
% 采用函数deconvwnr()进行图像的维纳滤波复原
% 调用格式为J=deconvwnr(I,PSF,NSR)其中PSF为点扩展函数,NSR为信噪比(电子设备中的信号与噪声的比例)
K=deconvwnr(J,PSF,NSR);

subplot(131),imshow(I);
title('原始图像');
subplot(132),imshow(J);
title('运动模糊图像');
subplot(133),imshow(K);
title('维纳滤波复原后的图像');

实现效果：

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

close all;
clear all;
clc;
I=imread('cameraman.tif');
I=im2double(I);
LEN=21;
THETA=11;
% h = fspecial('motion',len,theta)返回一个过滤器,使其在与图像卷积后近似相机的线性运动
% len指定运动的长度,theta指定逆时针方向的运动角度（以度为单位）
% 过滤器将成为水平和垂直运动的矢量。默认的len是9,默认的theta是0,其对应于9个像素的水平运动
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);% 给运动模糊图像添加均值为0,方差为0.0001的高斯噪声

figure(1);
subplot(121),imshow(I);
title('原始图像');
subplot(122),imshow(K);
title('运动模糊和添加噪声图像');

NSR1=0;
L1=deconvwnr(K,PSF,NSR1);% L1为NSR为0时得到的复原图像
NSR2=noise_var/var(I(:));% 加入的高斯噪声方差0.0001与原始灰度图像的方差0.0598之比
L2=deconvwnr(K,PSF,NSR2);% L2为采用真实NSR时得到的复原图像

figure(2);
subplot(121),imshow(L1);
title('NSR为0时的复原图像');
subplot(122),imshow(L2);
title('NSR为真实值时的复原图像');

实现效果：

（3）通过图像的自相关信息进行图像复原

close all;
clear all;
clc;
I=imread('rice.png');
I=im2double(I);
LEN=20;
THETA=10;
% h = fspecial('motion',len,theta)返回一个过滤器,使其在与图像卷积后近似相机的线性运动
% len指定运动的长度,theta指定逆时针方向的运动角度（以度为单位）
% 过滤器将成为水平和垂直运动的矢量。默认的len是9,默认的theta是0,其对应于9个像素的水平运动
PSF=fspecial('motion',LEN,THETA);
J=imfilter(I,PSF,'conv','circular');% 表示加入运动模糊

figure(1);
subplot(121),imshow(I);
title('原始图像');
subplot(122),imshow(J);
title('运动模糊图像');

noise=0.03*randn(size(I));% 返回一个n*n从标准正态分布中得到的随机项的矩阵
K=imadd(J,noise);% 可以使用函数imadd()为图像添加通用噪声(imnoise函数可以用来为图像添加特定噪声比如高斯/泊松/斑点/椒盐噪声)
% fft2-傅里叶变换、ifft2-傅里叶反变换、fftshift-将零频分量移到频谱中心、ifftshift-逆零频平移
% fft2和fftshift配合使用,ifft2和ifftshift配合使用,两个变换方向相反
% 一般都是先使用fft2或ifft2后再使用abs最后使用fftshift或ifftshift
NP=abs(fft2(noise)).^2;
NCORR=ifftshift(real(ifft2(NP)));% 噪声的自相关函数

IP=abs(fft2(I)).^2;
ICORR=ifftshift(real(ifft2(IP)));% 图像的自相关函数
% L=deconvwnr(K,PSF,NCORR,ICORR)该函数中参数NCORR为噪声的自相关函数,ICORR为原始图像的自相关函数
L=deconvwnr(K,PSF,NCORR,ICORR);

figure(2);
subplot(121),imshow(K);
title('运动模糊和添加噪声图像');
subplot(122),imshow(L);
title('维纳滤波复原后的图像');

实现效果：

四、约束最小二乘法复原

1、基本原理

在MATLAB软件中，采用函数deconvreg()进行图像的约束最小二乘法复原。该函数的调用格式如下：

2、matlab实现

（1）通过约束最下二乘法进行图像复原

close all;
clear all;
clc;
I=imread('rice.png');
I=im2double(I);
% h = fspecial('gaussian',hsize,sigma)返回大小为hsize且具有标准偏差sigma的旋转对称高斯低通滤波器。
% hsize表示过滤器大小,默认值为3*3(本例中过滤器大小为8*8),sigma为滤波器的标准值,单位为像素,默认值为 0.5
PSF=fspecial('gaussian',8,4);
% g=imfilter(f,w,filtering_mode,boundary_options,size_options)
% f为输入图像,w为滤波掩模,g为滤波后图像.filtering_mode用于指定在滤波过程中是使用“corr(相关-默认)”还是“conv(卷积)”.boundary_options用于处理边界充零问题,边界的大小由滤波器的大小确定
J=imfilter(I,PSF,'conv');

figure(1);
subplot(121),imshow(I);
title('原始图像');
subplot(122),imshow(J);% 显示退化后的图像
title('退化图像');

v=0.02;
K=imnoise(J,'gaussian',0,v);% 添加高斯噪声
NP=v*numel(I);% 返回数组元素的数目,等同于prod(size(I))(prod用于表示数组元素的乘积)
% J=deconvreg(I,PSF,NOISERPOWER)可用于图像的约束最小二乘法复原,其中PSF为扩展函数,NOISEPOWER为噪声强度,默认为0
L=deconvreg(K,PSF,NP);

figure(2);
subplot(121),imshow(K);
title('添加噪声的退化图像');
subplot(122),imshow(L);
title('约束最小二乘法复原后的图像');

实现效果：

（2）通过拉格朗日算子进行图像复原

close all;
clear all;
clc;
I=imread('rice.png');
I=im2double(I);
% h = fspecial('gaussian',hsize,sigma)返回大小为hsize且具有标准偏差sigma的旋转对称高斯低通滤波器。
% hsize表示过滤器大小,默认值为3*3(本例中过滤器大小为10*10),sigma为滤波器的标准值,单位为像素,默认值为 0.5
PSF=fspecial('gaussian',10,5);
% g=imfilter(f,w,filtering_mode,boundary_options,size_options)
% f为输入图像,w为滤波掩模,g为滤波后图像.filtering_mode用于指定在滤波过程中是使用“corr(相关-默认)”还是“conv(卷积)”.boundary_options用于处理边界充零问题,边界的大小由滤波器的大小确定
J=imfilter(I,PSF,'conv');
v=0.02;
K=imnoise(J,'gaussian',0,v);% 添加高斯噪声
NP=v*numel(I);% 返回数组元素的数目,等同于prod(size(I))(prod用于表示数组元素的乘积)
% [J,LAGRA]=deconvreg(I,PSF,...)该函数返回LAGRA为最终采用的拉格朗日算子(NP指的是噪声强度)
[L,LAGRA]=deconvreg(K,PSF,NP);
% edgetaper函数用于对图像边缘进行模糊处理
edged=edgetaper(K,PSF);

figure(1);
subplot(131),imshow(I);
title('原始图像');
subplot(132),imshow(K);% 显示退化图像
title('退化图像');
subplot(133),imshow(edged);% 显示图像边缘
title('原始图像边缘');
% J=deconvreg(I,PSF,NOISEPOWER,LRANGE)该函数中对参数LRANGE进行设置,其为拉格朗日算子搜索范围
M1=deconvreg(edged,PSF,[],LAGRA);    % 图像复原
M2=deconvreg(edged,PSF,[],LAGRA*30); % 增大拉格朗日算子
M3=deconvreg(edged,PSF,[],LAGRA/30); % 减小拉格朗日算子

figure(2);
subplot(131),imshow(M1);
title('采用拉格朗日算子复原后的图像');
subplot(132),imshow(M2);
title('增大拉格朗日算子复原后的图像');
subplot(133),imshow(M3);
title('缩小拉格朗日算子复原后的图像');

实现效果：

五、Lucy-Richardson复原

1、基本原理

在MATLAB软件中，函数deconvlucy0)采用加速收敛的Lucy-Richardson算法对图像进行复原。该函数的详细调用格式如下：

2、matlab实现

（1）对运动模糊图像采用Lucy-Richardson算法进行复原

close all;
clear all;
clc;
I=imread('rice.png');
I=im2double(I);
LEN=30;
THETA=20;
% h = fspecial('motion',len,theta)返回一个过滤器,使其在与图像卷积后近似相机的线性运动
% len指定运动的长度,theta指定逆时针方向的运动角度（以度为单位）
% 过滤器将成为水平和垂直运动的矢量。默认的len是9,默认的theta是0,其对应于9个像素的水平运动
PSF=fspecial('motion',LEN,THETA);
J=imfilter(I,PSF,'circular','conv');% 表示加入运动模糊

figure(1);
subplot(121),imshow(I);
title('原始图像');
subplot(122),imshow(J);% 显示退化图像
title('退化图像');
% 函数deconvlucy采用加速收敛的Lucy-Richardson算法对图像进行复原
% J=deconvlucy(I,PSF,NUMIT)该函数中参数NUMIT为算法的迭代次数,默认值为10
K=deconvlucy(J,PSF,5);% 复原,迭代次数为5
L=deconvlucy(J,PSF,15);% 复原,迭代次数为15(增加迭代次数后图像变得更加清晰)

figure(2);
subplot(121),imshow(K);
title('5次迭代后的图像');
subplot(122),imshow(L);
title('15次迭代后的图像');

实现效果：

（2）对高斯噪声采用Lucy-Richardson算法进行图像复原

close all;
clear all;
clc;
I=imread('cameraman.tif');
I=im2double(I);
PSF=fspecial('gaussian',7,10); % 产生PSF
v=0.0001;
J=imnoise(imfilter(I,PSF),'gaussian',0,v);% 添加高斯噪声

figure(1);
subplot(121),imshow(I);
title('原始图像');
subplot(122),imshow(J);% 显示退化图像
title('退化图像');

WT=zeros(size(I));
WT(5:end-4,5:end-4)=1;% 设置参数WEIGHT的数值
% J=deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT)中NUMIT为算法迭代次数(本例为20),DAMPAR为偏差阈值默认为0,WEIGHT为像素加权值默认为原始图像的数值
K=deconvlucy(J,PSF,20,sqrt(v));
L=deconvlucy(J,PSF,20,sqrt(v),WT);


figure(2);
subplot(121),imshow(K);
title('WEIGHT为默认值时复原后的图像');
subplot(122),imshow(L);
title('WEIGHT设置复原后的图像');

实现效果：

六、盲解卷积复原

1、基本原理

前面介绍的图像复原方法，需要预先知道退化图像的PSF。在实际应用中，经常在不知道PSF的情况下对图像进行复原。盲解卷积复原方法，不需要预先知道PSF，而且可以对PSF进行估计。盲解卷积复原算法的优点是在对退化图像无先验知识的情况下，仍然能够进行复原。

在MATLAB软件中，采用函数deconvblind()度退化的模糊图像进行盲解卷积复原。该函数的详细调用格式如下：

2、matlab实现

（1）对运动模糊图像采用盲解卷积算法进行复原

close all;
clear all;
clc;
I=imread('cameraman.tif');
I=im2double(I);
LEN=25;
THETA=20;
% h = fspecial('motion',len,theta)返回一个过滤器,使其在与图像卷积后近似相机的线性运动
% len指定运动的长度,theta指定逆时针方向的运动角度（以度为单位）
% 过滤器将成为水平和垂直运动的矢量。默认的len是9,默认的theta是0,其对应于9个像素的水平运动
PSF=fspecial('motion',LEN,THETA);
J=imfilter(I,PSF,'conv','circular');% 表示加入运动模糊
% [J,PSF]=deconvblind(I,INITPSF,NUMIT),其中输入参数INITPSF为PSF的估计值
% 输出参数PSF为算法实际采用的PSF值,NUMIT为算法的迭代次数默认为10
[K,PSF2]=deconvblind(J,PSF,30);

figure(1);
subplot(121),imshow(PSF,[]);% 显示原PSF
title('原PSF显示为图像');
subplot(122),imshow(PSF2,[]);% 显示估计的PSF
title('估计得到的PSF显示为图像');
axis auto;% axis auto表示根据x,y,z的范围自动确定坐标轴的显示范围

figure(2);
subplot(121),imshow(J);% 显示退化图像
title('退化图像');
subplot(122),imshow(K);% 显示复原图像
title('复原图像');

实现效果：

（2）% 对退化图像采用盲解卷积算法进行复原

close all;
clear all;
clc;
I=checkerboard(8);% checkerboard函数用于创建棋盘图像,I=checkerboard(n)指定棋盘图像中每个单元的边长为n个像素
PSF=fspecial('gaussian',7,10);
v=0.001;
J=imnoise(imfilter(I,PSF),'gaussian',0,v);% 添加高斯噪声
WT=zeros(size(I));
WT(5:end-4,5:end-4)=1;
% [J,PSF]=deconvblind(I,INITPSF,NUMIT,DAMPAR,WEIGHT),其中输入参数INITPSF为PSF的估计值
% 输出参数PSF为算法实际采用的PSF值,NUMIT为算法的迭代次数默认为10
% DAMPAR为偏移阈值默认值为0,WEIGHT为像素的加权值默认为原始图像的数值
[K,PFS2]=deconvblind(J,PSF,20,10*sqrt(v),WT);

subplot(131),imshow(I);
title('原始图像');
subplot(132),imshow(J);% 显示退化图像
title('退化图像');
subplot(133),imshow(K);% 显示复原图像
title('复原图像');

实现效果：