在图像分析领域中，二维图像目标物体的边界周长是一个十分重要的目标形态特征，并且在现实生活中具有广泛的应用。一方面，目标边界的模糊性和复杂性给计算带来了很大的困难；另一方面，噪声等外在因素也会干扰计算。因而，不少学者都围绕这个问题展开了研究，探讨如何精确地计算目标物体的边界长度。目前，大部分文献提出的方法都是基于二值图像，即图像中只存在表示目标和背景的两类值。这类方法较常见的有数字化直线段片断法(digital straight segment, DSS)^[1]、最小多边形法(minimum length polygon, MLP)^[2]和基于B样条曲线的曲线拟合方法^[3-6]。这些方法估算目标周长时，虽然最终能获得较精确的估算结果，但是需要先提取目标的边界信息，这个过程需要将一般图像转换成二值图像，在转换过程中非二值图像的目标边界信息不可避免会产生丢失，最终导致估算结果的不精确。相比二值图像，数字化灰度图像包含的信息更加丰富^[7-8]，因而直接用灰度图像进行能使周长估算的结果更符合实际。

Sladoje和Lindblad提出了一种基于灰度图像的高精度曲线长度估算方法，简称HPGL方法^[9-11]。HPGL方法本身对分辨率较低和目标存在模糊的图像的计算结果并不理想，其主要原因是对于边界厚度超过1个像素宽度的目标边界可能造成边界重复计算的问题。对于这两种情况，文献[9]提出了一种基于图像形态学的像素覆盖双阈值预处理方法，这种方法虽然避免了重复计算问题，但是没有合理利用边界像素信息，估算结果误差较大。文献[12-13]分别提出基于边界跟踪的BTHPGL方法和基于图像粒的系列方法，两类方法改进了HPGL方法在边界模糊情况下的不足，在简单的像素覆盖预处理的基础上，实现对目标周长较精确的估算。BTHPGL方法与基于图像粒的方法在像素覆盖数字化处理后，图像并不能在所有情况下都排除噪声的干扰，噪声可能导致方法精确度下降甚至失效。并且这些方法并没有考虑目标边界像素之间的关系，可能导致最终估算结果精确度不高。

为了进一步提高算法的适应性和准确性，本文提出了一种基于图像森林变换^[14]的周长估算方法，将数字图像看作包含节点和弧的图结构，充分利用像素之间的灰度关系，优化像素覆盖数字化后的图像的边界信息，结合经典的周长估算方法，得出最终的估算结果。本文设计的实验，验证了算法的有效性。

1 相关工作 1.1 像素覆盖数字化

像素覆盖数字化^[9]是由文献提出的高效的图像分割方法。其核心思想是待分割目标中的像素所占百分比和目标所覆盖面积通常成正比。

定义1   假定一个连续目标S⊂R²，投影到2维数字化网格平面上，每个网格对应一个像素p_{i, j}, 目标S的像素覆盖数字化定义为

$ D\left( S \right)\approx \left\{ \left( \left( i, \text{ }j \right), \frac{A({{p}_{i, j}}\cap S)}{A({{p}_{i, j}})} \right)\left| \left( i, \text{ }j \right)\in \text{ }{{Z}^{2}} \right.) \right\} $

(1)

式中：A(X)表示集合X的面积；Z²表示2维数字化网格平面。

上述定义描述的是分配给网格像素的真实覆盖值，而实际在计算机中处理的图像会受灰度级限制，一般图像处理系统中灰度共256级，因此引出下列定义。

定义2   假定一个连续目标S⊂ R²，投影到2维数字化网格平面上，并给定一个n级的量化标准。量化级为n的目标S的像素覆盖数字化定义为

$ D\left( S \right)\approx \left\{ \left( \left( i, \text{ }j \right), \frac{1}{n}\left\lfloor n\frac{A({{p}_{i, \text{ }j}}\cap S)}{A({{p}_{i, \text{ }j}})}+\frac{1}{2} \right\rfloor \right)\left| \left( i, \text{ }j \right)\in {{Z}^{2}}) \right. \right\} $

(2)

式中$\left\lfloor x \right\rfloor $表示小于等于x的最大整数。

定义2将分配给网格的像素值限制在一个n级的量化数字化集合中。特别地，当n=1时表示二值图像，n=255时表示8位的灰度图像。由于目标与每个像素网格相交面积无法计算，式(1)、式(2) 只能用于理论描述。因此本文使用文献[10]中的双阈值分割算法来获得目标图像的像素覆盖数字化。

1.2 DSS、MLP、HPGL和粒度化周长估算方法

