基于锥形束CT数据估算颌骨骨量变化的方法研究*

李伟志 王瑞永 刘 敏 顾晓宇 乔惠婷

颌骨的骨量变化，尤其是牙槽突骨量减少，是口腔医学临床工作中常常要面对的问题。牙周病的进展和治疗会使牙槽突产生不同程度的骨量变化[1]。牙拔除术不可避免地导致牙槽骨骨量减少[2]，即使在牙齿拔除（或脱落）后自然愈合期牙槽突骨量减少还会进行性发展[3，4]。有效定量定位地评估颌骨骨量变化程度对于评价牙拔除术对牙槽骨的损伤程度、牙周病进程及正畸正颌美容整形等手术效果有重要意义[5]。近年来微创拔牙相关手术技术不断发展，旨在减少对牙槽骨的损伤[6]；组织工程技术的应用（比如牙种植或牙周相关植骨手术）在一定程度上可以修复牙槽突骨量减少的情况[7]，牙槽突骨量变化的准确定位和对应体积的有效估算对于口腔临床工作中新技术及新方法有重要的支撑意义[8]。

借助锥形束CT（cone-beam computed tomography，CBCT）可以无创地进行颌骨（含牙槽突）的骨量探测[9]。目前对于颌骨牙槽突骨量变化的检测，主要通过人工设定标志点并进行测量完成[10，11]，这种方法对成像方位有严格要求[12]，得到的测量值仅为部分二维平面特征，无法对三维空间内颌骨牙槽骨的骨量变化进行全面评估，并且所获得颌骨牙槽骨特征值受阅片人为因素影响较大。Linderup 等人[13]提出基于两次CBCT 成像，分别建立颌面部三维模型，通过比较三维模型体积的变化来进行颌骨牙槽骨-骨量变化体积估算。由于涉及颌面部的三维建模工作量较大，不易在临床使用，因此临床急需自动或半自动的颌骨牙槽骨骨量变化的估算方法，来实现不同时期颌骨牙槽骨骨量变化情况的定量评估。

本文提出了一种基于两次CBCT 成像数据进行半自动颌骨牙槽骨骨量变化定量定位估算的方法，通过在不同时间点（如牙拔除术前后）分别进行CBCT 成像，以初始状态时刻的CBCT 图像作为基准，利用最大化互信息配准法，人工结合计算机程序化算法半自动将骨量变化估计时刻的CBCT 图像与基准进行配准，对比两幅已配准的CBCT 图像中颌骨差异，分离存在差异的体素，并根据成像参数估算两次成像之间的颌骨牙槽骨的骨量变化。本文利用平台实验验证该方法的可行性，并结合光学扫描实验对该方法进行了量化校正，证明该方法可以明确定位颌骨骨量变化位置并实现颌骨骨量变化体积的准确估算。

1.方法描述

1.1 颌骨牙槽骨骨量（体积）变化估算方法流程

估算从t1 到t2 一段时间内颌骨牙槽骨骨量变化的方法，需要实施2 次CBCT 成像，在配准两幅三维图像后，自动实现图像的差异提取及体素估算。首先，要在初始状态t1 时刻对颌面部进行CBCT 成像，得到参考图像，而后在拟进行颌骨牙槽骨骨量变化估算的t2 时刻再次对颌面部进行CBCT 成像，得到浮动图像。基于骨骼的刚性特点，配准算法选用刚体变换。估算刚体变换参数，把变换参数的初始值设为零。配准过程中利用估算参数对浮动图像进行变换。选用三线性插值方法进行体素插值，使参考图像与浮动图像具有相同的空间分辨率。定义互信息为两图像的相似性测度，用它衡量两幅图像配准后的匹配程度。采用方向加速法，迭代执行上述步骤，通过调整几何变换参数，使浮动图像和参考图像的互信息测度最大。配准完成后，参考图像和浮动图像做差获得差异图像。针对差异图像进行阈值分割获得颌骨牙槽骨变化量，完成一段时间内骨量实际变化的定量估算。

