2015-03-22 20:39:49 ispforfun 阅读数 10115

  维纳滤波器是一种自适应的滤波器,在数字信号处理中有着广泛的应用。ispforfun会在从今天开始定期给大家带来维纳滤波器在图像处理中应用。本节讲诉维纳滤波器在图像去噪中的简单应用。

       让我们从Matlab中的函数wiener2开始。Matlab的help中对wiener2的说明如下:

       wiener2估计图像中每个像素的局部均值和局部方差,


      其中,是图像A中每个像素的NXM邻域。维纳滤波器的去噪原理如下:


      其中是噪声的方差,如果没有给出来,那么会根据图像的局部方差进行估计。

     

     ispforfun根据这段原理,编写了如下的matlab代码,共享给大家。

function J = wiener_image_denoising_filter(I, Win, noise_var)
% J = wiener_image_denoising_filter(I, w)
% -----------------------------------------
% by: ispforfun
% I is the input image. Single channel image.
% Win is the window size for filtering.

[img_height, img_width] = size(I);
half_win_size = floor(Win/2.0);

[X, Y] = meshgrid(1:img_width, 1:img_height);
[dX, dY] = meshgrid(-half_win_size(2):half_win_size(2), -half_win_size(1):half_win_size(1));

dX = reshape(dX, [1 1 Win(1) Win(2)]);
dY = reshape(dY, [1 1 Win(1) Win(2)]);
X = repmat(X, [1 1 Win(1) Win(2)]) + repmat(dX, [img_height img_width 1 1]);
Y = repmat(Y, [1 1 Win(1) Win(2)]) + repmat(dY, [img_height img_width 1 1]);

X(X<1) = 2 - X(X<1);
X(X>img_width) = 2*img_width - X(X>img_width);
Y(Y<1) = 2 - Y(Y<1);
Y(Y>img_height) = 2*img_height - Y(Y>img_height);

patch = @(f)f(Y + (X-1)*img_height);

P = patch(I);

img_mean = local_image_mean(P);
img_var = local_image_var(P, img_mean, Win);

if nargin < 3
    estimated_noise_var = sum(sum(img_var))/(img_height * img_width);
else
    estimated_noise_var = noise_var;
end

gain = (estimated_noise_var == 0.0).*zeros(img_height, img_width) + ...
    (estimated_noise_var ~= 0.0).*((img_var - estimated_noise_var)./img_var);

J = img_mean + gain.*(double(I) - img_mean);

J = min(max(J, 0), 255);

J = uint8(J);

%gain = (img_var - estimated_noise_var)./img_var;



function img_mean = local_image_mean(P)
% img_mean = local_image_mean(P)

local_img_sum = sum(sum(double(P), 4), 3);
sum_size = size(P,3)*size(P,4);

img_mean = local_img_sum/sum_size;


function img_var = local_image_var(P, img_mean, Win)
% img_var = local_img_var(P, img_mean)

diff_img_patch = double(P) - repmat(img_mean, [1 1 Win(1) Win(2)]);
diff_img_patch_pow2 = diff_img_patch.^2;
local_diff_img_sum = sum(sum(double(diff_img_patch_pow2), 4), 3);
img_var = local_diff_img_sum/(Win(1)*Win(2) - 1);

实验如下,富含噪声的图像如下:


去噪后的图像如下:



       

2018-01-15 11:03:48 u013165921 阅读数 3693

实验要求

  (a) 编写一个给图像中添加高斯噪声的程序,程序的输入参数为噪声的均值与方差。

  (b) 编写程序实现公式(5.6-11)所示的污损滤波;

  (c) 如图5.26(b)所示,对图像5.26(a) 进行+45o 方向,T = 1 的污损滤波;

  (d) 对污损后的图像加入均值为0,方差为10 的高斯噪声;

  (e) 编写程序使用公式(5.8-6)所示的参数维纳滤波对图像进行恢复。


技术论述

1、高斯噪声

  高斯噪声是指概率密度函数服从高斯分布(即正态分布)的一类噪声。常见的高斯噪声包括起伏噪声、宇宙噪声、热噪声和散粒噪声等等。

  在空间域和频率域中,由于高斯噪声在数学上的易处理性,故实践中常用这种噪声模型。

这里写图片描述

  高斯随机变量z的PDF由上式给出,其中z表示灰度值,u表示z的均值,σ表示z的方差。当z服从高斯分布时,其值有大约70%落在范围[(u-σ),(u+σ)]内,有大约95%落在范围[(u-2σ),(u+2σ)]内。

这里写图片描述

  在MATLAB中,通常使用randn()函数来产生标准正态分布的随机数或矩阵。