DSS (digital straight segment)方法^[1]是一种将数字曲线划分为数字直线片段的几何方法。该方法用数字化直线片段取代数字化边界曲线，按一定规则使这些片段尽量长。DSS方法主要应用于具有封闭边界的目标区域的周长计算。对于一个目标区域，把像素看作网格正方形结构，当其边和顶点作为边界时，只要目标区域连通，那么所求的目标区域的边界就是闭合的。根据方向偏转判别公式确定DSS点，然后累加相邻DSS点的距离，得到最终周长。

MLP (minimum length polygon)方法^[2]与DSS方法相似，也是一种几何方法，其核心思想是找到最小多边形使其包含数字化目标区域，初始化条件与DSS方法基本一致，不同的是在确定每个特征点的同时需要确定该点与该点相邻的两点连线的相对位置，根据这个位置确定一个凸性。最后，根据MLP点的位置和凹凸属性累计局部线段长度，得到目标周长。

HPGL^[8]方法利用数字化图像中的灰度级信息进行边界周长估计。与上述两种求数字图像周长方法相比，该方法能达到更高的计算精度和准确度。HPGL方法以半平面上直线段的长度计算理论为基础，根据像素覆盖数字化图像中分配给每个像素的量化灰度值，进行局部周长估计，其局部周长计算依据公式：

$ l=\int_{a}^{b}{\sqrt{1+{{\left| {f}'\left( x \right) \right|}^{2}}}\text{d}x} $

(3)

式(3) 可计算曲线在区间[a, b]上的曲线长度。

最后，对每一次局部计算求和并进行量化误差的矫正，得到最终估算结果。HPGL方法在单像素宽度目标边界的周长计算应用上能取得很高精度的结果。

2 结合IFT的目标周长估算

第1.1节阐述了像素覆盖数字化的简易方法，并用图片直观展示了处理效果。由图 1可知，在图像覆盖分割之后多获得的图像仍可能存在模糊现象，虽然该方法对噪声不敏感，但噪声的存在还是会对分割结果产生影响。从1.2节对不同周长估算方法的阐述可知，一般方法对模糊图像不适用，而基于图像粒的方法在确定边界厚度时图像中的噪声因素会干扰算法确定的粒度大小，从而导致估算结果误差变大。

针对在对灰度图像像素覆盖分割之后获得的像素覆盖数字化图像可能仍存在模糊和噪声，导致周长估算不精确的现象(见图 1)，结合经典的周长估算方法，提出基于图像森林变换(IFT)的周长估算方法(IDSS、IMLP、IGL)。这类方法利用像素覆盖数字化之后的目标边界处覆盖值与节点间的相互关系，对这些覆盖值所代表的像素点进行优化分类，将模糊边界转换成清晰边界，获得理想的像素覆盖数字化图像。最后结合经典方法，估算出图像中目标的周长。

2.1 图像森林变换

图像森林变换(IFT)^[12]是一种基于图论的图像处理方法，该方法将图像处理问题简化为在图形中求最短路径森林的统一有效方法。

对于一幅给定的二维灰度数字图像I, 图像中的信息可表示为属性对(D_I, I)，其中D_I∈Z²表示图像域，即灰度图像在二维空间的位置集合。向量I表示为I(p)= I₁(p)∈ N，其中p表示图像中的像素，N表示非负数集合，对于每个p∈D_I，I(p)表示图像像素p所在的位置的灰度数值。

在上述定义的基础上，IFT将图像I映射为图结构G=(N, A)。N为图的节点集合且N⊆ D_I, 每个像素点对应一个节点。A为图节点之间的连接关系且A⊆ D_I×D_I，A(p)表示像素对应节点p的连接节点集合。为了利用IFT实现不同的图像处理目标，需要根据图像的局部特性为每对连接点之间的弧分配一个权值w, 通常选择的是定义一个基于灰度强度(像素值)或灰度梯度的代价函数f。在图结构中，一条路径定义为从起始节点按连接关系连接路径上各个节点，直至终止节点的有序序列，而带权路径长度为路径上每两个相邻节点之间路径代价之和。若将路径定义为P =〈p₁, p₂, p₃, …, p_i〉, 则路径长度P可通过式(4) 计算：

$ \left| P \right|={{w}_{{{p}_{1}}{{p}_{_{2}}}}}+{{w}_{{{p}_{1}}{{p}_{_{2}}}}}+\ldots +{{w}_{{{p}_{i-1}}{{p}_{~}}_{i}}} $

(4)

式中w_xy表示节点x和y之间的弧的权值。

IFT的一般过程：进行图像森林变换之前，将像素标签L初始化；在变换过程中，计算每个非种子节点到种子节点的最短路径，根据计算结果更新相应节点的L；最后，在变换结束之后，根据L的值后续处理图像。具体的变换过程根据实际应用需求确定，变换的输出是最短路径的生成森林与由所有节点L值所确定的与输入图像对应的标签字典，若输入图像是二维灰度图像，则标签字典是一个由L的数值填充的二维矩阵。图 2为文献[14]作者展示的一般的算法效果。

2.2 基于IFT的周长估算方法

由1.1节阐述的内容可知，在对灰度图像预处理过程中需要进行像素覆盖数字化。而对于低分辨率与目标模糊的图像，在进行覆盖数字化后，其目标与背景之间可能仍存在一个过渡的像素覆盖值过渡带，为了获得更加精确的周长估算结果，需要将这些覆盖值点重新分类成目标或背景，从而获得清晰的目标边界。根据覆盖像素值，计算边界每个过渡带的点到完全覆盖的点(覆盖值为1) 和完全未覆盖的点(覆盖值为0) 的最短路径长度。若过渡带的覆盖点到完全覆盖点的路径长度小，则将其分类在完全覆盖的点一类，否则分类给完全未覆盖的点一类。若像素覆盖数字化后图像清晰，则图像真实边界应在覆盖值为0和1的像素之间，否则真实边界应在覆盖值为0.5的像素附近。由此可以定义路径代价函数f为

$ f=\left| {{p}_{~}}_{i}-{{p}_{~}}_{j} \right|+\varepsilon $

(5)

式中：p_x表示像素覆盖值，即覆盖数字化之后的新灰度值；ε表示极小的正数，用于确保路径长度是递增序列。为每个点对应添加一个用于记录在变换过程中记录当前最短路径长度的属性c。

为了简化计算，记A为求像素的4邻接像素，则A(p)表示p的4邻接点的集合。根据以上描述，基于IFT的边界优化算法的实现步骤如下：

1) 输入一幅灰度图像I, 将图像I进行简单的像素覆盖数字化。

2) 将所有覆盖值为0和1的点作为种子点集合S的元素。为所有点初始化标签L, 将表示目标节点的标签初始化为0，表示背景节点的标签初始化为1，其他覆盖值点的标签初始化为-1。定义队列数据结构Q, 将S中的点按序加入队列Q, 并将它们的属性c初始化为0，将其余点的权值c属性初始化为+∞。

3) 从队列Q中获取队列头部的节点p，对它的所有邻接节点q∈A(p)按公式(5) 计算它们之间的路径长度w_pq，若w_pq+c_p小于q的权值c_q, 则将w_pq+c_p赋值给c_q，将p的分类标签L_p赋值给q的分类标签L_q，将节点q也加入队列Q。

4) 重复步骤3) 直到队列Q为空队列。根据标签变量L的值重新分配图像灰度值，得到一幅目标边缘清晰的二值图像。

