chap4 图像复原 part1

Reference：《数字图像处理MATLAB版》冈萨雷斯

1.噪声模型

常见的噪声模型（z是随机变量，灰度值； $\mu$ 表示均值， ${s}^{2}$ 表示方差）：
1、高斯噪声
2、瑞里噪声其中， $μ=a+\sqrt{(\frac{\pi*b}{4})},s^2=\frac{b(4-\pi)}{4}$
3、伽马噪声其中， $\mu=\frac{b}{a}$ , $\sigma^2=\frac{b}{a^2}$
4、指数分布噪声其中， $μ=\frac{1}{a},\sigma^2=\frac{1}{a^2}$
5、均匀分布噪声其中， $μ=\frac{(a+b)}{2},\sigma^2=\frac{(b-a)^2}{12}$
6、椒盐噪声其中盐表示亮点，椒表示暗点。
7、周期噪声，比如空间域图像受到正弦波信号干扰。一幅图像的周期噪声是在图像获取期间由电力或机电干扰产生的
纯正弦波的二维DFT： $sin(2\pi\mu_0x+2\pi\nu_0y) \rightleftharpoons \frac{1}{2}j[\delta(\mu + M\mu_0, \nu +N\nu_0)-\delta(\mu+M\mu_0,\nu+N\nu_0)]$
而， $\mu\in[0,M],\nu\in[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的左右两侧分别剪去了 $\frac{w-1}{2}$ 个元素。

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)));