CN107203998B - 一种对锥束ct图像进行牙列分割的方法 - Google Patents

一种对锥束ct图像进行牙列分割的方法 Download PDF

Info

Publication number
CN107203998B
CN107203998B CN201610157705.3A CN201610157705A CN107203998B CN 107203998 B CN107203998 B CN 107203998B CN 201610157705 A CN201610157705 A CN 201610157705A CN 107203998 B CN107203998 B CN 107203998B
Authority
CN
China
Prior art keywords
voxel
image
voxels
dentition
segmentation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201610157705.3A
Other languages
English (en)
Other versions
CN107203998A (zh
Inventor
裴玉茹
艾兴胜
查红彬
许天民
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peking University
Original Assignee
Peking University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peking University filed Critical Peking University
Priority to CN201610157705.3A priority Critical patent/CN107203998B/zh
Publication of CN107203998A publication Critical patent/CN107203998A/zh
Application granted granted Critical
Publication of CN107203998B publication Critical patent/CN107203998B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20092—Interactive image processing based on input by user
    • G06T2207/20104—Interactive definition of region of interest [ROI]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30036—Dental; Teeth

Abstract

本发明公布了一种对锥束CT(CBCT)图像进行牙列分割的方法，针对锥束CT图像中的感兴趣体图像区域定义图结构，基于半监督的随机游走算法和三维可变形模型配准定义的柔化约束，从CBCT图像中分割得到完整牙列。其中，利用三维可变形模型引入体图像的柔化约束，用于处理基于半监督标签扩散的体图像分割中的噪声；还采用迭代修正方法，对柔化约束下标记扩散以及三维模型对表面体素集的拟合问题，进行迭代求解，可有效地消除分割误差，改进单次标签扩散所获取的牙列分割，提高分割结果的精度，可满足口腔临床对于精度的要求。

Description

一种对锥束CT图像进行牙列分割的方法
技术领域
本发明涉及计算机视觉和图像处理技术领域，尤其涉及一种对锥束CT(CBCT)图像进行牙列分割的方法。
背景技术
锥束CT(CBCT)图像常用于颌面与口腔正畸外科辅助手术预测以及牙列对齐计划的制定。CBCT图像可以提供病人特定的解剖结构信息，其中包含完整的牙列信息。而传统的基于光学的方式，例如三维激光扫描以及立体视觉设备，仅仅能够获取牙冠的几何信息，其中由于牙根被埋藏在牙龈以及颌骨中无法获取其几何形态。在临床口腔外科中，CBCT图像由于其低放射剂量以及较低的采集成本，相对于传统的CT成像技术具有极大的优势。但是，较低的放射剂量以及信噪比一般会造成CBCT图像质量较差，常常出现图像退化现象，例如口腔治疗以及牙齿种植在牙颌结构中放置的金属物体造成的线束硬化(beam hardening)问题，有限视域的截断，以及由于牙齿与周围牙槽骨相似的灰度而造成的模糊的轮廓。同时在图像采集过程中，病人可能有微小的头部运动，也会造成牙齿轮廓的模糊。特别是为了诊断牙齿畸形，CBCT在采集过程中通常上下牙列处于咬合状态，这使得上下牙列在咬合部分的分割变得更加困难。
近年来有大量技术被应用于从CT以及CBCT图像中进行牙列分割，其中包括使用积分灰度投影、基于阈值和区域增长的方法、图割法、基于水平集活动轮廓的方法等。积分灰度投影方法仅仅能获取牙列的粗略分割，例如不同牙齿之间的分割平面。阈值与区域增长方法对于CBCT中可能出现的牙列与牙槽骨以及其它组织之间模糊的轮廓难以有效处理，其分割与真实的牙齿轮廓的一致性较差。基于交互的分割方法例如图割法非常依赖初始前景区域的定义，对于不同的初始前景定义所得到的分割常常存在差异，难以重复分割结果。为了获取可靠的分割，活动轮廓方法诸如水平集技术集成形状以及灰度先验、从切片图像分割获取的三维形状先验、以及隐参数表达模型，逐层对CBCT图像进行分割。在基于水平集的方法中，对于约束项需要进行精心设计以避免分割过程中可能出现的收缩(shrinkage)以及泄露(leakage)问题。此外，在基于活动轮廓方法进行逐层分割中还可能有误差累计的问题。此外，三维牙列模型被引入分割系统，包括三维形状图集以及三维形状统计模型。图集以及统计形状模型一般来自大量的CT图像，这势必会增加对数据的需求以及数据处理的负担。
发明内容
为了克服上述现有技术的不足，本发明提供一种对锥束CT(CBCT)图像进行牙列分割的方法，基于可变形模型的随机游走方法，从CBCT图像中分割得到完整牙列。
本发明的原理是：将体图像中的分割问题定义为图结构中的标签扩散问题，其中，预先对少量切片图像进行标注，并使用扩散机制对感兴趣体图像(volume of interest，VOI)中其它像素进行自动标注。首先，利用原始的随机游走方法对CBCT图像中牙列进行初始分割，基于原始的随机游走的初始分割通常会包含分割错误，例如由于CBCT图像退化出现的线束硬化、图像中模糊的牙列轮廓所造成的标记扩散误差。因此，引入三维可变形模型用于拟合由标记扩散得到的体图像分割，利用组合Logistic函数，并基于经验牙本质厚度定义柔化约束，将柔化约束引入随机游走算法中更新体素的类别标签的估计。此外，为了改善从单步柔化约束下的标记扩散所得到的图像分割精度，引入迭代修正算法，针对系统的两个方面，即柔化约束下标记扩散所对应的大规模稀疏线性系统和三维模型到图像分割所得到的牙列表面体素的拟合，进行迭代求解。在每步迭代中，三维可变形模型的配准可以看作是对体素标记的正则化，该过程可以有效地消除分割误差。
本发明提供的技术方案是：
一种对锥束CT图像进行牙列分割的方法，针对锥束CT图像中的感兴趣体图像进行牙列分割，得到牙列表面模型；包括如下步骤：
1)针对锥束CT图像中的感兴趣体图像区域定义图结构；所述图结构中，结点集合对应感兴趣体图像区域中所有的体素；边连接集合为感兴趣体图像区域中相邻体素之间的边连接；
2)根据牙列分布情况定义类别标签集合，针对步骤1)所述感兴趣体图像中所有的体素，采用随机游走方法进行牙列的初始分割，得到所述所有体素的类别，对所述所有体素进行类别标记；
3)根据步骤2)得到的所有体素的类别，利用其中属于牙列的表面体素定义牙列的表面形态，生成初始分割得到的表面体素集；进一步利用三维可变形模型，通过非刚性变形配准拟合所述牙列表面体素集合，得到三维可变形模型的非刚性变换参数，用于拟合所述初始分割得到的表面体素集；
4)基于三维可变形模型定义柔化约束，并基于柔化约束下的随机游走方法进行牙列分割，得到感兴趣体图像区域中所有体素的类别标签；
5)通过迭代修正过程对牙列分割进行迭代修正，结合步骤4)所述柔化约束下的类别标签扩散和步骤3)所述三维可变形模型的非刚性配准过程，改进单次柔化约束下随机游走方法所得到的体图像分割结果。
针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤1)所述图结构中，采用基于上下文的体素特征描述子，通过体素的表观计算相邻体素之间的相似度，用于表示体素的上下文特征差异；所述体素特征描述子向量中的元素对应当前体素的包围图像块与模式中的采样体素的包围图像块之间的像素灰度差异累计；图结构中边连接上的权重通过体素上下文特征描述差异和灰度差异的加权组合得到。
上述上下文的体素描述子在当前体素的一个半径为
Figure BDA0000944370090000031
的包围球体中通过随机采样生成一个模式P，所述模式P由采样得到的体素集合表示；所述体素特征描述子向量f(vi)中的元素对应当前体素的包围图像块与模式P中的采样体素的包围图像块之间的像素灰度差异累计，表示为式11：
f(vi)＝{dp(C(vi),C(vi+γ))|vi+γ∈P} (式11)
式11中，f(vi)表示当前体素vi的特征描述子向量；C(vi)表示体素vi的包围图像块，dp为图像块的灰度差异；γ为当前体素与模式P中采样体素之间的位移；C(vi+γ)表示模式中的采样体素的包围图像块；其中，图像块的灰度差异dp定义为式12：
Figure BDA0000944370090000032
式12中，dp为图像块的灰度差异；C(vi)与C(vk)分别表示体素vi与vk的包围图像块；δv为包围体素块中体素的位移；vi+δv属于图像块C(vi)；I(vi+δv)为体素vi+δv的灰度；I(vk+δv)为体素vk+δv的灰度；
所述图像块的灰度差异可通过图像的线性卷积算子计算得到：对差异图像V-V(γ)进行卷积，通过式13得到：
dp＝(V-V(γ))2*G (式13)
式13中，dp为图像块的灰度差异；γ为体素vi与模式中一个采样体素之间的位移向量；对体图像进行平移操作得到V(γ)；V-V(γ)为差异图像；G为卷积核；
图结构中边连接上的权重通过上下文特征描述差异与灰度差异的加权组合得到，如式14：
Figure BDA0000944370090000041
式14中，Wij为体素i与相邻体素j之间边连接的权重；α为常数系数用于调节上下文特征以及体素灰度对于权重计算的影响；Ii与Ij对应第i个与第j个体素的灰度，fi与fj对应第i个与第j个体素的上下文描述；Ψ为图结构的边连接集合。
针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤2)采用随机游走方法对感兴趣体图像区域中所有的体素的类别进行标记，得到牙列的初始分割结果；具体通过最小化式21的能量函数得到感兴趣体图像中体素类别标记：
Figure BDA0000944370090000042
式21中，Erw(X)表示基于随机游走算法的图像分割能量；|S|为系统类别标签的个数；
Figure BDA0000944370090000043
表示第i个结点被分为第s类的概率；系统中每个结点对均对应一个|S|维的概率向量，其中第i个结点所属的类别定义为
Figure BDA00009443700900000410
表示预先由用户交互定义的结点概率，当第i个结点被用户交互定义为第s类，则
Figure BDA0000944370090000046
取值为1，否则为0；nL是用户预先交互标注的体素个数；nV为感兴趣体图像中所有体素的个数；μ1与μ2为常数；L为拉普拉斯矩阵；H为对角线指示矩阵，当感兴趣体图像中的第i个结点被预先标注时，H的元素Hii＝1；X与
Figure BDA0000944370090000047
表示感兴趣体图像中所有体素对应的类别概率
Figure BDA0000944370090000048
以及
Figure BDA0000944370090000049
所组成的矩阵；。
其中，将所述式21的目标函数转化为矩阵形式，通过将所述式21的一阶导数设为0，将目标函数求解转化为求解式22所示的大规模稀疏的线性系统：
Figure BDA0000944370090000051
式22中，
Figure BDA0000944370090000052
表示单位矩阵；L为拉普拉斯矩阵；H为对角线指示矩阵，当感兴趣体图像中的第i个结点被预先标注时，H的元素Hii＝1；X与
Figure BDA0000944370090000053
表示感兴趣体图像中所有体素对应的类别概率
Figure BDA0000944370090000054
以及
Figure BDA0000944370090000055
所组成的矩阵；
通过求解式22的线性系统，得到感兴趣体图像中所有体素的类别，作为牙列的初始分割结果。
针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤3)所述三维可变形模型用于处理锥束CT图像自身的局部退化造成的边界混淆问题；所述通过非刚性变形配准拟合所述牙列表面体素集合，具体包括如下步骤：
31)利用Procrustes分析计算三维模型的全局变换使其最大程度地与表面体素集Yc一致；
32)随后利用非刚性配准获取特定的三维模型，使得所述特定的三维模型与初始分割的表面体素集一致；
33)为降低非刚性变换的参数空间，利用嵌入变形方法，将所述特定的三维模型表面的非刚性变形定义为一个稀疏基网格顶点变换的线性组合，三维模型非刚性配准的能量函数定义如式31：
Figure BDA0000944370090000056
式31中，
Figure BDA0000944370090000057
为变形前后表面顶点的位移；变形后的三维模型Ye的顶点为
Figure BDA0000944370090000058
其中，权重βij对应顶点
Figure BDA0000944370090000059
与其在基网格B上的近邻顶点ηj，并基于两点之间的欧式距离定义为
Figure BDA00009443700900000510
其中参数
Figure BDA00009443700900000511
Tj为第j个基网格顶点上的三维变换矩阵；能量函数Ere的第一项用于最小化变形后的三维模型表面与体图像分割获取的表面体素集之间的距离；dc表示变换后的顶点
Figure BDA00009443700900000512
与标签扩散得到的表面体素集Yc之间的距离,
Figure BDA00009443700900000513
其中
Figure BDA00009443700900000514
为表面体素；第二项用于保持三维模型表面平滑，该项最小化顶点
Figure BDA0000944370090000061
与其近邻
Figure BDA0000944370090000062
的空间变换之间的差异；ω为常数系数，ne为三维模型表面顶点个数。
34)通过最小化能量函数Ere，得到三维模型的非刚性变换参数，用于配准拟合所述牙列表面体素集合。
针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤4)所述基于柔化约束下的随机游走方法进行牙列分割，包括如下步骤：
41)通过设定与体图像分割结果拟合后的三维模型表面内部体素更有可能属于牙列，根据体素结点概率定义柔化约束；
42)设定柔化约束下的感兴趣体图像的标签扩散能量为式43：
Erwsc(X)＝λ1Erw(X)+λ2Esc(X) (式43)
其中Erw与Esc式21式42所定义的随机游走以及柔化约束的能量。λ1与λ2为常数系数。
将式43中的能量函数的一阶导数设为0，得到关于体素结点概率向量的线性系统如式44，
Figure BDA0000944370090000063
其中L为拉普拉斯矩阵；H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；X与
Figure BDA0000944370090000064
表示VOI中所有体素对应的类别概率
Figure BDA0000944370090000065
以及用户定义的先验概率
Figure BDA0000944370090000066
所组成的矩阵。
Figure BDA0000944370090000067
对应感兴趣区域中所有体素的柔化约束所构成的矩阵；λ1与λ2为常数系数,实验中分别设置为1、0.1；μ1与μ2为常数，本实验中分别设置为1与0.001；
Figure BDA0000944370090000068
表示单位矩阵。
43)利用共轭梯度方法求解线性系统，该线性系统需要求解|S|-1次，其中|S|为系统中类别标签的个数；通过求解线性系统得到感兴趣体图像中所有体素的类别标签。
进一步地，步骤41)采用以下步骤定义所述柔化约束：
410)根据牙齿表面内部的灰度分布，对牙本质赋以更大的属于牙齿的概率，定义基于三维可变形模型属于前景牙齿的概率为式41：
Figure BDA0000944370090000069
式41中，pd(r,θ)为基于三维可变形模型属于前景牙齿的概率；r1(θ)与r2(θ)表示牙齿的外轮廓以及牙髓腔轮廓，其中r2(θ)＝r1(θ)-a，a为经验的牙本质厚度；
Figure BDA0000944370090000071
参数η用于控制如上概率函数经过牙齿轮廓时函数的形状；
411)针对没有牙髓腔的分层图像，基于三维可变形模型属于前景牙齿的概率定义为式411：
Figure BDA0000944370090000072
式411中，
Figure BDA0000944370090000073
412)针对牙齿可能出现的多牙根以及多牙尖结构情况，牙齿在切片上呈多轮廓状态时基于三维可变形模型属于前景牙齿的概率表示为式412：
Figure BDA0000944370090000074
式412中，nc为当前牙齿在切片图像中轮廓的个数；pd为单外轮廓的基于三维可变形模型属于前景牙齿的概率(式41与式411)；r,θ为体素在极坐标系中的坐标；
413)利用卡方核的一对多支持向量机生成分类器，所述分类器的输出被归一化到0-1之间，表示基于体素表观属于前景(牙列)的概率pc，柔化约束定义为式413：
Figure BDA0000944370090000075
式413中，pc表示由体素vi表观估计的体素属于前景(牙列)的概率；
Figure BDA0000944370090000076
表示在体素vi基于三维可变形模型属于前景牙齿的概率；r,θ为体素在极坐标系中的坐标；
Figure BDA0000944370090000077
是对体素vi相对于类别s所定义的柔化约束。
针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤5)所述通过迭代修正过程对牙列分割进行迭代修正，具体包括如下步骤：
51)将柔化约束下的标签扩散和三维模型的非刚性配准能量结合，定义能量如式51：
E(X,T)＝λ1Erw+λ2Esc+λ3Ere (式51)
式51中，E(X,T)是基于三维可变形模型随机游走算法进行牙列分割的总能量函数；Erw与Esc为式21与式42所定义的随机游走以及柔化约束的能量函数；Ere为式31所定义三维模型非刚性配准的能量；λ1，λ2与λ3为常数，实验中分别设置为1、0.1、1；；
52)在第一个阶段，设定基于标签扩散得到的体图像分割获取牙列表面体素集，由于表面体素集由第i-1次迭代中的CBCT图像分割得到，所以将与柔化约束下的随机游走相关的项Erw和Esc从能量函数式51中去掉，通过最小化能量E(X(i-1),T)对三维模型进行变形，用于拟合牙列的表面体素集；此阶段对于三维模型变形参数的求解和步骤3)相同；
53)在第二个阶段，去掉能量函数中与三维模型配准相关的项Ere，最小化能量E(X,T(i-1))，得到感兴趣图像中所有体素的标签；此阶段通过求解一个大规模稀疏的线性系统求解所有体素的标签，和步骤4)所述柔化约束下的随机游走方法相同；
54)迭代优化进行步骤52)和步骤53)，当相邻两步迭代中三维模型更新||Ye(i)-Ye (i-1)||2小于预先给定的阈值时终止迭代，得到感兴趣图像中所有体素的标签，作为迭代修正后的牙列分割结果。其中Ye(i)与Ye(i-1)分别为第i步与第i-1步三维表面模型。
针对上述对锥束CT图像进行牙列分割的方法，进一步地，所述得到的牙列表面模型与手工分割得到的表面模型的平均距离误差在0.3mm以下。
与现有技术相比，本发明的有益效果是：
本发明提供一种对锥束CT(CBCT)图像进行牙列分割的方法，基于可变形模型的随机游走方法，从CBCT图像中分割得到完整牙列。本发明方法结合半监督的随机游走算法和三维可变形模型配准定义的柔化约束，获取牙列的可靠分割。其中，利用三维可变形模型引入体图像的柔化约束，可以看作为正则化项用于处理基于半监督标签扩散的体图像分割中的噪声。本发明还采用迭代修正方法，对柔化约束下标记扩散以及三维模型对表面体素集的拟合问题，进行迭代求解，可有效地消除分割误差，改进单次标签扩散所获取的牙列分割，提高分割结果的精度。
利用本发明提供的方法，可以对医学临床CBCT图像进行自动的牙列分割，分割得到的三维牙列模型可用于颌面与口腔正畸外科辅助手术预测以及牙列对齐计划的制定，该牙列分割可满足口腔临床对于精度的要求。度量自动分割所得到的牙列表面模型与手工分割对应的表面模型的平均距离，该距离误差在0.3mm以下。
附图说明
图1是本发明提供方法的流程框图。
图2是本发明中采用的基于上下文的体素描述子的结构示意图；
其中，vi为当前体素；
Figure BDA0000944370090000091
为当前体素的一个包围球体的半径；γ为当前体素与模式P中采样体素之间的位移；C(vi)为当前体素vi的包围图像块；C(vi+γ)为模式中的采样体素的包围图像块。
具体实施方式
下面结合附图，通过实施例进一步描述本发明，但不以任何方式限制本发明的范围。
CBCT图像是三维体图像，由一系列的二维图像组成，本说明将其中二维图像称作切片图像(分层图像)。本发明实施例针对医学临床CBCT图像，基于CBCT图像中感兴趣体图像区域中定义的图结构和少量用户交互标注的分层图像(切片图像)，对CBCT图像中的牙列进行分割。其中，采用三维可变形模型定义柔化约束，通过基于柔化约束下的随机游走算法更新体图像中的牙列分割，再以迭代修正方法获取可靠的分割结果。
图1是本发明提供方法的流程框图。具体地，本发明提供方法包括如下步骤：
步骤1：设定CBCT图像中的感兴趣体图像(VOI)，定义VOI中的图结构；
本方法以一种半监督的方式对CBCT图像进行牙列分割，其中，将用户所指定的少量切片图像中的标记扩散到感兴趣体图像的整个区域。该扩散过程在VOI中定义的图结构中进行。针对VOI定义的图结构中，结点集合对应VOI中所有的体素，边连接集合为相邻体素之间的边连接。
一般地，图结构的边连接集合包含两类，一类对应切片图像内部的边连接，该连接对应切片图像中规则的像素网格；另一类对应切片图像之间的边连接，依据由稠密光流算法获取的切片图像之间的对应体素建立边连接。值得注意的是，由稠密光流算法建立的对应并非一一对应，即第i层切片图像中的一个结点可能与第(i+1)层切片图像中不止一个结点对应。该情况常常发生在牙齿区域分叉的切片图像，例如磨牙的多牙根区域，或者牙冠中的多牙尖区域。
系统基于体素的表观计算相邻体素之间的相似度(为上下文特征描述差异与灰度差异的加权组合)。使用体素的灰度I计算体素的相似度是一种较直观的方式，但是仅仅利用体素灰度定义体素相似度对于低辐射剂量采集的CBCT图像并不适合。这是由于CBCT图像中可能存在退化区域，其中仅由灰度相似并不能说明对应体素相似。因此系统引入基于上下文的体素描述子。
图2是基于上下文的体素描述子；其中，当前体素的一个半径为
Figure BDA0000944370090000102
的包围球体中，γ为当前体素与模式P中采样体素之间的位移；通过随机采样生成模式P；特征描述子向量f(vi)中的元素对应当前体素vi的包围图像块C(vi)与模式中的采样体素的包围图像块C(vi+γ)之间的像素灰度差异累计。
为了降低计算的复杂度，该上下文的体素描述子并未对当前体素的所有近邻体素进行计算，而是在当前体素的一个半径为
Figure BDA0000944370090000103
的包围球体中通过随机采样生成一个模式P。该模式由采样得到的体素集合表示。特征描述子向量f(vi)中的元素对应当前体素的包围图像块与模式中的采样体素的包围图像块之间的像素灰度差异累计，表示为式11：
f(vi)＝{dp(C(vi),C(vi+γ))|vi+γ∈P} (式11)
式11中，f(vi)表示当前体素vi的特征描述子向量；C(vi)表示体素vi的包围图像块，dp为图像块的灰度差异；γ为当前体素与模式P中采样体素之间的位移；C(vi+γ)表示模式中的采样体素的包围图像块。其中，图像块的灰度差异dp定义为式12：
Figure BDA0000944370090000101
式12中，dp为图像块的灰度差异；C(vi)与C(vk)分别表示体素vi与vk的包围图像块；δv为包围体素块中体素的位移；vi+δv属于图像块C(vi)；I(vi+δv)为体素vi+δv的灰度；I(vk+δv)为体素vk+δv的灰度；
图像块的灰度差异可以看作图像块之间对应体素灰度差异的求和，该计算可以通过图像的线性卷积算子实现。如果体素vi与模式中一个采样体素之间的位移向量为γ，则可预先对体图像进行平移操作得到V(γ)。图像块的灰度差异dp则可以对差异图像V-V(γ)进行卷积得到，即通过式13得到：
dp＝(V-V(γ))2*G (式13)
式13中，dp为图像块的灰度差异；G为卷积核,V(γ)为对原始体图像V进行平移γ后得到的体图像。
结合体素的上下文特征描述，图结构中边连接上的权重定义为上下文特征描述差异与灰度差异的加权组合，如式14：
Figure BDA0000944370090000111
式14中，Wij为体素i与相邻体素j之间边连接的权重；α为常数系数用于调节上下文特征以及体素灰度对于权重计算的影响；Ii与Ij对应第i个与第j个体素的灰度，fi与fj对应第i个与第j个体素的上下文描述。Ψ为图结构的边连接集合；本实施例中，系统预先将CBCT图像的灰度范围变换到[0,255]区间，其中α被设置为0.1。通过式14计算得到体素的相似度矩阵。
步骤2：采用随机游走算法进行牙列的初始分割；
系统根据牙列分布情况定义标签集合S，其中对不同的牙齿分配不同的标签，正常的牙列包含32颗牙齿，其中第三磨牙，即智齿通常未发育完全，系统将其排除在外，还余下28颗牙齿，因而加上背景类别，系统的标签集合中将定义29类标签。但在实际CBCT数据中，由于牙齿缺失等问题，实际CBCT图像中包含的类别可能少于该类别个数。在初始牙列分割中，系统采用原始的随机游走算法对感兴趣体图像区域进行标记。VOI中体素标记可以通过最小化如下的能量函数(式21)得到。
Figure BDA0000944370090000112
式21中，Erw(X)表示基于随机游走算法的图像分割能量；|S|为系统类别标签的个数；
Figure BDA0000944370090000113
表示第i个结点被分为第s类的概率；系统中每个结点均对应一个|S|维的概率向量，其中第i个结点所属的类别定义为
Figure BDA0000944370090000114
表示预先由用户交互定义的结点概率，当第i个结点被用户交互定义为第s类，则
Figure BDA0000944370090000121
取值为1，否则为0；nL是用户预先交互标注的体素个数；nV为VOI中所有体素的个数；μ1与μ2为常数，本实验中分别设置为1与0.001。上述目标函数(式21)可以转化为矩阵形式，其中L为拉普拉斯矩阵，在给定步骤1)中通过公式14得到的相似度矩阵W后，拉普拉斯矩阵L＝D-W；其中，对角线矩阵D中元素定义为
Figure BDA0000944370090000122
H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；X与
Figure BDA0000944370090000123
表示VOI中所有体素对应的类别概率
Figure BDA0000944370090000124
以及
Figure BDA0000944370090000125
所组成的矩阵。通过将上面的能量函数的一阶导数设为0，该能量函数可以转化为求解如下大规模稀疏的线性系统，表示为式22：
Figure BDA0000944370090000126
式22中，
Figure BDA0000944370090000127
表示单位矩阵；L为拉普拉斯矩阵；H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；μ1与μ2为常数，本实验中分别设置为1与0.001。给定用户定义的形状先验
Figure BDA0000944370090000128
可通过求解上面的线性系统(式22)得到VOI中所有体素的类别，即得到牙列的初始分割结果。
步骤3：利用三维可变形模型，通过非刚性变形配准拟合步骤2进行图像分割所得到的牙列表面体素集合，得到三维可变形模型的非刚性变换参数；
步骤2得到VOI中所有体素的类别，利用其中属于牙列表面的体素，可以定义牙列的表面形态；考虑到CBCT图像中可能存在的退化现象，以及处于咬合状态上下颌接触的牙齿，基于原始的随机游走技术获取的初始分割与真实牙齿轮廓之间会存在较大的差异。系统引入三维可变形模型处理CBCT图像自身的局部退化问题所造成的边界混淆问题。三维可变形模型来自对实体牙齿模型的三维光学扫描，其具有良好的拓扑结构定义。从CBCT图像分割可以定义牙列表面的体素集合，即属于牙列的体素如果其近邻体素属于背景，则该体素被放入表面体素集Yc。对三维可变形模型施加非刚性变换，以拟合初始分割得到的表面体素集。与体图像分割结果拟合后的三维模型表面内部体素更有可能属于牙列，并基于该假设定义柔化约束。
三维模型非刚性拟合分为两个阶段，首先利用Procrustes分析计算三维模型的全局变换使其最大程度地与表面体素集Yc一致；随后利用非刚性配准获取特定病人的三维模型，使其与初始分割的表面体素集一致。为了降低非刚性变换的参数空间，系统利用嵌入变形技术，将三维模型表面的非刚性变形定义为一个稀疏基网格顶点变换的线性组合，其中变形后的三维模型顶点为
Figure BDA0000944370090000131
其中，权重βij对应顶点
Figure BDA0000944370090000132
与其在基网格B上的近邻顶点ηj，并基于两点之间的欧式距离定义为
Figure BDA0000944370090000133
其中参数
Figure BDA0000944370090000134
Tj为第j个基网格顶点上的三维变换矩阵，三维模型非刚性配准的能量函数定义如式31：
Figure BDA0000944370090000135
式31中，
Figure BDA0000944370090000136
即变形前后表面顶点的位移；能量函数Ere的第一项用于最小化变形后的三维模型表面与体图像分割获取的表面体素集之间的距离；dc表示变换后的顶点
Figure BDA0000944370090000137
与标签扩散得到的表面体素集Yc之间的距离,
Figure BDA0000944370090000138
其中
Figure BDA0000944370090000139
为表面体素；第二项用于保持三维模型表面平滑，该项最小化顶点
Figure BDA00009443700900001310
与其近邻
Figure BDA00009443700900001311
的空间变换之间的差异。ω为常数系数，ne为三维模型表面顶点个数。通过最小化能量函数Ere可以得到三维模型的非刚性变换参数。
步骤4：基于柔化约束下的随机游走算法进行牙列分割，得到VOI中所有体素的类别标签，即哪些体素属于牙列(前景)，哪些体素属于背景；
可以直观看到，与体图像分割得到的表面体素集配准后的三维模型表面内部的体素更有可能属于牙列。但是，在体图像中牙齿表面内部体素并不具有均匀的灰度，其中牙本质具有较亮的灰度，而牙髓腔灰度较低，考虑到牙齿表面内部的灰度分布，系统对牙本质赋以更大的属于牙齿的概率。由于三维模型来自对牙齿模型的光学扫描，该模型不包括内部的牙髓腔的表面，利用经验的牙本质厚度定义a，以及组合Logistic函数依据步骤3得到的变形后的三维模型表面定义柔化约束。极坐标系的中心定义为变形后三维模型与切片图像相交轮廓的中心，令r1(θ)与r2(θ)表示牙齿的外轮廓以及牙髓腔轮廓，其中r2(θ)＝r1(θ)-a。基于三维可变形模型属于前景牙齿的概率pd定义为式41：
Figure BDA00009443700900001312
式41中，pd(r,θ)为基于三维可变形模型属于前景牙齿的概率；r,θ为体素在极坐标系中的坐标；r1(θ)与r2(θ)表示牙齿的外轮廓以及牙髓腔轮廓，其中r2(θ)＝r1(θ)-a，a为经验的牙本质厚度；
Figure BDA0000944370090000141
参数η用于控制如上概率函数经过牙齿轮廓时函数的形状。在实验中，参数η的值设为0.99。
如果考虑没有牙髓腔的分层图像，基于三维可变形模型属于前景牙齿的概率定义为式411：
Figure BDA0000944370090000142
其中，pd(r,θ)为基于三维可变形模型属于前景牙齿的概率；r,θ为体素在极坐标系中的坐标；r1(θ)表示牙齿的外轮廓；
Figure BDA0000944370090000143
其中参数η用于控制如上概率函数经过牙齿轮廓时函数的形状。在实验中，参数η的值设为0.99。
上述的基于三维可变形模型属于前景牙齿的概率针对切片图像中的单轮廓定义，考虑到牙齿可能出现的多牙根以及多牙尖结构，该定义可以直接扩展到牙齿在切片上的多轮廓状态，表示为式412：
Figure BDA0000944370090000144
其中，nc为当前牙齿在切片图像中轮廓的个数；pd(r,θ)为式41与式411定义的在坐标(r,θ)处基于三维可变形模型属于前景牙齿的概率；
Figure BDA0000944370090000145
为扩展到牙齿在切片上的多轮廓状态时基于三维可变形模型属于前景牙齿的概率。
由于三维模型变形后与牙齿真实表面之间可能存在空隙，系统还利用体素的表观定义体素属于牙列的分类器。从交互标注的切片图像中学习该分类器，利用卡方(chi-square)核的一对多(one-vs-rest)支持向量机(SVM)生成分类器，该分类器的输出被归一化到0-1之间，表示体素属于前景(牙列)的概率pc。柔化约束定义为式413：
Figure BDA0000944370090000146
式413中，pc表示由体素表观估计的体素vi属于前景(牙列)的概率；
Figure BDA0000944370090000147
表示基于三维可变形模型vi属于前景牙齿的概率；r,θ为体素vi在极坐标系中的坐标。
Figure BDA00009443700900001513
是对体素vi相对于类别s所定义的柔化约束。
柔化约束的能量项定义为式42：
Figure BDA0000944370090000151
其中，
Figure BDA0000944370090000152
表示第i个结点被分为第s类的概率；X表示VOI中所有体素对应的类别概率
Figure BDA0000944370090000153
所组成的矩阵；|S|表示系统中类别个数；nV为感兴趣区域中体素的个数；
Figure BDA0000944370090000154
为式413所定义的柔化约束；
Figure BDA0000944370090000155
对应感兴趣区域中所有体素的柔化约束所构成的矩阵；给定柔化约束，组合式21与式42，将感兴趣体图像的标签扩散能量定义为式43：
Erwsc(X)＝λ1Erw+λ2Esc (式43)
其中Erw与Esc式21式42所定义的随机游走以及柔化约束的能量。λ1与λ2为常数系数。将该能量函数(式43)的一阶导数设为0，可得到关于体素结点概率向量的线性系统，利用共轭梯度方法求解如下线性系统(式44)。
Figure BDA0000944370090000156
其中L为拉普拉斯矩阵；H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；X与
Figure BDA0000944370090000157
表示VOI中所有体素对应的类别概率
Figure BDA0000944370090000158
以及用户定义的先验概率
Figure BDA0000944370090000159
所组成的矩阵。
Figure BDA00009443700900001510
对应感兴趣区域中所有体素的柔化约束所构成的矩阵；λ1与λ2为常数系数,实验中分别设置为1、0.1；μ1与μ2为常数，本实验中分别设置为1与0.001；
Figure BDA00009443700900001511
表示单位矩阵。由于xi是|S|维的向量，其中|S|为系统中类别标签的个数，该线性系统需要求解|S|-1次，其中
Figure BDA00009443700900001512
求解线性系统可以得到VOI中所有体素的类别标签。
步骤5：对牙列分割迭代修正；
为了改进单次柔化约束下随机游走算法所得到的体图像分割，系统引入迭代修正过程，对柔化约束下标记扩散，以及三维模型对表面体素集的拟合问题，进行迭代求解，以改进单次标签扩散所获取的牙列分割。
将柔化约束下的标签扩散以及三维模型的非刚性配准能量结合，定义能量如下：
E(X,T)＝λ1Erw+λ2Esc+λ3Ere (式51)
式51中，E(X,T)是基于三维可变形模型随机游走算法进行牙列分割的总能量函数；Erw与Esc为式21式42所定义的随机游走以及柔化约束的能量函数；Ere为式31所定义三维模型非刚性配准的能量；λ1，λ2与λ3为常数，实验中分别设置为1、0.1、1。
该能量函数的优化可以分解为两个子问题。在第一个阶段，给定基于标签扩散得到的体图像分割定义牙列表面体素集，对三维模型进行变形，用于拟合更新后的表面体素集。该变形通过最小化能量E(X(i-1),T)得到，其中表面体素集由第i-1次迭代中的CBCT图像分割得到，因而与柔化约束下的随机游走相关的项，即Erw与Esc可以从能量函数中去掉。对于三维模型变形参数的求解与步骤3)的描述一致。
在第二个阶段，最小化能量E(X,T(i-1))以得到VOI中所有体素的标签。由于在标签扩散过程中使用的柔化约束来自第i-1步迭代中的三维模型配准，因而能量函数中与三维模型配准相关的项，即Ere可被去掉，该阶段对应步骤4)中柔化约束下的随机游走算法，通过求解一个大规模稀疏的线性系统求解所有体素的标签。
随着两个子问题的优化迭代进行，变形后的三维物体表面与基于标签扩散所得到的表面体素集会越来越接近。当相邻两步迭代中三维模型更新||Ye(i)-Ye(i-1)||2(其中Ye(i)与Ye(i-1)分别为第i步与第i-1步三维表面模型)小于预先给定的阈值时，终止迭代，得到VOI中所有体素的标签，作为迭代修正后的牙列分割结果。
为了验证基于三维可变形模型随机游走算法对牙列进行分割的精度，利用Dice测度计算自动分割与手工分割的相似度，该相似度大于0.95。同时，度量自动分割所得到的牙列表面模型与手工分割对应的表面模型的平均距离，该距离误差在0.3mm以下，可满足口腔临床对于精度的要求。
利用本发明的方法，可以对临床采集CBCT图像进行牙列分割。该方法结合半监督的随机游走算法，以及由三维可变形模型配准定义的柔化约束获取牙列的可靠分割。其中利用三维可变形模型引入柔化约束，也可以看作为正则化项用于处理基于半监督标签扩散的体图像分割中的噪声。该发明引入迭代修正方法，其中对系统的两个子问题，即柔化约束下标记扩散，以及三维模型对表面体素集的拟合问题，进行迭代求解以改进单次标签扩散所获取的牙列分割。
需要注意的是，公布实施例的目的在于帮助进一步理解本发明，但是本领域的技术人员可以理解：在不脱离本发明及所附权利要求的精神和范围内，各种替换和修改都是可能的。因此，本发明不应局限于实施例所公开的内容，本发明要求保护的范围以权利要求书界定的范围为准。