在本文的IFT方法处理过程中，以256个灰度级的图像为例，获得标签集合L。按L中的标签值生成新的图像I′，I′中的像素值按式(6) 计算：

$ {p}'=\left\{ \begin{align} & 255, \ \ \ \ \ \ \ \ {{L}_{p}}=1 \\ & 0, \ \ \ \ \ \ \ \ {{L}_{p}}=0 \\ \end{align} \right. $

(6)

式中：p′表示节点p对应的像素的位置；L_p表示p节点位置的标签值；255和0表示最大灰度和最小灰度。

新的图像为目标边界清晰的二值图像。该方法能消除图像目标边界上的噪声对图像目标周长产生的影响，并且确定了目标边界的近似位置。本节算法的流程如图 3所示。

接下来就结合经典的方法(DSS、MLP、HPGL)，以IFT方法的处理结果作为后续算法输入，对最后的图像进行目标周长的估算。特别地, 以变换后的图片作为输入图像，结合HPGL估算周长需引入的误差矫正参数y_n，根据公式(7)^[8]确定：

$ {{y}_{n}}=\frac{2}{1+{{\left( \sqrt{{{n}^{2}}+1}-n \right)}^{2}}+1} $

(7)

其中n为数字量化的级数, 本文中n=1。

3 结合IFT的目标周长估算

为了验证本文所提出的方法的准确性，对模糊图像和含噪声图像目标周长估算的适应性，本节将该方法用于计算不同合成图像的周长，并将所计算的结果与标准周长、BTHPGL方法^[10]和对应的基于自适应图像粒的ADSS、AMLP、AGL^[11]方法实验所得的结果进行比较。

实验平台为Intel(R) Core(TM) i5-3470M CPU @3.20 GHz RAM 8GB的PC机，64位Win7系统和MATLAB 2012b。实验中的图像是通过选择适当的函数合成的简单目标图像，通过积分计算目标的标准周长。实验中还进行以下操作：对图像目标边界进行人工加厚模糊化，选取不同的模糊核对图像模糊，然后进行像素覆盖数字化，用于验证算法对模糊边界的适应性。对图像加入噪声，用于验证算法对含噪声图像的适应性。实验都是以单位像素宽为单位。

3.1 模糊图像周长实验

