逆滤波和维纳滤波

简介

在图像的获取、传输以及记录保存过程中，由于各种因素，如成像设备与目标物体的相对运动，大气的湍流效应，光学系统的相差，成像系统的非线性畸变，环境的随机噪声等原因都会使图像产生一定程度的退化，图像退化的典型表现是图像出现模糊、失真，出现附加噪声等。由于图像的退化，使得最终获取的图像不再是原始图像，图像效果明显变差。为此， 要较好地显示原始图像，必须对退化后的图像进行处理，恢复出真实的原始图像，这一过程就称为图像复原。

大气湍流退化

对经过大气湍流退化的图片实现全逆滤波，半径受限逆滤波以及维纳滤波，并对比。

三种滤波思想

直接全逆滤波、半径受限逆滤波和维纳滤波

直接全逆滤波

如果退化图像为 g(x, y)，原始图像为 f (x, y)，在不考虑噪声的情况下，由傅立叶变换的卷积定理可知有下式成立 G(u, v) = H(u, v) ∗ F (u, v)。由此可见，如果已知退化图像的傅里叶变换和系统的冲击响应函数 G(u, v)，则可以求出原图像的傅里叶变换 H(u, v)，之后使用傅里叶逆变换即可得到原图像。这就是全逆滤波的主要思想。

半径受限逆滤波

全逆滤波在有噪声的情况下：

F ( u , v ) = G ( u , v ) / H ( u , v ) F (u, v) = G(u, v) /H(u, v) F(u,v)=G(u,v)/H(u,v)

F ’ ( u , v ) = G ( u , v ) / H ( u , v ) − N ( u , v ) / H ( u , v ) F ’(u, v) =G(u, v)/H(u, v) −N (u, v) /H(u, v) F’(u,v)=G(u,v)/H(u,v)−N(u,v)/H(u,v)

N (u, v) 是噪声的傅里叶变换。如果此时在 N (u, v) 不为 0 的地方，H(u, v) 过零或者比较小， 则会导致噪声被放大。因此对 F (u, v) 首先使用低通滤波器过滤掉高频的噪声，之后再除以H(u, v)，这样虽然有可能会导致图像的高频细节被忽略，但是相对于全逆滤波更加稳定。

维纳滤波

维纳滤波就是最小二乘滤波，它是使原始图像 f(x,y) 与其恢复图像 f’(x,y) 之间的均方误差最小的复原方法。具有较好的抑制噪音的能力。由于在日常使用中，无法确定图片的信噪比，因此在公式中使用参数 k 来代替信噪比。其公式为：

F ’ ( u , v ) = 1 H ( u , v ) ∗ ∣ H ( u , v ) ∣ 2 ∣ H ( u , v ) ∣ 2 + k F’(u,v) = \frac{1}{H(u,v)} * \frac{|H(u,v)|^{2}}{|H(u,v)|^{2} + k} F’(u,v)=H(u,v)1​∗∣H(u,v)∣2+k∣H(u,v)∣2​

代码1

主函数

clc;clear;close all;%% 课本图 5.28% 读取图片im = imread('demo-1.jpg'); % 原始图像 480x480 uint8%% 图像退化（大气湍流模型）% Output（H：退化模型， im_f：退化后图片）[H, im_f,im_F] = atmosph(im); %% 全逆滤波，半径受限逆滤波D0 = 60;% Input（im_f：退化图片，H：退化模型，D0：半径）% Output（im_inverse：全逆滤波结果，im_inverse_b：半径受限逆滤波）[im_inverse, im_inverse_b] = my_inverse(im_f, H, D0); % [im_inverse, im_inverse_b] = my_inverse_F(im_f, H, D0,0.001); %% 维纳滤波K = 0.0001;% Input（im_f：退化图片，H：退化模型，K：维纳滤波常数）im_wiener = my_wiener(im_f, H, K);%% 显示结果figure(1); subplot(231); imshow(im); title('原图'); axis onsubplot(232); imshow(im_f); title('大气湍流(k=0.0025)'); axis onsubplot(233); imshow(im_inverse); title('全逆滤波'); axis onsubplot(234); imshow(im_inverse_b); title('半径受限的逆滤波'); axis onsubplot(235); imshow(im_wiener); title('维纳滤波'); axis on

全逆滤波和半径受限的全逆滤波

