基于MicroCT的骨小梁參數測量系統的應用效果分析

吳沛澤，羅守華，陳功，柳慧芬，陳斌

1.東南大學 生物科學與醫學工程學院，江蘇 南京 210018；2.南京中醫藥大學附屬醫院，江蘇 南京 210029；3.南京市口腔醫院 南京大學醫學院附屬口腔醫院，江蘇 南京 210001

基于MicroCT的骨小梁參數測量系統的應用效果分析

吳沛澤1，羅守華1，陳功2，柳慧芬3，陳斌3

1.東南大學 生物科學與醫學工程學院，江蘇 南京 210018；2.南京中醫藥大學附屬醫院，江蘇 南京 210029；3.南京市口腔醫院 南京大學醫學院附屬口腔醫院，江蘇 南京 210001

為了更全面便捷地衡量動物臨床試驗骨質疏松情況，開發了一套基于MircoCT圖像的用于測量骨小梁參數的軟件系統。在保留骨結構完好的前提下，使用12組成年豬額骨及頜骨樣本進行MicroCT掃描成像，并根據二維斷層圖像，經過導入數據、三維可視化、圈定感興趣區、二值化處理等步驟后進行骨密度、骨小梁厚度等參數測量。測量結果與BoneJ進行對比，結果驗證了本系統具有較高的測量準確性（各項參數相對誤差均＜5%），可視化效果直觀，在動物骨小梁參數測量方面具有較高的應用價值。

骨小梁；距離變換；三維可視化；骨密度；骨小梁厚度

骨小梁是松質骨的主要結構，其形態結構參數是研究松質骨病變，判斷骨質疏松情況的重要參考依據。研究動物骨小梁在試驗前后的變化情況對于臨床醫學具有重要意義。MircoCT（又稱微型CT或顯微CT），是指可測定體素空間分辨率小于100的X線三維成像系統，具有空間分辨率高、焦點小、輻射劑量少等特點，是研究動物骨小梁結構的有力輔助工具[1]。而對于動物骨小梁的結構參數測量，除了比利時SkyScan有配合測量其MircoCT圖像的對應骨分析軟件以外，目前國內對動物骨小梁結構多數還在使用顯微切片或掃描電鏡觀察等二維圖像分析方法，缺少對動物骨小梁完整三維結構測量和觀察的輔助工具[2]。因此，針對MircoCT掃描骨小梁圖像，開發了一套基于C++的骨小梁參數測量系統，能夠測得骨密度、骨小梁厚度、骨小梁間隙、骨小梁數量等多個參數，可以對動物骨病變情況有更全面的衡量。

1 測量原理

1.1 測量骨密度原理

骨密度，全稱骨骼礦物質密度，是反映骨質疏松程度，衡量骨骼健康情況最為重要的一個參數。使用定量CT（QCT）所測密度為該區域的體密度，單位為g/cm3。

CT成像中當X線束進入物體時，將受到物體對X線的吸收和散射，并且以吸收為主。X線束透過物體后剩余能量將被檢測器所接受作為投影數據，投影數據經過CT重建后獲得最終CT掃描圖像，因此掃描圖像的CT值可以直接反映物體對X線束的吸收情況。而物體對X線的吸收程度與物體密度、原子系數以及X線的吸收能量等參數密切相關[3]。在同一CT相同電壓電流的拍攝條件下，可以保證X線的吸收能量相同。使用專業制造與骨成分近似的模體對比拍攝，可以保證原子系數相同。在控制其他變量的情況下可視為物體對X線的吸收程度與物體的密度成正比。即得到以下公式：

其中，CTLevel表示當前測量區域平均CT值，BMDh表示模體高密度部分密度，BMD1表示模體低密度部分密度，CTh表示模體高密度部分測量CT值，CT1表示模體低密度部分測量CT值。計算結果的單位通常用mg/cm3表示。其中，對于比較好的骨密度模體，有不止兩個高低密度值的情況，也可以采用最小二乘法來求得目標區域骨密度。

1.2 測量骨小梁厚度原理

骨小梁厚度指的是骨小梁的平均厚度，是衡量骨小梁微觀形態最主要的參數之一，它與骨小梁間隙和骨小梁數量這兩個參數在計算上是緊密聯系的。結合骨小梁間隙（Th.Sp）與骨小梁數量（Th.b）這兩個參數可以判斷骨質疏松程度。當發生骨質疏松時，骨小梁厚度值會降低，而骨小梁間隙增大，骨小梁數量減小。

