从物体三维图像中获取四面体网格的方法

摘要

本发明涉及一种图像处理技术领域的从物体三维图像中获取四面体网格的方法。本发明直接在运行的计算机上进行，首先读入物体三维图像，然后依次通过图像前处理模块、采样模块以及网格重建模块处理，得到四面体网格。其中图像前处理模块对输入的二值图像反复进行二值腐蚀操作获得图像的距离变换结果，在距离变换的结果上计算出密度函数，并将其输出到采样模块。采样模块根据密度函数对图像进行采样，得到所需要的空间采样点，并将其输出到网格重建模块，网格重建模块由上述采样点构建最终的Delaunay网格。本发明易于实现，在网格质量上相比于现有方法，拥有良好的表现。

说明书

技术领域

本发明涉及的是一种图像处理技术领域的方法，具体是一种从物体三维图像中获取四面体网格的方法。

背景技术

随着科学技术的高速发展，高科技医疗设备的不断涌现为医疗的现代化提供了越来越多的帮助。在众多的计算机辅助医疗系统中，获取建立系统分析数据所需要的四面体网格，是系统必须具备的功能模块。该模块需要运用计算机图形学、计算几何、数字图像处理相关专业知识，将以像素为基本表示单位的图像，转换生成为以四面体为基本单位的网格数据。在计算机辅助医疗系统中，需要将真实病人的图像数据(如CT、MRI等)，转换为计算机图形学范畴内的四面体网格，再进行后继的分析处理。生成的网格的质量将严重影响后继处理的效果及正确性。因此，国际上众多研究机构在该技术上倾注了大量的资源和精力。

当前医学数据四面体网格重建，一般分为两个阶段。第一阶段为数据的表面网格重建，即抽取出图像的边缘信息，而后通过一系列处理，将该边缘信息转换为三角形或其他形式的表面网格。代表方法有Marching Cube、α-Shape、CRUST等。第二阶段为四面体网格重建阶段，即在前一阶段的基础上，用四面体网格填充表面网格，得到最终的四面体网格。代表方法很多，一般说来，可以分为基于栅格的方法、基于Delaunay网格的方法、前进波方法等。一般的医学仿真软件，通过以上两个阶段的处理，得到最终所需要的四面体网格。

然而，上述处理方法存在以下几个问题：首先，两个阶段的衔接存在问题。一般情况下，第二阶段对第一阶段的处理结果有一定的要求，如必须生成封闭的表面网格数据、质量要达到某些指标等等。在实际应用场景中，这些要求第一阶段的技术未必能够百分之一百的达到，可能导致最终无法正确地生成网格或者生成低质量的网格。其次，绝大多数的处理技术，输出的网格的正确性及质量，严重依赖于用户手工设定的参数值。这样一来，在实际应用中，用户往往需要凭经验、多次试验才能得到所需要的一定质量的网格，带来了相当的麻烦。再者，上述方法，基本上来源于计算机辅助设计制造领域，并未针对医学图像数据的特点，进行相应的优化和处理，导致难以生成高质量的网格。最后，由于实现上述技术的复杂性较高，在应用到实际的工程领域时，亦会存在不小的困难。

经对现有技术的文献检索，N.Amenta等在Proceedings of the 25th annualconference on Computer graphics and interactive techniques(98年出版的第25届计算机图形与交互技术年会论文集)，1998.pages 415-421，上发表的“A new Voronoi-based surface reconstruction algorithm”(《一种基于Voronoi图的表面重建方法》)中提出了一种利用插入某些特殊的空间点，通过两次构建Voronoi图的方式建立表面网格数据。该方法存在的问题是，最终输出的表面三角形网格，不能保证为流型网格，且可能出现空洞，不能保证为封闭网格。这对于后继的造型技术，如四面体划分等，存在很大的障碍。再如，Qiang Du等人在International Journal for Numerical Methods in Engineering杂志Volume 56 Issue 9，Pages 1355-1373中发表的“Tetrahedral mesh generationand optimization based on Centroidal Voronoi tessellations”(基于质心Voronoi图的四面体网格划分及优化)。该文主要针对从表面网格到四面体网格的划分与网格质量优化，其主要问题在于该方法单纯的追求四面体的质量，虽然对四面体的质量进行了局部优化，但是，在全局上并未能保证四面体的采样率与物体几何形态相关。此外，该文仅给出了简单数据的测试结果，并未给出对于复杂形态物的建模效果，其在复杂场景下的实用性值得怀疑。

发明内容

本发明的目的在于克服现有方法中的不足，提供一种从物体三维图像中获取四面体网格的方法，直接从图像数据中进行网格重建，同时，针对医学图像数据的特点，考虑复杂图像的优化，生成Delaunay网格，方便后继处理，保证网格质量。

