基于MicroCT的骨小梁參數(shù)測(cè)量系統(tǒng)的應(yīng)用效果分析

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

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

基于MicroCT的骨小梁參數(shù)測(cè)量系統(tǒng)的應(yīng)用效果分析

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

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

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

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

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

1 測(cè)量原理

1.1 測(cè)量骨密度原理

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

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

其中，CTLevel表示當(dāng)前測(cè)量區(qū)域平均CT值，BMDh表示模體高密度部分密度，BMD1表示模體低密度部分密度，CTh表示模體高密度部分測(cè)量CT值，CT1表示模體低密度部分測(cè)量CT值。計(jì)算結(jié)果的單位通常用mg/cm3表示。其中，對(duì)于比較好的骨密度模體，有不止兩個(gè)高低密度值的情況，也可以采用最小二乘法來(lái)求得目標(biāo)區(qū)域骨密度。

1.2 測(cè)量骨小梁厚度原理

骨小梁厚度指的是骨小梁的平均厚度，是衡量骨小梁微觀形態(tài)最主要的參數(shù)之一，它與骨小梁間隙和骨小梁數(shù)量這兩個(gè)參數(shù)在計(jì)算上是緊密聯(lián)系的。結(jié)合骨小梁間隙（Th.Sp）與骨小梁數(shù)量（Th.b）這兩個(gè)參數(shù)可以判斷骨質(zhì)疏松程度。當(dāng)發(fā)生骨質(zhì)疏松時(shí)，骨小梁厚度值會(huì)降低，而骨小梁間隙增大，骨小梁數(shù)量減小。

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

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

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

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

1.3 測(cè)量骨小梁結(jié)構(gòu)指數(shù)原理

結(jié)構(gòu)模型指數(shù)（Structure Model Index，SMI）表達(dá)的是骨小梁的三維形態(tài)相對(duì)而言接近圓盤狀或圓柱狀的程度。當(dāng)出現(xiàn)骨質(zhì)疏松時(shí)，骨小梁的形態(tài)結(jié)構(gòu)將會(huì)從盤狀向桿狀轉(zhuǎn)變，而這個(gè)參數(shù)就能用于衡量骨小梁的形變程度。一個(gè)理想的圓盤其SMI值為0，一個(gè)理想的圓柱其SMI值為3，一個(gè)理想的球體其SMI值為4。相應(yīng)地，一個(gè)理想的圓柱狀孔洞SMI值為-3，一個(gè)理想的圓柱狀孔洞SMI值為-4。因?yàn)镾MI涉及到表面凸面曲率的測(cè)量，如果目標(biāo)是凹面，則其測(cè)量值為負(fù)數(shù)。

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

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

2 測(cè)量方法

2.1 材料準(zhǔn)備及圖像獲取

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

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

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

2.2 測(cè)量流程

測(cè)量流程，見圖1。

圖1 測(cè)量流程示意圖

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

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

（2）三維可視化的操作界面，見圖2。為了直觀地進(jìn)行ROI圈定和觀察目標(biāo)松質(zhì)骨的骨質(zhì)疏松情況，需要先對(duì)圖像進(jìn)行三維可視化。使用OpenGL庫(kù)進(jìn)行體繪制來(lái)實(shí)現(xiàn)三維可視化，其具體方法是導(dǎo)入數(shù)據(jù)經(jīng)采樣和插值后，將每個(gè)體素構(gòu)造為獨(dú)立的理想化物理模型，考慮其介質(zhì)屬性，按照phone光照模型調(diào)節(jié)傳遞函數(shù)對(duì)每個(gè)體素分配光強(qiáng)和不透明度，然后沿視線觀察方向積分，最后形成半透明的三維投影圖像。

圖2 三維可視化示意圖

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

圖3 圈定ROI示意圖

（4）二值化處理的操作界面，見圖4。在計(jì)算骨結(jié)構(gòu)參數(shù)前，需要將骨組織和非骨組織區(qū)分，故先將圖像進(jìn)行二值化處理。具體過(guò)程為通過(guò)手動(dòng)對(duì)斷層圖像調(diào)窗，將非骨部分?jǐn)?shù)據(jù)置為0，將骨部分?jǐn)?shù)據(jù)置為1，同時(shí)將調(diào)窗結(jié)果同步至三維可視化窗口，直觀顯示調(diào)窗效果，便于判斷是否將骨組織和非骨組織區(qū)分開。圖像二值化數(shù)據(jù)可以保存至硬盤便于再次查看和測(cè)量。

