数字图像处理之图像复原（1）
说明：图像复原是一种客观的操作，通过使用退化现象的先验知识重建或恢复一副退化的图像；图像在形成、传输和记录的过程中，由于受多种原因的影响，图像的质量会有下降，典型表现为图像模糊、失真、噪声等，这一降质的过程称为图像的退化。而图像复原试图利用退化现象的某种先验知识（即退化模型），把已经退化了的图像加以重建和复原。其目的就是尽可能地减少或去除在获取图像过程中发的图像质量的下降（退化），恢复被退化图像的本来面目。
（一）图像退化、复原处理的模型
如图所示，用退化函数对退化过程建模，它和附加噪声选项一起，作用于输入图像f(x,y)，产生一幅退化的图像g(x,y)：
如果H是线性的、空间不变的过程，那么退化图像在空间域将由下面的式子表示：
等价的频域表达式为：
（二）噪声模型类型
说明：能够模拟噪声的行为和效果的能力是图像复原的核心。下面我们将对两种基本噪声模型进行演示：空域中的噪声（用噪声概率密度函数来描述）以及频域中的噪声（用噪声的各种傅里叶特性来描述）。
- imnoise 能够产生可用类型和参数的噪声。
- imnoise2能够产生一个M*N的噪声数组R，它不以任何方式缩放，另一个主要不同是imnoise输出一个有噪声的图像，而imnoise2产生噪声模式本身。
- imnoise3产生周期噪声模型。
（1）用imnoise函数为图像添加噪声
在matlab中采用imnoise函数使图像添加上噪声。其基本语法如下：
g = imnoise(f,type,parameters)
注意： 使用函数之前需要将图像转换成范围在【0，1】的double类型！
这个函数的语法形式如下：
g = imnoise(f,'gaussian',m,var) %将均值为m、方差为var的高斯噪声加到图像f上，默认值均值为0、方差为0.01的噪声
g = imnoise(f,'localvar',V) %将均值为0、局部方差为V的高斯噪声加到图像f上
g = imnoise(f,'localvar',image_intensity,var) %将均值为0的高斯噪声加到图像f上,噪声的局部方差var是图像f的灰度值函数，自变量image_intensity和var是大小相同的向量，且向量image_intensity必须包含归一化的灰度值（范围是【0,1】）
g = imnoise(f,'salt & pepper',d) %用椒盐噪声污染图像f，d是噪声密度（默认值为0.05）
g = imnoise(f,'speckle',var) %用方程式g=f+n.*f将乘性噪声添加到图像f上，var默认值是0.04
g = imnoise(f,'poisson') %生成泊松噪声
实验代码及结果如下：（分别为高斯噪声、椒盐噪声、乘性噪声、泊松噪声）
f = imread('C:\Users\Public\Pictures\Sample Pictures\horse.jpg');
f1 = rgb2gray(f);
g = imnoise(f1,'gaussian',0,0.01);
g1 = imnoise(f1,'gaussian',0.5,0.01);
g2 = imnoise(f1,'gaussian',0,0.05);
subplot(131),imshow(g),title('默认值的高斯噪声污染的图像');
subplot(132),imshow(g1),title('均值为0.5、方差为0.01的高斯噪声污染的图像');
subplot(133),imshow(g2),title('均值为0、方差为0.05的高斯噪声污染的图像');
结论：
图一与图二对比：均值增大，图像越明亮，大部分颜色细节丢失，只有原图像中颜色较深的部分显示出来，因此这种噪声污染比较严重；
图一与图三对比：方差变大，噪声越明显，同时也会丢失细节，虽然看起来没有图二损失的严重，但是会有点状噪声产生，影响辨识。
f = imread('C:\Users\Public\Pictures\Sample Pictures\horse.jpg');
f1 = rgb2gray(f);
g = imnoise(f,'salt & pepper',0.05);
g1 = imnoise(f,'salt & pepper',0.08);
g2 = imnoise(f,'salt & pepper',1);
subplot(131),imshow(g),title('噪声密度为默认值的椒盐噪声污染的图像');
subplot(132),imshow(g1),title('噪声密度为0.08的椒盐噪声污染的图像');
subplot(133),imshow(g2),title('噪声密度为1的椒盐噪声污染的图像');
结论： 噪声密度的增大使得图像的丢失细节越多，最后完全消除图像本来的模样。噪声密度达到0.08时还可以识别图像内容，但在噪声密度达到1时，图片中的细节及内容完全被遮盖住，无法识别，这种噪声污染是很严重的。
f = imread('C:\Users\Public\Pictures\Sample Pictures\horse.jpg');
f1 = rgb2gray(f);
g = imnoise(f,'speckle',0.01);
g1 = imnoise(f,'speckle',0.04);
g2 = imnoise(f,'speckle',0.08);
subplot(131),imshow(g),title('参数为0.01乘性噪声污染的图像');
subplot(132),imshow(g1),title('参数为0.04乘性噪声污染的图像');
subplot(133),imshow(g2),title('参数为0.08乘性噪声污染的图像');
结论：
如图所示，乘性噪声不仅会模糊图像，而且会像椒盐噪声一样，产生黑白交替的点状。随着参数的增大，图像被污染得越明显。乘性噪声的参数为0.01时，噪声对图像的影响较小，参数的增大使得图像中某些细节的丢失，而且造成的污染很严重。
f = imread('C:\Users\Public\Pictures\Sample Pictures\horse.jpg');
f1 = rgb2gray(f);
g = imnoise(f,'poisson');
imshow(g),title('泊松噪声污染的图像');
结论： 泊松噪声会模糊图像，并不像椒盐噪声会产生黑白交替的点状。
（2）用给定分布产生空间随机噪声（函数imnoise2）
说明：通常，在函数imnoise中能够产生可用类型和参数的噪声是很有必要的。空间噪声值是随机数，以概率密度函数（PDF）或是等价的、相应的累积分布函数（CDF）为特征。
基本语法如下：
A = rand(M,N) %数组元素在区间[0,1]内均匀分布的数，且产生单一的随机数
A = randn(M,N) %数组元素是零均值、单位方差的正态（高斯）数，若无参数也是产生单一随机数
I = find(A); % 返回I中A的所有非零元素的线性索引，没找到则返回一个空矩阵
[r,c] = find(A); % 返回矩阵A的非零元素的行和列索引
[r,c,v] = find(A); % 除了返回行索引和列索引外，还以列向量v返回A的非零值
I = find(A<128);
A(I) = 0; % 寻找图像中值小于128的像素并把它们设置为0
I = find(A >=64 & A <= 192);
A(I) = 128; % 将区间[64,192]的所有像素置为128
特别说明：
由椒盐噪声产生的噪声数组有三个值：分别对应胡椒噪声的0，盐粒噪声的1和无噪声的0.5。
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-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.5 $ $Date: 2003/10/12 23:37:29 $
% 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 = a*exp(b*randn(M, N));
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 = imnoise2('gaussian',100000,1,0,1);
hist(r,50); %高斯随机数的直方图
>> r1 = imnoise2('uniform',10000,1,0,1);
>> hist(r1,50) %均匀随机数的直方图
>> r5 = imnoise2('rayleigh',10000,1,0,1);
>> hist(r5,50) %瑞利分布直方图
>> r4 = imnoise2('lognormal',10000,1,0.25,1);
hist(r4,50) %得对数正态随机数直方图
>> r6 = imnoise2('exponential',10000,1);
>> hist(r6,50) %指数随机数直方图
（3）周期噪声（函数imnoise3）
说明：一幅图像的周期噪声典型地产生于图像获取过程中的电器和/或电动机械的干扰。这是本章节中唯一考虑的一种空间依赖型噪声。周期噪声的模型是2D正弦波，表达式如下：
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 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)));
实验代码及结果：
C = [0 64;0 128;32 32;64 0;128 0;-32 32];
[r,R,S] = imnoise3(512,512,C);
subplot(1,2,1);imshow(S,[ ]);title('指定冲击的频谱');
subplot(1,2,2);imshow(r,[ ]);title('相应的空间正弦噪声模式');
（4）估计噪声参数
说明：估计周期噪声参数的典型方法是分析图像的傅里叶谱。周期噪声往往产生频率尖峰，频率尖峰可以通过目测来检测。当噪声尖峰足够明显时，或在干扰频率的一些知识可用的情况下，自动分析是有可能的。对于空间域噪声，PDF的参数可以通过传感器的技术来说明，但是通过样本图像来估计它们是很有必要的。
函数statmoments计算平均值和n阶中心距：
[v,unv] = statmoments(p,n) %p是直方图向量，n是计算的矩的数量
对于uint8类图像，p的分量等于28；对于uint16类图像，等于216；对于single和double类图像，等于28或216。
使用函数roipoly选择感兴趣的区域：
B = roipoly(f,c,r) %f是感兴趣图像，c和r是相应多边形的顶点列坐标和行坐标
下面函数计算ROI的直方图：
function [p, npix] = histroi(f, c, r)
% 用于计算一幅图像在多边形区域内的直方图，多边形的顶点由c和r指定，通过函数roipoly复制所定义的多边形区域
% Generate the binary mask image.
B = roipoly(f, c, r);
% Compute the histogram of the pixels in the ROI.
p = imhist(f(B));
% Obtain the number of pixels in the ROI if requested in the output.
if nargout > 1
npix = sum(B(:));
end
实验代码及结果：
f = imread('C:\Users\Public\Pictures\Sample Pictures\flower.jpg');
g = rgb2gray(f); %将彩色图像转为灰度图像
figure,imshow(g);
c=[222 272 300 270 221 194];
r=[21 21 75 121 121 75];
[B,c,r] = roipoly(g); %产生交互式模板B
[h,npix] = histroi(g,c,r);
figure,bar(h,1);
模板B：
结果图：
结论： 用所选用的模板来估计原图像的灰度级分布，显示出图像大概的直方图。
X = imnoise2('gaussian',npix,1,147,20);
figure,hist(X,130);
axis([0 300 0 140]);
（三）仅有噪声的复原——空间滤波
（1）空间噪声滤波器
使用spfilt自定义函数：
function f = spfilt(g, type, m, n, parameter)
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
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'
f = medfilt2(g, [m n], 'symmetric');
case 'max'
f = ordfilt2(g, m*n, ones(m, n), 'symmetric');
case 'min'
f = ordfilt2(g, 1, ones(m, n), 'symmetric');
case 'midpoint'
f1 = ordfilt2(g, 1, ones(m, n), 'symmetric');
f2 = ordfilt2(g, m*n, ones(m, n), 'symmetric');
f = imlincomb(0.5, f1, 0.5, f2);
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
实验代码及结果如下：
f = imread('C:\Users\Public\Pictures\Sample Pictures\flower.jpg');
g = rgb2gray(f);
[M,N] = size(g);
R1 = imnoise2('salt & pepper',M,N,0.1,0);
gp = g;
gp(R1 == 0) = 0;
R2 = imnoise2('salt & pepper',M,N,0,0.1);
gs = g;
gs(R2 == 1) = 255;
fp = spfilt(gp,'chmean',3,3,1.5); %用反调和均值滤波
fs = spfilt(gs,'chmean',3,3,1.5); %用反调和均值滤波
fpmax = spfilt(gp,'max',3,3);
fsmin = spfilt(gs,'min',3,3);
subplot(321);imshow(R1,[]);title('被胡椒噪声污染的图像');
subplot(322);imshow(R2,[ ]);title('被盐粒噪声污染的图像');
subplot(323);imshow(fp,[]);title('滤除胡椒噪声后的图像');
subplot(324);imshow(fs,[ ]);title('滤除被盐粒噪声后的图像');
subplot(325);imshow(fpmax,[]);title('用最大滤波后的结果');
subplot(326);imshow(fsmin,[ ]);title('用最小滤波后的结果');
结论： 被胡椒噪声污染的图像会产生白色背景和黑色颗粒，被盐粒噪声污染的图像会产生黑色背景和白色颗粒。用反调和均值法滤除胡椒噪声，基本可以滤除噪声并还原图像，虽然会产生模糊；但用反调和均值法滤除盐粒噪声的效果较差，只能还原主体图像而背景图像仍然没有被还原。最大值滤波器和最小值滤波器能够比较好的滤除胡椒噪声和盐粒噪声，反调和均值则会模糊图像。
（2）自适应空间滤波器
用adpmedian函数可实现此算法，语法如下：
f = adpmedian(g,Smax) %g是被滤波的图像，Smax是允许的最大自适应滤波器窗口的尺寸
实验代码及结果如下所示：
f = imread('C:\Users\Public\Pictures\Sample Pictures\flower.jpg');
g1 = rgb2gray(f);
g = imnoise(g1,'salt & pepper',0.25);
f1 = medfilt2(g,[7 7],'symmetric');
f2 = adpmedian(g,7);
subplot(131);imshow(g);title('被椒盐噪声污染的图像');
subplot(132);imshow(f1);title('中值滤波器得到的图像');
subplot(133);imshow(f2);title('自适应中值滤波器得到的图像');
结论： 经过对比，自适应中值滤波器比中值滤波器复原效果更加好。中值滤波器会模糊图像，而自适应中值滤波器会更好的还原图像，细节被复原得很清晰，基本与原图像无差别。
（四）退化函数建模
在图像复原问题中，遇到的主要退化是图像模糊。图像模糊可以使用fspecial来建模：
PSF = fspecial('motion',len,theta) %len和theta的默认值分别是9和0
使用函数imfilter来创建已知PSF或用刚才描述的方法计算得到的PSF的退化图像：
g = imfilter(f,PSF,'circular'); %'circular'用来减少边缘效应
之后我们通过添加适当的噪声来完成退化图像模拟
g = g + noise; %noise和g是大小相同的随机噪声图像
用函数checkerboard产生的测试模式对于实现这个目的非常有用，因为大小可以缩放，但不会影响主要特征。语法为：
C = checkerboard(NP,M,N) %NP是每个正方形一边的像素数，M是行数，N是列数
用下面的命令可以产生亮正方形全是白色的棋盘板：函数checkerboard产生的图像属于double类，值在区间[0,1]内。
R = checkerboard(NP,M,N) > 0.5;
实验代码及结果：
f = checkerboard(8); %产生棋盘板图像
figure,subplot(221),imshow(f);title('原始图像');
PSF = fspecial('motion',7,45);
gb = imfilter(f,PSF,'circular'); %产生退化图像
subplot(222),imshow(gb);title('使用motion滤波器模糊后图像');
noise=imnoise(zeros(size(f)),'gaussian',0,0.002);%产生高斯噪声
subplot(223),imshow(noise,[]);title('高斯纯噪声');
g=gb+noise;
subplot(224),imshow(g,[]);%模糊加噪声后的图像
title('模糊加噪声后的图像');
（五）直接逆滤波
逆滤波，就是简单将退化函数去除，直接的逆滤波没有什么意义，只处理了靠近直流分量的部分，其他不做处理。由退化函数H退化的图像复原的最简单的方法是直接做逆滤波，设图像退化前的傅里叶变换为F(u,v)，退化后的傅里叶变换为G(u,v)，系统函数即退化函数的傅里叶变换为H(u,v)。
所谓直接逆滤波，就是用退化函数除退化图像的傅里叶变换，得到退化前图像的傅里叶变换的估计。公式如下所示：
由该式可知，即使知道退化函数，也不能准确的复原图像，因为N(u,v)未知，甚至有更糟的情况是如果退化函数是零或是非常小的值，则N(u,v)/H(u,v)的值比较大，很容易支配F(u,v)的估计值。在无噪声的情况下，逆滤波是完美的。但实际中，所有都会噪声，因此这个类型逆滤波的方法很少使用到。
（六）维纳滤波
维纳滤波也称为最小均方差误差滤波，它是一种最早、也是最为熟知的线性图像复原方法。这种滤波就是使图像尽可能的平滑，可以消除很严重的噪声，恢复图像。
维纳滤波是通过函数deconvwnr实现的，语法如下：（g是退化图像，frest是复原图像）
frest = deconvwnr(g,PSF); %信噪比为0
frest = deconvwnr(g,PSF,NAPR); %信噪比已知或是常量或是数组，NAPR是标量输入
frest = deconvwnr(g,PSF,NACORR,FACORR); %噪声和未退化的图像的自相关函数NACORR,FACORR是已知的
实验代码及结果：（对上面退化模型进行滤波）
subplot(221),imshow(g,[]);%模糊加噪声后的图像
title('模糊加噪声后的图像');
frest1=deconvwnr(g,PSF); %维纳滤波
subplot(222),imshow(frest1,[]),title('维纳滤波后的结果');
Sn = abs(fft2(noise)).^2; %噪声功率谱
nA=sum(Sn(:))/prod(size(noise)); %平均噪声功率谱
Sf = abs(fft2(f)).^2; %信号功率谱
fA=sum(Sf(:))/prod(size(f)); %平均信号功率谱
R = nA/fA; %求信噪比
frest2=deconvwnr(g,PSF,R);
subplot(223),imshow(frest2,[]),title('使用常数比例的维纳滤波');
NACORR = fftshift(real(ifft(Sn)));
FACORR = fftshift(real(ifft(Sf)));
frest3=deconvwnr(g,PSF,NACORR,FACORR); %自相关后的维纳滤波
subplot(224),imshow(frest3,[]),title('使用自相关函数的维纳滤波');
结论：自相关函数的维纳滤波效果比较好，简单的维纳滤波效果最差。自相关函数的维纳滤波并没有很好地完全还原图像，而且还有一些噪声存在；相比之下，简单的维纳滤波存在的噪声更加多，效果比较差，亮度也没有得到还原。
1098
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
