摘  要： 在岩心三维模型中，中轴是描述孔隙结构特征的一种重要表示方法。针对现有的拓扑细化和距离变换方法存在的中轴提取不准确和不连续的问题，提出了一种使用快速行进水平集方法进行距离变换的三维中轴提取算法。对比实验结果表明，该算法中提取的中轴在准确性和连续性上有着较好的保证，实际应用中效果良好。
　　关键词： 中轴提取；水平集；快速行进方法；距离场
0 引言
　　三维模型分析技术在分析岩心孔隙结构和其相应的统计特性中有着重要的应用，中轴骨架是孔隙三维模型分析中的一项重要参数。中轴骨架表现出了三维模型的拓扑结构，使用中轴描述三维模型不仅能很好地表示三维模型的结构信息，而且还能提高内存的使用率和数据压缩率。在岩心三维模型分析中，将中轴应用于模型的迂曲度计算，有着重要的实际意义和工程应用价值。
　　国内外已有许多学者对三维中轴提取进行了研究，目前的三维图像骨架化方法主要有拓扑细化方法、基于距离变换的方法、基于V图的方法等。
　　拓扑细化算法[1-2]能较好地保证拓扑结构，但对边界噪声较为敏感，不能保证中轴位置的准确性。距离变换算法[3-5]能够提取出不规则物体的骨架并且具有平移、旋转、缩放不变性，但是很难保证模型的连通性。V图方法[6-7]是利用V图思想以中轴面或中轴线为基础得到模型的骨架，一般用于生成多尺度的骨架。此外还有将距离变换和迭代并行细化相结合的方法[8]，得出的骨架具有良好的连通性和拓扑等价性，但对于三维体数据需要多次迭代，计算量较大。
　　本文针对岩心三维孔隙模型进行迂曲度计算等需求，提出了一种使用水平集进行距离变换的方法求解三维结构的骨架，使用水平集的距离变换算法与一般的距离变换方法相比具有更好的稳定性和拓扑无关性，解决了距离变换存在的连通性问题。该方法提取出的中轴骨架具有很好的连通性和拓扑等价性。
1 水平集方法
　　水平集方法是1988年由OSHER S和SETHIAN J首次提出的[9]，用于解决遵循热力学方程下火苗外形的变化过程。其基本思想是将曲面（用x表示）看成高一维空间中某一函数（称为水平集函数）的零水平集，同时曲面的演化也扩充到高一维的空间中。曲线的演化转化成一个纯粹的求偏微分方程数值解问题。在任何时间，波前的位置可以由水平集函数的零水平集来确定。将水平集函数按照它所满足的发展方程进行演化或迭代，由于水平集函数不断进行演化，所以对应的零水平集也在不断变化，当水平集演化趋于平稳时，演化停止，得到界面形状。
　　假设以速度F沿x（t）法线的方向进行移动，若波前的粒子在路径上，则粒子的水平集值必定为0：

　　其全微分方程为：

　　由于曲线沿法线方向移动，F可以表示为：

　　结合以上各式可以得出：

　　这就是参考文献[9]给出的水平集公式。
2 快速行进方法
　　考虑波前以速度F沿一个方向运动，根据速度函数的符号，波前单调递增或递减传播，计算每个点的到达时间。运动公式可以表示为：

　　其中，F是速度函数，T是到达时间。
　　在三维空间中，求解上述等式需要计算出每个点（x，y，z）的到达时间。三维空间中梯度的离散化形式为：

　　为了求解这个等式，SETHIAN J[10]提出了一种简单但是低精度的公式：

　　式中前向差分和后向差分的表现形式如下：

　　快速行进算法的更新过程如下：
　　（1）从初始点开始推进算法。将初始点加入到已知点集合，标记为Known。
　　（2）找出Known点的邻域，对邻域点计算到达时间，将求解点标记为trial。
　　（3）选取trial中到达时间最小的点，标记为Known。
　　（4）查找所有Known点邻域内所有的未标记点，对其计算到达时间，将求解点标记为trial。
　　（5）转至步骤（3），直到遍历所有点，退出循环，算法结束。
3 使用水平集方法进行距离变换求解骨架
　　在本文算法中，首先提取出三维模型中所有的连通目标，独立处理每个连通目标。对每个连通目标进行距离变换，求解距离场D（x，y，z）。距离场是目标中所有体素点到边界距离的最小值，可以通过计算欧氏距离得到。选取距离场中具有最大距离的体素点作为全局最大距离点。根据距离场计算速度函数：

