CN114241173A - 一种牙齿cbct图像三维分割方法及系统 - Google Patents

Abstract

本发明公开了一种牙齿CBCT图像三维分割方法及系统，首先采用各向异性扩散方法对牙齿CBCT图像进行预处理操作，在去除CBCT图像噪声干扰的情况下，还可以保留图像边缘细节特征，克服了高斯滤波预处理的缺陷；然后使用阈值大律法与形态学相结合的方法，得到牙齿的感兴趣区域，并由此产生初始水平集函数；在此基础上，采用三维可变尺度区域拟合分割模型提取牙齿区域，使得牙齿分割结果的鲁棒性和准确性都明显提升。

Description

一种牙齿CBCT图像三维分割方法及系统

技术领域

本发明属于图像分割技术领域，具体涉及一种牙齿CBCT图像三维分割方法及系统的设计。

背景技术

随着医学技术的不断发展与推进，锥形束CT(Cone Beam Computer Tomography，CBCT)在口腔疾病的诊断与治疗中获得了广泛的应用。与传统的计算机断层扫描图像(CT)相比，CBCT具有更短的扫描时间、更低的辐射剂量和更高的图像空间分辨率。因此，目前CBCT已成为诊断牙齿疾病与牙齿正畸最为重要和关键的检查工具，也是牙齿疾病诊疗前必做的检查项目。

对牙齿CBCT图像中的每个牙齿进行精确自动地分割，是牙齿疾病诊断与牙齿矫正等临床决策与制定治疗计划的关键。每个牙齿在结构上包括牙釉质、牙本质、牙髓等部分，然而，在CBCT图像中，牙齿与牙槽骨的灰度值相近，相邻牙齿相互粘连、边界不明显，而且牙齿又分为单牙根(如尖牙、切牙)以及多牙根(磨牙、第二磨牙等)等多种形状结构，这些问题对精确分割单颗牙齿带来很大困难。此外，牙釉质、牙本质以及牙髓的结构较小，且相互间的灰度不均匀性也给分割牙齿带来了极大的挑战。

牙齿CBCT图像分割算法利用人工智能技术在医学图像数据分析中的优势，能够自动或半自动的实现牙齿图像分割的过程，从而提高口腔疾病的诊疗效率、降低诊疗难度，在口腔正畸(各种牙列不齐的矫正等)、口腔外科(种植牙、拔牙等)和口腔修复(铸造金属全冠、烤瓷全冠等)在内的相关口腔疾病的临床诊疗过程中存在重要的临床需求，能够有效的缓解我国口腔医疗人力资源存在的持续紧张的态势。因此，牙齿CBCT图像分割算法的研究具有重要的科学意义和应用价值。精确且完整地分割三维牙齿图像数据，在牙齿正畸中起着至关重要的作用，尤其是在自动排牙以及牙齿生物力学仿生的口腔临床应用中尤为关键。

近年来，国内外学者提出了各种牙齿CBCT图像分割方法，包括对健康牙齿图像的自动分割，以及不健康牙齿和人工牙图像的分割。不健康的牙齿有错位牙、阻生牙和埋伏牙，人工牙有种植牙和活动假牙等。现有文献中的牙齿CBCT图像分割方法主要是水平集、深度学习、区域增长、图割、B样条曲线拟合等。目前基于水平集的方法大多使用二维逐层分割的方式，如在基于边缘信息的水平集方法中，引入了灰度先验能量项，使得每层牙齿图像的分割都可以使用上一层分割结果的灰度先验信息，这种方法有助于多牙根牙齿图像的分割，也可以引入保持凸性的能量项，利用单牙根牙齿的解剖学形状上凸的特点，实现了对切牙和尖牙的高精度分割。以上方法主要使用二维的水平集逐层分割牙齿图像，这种做法没有利用牙齿的三维几何性质，而且逐层分割容易在跟踪过程中产生误差。在牙齿CBCT图像分割的早期研究中，通常采用基于CV模型实现了健康牙齿图像的三维水平集分割，但该方法分割精度较低，鲁棒性较差，主要因为CV模型在对灰度不均匀且背景复杂的CBCT牙齿图像的分割上表现不佳。此外，还有区域增长与Canny算子相结合的方法，实现磨牙图像的高精度分割，但该方法容易出现过分割现象，且方法的鲁棒性较差。还有基于B样条曲线拟合的方法对牙齿图像进行逐层分割，但该方法对手动初始化的要求较高。

