1. 图像复原与重建
1.1. 噪声模型
常见的噪声包括胡椒噪声、盐粒噪声、高斯噪声、均匀噪声、瑞利噪声、指数噪声、伽马噪声
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
f = imread('../pic/1.jpg');
figure;
subplot(4,2,1);
imshow(f)
title('原始图像');
[M,N] = size(f);
R = imnoise2('salt & pepper',M,N,0.1,0);
c = R == 0;
gp = f;
gp(c) = 0;
subplot(4,2,2);
imshow(gp)
title('胡椒噪声图像(0.1)')
R = imnoise2('salt & pepper',M,N,0,0.1);
c = find(R == 1);
gs = f;
gs(c) = 255;
subplot(4,2,3);
imshow(gs)
title('盐粒噪声图像(0.1)')
gg = imnoise(f,'gaussian',0,1);
subplot(4,2,4);
imshow(gg);
title('高斯噪声图像(0,1)');
R = imnoise2('uniform',M,N,0,1);
R = uint8(R);
c = R ==0;
gu = f;
gu(c) = 0;
subplot(4,2,5);
imshow(gu);
title('均匀噪声图像');
R= imnoise2('rayleigh',M,N,0.1,1);
R =uint8(R);
c = R ==0;
gr = f;
gr(c)=0;
subplot(4,2,6);
imshow(gr);
title('瑞利噪声图像');
R= imnoise2('exponential',M,N,1,1);
R = uint8(R);
c = R ==0;
ge =f;
ge(c)=0;
subplot(4,2,7);
imshow(ge);
title('指数噪声图像');
R = imnoise2('erlang',M,N,2,5);
R = uint8(R);
c = R ==0;
ger =f;
ger(c)=0;
subplot(4,2,8);
imshow(ger);
title('伽马噪声图像');
1.2. 图像复原
1.2.1. 空间域滤波
function f = spfilt(g, type, m, n, parameter)
%SPFILT Performs linear and nonlinear spatial filtering.
% F = SPFILT(G, TYPE, M, N, PARAMETER) performs spatial filtering
% of image G using a TYPE filter of size M-by-N. Valid calls to
% SPFILT are as follows:
%
% F = SPFILT(G, 'amean', M, N) Arithmetic mean filtering.
% F = SPFILT(G, 'gmean', M, N) Geometric mean filtering.
% F = SPFILT(G, 'hmean', M, N) Harmonic mean filtering.
% F = SPFILT(G, 'chmean', M, N, Q) Contraharmonic mean
% filtering of order Q. The
% default is Q = 1.5.
% F = SPFILT(G, 'median', M, N) Median filtering.
% F = SPFILT(G, 'max', M, N) Max filtering.
% F = SPFILT(G, 'min', M, N) Min filtering.
% F = SPFILT(G, 'midpoint', M, N) Midpoint filtering.
% F = SPFILT(G, 'atrimmed', M, N, D) Alpha-trimmed mean filtering.
% Parameter D must be a
% nonnegative even integer;
% its default value is D = 2.
%
% The default values when only G and TYPE are input are M = N = 3,
% Q = 1.5, and D = 2.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.6 $ $Date: 2003/10/27 20:07:00 $
% Process inputs.
if nargin == 2
m = 3; n = 3; Q = 1.5; d = 2;
elseif nargin == 5
Q = parameter; d = parameter;
elseif nargin == 4
Q = 1.5; d = 2;
else
error('Wrong number of inputs.');
end
% Do the filtering.
switch type
case 'amean'
w = fspecial('average', [m n]);
f = imfilter(g, w, 'replicate');
case 'gmean'
f = gmean(g, m, n);
case 'hmean'
f = harmean(g, m, n);
case 'chmean'
f = charmean(g, m, n, Q);
case 'median'
if length(size(g))==3
f1 = medfilt2(g(:,:,1), [m n], 'symmetric');
f2 = medfilt2(g(:,:,2), [m n], 'symmetric');
f3 = medfilt2(g(:,:,3), [m n], 'symmetric');
f=cat(3,f1,f2,f3);
else
f = medfilt2(g, [m n], 'symmetric');
end
case 'max'
if length(size(g))==3
f1 = ordfilt2(g(:,:,1), m*n, ones(m, n), 'symmetric');
f2 = ordfilt2(g(:,:,2), m*n, ones(m, n), 'symmetric');
f3 = ordfilt2(g(:,:,3), m*n, ones(m, n), 'symmetric');
f=cat(3,f1,f2,f3);
else
f = ordfilt2(g, m*n, ones(m, n), 'symmetric');
end
case 'min'
if length(size(g))==3
f1 = ordfilt2(g(:,:,1), 1, ones(m, n), 'symmetric');
f2 = ordfilt2(g(:,:,2), 1, ones(m, n), 'symmetric');
f3 = ordfilt2(g(:,:,3), 1, ones(m, n), 'symmetric');
f=cat(3,f1,f2,f3);
else
f = ordfilt2(g,1, ones(m, n), 'symmetric');
end
case 'midpoint'
if length(size(g))==3
f1 = ordfilt2(g(:,:,1), 1, ones(m, n), 'symmetric');
f2 = ordfilt2(g(:,:,2), 1, ones(m, n), 'symmetric');
f3 = ordfilt2(g(:,:,3), 1, ones(m, n), 'symmetric');
g1=cat(3,f1,f2,f3);
else
g1 = ordfilt2(g,1, ones(m, n), 'symmetric');
end
if length(size(g))==3
f1 = ordfilt2(g(:,:,1), m*n, ones(m, n), 'symmetric');
f2 = ordfilt2(g(:,:,2), m*n, ones(m, n), 'symmetric');
f3 = ordfilt2(g(:,:,3), m*n, ones(m, n), 'symmetric');
g2=cat(3,f1,f2,f3);
else
g2 = ordfilt2(g, m*n, ones(m, n), 'symmetric');
end
f = imlincomb(0.5, g1, 0.5, g2);
case 'atrimmed'
if (d <= 0) | (d/2 ~= round(d/2))
error('d must be a positive, even integer.')
end
f = alphatrim(g, m, n, d);
otherwise
error('Unknown filter type.')
end
%-------------------------------------------------------------------%
function f = gmean(g, m, n)
% Implements a geometric mean filter.
inclass = class(g);
g = im2double(g);
% Disable log(0) warning.
warning off;
f = exp(imfilter(log(g), ones(m, n), 'replicate')).^(1 / m / n);
warning on;
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function f = harmean(g, m, n)
% Implements a harmonic mean filter.
inclass = class(g);
g = im2double(g);
f = m * n ./ imfilter(1./(g + eps),ones(m, n), 'replicate');
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function f = charmean(g, m, n, q)
% Implements a contraharmonic mean filter.
inclass = class(g);
g = im2double(g);
f = imfilter(g.^(q+1), ones(m, n), 'replicate');
f = f ./ (imfilter(g.^q, ones(m, n), 'replicate') + eps);
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function f = alphatrim(g, m, n, d)
% Implements an alpha-trimmed mean filter.
inclass = class(g);
g = im2double(g);
f = imfilter(g, ones(m, n), 'symmetric');
for k = 1:d/2
f = imsubtract(f, ordfilt2(g, k, ones(m, n), 'symmetric'));
end
for k = (m*n - (d/2) + 1):m*n
f = imsubtract(f, ordfilt2(g, k, ones(m, n), 'symmetric'));
end
f = f / (m*n - d);
f = changeclass(inclass, f);
1.2.2. 均值滤波器
算术均值滤波器
几何均值滤波器
谐波均值滤波器
逆谐波滤均值波器
Q称为滤波器的阶数。当Q为正数时，用于消除“胡椒”噪声；
当Q为负数时，用于消除“盐”噪声，但不能同时消除“椒盐”噪声
当Q = 0，逆谐波均值滤波器转变为算术均值滤波器
当Q = -1，逆谐波均值滤波器转变为谐波均值滤波器
f = imread('../pic/5.jpg');
[M,N] = size(f);
R =imnoise2('gaussian',M,N,1,20);
R = uint8(R);
idx = R==1;
gg=f;
gg(idx)=255;
figure;
subplot(231);
imshow(gg);
title('高斯噪声图像');
R1 = imnoise2('salt & pepper',M,N,0.1,0);
idx = R1==0;
gp= f;
gp(idx)=0;
R2 = imnoise2('salt & pepper',M,N,0,0.1);
idx = find(R2==1);
gs=f;
gs(idx)=255;
subplot(232);
imshow(gs);
title('盐粒噪声图像');
subplot(233);
imshow(gp);
title('胡椒噪声图像');
f1= spfilt(gs,'hmean',3,3);
subplot(234);
imshow(f1);
title('谐波均值滤波器处理后的盐粒噪声图像');
f2= spfilt(gp,'hmean',3,3);
subplot(235);
imshow(f2);
title('谐波均值滤波器处理后的胡椒噪声图像');
f3= spfilt(gg,'hmean',3,3);
subplot(236);
imshow(f3);
title('谐波均值滤波器处理后的高斯噪声图像');
1.2.3. 顺序统计滤波器
最大、最小滤波器
figure; subplot(2,2,1); imshow(gs); title('胡椒粒噪声图像(0.1)'); fpmax = spfilt(gp,'max'); subplot(2,2,2); imshow(fpmax); title('用3*3最大滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果'); subplot(2,2,3); imshow(gs); title('胡椒粒噪声图像(0.1)'); fsmin = spfilt(gp,'min',3,3); subplot(224); imshow(fsmin); title('用3*3最小滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果');
中值滤波器
fmean = spfilt(gp,'median',3,3); figure; subplot(221); imshow(gp);title('胡椒噪声图像(0.1)'); subplot(222); imshow(fmean);title('中值滤波后的图像'); favg = spfilt(gp,'amean',3,3); subplot(223); imshow(favg);title('算术均值滤波后的图像'); subplot(224); fgm = spfilt(gp,'gmean',3,3); imshow(fgm);title('几何均值滤波后的图像');
中点滤波器
结合了顺序统计和求平均
对于高斯和均匀随机分布这类噪声有最好的效果
figure;
subplot(221);
imshow(gg);
title('高斯噪声');
subplot(222);
imshow(gu);
title('均匀噪声');
f1 = spfilt(gg,'midpoint',3,3);
subplot(223);
imshow(f1);title('中点滤波后的高斯噪声污染图像');
subplot(224);
f2 = spfilt(gu,'midpoint',3,3);
imshow(f2);title('中点滤波后的均匀噪声污染图像');
阿尔法滤波器
- 在邻域内去掉g(s,t)最高灰度值的d/2和最低灰度值的d/2对于高斯和均匀随机分布这类噪声有最好的效果
- 代表剩余的个像素
- 当d=0，退变为算术均值滤波器
- 当d=(mn-1)/2，退变为中值滤波器
- 当d取其它值时，适用于包括多种噪声的情况下，例如高斯噪声和椒盐噪声混合的情况
f = imread('../pic/4.jpg');
figure;
subplot(221);
imshow(f);
title('原始图像');
f2 = imnoise(f,'salt & pepper',0.1);
subplot(222);
imshow(f2);
title('增加椒盐噪声后的图像（0.1）');
[M,N] =size(f);
R = imnoise2('gaussian',M,N,0,100);
R = uint8(R);
c = R==1;
gg = f;
gg(c) = 255;
subplot(223)
imshow(gg);
title('增加高斯噪声后的图像（0,100）');
sp = spfilt(f,'gmean',8,8);
subplot(224);
imshow(sp);
title('用8*8几何平均滤波器对[被高斯噪声污染的图像]处理后的结果');
imshow(f);
title('原图像');
f = imnoise(f,'salt & pepper',0.1);
imshow(f);
title('增加椒盐噪声后的图像');
sp =spfilt(f2,'median',3,3); %中值滤波器
figure;
imshow(sp);
title('用3*3的均值滤波器对[被概率为0.1的椒盐噪声污染的图像]滤波后的结果');
fmax = spfilt(f,'max',3,3);
fmin=spfilt(f,'min',3,3);
fmid = spfilt(f,'midpoint',3,3);
figure;
imshow(fmax);title('用3*3的最大滤波器对[被概率为0.1的椒盐噪声污染的图像]滤波后的结果');
figure;
imshow(fmax);title('用3*3的最小滤波器对[被概率为0.1的椒盐噪声污染的图像]滤波后的结果');
figure;
imshow(fmax);title('用3*3的中点滤波器对[被概率为0.1的椒盐噪声污染的图像]滤波后的结果');
sp =spfilt(f(:,:,1),'atrimmed',3,3);
sp2 =spfilt(f(:,:,2),'atrimmed',3,3);
sp3 =spfilt(f(:,:,3),'atrimmed',3,3);
sp = cat(3,sp,sp2,sp3);
figure;
imshow(sp);
title('用3*3的阿尔法滤波器对[被概率为0.1的椒盐噪声污染的图像]滤波后的结果');
1.2.4. 自适应滤波器
function f = adpmedian(g, Smax)
%ADPMEDIAN Perform adaptive median filtering.
% F = ADPMEDIAN(G, SMAX) performs adaptive median filtering of
% image G. The median filter starts at size 3-by-3 and iterates up
% to size SMAX-by-SMAX. SMAX must be an odd integer greater than 1.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.5 $ $Date: 2003/11/21 14:19:05 $
% SMAX must be an odd, positive integer greater than 1.
if (Smax <= 1) | (Smax/2 == round(Smax/2)) | (Smax ~= round(Smax))
error('SMAX must be an odd integer > 1.')
end
[M, N] = size(g);
% Initial setup.
f = g;
f(:) = 0;
alreadyProcessed = false(size(g));
% Begin filtering.
for k = 3:2:Smax
zmin = ordfilt2(g, 1, ones(k, k), 'symmetric');
zmax = ordfilt2(g, k * k, ones(k, k), 'symmetric');
zmed = medfilt2(g, [k k], 'symmetric');
processUsingLevelB = (zmed > zmin) & (zmax > zmed) & ...
~alreadyProcessed;
zB = (g > zmin) & (zmax > g);
outputZxy = processUsingLevelB & zB;
outputZmed = processUsingLevelB & ~zB;
f(outputZxy) = g(outputZxy);
f(outputZmed) = zmed(outputZmed);
alreadyProcessed = alreadyProcessed | processUsingLevelB;
if all(alreadyProcessed(:))
break;
end
end
% Output zmed for any remaining unprocessed pixels. Note that this
% zmed was computed using a window of size Smax-by-Smax, which is
% the final value of k in the loop.
f(~alreadyProcessed) = zmed(~alreadyProcessed);
figure;
subplot(221);
imshow(gg);
title('高斯噪声污染的图像(0.01)');
famean=spfilt(gg,'amean',8,8);
fgmean=spfilt(gg,'gmean',8,8);
subplot(222);
imshow(famean);title('8*8的算术均值滤波器处理后的图像');
subplot(223);
imshow(fgmean);title('8*8的几何均值滤波器处理后的图像');
am = adpmedian(gg(:,:,1),3);
am2 = adpmedian(gg(:,:,2),3);
am3 = adpmedian(gg(:,:,3),3);
am = cat(3,am,am2,am3);
subplot(224);
imshow(am);
title('用7*7的自适应滤波器对[被概率为0.1的高斯噪声污染的图像]滤波后的结果');
