去除高斯噪声的图像处理技术.doc

去除高斯噪声的图像处理技术
-------验证《image enhancement based on noise estimation》
一．概述
大家都知道这样的一个事实，对图象躁声的处理存在一个平滑和锐化的矛盾。躁声在图象中对于高频部分的贡献比正常数据要大，也就是说它与它周围的象素存在“突变”，这就是我们平滑的原因。但是“突变”也可能是边界，只进行平滑可能会模糊边界，得到不满意的结果。“锐化”就是为了突显边界，这两者之间存在一个权衡。我们以8位的BMP图，和未知的高斯噪声为例子。为了能有效的验证我们算法的正确性，我们先写给出制造高斯躁声的程序（可以调节方差），然后给出有效的消除高斯躁声的算法。我们程序的编译运行的环境是WindowsXP系统，VC++2005。
二. 设计思路
1.制造高斯躁声的设计：
要研究图像的增强与恢复，就必须先了解噪声，高斯噪声是自然界中存在最广泛的噪声，因此也成为我们本次实验的对象噪声，其振幅分布服从高斯分布即正态分布。含有高斯噪声的图像，可以看成是图像的每一个像素值加上一个高斯分布的随机数。因此，要生成含高斯噪声的图像，就必须先产生高斯分布的随机数，然后再将随机数加入图像。
根据box-muller算法，可以利用极坐标的原理，由均匀分布的随机数生成标准正态分布的随机数。但经测试，系统程序库中的随机函数rand生成的假随机数，在分布精度上并不满足均匀分布的统计特性。故需要重新设定随机种子与随机算法，以生成较为严格的均匀分布随机数。
2.消除高斯噪声的设计：
这个算法是用于图象增强，但其设计是源于图象恢复，其实本质上也是种图象恢复。我们首先要估计出在原图象上的高斯躁声，用较为准确的估计值确定参数。我觉得最为关键的点，就在于方差估计算法的设计。在确定好参数以后，我们通过两次不同的平滑和一次锐化，可以有效的消除高斯躁声。具体的设计可以从算法设计中看出。
三. 算法设计
1.制造高斯躁声的算法：
生成含高斯噪声的算法可以转化为两个算法的实现：第一，生成均匀分布随机数；第二，实现均匀分布随机数向高斯分布随机数转化的box-muller算法。
系统程序库中的随机数发生器rand，实际上是一个函数子过程，通过线形同余法，即采用递推关系，这样使得计算量变得较少，但也因此造成了产生的随机数的序列相关性，因此，要设计或者改进这种随机数发生器，就必须破坏这种序列相关性。在本实验中，我们采用混洗的方法破坏序列相关性，即设V1,V2,V3,…,Vn是由rand产生的n个随机数，现随机的取一正整数j（1<=j<=n），取Vj为一要求的随机数，而Vj再由rand生成的另一随机数替换，替换后再由V1,V2…,Vn中随机的取一个为下一次要求的随机数，依此重复。为进一步改进产生随机数的统计特性，本实验函数ran1(idum)采用三个线形同余发生器，公共组成，其中，第一个线形同余发生器用于产生随机数的最高有效位部分，第二个用于产生随机数的最低有效位部分，第三个用于控制混洗过程。经测试，ran1(idum)所产生的随机数较好的体现了平均分布的统计特性。具体实现方法，及代码，见代码部分。
用变换的方法可以实现均匀分布的随机数向正态分布的转换。根据box-muller算法，设是的均匀分布随机函数，则可作变换：

等价的，可得：

因此，是两个独立的标准正态分布随机数。
进一步，考虑到三角函数的计算量大，因此，把取换为在单位圆内取随机点，由代替，而该点与轴的角代替2，则：

现设是区间上的独立的均匀分布随机数，且：

显然，是上的均匀随机数，若，则：