圖4 二值化處理示意圖

（5）參數(shù)測(cè)量的操作界面，見圖5。參數(shù)測(cè)量分為“骨密度測(cè)量”和“骨小梁參數(shù)測(cè)量”兩個(gè)部分，其中“骨密度測(cè)量”使用未二值化的原始圖像數(shù)據(jù)，通過(guò)計(jì)算區(qū)域灰度平均值，再導(dǎo)入模體圖像數(shù)據(jù)，計(jì)算模體圖像灰度值并輸入模體密度，按照第一節(jié)所述原理進(jìn)行密度計(jì)算?！肮切×簠?shù)測(cè)量”使用原始圖像二值化后數(shù)據(jù)，經(jīng)過(guò)去噪（消除圖像中由于CT掃描造成的非骨部分的噪點(diǎn)），按照第一節(jié)所述原理對(duì)所有骨小梁參數(shù)進(jìn)行測(cè)量。最后輸出測(cè)量結(jié)果并生成測(cè)量報(bào)告。

圖5 測(cè)量結(jié)果示意圖

3 測(cè)量結(jié)果

3.1 可視化效果

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

3.2 軟件測(cè)量結(jié)果與BoneJ誤差對(duì)比

ImageJ是一個(gè)基于Java的公共的圖像處理軟件[8]，它是由美國(guó)國(guó)立衛(wèi)生研究院（National Institutes of Health，NIH）的Wayne Rasband等人開發(fā)，可運(yùn)行于Microsoft Windows、Mac OS、Mac OS X、Linux和Sharp Zaurus等多種平臺(tái)，在世界范圍的科學(xué)研究領(lǐng)域中廣泛應(yīng)用。其中BoneJ屬于ImageJ中專門測(cè)量骨結(jié)構(gòu)參數(shù)的插件[9]，由 Doube M、K?osowski MM、Arganda-Carreras I等人研發(fā)。其測(cè)量結(jié)果具有一定的權(quán)威性和參考價(jià)值，本文采用BoneJ測(cè)量來(lái)驗(yàn)證本系統(tǒng)的測(cè)量準(zhǔn)確度。

測(cè)量結(jié)果見表1，與BoneJ測(cè)量結(jié)果對(duì)比，12組標(biāo)本骨密度的平均相對(duì)誤差為0.0%，骨小梁厚度的平均相對(duì)誤差為4.275%，骨小梁間隙的平均相對(duì)誤差為1.127%，骨小梁數(shù)量的平均相對(duì)誤差為1.952%，骨小梁結(jié)構(gòu)指數(shù)的平均相對(duì)誤差為3.954%。

分析誤差原因，主要有如下幾個(gè)方面：① 計(jì)算骨表面積方法不同，BoneJ中使用三角面片法計(jì)算骨表面積，而本系統(tǒng)使用體素統(tǒng)計(jì)計(jì)算表面積。表面積計(jì)算的誤差直接影響到參數(shù)計(jì)算結(jié)果；② 圈定ROI時(shí)有部分骨組織被截?cái)?，形成邊際效應(yīng)影響測(cè)量結(jié)果。

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

4 結(jié)論

本測(cè)量系統(tǒng)能夠?qū)icroCT斷層掃描獲得的動(dòng)物骨小梁圖像進(jìn)行骨密度、骨小梁厚度、骨結(jié)構(gòu)指數(shù)等參數(shù)進(jìn)行全面準(zhǔn)確的測(cè)量，并提供骨小梁結(jié)構(gòu)的三維可視化效果，為臨床研究動(dòng)物骨質(zhì)疏松程度、骨健康情況提供了一個(gè)更全面的參考工具。

表1 豬頜骨各項(xiàng)參數(shù)測(cè)量結(jié)果

[1]高鵬,戎軍艷,廖琪梅,等.MicroCT系統(tǒng)軟件平臺(tái)的設(shè)計(jì)與實(shí)現(xiàn)[J].醫(yī)療衛(wèi)生裝備,2014,35(12):22-24.

[2]付鑫,馬劍雄,董寶康.骨微結(jié)構(gòu)檢測(cè)方法研究進(jìn)展[J].中國(guó)骨與關(guān)節(jié)外科,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

國(guó)家自然科學(xué)基金（61127002，61179035）；中央高?；緲I(yè)務(wù)費(fèi)重點(diǎn)項(xiàng)目培育計(jì)劃（021414310017）。

作者郵箱：owenwu_mx@163.com