基于XCT序列圖像的植物根系三維矢量模型構建方法

李一海　周學成　陳富強

摘要：為實現植物根系三維構型參數的原位、無損、準確、快速和自動化測量，提出一種基于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.