是两个独立的标准正态分布随机数。具体实现方法，及代码，见代码部分。
基于8位图像的特点：每个像素由256个灰度等级表示，故需要按照特定的方差，将产生的正态分布随机数分进512个区段中。因此，在程序中添加了cint(double x)函数和方差的相关常量0.05后，完成了基于box-muller和极坐标理论生成的正态分布随机数在8位图像中的高斯噪声添加工作。

2.消除高斯噪声的算法：
在我们的算法里面，主体分为两大块：
我们先设定象素值L=256。
a.高斯躁声参数的估计：
平均值u：常被认为是0。、
方差σ：所要估计的。先给出几个涉及到的方程：
算法的主要思路：
处理的象素集合：
公式1：
公式1为一个滤波器，其中，b,c取经验值（b=1,c=3），选择滤波器中的参数a（我们选定的范围是从1到100）用于求解噪声方差。
公式2：
公式3：
公式4：
当选择一个a的值后，式1滤波结果图像（记为G(a)），按式2计算G(a)的边界梯度，统计边界梯度图的直方图（H（g,a）中g指梯度（由于求梯度前做过滤波，就与a有关了），Hmax(a)为该直方图的最大值（公式3），当然不同的a选择就会有不同的G(a)，当然也就会有不同的梯度图，直方图也随之不同，其中获得Hmax最大值的那个a选择，就是公式4中的a0。得到了a0以后，除以2就得到我们估计的方差。（在我们后面的验证中可以看到该方法很有效）。
在我们程序中函数是：int estimate()和long H(BYTE* y, int a0)
主要的过程：
while(a<100)
{
for(i=1; i<m_ImageHeight-1; i++)
{
for(j=1; j<m_ImageWidth-1; j++)
{
temp[0]=Bt(m_imagedata[i*m_ImageWidth+j],m_imagedata[i*m_ImageWidth+j-1],a,1,3);
temp[1]=Bt(m_imagedata[i*m_ImageWidth+j],m_imagedata[i*m_ImageWidth+j+1],a,1,3);
temp[2]=Bt(m_imagedata[i*m_ImageWidth+j],m_imagedata[(i-1)*m_ImageWidth+j],a,1,3);
temp[3]=Bt(m_imagedata[i*m_ImageWidth+j],m_imagedata[(i+1)*m_ImageWidth+j],a,1,3);

y[i*m_ImageWidth+j]=m_imagedata[i*m_ImageWidth+j]+( temp[0]+temp[1]+temp[2]+temp[3])/4;
}
}
a++;
long q= H(y,a);
if( max<q )
{
max=q;
a0=a;
}
}
估计出方差σ后，可以根据以下的式子计算出参数（由经验得到的）：
a1=(19*σ-29)/6；
a2=σ/2；
d=(5/4*σ+25/4)；
b.主体计算：
在程序中的函数是：void Process( float a1, float a2, float d )
涉及到共三个过程（按顺序执行）：
1.平滑：
这里提到了Bt函数：float Bt(BYTE soursData, BYTE aroundData, float a, float b, float c)，这个函数实际是个二维图：
float Bt(BYTE soursData, BYTE aroundData, float a, float b, float c)
{
int i=-soursData+aroundData;
if( i<=-c*a || i>=c*a )
{
return 0;
}
else if( i>-c*a && i<=-a*b )// || ( i>a && i<b*a )
{
return (i+c*a)/(b-c);
}
else if( i>-b*a && i<=-a )
{
return -a;
}
else if( i>-a && i<=a )
{
return i;
}
else if( i>a && i<=b*a )
{
return a;
}
else if( i>b*a && i<c*a )
{
return (i-c*a)/(b-c);
}
return 0;
}
有了这个函数以后，我们可以根据以下的方程处理：
先给出四个处理的象素集合：

，
其中：，y是处理后的象素值（其对应的位置是I,j），x是处理前的象素值（其对应的位置是I,j，原图象的象素），a1可以根据预估计的σ（公式是a1=(19*σ-29)/6）计算得到，可以根据经验设定b1=1.5, c1=3。
经过计算后，我们进入下一步处理。
2.平滑：
用的是和以上一样的函数float Bt(BYTE soursData, BYTE aroundData, float a, float b, float c)。
象素计算的集合和第一步平滑一样。
可以根据这个方程进行我们的第2步平滑：

