基于錐形束CT數(shù)據(jù)估算頜骨骨量變化的方法研究*

李偉志 王瑞永 劉 敏 顧曉宇 喬惠婷

頜骨的骨量變化，尤其是牙槽突骨量減少，是口腔醫(yī)學(xué)臨床工作中常常要面對的問題。牙周病的進(jìn)展和治療會使牙槽突產(chǎn)生不同程度的骨量變化[1]。牙拔除術(shù)不可避免地導(dǎo)致牙槽骨骨量減少[2]，即使在牙齒拔除（或脫落）后自然愈合期牙槽突骨量減少還會進(jìn)行性發(fā)展[3，4]。有效定量定位地評估頜骨骨量變化程度對于評價牙拔除術(shù)對牙槽骨的損傷程度、牙周病進(jìn)程及正畸正頜美容整形等手術(shù)效果有重要意義[5]。近年來微創(chuàng)拔牙相關(guān)手術(shù)技術(shù)不斷發(fā)展，旨在減少對牙槽骨的損傷[6]；組織工程技術(shù)的應(yīng)用（比如牙種植或牙周相關(guān)植骨手術(shù)）在一定程度上可以修復(fù)牙槽突骨量減少的情況[7]，牙槽突骨量變化的準(zhǔn)確定位和對應(yīng)體積的有效估算對于口腔臨床工作中新技術(shù)及新方法有重要的支撐意義[8]。

借助錐形束CT（cone-beam computed tomography，CBCT）可以無創(chuàng)地進(jìn)行頜骨（含牙槽突）的骨量探測[9]。目前對于頜骨牙槽突骨量變化的檢測，主要通過人工設(shè)定標(biāo)志點(diǎn)并進(jìn)行測量完成[10，11]，這種方法對成像方位有嚴(yán)格要求[12]，得到的測量值僅為部分二維平面特征，無法對三維空間內(nèi)頜骨牙槽骨的骨量變化進(jìn)行全面評估，并且所獲得頜骨牙槽骨特征值受閱片人為因素影響較大。Linderup 等人[13]提出基于兩次CBCT 成像，分別建立頜面部三維模型，通過比較三維模型體積的變化來進(jìn)行頜骨牙槽骨-骨量變化體積估算。由于涉及頜面部的三維建模工作量較大，不易在臨床使用，因此臨床急需自動或半自動的頜骨牙槽骨骨量變化的估算方法，來實(shí)現(xiàn)不同時期頜骨牙槽骨骨量變化情況的定量評估。

本文提出了一種基于兩次CBCT 成像數(shù)據(jù)進(jìn)行半自動頜骨牙槽骨骨量變化定量定位估算的方法，通過在不同時間點(diǎn)（如牙拔除術(shù)前后）分別進(jìn)行CBCT 成像，以初始狀態(tài)時刻的CBCT 圖像作為基準(zhǔn)，利用最大化互信息配準(zhǔn)法，人工結(jié)合計算機(jī)程序化算法半自動將骨量變化估計時刻的CBCT 圖像與基準(zhǔn)進(jìn)行配準(zhǔn)，對比兩幅已配準(zhǔn)的CBCT 圖像中頜骨差異，分離存在差異的體素，并根據(jù)成像參數(shù)估算兩次成像之間的頜骨牙槽骨的骨量變化。本文利用平臺實(shí)驗(yàn)驗(yàn)證該方法的可行性，并結(jié)合光學(xué)掃描實(shí)驗(yàn)對該方法進(jìn)行了量化校正，證明該方法可以明確定位頜骨骨量變化位置并實(shí)現(xiàn)頜骨骨量變化體積的準(zhǔn)確估算。

1.方法描述

1.1 頜骨牙槽骨骨量（體積）變化估算方法流程

