技术领域
[0001] 本
发明涉及
计算机视觉和
图像处理技术领域，尤其涉及一种对锥束CT(CBCT)图像进行牙列分割的方法。
背景技术
[0002] 锥束CT(CBCT)图像常用于颌面与
口腔正畸外科辅助手术预测以及牙列对齐计划的制定。CBCT图像可以提供病人特定的解剖结构信息，其中包含完整的牙列信息。而传统的基于光学的方式，例如三维激光扫描以及
立体视觉设备，仅仅能够获取
牙冠的几何信息，其中由于
牙根被埋藏在
牙龈以及颌骨中无法获取其几何形态。在临床口腔外科中，CBCT图像由于其低放射剂量以及较低的采集成本，相对于传统的CT成像技术具有极大的优势。但是，较低的放射剂量以及
信噪比一般会造成CBCT图像
质量较差，常常出现图像退化现象，例如口腔
治疗以及
牙齿种植在牙颌结构中放置的金属物体造成的
线束硬化(beam hardening)问题，有限视域的截断，以及由于牙齿与周围
牙槽骨相似的灰度而造成的模糊的轮廓。同时在
图像采集过程中，病人可能有微小的头部运动，也会造成牙
齿轮廓的模糊。特别是为了诊断牙齿畸形，CBCT在采集过程中通常上下牙列处于咬合状态，这使得上下牙列在咬合部分的分割变得更加困难。
[0003] 近年来有大量技术被应用于从CT以及CBCT图像中进行牙列分割，其中包括使用积分灰度投影、基于
阈值和区域增长的方法、图割法、基于
水平集活动轮廓的方法等。积分灰度投影方法仅仅能获取牙列的粗略分割，例如不同牙齿之间的分割平面。阈值与区域增长方法对于CBCT中可能出现的牙列与牙槽骨以及其它组织之间模糊的轮廓难以有效处理，其分割与真实的牙齿轮廓的一致性较差。基于交互的分割方法例如图割法非常依赖初始前景区域的定义，对于不同的初始前景定义所得到的分割常常存在差异，难以重复分割结果。为了获取可靠的分割，活动轮廓方法诸如水平集技术集成形状以及灰度先验、从切片
图像分割获取的三维形状先验、以及隐参数表达模型，逐层对CBCT图像进行分割。在基于水平集的方法中，对于约束项需要进行精心设计以避免分割过程中可能出现的收缩(shrinkage)以及泄露(leakage)问题。此外，在基于活动轮廓方法进行逐层分割中还可能有误差累计的问题。此外，三维牙列模型被引入分割系统，包括三维形状图集以及三维形状统计模型。图集以及统计形状模型一般来自大量的CT图像，这势必会增加对数据的需求以及
数据处理的负担。
发明内容
[0004] 为了克服上述
现有技术的不足，本发明提供一种对锥束CT(CBCT)图像进行牙列分割的方法，基于可
变形模型的随机游走方法，从CBCT图像中分割得到完整牙列。
[0005] 本发明的原理是：将体图像中的分割问题定义为图结构中的标签扩散问题，其中，预先对少量切片图像进行标注，并使用扩散机制对感兴趣体图像(volume of interest，VOI)中其它
像素进行自动标注。首先，利用原始的随机游走方法对CBCT图像中牙列进行初始分割，基于原始的随机游走的初始分割通常会包含分割错误，例如由于CBCT图像退化出现的线束硬化、图像中模糊的牙列轮廓所造成的标记扩散误差。因此，引入三维可变形模型用于拟合由标记扩散得到的体图像分割，利用组合Logistic函数，并基于经验牙本质厚度定义柔化约束，将柔化约束引入随机游走
算法中更新
体素的类别标签的估计。此外，为了改善从单步柔化约束下的标记扩散所得到的图像分割
精度，引入
迭代修正算法，针对系统的两个方面，即柔化约束下标记扩散所对应的大规模稀疏线性系统和三维模型到图像分割所得到的牙列表面体素的拟合，进行迭代求解。在每步迭代中，三维可变形模型的配准可以看作是对体素标记的正则化，该过程可以有效地消除分割误差。
[0006] 本发明提供的技术方案是：
[0007] 一种对锥束CT图像进行牙列分割的方法，针对锥束CT图像中的感兴趣体图像进行牙列分割，得到牙列表面模型；包括如下步骤：
[0008] 1)针对锥束CT图像中的感兴趣体图像区域定义图结构；所述图结构中，结点集合对应感兴趣体图像区域中所有的体素；边连接集合为感兴趣体图像区域中相邻体素之间的边连接；
[0009] 2)根据牙列分布情况定义类别标签集合，针对步骤1)所述感兴趣体图像中所有的体素，采用随机游走方法进行牙列的初始分割，得到所述所有体素的类别，对所述所有体素进行类别标记；
[0010] 3)根据步骤2)得到的所有体素的类别，利用其中属于牙列的表面体素定义牙列的表面形态，生成初始分割得到的表面体素集；进一步利用三维可变形模型，通过非刚性变形配准拟合所述牙列表面体素集合，得到三维可变形模型的非刚性变换参数，用于拟合所述初始分割得到的表面体素集；
[0011] 4)基于三维可变形模型定义柔化约束，并基于柔化约束下的随机游走方法进行牙列分割，得到感兴趣体图像区域中所有体素的类别标签；
[0012] 5)通过迭代修正过程对牙列分割进行迭代修正，结合步骤4)所述柔化约束下的类别标签扩散和步骤3)所述三维可变形模型的非刚性配准过程，改进单次柔化约束下随机游走方法所得到的体图像分割结果。
[0013] 针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤1)所述图结构中，采用基于上下文的体素特征描述子，通过体素的表观计算相邻体素之间的相似度，用于表示体素的上下文特征差异；所述体素特征描述子向量中的元素对应当前体素的包围图像
块与模式中的
采样体素的包围图像块之间的像素灰度差异累计；图结构中边连接上的权重通过体素上下文特征描述差异和灰度差异的加权组合得到。
[0014] 上述上下文的体素描述子在当前体素的一个半径为 的包围球体中通过随机采样生成一个模式P，所述模式P由采样得到的体素集合表示；所述体素特征描述子向量f(vi)中的元素对应当前体素的包围图像块与模式P中的采样体素的包围图像块之间的像素灰度差异累计，表示为式11：
[0015] f(vi)＝{dp(C(vi),C(vi+γ))|vi+γ∈P} (式11)
[0016] 式11中，f(vi)表示当前体素vi的特征描述子向量；C(vi)表示体素vi的包围图像块，dp为图像块的灰度差异；γ为当前体素与模式P中采样体素之间的位移；C(vi+γ)表示模式中的采样体素的包围图像块；其中，图像块的灰度差异dp定义为式12：
[0017]
[0018] 式12中，dp为图像块的灰度差异；C(vi)与C(vk)分别表示体素vi与vk的包围图像块；δv为
包围体素块中体素的位移；vi+δv属于图像块C(vi)；I(vi+δv)为体素vi+δv的灰度；I(vk+δv)为体素vk+δv的灰度；
[0019] 所述图像块的灰度差异可通过图像的线性卷积算子计算得到：对差异图像V-V(γ)进行卷积，通过式13得到：
[0020] dp＝(V-V(γ))2*G (式13)
[0021] 式13中，dp为图像块的灰度差异；γ为体素vi与模式中一个采样体素之间的位移向量；对体图像进行平移操作得到V(γ)；V-V(γ)为差异图像；G为卷积核；
[0022] 图结构中边连接上的权重通过上下文特征描述差异与灰度差异的加权组合得到，如式14：
[0023]
[0024] 式14中，Wij为体素i与相邻体素j之间边连接的权重；α为常数系数用于调节上下文特征以及体素灰度对于权重计算的影响；Ii与Ij对应第i个与第j个体素的灰度，fi与fj对应第i个与第j个体素的上下文描述；Ψ为图结构的边连接集合。
[0025] 针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤2)采用随机游走方法对感兴趣体图像区域中所有的体素的类别进行标记，得到牙列的初始分割结果；具体通过最小化式21的
能量函数得到感兴趣体图像中体素类别标记：
[0026]
[0027] 式21中，Erw(X)表示基于随机游走算法的图像分割能量；|S|为系统类别标签的个数； 表示第i个结点被分为第s类的概率；系统中每个结点对均对应一个|S|维的概率向量，其中第i个结点所属的类别定义为 表示预先由用户交互定义的结点概率，当第i个结点被用户交互定义为第s类，则 取值为1，否则为0；nL是用户预先交互标注的体素个数；nV为感兴趣体图像中所有体素的个数；μ1与μ2为常数；L为拉普拉斯矩阵；H为对
角线指示矩阵，当感兴趣体图像中的第i个结点被预先标注时，H的元素Hii＝1；X与 表示感兴趣体图像中所有体素对应的类别概率 以及 所组成的矩阵；。
[0028] 其中，将所述式21的目标函数转化为矩阵形式，通过将所述式21的一阶导数设为0，将目标函数求解转化为求解式22所示的大规模稀疏的线性系统：
[0029]
[0030] 式22中，表示单位矩阵；L为拉普拉斯矩阵；H为对角线指示矩阵，当感兴趣体图像中的第i个结点被预先标注时，H的元素Hii＝1；X与 表示感兴趣体图像中所有体素对应的类别概率 以及 所组成的矩阵；
[0031] 通过求解式22的线性系统，得到感兴趣体图像中所有体素的类别，作为牙列的初始分割结果。
[0032] 针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤3)所述三维可变形模型用于处理锥束CT图像自身的局部退化造成的边界混淆问题；所述通过非刚性变形配准拟合所述牙列表面体素集合，具体包括如下步骤：
[0033] 31)利用Procrustes分析计算三维模型的全局变换使其最大程度地与表面体素集Yc一致；
[0034] 32)随后利用非刚性配准获取特定的三维模型，使得所述特定的三维模型与初始分割的表面体素集一致；
[0035] 33)为降低非刚性变换的参数空间，利用嵌入变形方法，将所述特定的三维模型表面的非刚性变形定义为一个稀疏基网格
顶点变换的线性组合，三维模型非刚性配准的能量函数定义如式31：
[0036]
[0037] 式31中， 为变形前后表面顶点的位移；变形后的三维模型Ye的顶点为 其中，权重βij对应顶点 与其在基网格B上的近邻顶点ηj，并基于两点
之间的欧式距离定义为 其中参数 Tj为第j个基网格
顶点上的三维变换矩阵；能量函数Ere的第一项用于最小化变形后的三维模型表面与体图像分割获取的表面体素集之间的距离；dc表示变换后的顶点 与标签扩散得到的表面体素集Yc之间的距离, 其中 为表面体素；第二项用于保持三维模型表面平
滑，该项最小化顶点 与其近邻 的空间变换之间的差异；ω为常数系数，ne为三维模型表面顶点个数。
[0038] 34)通过最小化能量函数Ere，得到三维模型的非刚性变换参数，用于配准拟合所述牙列表面体素集合。
[0039] 针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤4)所述基于柔化约束下的随机游走方法进行牙列分割，包括如下步骤：
[0040] 41)通过设定与体图像分割结果拟合后的三维模型表面内部体素更有可能属于牙列，根据体素结点概率定义柔化约束；
[0041] 42)设定柔化约束下的感兴趣体图像的标签扩散能量为式43：
[0042] Erwsc(X)＝λ1Erw(X)+λ2Esc(X) (式43)
[0043] 其中Erw与Esc式21式42所定义的随机游走以及柔化约束的能量。λ1与λ2为常数系数。
[0044] 将式43中的能量函数的一阶导数设为0，得到关于体素结点概率向量的线性系统如式44，
[0045]
[0046] 其中L为拉普拉斯矩阵；H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；X与 表示VOI中所有体素对应的类别概率 以及用户定义的先验概率 所组成的矩阵。 对应感兴趣区域中所有体素的柔化约束所构成的矩阵；λ1与λ2为常数系数,实验中分别设置为1、0.1；μ1与μ2为常数，本实验中分别设置为1与0.001；表示单位矩阵。
[0047] 43)利用共轭梯度方法求解线性系统，该线性系统需要求解|S|-1次，其中|S|为系统中类别标签的个数；通过求解线性系统得到感兴趣体图像中所有体素的类别标签。
[0048] 进一步地，步骤41)采用以下步骤定义所述柔化约束：
[0049] 410)根据牙齿表面内部的灰度分布，对牙本质赋以更大的属于牙齿的概率，定义基于三维可变形模型属于前景牙齿的概率为式41：
[0050]
[0051] 式41中，pd(r,θ)为基于三维可变形模型属于前景牙齿的概率；r1(θ)与r2(θ)表示牙齿的外轮廓以及
牙髓腔轮廓，其中r2(θ)＝r1(θ)-a，a为经验的牙本质厚度；参数η用于控制如上概率函数经过牙齿轮廓时函数的形状；
[0052] 411)针对没有牙髓腔的分层图像，基于三维可变形模型属于前景牙齿的概率定义为式411：
[0053]
[0054] 式411中，
[0055] 412)针对牙齿可能出现的多牙根以及多
牙尖结构情况，牙齿在切片上呈多轮廓状态时基于三维可变形模型属于前景牙齿的概率表示为式412：
[0056]
[0057] 式412中，nc为当前牙齿在切片图像中轮廓的个数；pd为单外轮廓的基于三维可变形模型属于前景牙齿的概率(式41与式411)；r,θ为体素在极
坐标系中的坐标；
[0058] 413)利用卡方核的一对多
支持向量机生成分类器，所述分类器的输出被归一化到0-1之间，表示基于体素表观属于前景(牙列)的概率pc，柔化约束定义为式413：
[0059]
[0060] 式413中，pc表示由体素vi表观估计的体素属于前景(牙列)的概率； 表示在体素vi基于三维可变形模型属于前景牙齿的概率；r,θ为体素在极坐标系中的坐标； 是对体素vi相对于类别s所定义的柔化约束。
[0061] 针对上述对锥束CT图像进行牙列分割的方法，进一步地，步骤5)所述通过迭代修正过程对牙列分割进行迭代修正，具体包括如下步骤：
[0062] 51)将柔化约束下的标签扩散和三维模型的非刚性配准能量结合，定义能量如式51：
[0063] E(X,T)＝λ1Erw+λ2Esc+λ3Ere (式51)
[0064] 式51中，E(X,T)是基于三维可变形模型随机游走算法进行牙列分割的总能量函数；Erw与Esc为式21与式42所定义的随机游走以及柔化约束的能量函数；Ere为式31所定义三维模型非刚性配准的能量；λ1，λ2与λ3为常数，实验中分别设置为1、0.1、1；；
[0065] 52)在第一个阶段，设定基于标签扩散得到的体图像分割获取牙列表面体素集，由于表面体素集由第i-1次迭代中的CBCT图像分割得到，所以将与柔化约束下的随机游走相关的项Erw和Esc从能量函数式51中去掉，通过最小化能量E(X(i-1),T)对三维模型进行变形，用于拟合牙列的表面体素集；此阶段对于三维模型变形参数的求解和步骤3)相同；
[0066] 53)在第二个阶段，去掉能量函数中与三维模型配准相关的项Ere，最小化能量E(X,T(i-1))，得到感兴趣图像中所有体素的标签；此阶段通过求解一个大规模稀疏的线性系统求解所有体素的标签，和步骤4)所述柔化约束下的随机游走方法相同；
[0067] 54)迭代优化进行步骤52)和步骤53)，当相邻两步迭代中三维模型更新||Ye(i)-Ye(i-1)||2小于预先给定的阈值时终止迭代，得到感兴趣图像中所有体素的标签，作为迭代修正后的牙列分割结果。其中Ye(i)与Ye(i-1)分别为第i步与第i-1步三维表面模型。
[0068] 针对上述对锥束CT图像进行牙列分割的方法，进一步地，所述得到的牙列表面模型与手工分割得到的表面模型的平均距离误差在0.3mm以下。
[0069] 与现有技术相比，本发明的有益效果是：
[0070] 本发明提供一种对锥束CT(CBCT)图像进行牙列分割的方法，基于可变形模型的随机游走方法，从CBCT图像中分割得到完整牙列。本发明方法结合半监督的随机游走算法和三维可变形模型配准定义的柔化约束，获取牙列的可靠分割。其中，利用三维可变形模型引入体图像的柔化约束，可以看作为正则化项用于处理基于半监督标签扩散的体图像分割中的噪声。本发明还采用迭代修正方法，对柔化约束下标记扩散以及三维模型对表面体素集的拟合问题，进行迭代求解，可有效地消除分割误差，改进单次标签扩散所获取的牙列分割，提高分割结果的精度。
[0071] 利用本发明提供的方法，可以对医学临床CBCT图像进行自动的牙列分割，分割得到的三维牙列模型可用于颌面与
口腔正畸外科辅助手术预测以及牙列对齐计划的制定，该牙列分割可满足口腔临床对于精度的要求。度量自动分割所得到的牙列表面模型与手工分割对应的表面模型的平均距离，该距离误差在0.3mm以下。
附图说明
[0073] 图2是本发明中采用的基于上下文的体素描述子的结构示意图；
[0074] 其中，vi为当前体素；为当前体素的一个包围球体的半径；γ为当前体素与模式P中采样体素之间的位移；C(vi)为当前体素vi的包围图像块；C(vi+γ)为模式中的采样体素的包围图像块。
具体实施方式
[0075] 下面结合附图，通过
实施例进一步描述本发明，但不以任何方式限制本发明的范围。
[0076] CBCT图像是三维体图像，由一系列的二维图像组成，本说明将其中二维图像称作切片图像(分层图像)。本发明实施例针对医学临床CBCT图像，基于CBCT图像中感兴趣体图像区域中定义的图结构和少量用户交互标注的分层图像(切片图像)，对CBCT图像中的牙列进行分割。其中，采用三维可变形模型定义柔化约束，通过基于柔化约束下的随机游走算法更新体图像中的牙列分割，再以迭代修正方法获取可靠的分割结果。
[0077] 图1是本发明提供方法的流程框图。具体地，本发明提供方法包括如下步骤：
[0078] 步骤1：设定CBCT图像中的感兴趣体图像(VOI)，定义VOI中的图结构；
[0079] 本方法以一种半监督的方式对CBCT图像进行牙列分割，其中，将用户所
指定的少量切片图像中的标记扩散到感兴趣体图像的整个区域。该扩散过程在VOI中定义的图结构中进行。针对VOI定义的图结构中，结点集合对应VOI中所有的体素，边连接集合为相邻体素之间的边连接。
[0080] 一般地，图结构的边连接集合包含两类，一类对应切片图像内部的边连接，该连接对应切片图像中规则的像素网格；另一类对应切片图像之间的边连接，依据由稠密光流算法获取的切片图像之间的对应体素建立边连接。值得注意的是，由稠密光流算法建立的对应并非一一对应，即第i层切片图像中的一个结点可能与第(i+1)层切片图像中不止一个结点对应。该情况常常发生在牙齿区域分叉的切片图像，例如磨牙的多牙根区域，或者牙冠中的多牙尖区域。
[0081] 系统基于体素的表观计算相邻体素之间的相似度(为上下文特征描述差异与灰度差异的加权组合)。使用体素的灰度I计算体素的相似度是一种较直观的方式，但是仅仅利用体素灰度定义体素相似度对于低
辐射剂量采集的CBCT图像并不适合。这是由于CBCT图像中可能存在退化区域，其中仅由灰度相似并不能说明对应体素相似。因此系统引入基于上下文的体素描述子。
[0082] 图2是基于上下文的体素描述子；其中，当前体素的一个半径为 的包围球体中，γ为当前体素与模式P中采样体素之间的位移；通过随机采样生成模式P；特征描述子向量f(vi)中的元素对应当前体素vi的包围图像块C(vi)与模式中的采样体素的包围图像块C(vi+γ)之间的像素灰度差异累计。
[0083] 为了降低计算的复杂度，该上下文的体素描述子并未对当前体素的所有近邻体素进行计算，而是在当前体素的一个半径为 的包围球体中通过随机采样生成一个模式P。该模式由采样得到的体素集合表示。特征描述子向量f(vi)中的元素对应当前体素的包围图像块与模式中的采样体素的包围图像块之间的像素灰度差异累计，表示为式11：
[0084] f(vi)＝{dp(C(vi),C(vi+γ))|vi+γ∈P} (式11)
[0085] 式11中，f(vi)表示当前体素vi的特征描述子向量；C(vi)表示体素vi的包围图像块，dp为图像块的灰度差异；γ为当前体素与模式P中采样体素之间的位移；C(vi+γ)表示模式中的采样体素的包围图像块。其中，图像块的灰度差异dp定义为式12：
[0086]
[0087] 式12中，dp为图像块的灰度差异；C(vi)与C(vk)分别表示体素vi与vk的包围图像块；δv为包围体素块中体素的位移；vi+δv属于图像块C(vi)；I(vi+δv)为体素vi+δv的灰度；I(vk+δv)为体素vk+δv的灰度；
[0088] 图像块的灰度差异可以看作图像块之间对应体素灰度差异的求和，该计算可以通过图像的线性卷积算子实现。如果体素vi与模式中一个采样体素之间的位移向量为γ，则可预先对体图像进行平移操作得到V(γ)。图像块的灰度差异dp则可以对差异图像V-V(γ)进行卷积得到，即通过式13得到：
[0089] dp＝(V-V(γ))2*G (式13)
[0090] 式13中，dp为图像块的灰度差异；G为卷积核,V(γ)为对原始体图像V进行平移γ后得到的体图像。
[0091] 结合体素的上下文特征描述，图结构中边连接上的权重定义为上下文特征描述差异与灰度差异的加权组合，如式14：
[0092]
[0093] 式14中，Wij为体素i与相邻体素j之间边连接的权重；α为常数系数用于调节上下文特征以及体素灰度对于权重计算的影响；Ii与Ij对应第i个与第j个体素的灰度，fi与fj对应第i个与第j个体素的上下文描述。Ψ为图结构的边连接集合；本实施例中，系统预先将CBCT图像的灰度范围变换到[0,255]区间，其中α被设置为0.1。通过式14计算得到体素的相似度矩阵。
[0094] 步骤2：采用随机游走算法进行牙列的初始分割；
[0095] 系统根据牙列分布情况定义标签集合S，其中对不同的牙齿分配不同的标签，正常的牙列包含32颗牙齿，其中第三磨牙，即
智齿通常未发育完全，系统将其排除在外，还余下28颗牙齿，因而加上背景类别，系统的标签集合中将定义29类标签。但在实际CBCT数据中，由于牙齿缺失等问题，实际CBCT图像中包含的类别可能少于该类别个数。在初始牙列分割中，系统采用原始的随机游走算法对感兴趣体图像区域进行标记。VOI中体素标记可以通过最小化如下的能量函数(式21)得到。
[0096]
[0097] 式21中，Erw(X)表示基于随机游走算法的图像分割能量；|S|为系统类别标签的个数； 表示第i个结点被分为第s类的概率；系统中每个结点均对应一个|S|维的概率向量，其中第i个结点所属的类别定义为 表示预先由用户交互定义的结点概率，当第i个结点被用户交互定义为第s类，则 取值为1，否则为0；nL是用户预先交互标注的体素个数；nV为VOI中所有体素的个数；μ1与μ2为常数，本实验中分别设置为1与0.001。上述目标函数(式21)可以转化为矩阵形式，其中L为拉普拉斯矩阵，在给定步骤1)中通过公式14得到的相似度矩阵W后，拉普拉斯矩阵L＝D-W；其中，对角线矩阵D中元素定义为 H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；X与 表示VOI中所有体素对应的类别概率 以及 所组成的矩阵。通过将上面的能量函数的一阶导数设为0，该能量函数可以转化为求解如下大规模稀疏的线性系统，表示为式22：
[0098]
[0099] 式22中， 表示单位矩阵；L为拉普拉斯矩阵；H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；μ1与μ2为常数，本实验中分别设置为1与0.001。给定用户定义的形状先验 可通过求解上面的线性系统(式22)得到VOI中所有体素的类别，即得到牙列的初始分割结果。
[0100] 步骤3：利用三维可变形模型，通过非刚性变形配准拟合步骤2进行图像分割所得到的牙列表面体素集合，得到三维可变形模型的非刚性变换参数；
[0101] 步骤2得到VOI中所有体素的类别，利用其中属于牙列表面的体素，可以定义牙列的表面形态；考虑到CBCT图像中可能存在的退化现象，以及处于咬合状态上下颌
接触的牙齿，基于原始的随机游走技术获取的初始分割与真实牙齿轮廓之间会存在较大的差异。系统引入三维可变形模型处理CBCT图像自身的局部退化问题所造成的边界混淆问题。三维可变形模型来自对实体牙齿模型的三维光学扫描，其具有良好的拓扑结构定义。从CBCT图像分割可以定义牙列表面的体素集合，即属于牙列的体素如果其近邻体素属于背景，则该体素被放入表面体素集Yc。对三维可变形模型施加非刚性变换，以拟合初始分割得到的表面体素集。与体图像分割结果拟合后的三维模型表面内部体素更有可能属于牙列，并基于该假设定义柔化约束。
[0102] 三维模型非刚性拟合分为两个阶段，首先利用Procrustes分析计算三维模型的全局变换使其最大程度地与表面体素集Yc一致；随后利用非刚性配准获取特定病人的三维模型，使其与初始分割的表面体素集一致。为了降低非刚性变换的参数空间，系统利用嵌入变形技术，将三维模型表面的非刚性变形定义为一个稀疏基网格顶点变换的线性组合，其中变形后的三维模型顶点为 其中，权重βij对应顶点 与其在基网格B上的近邻顶点ηj ，并基于两点之间的欧式距离定义为 其中参数Tj为第j个基网格顶点上的三维变换矩阵，三维模型非刚性配准的能量
函数定义如式31：
[0103]
[0104] 式31中， 即变形前后表面顶点的位移；能量函数Ere的第一项用于最小化变形后的三维模型表面与体图像分割获取的表面体素集之间的距离；dc表示变换后的顶点 与标签扩散得到的表面体素集Yc之间的距离, 其中 为表
面体素；第二项用于保持三维模型表面平滑，该项最小化顶点 与其近邻 的空间变换之间的差异。ω为常数系数，ne为三维模型表面顶点个数。通过最小化能量函数Ere可以得到三维模型的非刚性变换参数。
[0105] 步骤4：基于柔化约束下的随机游走算法进行牙列分割，得到VOI中所有体素的类别标签，即哪些体素属于牙列(前景)，哪些体素属于背景；
[0106] 可以直观看到，与体图像分割得到的表面体素集配准后的三维模型表面内部的体素更有可能属于牙列。但是，在体图像中牙齿表面内部体素并不具有均匀的灰度，其中牙本质具有较亮的灰度，而牙髓腔灰度较低，考虑到牙齿表面内部的灰度分布，系统对牙本质赋以更大的属于牙齿的概率。由于三维模型来自对牙齿模型的光学扫描，该模型不包括内部的牙髓腔的表面，利用经验的牙本质厚度定义a，以及组合Logistic函数依据步骤3得到的变形后的三维模型表面定义柔化约束。极坐标系的中心定义为变形后三维模型与切片图像相交轮廓的中心，令r1(θ)与r2(θ)表示牙齿的外轮廓以及牙髓腔轮廓，其中r2(θ)＝r1(θ)-a。基于三维可变形模型属于前景牙齿的概率pd定义为式41：
[0107]
[0108] 式41中，pd(r,θ)为基于三维可变形模型属于前景牙齿的概率；r,θ为体素在极坐标系中的坐标；r1(θ)与r2(θ)表示牙齿的外轮廓以及牙髓腔轮廓，其中r2(θ)＝r1(θ)-a，a为经验的牙本质厚度； 参数η用于控制如上概率函数经过牙齿轮廓时函数的形状。在实验中，参数η的值设为0.99。
[0109] 如果考虑没有牙髓腔的分层图像，基于三维可变形模型属于前景牙齿的概率定义为式411：
[0110]
[0111] 其中，pd(r,θ)为基于三维可变形模型属于前景牙齿的概率；r,θ为体素在极坐标系中的坐标；r1(θ)表示牙齿的外轮廓； 其中参数η用于控制如上概率函数经过牙齿轮廓时函数的形状。在实验中，参数η的值设为0.99。
[0112] 上述的基于三维可变形模型属于前景牙齿的概率针对切片图像中的单轮廓定义，考虑到牙齿可能出现的多牙根以及多牙尖结构，该定义可以直接扩展到牙齿在切片上的多轮廓状态，表示为式412：
[0113]
[0114] 其中，nc为当前牙齿在切片图像中轮廓的个数；pd(r,θ)为式41与式411定义的在坐标(r,θ)处基于三维可变形模型属于前景牙齿的概率； 为扩展到牙齿在切片上的多轮廓状态时基于三维可变形模型属于前景牙齿的概率。
[0115] 由于三维模型变形后与牙齿真实表面之间可能存在空隙，系统还利用体素的表观定义体素属于牙列的分类器。从交互标注的切片图像中学习该分类器，利用卡方(chi-square)核的一对多(one-vs-rest)支持向量机(SVM)生成分类器，该分类器的输出被归一化到0-1之间，表示体素属于前景(牙列)的概率pc。柔化约束定义为式413：
[0116]
[0117] 式413中，pc表示由体素表观估计的体素vi属于前景(牙列)的概率； 表示基于三维可变形模型vi属于前景牙齿的概率；r,θ为体素vi在极坐标系中的坐标。 是对体素vi相对于类别s所定义的柔化约束。
[0118] 柔化约束的能量项定义为式42：
[0119]
[0120] 其中， 表示第i个结点被分为第s类的概率；X表示VOI中所有体素对应的类别概率 所组成的矩阵；|S|表示系统中类别个数；nV为感兴趣区域中体素的个数； 为式413所定义的柔化约束； 对应感兴趣区域中所有体素的柔化约束所构成的矩阵；给定柔化约束，组合式21与式42，将感兴趣体图像的标签扩散能量定义为式43：
[0121] Erwsc(X)＝λ1Erw+λ2Esc (式43)
[0122] 其中Erw与Esc式21式42所定义的随机游走以及柔化约束的能量。λ1与λ2为常数系数。将该能量函数(式43)的一阶导数设为0，可得到关于体素结点概率向量的线性系统，利用共轭梯度方法求解如下线性系统(式44)。
[0123]
[0124] 其中L为拉普拉斯矩阵；H为对角线指示矩阵，当VOI中的第i个结点被预先标注时，H的元素Hii＝1；X与 表示VOI中所有体素对应的类别概率 以及用户定义的先验概率所组成的矩阵。 对应感兴趣区域中所有体素的柔化约束所构成的矩阵；λ1与λ2为常数系数,实验中分别设置为1、0.1；μ1与μ2为常数，本实验中分别设置为1与0.001；表示单位矩阵。由于xi是|S|维的向量，其中|S|为系统中类别标签的个数，该线性系统需要求解|S|-1次，其中 求解线性系统可以得到VOI中所有体素的类别标签。
[0125] 步骤5：对牙列分割迭代修正；
[0126] 为了改进单次柔化约束下随机游走算法所得到的体图像分割，系统引入迭代修正过程，对柔化约束下标记扩散，以及三维模型对表面体素集的拟合问题，进行迭代求解，以改进单次标签扩散所获取的牙列分割。
[0127] 将柔化约束下的标签扩散以及三维模型的非刚性配准能量结合，定义能量如下：
[0128] E(X,T)＝λ1Erw+λ2Esc+λ3Ere (式51)
[0129] 式51中，E(X,T)是基于三维可变形模型随机游走算法进行牙列分割的总能量函数；Erw与Esc为式21式42所定义的随机游走以及柔化约束的能量函数；Ere为式31所定义三维模型非刚性配准的能量；λ1，λ2与λ3为常数，实验中分别设置为1、0.1、1。
[0130] 该能量函数的优化可以分解为两个子问题。在第一个阶段，给定基于标签扩散得到的体图像分割定义牙列表面体素集，对三维模型进行变形，用于拟合更新后的表面体素集。该变形通过最小化能量E(X(i-1),T)得到，其中表面体素集由第i-1次迭代中的CBCT图像分割得到，因而与柔化约束下的随机游走相关的项，即Erw与Esc可以从能量函数中去掉。对于三维模型变形参数的求解与步骤3)的描述一致。
[0131] 在第二个阶段，最小化能量E(X,T(i-1))以得到VOI中所有体素的标签。由于在标签扩散过程中使用的柔化约束来自第i-1步迭代中的三维模型配准，因而能量函数中与三维模型配准相关的项，即Ere可被去掉，该阶段对应步骤4)中柔化约束下的随机游走算法，通过求解一个大规模稀疏的线性系统求解所有体素的标签。
[0132] 随着两个子问题的优化迭代进行，变形后的三维物体表面与基于标签扩散所得到的表面体素集会越来越接近。当相邻两步迭代中三维模型更新||Ye(i)-Ye(i-1)||2(其中Ye(i)与Ye(i-1)分别为第i步与第i-1步三维表面模型)小于预先给定的阈值时，终止迭代，得到VOI中所有体素的标签，作为迭代修正后的牙列分割结果。
[0133] 为了验证基于三维可变形模型随机游走算法对牙列进行分割的精度，利用Dice测度计算自动分割与手工分割的相似度，该相似度大于0.95。同时，度量自动分割所得到的牙列表面模型与手工分割对应的表面模型的平均距离，该距离误差在0.3mm以下，可满足口腔临床对于精度的要求。
[0134] 利用本发明的方法，可以对临床采集CBCT图像进行牙列分割。该方法结合半监督的随机游走算法，以及由三维可变形模型配准定义的柔化约束获取牙列的可靠分割。其中利用三维可变形模型引入柔化约束，也可以看作为正则化项用于处理基于半监督标签扩散的体图像分割中的噪声。该发明引入迭代修正方法，其中对系统的两个子问题，即柔化约束下标记扩散，以及三维模型对表面体素集的拟合问题，进行迭代求解以改进单次标签扩散所获取的牙列分割。
[0135] 需要注意的是，公布实施例的目的在于帮助进一步理解本发明，但是本领域的技术人员可以理解：在不脱离本发明及所附
权利要求的精神和范围内，各种替换和
修改都是可能的。因此，本发明不应局限于实施例所公开的内容，本发明要求保护的范围以权利要求书界定的范围为准。
