【图像复原】逆滤波、维纳滤波、最小二乘方滤波（代码实现与分析）

Waygoer

已于 2024-03-04 21:31:48 修改

阅读量3k

点赞数 18

分类专栏：数字图像处理文章标签：算法人工智能图像处理计算机视觉

于 2024-03-04 15:28:15 首次发布

本文链接：https://blog.csdn.net/wyr1849089774/article/details/136452105

版权

数字图像处理专栏收录该内容

5 篇文章 1 订阅

订阅专栏

本篇博客聚焦用于图像复原技术的逆滤波、维纳滤波、最小二乘方滤波方法进行原理剖析、代码实现和结果总结，代码含有详细注释，希望帮助大家理解。

以下将从图像复原概念、逆滤波、维纳滤波、最小二乘方滤波三个块题进行阐述，每一块题将分为实验原理、实验代码展示、实验结果分析三部分进行讲解。

图像复原

图像复原的基本概念

图像复原是图像处理重要的研究领域。在成像过程中，由于成像系统各种因素的影响，可能使获得的图像不是真实景物的完善影像。图像在形成、传播和保存过程中使图像质量下降的过程，称为图像退化。图像复原就是重建退化的图像，使其最大限度恢复景物原貌的处理。

图像复原的概念与图像增强相似。但图像增强可以针对本来完善的图像，经过某一处理，使其适合于某种特定的应用，是一个主观的过程。图像复原的目的也是改善图像质量，但图像复原更偏向于利用退化过程的先验知识使已被退化的图像恢复本来面目，更多的是一个客观过程。

引起图像退化的因素包括由光学系统、运动等造成的图像模糊，以及源自电路和光学因素的噪声等。

图像复原的一般模型

一般将图像的退化过程模型化为一个退化函数和一个加性噪声项.设原始输入图像为$f(x,y)$，退化函数为 $h(x,y)$ ，加性噪声为 $n(x,y)$ ，产生的退化图像为 $g(x,y)$ ，复原滤波后重建的复原图像为 $f(x,y)$ 。退化和复原过程如下图所示。

如果系统H是一个线性、位置不变性的过程，那么在空间域中给出的退化图像可由下式给出:

$g(x, y)=h(x, y)^{*} f(x, y)+n(x, y)$

其中(*)表示空间卷积。由数字信号处理的知识可知，空域卷积在频域上可用乘积表示，因此上式等价于以下的频域表达式:

$G(u, v)=H(u, v) F(u, v)+N(u, v)$

逆滤波

实验原理

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

$\widehat{F}(u, v)=\frac{G(u, v)}{H(u, v)}$

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

$\widehat{F}(u, v)=G(u, v) + \frac{N(u, v)}{H(u, v)}$

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

实验代码实现

使用逆滤波进行图像复原时，构建逆滤波器有两种不同的方案。

方案一

该种实现方案的思路在于——对于图像而言，一般高频分量的 $H(u, v)$ 比较小，故给定截止频率 $\omega_{0}$ ，对于小于这一截止频率的低频分量进行正常的计算，而对于高于这一截止频率的高频分量，则直接赋值为1，不进行计算，以避免出现较大误差。

以下给出该部分的实现代码，该部分代码选取128，108，98，88，78，68，48七个截止频率进行处理，使用MATLAB完成实现。

I = imread("blurred_image.jpg");
subplot(2, 4, 1);
imshow(I);

radius = [128, 108, 98, 88, 78, 68, 48];
for i = 1 : length(radius)
    I = inverse_recovery(radius(i));
    subplot(2, 4, i + 1);
    imshow(I);
    title("滤波半径为" + radius(i) +"的逆滤波复原图像");
%     imwrite(I, "滤波半径为" + radius(i)+".jpg");
end


function result_image = inverse_recovery(radius)
    %读入待处理的图像
    image_d = imread("blurred_image.jpg");
    %将图像进行double类型归一化
    image_d = im2double(image_d);
    %将图像从时域转到频域
    image_d = fftshift(fft2(image_d));
    %获取频域的尺寸
    [M, N] = size(image_d);
    %预存大气湍流模型的退化函数
    [u,v]=meshgrid(1:M,1:N);%生成二维坐标系
    H=exp(-0.0025* ( (u'-M/2).^2+(v'-N/2).^2).^(5/6) );
    %定义逆滤波的滤波半径
    threshold = radius;
    if threshold > M / 2
        image_d = image_d / (H + eps);
    else 
        for i = 1 : M
            for j = 1 : N
                if sqrt((i - M / 2).^2 + (j - N / 2).^2) < threshold
                    image_d(i, j) = image_d(i, j) ./ (H(i, j) + eps);
                end
            end
        end
    end
    
    result_image = ifft2(ifftshift(image_d));
    result_image = uint8(abs(result_image)*255);