　　选取全局最大距离点作为目标的中心点，将目标的边界体素点作为水平集曲线，使用快速行进方法进行演化，演化速度由速度函数决定。在每个体素更新的过程中，快速行进算法计算出了当前体素点到起点的到达时间。
　　对于各向同性的快速行进方法，沿着波前的垂直方向波传播得最快。由于每个中轴点梯度的方向始终垂直于波前，因此中轴体素点具有最快的传播速度，可以根据这一特性来提取出中轴体素点，三维模型的中轴即可通过选定的初始体素沿着到达时间的梯度方向回溯求得。
　　三维模型单个孔隙目标中轴提取流程如下：
　　（1）求解三维模型的距离场，根据距离场求解出速度函数和全局最大距离点。
　　（2）提取三维模型中所有独立的连通目标，对每个目标，将全局最大距离点加入点源集合，利用速度图像，用快速行进方法求解水平集，得到到达时间，将距离点源集合到达时间最大的点作为回溯起始点。
　　（3）从回溯起始点开始，沿着到达时间的梯度方向开始回溯，一直回溯到点源集合为止。回溯完成后得到一条由最远点到点源集合的一条分支。
　　（4）判断分支长度。若分支长度大于距离场中的最大值，则认为这条分支是中轴的一部分，把分支加入中轴中，并且将该分支中的所有点加入点源集合；若分支长度小于距离场中的最大值，则剔除掉该分支，并且该目标的中轴提取完成。
4 实验结果分析
　　本文以下的测试使用图1所示的岩心三维模型，针对拓扑细化方法和本文提出的基于快速行进水平集方法两种不同的中轴提取算法进行了对比。两种方法提取中的中轴结果如图2所示。

　　从图2可以看出，拓扑细化方法保持了中轴的连通性，但是逐步剔除边界点时，为了保证单像素性，产生了很多折线。使用这种中轴计算三维孔隙图像的迂曲度时会产生较大的误差。从图3可以看出，本文使用的算法提取出的中轴较为平滑，在保持连通性的同时，很好地解决了折线的问题。图4为本文提取出的中轴和三维模型的叠加显示，可以看出本文算法提取出的中轴在连通性和正确性上得到了保证。

5 结论
　　本文使用快速行进的水平集方法实现了岩心三维孔隙模型的中轴提取，通过对比实验，可以看出本文提出的算法有着明显的优势，提取出的骨架有着很好的准确性和连通性，有实际意义和应用价值。
　　参考文献
　　[1] 王广垒，张维忠，宋明玉，等.基于数学形态学的鞋楦特征曲线骨架的提取方法[J].青岛大学学报（自然科学版），2012，25（2）：43-46.
　　[2] 徐莹.基于数学形态学的图像骨架提取和复原的改进算法[J].成都信息工程学院学报，2009，24（3）：259-263.
　　[3] 张国栋，韩佳池.基于模糊距离变换的骨架剪枝算法[J].沈阳航空航天大学学报，2012，29（1）：64-69.
　　[4] GAGVANI N， KENCHAMMANA H D， SILVER D. Volume animation using the skeleton tree[C]. Proceedings of IEEE Volume Visualization， 1998：47-53.
　　[5] DEY T K， SUN J. Defining and computing curve-skeletons with medial geodesic function[C]. Proceedings of the fourth Eurographics Symposium on Geometry processing， AirelaVille， Switzerland， Eurographics Association， 2006： 143-152.
　　[6] 刘辉，秦茂玲，徐海峰.基于Reeb图的三维网格模型骨架提取算法[J].信息技术与信息化，2012，24（5）：1672-9528.
　　[7] 吴艳花.三维模型骨架提取算法及其在检索中的应用[D].广州：中山大学，2013.
　　[8] 滕奇志，康瑕，唐棠，等.基于升序复核的并行三维图像骨架化算法[J].光学精密工程，2009，17（10）：2528-2534.
　　[9] OSHER S， SETHIAN J A. Fronts propagating with curvature-dependent speed： algorithms based on Hamilton-Jacobi formulations[J]. Journal of Computational Physics， 1988，79：12-49.
　　[10] SETHIAN J A. Level sets methods and fast marching methods（2nd edition）[M].Cambridge University Press，1999.