土壤孔隙结构三维重建软件系统的设计与实现

吕菲

土壤结构尤其是土壤孔隙结构特征不仅决定了许多关键的土壤物理性质(水分持留、传导能力等)，而且也控制着发生在土壤中的诸多重要生态过程(如土壤水分运动、养分迁移等)。快速获取原状孔隙结构的定量信息是从根本上认识土壤结构与功能间关系及有关过程内在机制的必要前提，也是土壤学工作者长期关注的热点和难点问题之一。

土壤微形态学研究中通常采用的切片方法不仅制备过程费时费力，而且根据二维的切片信息推断三维的土壤结构容易出现偏差甚至错误。近20多年来，CT(Computerized Tomography)、核磁共振(NMR)等无损、显微探测技术的长足进步使得高精度数字图像的快速获取成为可能，其在医学、工程技术等领域的成功应用、推广为我们提供了可供借鉴的新思路。在土壤学研究领域，国内近年也逐渐有学者开展了此方面的研究，如冯杰和郝振纯(2002年)、贺秀斌等(2005年)、吕菲等(2007年)等。但就目前国内的研究现状来看，由于缺少通用的或有针对性的数字图像分析、三维重建及可视化软件工具，已经大大制约了研究工作的深入、全面开展。针对这一问题，本文尝试开发基于计算机图形学的、集图像数据管理、数字图像分析处理及三维可视化为一体的土壤孔隙结构分析与三维重建软件系统，以期为当前有关领域的研究提供一个功能强大、简单易用的操作平台。

1 总体设计思路

由于实际土壤在孔隙形态学上的复杂性，为了便于三维网络模型的建立与求解，孔隙体和孔喉通常被简化成简单的几何形状:孔隙体多以球形或立方体表示，孔喉则以圆柱形或其他具有代表性的非环形管道表示。虽然这样形成的网络模型是三维的，但是在视觉上与实际土壤内部结构还有较大差距。土壤孔隙结构三维重建系统是将一系列二维切片信息合成为三维立体影像的软件，可以直接显示出一个生动的虚拟实体，更直观、更准确的观察样本土壤的内部结构，有助于了解大孔隙对土壤有效水力学性质的影响。

软件开发环境包括:采用Delphi7.0作为软件包的开发工具，OpenGL作为三维图形显示的开发工具。硬件环境:Intel奔腾4 CPU 2.4 GHz，512 M内存，GeForce4显卡。土壤孔隙结构三维重建系统主要由主模块、读取图片模块、图片观察模块、三维重构模块以及数据存储模块组成。主模块是系统的主体，负责完成系统的所有操作和各个模块的激活;读取图片模块将由X射线断层摄影术(CT)、核磁共振(NMR)、扫描电镜等无损技术获得的二维连续切片灰度图像数据读取到系统中;在图片观察模块中可以对这些二维切片进行缩放观察、开运算、二值化处理等;系统在三维重构模块中将处理后的二维连续切片按照设定的灰度值范围以及阈值进行三维重建，重建出的三维立体影像可以实现旋转、拖动、平移、光照控制等功能;整个过程中用到的数据都统一在数据储存模块中进行存放或读取。整个系统的操作流程见图1。

2 三维重建算法

三维重构是本软件系统的核心功能，其算法设计将直接影响系统的运行效率、计算精度和可视化效果。用于三维重建的算法和技术有多种，如移动立方体算法(Marching Cubes，简称为MC算法)、剖分立方体法(Dividing Cubes)、射线投射法(Ray Casting)、剪切—曲变法(Shear-Wrap)等(张季和王宜杰，2006年)。其中，由Lorensen和Cline(1987年)提出的MC算法是发展较为成熟且得到最广泛应用的一种表面重建算法。其基本思路是通过在三维网格顶点上的采样构建出平面的组合来近似表达空间曲面，即将空间分割为一系列的立方体，然后将每一立方体在三维数据场中移动，并对立方体的顶点与等值曲面进行判断比较，用合适数量的多边形取代立方体，所有产生的多边形将形成许多个等值表面，这些表面就是对输入数据的近似表达(蒋先刚，2006年)。虽然MC算法的内存开销较大，但其在编程上容易实现，且得到的三维图像品质较高。本文中采用的MC算法的程序框图如图2所示。

算法首先建立单个立方体与等值曲面相交状态的索引值表。对于一个已知的等值面的值，立方体每个定点处的值或大于或小于或等于该值。若顶点数据值大于或等于等值面的值，则该顶点位于等值面之上，赋值为1;否则该顶点位于等值面之下，赋值为0。等值面只与那些两端点取不同值的边相交。因为立方体有8个顶点，每个顶点有两种可能的状态值，故立方体与等值面相交的情况最多有256种(28)。将各顶点的状态值拼成一个字节进行编码即可构成一张索引值表(见图3)，对应256种可能的相交情况。如索引值为二进制数00001000，指出等值面与立方体的e3，e4，e11三条边相交。顶点v4高于等值面，其余顶点都低于等值面。