end

方案二

该种实现方法的思路在于——直接针对 $H(u, v)$ 的大小，对于较小的值则直接赋值为一个合适的 $k$ 值，对于其他的值则正常进行运算，从而避免出现较大的误差。在该式中， $k$ 和 $d$ 的值均为小于1的常数。

以下给出该部分的实现代码，该部分代码选取0.5, 0.1, 0.05, 0.03, 0.01, 0.005, 0.001七个频率值（ $d$ ）进行处理， $k$ 设置为0，使用MATLAB完成实现。

I = imread("blurred_image.jpg");
subplot(2, 4, 1);
imshow(I);

d = [0.5, 0.1, 0.05, 0.03, 0.01, 0.005, 0.001];
for i = 1 : length(d)
    I = inverse_recovery(d(i));
    subplot(2,4,i+1);
    imshow(I);
    title("设置频率范围d=" + d(i) +" k=0的逆滤波复原图像");
%     imwrite(I, "频率范围d=" + d(i) +"_k=0.jpg");
end


function result_image = inverse_recovery(d)
    %读入待处理的图像
    image_d = imread("blurred_image.jpg");
    %将图像进行double类型归一化
    image_d = im2double(image_d);
    %将图像从时域转到频域
    image_d = fftshift(fft2(image_d));
    %获取频域的尺寸
    [M, N] = size(image_d);
    %预存大气湍流模型的退化函数
    [u,v]=meshgrid(1:M,1:N);%生成二维坐标系
    H=exp(-0.0025* ( (u'-M/2).^2+(v'-N/2).^2).^(5/6) );
    %定义逆滤波的改进方式
    HC=zeros(M, N);
    M1=H>d;
    HC(M1)=1./H(M1);
    image_d = image_d.*HC;
    
    result_image = ifft2(ifftshift(image_d));
    result_image = uint8(abs(result_image)*255);
end

实验结果展示与分析

使用逆滤波进行图像处理时，对逆滤波器的构建可以有两种方式，在上一部分“代码具体实现”中已进行详细说明。

方案一

对于第一种实现方案，即

来说，处理后的实现结果如下图所示。

图像本身的高频分量较少，因此高频分量对应的 $H(u,v)$ 也较小。

可以看到，当滤波半径较大时，图像复原效果很差，甚至看不清原有的图像轮廓，这是由于滤波半径选取过大导致部分较小的H(u,v)量被纳入计算，从而使得误差加大，图像基本无法得到复原；随着半径逐渐缩小，半径选取逐渐合理，图像复原效果逐渐变优；然而，随着滤波半径的进一步缩小，信息量越来越少，导致图像逐渐模糊，图像复原效果逐渐变差。

综合选取的七种图像观察可得，滤波半径为88的逆滤波复原图像效果最佳。

方案二

对于第二种实现方案，即

来说，选取频率范围d分别为[0.5, 0.1, 0.05, 0.03, 0.01, 0.005, 0.001],替换数值k=0，所得结果如下图所示。

该方案的改进思路在于，当 $H(u,v)$ 本身数值过小时，直接将其数值替换为较小值（或0），而避免其参与 $1/ H(u,v)$ 的计算导致产生较大的误差。

可以看到，当频率范围 $d$ 设置的过小时，有些较小的数值没有被过滤掉，参与到 $1/ H(u,v)$ 的计算后会产生较大的误差，导致图像复原效果极差；当频率范围 $d$ 设置的过大时，有些反应关键信息的 $H(u,v)$ 被过滤掉，此时图像复原效果也并不理想，由图中可以看出图像比较模糊；当滤波半径选取为合适的值时，较为理想的图像复原效果可以比较直观地展现。综合选取的七种处理图像观察可得，频率范围 $d=0.03$ ， $k=0$ 的逆滤波复原图像效果最佳。

维纳滤波

实验原理

假设图像和噪声函数都是随机变量，维纳滤波方法以均分误差为评估标准，找到最小方误差下的未污染图像 $f(x,y)$ 的一个估计 $\hat{f}(x,y)$ ，这种误差度量为：