本发明是通过以下技术方案实现的。本发明直接在运行的计算机上进行，首先读入物体三维图像，然后依次通过图像前处理模块、采样模块以及网格重建模块处理，得到物体三维图像四面体网格。其中图像前处理模块对输入的二值图像反复进行二值腐蚀操作获得图像的距离变换结果，在距离变换的结果上计算出密度函数，并将其输出到采样模块。采样模块根据密度函数对图像进行采样，得到所需要的空间采样点，并将其输出到网格重建模块，网格重建模块由上述采样点构建最终的Delaunay网格。该方案稳定可靠，自动化程度高，且实现方便，非常容易在实际的项目中加以应用。

其中：

所述的图像距离变换，指的是求取图像每个像素点到物体边界的最短距离。该距离可以直观的理解为，以相应像素位置为圆心，能够在物体内部放置的最大的球体的半径的大小。距离变换所得到的数值，反映了在不同位置，物体的几何形态特征。半径越大，说明物体在该像素位置的几何形态变化越平缓，反之，说明该像素位置的物体几何形态变化越剧烈。

所述的质心Voronoi图，是一种特殊的Voronoi图。输入一组采样点，每个采样点可以得到一个对应的Voronoi单元，单元内的所有空间点，到该采样点的距离，小于到其他采样点的距离。所有的Voronoi单元，形成了对空间的一个划分，该划分被成为对空间的Voronoi图。质心Voronoi图是Voronoi图的一个特例，即，给定空间内的密度函数，用于划分Voronoi图的所有采样点，恰好位于其对应的Voronoi单元的质心位置。

所述的Delaunay网格，是一种特殊的网格结构。三维Delaunay网格由四面体组成。每个四面体的外接球内，不包含除该四面体的四个顶点之外的其他采样点。在数学上，Delaunay网格与Voronoi图呈对偶关系。用于生成Voronoi图的采样点，即为Delaunay网格中四面体的顶点。当两个Voronoi单元相接壤，则连接该区域对应的采样点，将得到对应的Delaunay网格内的一条边。

本发明包括如下步骤：

第一步，将物体的三维图像数据，以二值的形式读入系统；

第二步，图像前处理模块接收输入的二值图像，并对其反复进行二值腐蚀操作，将第k次图像二值腐蚀中被腐蚀的像素置为k，直至所有有值的像素点都被赋值为止，得到图像距离变换的结果，k为自然数；

第三步，图像前处理模块在获得图像距离变换的结果基础上，遍历所有像素点，获得每个像素点的密度值，并将其输出到采样模块；

第四步，采样模块接收图像前处理模块输出的像素点的密度值，在图像内随机选取初始采样点，然后进行迭代操作；

第五步，采样模块每次迭代获得当前采样点集的Voronoi图，从而得到图像中每个Voronoi单元的质心及质量，将获得的Voronoi单元中对应的采样点位移到该质心位置；

第六步，检测Voronoi单元的质量，根据其质量控制在图像不同区域的采样率，如果质量小于用户设定的参数值，则删除该采样点；如果质量大于用户设定的参数值，则在该采样点附近随机插入新的样点，每次迭代得到的采样点趋近于Voronoi单元的质心位置，采样点的变化将会逐步减小，直至最终几乎没有发生位移。如果插入或者删除的采样点数量小于给定的某个值，且采样点的位移和小于指定值，则采样模块结束处理，将处理的结果输出到网格重建模块；

第七步，网格重建模块接收采样模块输出的采样点，由边界采样点构造对应的Voronoi图，以此得到所有的Voronoi单元，并找到每个Voronoi单元的两个极点，将找到的所有的位于物体外部的极点，记为集合P，然后根据所有的采样点以及极点构建Delaunay网格；

第八步，检查构造出来的Delaunay网格内的所有四面体，如果四面体至少有一个顶点在集合P中，则该四面体在物体的外部，否则就在物体内部，删除外部四面体，保留内部四面体，得到最终需要的四面体网格数据并输出。

所述极点，即该边界采样点所在的Voronoi单元中(呈细长状)，在长轴方向上，距离最远的两个顶点。

所述每个像素点的密度值通过物体密度函数得到，当当前像素点位于物体外部时，密度函数值为0。否则，根据下述公式进行计算：

$\rho (x,y,z)=\frac{\tan \left[\alpha (1.0-\frac{d(x,y,z)-1}{\max \{d(x,y,z)-1\}})\right]}{\tan \alpha }(1.0\le \alpha <\frac{\Pi }{2})$

其中，ρ(x，y，z)为密度函数，d(x，y，z)为前述的距离变换函数，α为用户输入的参数。