估算從t1 到t2 一段時間內(nèi)頜骨牙槽骨骨量變化的方法，需要實(shí)施2 次CBCT 成像，在配準(zhǔn)兩幅三維圖像后，自動實(shí)現(xiàn)圖像的差異提取及體素估算。首先，要在初始狀態(tài)t1 時刻對頜面部進(jìn)行CBCT 成像，得到參考圖像，而后在擬進(jìn)行頜骨牙槽骨骨量變化估算的t2 時刻再次對頜面部進(jìn)行CBCT 成像，得到浮動圖像。基于骨骼的剛性特點(diǎn)，配準(zhǔn)算法選用剛體變換。估算剛體變換參數(shù)，把變換參數(shù)的初始值設(shè)為零。配準(zhǔn)過程中利用估算參數(shù)對浮動圖像進(jìn)行變換。選用三線性插值方法進(jìn)行體素插值，使參考圖像與浮動圖像具有相同的空間分辨率。定義互信息為兩圖像的相似性測度，用它衡量兩幅圖像配準(zhǔn)后的匹配程度。采用方向加速法，迭代執(zhí)行上述步驟，通過調(diào)整幾何變換參數(shù)，使浮動圖像和參考圖像的互信息測度最大。配準(zhǔn)完成后，參考圖像和浮動圖像做差獲得差異圖像。針對差異圖像進(jìn)行閾值分割獲得頜骨牙槽骨變化量，完成一段時間內(nèi)骨量實(shí)際變化的定量估算。

圖1 估算頜骨牙槽骨骨量變化的方法流程圖

1.2 基于最大互信息的圖像配準(zhǔn) 參考圖像與浮動圖像的配準(zhǔn)是自動或半自動實(shí)現(xiàn)缺損體積估計的重要基礎(chǔ)，可以克服兩次獨(dú)立成像體位不同等成像因素的影響，實(shí)現(xiàn)影像數(shù)據(jù)的有效對比，具體算法流程如圖2 所示?；バ畔⑹且粋€統(tǒng)計相關(guān)性的測度，通常用于表示兩個隨機(jī)變量的依賴程度，也可以表示兩幅圖像的匹配程度，兩幅圖像的互信息值最大時，配準(zhǔn)精度最高[14]。參考圖像A 與浮動圖像B的互信息MI(A，B)定義如下：

其中，PA(a)和PB(b)分別是參考圖像A 和浮動圖像B 的灰度邊緣概率分布，PAB(a，b)為參考圖像A 和浮動圖像B 的聯(lián)合概率分布。直接計算兩幅圖像的聯(lián)合概率分布比較困難，可以借助于兩幅圖像的聯(lián)合直方圖h(a，b)獲得[15]。

圖像的精細(xì)配準(zhǔn)過程是基于圖像互信息的目標(biāo)函數(shù)優(yōu)化問題，其目標(biāo)是找到使浮動圖與參考圖像對齊的空間映射，即變換矩陣T（包含旋轉(zhuǎn)信息與位移信息）。變換矩陣可以實(shí)現(xiàn)圖像在三維空間的方位變換，本算法采用Powell 搜索求解最優(yōu)化變換矩陣。三維空間優(yōu)化過程是由4 次一維搜索組成，分別沿x，y，z 三個正交方向分別搜索得到單向最優(yōu)值x’，y’和z’，再沿初始點(diǎn)與單向最優(yōu)值對應(yīng)點(diǎn)（x’，y’，z’）連線方法進(jìn)行一維搜索得到最優(yōu)值，多次迭代后可得到最優(yōu)解。

圖2 基于最大互信息的配準(zhǔn)算法流程圖

1.3 差異圖像的分割 將配準(zhǔn)后的浮動圖像和參考圖像的對應(yīng)體素作差后，獲得差異圖像。差異圖像可以去除牙槽骨周圍組織的影響，突出牙槽骨骨量變化。針對差異圖像再進(jìn)行閾值分割，獲得兩次成像牙槽骨變化的區(qū)域。閾值分割是一種直接檢測區(qū)域的分割方法，它假設(shè)目標(biāo)內(nèi)相鄰像素間的灰度值是相似的，但目標(biāo)與背景的像素在灰度值上有差異。在灰度直方圖上，目標(biāo)和背景對應(yīng)不同的峰值，選取的閾值位于兩個峰谷間，從而將目標(biāo)和背景分開。由于受不同成像因素的影響兩次成像配準(zhǔn)后并非完全相同，做差后得到的差異圖像會體現(xiàn)出成像差異，不同于組織無變化區(qū)域的噪聲數(shù)據(jù)，牙槽骨骨量變化區(qū)域會呈現(xiàn)明顯數(shù)值，因此對于配準(zhǔn)后作差得到的差異圖像進(jìn)行閾值分割可以選出三維圖像中骨量變化的體素。

