五、图像的复原与重建（MATLAB）上

一、图像退化/复原处理模型

原始图像为f(x,y),退化过程建模为一个退化函数H和一个加性噪声k(x,y)，退化之后得到的图像为g(x,y),数学模型表达式为：
g(x,y)=H[f(x,y)]+k(x,y)
图像复原的目的就是通过退化函数，以及加性噪声来得到初始图像的估计值f’(x,y)。
如果 H 是一个线性的，空间不变的过程，那么空域中的表达式可以写成
g(x,y)=h(x,y)卷积f(x,y)+k(x,y)
等价的频域表达式为：
G(u,v)=H(u,v)F(u,v)+K(u,v)

二、噪声模型

模拟噪声特性和影响的能力是图像复原的核心。
1，使用函数imnoise对图像添加噪声
函数语法为：
g=imnoise(f,type,parameters)
其中，f是输入图像，type是噪声类型，parameters指一些参量。函数 imnoise在给图像添加噪声之前，将图像转换为 [0，1] 中的double类。指定噪声参数时，必须考虑这一点。例如，要将均值为64，方差为400的高斯噪声添加到一幅uint8类图像上，需要将均值标度为64/255并将方差标度为400/（255^2）后，作为imnoise的输入。
g=imnoise(f,‘gaussian’,m,var) 将均值m,方差为var的高斯噪声添加到原图像上。
g=imnoise(f,‘localvar’,v) 将均值为0，局部方差为v的高斯噪声添加到源图像上。其中V是一个大小与f相同的数组，它在每个点包含了期望的方差值。
g=imnoise(f,‘localvar’,image_intensity,var) ,image_intensity是自变量，局部方差var是图像灰度值的函数。
g=imnoise(f,‘salt&pepper’,d) 用椒盐噪声污染图像f,其中d是噪声密度，默认值为0.05。
g=imnoise(f,‘speckle’,var) 使用方程g=f+n.*f将乘性噪声添加到f上,其中n是均值为0，方差为var(默认为0.04)的均匀分布的随机噪声。

subplot(121),imshow(f)
>> g=imnoise(f,'gaussian');
>> subplot(122),imshow(g)

添加高斯噪声到原图像上，参数默认值是均值为0，方差为0.01。
再保持均值不变，方差为0.1；以及保持方差不变，均值为0.5。

 f=imread('C:/experiment/test1.jpg');
>> f=rgb2gray(f);
>> subplot(221),imshow(f)
>> g1=imnoise(f,'gaussian');
>> subplot(222),imshow(g1),title('噪声均值为0，方差为0.01')
>> g2=imnoise(f,'gaussian',0,0.1);
>> subplot(223),imshow(g2),title('噪声均值为0，方差为0.1')
>> g3=imnoise(f,'gaussian',0.5,0.01);
>> subplot(224),imshow(g3),title('噪声均值为0.5，方差为0.01')

观察图像：


观察对比图像，保持均值不变，方差为0.1（方差增加），噪声干扰更大，分辨率降低，图像更加模糊。保持方差不变，均值为0.5（均值增加），图像更白。
换一副图像验证结果。

观察对比可知，上述结论依然适用。

下面通过rand函数生成一个随机数组最为V，添加噪声到原图像

 V=rand(size(f));
 g1=imnoise(f,'localvar',V);

2，使用规定分布生成空间随机噪声
由概率论知识可知，如果w是在区间（0，1）内均匀分布的一个随机变量，通过求解下面的方程，就可以得到一个具有规定（累积分布函数）CDF和F的随机变量z：

假设有一个在(0,1)内的均匀随机数w的生成器，并假设我们要用该生成器来生成具有瑞丽CDF的随机数z,它有如下形式：

式中，b>0。为求解z,解方程

等价为

运用这个结果，在MATLAB中，很容易生成一个随机数的数组R

 R = a + sqrt (b*log(1-rand(M,N)));

z的表达式有时称为随机数生成器方程，它确定了生成期望随机数的方式。下表列出了当前感兴趣的随机变量，及其PDF,CDF，随机数生成器方程。

在像以低照明度运行的成像传感器中，高斯噪声看作其中的近似；在不完善的开关设备中，会出现椒盐噪声；感光乳剂中的银粒大小是由对数正态分布描述的一个随机变量；深度成像中会出现瑞丽噪声，而在描述激光成像中的噪声时，指数和爱尔兰噪声有用。