本发明利用本发明进行物体三维图像数据的四面体网格重建，具备如下优点：(1)在工程应用中便于操作和实现。在实际应用中，可以在非常短的时间内对该技术进行实现并集成至系统中。(2)网格质量良好，非常适合于进行几何造型及其后继处理。该方法生成的网格，是标准的Delaunay网格，具备Delaunay网格应有的所有的数学特性。且，网格的分布合理均匀，在物体中心附近，网格体积较大，越接近边界，网格体积越小。(3)一步构建网格，省去了用户需要多步进行构建、反复调节参数的麻烦，避免了原有表面重建、四面体网格划分两步走的数据无法交接的问题。(4)该方法针对医学图像数据的特点，非常适合于复杂形态的物体的网格获取，尤其是人体各个组织器官。

附图说明

图1本发明方法流程图

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

本实施例在CPU为Pentuim M 1.8GHz，显卡为NVIDIA Geforce 9800GTS，内存为2.0GB的计算机中实现，编程语言为C++。

本实施例的实施流程如图1所示。

测试用两个图像具体参数情况如表1所示：

表1 实施例中各种参数情况

 

情况图像大小四面体个数顶点数图像A168x142x22113702525905图像B168x142x22138434869050

第一步，将图像A和图像B的上述三维图像数据，以二值的形式读入系统；

第二步，图像前处理模块接收输入的二值图像，并对其反复进行二值腐蚀操作，将第k次图像二值腐蚀中被腐蚀的像素置为k，直至所有有值的像素点都被赋值为止，得到图像距离变换的结果，k为自然数；

第三步，图像前处理模块在获得图像距离变换的结果基础上，遍历所有像素点，获得每个像素点的密度值，并将其输出到采样模块；

第四步，采样模块接收图像前处理模块输出的像素点的密度值，在图像内随机选取初始采样点，然后进行迭代操作；

第五步，采样模块每次迭代获得当前采样点集的Voronoi图，从而得到图像中每个Voronoi单元的质心及质量，将获得的Voronoi单元中对应的采样点位移到该质心位置；

第六步，检测Voronoi单元的质量，根据其质量控制在图像不同区域的采样率，如果质量小于用户设定的参数值，则删除该采样点；如果质量大于用户设定的参数值，则在该采样点附近随机插入新的样点，每次迭代得到的采样点趋近于Voronoi单元的质心位置，采样点的变化将会逐步减小，直至最终几乎没有发生位移。如果插入或者删除的采样点数量小于给定的某个值，且采样点的位移和小于指定值，则采样模块结束处理，将处理的结果输出到网格重建模块；

第七步，网格重建模块接收采样模块输出的采样点，由边界采样点构造对应的Voronoi图，以此得到所有的Voronoi单元，并找到每个Voronoi单元的两个极点，将找到的所有的物体位于外部的极点，记为集合P，然后根据所有的采样点以及极点构建Delaunay网格；

第八步，检查构造出来的Delaunay网格内的所有四面体，如果四面体至少有一个顶点在集合P中，则该四面体在物体的外部，否则就在物体内部，删除外部四面体，保留内部四面体，得到最终需要的四面体网格数据并输出。

所述极点，即该边界采样点所在的Voronoi单元中(呈细长状)，在长轴方向上，距离最远的两个顶点。

所述每个像素点的密度值通过物体密度函数得到，当当前像素点位于物体外部时，密度函数值为0。否则，根据下述公式进行计算：

$\rho (x,y,z)=\frac{\tan \left[\alpha (1.0-\frac{d(x,y,z)-1}{\max \{d(x,y,z)-1\}})\right]}{\tan \alpha }(1.0\le \alpha <\frac{\Pi }{2})$

其中，ρ(x，y，z)为密度函数，d(x，y，z)为前述的距离变换函数，α为用户输入的参数。

图像A的实例与已有方法的比较。方法一为栅格化网格方法，即先将物体所在包围盒划分为相同大小的立方体，而后再将每个小立方体划分为四面体。方法二使用原有的两步走的策略。使用TightCocone方法进行表面重建后，再使用NETGEN方法进行四面体网格划分。

表2中，计算、统计了以下几个质量衡量指标。这些指标，数值越高，表明网格质量越高。

参数1：，衡量基于该网格进行线性插值后的函数，与真实函数的误差限大小。其中，V为四面体的体积，r_mc为包围该四面体的最小包围球的半径大小。

参数2：，衡量及于该网格进行线性插值后的函数的一阶导，与真实函数的一阶导之间的误差限。其中A_m为四面体的三角形面的面积，1_ij为四面体的边长。

参数3：，衡量该网格对有限元中刚度矩阵的条件数的影响。

表2 图像A实例实施结果及与传统方法的比较情况

从表2中可以知道，不论从网格的何种指标上来看，本实施例所得的网格，都具备更加合理优异的全局分布。即，高质量的网格所占的比例，远远高于其他方法。如，参数3中，本实施例高质量网格的比例为37.033％，大于其他两方法的23.732％和14.230％。这表明，用此方法得到的网格，在后继的其他模块(如数值计算)的相关任务中，将比其他方法更具优势。