$e^{2} =E\left \{ \left ( f-\hat{f} \right )^{2} \right \}$

通常情况下，假设图像和噪声是相互独立的且至少存在一个均值为零，恢复图像和模糊图像两者灰度级呈线性关系时，维纳滤波可表示为：

$\hat{F} (u,v)=\left ( \frac{1}{H(u,v)} \frac{|H(u,v)|^{2} }{|H(u,v)|^{2}+K} \right )G(u,v)$

$|H(u,v)|^2=\bar{H} (u,v)H(u,v)$

$K=\frac{S_{n} (u,v)}{S_{f} (u,v)}$

其中， $S_{n} (u,v)$ 为噪声功率谱， $S_{f} (u,v)$ 为原图像功率谱。当信噪比较大时， $K\rightarrow 0$ ，此时维纳滤波可转换为逆滤波；当信噪比较小时， $K\rightarrow+\infty$ ，图像存在难以复原的问题。

实验代码实现

维纳滤波的图像复原操作相关代码如下所示，该部分代码由MATLAB实现。

% 读入模糊图像
blurred_image = imread('blurred_kkimage.jpg');

% 将模糊图像转换为双精度数据类型
blurred_image_double = im2double(blurred_image);

% 获取图像的尺寸
[a, b] = size(blurred_image_double);

% 将模糊图像从时域转换为频域
fft_blurred_image = fft2(blurred_image_double);
fft_blurred_image = fftshift(fft_blurred_image);

% 构建退化函数（模糊核）的频域表示
[u, v] = meshgrid(1:b, 1:a); % 注意 u 和 v 的顺序与之前不同
H = exp(-0.0025 * ((u - b/2).^2 + (v - a/2).^2).^(5/6));

% 定义维纳滤波器的参数
K = 0.01; % 维纳滤波器的常数
gamma = 0.0001; % 正则化参数

% 计算维纳滤波器的频域表示
weiner_filter = conj(H) ./ (abs(H).^2 + K);

% 应用维纳滤波器
fft_restored_image = fft_blurred_image .* weiner_filter;

% 将图像从频域转换回时域
restored_image = ifftshift(fft_restored_image);
restored_image = ifft2(restored_image);
restored_image = uint8(abs(restored_image)*256);

imshow(restored_image);
imwrite(restored_image,"restored_kkimage.jpg")

实验结果展示与分析

使用维纳滤波处理图像所得结果如下图所示。

观察图像的特性以及存在的噪声水平。如果图像受到较强的噪声干扰或模糊，较高的 $K$ 值可能会产生更好的恢复效果。而对于较低噪声水平或轻微模糊的图像，较低的 $K$ 值可能更适合。经过参数修改和对照，我们选取 $K=0.01$ 作为最终的K值，获得处理效果良好的图像。

约束最小二乘方滤波

实验原理

一个图像采集系统输入输出的关系可以表示为 $g(x,y)=H[f(x,y)]+\eta(x,y)$ ，表达为向量-矩阵形式为 $g=Hf+\eta$ ，明确地以矩阵形式来表达问题可以简化复原技术的推导。约束最小二乘方滤波的核心是解决退化函数H对噪声的敏感性问题，而减少噪声敏感性问题的一种方法是以平滑度量的最佳复原为基础的，如图像的二阶导数即拉普拉斯变换。于是，我们期望找到一个带约束条件的最小准则函数 $C$ ，定义如下：

其约束条件为：

$\|\mathbf{g}-\mathbf{H} \hat{\mathbf{f}}\|^{2}=\|\boldsymbol{\eta}\|^{2}$

其中 $\|\mathbf{w}\|^{2} \triangleq \mathbf{w}^{T} \mathbf{w}$ 是欧几里德向量范数， $\hat{f}$ 是未退化图像的估计， $\nabla^{2}$ 是拉普拉斯算子。

采用拉格朗日乘数法，在频率域中，函数 $C$ 可表示为 $C=\|P \hat{F}\|_{2}^{2}$ ，其中 $P$ 是拉普拉斯算子的傅里叶变换。则频率域中，拉格朗日函数为

$\lambda\left(\|G-H \hat{F}\|_{2}^{2}-\|N\|_{2}^{2}\right)$

其中 $N$ 是加性噪声 $\eta$ 的傅里叶变换，对 $\hat{F}$ 求导，得到 $\hat{F}$ 的最小值表达式如下：

$\hat{F}(u, v)=\left[\frac{H^{*}(u, v)}{|H(u, v)|^{2}+\gamma|P(u, v)|^{2}}\right] G(u, v)$