定义一个M函数 imnoise2，生成一个M*N的噪声数组，生成的是噪声模式本身。

function R=imnoise2(type,varargin)
[M,N,a,b]=setDefaults(type,varargin{:});
switch lower(type)
    case 'uniform'
        R=a+(b-a)*rand(M,N);
    case 'gaussian'
        R=a+b*randn(M,N);
    case 'salt&pepper'
        R=saltpepper(M,N,a,b);
    case 'lognormal'
        R=exp(b*randn(M,N)+a);
    case 'rayleigh'
        R=a+(-b*log(1-rand(M,N))).^0.5
    case 'exponential'
        R=exponential(M,N,a,b);
    case 'erlang'
        R=erlang(M,N,a,b);
    otherwise
        error('unkonow distribution type')
end

查看各种噪声直方图（其中参数是函数中的默认值）

注意，由椒盐噪声生成的噪声数组有三个值：对应于胡椒噪声的0，对应于盐粒噪声的1，以及对应于无噪声的0.5。为了使这个数组有用，还需要对他进一步的处理。椒盐噪声污染图像时，需要找到R中所有的0值坐标，把相应原图像的坐标位置置为最小灰度值，找到R中值为1的坐标，把相应原图像的坐标位置置为最大灰度值。代码实现为：

function R=saltpepper(M,N,a,b)
if (a+b)>1
    error('the sum must not exceed 1')
end
R(1:M,1:N)=0.5;
X=rand(M,N);
R(X<=a)=0;
u=a+b;
R(X>a&X<=u)=1;

运用椒盐噪声污染图像：

[M,N]=size(f);
>> R=imnoise2('salt&pepper',M,N,0.1,0);
>> gp=f;
>> gp(R==0)=0;
>> subplot(121),imshow(f);
>> subplot(122),imshow(gp);

图像是被概率只有0.1的胡椒噪声污染的uint8类图像，图像上出现大量黑点（胡椒状）。

[M,N]=size(f);
R=imnoise2('salt&pepper',M,N,0,0.1);
gp=f;
gp(R==1)=255;
subplot(121),imshow(f);
subplot(122),imshow(gp);

图像是被概率只有0.1的盐粒噪声污染的uint8类图像，图像上出现大量白点（盐粒状）。
3，周期噪声
图像出现周期噪声通常源于电气，电机的干扰，通常通过在频域滤波来处理。周期噪声的噪声模型是一个离散的二维正弦波，方程为：

其中，x=0,1,2，…，M-1且y=0,1,2,…,N-1；A 是振幅，u0和 v0 分别确定了关于x 轴和y 轴的正弦频率。Bx 和 By分别是关于原点的相移。
定义一个M函数接受任意数量的脉冲位置（频率坐标），每个脉冲位置都有自己的振幅，频率和相移参数。该函数还能够输出各个正弦波之和的傅里叶变换，以及谱。

function[r,R,S]=imnoise3(M,N,C,A,B)
K=size(C,1);
if nargin<4
    A=ones(1,K);
end
if nargin<5
    B=zeros(K,2);
end
R=zeros(M,N);
for j=1:K
    u1=floor(M/2)+1-C(j,1);
    v1=floor(N/2)+1-C(j,2);
    R(u1,v1)=i*M*N*(A(j)/2)*exp(-i*2*pi*(C(j,1)*B(j,1)/M)+C(j,2)*B(j,2)/N);
    u2=floor(M/2)+1+C(j,1);
    v2=floor(N/2)+1+C(j,2);
    R(u2,v2)=-i*M*N*(A(j)/2)*exp(i*2*pi*(C(j,1)*B(j,1)/M)+C(j,2)*B(j,2)/N);
end
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(321),imshow(S,[]),title('规定脉冲的谱')
>> subplot(322),imshow(r,[]),title('空域中相应的正弦噪声模型')
>>  C=[0 32;0 64;16 16;32 0;64 0;-16 16];
>> [r,R,S]=imnoise3(512,512,C);
>> subplot(323),imshow(S,[])
>> subplot(324),imshow(r,[])
>> C=[6,32;-2 2];
>> [r,R,S]=imnoise3(512,512,C);
>> subplot(325),imshow(r,[])
>> A=[1,5];
>> [r,R,S]=imnoise3(512,512,C,A);
>> subplot(326),imshow(r,[])