以上二维逐层分割的方法的鲁棒性通常表现较低，因为每颗牙齿切面图像的分割中均使用了待分割牙齿和相邻牙齿在其它层的分割结果，这会导致分割中出现的任何误差都有可能对所有牙齿的精确分割造成严重地影响。

发明内容

本发明的目的是为了解决现有二维逐层分割的牙齿CBCT图像分割方法分割精度较低且鲁棒性较差的问题，提出了一种牙齿CBCT图像三维分割方法及系统。

本发明的技术方案为：一种牙齿CBCT图像三维分割方法，包括以下步骤：

S1、获取待分割的牙齿CBCT图像I₀。

S2、采用各向异性扩散方法对牙齿CBCT图像I₀进行预处理，得到预处理图像。

S3、采用自动阈值法从预处理图像中提取出牙齿区域I。

S4、根据牙齿区域I，采用三维可变尺度区域拟合分割模型得到牙齿三维分割结果。

进一步地，步骤S3包括以下分步骤：

S31、设置阈值t₁去除预处理图像中的背景区域，得到去背景图像。

S32、计算去背景图像的灰度直方图，通过阈值大律法，在灰度直方图的最大化类间方差和最小化类间方差之间计算得到一个最佳阈值t₂，通过最佳阈值t₂分割去背景图像中的牙齿区域I₁和软组织区域I₂。

S33、采用形态学方法中的闭操作对牙齿区域I₁进行补洞操作，得到牙齿区域I。

进一步地，步骤S31中阈值t₁的计算公式为：

t₁＝maxT×2％

其中maxT为预处理图像的最大灰度值。

进一步地，步骤S4包括以下分步骤：

S41、假设牙齿区域I表示为I:Ω→R^d，其中Ω表示图像域，d表示图像维数且d≥1。

S42、设曲线C为图像域Ω中的一条闭合曲线，将整个牙齿区域I分为内部区域Ω₁＝inside(C)和外部区域Ω₂＝outside(C)。

S43、定义三维空间中任意一点x∈Ω满足局部拟合能量函数E_x(C,f₁(x),f₂(x))，其中f₁(x)和f₂(x)分别表示在以x为圆心的小区域内Ω₁＝inside(C)和Ω₂＝outside(C)的各种局部区域内的近似拟合值。

S44、对局部拟合能量函数E_x(C,f₁(x),f₂(x))进行积分，得到能量函数E(C,f₁(x),f₂(x))。

S45、将能量函数E(C,f₁(x),f₂(x))转化为等价水平集函数

其中

表示水平集函数。

S46、在等价水平集函数

中引入距离正则项

和长度项

得到完整水平集形式的能量泛函

S47、将能量泛函

的极小化形式通过梯度下降流表示，得到水平集函数

的0-水平集曲线。

S48、不断迭代更新能量泛函

的极小化形式，直到水平集函数

的0-水平集曲线到达真实的牙齿边界处，得到牙齿三维分割结果。

进一步地，步骤S43中的局部拟合能量函数E_x(C,f₁(x),f₂(x))具体为：

E_x(C,f₁(x),f₂(x))＝λ₁∫_inside(C)K(x-y)|I(y)-f₁(x)|²dy+λ₂∫_outside(C)K(x-y)|I(y)-f₂(x)|²dy

其中λ₁,λ₂均为正常数，K(·)表示核函数，I(·)表示牙齿区域，y表示图像域Ω中的所有像素点。

进一步地，步骤S45中的等价水平集函数

具体为：

其中K_σ(·)表示高斯核函数，H(·)表示赫维赛德函数，y表示图像域Ω中的所有像素点。

进一步地，步骤S46中的距离正则项

和长度项

具体为：

其中δ(·)表示狄拉克函数，

表示水平集函数的梯度。

进一步地，步骤S46中的能量泛函

具体为：

其中μ,v均为非负常数。

进一步地，步骤S47中水平集函数

的0-水平集曲线表示为：

其中中间函数e₁(x)和e₂(x)分别表示为：

e₁(x)＝∫_ΩK_σ(y-x)|I(x)-f₁(y)|²dy

e₂(x)＝∫_ΩK_σ(y-x)|I(x)-f₂(y)|²dy

其中δ(·)表示狄拉克函数，

表示水平集函数的梯度，K_σ(·)表示高斯核函数，I(·)表示牙齿区域，y表示图像域Ω中的所有像素点，div(·)表示散度函数。

本发明还提供了一种牙齿CBCT图像三维分割系统，用于配置执行上述牙齿CBCT图像三维分割方法。

本发明的有益效果是：