图1 估算颌骨牙槽骨骨量变化的方法流程图

1.2 基于最大互信息的图像配准 参考图像与浮动图像的配准是自动或半自动实现缺损体积估计的重要基础，可以克服两次独立成像体位不同等成像因素的影响，实现影像数据的有效对比，具体算法流程如图2 所示。互信息是一个统计相关性的测度，通常用于表示两个随机变量的依赖程度，也可以表示两幅图像的匹配程度，两幅图像的互信息值最大时，配准精度最高[14]。参考图像A 与浮动图像B的互信息MI(A，B)定义如下：

其中，PA(a)和PB(b)分别是参考图像A 和浮动图像B 的灰度边缘概率分布，PAB(a，b)为参考图像A 和浮动图像B 的联合概率分布。直接计算两幅图像的联合概率分布比较困难，可以借助于两幅图像的联合直方图h(a，b)获得[15]。

图像的精细配准过程是基于图像互信息的目标函数优化问题，其目标是找到使浮动图与参考图像对齐的空间映射，即变换矩阵T（包含旋转信息与位移信息）。变换矩阵可以实现图像在三维空间的方位变换，本算法采用Powell 搜索求解最优化变换矩阵。三维空间优化过程是由4 次一维搜索组成，分别沿x，y，z 三个正交方向分别搜索得到单向最优值x’，y’和z’，再沿初始点与单向最优值对应点（x’，y’，z’）连线方法进行一维搜索得到最优值，多次迭代后可得到最优解。

图2 基于最大互信息的配准算法流程图

1.3 差异图像的分割 将配准后的浮动图像和参考图像的对应体素作差后，获得差异图像。差异图像可以去除牙槽骨周围组织的影响，突出牙槽骨骨量变化。针对差异图像再进行阈值分割，获得两次成像牙槽骨变化的区域。阈值分割是一种直接检测区域的分割方法，它假设目标内相邻像素间的灰度值是相似的，但目标与背景的像素在灰度值上有差异。在灰度直方图上，目标和背景对应不同的峰值，选取的阈值位于两个峰谷间，从而将目标和背景分开。由于受不同成像因素的影响两次成像配准后并非完全相同，做差后得到的差异图像会体现出成像差异，不同于组织无变化区域的噪声数据，牙槽骨骨量变化区域会呈现明显数值，因此对于配准后作差得到的差异图像进行阈值分割可以选出三维图像中骨量变化的体素。

1.4 颌骨牙槽骨变化量估算 在实现骨量变化区域进行精确分割后，针对已标记的变化体素，使用如下公式进行计算可以得到骨量变化的实际体积。颌骨牙槽骨的骨量变化体积的计算公式为：

其中：V 为估算出的骨量变化体积，N 表示该区域的断层数目；Ai表示第i 断层内-发生骨量变化的区域的面积；Δh表示CBCT 重建数据的断层层厚，Mi为对差异图像分割后第i 层所分割出的像素数，a0为每一断层内单位像素所对应的面积，M 为对差异图像分割后得到骨量变化区域的体素数量，V0为三维影像中单位体素所对应的面积。该实验中所用CBCT 重建层厚0.25 mm，每一断层内像素对应面积0.0625 mm2，每个体素体积为0.0156 mm3。

2.方法的实验验证

2.1 实验对象 本文利用人类湿下颌骨标本及密度与人体牙槽骨接近的猪肋骨块进行实验，其中涉及1套人类湿下颌骨标本，该下颌标本存在10颗完整牙齿，5 颗残根及残冠，如图3（1）所示，猪骨块共8 块，取自成年猪第4 肋，骨块由小到大排列，如图3（2）所示。猪肋骨块与人体骨松质具有相近的成像特性，而发生在骨松质的下颌骨骨量改变非常不容易被察觉、辨别，因此本文通过在下颌骨标本右下颌第二磨牙远中的磨牙后垫所对应的骨组织表面放置猪肋骨块模拟多种下颌骨骨量变化情况。本研究获得北京市海淀医院生物伦理委员会审查批准（批准号2019022）。