1.4 頜骨牙槽骨變化量估算 在實(shí)現(xiàn)骨量變化區(qū)域進(jìn)行精確分割后，針對已標(biāo)記的變化體素，使用如下公式進(jìn)行計算可以得到骨量變化的實(shí)際體積。頜骨牙槽骨的骨量變化體積的計算公式為：

其中：V 為估算出的骨量變化體積，N 表示該區(qū)域的斷層數(shù)目；Ai表示第i 斷層內(nèi)-發(fā)生骨量變化的區(qū)域的面積；Δh表示CBCT 重建數(shù)據(jù)的斷層層厚，Mi為對差異圖像分割后第i 層所分割出的像素數(shù)，a0為每一斷層內(nèi)單位像素所對應(yīng)的面積，M 為對差異圖像分割后得到骨量變化區(qū)域的體素數(shù)量，V0為三維影像中單位體素所對應(yīng)的面積。該實(shí)驗(yàn)中所用CBCT 重建層厚0.25 mm，每一斷層內(nèi)像素對應(yīng)面積0.0625 mm2，每個體素體積為0.0156 mm3。

2.方法的實(shí)驗(yàn)驗(yàn)證

2.1 實(shí)驗(yàn)對象 本文利用人類濕下頜骨標(biāo)本及密度與人體牙槽骨接近的豬肋骨塊進(jìn)行實(shí)驗(yàn)，其中涉及1套人類濕下頜骨標(biāo)本，該下頜標(biāo)本存在10顆完整牙齒，5 顆殘根及殘冠，如圖3（1）所示，豬骨塊共8 塊，取自成年豬第4 肋，骨塊由小到大排列，如圖3（2）所示。豬肋骨塊與人體骨松質(zhì)具有相近的成像特性，而發(fā)生在骨松質(zhì)的下頜骨骨量改變非常不容易被察覺、辨別，因此本文通過在下頜骨標(biāo)本右下頜第二磨牙遠(yuǎn)中的磨牙后墊所對應(yīng)的骨組織表面放置豬肋骨塊模擬多種下頜骨骨量變化情況。本研究獲得北京市海淀醫(yī)院生物倫理委員會審查批準(zhǔn)（批準(zhǔn)號2019022）。

圖3 CBCT成像實(shí)驗(yàn)所用材料：（1）人體濕下頜骨，（2）豬骨塊

2.2 實(shí)驗(yàn)過程 使用口腔錐形束CT（北京朗視儀器有限公司），對人類濕下頜骨標(biāo)本及豬肋骨塊進(jìn)行成像，其中對人類濕下頜骨標(biāo)本單獨(dú)進(jìn)行成像，并將8塊豬骨塊分別放置在人類濕下頜骨標(biāo)本上，分別進(jìn)行全方位掃描成像。所有研究對象掃描參數(shù)設(shè)定為：電流4 mA，電壓100 kV，將掃描數(shù)據(jù)用錐形束CT自帶的軟件進(jìn)行三維圖像重建，重建后圖像以DICOM格式存儲，共得到三維錐形束CT圖像9套，其中1套為獨(dú)立人類濕下頜骨三維影像，另8套分別為帶有不同體積的豬骨塊的人類濕下頜骨三維影像，各體素均為0.25 mm×0.25 mm×0.25 mm。

獲得錐形束CT三維圖像后，基于VS2017平臺C++語言，在主頻2.2 GHz，Intel-Xeon E5-2650 v4處理器，64GByte內(nèi)存環(huán)境下運(yùn)行算法程序，分別對8組對照圖像進(jìn)行配準(zhǔn)、做差、分割后實(shí)現(xiàn)下頜骨骨量變化體積的估算。