計算骨小梁厚度可采用Hildebrand和Ruegsegger[4]（1997）所提出的算法。這種算法基于MicroCT拍攝的3D圖像分析，不同于過去借助形狀建模的算法，這種獨立于模型的算法更加精確。根據該算法，對物體上一個點的局部厚度可以定義為滿足以下兩個條件的最大球體的直徑：① 球體包含該點（該點不一定是球心）；② 球體邊界在物體內。

基于這個定義，一段骨小梁結構的厚度即骨小梁內骨架上所有點局部厚度的平均值。而圖像骨架可以使用距離變換法（Distance Transform）求出[5]。距離變換法指的是對二值圖像，計算每個前景像素最近的背景像素距離，將其變換為灰度圖像，即距離圖像的一種算法。

通過距離變換求得灰度圖像的峰值即原二值圖像的骨架，對骨架上的點遍歷計算局部最大球體直徑，求和平均即得到骨小梁厚度。計算骨小梁間隙的過程與骨小梁厚度相同，只是需要先將圖像取反再進行運算。而骨小梁平均數量則是由以下公式得出[6]：

其中，Tb.Th表示骨小梁厚度，Tb.Sp表示骨小梁間隙，因此骨小梁數量常用單位為 1/mm。

1.3 測量骨小梁結構指數原理

結構模型指數（Structure Model Index，SMI）表達的是骨小梁的三維形態相對而言接近圓盤狀或圓柱狀的程度。當出現骨質疏松時，骨小梁的形態結構將會從盤狀向桿狀轉變，而這個參數就能用于衡量骨小梁的形變程度。一個理想的圓盤其SMI值為0，一個理想的圓柱其SMI值為3，一個理想的球體其SMI值為4。相應地，一個理想的圓柱狀孔洞SMI值為-3，一個理想的圓柱狀孔洞SMI值為-4。因為SMI涉及到表面凸面曲率的測量，如果目標是凹面，則其測量值為負數。

計算SMI是基于三維體素模型的膨脹求得的。首先將目標松質骨進行二值化，計算其體積與表面積，然后對二值化后物體進行三維膨脹，再次計算體積與表面積（這個步驟與計算形狀指數Tb.Pf的過程相同），然后SMI可用如下公式求出[7]：

需要說明的是對于封閉孔洞而言，其表面是凹陷的，因此求得的SMI也是負數，由于對這樣的封閉空間膨脹后物體的表面積變小了，所求得的表面積差S'也就是負數。因此對于一個包含封閉孔洞超過50%以上空間的物體（如骨小梁），最終SMI參數將是一個負數。也可以說，參數SMI與體積百分比關系緊密。由于人為圈定目標區域產生的邊和角改變了測量物體的體積，也會增大測量的SMI值。

2 測量方法

2.1 材料準備及圖像獲取

實驗儀器：本次實驗使用的MircoCT為東南大學生物科學與醫學工程學院利用國家自然科學基金科學儀器研究專項資金研發的高分辨顯微CT（Hiscan M1000），最大分辨率為20 μm，最大掃描范圍可達到60 mm×200 mm（多次掃描）。圖像重建及后處理軟件為蘇州海斯菲德信息技術有限公司提供的MCT_Recon_Sever（版本號1.0.10）。

材料制備：從正常豬體內使用鋸子分離獲得頜骨12組，然后采用濃度4%的福爾馬林浸泡，保持其骨結構不變，骨礦含量不流失（由于骨組織與非骨組織在圖像灰度值上有很大差異，因此在掃描前不需要特意剔除骨表面的其他組織）。用鋸子將頜骨樣本切割為長度不超過15 cm的骨塊，使其能夠完全放入CT的聚丙烯樣本槽。

掃描過程：將12組豬頜骨樣本依次置于聚丙烯樣本槽中，將掃描X 射線的能量設置為400 kV，電流設置為200 μA，掃描旋轉360°，采集720張投影圖，每張投影圖的曝光時間為50 ms，每個樣本平均掃描時間3min。使用MCT_Recon_ Sever對每個樣本重建，平均獲得800張大小為1874×1874像素的斷層重建數據。重建圖像分辨率為25 μm，Z方向間隔25 μm。

2.2 測量流程

測量流程，見圖1。

圖1 測量流程示意圖