function [im_inverse, im_inverse_b] = my_inverse(img, H, D0)[M,N] = size(img);%图像进行二位傅里叶变换，之后移动到中心点im_double = mat2gray(img,[0 255]);im_F = fftshift(fft2(im_double));% 空域 > 频域% 10阶巴特沃斯低通滤波器H2 = zeros(M,N);for x = (-M/2):1:(M/2)-1for y = (-N/2):1:(N/2)-1D2 = (x^2 + y^2)^(1/2);H2(x+(M/2)+1,y+(N/2)+1) = 1/(1+(D2/D0)^20);%十阶 endend%滤波后图像im_H2 = im_F .* H2; % 全逆滤波% 抑制低幅值的频率的全逆滤波% im_inverse = zeros(M,N);% num_F = abs(im_F);% for x = 1:M%for y = 1:N%if(num_F(x,y)>0.78)% im_inverse(x,y) = im_F(x,y) / H(x,y);%else% im_inverse(x,y) = im_F(x,y);%end %end% end%直接全逆滤波im_inverse = im_F ./ H; im_inverse_double = real(ifft2(ifftshift(im_inverse))); % 频域 > 空域im_inverse = im2uint8(mat2gray(im_inverse_double));% 半径受限逆滤波% 抑制低幅值的频率的全逆滤波% im_inverse_b = zeros(M,N);% num_H2 = abs(im_H2);% for x = 1:M%for y = 1:N%if(num_H2(x,y)>0.78)% im_inverse_b(x,y) = im_H2(x,y) / H(x,y);%else% im_inverse_b(x,y) = im_H2(x,y);%end %end% end% 直接全逆滤波im_inverse_b = im_H2 ./ H; im_inverse_b_double = real(ifft2(ifftshift(im_inverse_b))); % 频域 > 空域im_inverse_b = im2uint8(mat2gray(im_inverse_b_double));end

实验结果1

在本次实验中，首先对原图进行大气湍流退化（退化系数 k = 0.0025），之后分别对退化的图像进行全逆滤波，半径受限逆滤波（取半径 D0 = 60），维纳滤波（取信噪比 k = 0.0001）。

实验结果如下：

如上图所示，实验中全逆滤波的结果是一团噪声，但是在图像退化过程中未引入噪声。在 实验中，经过反复调试，终于发现问题出在图像退化时的傅里叶变换的时候。在图像退化的反傅里叶变换时，存在一次对图像的取实部，这个操作没有问题，但是由于图像是离散的，当图像再次傅里叶变换，会在各个频段出现看似随机的噪声。这直接导致了由于 H(u,v) 过零导致的噪声放大。

为了说明全逆滤波对无噪声图像有着极好的恢复能力，重写给定退化函数，增加返还量im_F，不做反变换，用来使用全逆滤波。新的实验结果如下图：

可以看出此时全逆滤波的效果很明显。

实验分析1

从上述实验结果可以看出，全逆滤波对无噪声的图像，具有极好的恢复能力，但是及其不稳定，可能由于各种各样的噪声，导致图像不仅没有变得清晰，反而变的更加模糊。而半径受限的全逆滤波，尽管理想的效果不如全逆滤波，会丢失很多细节，但是对于噪声的忍耐能力更强，相对也更加稳定。从上面比较来看，维纳滤波是效果最好的，并且其对噪声的抑制能力最强，但是需要知道图像的信噪比，才能最好的滤波。

运动模糊和高斯噪声污染

运动模糊是图像中的对象由于平移运动，导致的多幅图像的线性堆积。高斯噪声是服从高斯分布的亮度随机变化，是图像中常见的噪声污染。

运动模糊和高斯噪声设计思路

运动模糊是由于物体运动导致的，高斯噪声是一种由于设备、环境导致的随机性误差。

代码2

主函数