(1)本发明采用各向异性扩散方法对牙齿CBCT图像进行预处理操作，在去除CBCT图像噪声干扰的情况下，还可以保留图像边缘细节特征，克服了高斯滤波预处理的缺陷。

(2)本发明使用阈值大律法与形态学相结合的方法，得到牙齿的感兴趣区域，并由此产生初始水平集函数，用于后续牙齿区域的准确分割。

(3)现有的二维逐层分割方法鲁棒性较低，最终的分割结果的准确性会由于每一层的误差累加而降低，本发明采用三维可变尺度区域拟合分割模型提取牙齿区域，可以避免上述问题的发生，牙齿分割结果的鲁棒性和准确性都明显提升。

附图说明

图1所示为本发明实施例提供的一种牙齿CBCT图像三维分割方法流程图。

具体实施方式

现在将参考附图来详细描述本发明的示例性实施方式。应当理解，附图中示出和描述的实施方式仅仅是示例性的，意在阐释本发明的原理和精神，而并非限制本发明的范围。

实施例一：

本发明实施例提供了一种牙齿CBCT图像三维分割方法，如图1所示，包括以下步骤S1～S4：

S1、获取待分割的牙齿CBCT图像I₀。

S2、采用各向异性扩散方法对牙齿CBCT图像I₀进行预处理，得到预处理图像。

各向异性扩散方法克服了高斯滤波的缺陷，可以在保留图像边缘细节特征的同时去除噪声的干扰，平滑图像，增强图像边缘的对比度。

S3、采用自动阈值法从预处理图像中提取出牙齿区域I。

步骤S3包括以下分步骤S31～S33：

S31、设置阈值t₁去除预处理图像中的背景区域，得到去背景图像。

牙齿CBCT图像主要由牙齿区域、软组织区域以及背景区域构成，背景区域一般对应空气或水，本发明实施例中，设置阈值t₁＝maxT×2％用于去除背景区域，其中maxT为预处理图像的最大灰度值。

S32、计算去背景图像的灰度直方图，通过阈值大律法，在灰度直方图的最大化类间方差和最小化类间方差之间计算得到一个最佳阈值t₂，通过最佳阈值t₂分割去背景图像中的牙齿区域I₁和软组织区域I₂。

S33、采用形态学方法中的闭操作对牙齿区域I₁进行补洞操作，得到牙齿区域I。

由于牙齿内部的牙髓区域的灰度值较低，使用步骤S31～S32时，这一部分会被视为背景区域，因此会造成牙齿区域出现空洞，在这种情况下，需要进行补洞操作。

S4、根据牙齿区域I，采用三维可变尺度区域拟合分割模型得到牙齿三维分割结果。

步骤S4包括以下分步骤S41～S48：

S41、假设牙齿区域I表示为I:Ω→R^d，其中Ω表示图像域，d表示图像维数且d≥1。

S42、设曲线C为图像域Ω中的一条闭合曲线，将整个牙齿区域I分为内部区域Ω₁＝inside(C)和外部区域Ω₂＝outside(C)。

S43、定义三维空间中任意一点x∈Ω满足局部拟合能量函数E_x(C,f₁(x),f₂(x))：

E_x(C,f₁(x),f₂(x))＝λ₁∫_inside(C)K(x-y)|I(y)-f₁(x)|²dy+λ₂∫_outside(C)K(x-y)|I(y)-f₂(x)|²dy

其中λ₁,λ₂均为正常数，K(·)表示核函数，I(·)表示牙齿区域，y表示图像域Ω中的所有像素点，f₁(x)和f₂(x)分别表示在以x为圆心的小区域内Ω₁＝inside(C)和Ω₂＝outside(C)的各种局部区域内的近似拟合值。

三维可变尺度区域拟合分割模型中的一个关键因素是具有局部化特征的核函数K(·)，该函数需满足如下条件：

(1)K(-u)＝K(u)；

(2)如果|u|<|v|，则K(u)≥K(v)，而且lim_u→+∞K(u)＝0；

(3)∫K(x)dx＝1。

核函数K(·)可以用于调节和控制f₁(x)和f₂(x)的大小，本发明实施例中使用高斯核函数。局部拟合能量函数E_x(C,f₁(x),f₂(x))的范围只是以x为中心的一个小邻域，因为f₁(x)和f₂(x)的值只是拟合x附近的图像灰度值，由于高斯核函数具有局部特性，服从正态分布，因此对f₁(x)和f₂(x)的值产生主要影响的是当y在x附近时，而当y远离x时，高斯核函数对f₁(x)和f₂(x)的值几乎不产生影响。

