锥形束计算机断层扫描（Cone Beam Computed Tomography，CBCT）即锥形束CT，因具有扫描时间短、辐射剂量小、设备轻便等优点而被广泛应用于口腔诊断领域.通过对CBCT所得二维图像进行重建，实现牙齿及牙槽骨等骨性结构的三维显示，这将大大减少医生的诊断难度，同时为治疗方案的确定和虚拟手术的实现提供便利.但由于CBCT具有辐射性，磁共振成像（Magnetic Resonance Imaging，MRI）逐渐成为了更安全的口腔检测成像手段.因为成像原理的不同，牙齿及牙槽骨在CBCT图像中为高密度影像，而在MRI中为低密度影像，且由于MRI主要适用于软组织成像，骨性组织分辨率较差，因此虽能基本完整地展现牙齿与牙槽骨的形态与结构，但依然存在相邻牙齿间黏连、牙根处拓扑结构改变、牙齿与牙槽骨皮质连接处等复杂结构无法清晰分辨的问题.利用MRI检测下颌骨病，并将牙齿从图像分割出来后，可以对牙槽骨和下颌骨病变作出更准确地诊断.

把牙齿和牙槽骨分别从二维图像中分割出来，是实现骨性结构三维重建的基本前提；并且分割的精确度将直接影响重建的效果.因此，如何实现二维图像中牙齿及牙槽骨的精确分割就成为一个关键问题.水平集方法由Osher和Sethian^[1]在1988年提出，因其在描述复杂拓扑结构上的独特优势，水平集方法近年来被广泛应用于图像识别、检测等领域.国内外研究人员对于改进水平集算法有过许多研究.Gao等^[2]首次将基于边界、基于区域及基于先验形状约束能量项加入到水平集中，提出了一种混合的水平集模型.Zhang等^[3]基于反应扩散正则化提出了一种新的水平集演化方法，佐以两步裂项法，可大大降低计算的复杂性.由于单一的能量模型很难准确分割目标，Wang等^[4]提出了一种由多尺度局部似然图像拟合（Local Likelihood Image Fitting，LLIF）能量项、自适应先验形状约束能量项和反应扩散（Reaction Diffusion，RD）正则化能量项组成的混合水平集模型，能够有效提高分割精度.

Li等^[5]于2010年提出一种避免初始化水平集函数的距离正则化水平集演化（Distance Regularized Level Set Evolution，DRLSE）模型，其通过加入距离正则项（Distance Regularized Term）使零水平集函数在向目标物体边缘移动的过程中，保持了其符号距离特性.但通过观察其距离正则项所对应的势阱函数，发现该势阱函数在水平集函数梯度$\nabla \phi $为0时会出现演化速度变至无穷的现象，这将直接导致水平集函数越过目标物体的某些较弱的边缘，出现曲线侵入.而牙齿的牙冠部分通常情况下存在紧贴或挤压，这种情况在CBCT与磁共振图像上即表现为弱边缘，故无法将该种势阱函数所得到的距离正则项直接应用于对CBCT与磁共振图像中牙齿的分割.基于此，本文在该势阱函数的基础上进行了改进，提出了一种新的势阱函数.该势阱函数所产生的正则项克服了原势阱函数对弱边缘不敏感的缺点、并加快水平集函数演化速度，更适合用于对牙齿牙槽骨图像的分割.

水平集方法能够简明地表示复杂拓扑结构的基础思想在于其将n维曲面的演化问题转化为n+1维空间中隐式方程解集合，即水平集的演化.因此二维曲线运动将转化为三维曲面运动，即该曲线每一个时刻的变化均可利用三维曲面的零水平面来表达.图 1为水平集方法示意图.

但利用水平集进行曲线演化的缺点是要进行周期性初始化，为避免复杂的初始化过程，Li等^[5]提出了归一化符号距离函数.通过在能量泛函中构造势阱函数决定的正则项来保证每一次演化后符号距离函数的值满足一定条件，从而避免了对水平集函数即该符号距离函数的重新初始化.水平集函数定义为：