实验选取圆、弦月形、阴阳图形的一半、正方形旋转45 °与正方形旋转22.5°的5种形状合成的正方形图像，其中圆的直径为图像边长的4/5，正方形边长为图像边长的3/5。对每个图形处理并分6类，名称后缀为1的表示理想状态的图像，后缀为2的表示以散焦半径为2的散焦模糊的图像，后缀为3的表示图像中目标物体按与水平方向成15°角以匀速直线运动方式产生5个像素宽模糊的图像；后缀为4~6的分别表示以图 4、5所示的模糊核对原图像进行模糊的图像。由于DSS、MLP和HPGL方法对模糊图像不适用，所以这里选取BTHPGL、ADSS、AMLP和AGL方法作为对比方法，本文所提出的改进方法为IDSS、IMLP和IGL。测试图像和模糊核数据采用文献[15-16]中的数据集，表 1~5是部分实验结果，选择图像规格为512×512个像素, 统计的是与标准周长相对误差的绝对值。

3.2 模糊图像实验结果分析

图 6是实验结果的平均相对误差的折线图，分析图 6和上述表格中的实验数据不难发现，对于不同种类的模糊图像目标周长估算使用ADSS、AMLP、AGL等方法产生的误差很不稳定，这是由于基于图像粒的方法适合边界像素分布均匀的情况且经过粒度处理的图像本身也会增加估算误差。而基于IFT的方法总体上对周长的估算精确度更高：IDSS、IMLP方法相比ADSS和AMLP方法具有更好的稳定性，并且大部分情况下精确度更高，而IGL也比BTHPGL和AGL方法估算更加稳定，并且在对曲线的估计上精确度也更高。

当图像完全符合像素覆盖数字化定义(本文中这类图像名称后缀为1) 时，基于IFT的方法会将图像先转换为二值图像，这会增加周长估算的误差，但从表 1~5中数据可以看出，误差在可承受的范围之内并且现实图像经过像素覆盖数字化后大部分情况下很难达到理想的效果。因此, 可以认为基于IFT改进的周长估算方法优于基于图像粒的方法。

3.3 复杂目标含噪图像周长实验

在平面坐标系中，选取直线y=820-2x与y=$\frac{10}{7}$x+160，与曲线(x-300)²+(y-300)²=200²，三角函数$y=5\cos \frac{\text{ }\!\!\pi\!\!\text{ x}}{50}+{{\sin }^{2}}\frac{\text{ }\!\!\pi\!\!\text{ x}}{50}+150$与$y=10\cos \frac{\text{ }\!\!\pi\!\!\text{ x}}{50}+10\sin \frac{\text{ }\!\!\pi\!\!\text{ x}}{50}+450$围成的边界图像，用MATLAB软件工具合成图像，然后增加5组均值为0方差(图像归一化后)为0.01~0.05的随机噪声的图像。实验合成的曲线标准周长根据积分计算结果为1 269。图 7列出了基于散焦模糊的无噪声和含噪声均值为0方差为0.05的两幅图像。

实验方法与3.1节中的实验方法相同，图片编号Pic1、Pic2、Pic3、Pic4、Pic5和Pic6对应无噪声、方差为0.01、0.02、0.03、0.04和0.05的图像。表 6为含噪声的模糊图像应用相应周长估算方法的实验结果。

3.4 含噪图像实验结果分析

从表 6的实验数据可知，基于图像粒的方法精确度明显下降，分析其原理，这类方法会因噪声的干扰，导致所求得的边界宽度小于实际宽度，从而导致重复计算问题，求得的结果比实际大很多。基于IFT改进的方法对噪声表现出了良好的适应性与稳定性，并且具有较高的精确度。从表 7的实验数据可以看出，基于图像粒的方法需要根据其概率模型计算出估计的边界宽度，从而确定粒度尺寸，计算过程耗时较多。然而基于IFT的方法虽然过程较复杂，但主要是需要辅助的数据结构，其运行速度相对很快。而BTHPGL、ADSS、AMLP等方法需要追踪边界，图像中的噪声会干扰边界的获取，导致估算失败，因此未在表中列出。通过这些实验可得出结论，本文算法相比现有的算法具有更好的抗噪性和更快的运行时间。

4 结束语

本文基于IFT理论，结合经典估算方法，提出具有高精确度和广泛适应性的灰度目标边界周长估算方法。该方法继承了经典方法的高精确度，并且对模糊图像具有适应性，解决了图像在模糊情况下的周长难以估算问题。因此本文方法能达到良好的精确度、适应性和稳定性，在实际图像特征提取中能得到广泛应用。另外, 本文所设计的方法考虑的情况并不全面，算法精确度仍有提升空间。而在更复杂的情况下估算出更精确的灰度图像目标周长，这将是今后的研究方向。