S44、对局部拟合能量函数E_x(C,f₁(x),f₂(x))进行积分，得到能量函数E(C,f₁(x),f₂(x))。

为了找到完整图像上的目标边界，必须找到一个轮廓曲线C使得能量函数在整个图像域Ω上的所有点x均为最小值，这个目标可以通过对局部拟合能量函数E_x(C,f₁(x),f₂(x))进行积分来实现，最终得到能量函数E(C,f₁(x),f₂(x))：

E(C,f₁(x),f₂(x))＝∫_ΩE_x(C,f₁(x),f₂(x))d_x

S45、为了方便处理CBCT图像中牙齿的拓扑结构变化，将能量函数E(C,f₁(x),f₂(x))转化为等价水平集函数

其中

表示水平集函数，K_σ(·)表示高斯核函数，H(·)表示赫维赛德函数，y表示图像域Ω中的所有像素点。

S46、为了保证水平集函数

可以稳定地演化，在等价水平集函数

中引入距离正则项

和长度项

得到完整水平集形式的能量泛函

本发明实施例中，距离正则项

和长度项

具体为：

其中δ(·)表示狄拉克函数，

表示水平集函数的梯度。

能量泛函

具体为：

其中μ,v均为非负常数。

S47、将能量泛函

的极小化形式通过梯度下降流表示，得到水平集函数

的0-水平集曲线。

本发明实施例中，通过对步骤S46中的能量泛函极小化来使得水平集函数

在目标边界处收敛停止迭代。对能量泛函的极小化的一般做法是使用标准下降法(或最速下降法)。对于一个给定的水平集函数

对f₁(x)和f₂(x)相对应的能量泛函进行极小化，通过变分法可得f₁(x)和f₂(x)的表达式为：

其中K_σ(·)表示高斯核函数，H_ε(·)表示隶属度函数。

保持f₁(x)和f₂(x)固定不变，能量泛函

的极小化形式通过梯度下降流表示为：

其中中间函数e₁(x)和e₂(x)分别表示为：

e₁(x)＝∫_ΩK_σ(y-x)|I(x)-f₁(y)|²dy

e₂(x)＝∫_ΩK_σ(y-x)|I(x)-f₂(y)|²dy

其中δ(·)表示狄拉克函数，

表示水平集函数的梯度，K_σ(·)表示高斯核函数，I(·)表示牙齿区域，y表示图像域Ω中的所有像素点，div(·)表示散度函数。

S48、不断迭代更新能量泛函

的极小化形式，直到水平集函数

的0-水平集曲线到达真实的牙齿边界处，得到牙齿三维分割结果。

实施例二：

本发明实施例提供了一种牙齿CBCT图像三维分割系统，用于配置执行实施例一中的牙齿CBCT图像三维分割方法。

本发明实施例中的系统可以为一种电子设备，包括存储器、处理器及存储在所述存储器上并在所述处理器上运行的计算机程序，所述处理器执行所述程序，以实现如实施例1中所述的牙齿CBCT图像三维分割方法的部分或全部步骤。

本发明实施例中，电子设备可以包括：处理器，存储器，总线和通信接口，处理器、通信接口和存储器通过总线连接，存储器中存储有可在处理器上运行的计算机程序，处理器运行该计算机程序时执行本申请前述实施例1所提供的牙齿CBCT图像三维分割方法的部分或全部步骤。

本发明实施例中的系统还可以为一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，执行所述计算机程序时实现如实施例1所述的牙齿CBCT图像三维分割方法的部分或全部步骤。

上述计算机可读存储介质可以中由任何类型的易失性或非易失性存储设备或它们的组合实现，如静态随机存取存储器(SRAM)，可擦除可维和只读存储器(EEPROM)，可擦除可编程只读存储器(EPROM)，可编程只读存储器(PROM)，只读存储器(ROM)，磁存储器，快闪存储器，磁盘或光盘，可读存储介质可以是通用或专用计算机能够存取的任何可用介质。可读存储介质耦合至处理器，从而使处理器能够从该可读存储介质读取信息，且可向该可读存储介质写入信息，可读存储介质也可以是处理器组成部分，处理器和可读存储介质可以位于专用集成电路(ASIC)中，处理器和可读存储介质也可以作为分立组件存在于点去配准系统中。

本发明实施例可提供为方法、装置或计算机程序产品，因此，本发明可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且本发明可采用一个或多个其中包括有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。参照根据本发明实施例的方法、设备(装置)、和计算机程序产品的流程和/或方框图来描述的，应理解可由计算机程序指令实现流程图和/或方框图中每一流程和/或方框图、以及流程图和/或方框图的结合，可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。这些计算机程序指令也可装到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图中一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (10)

