Reference:《数字图像处理MATLAB版》冈萨雷斯
1.噪声模型
常见的噪声模型(z是随机变量,灰度值;
μ
表示均值,
s2
表示方差):
1、高斯噪声
2、瑞里噪声 其中,
μ=a+(π∗b4)−−−−√,s2=b(4−π)4
3、伽马噪声 其中,
μ=ba
,
σ2=ba2
4、指数分布噪声 其中,
μ=1a,σ2=1a2
5、均匀分布噪声 其中,
μ=(a+b)2,σ2=(b−a)212
6、椒盐噪声其中盐表示亮点,椒表示暗点。
7、周期噪声,比如空间域图像受到正弦波信号干扰。一幅图像的周期噪声是在图像获取期间由电力或机电干扰产生的
纯正弦波的二维DFT:
sin(2πμ0x+2πν0y)⇌12j[δ(μ+Mμ0,ν+Nν0)−δ(μ+Mμ0,ν+Nν0)]
而,
μ∈[0,M],ν∈[0,N]
2. Matalb的噪声生产函数imnoise
g = imnoise(f,type,parameters);
%Matlab自带函数
g=imnoise(f,'gaussian',m, var);%默认均值0,方差0.01
g=imnoise(f,'gaussian', var);%均值0,局部方差V.V与f大小相同,每一个点是包含了期望的方差
g=imnoise(f,'localvar',image_intensity,var);%均值为0的高斯噪声添加到f上,其中噪声的局部方差var是f的灰度值的函数。image_intensity和var是相同大小的向量。plot(image_intensity,var)绘制出噪声方差和图像灰度间的函数关系。向量image_intensity是[0,1]内的归一化灰度值
g=imnoise(f,'salt&pepper',d);%d是噪声密度,大约都有d*numel(f)个像素收到了影响
g=imnoise(f,'speckle',var);%使用方程g=f+n.*f将乘性噪声添加到f上,其中n是均值为0,方差为var的均匀分布的随机噪声。
g=imnoise(f,'poisson');%从数据中生成泊松噪声,而不是将人为噪声添加到数据中。
3. 使用规定分布生成空间随机噪声
方法:F(z) = w;反解出 z=F−1(w)
具体实现见imnoise2。
4. 估计噪声参数
(1) 画出ROI区域;
(2) 计算这个区域内的直方图
(3) 画出直方图,计算均值和方差(用中心矩的计算函数)
(4) 用计算出来的均值和方差,带入到imnoise2函数中,比如高斯,瑞利分布等,看哪一种情况最为匹配。
其中,imnoise2的M要与ROI的像素总数一样,这样才有可比较性。
其中要用到的函数如下:
[B,c,r] = roipoly(f);%在f的图片中选取一个感兴趣区域。交互式地选择感兴趣区域
[h,npix] = histroi(f,c,r);%计算直方图和总像素个数
figure,bar(h,1)
[v,unv] = statmoments(h,2);%计算2阶中心矩,其中h为直方图的一列,n为计算n阶中心矩
X = imnoise2('gaussian', npix,1,147,20);%上面计算出均值为147,方差选为400(标准差约为200)
figure, hist(X,130)
axis([0 300 0 140])
5. 仅有噪声的复原
工具箱使用imfilter来实现线性空间滤波,
g = imfilter(f,w,filtering_mode, boundary_options, size_options);
filtering_mode 滤波模式 ‘corr’相关;’conv’卷积; 默认是执行相关操作
boundary_options 边界选项:p 边界填充p;默认选项,P=0;
‘replicate’ 复制边界像素来填充;
‘symmetric’ 图像的大小通过边界镜像反射来扩展;’circular’ 将图像处理为二维周期函数的一个周期来扩展size_options 大小选项:
‘full’输出图像与扩展后的图像大小相同;
‘same’,输出图像与输入图像的大小相同;默认选项为same.
(1)imfilter使用时的注意事项
(1)空域滤波函数imfilter的使用
卷积与相关
f corr w != w corr f;卷积的顺序对结果有影响,两者正好相反,如g与g1;
当滤波器是对称时,imfilter采用corr 还是 conv,结果都是一样的。因为滤波器翻转180后,还是等于它自身。full 与 same的区别
full是完全卷积,g = f+w-1 (元素个数);same 是w的中心点与f的第一个对齐开始计算g(1),到w的中心点与f的最后一个点对齐为止。即g的左右两侧分别剪去了 w−12 个元素。f = [0 0 0 1 0 0 0 0]; w = [1 2 3 2 0]; g = imfilter(f, w, 'corr','full') g = 0 0 0 0 2 3 2 1 0 0 0 0 g1 = imfilter(w, f, 'corr', 'full') g1 = 0 0 0 0 1 2 3 2 0 0 0 0
比较各个参数的结果
这个地方有一个重要问题,开始运行出来效果都是一样的,没有区别。原因是:f是uint8,不是doulbe类型,转换成double类型,就会有结果。引起意外的原因:模板的系数不在[0,1]内求和,引起滤波后的结果超出[0,255]. imfilter将输出转换成与输入相同的类型。
解决办法:(1)归一化滤波器系数,使得系数的和在[0,1]内;(2)f转换为double或者single输入
f = imread('Fig0315(a)(original_test_pattern).tif'); f = double(f); w = ones(31); subplot(2,3,1);imshow(f,[]);title('原始图像') gd = imfilter(f, w);%默认相关,默认边界扩展0,默认'same' subplot(2,3,2);imshow(gd,[]);title('默认参数的结果 p=0 same') gr = imfilter(f, w, 'replicate'); subplot(2,3,3);imshow(gr,[]);title('复制边界像素') gc = imfilter(f, w, 'circular'); subplot(2,3,4);imshow(gc, []);title('边界当成周期性函数'); gs = imfilter(f, w, 'symmetric'); subplot(2, 3, 5);imshow(gs,[]);title('边界为镜像像素') f8 = im2uint8(f); gr8 = imfilter(f8, w, 'replicate'); subplot(2,3,6);imshow(gr8,[]);title('将原始图像转换为uint8 ,replicate')
附录1 如何将一个矩阵归一化到 [0,1]
很简单,用函数mapminmax,文档太长我就不翻译了,只提醒几个关键
1: 默认的map范围是[-1, 1],所以如果需要[0, 1],则按这样的格式提供参数:
MappedData = mapminmax(OriginalData, 0, 1);
2 只按行归一化,如果是矩阵,则每行各自归一化,如果需要对整个矩阵归一化,用如下方法:
FlattenedData = OriginalData(:)’; % 展开矩阵为一列,然后转置为一行。
MappedFlattened = mapminmax(FlattenedData, 0, 1); % 归一化。
MappedData = reshape(MappedFlattened, size(OriginalData)); % 还原为原始矩阵形式。此处不需转置回去,因为reshape恰好是按列重新排序法二:X-X平均/(Xmax-Xmin);但是这样的结果在[-1,1]之间
X -XMin/(Xmax-Xmin);这个就在[0,1]之间了。
附录2 周期噪声——正弦噪声 imnoise3也有详细介绍
%sin(2pi*u0*x + 2pi*v0*y)
%%备注效果:F=fft2(f)还有效果,但是一圈有很多个点,跟以为的只有两个点不太符合,也不知正确与否
%%F=fft2(MappedData)反而更糟糕。
u0 = 30; v0 =30;
% f = imread('lena.bmp');
% f = f(:,:,1);
for i =1:256
for j = 1:256
f(i,j) = sin(2*pi*u0*i+2*pi*v0*j);
end
end
%测试归一化函数
OriginalData =f;
FlattenedData = OriginalData(:);%展开成一列,做归一化要转置一下
MappedFlattened = mapminmax(FlattenedData', 0, 1); %归一化
MappedData = reshape(MappedFlattened, size(OriginalData)); %还原为原始矩阵形式。此处不需转置回去,因为reshape恰好是按列重新排序
MappedData = MappedData*255;
F = fft2(f);
S = log(1+fftshift(abs(F)));
imshow(S,[]);
附录3 imnoise2 imnoise3具体代码——详细研究,才能明白整个过程
%%产生六种噪声,根据其分布来计算的,见上图表格,M*N矩阵大小,a,b为参数
function R = imnoise2(type, M, N, a, b)
%IMNOISE2 Generates an array of random numbers with specified PDF.
% R = IMNOISE2(TYPE, M, N, A, B) generates an array, R, of size
% M-by-N, whose elements are random numbers of the specified TYPE
% with parameters A and B. If only TYPE is included in the
% input argument list, a single random number of the specified
% TYPE and default parameters shown below is generated. If only
% TYPE, M, and N are provided, the default parameters shown below
% are used. If M = N = 1, IMNOISE2 generates a single random
% number of the specified TYPE and parameters A and B.
%
% Valid values for TYPE and parameters A and B are:
%
% 'uniform' Uniform random numbers in the interval (A, B).
% The default values are (0, 1).
% 'gaussian' Gaussian random numbers with mean A and standard
% deviation B. The default values are A = 0, B = 1.
% 'salt & pepper' Salt and pepper numbers of amplitude 0 with
% probability Pa = A, and amplitude 1 with
% probability Pb = B. The default values are Pa =
% Pb = A = B = 0.05. Note that the noise has
% values 0 (with probability Pa = A) and 1 (with
% probability Pb = B), so scaling is necessary if
% values other than 0 and 1 are required. The noise
% matrix R is assigned three values. If R(x, y) =
% 0, the noise at (x, y) is pepper (black). If
% R(x, y) = 1, the noise at (x, y) is salt
% (white). If R(x, y) = 0.5, there is no noise
% assigned to coordinates (x, y).
% 'lognormal' Lognormal numbers with offset A and shape
% parameter B. The defaults are A = 1 and B =
% 0.25.
% 'rayleigh' Rayleigh noise with parameters A and B. The
% default values are A = 0 and B = 1.
% 'exponential' Exponential random numbers with parameter A. The
% default is A = 1.
% 'erlang' Erlang (gamma) random numbers with parameters A
% and B. B must be a positive integer. The
% defaults are A = 2 and B = 5. Erlang random
% numbers are approximated as the sum of B
% exponential random numbers.
% Copyright 2002-2006 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.6 $ $Date: 2006/07/15 20:44:52 $
% Set default values.
if nargin == 1
a = 0; b = 1;
M = 1; N = 1;
elseif nargin == 3
a = 0; b = 1;
end
% Begin processing. Use lower(type) to protect against input being
% capitalized.
switch lower(type)
case 'uniform'
R = a + (b - a)*rand(M, N);
case 'gaussian'
R = a + b*randn(M, N);
case 'salt & pepper'
if nargin <= 3
a = 0.05; b = 0.05;
end
% Check to make sure that Pa + Pb is not > 1.
if (a + b) > 1
error('The sum Pa + Pb must not exceed 1.')
end
R(1:M, 1:N) = 0.5;
% Generate an M-by-N array of uniformly-distributed random numbers
% in the range (0, 1). Then, Pa*(M*N) of them will have values <=
% a. The coordinates of these points we call 0 (pepper
% noise). Similarly, Pb*(M*N) points will have values in the range
% > a & <= (a + b). These we call 1 (salt noise).
X = rand(M, N);
c = find(X <= a);
R(c) = 0;
u = a + b;
c = find(X > a & X <= u);
R(c) = 1;
case 'lognormal'
if nargin <= 3
a = 1; b = 0.25;
end
R = exp(b*randn(M, N) + a);
case 'rayleigh'
R = a + (-b*log(1 - rand(M, N))).^0.5;
case 'exponential'
if nargin <= 3
a = 1;
end
if a <= 0
error('Parameter a must be positive for exponential type.')
end
k = -1/a;
R = k*log(1 - rand(M, N));
case 'erlang'
if nargin <= 3
a = 2; b = 5;
end
if (b ~= round(b) | b <= 0)
error('Param b must be a positive integer for Erlang.')
end
k = -1/a;
R = zeros(M, N);
for j = 1:b
R = R + k*log(1 - rand(M, N));
end
otherwise
error('Unknown distribution type.')
end
%%产生周期正弦噪声 指定频率 R是其傅里叶变换,可以看下书上的公式推导
function [r, R, S] = imnoise3(M, N, C, A, B)
%IMNOISE3 Generates periodic noise.
% [r, R, S] = IMNOISE3(M, N, C, A, B), generates a spatial
% sinusoidal noise pattern, r, of size M-by-N, its Fourier
% transform, R, and spectrum, S. The remaining parameters are:
% C是K对频率值
% C is a K-by-2 matrix with K pairs of freq. domain coordinates (u,
% v) that define the locations of impulses in the freq. domain. The
% locations are with respect to the frequency rectangle center at
% (floor(M/2) + 1, floor(N/2) + 1). The impulse locations are spe-
% cified as increments with respect to the center. For ex, if M =
% N = 512, then the center is at (257, 257). To specify an impulse
% at (280, 300) we specify the pair (23, 43); i.e., 257 + 23 = 280,
% and 257 + 43 = 300. Only one pair of coordinates is required for
% each impulse. The conjugate pairs are generated automatically.
%
% A is a 1-by-K vector that contains the amplitude of each of the
% K impulse pairs. If A is not included in the argument, the
% default used is A = ONES(1, K). B is then automatically set to
% its default values (see next paragraph). The value specified
% for A(j) is associated with the coordinates in C(j, 1:2).
%
% B is a K-by-2 matrix containing the Bx and By phase components
% for each impulse pair. The default values for B are B(1:K, 1:2)
% = 0.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.5 $ $Date: 2004/11/04 22:32:42 $
% Process input parameters.
[K, n] = size(C);
if nargin == 3
A(1:K) = 1.0;
B(1:K, 1:2) = 0;
elseif nargin == 4
B(1:K, 1:2) = 0;
end
% Generate R.
R = zeros(M, N);
for j = 1:K
u1 = M/2 + 1 + C(j, 1); v1 = N/2 + 1 + C(j, 2);
R(u1, v1) = i * (A(j)/2) * exp(i*2*pi*C(j, 1) * B(j, 1)/M);
% Complex conjugate.
u2 = M/2 + 1 - C(j, 1); v2 = N/2 + 1 - C(j, 2);
R(u2, v2) = -i * (A(j)/2) * exp(i*2*pi*C(j, 2) * B(j, 2)/N);
end
% Compute spectrum and spatial sinusoidal pattern.
S = abs(R);
r = real(ifft2(ifftshift(R)));