建立相交状态的索引值表后对每个立方体内的等值面进行三角化。通过以下三个条件:1)绕轴旋转一定角度;2)图形相对轴面镜面对称;3)顶点的位置状态取反，将256种位置状态减少到不同构的15种(见图4)。这样就更加容易建立起三角面集合近似表达物体表面。

在三角化过程中物体表面与立方体边的交点可以用线性插值的方法得到，精确的插值点取决于插值过程中用到的量值(对应二维图像中的阈值)。例如，P1和P2是邻接边的两个端点，V1和V2是这两个点对应的采样值(对应二维图像中的灰度值)，物体表面与边的交点P的坐标可由以下公式计算得到:

用一定规格的立方体分割空间，确定网格点的坐标及其对应的采样值(对应二维图像中的灰度值)。按顺序选择一个立方体，根据其8个顶点的不同标识获得索引值并查表获得当前立方体各相交边所在，依次在各交边上找到交点坐标，根据索引值和交点坐标形成当前立方体内三角形片的构成，选择下一个立方体直到所有立方体都找到与之对应的三角形片构成。这样构成的三角形片整体就是三维表面的近似表示。

从二维图像提取出三维空间信息后，需要采用三维图形开发工具将模型数据进行可视化。本文中的重建算法和三维可视化通过Borland Delphi来编程实现，其中使用了与Delphi7.0兼容的OpenGL 1.1。

实现的系统功能主要包括:1)绘制三维物体:通过读取三角形数据链表绘制需要的三维物体;2)查看物体:建立视点，指定观察角度、方向和观察的范围，实现对象的旋转、缩放;3)指定颜色模式和光照:对绘制的三维物体加以着色、添加材质、光照，使其逼真于真实环境下的物体特性。

3 实例分析

供试土壤取自黄淮海冲积平原的中国科学院河南封丘生态实验站内，土壤质地为壤质粘土(国际制)。按常规分析方法测定土壤容重为1.54 g/cm3，孔隙度为0.42。

CT扫描样本用高20 cm、内径10 cm的PVC管采集，采样深度35 cm～70 cm。样本运至中国科学院兰州寒区旱区研究所进行测定，采用Siemens公司生产的SOMATOM-PLUS螺旋式CT机，综合图像处理系统为FLS-7。最大扫描分辨率为62.5μm，相邻两层切片间距为1mm。选取样本中形状保留完好的2 cm厚土柱进行扫描，共获得连续切片数字图像20张。

首先打开土壤三维重构系统并读取切片图像集，进入二维切片图像操作界面，选取原状土柱CT扫描获得的8张连续图片，在界面中可设置三维重构区域的灰度值边界以及二值化阈值。阈值的选择可根据实际土壤孔隙度按照单阈值法确定。

生成的土壤孔隙结构三维重建图像如图5所示，网格体数目设置为128×128×7，图5中黑色实体部分为土壤孔隙，图5中显示出较大孔隙(＞750μm)的位置也与试验土样相一致。在此界面中可以对三维重构模型进行旋转、平移或选取特定部分显示等功能，更好的了解土壤的内部微观结构。

4 结语

土壤孔隙结构三维重建系统可以根据CT，NMR等获得的二维序列灰度图像，通过计算机图像处理与可视化技术，构建出三维土壤孔隙的虚拟实体，逼真、直观、准确的显示出土壤微观结构，为土壤微观模型的构建提供了一个可视化的软件平台。

在今后的研究中，需在以下几个方面对该系统作进一步完善和改进:1)土壤孔隙结构三维重构系统仅生成了土壤孔隙的虚拟实体，缺乏对二维图片以及三维重构模型分析计算的功能，不能得到实际土壤有效水力学性质方面的结果。2)对每个立方体进行三角形构造过程中，并没有考虑其周围立方体的情况，造成部分计算的冗余，应该进一步改进重建算法提高重建速度与精度。3)数据存储功能尚不完善，一次重建后得到的图形数据应当以标准格式储存起来，以便下次直接调用。储存格式可以与其他三维绘图工具软件(AutoCAD，3DS等)相一致，方便相互调用。

［1］ W.E.Lorensen，H.E.Cline.Marching Cubes：A High Resolution 3D Surface Construction Algorithm［J］.Computer Graphics，1987，21(4)：163-169.

［2］ 冯 杰，郝振纯.CT扫描确定土壤大孔隙分布［J］.水科学进展，2002，13(5)：611-617.

［3］ 吕 菲，刘建立，何 娟.利用CT数字图像和网络模型预测近饱和土壤水力学性质［J］.农业工程学报，2008，24(5)：10-14.

［4］ 曾 筝，董芳华，陈 晓，等.利用MATLAB实现CT断层图像的三维重建［J］.CT理论与应用研究，2004，13(2)：24-29.

［5］ 蒋先刚.基于Delphi的数字图像处理工程软件设计［M］.北京：中国水利水电出版社，2006：291.

［6］ 张 季，王宜杰.医学图像三维重建方法的比较研究［J］.医学信息，2006，19(5)：948-950.