2、污损滤波

  通过均匀线性运动模糊对图像进行污损滤波。假定图像只在x方向以给定的速度x0(t)=at/T 做匀速直线运动。当 t=T 时,图像位移的总距离为a,则有

这里写图片描述

  若允许y分量也随着变化,按 y0=bt/T 给出的运动,则退化函数变为

这里写图片描述

3、参数维纳滤波

  维纳滤波又称最小均方滤波,是一种综合了退化函数和噪声统计特征进行复原处理的方法。它是一种基于最小均方误差准则、对平稳过程的最优估计器。这种滤波器的输出与期望输出之间的均方误差为最小,因此,维纳滤波是一个最佳滤波系统,可用于提取被平稳噪声所污染的信号。

  该方法建立在图像和噪声都是随机变量的基础上,目标是找到未污染图像f的一个估计,使它们之间的均方误差最小。这种误差度量为:

这里写图片描述

  当未退化图像的功率谱未知或者是不能估计时,通常使用的一种方法是由下面的表达式来近似:

这里写图片描述

  其中,K是一个加到|H(u,v)|2的所有项上的特定常数。


实验结果

这里写图片描述
这里写图片描述
这里写图片描述
这里写图片描述
这里写图片描述


实验程序

% --------------------------实验要求------------------------------------------

img = imread('Fig5.26(a).jpg');                             % 原图像
img_fouling = fouling_filter(img,0.1,0.1,1);                % 污损滤波
img_gaussian_10 = gaussian_noise(img_fouling,0,10);         % 高斯噪声
img_wiener = wiener_filter(img,img_gaussian_10,0.1);        % 维纳滤波

subplot(3,1,1);imshow(img);title('原图像');
subplot(3,1,2);imshow(img_fouling);title('污损滤波: +45°,T=1');
subplot(3,1,3);imshow(img_gaussian_10);  title('高斯噪声: 均值0,方差10');
figure;
imshow(img_wiener);title('维纳滤波图像');


% ------------------------ 高斯噪声比较-----------------------------------------

img_gaussian_100 = gaussian_noise(img_fouling,0,100);       % 均值0,方差100
img_gaussian_1000 = gaussian_noise(img_fouling,0,1000);     % 均值0,方差1000

figure;
subplot(2,1,1);imshow(img_gaussian_100); title('高斯噪声: 均值0,方差100(仅作对比)');
subplot(2,1,2);imshow(img_gaussian_1000);title('高斯噪声: 均值0,方差1000(仅作对比)');


% ------------------------ 调整合适参数-----------------------------------------

img_fouling = fouling_filter(img,0.01,0.01,1);              % 污损滤波
img_gaussian_10 = gaussian_noise(img_fouling,0,10);       % 高斯噪声
img_wiener = wiener_filter(img,img_gaussian_10,0.1);        % 维纳滤波

figure;
subplot(3,1,1);imshow(img_fouling);title('污损滤波: +45°,T=1');
subplot(3,1,2);imshow(img_gaussian_10);title('高斯噪声: 均值0,方差10');
subplot(3,1,3);imshow(img_wiener);title('调整参数后的维纳滤波图像');


% ------------------------MATLAB自带函数------------------------------------

psf = fspecial('motion',25,45);
mf = imfilter(img,psf,'circular','conv');
noise = imnoise(img,'gaussian',0,10);
mfn = mf + noise;
nsr = sum(noise(:).^2)/sum(mf(:).^2);
img_out = deconvwnr(mfn,psf,nsr);

figure;
subplot(3,1,1);imshow(mf);title('Matlab自带污损滤波');
subplot(3,1,2);imshow(mfn);title('Matlab添加高斯噪声');
subplot(3,1,3);imshow(img_out);title('Matlab自带维纳滤波');
% -------------------------------------------END-----------------------------------------------
% --------------------------------------污损滤波-----------------------------------------------

function img_fouling = fouling_filter(img,a,b,T) 

[M,N] = size(img);
F = fft2(img);

for u = 1:M
    for v = 1:N
        K(u,v) = pi * (u * a + v * b);
        H(u,v) = T * sin(K(u,v)) * exp(-1j * K(u,v))/K(u,v);
        G(u,v) = H(u,v) * F(u,v);
    end
end

img_fouling = ifft2(G);
img_fouling = uint8(abs(img_fouling));

end
% ----------------------------------给图像添加高斯噪声-----------------------------------

function img_noise = gaussian_noise(img,mean,var) 

[M,N] = size(img);
add_noise = mean + randn(M,N)*sqrt(var);
img_noise = img + uint8(add_noise);

end
% --------------------------------------维纳滤波-----------------------------------------------

function img_wiener = wiener_filter(img_src,img_degradation,K) 

[M,N] = size(img_src);
S = fft2(img_src);
G = fft2(img_degradation);

for u = 1:M
    for v = 1:N
        H(u,v) = G(u,v)/S(u,v);
        F(u,v) = 1/H(u,v)*(abs(H(u,v)))^2/((abs(H(u,v)))^2+K)*G(u,v);
    end
end

img_wiener = ifft2((F));
img_wiener = uint8(abs(img_wiener));

end
2014-07-26 15:41:00 EbowTang 阅读数 7727

以下为来自Matlab的基本原理

wiener2

2-D adaptive noise-removal filtering

The syntax wiener2(I,[m n],[mblock nblock],noise) has been removed. Use the wiener2(I,[m n],noise) syntax instead.

Syntax

J = wiener2(I,[m n],noise)
[J,noise] = wiener2(I,[m n])

Description

wiener2 lowpass-filters a grayscale image that has been degraded by constant power additive noise. wiener2 uses a pixelwise adaptive Wiener method based on statistics estimated from a local neighborhood of each pixel.

J = wiener2(I,[m n],noise) filters the image I using pixelwise adaptive Wiener filtering, using neighborhoods of size m-by-n to estimate the local image mean and standard deviation. If you omit the [m n] argument, m and n default to 3. The additive noise (Gaussian white noise) power is assumed to be noise.

[J,noise] = wiener2(I,[m n]) also estimates the additive noise power before doing the filtering. wiener2 returns this estimate in noise.

Class Support

The input image I is a two-dimensional image of class uint8uint16int16single, or double. The output image J is of the same size and class as I.

Examples

For an example, see Remove Noise By Adaptive Filtering.

More About

collapse all

Algorithms

wiener2 estimates the local mean and variance around each pixel.

and

where  is the N-by-M local neighborhood of each pixel in the image Awiener2 then creates a pixelwise Wiener filter using these estimates,

where ν2 is the noise variance. If the noise variance is not given, wiener2 uses the average of all the local estimated variances.

References

[1] Lim, Jae S., Two-Dimensional Signal and Image Processing, Englewood Cliffs, NJ, Prentice Hall, 1990, p. 548, equations 9.26, 9.27, and 9.29.

See Also

filter2 | medfilt2


代码实现

1,算法流程:


2,Matlab代码实现:

简化版本:
function [f,noise] = adpwiener2(g,nhood)
noise = [];
g = im2double(g);

% 计算均值
localMean = filter2(ones(nhood), g) / prod(nhood);

% 计算方差
localVar = filter2(ones(nhood), g.^2) / prod(nhood) - localMean.^2;

% 如果需要,计算噪声能量
if (isempty(noise))
  noise = mean2(localVar);
end

% 计算结果
% f = localMean + (max(0, localVar - noise) ./ ...
%           max(localVar, noise)) .* (g - localMean); 
f = g - localMean;
g = localVar - noise; 
g = max(g, 0);
localVar = max(localVar, noise);
f = f ./ localVar;
f = f .* g;
f = f + localMean;


3,C++代码实现

处理结果与matlab库函数一模一样。
//对应matlab中的filter2,但是此处功能多了一些
double filter2(
	double* srcbuffer,
	int width,
	int height,
	int kernelsize,
	double* dstbuffer1,
	double* dstbuffer2
	)
{
	//kernel的中心位置
	int center_pos = kernelsize / 2;
	int i, j;
	double sumvar = 0.0;
	for (i = 0; i < height; ++i)    // 行
	{
		for (j = 0; j < width; ++j)    // 列
		{
			for (int m = 0; m < kernelsize; ++m)     // kernel行
			{
				int mm = kernelsize - 1 - m;      // 核的行索引

				for (int n = 0; n < kernelsize; ++n) // kernel列
				{
					int nn = kernelsize - 1 - n;  // 核的列索引
					// 输入图像信号的索引,用于检查边界
					int ii = i + (m - center_pos);
					int jj = j + (n - center_pos);
					// 忽视越界
					if (ii >= 0 && ii < height && jj >= 0 && jj < width)
					{
						dstbuffer1[i*width + j] += srcbuffer[ii*width + jj] / (kernelsize*kernelsize);
						dstbuffer2[i*width + j] += srcbuffer[ii*width + jj] * srcbuffer[ii*width + jj] / (kernelsize*kernelsize);
					}
				}
			}
			dstbuffer2[i*width + j] -= dstbuffer1[i*width + j] * dstbuffer1[i*width + j];
			sumvar += dstbuffer2[i*width + j];
		}
	}
	return sumvar / (width*height);//求取noise
}


void adpwiener2(
	unsigned char* inbuffer,
	int kernelsize,
	int width,
	int height,
	unsigned char* outbuffer
	)
{
	double len = 1.0*kernelsize*kernelsize;
	double* srcbuffer = new double[width*height];//对应matlab中的g
	double* localMeanArr = new double[width*height];//对应matlab中的localMean,以下同理
	double* localVarArr = new double[width*height];
	double* dstbuffer = new double[width*height];
	//初始化
	for (int i = 0; i < width*height; i++)
		localVarArr[i] = localMeanArr[i] = 0.0;
	for (int i = 0; i < width*height; i++)
		srcbuffer[i] = double(inbuffer[i]);

	double noise = filter2(srcbuffer, width, height, kernelsize, localMeanArr, localVarArr);

	for (int i = 0; i < width*height; i++)
	{
		dstbuffer[i] = srcbuffer[i] - localMeanArr[i];//f = g - localMean;
		srcbuffer[i] = localVarArr[i] - noise;//g = localVar - noise; 
		if (srcbuffer[i] < 0.0)//g = max(g, 0);
			srcbuffer[i] = 0.0;
		localVarArr[i] = max(localVarArr[i], noise);//localVar = max(localVar, noise);
		dstbuffer[i] = dstbuffer[i] / localVarArr[i] * srcbuffer[i] + localMeanArr[i];//f = f ./ localVar.*+localMean;
		outbuffer[i] = unsigned char(dstbuffer[i]);//转换成图像数据
	}

	delete[] srcbuffer;
	srcbuffer = NULL;
	delete[] localMeanArr;
	localMeanArr = NULL;
	delete[] localVarArr;
	localVarArr = NULL;
	delete[] dstbuffer;
	dstbuffer = NULL;
}




4,实验效果:

对原图处理结果,注意是5*5的窗口




对椒盐噪声处理结果,注意是5*5的窗口,结果较差


对高斯噪声的处理结果,注意是5*5的窗口,结果较优


参考文献:

【1】Lim, Jae S. Two-Dimensional Signal and Image Processing. Englewood Cliffs, NJ: Prentice Hall, 1990. pp. 536-540.


2019-01-14 16:38:19 baidu_37973494 阅读数 156

实验要求:

(a) 编写一个给图像中添加高斯噪声的程序,程序的输入参数为噪声的均值与方差。
(b) 编写程序实现公式(5.6-11)所示的污损滤波;
(c) 如图 5.26(b)所示,对图像 5.26(a) 进行+45o 方向,T = 1 的污损滤波;
(d) 对污损后的图像加入均值为0,方差为 10 的高斯噪声;
(d) 编写程序使用公式(5.8-6)

技术论述:

 

实验代码:

% ----------------------Parametric Wiener Filter---------------------------%
function main_test
clc;clear all;close all;
 
img = imread('Fig5.26(a).jpg');                             % 原图像
img_fouling = fouling_filter(img,0.1,0.1,1);                % 污损滤波
img_fouling1 = fouling_filter(img,0.01,0.01,1);             % 污损滤波
img_gaussian_10 = gaussian_noise(img_fouling,0,10);         % 添加高斯噪声
img_gaussian_10_1 = gaussian_noise(img_fouling1,0,10);      % 添加高斯噪声
 
figure;
subplot(2,3,1);imshow(img);title('图1(a):原图像');
subplot(2,3,2);imshow(img_fouling);title('图1(b):污损滤波: +45°,T=1,a=b=0.1');
subplot(2,3,3);imshow(img_fouling1);title('图1(c):污损滤波: +45°,T=1,a=b=0.01');
subplot(2,3,4);imshow(img_gaussian_10);title('图1(d):对图(b)加入均值0,方差10的高斯噪声');
subplot(2,3,5);imshow(img_gaussian_10_1);title('图1(e):对图(c)加入均值0,方差10的高斯噪声');
 
%% -------------------------维纳滤波函数-------------------------%
img_wiener = wiener_filter(img,img_gaussian_10,0.1);        % 维纳滤波
img_wiener1 = wiener_filter(img,img_gaussian_10_1,0.1);        % 维纳滤波
figure;
subplot(1,2,1);imshow(img_wiener);title('图2(a):对图1(d)进行维纳滤波图像');
subplot(1,2,2);imshow(img_wiener1);title('图2(b):对图1(e)进行维纳滤波图像');
 
%% ------------------------MATLAB自带函数-----------------------%
psf = fspecial('motion',25,45);
mf = imfilter(img,psf,'circular','conv');
noise = imnoise(img,'gaussian',0,10);
mfn = mf + noise;
nsr = sum(noise(:).^2)/sum(mf(:).^2);
img_out = deconvwnr(mfn,psf,nsr);
 