Claims (6)

1.一种对锥束CT图像进行牙列分割的方法，针对锥束CT图像中的感兴趣体图像进行牙列分割，得到牙列表面模型；包括如下步骤：
1)针对锥束CT图像中的感兴趣体图像区域定义图结构；所述图结构中，结点集合对应感兴趣体图像区域中所有的体素；边连接集合为感兴趣体图像区域中相邻体素之间的边连接；所述图结构中，采用基于上下文的体素特征描述子，通过体素的表观计算相邻体素之间的相似度，用于表示体素的上下文特征差异；所述体素特征描述子向量中的元素对应当前体素的包围图像块与模式中的采样体素的包围图像块之间的像素灰度差异累计；图结构中边连接上的权重通过体素上下文特征描述差异和灰度差异的加权组合得到；
所述上下文的体素描述子在当前体素的一个半径为
Figure FDA0002152374200000011
的包围球体中通过随机采样生成一个模式P，所述模式P由采样得到的体素集合表示；所述体素特征描述子向量f(vi)中的元素对应当前体素的包围图像块与模式P中的采样体素的包围图像块之间的像素灰度差异累计，表示为式11：
f(vi)＝{dp(C(vi),C(vi+γ))|vi+γ∈P} (式11)
式11中，f(vi)表示当前体素vi的特征描述子向量；C(vi)表示体素vi的包围图像块，dp为图像块的灰度差异；γ为当前体素与模式P中采样体素之间的位移；C(vi+γ)表示模式中的采样体素的包围图像块；其中，图像块的灰度差异dp定义为式12：
Figure FDA0002152374200000012
式12中，dp为图像块的灰度差异；C(vi)与C(vk)分别表示体素vi与vk的包围图像块；δv为包围体素块中体素的位移；vi+δv属于图像块C(vi)；I(vi+δv)为体素vi+δv的灰度；I(vk+δv)为体素vk+δv的灰度；
所述图像块的灰度差异再通过图像的线性卷积算子计算得到：对差异图像V-V(γ)进行卷积，通过式13得到：
dp＝(V-V(γ))2*G (式13)
式13中，dp为图像块的灰度差异；γ为体素vi与模式中一个采样体素之间的位移向量；对体图像进行平移操作得到V(γ)；V-V(γ)为差异图像；G为卷积核；
图结构中边连接上的权重通过上下文特征描述差异与灰度差异的加权组合得到，如式14：
Figure FDA0002152374200000021
式14中，Wij为体素i与相邻体素j之间边连接的权重；α为常数系数用于调节上下文特征以及体素灰度对于权重计算的影响；Ii与Ij对应第i个与第j个体素的灰度，fi与fj对应第i个与第j个体素的上下文描述；Ψ为图结构的边连接集合；
2)根据牙列分布情况定义类别标签集合，针对步骤1)所述感兴趣体图像中所有的体素，采用随机游走方法进行牙列的初始分割，得到所述所有体素的类别，对所述所有体素进行类别标记；
步骤2)采用随机游走方法对感兴趣体图像区域中所有的体素的类别进行标记，得到牙列的初始分割结果；具体通过最小化式21的能量函数得到感兴趣体图像中体素类别标记：
Figure FDA0002152374200000022
式21中，Erw(X)表示基于随机游走算法的图像分割能量；|S|为系统类别标签的个数；
Figure FDA0002152374200000023
表示第i个结点被分为第s类的概率；系统中每个结点对均对应一个|S|维的概率向量，其中第i个结点所属的类别定义为
Figure FDA0002152374200000024
Figure FDA0002152374200000025
表示预先由用户交互定义的结点概率，当第i个结点被用户交互定义为第s类，则
Figure FDA0002152374200000026
取值为1，否则为0；nL是用户预先交互标注的体素个数；nV为感兴趣体图像中所有体素的个数；μ1与μ2为常数；上述目标函数式21可转化为矩阵形式，其中L为拉普拉斯矩阵，在给定步骤1)中通过式14得到的相似度矩阵W后，拉普拉斯矩阵L＝D-W；其中，对角线矩阵D中元素定义为Dii＝∑jWij；H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；X与
Figure FDA0002152374200000027
表示VOI中所有体素对应的类别概率
Figure FDA0002152374200000031
以及
Figure FDA0002152374200000032
所组成的矩阵；
3)根据步骤2)得到的所有体素的类别，利用其中属于牙列的表面体素定义牙列的表面形态，生成初始分割得到的表面体素集；进一步利用三维可变形模型，通过非刚性变形配准拟合所述牙列表面体素集合，得到三维可变形模型的非刚性变换参数，用于拟合所述初始分割得到的表面体素集；
4)基于三维可变形模型定义柔化约束，再基于柔化约束下的随机游走方法进行牙列分割，得到感兴趣体图像区域中所有体素的类别标签；
所述基于柔化约束下的随机游走方法进行牙列分割，包括如下步骤：
41)通过设定与体图像分割结果拟合后的三维模型表面内部体素更有可能属于牙列，根据体素结点概率定义柔化约束；柔化约束的能量项定义为式42：
Figure FDA0002152374200000033
其中，
Figure FDA0002152374200000034
表示第i个结点被分为第s类的概率；X表示VOI中所有体素对应的类别概率
Figure FDA0002152374200000035
所组成的矩阵；|S|表示系统中类别个数；nV为感兴趣区域中体素的个数；
Figure FDA0002152374200000036
为定义的柔化约束；
Figure FDA0002152374200000037
对应感兴趣区域中所有体素的柔化约束所构成的矩阵；
42)设定柔化约束下的感兴趣体图像的标签扩散能量为式43：
Erwsc(X)＝λ1Erw+λ2Esc (式43)
其中，Erw与Esc分别为随机游走的能量项和柔化约束的能量项；
将式43中的能量函数的一阶导数设为0，得到关于体素结点概率向量的线性系统如式44，
Figure FDA0002152374200000038
其中L为拉普拉斯矩阵；H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；X与
Figure FDA0002152374200000039
表示VOI中所有体素对应的类别概率
Figure FDA00021523742000000310
以及用户定义的先验概率
Figure FDA00021523742000000311
所组成的矩阵；
Figure FDA00021523742000000312
对应感兴趣区域中所有体素的柔化约束所构成的矩阵；λ1与λ2为常数系数；μ1与μ2为常数；
Figure FDA00021523742000000313
表示单位矩阵；
43)利用共轭梯度方法求解线性系统，该线性系统需要求解|S|-1次，其中
Figure FDA00021523742000000314
|S|为系统类别标签的个数；通过求解线性系统得到感兴趣体图像中所有体素的类别标签；
5)通过迭代修正过程对牙列分割进行迭代修正，结合步骤4)所述柔化约束下的类别标签扩散和步骤3)所述三维可变形模型的非刚性配准过程，改进单次柔化约束下随机游走方法所得到的体图像分割结果。
2.如权利要求1所述对锥束CT图像进行牙列分割的方法，其特征是，将所述式21的目标函数转化为矩阵形式，通过将所述式21的一阶导数设为0，将目标函数求解转化为求解式22所示的大规模稀疏的线性系统：
Figure FDA0002152374200000041
式22中，
Figure FDA0002152374200000042
表示单位矩阵；L为拉普拉斯矩阵；H为对角线指示矩阵，当感兴趣体图像中的第i个结点被预先标注时，H的元素Hii＝1；X与
Figure FDA0002152374200000043
表示感兴趣体图像中所有体素对应的类别概率
Figure FDA0002152374200000044
以及
Figure FDA0002152374200000045
所组成的矩阵；μ1与μ2为常数，可分别设置为1与0.001；
通过求解式22的线性系统，得到感兴趣体图像中所有体素的类别，作为牙列的初始分割结果。
3.如权利要求1所述对锥束CT图像进行牙列分割的方法，其特征是，步骤3)所述三维可变形模型用于处理锥束CT图像自身的局部退化造成的边界混淆问题；所述通过非刚性变形配准拟合所述牙列表面体素集合，具体包括如下步骤：
31)利用Procrustes分析计算三维模型的全局变换使其最大程度地与表面体素集Yc一致；
32)随后利用非刚性配准获取特定的三维模型，使得所述特定的三维模型与初始分割的表面体素集一致；
33)为降低非刚性变换的参数空间，利用嵌入变形方法，将所述特定的三维模型表面的非刚性变形定义为一个稀疏基网格顶点变换的线性组合，三维模型非刚性配准的能量函数定义如式31：
Figure FDA0002152374200000046
式31中，
Figure FDA0002152374200000047
为变形前后表面顶点的位移；变形后的三维模型顶点为
Figure FDA0002152374200000048
其中，权重βij对应顶点
Figure FDA0002152374200000049
与其在基网格上的近邻顶点ηj，并基于两点之间的欧式距离定义为
Figure FDA0002152374200000051
其中参数
Figure FDA0002152374200000052
Tj为第j个基网格顶点上的三维变换矩阵；上述能量函数Ere的第一项用于最小化变形后的三维模型表面与体图像分割获取的表面体素集之间的距离；dc表示变换后的顶点
Figure FDA0002152374200000053
与标签扩散得到的表面体素集Yc之间的距离，
Figure FDA0002152374200000054
其中
Figure FDA0002152374200000055
为表面体素；第二项用于保持三维模型表面平滑，该项最小化顶点
Figure FDA0002152374200000056
与其近邻
Figure FDA0002152374200000057
的空间变换之间的差异；ω为常数系数；ne为三维模型表面顶点个数；
34)通过最小化能量函数Ere，得到三维模型的非刚性变换参数，用于配准拟合所述牙列表面体素集合。
4.如权利要求1所述对锥束CT图像进行牙列分割的方法，其特征是，步骤41)所述柔化约束通过如下过程定义：
410)根据牙齿表面内部的灰度分布，对牙本质赋以属于牙齿的概率，定义基于三维可变形模型属于前景牙齿的概率为式41：
Figure FDA0002152374200000058
式41中，pd(r,θ)为基于三维可变形模型属于前景牙齿的概率；r,θ为体素在极坐标系中的坐标；r1(θ)与r2(θ)表示牙齿的外轮廓以及牙髓腔轮廓，其中r2(θ)＝r1(θ)-a，a为经验的牙本质厚度；
Figure FDA0002152374200000059
参数η用于控制如上概率函数经过牙齿轮廓时函数的形状；
411)针对没有牙髓腔的分层图像，基于三维可变形模型属于前景牙齿的概率定义为式411：
Figure FDA00021523742000000510
式411中，pd(r,θ)为基于三维可变形模型属于前景牙齿的概率函数；r,θ为体素在极坐标系中的坐标；r1(θ)表示牙齿的外轮廓；
Figure FDA0002152374200000061
412)针对牙齿可能出现的多牙根以及多牙尖结构情况，当牙齿在切片上的多轮廓状态时基于三维可变形模型属于前景牙齿的概率表示为式412：
Figure FDA0002152374200000062
式412中，nc为当前牙齿在切片图像中轮廓的个数；pd(r,θ)为式41与式411定义的在坐标(r,θ)处基于三维可变形模型属于前景牙齿的概率；r,θ为体素在极坐标系中的坐标；
Figure FDA0002152374200000063
为扩展到牙齿在切片上的多轮廓状态时基于三维可变形模型属于前景牙齿的概率；
413)利用卡方核的一对多支持向量机生成分类器，所述分类器的输出被归一化到0-1之间，表示体素属于前景牙列的概率pc，柔化约束为式413：
Figure FDA0002152374200000064
式413中，pc表示由体素表观估计的体素属于前景牙列的概率；
Figure FDA0002152374200000065
表示基于三维可变形模型vi属于前景牙齿的概率；r,θ为体素在极坐标系中的坐标；
Figure FDA0002152374200000066
是对体素vi相对于类别s所定义的柔化约束。
5.如权利要求1所述对锥束CT图像进行牙列分割的方法，其特征是，步骤5)所述通过迭代修正过程对牙列分割进行迭代修正，具体包括如下步骤：
51)将柔化约束下的标签扩散和三维模型的非刚性配准能量结合，定义能量如式51：
E(X,T)＝λ1Erw+λ2Esc+λ3Ere (式51)
式51中，E(X,T)是基于三维可变形模型随机游走算法进行牙列分割的总能量函数；Erw与Esc为随机游走以及柔化约束的能量函数；Ere为式31所定义三维模型非刚性配准的能量；λ1，λ2与λ3为常数；
52)在第一个阶段，设定基于标签扩散得到的体图像分割获取牙列表面体素集，将与柔化约束下的随机游走相关的项Erw和Esc从能量函数式51中去掉，通过最小化能量E(X(i-1),T)对三维模型进行变形，用于拟合牙列的表面体素集；此阶段对于三维模型变形参数的求解和步骤3)相同；
53)在第二个阶段，去掉能量函数中与三维模型配准相关的项Ere，最小化能量E(X,T(i -1))，得到感兴趣图像中所有体素的标签；此阶段通过求解一个大规模稀疏的线性系统求解所有体素的标签，和步骤4)所述柔化约束下的随机游走方法相同；
54)迭代优化进行步骤52)和步骤53)，当相邻两步迭代中三维模型更新||Ye(i)-Ye(i-1)||2小于预先给定的阈值时终止迭代，其中Ye(i)与Ye(i-1)分别为第i步与第i-1步三维表面模型，得到感兴趣图像中所有体素的标签，作为迭代修正后的牙列分割结果。
6.如权利要求1～5任一项所述对锥束CT图像进行牙列分割的方法，其特征是，所述得到的牙列表面模型与手工分割得到的表面模型的平均距离误差在0.3mm以下。
CN201610157705.3A 2016-03-18 2016-03-18 一种对锥束ct图像进行牙列分割的方法 Active CN107203998B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610157705.3A CN107203998B (zh) 2016-03-18 2016-03-18 一种对锥束ct图像进行牙列分割的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610157705.3A CN107203998B (zh) 2016-03-18 2016-03-18 一种对锥束ct图像进行牙列分割的方法