1.一种牙齿CBCT图像三维分割方法，其特征在于，包括以下步骤：

S1、获取待分割的牙齿CBCT图像I₀；

S2、采用各向异性扩散方法对牙齿CBCT图像I₀进行预处理，得到预处理图像；

S3、采用自动阈值法从预处理图像中提取出牙齿区域I；

S4、根据牙齿区域I，采用三维可变尺度区域拟合分割模型得到牙齿三维分割结果。

2.根据权利要求1所述的牙齿CBCT图像三维分割方法，其特征在于，所述步骤S3包括以下分步骤：

S31、设置阈值t₁去除预处理图像中的背景区域，得到去背景图像；

S32、计算去背景图像的灰度直方图，通过阈值大律法，在灰度直方图的最大化类间方差和最小化类间方差之间计算得到一个最佳阈值t₂，通过最佳阈值t₂分割去背景图像中的牙齿区域I₁和软组织区域I₂；

S33、采用形态学方法中的闭操作对牙齿区域I₁进行补洞操作，得到牙齿区域I。

3.根据权利要求2所述的牙齿CBCT图像三维分割方法，其特征在于，所述步骤S31中阈值t₁的计算公式为：

t₁＝maxT×2％

其中maxT为预处理图像的最大灰度值。

4.根据权利要求1所述的牙齿CBCT图像三维分割方法，其特征在于，所述步骤S4包括以下分步骤：

S41、假设牙齿区域I表示为I:Ω→R^d，其中Ω表示图像域，d表示图像维数且d≥1；

S42、设曲线C为图像域Ω中的一条闭合曲线，将整个牙齿区域I分为内部区域Ω₁＝inside(C)和外部区域Ω₂＝outside(C)；

S43、定义三维空间中任意一点x∈Ω满足局部拟合能量函数E_x(C,f₁(x),f₂(x))，其中f₁(x)和f₂(x)分别表示在以x为圆心的小区域内Ω₁＝inside(C)和Ω₂＝outside(C)的各种局部区域内的近似拟合值；

S44、对局部拟合能量函数E_x(C,f₁(x),f₂(x))进行积分，得到能量函数E(C,f₁(x),f₂(x))；

S45、将能量函数E(C,f₁(x),f₂(x))转化为等价水平集函数

其中

表示水平集函数；

S46、在等价水平集函数

中引入距离正则项

和长度项

得到完整水平集形式的能量泛函

S47、将能量泛函

的极小化形式通过梯度下降流表示，得到水平集函数

的0-水平集曲线；

S48、不断迭代更新能量泛函

的极小化形式，直到水平集函数

的0-水平集曲线到达真实的牙齿边界处，得到牙齿三维分割结果。

5.根据权利要求4所述的牙齿CBCT图像三维分割方法，其特征在于，所述步骤S43中的局部拟合能量函数E_x(C,f₁(x),f₂(x))具体为：

其中λ₁,λ₂均为正常数，K(·)表示核函数，I(·)表示牙齿区域，y表示图像域Ω中的所有像素点。

6.根据权利要求4所述的牙齿CBCT图像三维分割方法，其特征在于，所述步骤S45中的等价水平集函数

具体为：

其中K_σ(·)表示高斯核函数，H(·)表示赫维赛德函数，y表示图像域Ω中的所有像素点。

7.根据权利要求4所述的牙齿CBCT图像三维分割方法，其特征在于，所述步骤S46中的距离正则项

和长度项

具体为：

其中δ(·)表示狄拉克函数，

表示水平集函数的梯度。

8.根据权利要求4所述的牙齿CBCT图像三维分割方法，其特征在于，所述步骤S46中的能量泛函

具体为：

其中μ,v均为非负常数。

9.根据权利要求4所述的牙齿CBCT图像三维分割方法，其特征在于，所述步骤S47中水平集函数

的0-水平集曲线表示为：

其中中间函数e₁(x)和e₂(x)分别表示为：

其中δ(·)表示狄拉克函数，

表示水平集函数的梯度，K_σ(·)表示高斯核函数，I(·)表示牙齿区域，y表示图像域Ω中的所有像素点，div(·)表示散度函数。

10.一种牙齿CBCT图像三维分割系统，其特征在于，所述牙齿CBCT图像三维分割系统用于配置执行如权利要求1-9任一所述的牙齿CBCT图像三维分割方法。

Images