图3 CBCT成像实验所用材料：（1）人体湿下颌骨，（2）猪骨块

2.2 实验过程 使用口腔锥形束CT（北京朗视仪器有限公司），对人类湿下颌骨标本及猪肋骨块进行成像，其中对人类湿下颌骨标本单独进行成像，并将8块猪骨块分别放置在人类湿下颌骨标本上，分别进行全方位扫描成像。所有研究对象扫描参数设定为：电流4 mA，电压100 kV，将扫描数据用锥形束CT自带的软件进行三维图像重建，重建后图像以DICOM格式存储，共得到三维锥形束CT图像9套，其中1套为独立人类湿下颌骨三维影像，另8套分别为带有不同体积的猪骨块的人类湿下颌骨三维影像，各体素均为0.25 mm×0.25 mm×0.25 mm。

获得锥形束CT三维图像后，基于VS2017平台C++语言，在主频2.2 GHz，Intel-Xeon E5-2650 v4处理器，64GByte内存环境下运行算法程序，分别对8组对照图像进行配准、做差、分割后实现下颌骨骨量变化体积的估算。

鉴于CBCT存在容积效应，基于CBCT图像通过算法的颌骨骨量缺损的估算值可能会出现少量偏差。本文使用光学三维扫描仪（Breuckmann precision in 3D），逐一对猪肋骨进行光学扫描得到猪肋骨块的实际体积。并利用MATLAB 软件，以光学扫描得到的体积为真值，通过一元线性回归的方法对基于CBCT 的颌骨骨量缺损算法进行定量线性校正，为确保在颌骨骨量无变化时能准确估算，回归方程常数项设为0，回归系数即为算法的校正系数。

3.结果

本文分别对人体湿下颌骨标本以及放置体积不同的猪肋骨块后的下颌骨标本进行CBCT 成像，共得到9套CBCT 图像，形成8组对照图像（放置猪骨块的下颌骨标本与无猪骨块下颌骨标本的2 次对照成像），图4 是在使用CBCT 常用设置，未进行人工处理条件下得到CBCT 断层图像，这是一组对照图像中的某一断层，其中红色圈出猪肋骨块放置的位置。

图4 人体湿下颌骨及猪骨块CBCT成像

表1 基于图像的骨量变化估算值与光学扫描值（mm3）

成像后，分别对8 组对照图像进行配准、做差、分割后实现下颌骨骨量变化体积的估算，得到的骨量变化数据（见表1）。算法能将人眼不易察觉的骨量变化情况清晰地显示出来，如图5 所示，图5（1）、（2）和（3）是利用CBCT（朗视）采用Ray-casting方法进行体绘制所得到的三维影像，其中图5（1）、（2）均是下颌骨上放置猪肋骨块后所得到的三维影像，在口腔科CBCT 常规设置下直接得到的三维图像图5（1）中并不能清晰地观测猪肋骨块，在人工手动调整显示窗位、窗宽后图5（2）中隐约可见猪肋骨块，图5（3）是无猪肋骨块的下颌骨独立成像所得到的对照三维影像，而图5（4）是利用算法估算出颌骨骨量变化情况的可视化显示，红色区域为骨量变化区域，与成像过程中所放置的猪骨块位置一致。

图5 CBCT图像三维可视化结果及算法标注骨量变化区域

通过光学扫描分别获得8 个大小不同猪骨块的实际体积（表1）作为骨量变化的真值对算法进行校正。所得到的校正系数为0.8044，校正后的确定系数R-square为0.988，校正结果如图6所示。

图6 骨量变化估算的定量校正结果

4.讨论