clc;clear;close all;%% 课本图 5.29% 读取图片im = imread('demo-2.jpg'); % 原始图像 688x688 uint8%% 图像退化（运动模糊+高斯噪声）% [~, im1_f] = motionblur(im, 0.01); % 噪声方差=0.01% [~, im2_f] = motionblur(im, 0.001); % [H, im3_f] = motionblur(im, 0.0000001); % H：退化模型[~, im1_f,im1_F] = motionblur(im, 0.01); % 噪声方差=0.01[~, im2_f,im2_F] = motionblur(im, 0.001); [H, im3_f,im3_F] = motionblur(im, 0.0000001); % H：退化模型%% 全逆滤波，半径受限逆滤波D0 = 33;[~, im1_inverse] = my_inverse(im1_f, H, D0);[~, im2_inverse] = my_inverse(im2_f, H, D0);[~, im3_inverse] = my_inverse(im3_f, H, D0);% [im1_inverse, ~] = my_inverse(im1_f, H, D0);% [im2_inverse, ~] = my_inverse(im2_f, H, D0);% [im3_inverse, ~] = my_inverse(im3_f, H, D0);% [~, im1_inverse] = my_inverse_F(im1_f, H, 33);% [~, im2_inverse] = my_inverse_F(im2_f, H, 33);% [~, im3_inverse] = my_inverse_F(im3_f, H, 33);% [im1_inverse, ~] = my_inverse_F(im1_f, H, D0 ,0.95);% [im2_inverse, ~] = my_inverse_F(im2_f, H, D0 ,0.31);% [im3_inverse, ~] = my_inverse_F(im3_f, H, D0 ,0.005);%% 维纳滤波%K = 0.0001;%0.1 0.008 0.0001im1_wiener = my_wiener(im1_f, H, 0.1);im2_wiener = my_wiener(im2_f, H, 0.008);im3_wiener = my_wiener(im3_f, H, 0.00001);%% 显示结果figure(1); subplot(331); imshow(im1_f); title('运动模糊+加性噪声(sigma)'); axis onsubplot(332); imshow(im1_inverse); title('逆滤波结果'); axis onsubplot(333); imshow(im1_wiener); title('维纳滤波结果'); axis onsubplot(334); imshow(im2_f); title('运动模糊+加性噪声(sigma*0.1)'); axis onsubplot(335); imshow(im2_inverse); title('逆滤波结果'); axis onsubplot(336); imshow(im2_wiener); title('维纳滤波结果'); axis onsubplot(337); imshow(im3_f); title('运动模糊+加性噪声(sigma*0.00001)'); axis onsubplot(338); imshow(im3_inverse); title('逆滤波结果'); axis onsubplot(339); imshow(im3_wiener); title('维纳滤波结果'); axis on

运动模糊

function [H,im_blured, im_blured_F] = motionblur(img, sigma)[M,N] = size(img);%图像进行二位傅里叶变换，之后移动到中心点im_double = mat2gray(img,[0 255]);im_F = fftshift(fft2(im_double));% 空域 > 频域H = zeros(M,N);gauss = 0 + sqrt(sigma) * randn(M,N);%二维高斯分布矩阵 0是均值% 运动模糊a = 0.1;b = 0.1;T = 1;for u = -M/2:1:M/2-1for v = -N/2:1:N/2-1L = ((a*u)+(b*v))*pi;if (L == 0)H(u + M/2 + 1,v + N/2 + 1) = 1;elseH(u + M/2 + 1,v + N/2 + 1) = T * sin(L) * exp( -L * 1i )/L;endendendim_G = im_F .* H;%im_double = real(ifft2(ifftshift(im_G))); % 频域 > 空域% noiseim_gauss = fftshift(fft2(gauss));% 空域 > 频域im_blured_F = im_G + im_gauss;%im_blured_F = im_G;im_blured = real(ifft2(ifftshift(im_blured_F))); % 频域 > 空域end

维纳滤波

function im_wiener = my_wiener(img, H, K)[M,N] = size(img);H2 = zeros(M,N);Mo = zeros(M,N);Me = zeros(M,N);% 空域 > 频域im_double = mat2gray(img,[0 255]);im_G = fftshift(fft2(im_double));%得到系数H2for u = 1 : Mfor v = 1 : NMo(u,v) = (abs(H(u,v)))^2;Me(u,v) = (H(u,v)*(Mo(u,v) + K));H2(u,v) = Mo(u,v)/Me(u,v);endendim_F = im_G.* H2;im_wiener_double = real(ifft2(ifftshift(im_F))); % 频域 > 空域im_wiener = im2uint8(mat2gray(im_wiener_double));end

实验结果2

首先是按照实验模板做了第一次实验，在实验中给三种情况下的维纳滤波分别调整的不同的信噪比系数 k，分别为 0.1、0.008、0.00001。但是针对逆滤波，无论是半径受限，还是未使用半径受限，都是一团噪声，无法产生需要的结果。如下图所示，（图 1 是全逆滤波，图 2是半径受限滤波）

实验思考

如上面所示，全滤波存在问题，但是按照理论知识，在噪声很小的情况下，逆滤波应该表现的相对清晰，因此，采用和大气湍流退化相同的处理方法，即对退化后的图像不进行反变换，直接全逆滤波处理，最终在半径受限的时候，得到了理想的处理结果，如下图所示：

此时，笔者在反复分析图像处理过程时，发现对于这种程度的高斯噪声，其影响的频率点的幅值不会特别高。因此，笔者写了一个低幅值抑制的函数先对图像进行处理，然后对图像使用逆滤波，最后得到了细节保留优于维纳滤波的处理算法。下面是处理的函数和算法。本次实验的下限 k 分别是0.95，0.31，0.005。

function [im_inverse, im_inverse_b] = my_inverse_F(img, H, D0, k)[M,N] = size(img);%图像进行二位傅里叶变换，之后移动到中心点im_double = mat2gray(img,[0 255]);im_F = fftshift(fft2(im_double));% 空域 > 频域% im_F = img;% 10阶巴特沃斯低通滤波器H2 = zeros(M,N);for x = (-M/2):1:(M/2)-1for y = (-N/2):1:(N/2)-1D2 = (x^2 + y^2)^(1/2);H2(x+(M/2)+1,y+(N/2)+1) = 1/(1+(D2/D0)^20);%十阶 endend%滤波后图像im_H2 = im_F .* H2; % 全逆滤波% 抑制低幅值的频率的全逆滤波im_inverse = zeros(M,N);num_F = abs(im_F);for x = 1:Mfor y = 1:Nif(num_F(x,y)>k)im_inverse(x,y) = im_F(x,y) / H(x,y);elseim_inverse(x,y) = im_F(x,y);end endend%直接全逆滤波%im_inverse = im_F ./ H; im_inverse_double = real(ifft2(ifftshift(im_inverse))); % 频域 > 空域im_inverse = im2uint8(mat2gray(im_inverse_double));% 半径受限逆滤波% 抑制低幅值的频率的全逆滤波im_inverse_b = zeros(M,N);num_H2 = abs(im_H2);for x = 1:Mfor y = 1:Nif(num_H2(x,y)>k)im_inverse_b(x,y) = im_H2(x,y) / H(x,y);elseim_inverse_b(x,y) = im_H2(x,y);end endend% 直接全逆滤波% im_inverse_b = im_H2 ./ H; im_inverse_b_double = real(ifft2(ifftshift(im_inverse_b))); % 频域 > 空域im_inverse_b = im2uint8(mat2gray(im_inverse_b_double));end

i f ( n u m F ( x , y ) > k ) 如 果 幅 值 大 于 下 限 k if (num_F (x, y) > k) 如果幅值大于下限 k if(numF​(x,y)>k)如果幅值大于下限k i n v e r s e ( x , y ) = i m F ( x , y ) / H ( x , y ) ; 正 常 使 用 逆 滤 波 inverse(x, y) = im_F (x, y)/H(x, y);正常使用逆滤波 inverse(x,y)=imF​(x,y)/H(x,y);正常使用逆滤波 e l s e 否 则 ， 说 明 其 信 息 重 要 度 低 else 否则，说明其信息重要度低 else否则，说明其信息重要度低 i n v e r s e ( x , y ) = i m F ( x , y ) ; 保 持 原 值 inverse(x, y) = im_F (x, y); 保持原值 inverse(x,y)=imF​(x,y);保持原值

实验分析2

这次实验发现，使用维纳滤波交互性的选择信噪比 K，可以得到较好的图像显示效果。如期望的，全逆滤波产生的图像质量非常差，即使使用了半径限制的逆滤波依旧会产生较大的噪声干扰。这是因为 H 的小值较多，会使噪声放大。

在实验中，通过限制低幅度值的频谱值，可以得到一定程度上优于维纳滤波的限制性逆滤波算法。但是这些滤波算法都没办法完全克服噪声的污染，随着加性的高斯噪声越来越大， 结果也越来越差。

实验小结

本次实验，研究了全逆滤波，半径受限的逆滤波，维纳滤波，以及最后通过图像的频谱图 设计出的频域幅度限制逆滤波算法。全逆滤波的抗干扰能力很差，但是如果是完全没噪声的图像，可以给出最好的修正效果。半径受限的逆滤波会损失高频的细节，但是可以有更好的抑制干扰的能力。维纳滤波要明显优于逆滤波，其通过修正信噪比 k 可以给出图像极好的恢复结果。

最后，通过图像的频谱图设计出的频域幅度限制逆滤波算法，可以通过对频谱幅值小的噪声有很好的过滤效果，尽管也会损失细节，但是可以通过设置噪声的下限来进行改善。在一定程度上，获得了比维纳滤波更好的结果。更换图片，重试发现算法有一定的普适性，具有不错的鲁棒性。