鑒于CBCT存在容積效應(yīng)，基于CBCT圖像通過算法的頜骨骨量缺損的估算值可能會出現(xiàn)少量偏差。本文使用光學(xué)三維掃描儀（Breuckmann precision in 3D），逐一對豬肋骨進(jìn)行光學(xué)掃描得到豬肋骨塊的實(shí)際體積。并利用MATLAB 軟件，以光學(xué)掃描得到的體積為真值，通過一元線性回歸的方法對基于CBCT 的頜骨骨量缺損算法進(jìn)行定量線性校正，為確保在頜骨骨量無變化時能準(zhǔn)確估算，回歸方程常數(shù)項(xiàng)設(shè)為0，回歸系數(shù)即為算法的校正系數(shù)。

3.結(jié)果

本文分別對人體濕下頜骨標(biāo)本以及放置體積不同的豬肋骨塊后的下頜骨標(biāo)本進(jìn)行CBCT 成像，共得到9套CBCT 圖像，形成8組對照圖像（放置豬骨塊的下頜骨標(biāo)本與無豬骨塊下頜骨標(biāo)本的2 次對照成像），圖4 是在使用CBCT 常用設(shè)置，未進(jìn)行人工處理?xiàng)l件下得到CBCT 斷層圖像，這是一組對照圖像中的某一斷層，其中紅色圈出豬肋骨塊放置的位置。

圖4 人體濕下頜骨及豬骨塊CBCT成像

表1 基于圖像的骨量變化估算值與光學(xué)掃描值（mm3）

成像后，分別對8 組對照圖像進(jìn)行配準(zhǔn)、做差、分割后實(shí)現(xiàn)下頜骨骨量變化體積的估算，得到的骨量變化數(shù)據(jù)（見表1）。算法能將人眼不易察覺的骨量變化情況清晰地顯示出來，如圖5 所示，圖5（1）、（2）和（3）是利用CBCT（朗視）采用Ray-casting方法進(jìn)行體繪制所得到的三維影像，其中圖5（1）、（2）均是下頜骨上放置豬肋骨塊后所得到的三維影像，在口腔科CBCT 常規(guī)設(shè)置下直接得到的三維圖像圖5（1）中并不能清晰地觀測豬肋骨塊，在人工手動調(diào)整顯示窗位、窗寬后圖5（2）中隱約可見豬肋骨塊，圖5（3）是無豬肋骨塊的下頜骨獨(dú)立成像所得到的對照三維影像，而圖5（4）是利用算法估算出頜骨骨量變化情況的可視化顯示，紅色區(qū)域?yàn)楣橇孔兓瘏^(qū)域，與成像過程中所放置的豬骨塊位置一致。

圖5 CBCT圖像三維可視化結(jié)果及算法標(biāo)注骨量變化區(qū)域

通過光學(xué)掃描分別獲得8 個大小不同豬骨塊的實(shí)際體積（表1）作為骨量變化的真值對算法進(jìn)行校正。所得到的校正系數(shù)為0.8044，校正后的確定系數(shù)R-square為0.988，校正結(jié)果如圖6所示。

圖6 骨量變化估算的定量校正結(jié)果

4.討論

在臨床中頜骨的骨量變化是難以直接測量的，因此頜骨的缺損與修復(fù)都難以進(jìn)行定位定量分析。將同一對象的兩次錐形束CT成像與圖像處理算法相結(jié)合，應(yīng)用于頜骨骨量變化的估算給口腔手術(shù)及治療技術(shù)的評價帶來了新的思路。本文提出基于兩次錐形束CT 成像數(shù)據(jù)估算頜骨骨量的變化值，并形成算法程序，可以通過牙拔除術(shù)前后兩次成像判斷牙拔除術(shù)對下頜牙槽骨的即刻損傷，也可以對牙周、正畸或正頜手術(shù)等治療后患者進(jìn)行定期成像，跟蹤治療進(jìn)展及此過程中頜骨骨量的長期變化，對于口腔科治療技術(shù)效果的定量評價有重要意義，可以證明并推動新治療技術(shù)的發(fā)展。

