3.2.2　关于图像平滑的最优停止时间

在图像处理过程中，噪声去除和边缘保护一直以来都是图像平滑过程中的一对矛盾．无论是直接的傅里叶（Fourier）变换或者其他滤波变换，还是基于偏微分方程或其他几何等形式的滤波方法，去噪过程都存在数据迭代的过程，即每一次在上一次平滑的基础上继续平滑．这种迭代需要处理的一个重要问题就是何时停止迭代．当然，有的去噪程序或去噪模块会设计一些自适应停止，但是这种自适应算法本身也就包含对停止时间的判断，如上一节的滤波过程就需要考虑停止时间．

1）判断平滑停止时间几种方法

这里主要介绍基于偏微分方程的图像平滑，其最优停止方法主要有两大类：一类必须知道无噪声理想图像或者图像噪声的相关信息，这一类通常用于理论分析；另一类用于解决实际问题，不必也无法知道图像噪声信息，所以实际应用性比较强．

首先，回顾一种基于计算相对方差的停止方法[41]

式中，f是无噪声图像，un是迭代n次后的滤波图像，var (u)表示图像函数u的方差，T是最佳迭代次数．而最小值点函数

其次，Weickert[41]介绍了噪声信噪比方法，提到如果噪声图像u0是由无噪声图像f加上加性噪声ne构成的，即u0=f+ne，那么，最佳停止次数可以利用信噪比SNR由下式计算：

然后，Solo在2001年[42]建议在迭代次数k处定义一个风险函数

其离散形式为

式中，E为均值，也即数学期望；k为平滑去噪迭代次数．理论上，这里应通过Rk取得最小值而得出停止时间k∗（k∗=k）．不过这里的困难在于f存在未知性．

Barcelos[43]在2003年提出了一种基于高斯（Gauss）分布的最佳停止时间判定方法．他在经过实验和数据分析后认为，最佳迭代次数N=σ2/（a·Δt），其中σ是零均值高斯分布的标准差，Δt是离散化的时间步长，a是高斯核

中的一个常数，一般取a=2或100（实际上，a只在使用热传导方程平滑时有具体的物理意义）．这种方法的缺点在于，对于真实图像而言，图像噪声信息难以被获取．

这里不再一一列举更多确定停止时间的理论分析方法．

2）平滑停止时间的计算

平滑停止时间的计算主要存在如下方法：

（1）基于相关系数法．

Mrázek和Navara[44]于2003年建议一种基于相关系数最小时的迭代停止准则，李世飞等[45-46]在这方面有详细分析．基本思想如下：

在图像噪声为加性噪声（即噪声与图像不相关）的前提下，图像被噪声污染的过程就可以用f+ne=u0表示．这里，u0为初始图像，f为理想的无噪声图像，ne为噪声．扩散滤波过程是从u0开始，依次迭代，第一次得到滤波图像u1，第二次得到u2，…，这样得出一列滤波后的图像列{un}．

再设ut为t次迭代获得的最佳图像，用u0-ut表示图像噪声，这样u0-ut和ut应存在最小的相关系数，此时t所对应的迭代次数即为所求．最小的相关系数计算式如下：

其中

这里，协方差，方差，均值．

从计算上看，这里不涉及无噪声图像，因此不影响计算．

（2）基于互信息的方法．

互信息是两个随机变量共同信息量的一种度量．互信息刻画最佳迭代次数的方法被郭圣文等[47-48]提及，其基本思想同用相关系数来刻画最佳迭代次数相似，最佳迭代次数对应互信息量最小的点．互信息的计算公式如下：

式中，p（x, y）是随机变量u0-ut和ut的联合分布概率密度函数；p（x）和p（y）分别是u0-ut和ut的边际分布概率密度函数．p（x）和p（y）可以通过计算p（x, y）得到．特别地，p（x, y）可以通过计算的联合直方图得到．

用二维随机变量（X, Y）的样本（x1，y1），（x2，y2），…，（xN, yN）构造估计联合概率密度的联合两维直方图[49]．二维联合直方图的X、Y轴代表随机变量X、Y的取值．对于图像而言，随机变量X、Y是图像的灰度，二维联合直方图是X、Y离散的图像灰度值，取值范围为[0, 255]，h（x, y）是两幅图像重叠部分的灰度值为（x, y）的像素的总个数．

同利用相关系数刻画最佳迭代次数相比，互信息方法不必事先对图像或噪声的性质做先验假设．具体的数值实验可以显示，基于相关系数或互信息这两种方法均能够获得良好的滤波效果．尽管这两种方法对各种不同的噪声均适用，但是由于两种方法在具体离散数值计算时采用的算法不同（基于相关系数的方法对图像的运算操作与图像尺寸有关，图像尺寸越大则计算量越大；互信息方法由于采用联合直方图来计算，运算操作与图像灰度级有关），在对图像进行同样扩散滤波时，互信息在较大尺寸图像计算方面用时较短．

3）数值实验

前面介绍了相关系数和基于互信息的最佳迭代停止方法，可以用于线性模型或者非线性各向异性扩散模型的平滑．这里采用通用性较强的热传导方程作为扩散模型来验证两种方法的有效性．在具体计算中，热传导方程采用中心差分格式进行离散[50]．

如图3-5所示为一个典型的使用热传导方程平滑去噪过程．图3-5（a）是被高斯噪声污染的图像；图3-5（b）是利用互信息判断最佳迭代次数时输出的图像（60次迭代）；图3-5（c）是热传导方程迭代作用于图像100次的图像，出现了过度的平滑与边缘消失．从图3-5（b）直观看，互信息方法在滤除噪声的同时基本上保护了边缘不被模糊．如图3-6所示为另一个算例．图3-6（a）是被椒盐噪声污染的图像；图3-6（b）是利用互信息最小时停止迭代后输出的图像（37次迭代）；图3-6（c）是热传导方程迭代作用于图像100次的图像．

图3-7中的曲线表示互信息和迭代次数之间的关系，使用的是热传导方程去噪方法，经过200次迭代．图3-7（a）与图3-5（b）对应，图3-7（b）与图3-6（b）对应．由图3-7可以看出曲线显示一个互信息量最小点，此点即为所求．曲线所显示的这种变化也与图3-5和图3-6中由（a）到（c）的变化一致．

如图3-8所示为使用互信息方法与相关系数方法进行噪声去除的效果比较．图3-8（a）是高斯噪声污染后的图像；图3-8（b）是利用互信息方法确定的最佳滤波图像，此种方法确定的最佳迭代次数为41次；图3-8（c）是利用相关系数方法确定的最佳滤波图像，此种方法确定最佳迭代次数为38次．

如图3-9所示为互信息量数值和相关系数随着迭代次数的变化情况，整个过程经过200次迭代．图3-9（a）对应图3-8（a）在去噪过程中互信息量与迭代次数的变化情况；图3-9（b）是与图3-8（a）对应的相关系数与迭代次数的变化情况．从图3-9可以看出，滤波图像效果相当，最佳迭代次数基本相同，但就所用时间来说互信息方法更优．这是因为计算互信息算法中的循环迭代次数与图像256个灰度级有关，与图像尺寸无关；而相关系数算法中的循环迭代次数与图像尺寸有关，图像尺度越大，循环迭代次数越多．另外，从与迭代次数的关系图上看，两种方法在最佳迭代次数与图形曲线变化上也是一致的．

进一步对信噪比SNR和峰值信噪比PSNR的分析，以及其他分析说明，可以参见李毅等[50]的工作，这里不再一一列举和深入说明．