在临床中颌骨的骨量变化是难以直接测量的，因此颌骨的缺损与修复都难以进行定位定量分析。将同一对象的两次锥形束CT成像与图像处理算法相结合，应用于颌骨骨量变化的估算给口腔手术及治疗技术的评价带来了新的思路。本文提出基于两次锥形束CT 成像数据估算颌骨骨量的变化值，并形成算法程序，可以通过牙拔除术前后两次成像判断牙拔除术对下颌牙槽骨的即刻损伤，也可以对牙周、正畸或正颌手术等治疗后患者进行定期成像，跟踪治疗进展及此过程中颌骨骨量的长期变化，对于口腔科治疗技术效果的定量评价有重要意义，可以证明并推动新治疗技术的发展。

本文基于两次锥形束CT成像，提出了半自动化估算颌骨变化的方法的具体流程，其中两次独立成像的锥形束CT图像的配准是重要基础，本文利用最大互信息原理实现了图像的自动配准算法，并基于配准结果实现了做差、分割及骨量变化的定量定位估算算法。尽管锥形束CT成像受体位、具体操作设置及多种因素影响，任意两次成像都存在差异，但基于互信息的配准算法有效避免了成像体位差异的影响因素。配准后做差得到的差异图像大大降低了周围组织的影响，将骨量变化凸显出来。后续在差异图像的感兴趣区域进行阈值分割，通过阈值筛选去掉无变化区域的影响，将骨量变化的区域分割出来，通过在计算机上运行程序，导入两次成像得到的CBCT图像数据，可同时实现骨量变化的定位标记与定量估计。

为验证方法的有效性，本文设计了利用人体湿下颌骨及猪肋骨块的下颌骨改变仿真成像实验，其中利用猪骨块模拟下颌骨松质骨骨量变化。本文选取了下颌牙槽骨骨松质区域，与猪肋骨块进行区域体素CT灰度值比较，人工选取的牙槽骨骨松质区域与猪肋骨块的CT灰度值非常接近，平均相对差异为1.1%。在此方法中没有其他人体组织CT灰度值介于猪骨和人骨之间对计算产生干扰，因此利用猪肋骨块进行下颌骨骨量变化成像实验是有效并且可行的。

基于8 组成像实验，对所提出方法进行了校验。对于多次独立成像所得到了锥形束CT颌骨图像，本文所提出的基于最大互信息的配准方法可以在不做特殊标记的前提下自动实现图像的配准，并估算出颌骨骨量变化并显示出相关变化的空间对应位置。随着成像所用猪肋骨大小的变化，算法能很好地估算出下颌骨骨量变化情况。利用光学扫描猪肋骨块得到每组成像对象的真实骨量差异，对该方法的估算结果进行比对，发现一致性很好，随着真实骨量差异增大，算法估算结果随之增大，并且算法估算结果与真实值具有很好的线性关系。这主要是由于CT 成像存在容积效应，基于CT 图像所估算的体积值会与真实体积有偏差。本文利用光学扫描结果对算法估算值进行修正，得到修正系数0.8044，修正后确定系数为0.988。确定系数r-square 的正常取值范围为[0，1]，其越接近1 表面估算结果与真值越接近，因此校正后的算法可以对0-500 mm3间骨量变化进行有效估算。

本文提出的方法是基于锥形束CT图像进行骨量变化的估算，估算方法的灵敏度受限于成像体素大小及算法重采样体素大小，该方法对于骨量变化估算的灵敏度为0.0156 mm3，高于人眼的灵敏度。该方法不仅可以定量估算出颌骨骨量变化，还可以标记出骨量变化的位置。如图4中红色显示区域，在算法中具有做差分割功能，可以将骨量变化区域单独分割出来，若与原图像融合即可清晰地显示颌骨骨量变化的区域，这对于口腔科手术及治疗的相关研究有重要价值。定量的骨量变化估算是评价手术等治疗手段效果的基本指标之一，而同时明确骨量变化的位置对于技术手段的改进有积极的指导意义。