（1）首先導入CT斷層掃描圖像數據，導入數據將以三視圖方式顯示。在三視圖中選擇目標區域，可以對該區域數據進行三維體繪制顯示。如果數據是已經分割完成的感興趣區（Region of Interest，ROI）數據，則可以直接進入二值化處理環節。在三維可視化環節，通過旋轉位移調整物體，可以直觀地觀察骨小梁形態，選擇最合適的角度獲得ROI的二維斷層圖像。然后在二維斷層圖像上手動圈定ROI區域。圈定時程序會根據手動圈定結果對中間所有層進行插值運算獲得所有層的圈定結果。圈定ROI完成后，將灰度圖像進行二值化并保存數據（灰度圖像數據保留在內存中不會刪除），然后按照上文所述算法進行不同的參數測量。

導入數據：由于圖像分辨率較高，單張DCM圖像大小可達到5～6 MB，全部數據導入將占用1 GB以上內存導致程序卡頓，因此在導入數據過程中并非將圖像數據全部讀入內存，而是使用了“預加載”技術，只讀入DCM圖像中的tag信息，在實際需要該層圖像數據時再讀取真正的圖像數據，極大提升了加載速度。同時考慮到內存空間限制，對于過大的數據會進行自動降采樣減小內存占用。

（2）三維可視化的操作界面，見圖2。為了直觀地進行ROI圈定和觀察目標松質骨的骨質疏松情況，需要先對圖像進行三維可視化。使用OpenGL庫進行體繪制來實現三維可視化，其具體方法是導入數據經采樣和插值后，將每個體素構造為獨立的理想化物理模型，考慮其介質屬性，按照phone光照模型調節傳遞函數對每個體素分配光強和不透明度，然后沿視線觀察方向積分，最后形成半透明的三維投影圖像。

圖2 三維可視化示意圖

（3）圈定ROI的操作界面，見圖3。計算骨密度、骨小梁厚度等參數，需要保證計算的圖像范圍只包含目標松質骨而不含有其他結構或空腔，否則會導致計算結果偏小。因此需要對ROI進行準確地圈定。圈定ROI時，首先將三維顯示的松質骨結構按照垂直屏幕方向獲得二維斷層圖像，瀏覽斷層圖像選擇ROI起始層和終止層，在起始層和終止層上再次手動圈定區域，將會在中間每一層按照漸變插值計算出該層的圈定區域，若中間層的圈定效果不滿意，可以在該層手動圈定，然后按照新的圈定結果重新計算每一層的圈定區域。圈定完成后，程序將ROI區域以外的數據清零，保證進行計算時不受到其他結構影響，并自動將當前圖像保存為新的DCM序列便于再次測量和查看結果。

圖3 圈定ROI示意圖

（4）二值化處理的操作界面，見圖4。在計算骨結構參數前，需要將骨組織和非骨組織區分，故先將圖像進行二值化處理。具體過程為通過手動對斷層圖像調窗，將非骨部分數據置為0，將骨部分數據置為1，同時將調窗結果同步至三維可視化窗口，直觀顯示調窗效果，便于判斷是否將骨組織和非骨組織區分開。圖像二值化數據可以保存至硬盤便于再次查看和測量。

圖4 二值化處理示意圖

（5）參數測量的操作界面，見圖5。參數測量分為“骨密度測量”和“骨小梁參數測量”兩個部分，其中“骨密度測量”使用未二值化的原始圖像數據，通過計算區域灰度平均值，再導入模體圖像數據，計算模體圖像灰度值并輸入模體密度，按照第一節所述原理進行密度計算。“骨小梁參數測量”使用原始圖像二值化后數據，經過去噪（消除圖像中由于CT掃描造成的非骨部分的噪點），按照第一節所述原理對所有骨小梁參數進行測量。最后輸出測量結果并生成測量報告。

圖5 測量結果示意圖

3 測量結果

3.1 可視化效果

豬頜骨三維效果示意圖，見圖6。在右側三視圖中框選出松質骨的大致區域后，可以看到完整松質骨的三維可視化效果，左下角標記了當前物體的旋轉情況，確定好旋轉角度后將按照垂直屏幕方向生成新的斷層圖像。相比在CT斷層圖像上選擇ROI的方法，使用三維可視化進行ROI選擇效果更加清楚直觀，操作更加方便準確。

3.2 軟件測量結果與BoneJ誤差對比

ImageJ是一個基于Java的公共的圖像處理軟件[8]，它是由美國國立衛生研究院（National Institutes of Health，NIH）的Wayne Rasband等人開發，可運行于Microsoft Windows、Mac OS、Mac OS X、Linux和Sharp Zaurus等多種平臺，在世界范圍的科學研究領域中廣泛應用。其中BoneJ屬于ImageJ中專門測量骨結構參數的插件[9]，由 Doube M、K?osowski MM、Arganda-Carreras I等人研發。其測量結果具有一定的權威性和參考價值，本文采用BoneJ測量來驗證本系統的測量準確度。