本文基于兩次錐形束CT成像，提出了半自動化估算頜骨變化的方法的具體流程，其中兩次獨(dú)立成像的錐形束CT圖像的配準(zhǔn)是重要基礎(chǔ)，本文利用最大互信息原理實(shí)現(xiàn)了圖像的自動配準(zhǔn)算法，并基于配準(zhǔn)結(jié)果實(shí)現(xiàn)了做差、分割及骨量變化的定量定位估算算法。盡管錐形束CT成像受體位、具體操作設(shè)置及多種因素影響，任意兩次成像都存在差異，但基于互信息的配準(zhǔn)算法有效避免了成像體位差異的影響因素。配準(zhǔn)后做差得到的差異圖像大大降低了周圍組織的影響，將骨量變化凸顯出來。后續(xù)在差異圖像的感興趣區(qū)域進(jìn)行閾值分割，通過閾值篩選去掉無變化區(qū)域的影響，將骨量變化的區(qū)域分割出來，通過在計算機(jī)上運(yùn)行程序，導(dǎo)入兩次成像得到的CBCT圖像數(shù)據(jù)，可同時實(shí)現(xiàn)骨量變化的定位標(biāo)記與定量估計。

為驗(yàn)證方法的有效性，本文設(shè)計了利用人體濕下頜骨及豬肋骨塊的下頜骨改變仿真成像實(shí)驗(yàn)，其中利用豬骨塊模擬下頜骨松質(zhì)骨骨量變化。本文選取了下頜牙槽骨骨松質(zhì)區(qū)域，與豬肋骨塊進(jìn)行區(qū)域體素CT灰度值比較，人工選取的牙槽骨骨松質(zhì)區(qū)域與豬肋骨塊的CT灰度值非常接近，平均相對差異為1.1%。在此方法中沒有其他人體組織CT灰度值介于豬骨和人骨之間對計算產(chǎn)生干擾，因此利用豬肋骨塊進(jìn)行下頜骨骨量變化成像實(shí)驗(yàn)是有效并且可行的。

基于8 組成像實(shí)驗(yàn)，對所提出方法進(jìn)行了校驗(yàn)。對于多次獨(dú)立成像所得到了錐形束CT頜骨圖像，本文所提出的基于最大互信息的配準(zhǔn)方法可以在不做特殊標(biāo)記的前提下自動實(shí)現(xiàn)圖像的配準(zhǔn)，并估算出頜骨骨量變化并顯示出相關(guān)變化的空間對應(yīng)位置。隨著成像所用豬肋骨大小的變化，算法能很好地估算出下頜骨骨量變化情況。利用光學(xué)掃描豬肋骨塊得到每組成像對象的真實(shí)骨量差異，對該方法的估算結(jié)果進(jìn)行比對，發(fā)現(xiàn)一致性很好，隨著真實(shí)骨量差異增大，算法估算結(jié)果隨之增大，并且算法估算結(jié)果與真實(shí)值具有很好的線性關(guān)系。這主要是由于CT 成像存在容積效應(yīng)，基于CT 圖像所估算的體積值會與真實(shí)體積有偏差。本文利用光學(xué)掃描結(jié)果對算法估算值進(jìn)行修正，得到修正系數(shù)0.8044，修正后確定系數(shù)為0.988。確定系數(shù)r-square 的正常取值范圍為[0，1]，其越接近1 表面估算結(jié)果與真值越接近，因此校正后的算法可以對0-500 mm3間骨量變化進(jìn)行有效估算。

本文提出的方法是基于錐形束CT圖像進(jìn)行骨量變化的估算，估算方法的靈敏度受限于成像體素大小及算法重采樣體素大小，該方法對于骨量變化估算的靈敏度為0.0156 mm3，高于人眼的靈敏度。該方法不僅可以定量估算出頜骨骨量變化，還可以標(biāo)記出骨量變化的位置。如圖4中紅色顯示區(qū)域，在算法中具有做差分割功能，可以將骨量變化區(qū)域單獨(dú)分割出來，若與原圖像融合即可清晰地顯示頜骨骨量變化的區(qū)域，這對于口腔科手術(shù)及治療的相關(guān)研究有重要價值。定量的骨量變化估算是評價手術(shù)等治療手段效果的基本指標(biāo)之一，而同時明確骨量變化的位置對于技術(shù)手段的改進(jìn)有積極的指導(dǎo)意義。