其中， $P(u, v)$ 是函数 $P(x, y) = \begin{bmatrix} 0 & -1 & 0\\ -1 & 4 & -1\\ 0 & -1 & 0\end{bmatrix}$ 的傅里叶变换。

实验代码实现

约束最小二乘方滤波的图像复原操作相关代码如下所示，该部分代码由MATLAB实现。

%读入图像为灰度图像
I = imread('I0.jpg');
%获取图像的长和宽
[h, w] = size(I);
figure(1);
imshow(I);
title('原始图像');
%将灰度图像从时域转为频域
I = im2double(I);%将图像的像素值转为double数据类型
fft_image = fft2(I);
fft_image = fftshift(fft_image);
%构建大气湍流模型的退化函数
[u,v]=meshgrid(1:h,1:w);%生成二维坐标系
 H=exp(-0.0025* ( (u'-h/2).^2+(v'-w/2).^2).^(5/6) );
 fft_image = fft_image.*H;
 trans_image = ifftshift(fft_image);
 trans_image = ifft2(trans_image);
 trans_image = uint8(abs(trans_image)*256);
 g=im2double(trans_image);
 figure(2);imshow(g);
 title('大气湍流模糊图像');
 % 添加加性噪声
 noise_mean = 0;
 noise_var = 0.00001;
 [h,w]=size(g);
 len = 21;
 theta = 11;
 PSF = fspecial('motion', len, theta);%产生点扩展函数
 g = imfilter(I, PSF,'conv', 'circular');
 g = imnoise(g,'gaussian',noise_mean, noise_var);
 H = psf2otf(PSF,[h,w]);
 p = [0 -1 0;-1 4 -1;0 -1 0];%拉普拉斯模板
 G=fft2(g);
 P = psf2otf(p,[h,w]);
 gamma=0.001;
 numerator = conj(H)
 denominator = H.^2+gamma*(P.^2);
 deblurred2 = ifft2(numerator.*G./denominator);%约束最小二乘方滤波在频率域中的表达式
 figure(3);
 imshow(deblurred2);
 title('约束最小二乘方滤波复原的图像')

实验结果展示与分析

使用约束最小二乘方滤波处理图像所得结果如下图所示。

在已知退化函数的情况下，对于含有不同水平噪声的运动模糊图像，最小二乘方滤波的图像复原效果应该来说都是比较好的，特别是对于高噪声和中等噪声的处理效果更好。但这一定程度上取决于 $\gamma$ 值的选择，如果 $\gamma$ （即 $\hat{F}(u, v)=\left[\frac{H^{*}(u, v)}{|H(u, v)|^{2}+\gamma|P(u, v)|^{2}}\right] G(u, v)$ 中的 $\gamma$ ）值不当，则复原的图像可能出现波纹或者较大模糊，若 $\gamma$ 值为0则简化为逆滤波。

三种滤波方法的横向比较

以下采用同一张模糊处理的图片，分别经逆滤波、维纳滤波、约束最小二乘方滤波后，将所得处理结果进行横向比较。

待处理图片如下：

处理后的结果对比：

从横向比较的图片中可以看到，三种滤波方法均能实现对模糊图像的有效复原，对于本次实验中选取的模糊处理图像来说，三种方法复原的效果相差不大。

但是综合来看，逆滤波的使用方法比较简单粗暴，需要不断调整以选取到合适的阈值进行逆滤波器的改进后才可达到比较理想的效果；对于维纳滤波和约束最小二乘方滤波来说，从实现原理上讲，两者都是有约束的，且不易受到较小量的极端影响，因此理论上的实现效果会相对更好一些，当然具体情况还需从实践中具体观察，不可一概而论。

以上为本篇博客的全部内容，希望对大家有所帮助🌹

Waygoer

关注

18
点赞
踩
59

收藏

觉得还不错? 一键收藏
0
评论
【图像复原】逆滤波、维纳滤波、最小二乘方滤波（代码实现与分析）

本篇博客聚焦用于图像复原技术的逆滤波、维纳滤波、最小二乘方滤波方法进行原理剖析、代码实现和结果总结，代码含有详细注释，希望帮助大家理解。以下将从图像复原概念、逆滤波、维纳滤波、最小二乘方滤波三个块题进行阐述，每一块题将分为实验原理、实验代码展示、实验结果分析三部分进行讲解。
复制链接

扫一扫