Publications (2)

Publication Number Publication Date
CN107203998A CN107203998A (zh) 2017-09-26
CN107203998B true CN107203998B (zh) 2020-04-03

Family

ID=59904364

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610157705.3A Active CN107203998B (zh) 2016-03-18 2016-03-18 一种对锥束ct图像进行牙列分割的方法

Country Status (1)

Country Link
CN (1) CN107203998B (zh)

Families Citing this family (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3462373A1 (en) * 2017-10-02 2019-04-03 Promaton Holding B.V. Automated classification and taxonomy of 3d teeth data using deep learning methods
CN107703029B (zh) * 2017-11-07 2019-05-10 大连理工大学 一种结合ct与pvt计算co2盐水扩散系数的方法
CN108257118B (zh) * 2018-01-08 2020-07-24 浙江大学 一种基于法向腐蚀和随机游走的骨折粘连分割方法
CN108765474A (zh) * 2018-04-17 2018-11-06 天津工业大学 一种针对ct与光学扫描牙齿模型的高效配准方法
CN108670451B (zh) * 2018-05-04 2020-10-13 上海正雅齿科科技股份有限公司 牙齿邻接面修补方法、装置、用户终端及存储介质
CN109978841B (zh) * 2019-03-12 2021-08-03 北京羽医甘蓝信息技术有限公司 基于深度学习的全景片阻生牙识别的方法和装置
CN109903396B (zh) * 2019-03-20 2022-12-16 洛阳中科信息产业研究院 一种基于曲面参数化的牙齿三维模型自动分割方法
CN110232684B (zh) * 2019-06-13 2023-05-23 大连理工大学 一种基于谱分析的三维医学图像自动分割方法
CN113139908B (zh) * 2020-01-17 2022-08-26 北京大学 一种三维牙列分割与标注方法
CN111265317B (zh) * 2020-02-10 2022-06-17 上海牙典医疗器械有限公司 一种牙齿正畸过程预测方法
CN111462159A (zh) * 2020-04-03 2020-07-28 哈尔滨理工大学 一种基于ct图像的变形模型的肝脏自动分割方法
CN112614127A (zh) * 2020-12-31 2021-04-06 北京朗视仪器有限公司 一种基于端到端的交互式三维cbct牙齿图像分割算法
CN113506303B (zh) * 2021-07-27 2024-01-30 四川九洲电器集团有限责任公司 一种交互式牙齿分割方法、装置和处理系统
CN113344950A (zh) * 2021-07-28 2021-09-03 北京朗视仪器股份有限公司 一种深度学习结合点云语义的cbct图像牙齿分割方法
CN113628222A (zh) * 2021-08-05 2021-11-09 杭州隐捷适生物科技有限公司 一种基于深度学习的3d牙齿分割和分类方法
CN114241173B (zh) * 2021-12-09 2023-03-21 电子科技大学 一种牙齿cbct图像三维分割方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104835153A (zh) * 2015-04-30 2015-08-12 天津大学 基于稀疏表示的非刚性表面对齐方法
CN104881681A (zh) * 2015-05-22 2015-09-02 浙江大学 基于混合图模型的图像序列类别标注方法
US9129363B2 (en) * 2011-07-21 2015-09-08 Carestream Health, Inc. Method for teeth segmentation and alignment detection in CBCT volume
CN104933709A (zh) * 2015-06-04 2015-09-23 西安理工大学 基于先验信息的随机游走ct肺组织图像自动分割方法
CN105279762A (zh) * 2015-11-20 2016-01-27 北京航空航天大学 一种口腔软硬组织ct序列与三维网格模型配准方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9129363B2 (en) * 2011-07-21 2015-09-08 Carestream Health, Inc. Method for teeth segmentation and alignment detection in CBCT volume
CN104835153A (zh) * 2015-04-30 2015-08-12 天津大学 基于稀疏表示的非刚性表面对齐方法
CN104881681A (zh) * 2015-05-22 2015-09-02 浙江大学 基于混合图模型的图像序列类别标注方法
CN104933709A (zh) * 2015-06-04 2015-09-23 西安理工大学 基于先验信息的随机游走ct肺组织图像自动分割方法
CN105279762A (zh) * 2015-11-20 2016-01-27 北京航空航天大学 一种口腔软硬组织ct序列与三维网格模型配准方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Quantification of Teeth in CT Images Using Statistcal Shape Model Based on Geometrical Complexity;Mohanmmadraza Soltaninejad等;《RearchGate》;20141128;第1-6页 *

Also Published As

Publication number Publication date
CN107203998A (zh) 2017-09-26

Similar Documents

Publication Publication Date Title
CN107203998B (zh) 一种对锥束ct图像进行牙列分割的方法
Chen et al. Automatic segmentation of individual tooth in dental CBCT images from tooth surface map by a multi-task FCN
US20200350059A1 (en) Method and system of teeth alignment based on simulating of crown and root movement
Montúfar et al. Automatic 3-dimensional cephalometric landmarking based on active shape models in related projections
US20210174543A1 (en) Automated determination of a canonical pose of a 3d objects and superimposition of 3d objects using deep learning
Qian et al. CephaNet: An improved faster R-CNN for cephalometric landmark detection
CN109903396B (zh) 一种基于曲面参数化的牙齿三维模型自动分割方法
CN108205806B (zh) 一种锥束ct图像三维颅面结构的自动解析方法
Zanjani et al. Mask-MCNet: tooth instance segmentation in 3D point clouds of intra-oral scans
CN110246580B (zh) 基于神经网络和随机森林的颅侧面影像分析方法和系统
CN106846346B (zh) 基于关键帧标记的序列ct图像骨盆轮廓快速提取方法
CN107680110B (zh) 基于统计形状模型的内耳三维水平集分割方法
Baka et al. Statistical shape model-based femur kinematics from biplane fluoroscopy
Chung et al. Automatic registration between dental cone-beam CT and scanned surface via deep pose regression neural networks and clustered similarities
CN106780491B (zh) Gvf法分割ct骨盆图像中采用的初始轮廓生成方法
US20080285822A1 (en) Automated Stool Removal Method For Medical Imaging
Cunha et al. A method for segmentation of dental implants and crestal bone
CN114241173B (zh) 一种牙齿cbct图像三维分割方法及系统
Qian et al. An automatic tooth reconstruction method based on multimodal data
CN114187293B (zh) 基于注意力机制和集成配准的口腔腭部软硬组织分割方法
KR101953629B1 (ko) 두개악안면 cbct 영상에서 형상제약 정보를 사용한 하악골 자동 분할 방법 및 장치
Woo et al. Automatic matching of computed tomography and stereolithography data
Antila et al. Segmentation of facial bone surfaces by patch growing from cone beam CT volumes
Chen et al. Detection of Various Dental Conditions on Dental Panoramic Radiography Using Faster R-CNN
CN114972360A (zh) 牙齿的计算机断层扫描图像的分割方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant