摘要：为实现植物根系三维构型参数的原位、无损、准确、快速和自动化测量，提出一种基于X射线计算机断层摄影术（X-ray computed tomography，简称XCT）序列图像的植物根系三维矢量模型构建方法。首先，根据根系的形态特征设计了根结点、根分枝、根系统3个基本结构；其次，对植物根系三维体数据（XCT序列图像）重切以提取根截面，并根据根截面质心和面积构建初始根结点，继而对根结点分组以构建初始根分枝，再通过重构根结点、重获根结点面积和根分枝关系判定构建完整的根系统（矢量模型）；再次，基于矢量模型设计根数、根长度、根体积、根表面积以及根夹角的计算方法；最后，在Windows平台上利用图形用户界面开发工具Qt、可视化工具包VTK和医学图像分割与配准工具包ITK实现上述结构和算法，并将上述参数的计算结果与手工测量值进行对比，验证所提出方法的可行性与准确性。
关键词：植物根系；XCT；矢量模型；三维构型
中图分类号： TP391文献标志码： A
文章编号：1002-1302（2017）14-0179-05
根系是植物赖以吸收養分、水分的重要器官，其构型描述了根系在土壤空间中的分布性状，对植物的吸收能力有着决定性的影响。根系构型的检测与分析方法一直是植物表型检测的重要研究内容，而三维原位构型的自动检测和分析方法是该领域目前亟待解决的一个技术瓶颈。因此，对植物根系三维构型参数进行原位、无损、准确、自动和快速地测量具有重要的意义。挖掘法、钉板法、土钻法以及土柱法等传统根系研究方法均无法满足要求，而借助X射线计算机断层摄影术（X-ray computed tomography，简称XCT）技术则有望实现上述目标[1]。XCT技术被应用于植物根系研究已有30多年历史[2-3]，但由于成像分辨率低等因素的限制，过去大多数研究都偏向图像分割与可视化方面[4-6]，虽然曾有研究人员基于XCT技术测量植物根系构型参数，但效果并不理想[7]。随着XCT技术的发展，其成像分辨率已大大提高，用其对植物根系三维构型参数进行准确的测量已经成为可能，因此近年来基于XCT技术测量植物根系构型参数的研究逐渐增多[8-9]。然而，相关软件的发展并没有跟上XCT硬件的发展，目前依然缺乏稳定可靠的基于XCT技术的植物根系三维构型参数自动测量系统[10]。
针对上述背景，并借鉴田绪红等提出的基于横截面算法的三维植物根系图像骨架生成方法[11]和李骈臻等提出的基于骨架模型的植物根系三维构型可视化方法[12]，本研究提出基于XCT序列图像的植物根系三维矢量模型构建方法，并基于所构建的矢量模型实现根数、根长度、根体积、根表面积以及根夹角的自动化测量，该方法的实现将为植物根系三维构型参数的原位、无损、准确、自动和快速测量提供一个全新的工具，对植物根系的定量研究具有重要意义。
1植物根系三维矢量模型构建的基本原理
首先，用160 kV XCT系统采集根系样本断层序列图像，然后对其进行滤波和分割；其次，将已分割的序列图像读进计算机内存，组成三维体数据，对其重切以提取重切片[13]，根据重切片中各根截面的质心和面积构建初始根结点，并根据根系的形态特征设计算法对上述根结点分组以构建根分枝；再次，对各根分枝进行曲线拟合，以固定间隔采样拟合曲线上的点作为新根心点，并在各新根心点处再次重切（切面与拟合曲线垂直且经过新根心点）三维体数据以获取新根截面，根据新根截面的面积、等效圆半径以及等效圆周长等信息构建新根结点以取代各根分枝的初始根结点；另外，通过判定各根分枝的连通性建立它们的拓扑关系以构建完整的根系统（植物根系三维矢量模型）；最后，在上述矢量模型的基础上实现三维构型参数的计算。
2植物根系三维矢量模型的数据结构设计
为实现上述算法步骤，首先须要设计相应的数据结构。根据植物根系的形态特征，本研究设计了RootNode、RootBranch、RootSystem等3个基本结构，分别代表根结点、根分枝、根系统。根系统包含根分枝，根分枝包含属于它的根结点列表。
2.1RootNode
RootNode包含根结点号（NodeId）、所属根分枝号（BranchId）、根心坐标（Centroid）、法向量（Normal）、面积（Area）、半径（Radius）、周长（Perimeter）、根截面圆率（Roundness）以及根心线在该根结点处的曲率（Curvature）等信息。
在三维空间中用垂直于根分枝生长方向的平面（切面）在特定位置截断该根分枝，将得到1个根截面（图1），然后可根据此根截面的质心、面积以及法向量等信息构建根结点，具体方法为将根截面质心进行坐标变换[14]得到的三维坐标作为根结点中心（根心点）；根结点的法向量取值为该根截面的法向量；根结点面积取值为该根截面的面积。
2.2RootBranch
RootBranch包含根分枝号（BranchId）、父根分枝号（ParentBranchId）、根结点号列表（NodeIds）、子根分枝号列表（ChildBranchIds）等信息。根分枝的构建，主要通过对上述根结点进行分组实现。
2.3RootSystem
RootSystem包含所有根结点、根分枝、计算特定根分枝参数（长度、体积以及表面积等）的方法以及计算根系统整体参数（根总长度、根总表面积以及根总体积等）的方法等。在根分枝的基础上，进一步建立它们的拓扑关系，即可构建完整的根系统。
3植物根系三维矢量模型构建方法
主要算法流程如图2所示。
3.1构建初始根结点
原XCT序列图像中各根截面是植物根系的横断面[16]，其质心接近于纵向生长的根分枝中心，但与横向生长的根分枝中心之间存在较大误差。因此，本研究采用基于平行轮廓线的重建方法[17]对植物根系XCT序列图像进行三维重建，然后分别沿Z轴、X轴和Y轴3个方向对其重切以获取重切片，再根据重切片中各根截面的质心和面积构建初始根结点，重切过程中若刚好切中侧根，也会导致根截面质心和实际根心[CM（25]之间存在较大误差，从而产生根心偏移现象（图3）。经过分析发现，根截面圆率越接近1.0，其质心与实际根心之间的误差就越小。因此可通过舍去圆率小于特定阈值的根截面以抑制根心偏移现象。根结点的根心效果如图4所示。
通过上述步骤构建了相互独立的根结点，这里进一步根据最短路径算法的思想[18]，加上连通性、面积变化率以及偏转角等约束条件，对它们进行分组以构建初始根分枝。
连通性：在2根结点中心之间作1直线段Lstraight，若该线段上的任意坐标点都落在根系体内，则认为这2个根结点是直连通的；在2根结点中心点之间作满足特定条件的任意抛物线Lstraight（抛物线经过这2根结点，它的顶点在垂直于Lstraight的平面上且与Lstraight之间的距离不能大于Lstraight长度的50%），只要其中1条抛物线上的任意坐标点都落在根系体内，则认为这2个根结点是曲连通。
面积变化率：如式（4）、式（5）、式（6）所示，其中areai和areai+1分别是根结点i和根结点（i+1）的面积，ratearea为面积变化率。
偏转角：如图5所示，nodei-1、nodei、nodei+1为相邻的3个根结点。依次在相邻2个根结点之间连接直线所形成的夹角θ即为偏转角。
有了上述定义，可将根结点分组的主要步骤描述为（1）若待分组根结点数小于1个，则停止分组过程；否则创建1个新根分枝，并在待分组根结点中随机选择1个作为新根分枝起点，然后将待分组根结点数减1，记录当前搜索方向为向后的；（2）若当前的搜索方向为向后的，将新根分枝最后1个根结点作为当前根结点，否则将新根分枝第1个根结点作为当前[CM（25]根结点；（3）在待分组根结点中搜索满足以下条件的根结
[FK（W9][TPLYH555.tif]
点作为新根分枝下一个候选根结点：（a）候选根结点与当前根结点之间的距离小于指定阈值；（b）候选根结点与当前根结点是连通的（直连通或者曲连通）；（c）候选根结点与当前根结点的面积变化率小于指定阈值；（d）当新根分枝的根结点数≥2个时，候选根结点与新根分枝前2个根结点所形成的偏转角小于指定阈值；（4）若候选根结点数≥1个，则将与当前根结点距离最小的候选根结點作为新根分枝的下一根结点（若搜索方向为向后的，则将新根结点添加到新根分枝最后的位置，否则将其添加到新根分枝最前的位置），然后跳回步骤（2）；（5）若候选根结点数等于0，则分2种情况处理，如果当前搜索方向为向后的，则改变为向前的，然后跳回步骤（2）；若当前搜索方向为向前的，则结束本根分枝的搜索，并跳回步骤（1）。
分组效果如图6所示，其中棕色皮肤为三维重建效果，蓝色点、红色点都是根结点中心（根心点），蓝色根心点被上述算法划分为同一组。由图6可见，分组结果与实际根分枝情况一致。
3.3重构根结点
虽然在构建初始根结点的步骤中剔除了误差较大的根心点，但是为了保证属于同一根分枝的相邻根结点间的距离不至于太大，以便于后续对它们进行分组，保留了误差相对较小的根心点，因此所得的根心点还是会存在一定程度的偏移（图7-a）。 上述根心偏移现象可以通过3次B样条拟合[14]来校正。具体方法是将各根分枝中除了首末根结点外圆率小于特定阈值的根结点都舍弃，然后用剩余的根结点作为控制点进行3次B样条拟合，并以特定间隔采样拟合曲线上的点作为新根心点以取代初始根心点（图7-b），从而构建没有面积的新根结点。
由于3次B样条曲线具有二阶连续性，可以计算曲线在各根结点处的一阶导数和二阶导数[19]，继而计算其曲率[CM（25][20]，即可分别得到新根结点的法向量（Normal）、曲率值。
3.4重获根结点面积
上述新根结点的面积须要通过对三维体数据（XCT序列图像）再次重切来获得。重切时将切面中心设为根结点的中心，并将切面的法向量设为根结点的法向量，从而使得切面垂直于根分枝生长方向。为了避免侧根对根截面积的影响，重切过程中须要舍弃圆率小于特定阈值的根截面。
由于舍去了部分根截面，导致某些根结点无法通过重切获得面积，对于这些根结点，可根据相邻根结点的半径进行线性插值来获得相应的半径，然后用该半径计算其面积和周长。
3.5判定根分枝关系
通过上述步骤构建了相互独立的根分枝，这里根据植物根系的形态特征设计算法进一步建立根分枝之间的关系，主要有确定分叉根结点和确定根分枝关系2个任务。
3.5.1确定分叉根结点的算法步骤
在父根分枝中找出所有与子根分枝首根结点具有连通性（直连通或曲连通）的根结点作为候选分叉根结点；设置由分叉根结点指向子根分枝首根结点的方向向量为Vd，子根分枝首根结点法向量为Vf（图8），它们之间的夹角为φ（图8）。依次计算上述各候选分叉根结点的φ值，然后选择φ值最小的候选分叉根结点作为最终分叉根结点。
3.5.2判定根分枝关系的算法步骤
（1）将第i（初始值为0）条根分枝作为当前根分枝；（2）调整当前根分枝中根结点的存储顺序，使得半径较大的根结点存储在较前的位置；（3）遍历其余根分枝，找出与当前根分枝首根结点具有指定连通性（直连通或曲连通）的根分枝添加到候选父根分枝列表中；（4）在当前根分枝首根结点与候选父根分枝的分叉根结点间作1直线段，若该直线段所经过的根分枝数大于2个，则将此根分枝从候选父根分枝列表中剔除；（5）比较当前根分枝首根结点的面积与候选根分枝分叉根结点面积，若分叉根结点面积小于当前根分枝首根结点面积，则将该根分枝从候选父根分枝列表中剔除；（6）在候选父根分枝列表中，选择其分叉根结点面积最大的根分枝作为当前根分枝的最终父根分枝；（7）递增i，若i≥根分枝数则结束搜索，否则返回步骤（1）。
3.5.3判定效果
如图9所示，蓝色根分枝为当前选中的根分枝，红色圈内蓝色线与红色线相交位置即为其父根分枝分叉根结点位置。由图9的判定效果可见，由上述算法建立的根分枝关系及所确定的分叉根结点位置与实际根系形态一致，取得较好效果。
4植物根系三维矢量模型基本构型参数计算方法
基于上述步骤构建的植物根系三维矢量模型，可计算根分枝数、根长度、根体积、根表面积以及根夹角。
4.1根分枝数
即植物根系三维矢量模型的根分枝数。
4.2根长度
同一根分枝相邻根结点的距离如式（7）所示。式（7）中（xi，yi，zi）、（xi+1，yi+1，zi+1）分别为第i根结点、第（i+1）根结点坐标，i为根结点在根分枝中的位置序号，取值为0，1，2，…，（n-1），其中n为根分枝的根结点数；D（i，i+1）为第i根结点与第（i+1）根结点之间的距离。
6结论
本研究提出并实现了基于XCT序列图像的植物根系三维矢量模型构建方法，在所构建矢量模型的基础上实现根数、根长度、根体积、根表面积以及根夹角的自动化测量，将上述自[CM（25]动测量值与手工测量值进行对比得知误差绝对值小于7%，从而验证了本研究所提出的方法的可行性与准确性。
参考文献：
[1]孙曰波，赵从凯. 根系研究方法进展[J]. 潍坊高等职业教育，2009，5（1）：52-55
[2]罗锡文，周学成，严小龙，等. 基于XCT技术的植物根系原位形态可视化研究[J]. 农业机械学报，2004，35（2）：104-106.
[3]Mooney S J，Pridmore T P，Helliwell J，et al. Developing X-ray computed tomography to non-invasively image 3-D root systems architecture in soil[J]. Plant Soil，2012，352（1/2）：1-22.
[4]周学成，罗锡文. 基于遗传算法的原位根系CT图像的模糊阈值分割[C]//中国农业工程学会学术年会论文摘要集. 大庆，2007：110.
[5]陈郁淦，周学成，罗锡文，等. 基于ITK和VTK的原位根系三维可视化研究[C]//中国农业工程学会学术年会论文集. 重庆，2011：1164-1168.
[6]李克新，李沐阳，薛瑞，等. 林木幼苗根系CT序列图像分割[J]. 森林工程，2014，30（1）：25-29.
[7]Heeraman D A，Hopmans J W，Clausnitzer V. Three dimensional imaging of plant roots in situ with X-Ray computed tomography[J]. Plant and Soil，1997，189（2）：167-179.
[8]Koebernick N，Weller U，Huber K，et al. In situ visualization and quantification of three dimensional root system architecture and growth using X-ray computed tomography[J]. Vadose Zone Journal，2014，13（8）：1-10.
[9]Metzner R，Eggert A，Dusschoten D V，et al. Direct comparison of MRI and X-ray CT technologies for 3D imaging of root systems in soil：potential and challenges for root trait quantification[J]. Plant Methods，2015，11：1-11.
[10]Kuijken R C P，van Eeuwijk F A，Marcelis L F M，et al. Root phenotyping：from component trait in the lab to breeding[J]. Journal of Experimental Botany，2015，66（18）：5389-5401.
[11]田绪红，李志垣，韩国强，等. 基于横截面算法的三维植物根系图象骨架生成方法[C]//第十二届全国图象图形学学术会议论文集. 北京，2005：655-658.
[12]李骈臻，周学成，张常玲，等. 基于骨架模型的植物根系三维构型可视化方法[J]. 计算机工程与设计，2014，35（11）：3913-3917.
[13]周振环，伍云智，赵明. 医学图像编程技术[M]. 北京：电子工业出版，2010.
[14]Hearn D，Baker M P. 计算机图形学[M].蔡士杰，译.3版. 北京：电子工业出版社，2013.
[15]Lehmann G. Label object representation and manipulation with ITK[J/OL]. Insight Journal，2008，1-34. http：//www.insight-journal.org/browse/publication/176
[16]周振环，郑小中，赵明. 三维图像编程实验[M]. 北京：电子工业出版社，2011.
[17]张俊华. 医学图像三维重建和可视化——VC++实现实例[M]. 北京：科学出版社，2014.
[18]程杰. 大话数据结构[M]. 北京：清华大学出版社，2011.
[19]江本赤. B样条曲线曲率简易求解算法[J]. 制造技术与机床，2014（10）：78-79.
[20]张学东. 空间曲线的曲率计算方法[J]. 塔里木农垦大学学报，2002，14（2）：37.
[21]方明亮，郭正光. 高等數学[M]. 广州：广东科技出版社，2008.