figure;
subplot(1,3,1);imshow(mf);title('图3(a):Matlab自带污损滤波');
subplot(1,3,2);imshow(mfn);title('图3(b):Matlab添加高斯噪声');
subplot(1,3,3);imshow(img_out);title('图3(c):Matlab自带维纳滤波');
% ---------------------------------END---------------------------------%

(a) 编写一个给图像中添加高斯噪声的程序,程序的输入参数为噪声的均值与方差。 
% --------------------------给图像添加高斯噪声---------------------------%
function img_noise = gaussian_noise(img,mean,var) 
% 程序的输入为噪声的均值mean和方差var,输出为噪声图像
 
[M,N] = size(img);
add_noise = mean + randn(M,N)*sqrt(var);
img_noise = img + uint8(add_noise);
end

(b) 编写程序实现公式(5.6-11)所示的污损滤波; 
% ------------------------污损滤波函数程序------------------------%
function img_fouling = fouling_filter(img,a,b,T) 
% 输入参数:img为输入图像,a,b,T为别为技术论述中的公式(9)的三个参数
 
[M,N] = size(img);
F = fft2(img);
 
for u = 1:M
    for v = 1:N
        K(u,v) = pi * (u * a + v * b);
        H(u,v) = T * sin(K(u,v)) * exp(-1j * K(u,v))/K(u,v);
        G(u,v) = H(u,v) * F(u,v);
    end
end
 
img_fouling = ifft2(G);
img_fouling = uint8(abs(img_fouling));
end

(e) 编写程序使用公式(5.8-6)所示的参数维纳滤波对图像进行恢复。 
% ----------------------------维纳滤波程序-----------------------------%
function img_wiener = wiener_filter(img_src,img_degradation,K) 
 
[M,N] = size(img_src);
S = fft2(img_src);
G = fft2(img_degradation);
 
for u = 1:M
    for v = 1:N
        H(u,v) = G(u,v)/S(u,v);
        F(u,v) = 1/H(u,v)*(abs(H(u,v)))^2/((abs(H(u,v)))^2+K)*G(u,v);
    end
end
 
img_wiener = ifft2((F));
img_wiener = uint8(abs(img_wiener));
end

 

2019-01-05 23:47:53 LVJINYANJY 阅读数 116
% 将一幅256×256的灰度图像用3×3平均滤波器进行模糊,分别再加上一定的高斯噪声和均匀噪声。
% 然后,设计一个维纳滤波器对这两幅图像进行复原,分别计算这两幅图像复原前后的PSNR。

% 维纳滤波函数wiener2(·)
% 均匀噪声的生成方法
% PSNR计算

clc;
clear;
I=imread('2-lena.tif'); 
[m,n]=size(I); 
%均值滤波 3*3 
J=imfilter(I,fspecial('average',3));   
PSF=fspecial('average',[3 3]);
%加高斯噪声和均匀噪声 
K1=imnoise(J,'gaussian',0,0.01); 
%产生随机的均匀噪声 
K2=imnoise(J,'speckle',0.04);  
%维纳滤波器复原 
noise=K1-J; 
NSR=sum(noise(:).^2)/sum(K1(:).^2); 
M1=deconvwnr(K1,PSF,NSR); 
noise=K2-J; 
NSR=sum(noise(:).^2)/sum(K2(:).^2); 
M2=deconvwnr(K2,PSF,NSR);

psnr1=PSNR(I,M1);
psnr2=PSNR(I,M2);

figure; 
subplot(321),imshow(I);
xname = sprintf('原图像');
xlabel(xname);
subplot(322),imshow(J);
xname = sprintf('均值滤波后图像');
xlabel(xname);
subplot(323),imshow(K1);
xname = sprintf('加高斯噪声后图像');
xlabel(xname);
subplot(324),imshow(K2); 
xname = sprintf('加均匀噪声后图像');
xlabel(xname);    
subplot(325),imshow(M1);
xname = sprintf('高斯噪声图像的维纳滤波复原 PSNR = %.2f dB\n',psnr1);
xlabel(xname);
subplot(326),imshow(M2);
xname = sprintf('均匀噪声图像的维纳滤波复原 PSNR = %.2f dB\n',psnr2);
xlabel(xname);

function psnr = PSNR(x,y)
x = double(x);
y = double(y);
cnt = length(x(:));
mse = sum((x(:)-y(:)).^2)/cnt;
psnr = 10*log10(255^2/mse);
end

 

没有更多推荐了,返回首页