(1)$\phi (x, y) = \left\{ \begin{gathered} - d, {\rm{ }}(x, y) \in {{\mathit{\Omega}} _{inside}} \\ d, {\rm{ }}(x, y) \in {{\mathit{\Omega}} _{outside}} \\ \end{gathered} \right.$

其中，d为正常数，Ω表示初始轮廓区域.初始轮廓内部设定为负值，初始轮廓外部设定为正值.为驱动该水平集函数表示的曲线向着感兴趣区域的边界移动，定义水平集能量泛函$E{\rm{(}}\phi {\rm{)}}$为：

(2)$E{\rm{(}}\phi {\rm{) = }}\mu {R_p}{\rm{(}}\phi {\rm{) + }}\lambda {L_g}{\rm{(}}\phi {\rm{)}} + \alpha {A_g}{\rm{(}}\phi {\rm{)}}$

(2) 式第一项中，${R_p}(\phi )$为距离正则项，由(3)式表示：

(3)${R_p}(\phi ) \triangleq \int_ {\mathit{\Omega }} {P(\left| {\nabla \phi } \right|){\rm{d}}x} $

其中P(·)为势阱函数，$\nabla $表示梯度算子，$\nabla \phi $为水平集函数梯度，根据势阱函数最小点处$\nabla \phi $的取值，分为单势阱函数和双势阱函数，在DRLSE模型中，选用了双势阱函数.由于构造势函数的最终目的是保证符号距离的值，因此采用这种势函数的正则项被称为距离正则项.由于对势阱函数取积分后所得的距离正则项可驱动曲线演化，所以针对不同的图像情况，选择恰当的势阱函数是得到准确结果的关键.

(2) 式第二项中，${L_g}(\phi )$为长度项，由(4)式表示：

(4)${L_g}(\phi ) \triangleq \int_\Omega {g\delta (\phi )\left| {\nabla \phi } \right|{\rm{d}}x} $

其中$\delta ( \cdot )$是Dirac函数，此处采用(5)式${\delta _\varepsilon }( \cdot )$逼近$\delta ( \cdot )$，$\varepsilon $常取1.5^[5]：

(5)${\delta _\varepsilon }(x) = \left\{ \begin{array}{l} \frac{1}{{2\varepsilon }}\left[ {1 + \cos (\frac{{\pi x}}{\varepsilon })} \right]{\rm{, }}\left| x \right| \leqslant \varepsilon \\ {\rm{0, }}\quad\quad\quad\quad\quad\quad \left| x \right| > \varepsilon \\ \end{array} \right.$

边缘指示函数g用(6)式表示，其中${I_0}$为待分割图像，${G_\sigma }$为标准偏差为σ的高斯核，${G_\sigma }*{I_0}$为高斯核与待分割图像作卷积平滑后的图像，$\nabla $为微分算子符号：

(6)$g = \frac{1}{{1 + {{\left| {\nabla {G_\sigma } * {I_0}} \right|}^2}}}$

在图像真实的边界处，梯度指向牙齿内部，然而相邻牙齿、边界的梯度和牙齿内部边界（牙釉质和牙本质的边界）的梯度都是指向牙齿外部.水平集符号距离函数的梯度方向总是指向轮廓的内部，即与真实边界的梯度方向夹角成锐角.因此，只有当轮廓边界的梯度方向和水平集符号距离函数的梯度方向指向同侧时，得到的轮廓为目标边界^[6].根据文献[7]的分析，g的范围在0~1之间，当接近目标区域边缘时取值较小，接近于0；当处于相同背景的区域时，取值较大，接近于1.

(2) 式第三项中，${A_g}(\phi )$为面积项，由(7)式表示：

(7)${A_g}(\phi ) \triangleq \int_\Omega {gH( - \phi ){\rm{d}}x} $

其中g为边缘指示函数，$H( \cdot )$为Heaviside函数，此处采用(8)式${H_\varepsilon }( \cdot )$逼近$H( \cdot )$，$\varepsilon $常取1.5，${H_\varepsilon }( \cdot )$与${\delta _\varepsilon }( \cdot )$的关系为${H_\varepsilon }^\prime (x) = {\delta _\varepsilon }(x)$：

(8) $H_{\varepsilon}(x)=\left\{\begin{array}{ll}\frac{1}{2}\left[1+\frac{x}{\varepsilon}+\frac{1}{\pi} \sin \left(\frac{\pi x}{\varepsilon}\right)\right], & |x| \leq \varepsilon \\1, & x>\varepsilon \\0, & x<-\varepsilon\end{array}\right.$

(2) 式中，$\mu $为正则权重系数；$\lambda $为长度权重系数；$\alpha $为面积权重系数，当$\alpha $为正时，轮廓向内收缩，$\alpha $为负时，轮廓向外扩张，从而驱动零水平集曲线向边界演化^[8].在曲线演化的过程中，若曲线尚未到达待分割图像目标区域的轮廓边界，能量泛函$E(\phi )$将驱动曲线向边界靠拢；当水平集函数位于目标位置时，能量泛函$E(\phi )$取得最小值，曲线停止演化，得到分割结果.

正则项的作用不仅是平滑水平集函数，更重要的一点是令至少在零水平集附近区域内满足$|\nabla \phi | = 1$，这样才可以保证曲线演化过程中相关计算的准确性.在势阱方程中用s表示$\nabla \phi $，可以通过构造在$s = 1$处取得最小值的势函数$p(s)$来实现以上目的.根据势阱函数最小点处s的取值，分为单势阱函数和双势阱函数.单势阱函数表达式简单，但在某些环境下，为得到最小化能量而衍生出的水平集演化会产生很严重的副作用，故多采用双势阱函数构造正则项.

Li等^[5]在DRLSE模型中采用了一种双势阱函数来保持水平集函数的符号距离特性即$|\nabla \phi | = 1$，在该模型中，势阱函数$p(s)$以及扩散速率函数${d_p}(s)$表达如下：

(9)$p(s) = \left\{ \begin{array}{l} \frac{1}{{{{(2{\rm{ \mathsf{ π} }})}^2}}}[1 - \cos (2{\rm{ \mathsf{ π} }}s)], {\rm{ }}s \leqslant 1 \\ \frac{1}{2}{(s - 1)^2}{\rm{, }}\quad\quad\quad\quad s > 1 \\ \end{array} \right.$

(10)${d_p}(s) = \frac{{p'(s)}}{s} = \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{2{\rm{ \mathsf{ π} }}s}}{\rm{sin(2 \mathsf{ π} }}s), {\rm{ }}s \leqslant 1} \\ {1 - \frac{1}{s}{\rm{, }}\quad\quad s > 1} \end{array}} \right.$

其中，s表示水平集梯度$\nabla \phi $.不难发现，该双势阱函数在零势阱附近的演化速度接近于1，达到演化速度的最大值，使得水平集急速向前演化，最终使弱边缘处的s减小至0，导致零水平集侵入被分割目标内部，无法实现对弱边缘的准确检测；这种不足使得将其运用在图像分割的算法中时条件受限.

针对Li提出的双势阱的不足，孙超等^[9]提出了V势阱函数，该函数统一了演化速度的数学表达形式，在避免零势阱修正的同时，降低了零势阱附近的演化速度.本文在双势阱的基础上提出了一种新的势阱函数—单勾势阱函数.势阱函数$p(s)$以及扩散速率函数${d_p}(s)$表达如下：

(11)$p(s) = \left\{ {\begin{array}{*{20}{l}} { - \frac{1}{4}{s^4}{\rm{ + }}\frac{1}{8}{s^2}{\rm{ + }}\frac{1}{{25}}{\rm{, }}s \leqslant 0.756} \\ {\frac{1}{2}{{(s - 1)}^2}{\rm{, }}\quad\quad s > 0.756} \end{array}} \right.$

(12)${d_p}(s) = \frac{{p'(s)}}{s} = \left\{ {\begin{array}{*{20}{l}} { - {s^2}{\rm{ + }}\frac{1}{4}{\rm{, }}\quad s \leqslant 0.756} \\ {1 - \frac{1}{s}{\rm{, }}\quad\quad s > 0.756} \end{array}} \right.$

如图 2(a)所示，本文的势阱函数同属于双势阱函数，由零势阱和壹势阱构成.扩散速率函数为水平集梯度$\nabla \phi $的势阱函数的导数除以$\nabla \phi $，它能表示$\nabla \phi $的变化方向和速度，为方便表示，用s代替水平集梯度$\nabla \phi $.如图 2(b)所示，当$0 \leqslant s \leqslant 0.5$时，势阱函数的演化速度大于零，水平集向前演化，s向0靠近，保持水平集为常数；当$0.5 < s < 1$，势阱函数的演化速度小于零，水平集向后演化，s增大，向1靠近；$s \geqslant 1$时，势阱函数的演化速度大于零，水平集向前演化，s减小，同样使s向1靠近.这使得水平集函数保持符号距离特性.

与Li等提出的DRLSE模型中的双势阱函数相比，本文提出的单勾势阱函数有如下优势：（1）扩散速率函数在$s = 0$处时，不存在分数计算，因此不需要对分母取0时的函数进行修正，避免了算法的下降；（2）如图 3所示，扩散速率函数在$s = 0$处时的值较小，即势阱函数在零势阱附近的演化速度偏小，不会导致s急速减小，避免了曲线前进速度太快导致水平集侵入被分割目标内部，因此可以增强对弱边缘的检测准确性.

在以往的水平集模型中，为了抑制由噪声引起的非均匀性，常常使用大尺度高斯滤波对图像进行预处理^[10].而图像中的纹理变化常常是非均匀性的，恒定方差的高斯滤波不能平滑图像区域的非均匀性.若方差较大，对应的高斯滤波器模板尺度也较大，图像边缘损失较多，容易导致目标区域出现过分割；反之，若方差较小，高斯滤波器模板尺度较小，图像非均匀性程度过大，容易导致目标区域出现欠分割^[10].另外，DRLSE模型是基于边缘的分割方法，此类方法的主要缺点有：（1）初始轮廓会影响分割的精确度；（2）对图像噪声很敏感；（3）在弱边缘情况下演化结果不好^[11].而边缘指示函数g依赖光滑图像的梯度，所以在分割含噪图像时，演化曲线受噪声干扰容易陷入虚假边缘^[12].因此，当目标之间的距离很近时，此时若使用大高斯核滤波，会使边界处的灰度值被其邻域的加权灰度值取代，而其邻域中像素多为目标像素，故原先低灰度边界灰度上升，到一定程度后，边界消失，目标边界可能相互连接，导致分割失败^[13].

为了得到较好的分割效果，将序贯滤波应用于本文的水平集模型，利用小方差高斯滤波的优势，对图像进行多次小方差滤波的叠加，平滑过程可抽象为：

(13)$I = {I_0} * g({\sigma _1}) * g({\sigma _2}) = {I_0} * [g({\sigma _1}) * g({\sigma _2})]$

其中，${I_0}$为输入图像，I为高斯滤波后的图像，$g({\sigma _1})$、$g({\sigma _2})$是方差分别为${\sigma _1}$、${\sigma _2}$的高斯核，根据高斯滤波函数的数学表达式，可将该滤波过程表达为：

(14)$g({\sigma _3}){\rm{ = }}g({\sigma _1}) * g({\sigma _2}){\rm{ = }}{\int_{ - \infty }^{ + \infty } {\rm{e}} ^{ - [\frac{{{x^2}}}{{2\sigma _1^2}} + \frac{{{{(t - x)}^2}}}{{2\sigma _2^2}}]}}{\rm{d}}x = \frac{{\sqrt {2\pi } {\sigma _1}{\sigma _2}}}{{\sqrt {\sigma _1^2 + \sigma _2^2} }}\exp \left[ { - \frac{{{t^2}}}{{2(\sigma _1^2 + \sigma _2^2)}}} \right]$

(14) 式表明，两次高斯函数方差${\sigma _1}$、${\sigma _2}$的序贯滤波可等效为方差为${\sigma _3}$的一次高斯滤波.则可理解为，两次较小方差的高斯滤波叠加产生的平滑效果可等效为一次较大方差的高斯滤波.基于此，将叠加较小方差高斯滤波的次数与水平集迭代次数相关联，在设定好小方差的初始值${\sigma _0}$后，可由前面的公式推得第k次序贯滤波的方差${\sigma _k}$为：

(15)${\sigma _k} = {\sigma _0}\sqrt k $

为避免水平集迭代次数过多导致目标区域过度平滑失去边缘特征而出现曲线侵入，文献[10]提出用区域置信度作为判定的标准之一，根据不同等效方差的高斯函数分割结果的区域置信度决定序贯滤波次数，其区域置信度Pr定义为：

(16)$Pr = \frac{{S({A^{k - 1}} \cap {A^k})}}{{\max \{ S({A^{k - 1}}), S({A^k})\} }}$

其中，${A^{k - 1}}$为第$(k - 1)$次小方差高斯滤波后，即第$(k - 1)$次水平集演化的目标分割区域，$S({A^{k - 1}})$为其对应区域的面积.同理，${A^k}$为第k次迭代得出的目标区域，$S({A^k})$为其对应的面积.(16)式的分子为两次分割区域的交集，在算法的实现上，判定区域的交集会稍微复杂一些.本文结合DRLSE算法的优势，在对牙齿进行分割时，常常将初始轮廓设定在目标区域外部，故水平集多保持向内收敛，故上式中的$S({A^{k - 1}} \cap {A^k})$即为$\min \{ S({A^{k - 1}}), S({A^k})\} $，意为$S({A^k})$与$S({A^{k - 1}})$中的最小值.这样，不需要在引入判定交集增加算法的运行时间，只需简单的选择结构语句就可实现区域置信度的判定.故可得区域置信度为：

(17)$Pr = \frac{{\min \{ S({A^{k - 1}}), S({A^k})\} }}{{\max \{ S({A^{k - 1}}), S({A^k})\} }}$

在以往的研究中，对于同一患者的不同牙齿多采用相同尺寸的牙盒作为初始轮廓，这样虽然在选取时较为方便，但却是以增加水平集迭代次数、增加数据提取时间作为代价的.本文将针对不同牙齿的特点，设置不同尺寸、不同迭代次数，以提高分割效率.为了进一步利用相邻序列图片中同一结构的相似性，本文将先验信息^[14]加入对下层图像的分割中，以当前层对单一牙齿的分割结果作为下层图像中同一牙齿的初始轮廓，以求在相对少的迭代次数下得到准确的结果.考虑到牙齿的部分图像会出现当前层轮廓小于下层轮廓的现象，故在定义下层初始轮廓时先对当前层最终轮廓做形态学膨胀处理，以避免出现曲线侵入而使得结果不准确.

采用阈值法^[15]对牙槽骨分割时，由于牙髓区域位于牙齿内部，并且灰度值较低，通常被视为背景区域，导致分割结果产生孔洞.Wang等^[16]通过孔洞填充校正阈值分割得到的掩模，再利用自适应扩散流（Active Diffusion Flow，ADF）主动轮廓模型对单颗牙齿进行分割.本文提出一种基于边缘信息检测的填充算法，如图 4所示，其主要思想为：对阈值分割后已被识别为背景的像素点进行十字形扩散检测，若以该像素点为中心向外扩散的同半径十字上均存在非背景点，则认为此处为孔洞即中心像素点为非背景点，依次检测非背景点后即可将牙槽骨分割过程中出现的孔洞全部填充，并避免出现边界泄露的情况.

牙齿牙槽骨分割算法步骤如图 5所示，具体步骤如下：

Step1. 参数初始化：设定步长、正则项系数$\mu $、长度项系数$\lambda $、面积项系数$\alpha $、内外迭代次数

Step2. 读取DICOM图像，对其进行灰度均一化处理和对比度增强

Step3. 牙齿分割

a. 水平集初始化，对第一层图像手动取框作为初始水平集曲线

b. 在区域置信度内进行小尺度高斯滤波，并计算图像梯度，生成边缘停止函数

c. 计算单勾势阱函数演化速度

d. 计算水平集能量项最小化迭代公式，使其驱使轮廓线停止在牙齿边缘处得到分割结果

e. 读取下一层图像，对上一层的分割结果小尺度膨胀后作为下一层的初始水平集曲线

f. 重复步骤b~e，直到分割完全部图层，完成此颗牙齿的序列分割

g. 重复步骤a~f，进行下一颗牙齿的分割，直到完成所有牙齿的序列分割

Step4. 牙槽骨分割

a. 通过阈值法对牙槽骨进行分割，并对孔洞进行填充

b. 读取下一层图像

c. 重复步骤a~c，直到分割完全部图层，完成牙槽骨分割

本文实验数据为10组CBCT口腔图像以及3组口腔磁共振图像，采集于上海市第九人民医院，其中CBCT扫描电压为90 kV，曝光时间为10.8 s，获得数据矩阵大小为595×595×358，像素分辨率为0.25×0.25×0.25 mm³；MRI扫描序列为T1W_TSE，TR为507 ms，层厚为4 mm，图片尺寸为336×336.图像的数据类型为16位无符号整型，保存的格式为DICOM文件.

本文的实验环境为：CPU处理器为Inter(R) Core(TM) i5 CPU(2.30GHz)，内存8.0GB，MATLAB R2016b软件下编写程序完成实验，实验参数为$\mu = 0.04$、$\lambda = 5$、$\alpha = 1.5$、内迭代8次、外迭代16次.

将16位数据转换为8位的灰度图像，手动画取合适大小的相同初始框，并选择相同的实验参数后，分别采用DRLSE模型和本文提出的水平集模型对CBCT图像及磁共振图像中的牙齿进行分割.图 6为CBCT图像分割结果，可以发现Li等的方法可能会把相邻牙齿的轮廓错分为目标牙齿的轮廓，导致欠分割.而本文提出的算法模型在相邻牙齿相互粘连、牙齿和牙槽骨灰度值相近的情况下，能够驱使水平集函数表示的曲线到达目标牙齿的边界处，而不会停留在相邻牙齿和牙槽骨的分界处.由此可见，本文的算法对牙齿有精确的分割能力，可以很好地改善DRLSE模型的欠分割的情况.图 7为文本算法运用于磁共振图像的实验结果，亦可准确有效地分割出磁共振图像中的单颗牙齿.

本文从分割的精确性及运算时间两个角度来评价算法的优劣.采用积重叠误差（Volumetric Overlap Error，VOE）^[17]来进行基于区域的精确性描述，VOE定义如下：

(18)$VOE = 1 - \frac{{\left| {A \cap B} \right|}}{{\left| {A \cup B} \right|}}$

其中A为算法自动分割的结果，B为手动勾画结果并作为标准.通过计算水平集分割算法和手动分割结果的交集和并集的体积比得到的就是两者的真实重叠度，VOE越小，说明分割的结果与标准的重叠度越高，即分割结果越准确.基于CBCT图像，本文算法与DRLSE分割结果的定量比较如表 1所示，可见本文算法在分割精确度上相比DRLSE有明显的提升，算法效率上也有一些提高，但由于序贯滤波增加了算法的时间，因此该优势并不明显.

利用本文提出的水平集算法对同一样本中同一颗牙齿的CBCT多层序列图像进行分割，结果如图 8所示.对第一层图像进行手动画框后，即可通过本文的水平集算法对目标牙齿进行准确的分割，将此层分割结果曲线进行小尺度膨胀后作为下一层的水平集初始轮廓，即可无中断地进行序列图像的牙齿分割，避免了重复初始化.实验结果表明，利用本文的方法能过对CBCT序列图像完成准确的牙齿分割，基本不会出现欠分割或是过分割的情况.

采用本文提出的水平集算法可对一层图像中的所有牙齿依次完成分割，得到的牙齿轮廓即为牙槽骨的内轮廓.再采用阈值分割法可得到牙槽骨的外轮廓，由于牙槽骨和牙齿内部存在灰度值较小的区域，阈值法得到的结果会存在孔洞，通过采用前文所述的孔洞填充的方法，选择尺度为15个像素点对其进行填充，即可分割出完整的牙槽骨，CBCT图像牙槽骨分割结果如图 9所示，红色曲线为各颗牙齿的轮廓，绿色曲线即为牙槽骨的外轮廓，红、绿曲线之间的区域即为分割出的牙槽骨.对下牙槽骨进行序列分割后，将结果进行三维重建，结果如图 10所示.

本文针对CBCT及磁共振图像中牙齿的特点，对以往水平集模型中不适合应用于此类图像的部分做出了改进，并将其应用在牙齿CBCT及磁共振图像分割中，得到了更为准确的分割结果.其中创新的单钩势阱函数克服了现有势阱函数使用限制多、弱边缘难以检测等缺点；将序贯滤波与水平集模型结合改善了原有水平集模型中单次使用大尺度高斯核而使图像信息丢失的情况；根据断层序列图像特点，将当前层分割结果作为下层图像初始框，实现先验信息的充分利用.