測量結果見表1，與BoneJ測量結果對比，12組標本骨密度的平均相對誤差為0.0%，骨小梁厚度的平均相對誤差為4.275%，骨小梁間隙的平均相對誤差為1.127%，骨小梁數量的平均相對誤差為1.952%，骨小梁結構指數的平均相對誤差為3.954%。

分析誤差原因，主要有如下幾個方面：① 計算骨表面積方法不同，BoneJ中使用三角面片法計算骨表面積，而本系統使用體素統計計算表面積。表面積計算的誤差直接影響到參數計算結果；② 圈定ROI時有部分骨組織被截斷，形成邊際效應影響測量結果。

圖6 豬頜骨三維效果示意圖

4 結論

本測量系統能夠對MicroCT斷層掃描獲得的動物骨小梁圖像進行骨密度、骨小梁厚度、骨結構指數等參數進行全面準確的測量，并提供骨小梁結構的三維可視化效果，為臨床研究動物骨質疏松程度、骨健康情況提供了一個更全面的參考工具。

表1 豬頜骨各項參數測量結果

[1]高鵬,戎軍艷,廖琪梅,等.MicroCT系統軟件平臺的設計與實現[J].醫療衛生裝備,2014,35(12):22-24.

[2]付鑫,馬劍雄,董寶康.骨微結構檢測方法研究進展[J].中國骨與關節外科,2009,2(6):509-515.

[3]Ma XQ,Overton TR.Automated Image Analysis for Bone Density Measurement Using Computed Tomography[J].IEEE Trans Med Imaging,1991,10(4):611-616.

[4]Hildebrand T,Rüegsegger P.A new method for the modelindependent assessment of thickness in three-dimensional images[J].J Microsc,1997,185(1):67-75.

[5]Borgefors G.Distance transformation in arbitrary dimensions[J].Comput Vis Graph Image Process,1984,27(3):321-345.

[6]Hahn M,Vogel M,Pompesius-Kempa M,et al.Trabecular bone pattern factor-a new parameter for simple quantification of bone microarchitecture[J].Bone,1992,13(4):327-330.

[7]Hildebrand T,Rüegsegger P.Quantifcation of Bone Microarchirecture with the Structure Model Index[J].Comput Methods Biomech Biomed Engin,1997,1(1):15-23.

[8]Schneider CA,Rasband WS,Eliceiri KW.NIH Image to ImageJ:25 years of image analysis[J].Nat Methods,2012,9(7):671-675.

[9]Doube M,K?osowski MM,Arganda-Carreras I,et al.BoneJ:Free and extensible bone image analysis in ImageJ[J].Bone,2010,47 (6):1076-1079.

Analysis of the Application Effects of Trabecular Bone Measurement System Development Based on MircoCT

In order to measure osteoporosis cases in animal clinical trials,the research developed a trabecular bone parameters measurement software system based on MircoCT images.On the premise of keeping skeletal structure intact,the paper carried on MicroCT scan on 12 samples of adult pig frontal bones and maxillary bones.Based on the two-dimensional sectional images,measurement of parameters such as bone density and trabecular thickness was conducted through the process of data import,threedimensional visualization,and delineation of region of interests,and binary image processing.The results of the measurement was compared with BoneJ,which confrmed the accuracy of measurement system (relative error ＜5% for each parameter).The system also allowed intuitive visualization and had a high application value in terms of animal trabecular parameter measurement.

trabecular;distance transform;three-dimensional visualization;bone density;trabecular thickness

WU Pei-ze1,LUO Shou-hua1,CHEN Gong2,LIU Hui-fen3,CHEN Bin3
1.Department of Biological Science and Medical Engineering,Southeast University,Nanjing Jiangsu 210018,China;2.The Affiliated Hospital of Nanjing University of Traditional Chinese Medicine,Nanjing Jiangsu 210029,China;3.Nanjing Stomatological Hospital,Medical School of Nanjing University,Nanjing Jiangsu 210001,China

R681；R814.42

A

10.3969/j.issn.1674-1633.2016.04.009

1674-1633(2016)04-0045-04

2016-01-12

2016-02-18

國家自然科學基金（61127002，61179035）；中央高校基本業務費重點項目培育計劃（021414310017）。

作者郵箱：owenwu_mx@163.com