其中，y是处理后的象素值（其对应的位置是I,j），x是处理前的象素值（其对应的位置是I,j，上次处理的结果），a2可以根据预估计的σ（公式是a2=σ/2）计算得到，可以根据经验设定b2=1.5, c2=3。
经过计算后，我们进入下一步处理。
3.锐化：
这里提到了Bt函数：float Gm(BYTE soursData, BYTE aroundData, float d)，这个函数实际是也个二维图：
float Gm(BYTE soursData, BYTE aroundData, float d)
{
int i=-soursData+aroundData;
int L=256;
int e=128;
float m1=2;
float m2=0.5;
float g=2L-2+4*d-3*e;

if( i<=-g )
{
return L-1;
}

else if( i>-g && i<=-e )
{
return (-0.5*i+1.5*e-2*d);
}

else if( i>-e && i<=-d)
{
return (-2*i-2*d);
}

else if( i>-d && i<=d)
{
return 0;
}

else if( i>d && i<=e)
{
return (-2*i+2*d);
}

else if( i>e && i<g)
{
return (-0.5*i-1.5*e+2*d);
}

else if( i>=g )
{
return (-L+1);
}
return 0;
}
得到了Gm函数以后，可以根据以下的函数：
处理的象素集合是：
在这个集合上，运用公式；

其中y是处理后的象素值（其对应的位置是I,j），x是处理前的象素值（其对应的位置是I,j，上次处理的结果），d可以根据预估计的σ（公式是d=(5/4*σ+25/4)）计算得到，其他的参数可以根据e=128，m1=f/(e-d)=2,m2=(L-1-f)/(g-e)=0.5得到，m1是处理前的象素值和处理后的象素值差的绝对值在d和e 之间时的直线斜率；m2是处理前的象素值和处理后的象素值差的绝对值在e和g 之间时的直线斜率，可以从我们的float Gm(BYTE soursData, BYTE aroundData, float d)函数中看出来。

经过以上处理以后，我们可以得到较好的消除高斯躁声的图象。
四. 实验结果及其说明：
测试1：
图片1（没有躁声，也就是说σ=0）

经过我们的填加躁声处理后，得到图片2（设定的σ=16）：

经过我们的预测函数，预测得到的σ=16，可见该函数还是比较准确的：

经过处理后的图象,σ=4：

注：
可以得到两点：
1.我们的填加高斯躁声的程序是正确的，欲估σ的函数也是比较有效的；
2.针对高斯躁声的图象加强的方法虽然不能将图片恢复到原样，但也取得了比较好的结果，证明了算法的有效。
测试2：
选用了原论文上的图片，其原文上给出的σ=10：

我们的预估函数得到的σ=10，再次验证该函数还是比较准确的：

经过处理后的图象,σ=3：

再次验证算法的有效
五. 总结
在开始阅读国外论文的时候，有些地方不是很了解，经过思考也还存在些问题，尤其是在预估计这个地方，最后请教老师，感到有豁然开朗的感觉。做这个程序，我一直很感兴趣，因为感到一种直观的结果，很有意思。我试图给出另一些算法，用于和该算法进行对比，但在进行滤波处理后，用一次微分锐化总觉得很过了，反观这个算法，它并没有单一的处理，不是一条直线，所以能很好的避免锐化过分的问题，我想这也就是float Gm(BYTE soursData, BYTE aroundData, float d)函数构造的思路。想想自己滤波的处理也并不好，这就是float Bt(BYTE soursData, BYTE aroundData, float a, float b, float c)和两次不同角度的平滑的原因，这样也有利于后面的锐化，这只是我的想法。
六. 参考资料
《image enhancement based on noise estimation》 ----------------------------------Fabrizio Russo
《数字图像处理第二版》---------------------------------------------------------------------冈萨雷斯