CN107564023B - 一种cbct牙齿分割及建模算法 - Google Patents

Abstract

本发明公开了一种CBCT牙齿分割及建模算法。采集CBCT二维图像序列，对图像序列中各层的图像进行平滑处理；构建所需要分割牙齿的初始轮廓，构建水平集函数；利用底层特征对水平集函数进行迭代计算直至收敛，获得分割轮廓结果；重复步骤对每层图像进行处理，第一次以牙齿轮廓清晰的一层作为初始层，并将当前层图像的分割轮廓结果进行扩张作为下一相邻层图像的初始轮廓，使得层之间传递迭代，获得所有层图像分割轮廓结果；由所有层的分割轮廓结果组成三维图像，重建三维模型。本发明能有效地从CBCT图像中抽取并重建整颗牙齿的三角网格模型，可辅助牙齿矫正过程的评估过程。

Description

一种CBCT牙齿分割及建模算法

技术领域

本发明涉及了数字口腔领域，具体地说是涉及了一种CBCT牙齿分割及建模算法。

背景技术

近年来，数字口腔技术发展迅速，牙齿隐形正畸也逐渐普及。在对牙齿正畸过程中，需要建立牙齿的三维网格模型，用于虚拟矫治。牙根的信息只能通过CBCT扫描得到，故提出一种CBCT牙齿分割及建模算法，用于从CBCT数据中抽取并重建牙齿的三维网格模型。

发明内容

为了解决背景技术中存在的问题，本发明提供了一种CBCT牙齿分割及建模算法。

本发明所采用的技术方案如下：

1)采集CBCT二维图像序列，对图像序列中各层的图像进行平滑处理；

2)构建所需要分割牙齿的初始轮廓，构建水平集函数；

3)利用底层特征对水平集函数进行迭代计算直至收敛，获得分割轮廓结果；

所述的底层特征包括梯度、灰度值和水平集函数。

4)重复上述步骤2)～3)对每层图像进行处理，第一次处理时选取CBCT二维图像序列中牙齿轮廓清晰的一层作为初始层，每次重复步骤2)～3)处理时将当前层图像的分割轮廓结果进行扩张作为下一相邻层图像的初始轮廓，使得层之间传递迭代，获得所有层图像针对所需要分割牙齿的分割轮廓结果；

具体实施中采用牙冠中点所在层作为初始层。

下一层是指相邻的上层或者下层，具体实施中采用牙冠中点所在层作为初始层后，向上和向下的相邻层图像分别进行传递迭代。

5)由所有层的分割轮廓结果组成三维图像，重建三维模型。

所述步骤B具体是包括：

2-1)建立在牙齿轮廓附近的凸多边形作为初始轮廓，也作为感兴趣区域(Regionof Interests，简称ROI)，并使得凸多边形仅包围该颗牙齿，而不包围其他牙齿；

2-2)构造一个以图像长宽为维度的矩阵，矩阵的大小和图像大小相同；

2-2)对矩阵赋值，使得初始轮廓外的元素为正值，初始轮廓内的元素为负值，以矩阵作为水平集函数Phi。此矩阵值的改变即代表了初始轮廓形状的改变。

所述步骤3)具体是包括：

3-1)根据图像分别沿X和Y方向的梯度值，X和Y方向是指图像的水平和竖直方向，采用以下公式计算梯度算子g：

其中，I表示灰度值矩阵，灰度值矩阵是由CBCT二维图像的所有像素值经归一化后组成的矩阵，G_X和G_Y表示对灰度值矩阵在X和Y方向上求导的结果，G_X ²表示矩阵G_X中的每个元素分别取平方；

3-2)利用梯度算子g采用以下公式计算图像的区域算子Area：

其中，矩阵的

符号是矩阵的逐元素乘法，表示两个矩阵同样位置的元素相乘，Phi表示图像的水平集函数，Phi_i,j表示水平集函数中第i列第j行的元素，i表示图像X方向像素点的序号，j表示图像Y方向像素点的序号，Del表示水平集函数的筛选函数，Del_i,j表示筛选函数Del中第i列第j行的元素，ε为预先设定的水平集函数阈值；

3-3)利用梯度算子g采用以下公式计算图像的边界算子Edge：

其中，

运算表示对矩阵的每个元素值计算开方，small为边界参数矩阵，其行列数与图像的行列数一致，所有元素取同一个值，具体实施中取一个很小的正数值如0.0001；Phi_X表示水平集函数在X方向的偏导，Phi_Y表示水平集函数在Y方向的偏导，v_x表示梯度算子在X方向的偏导，v_Y表示梯度算子在Y方向的偏导，N_x表示经过归一化的Phi_X函数，N_Y表示经过归一化的Phi_Y函数，curv表示函数N_x的曲率函数；

3-4)利用梯度算子g和边界算子Edge采用以下公式计算距离归一化算子Reg：

式中，lap表示水平集函数的拉普拉斯算子，S表示水平集函数的归一化矩阵，S'表示经过布尔筛选函数筛选后的水平集函数的归一化矩阵，S'_i,j表示S’矩阵的第i列第j行的元素，S_i,j表示S矩阵的第i列第j行的元素，A和B分别表示水平集函数的两个筛选矩阵，p表示双阱势函数，p’表示经过布尔筛选函数筛选后的p，p'_i,j表示p’矩阵的第i列第j行的元素，p_i,j表示p矩阵的第i列第j行的元素，de表示矩阵的布尔筛选函数，dp表示双阱势函数的归一化导函数；

3-5)如果当前层是第一层，则以所有项均为0的零矩阵作为形态学算子，直接进行步骤3-6)；

如果当前层不是第一层，则根据上一层的分割结果和当前的水平集函数计算形态学算子Shape：

Shape＝Phi-Phi₀

其中，Phi₀是上一层图像的分割结果，Phi是当前层图像的分割结果；

3-6)利用归一化算子、区域算子、边界算子和形态学算子，采用以下公式计算水平集函数进行迭代计算：

其中，Phi’表示迭代计算前的水平集函数，step表示水平集函数演化速度参数，α、β和γ分别表示区域算子的权重、边界算子的权重和形态学算子的权重；

迭代计算分为先后两个阶段：

第一阶段，取α、β和γ的值为1.5、5和0.25，迭代40次，使得水平集函数所表示的形态持续改变至收敛；

第一阶段，取α和β和γ的值为0、5和0.25，迭代10次，使得水平集函数所表示的形态不在改变，水平集函数的数值趋于稳定，完成水平集函数的演化；

3-7)然后对步骤3-6)最终获得的水平集函数进行处理，将其中正数的元素值置位0，将负数的元素值置位255，以255所包含的元素值对应到图像中的所有像素点作为分割轮廓内部区域，以0所包含的元素值对应到图像中的所有像素点作为分割轮廓外部区域，分割轮廓内部区域和分割轮廓外部区域之间的闭合路径作为分割轮廓，获得分割轮廓结果的二维位图图像。

所述步骤4)中，将当前层图像的分割轮廓结果进行扩张作为下一相邻层图像的初始轮廓具体为：将当前层的分割轮廓向外扩张形成新轮廓，扩张是指分割轮廓向外均匀偏置若干像素，新轮廓内所对应的水平集函数中的各个元素值均置为负值，新的轮廓外所对应的水平集函数中的各个元素值均置为正值。

具体实施中可将新轮廓内元素值置为-2，新轮廓外元素值置位2。

所述步骤5)具体是包括：

5-1)将分割轮廓结果的二维位图图像序列按照顺序依次叠加形成三维图像；

5-2)使用Marching Cubes算法(Marching cubes:A high resolution 3Dsurface construction algorithm)，以三维图像作为输入，以255作为Marching Cubes算法中的阈值参数，计算得到三维网格模型。

本发明的有益效果是：

本发明能够精确地从CBCT数据中抽取包括牙根在内的整颗牙齿的三角网格模型，有利于后期的矫正实现。由于本发明在图像之间的迭代和分割是自动的，故可以大大提升效率。

附图说明

图1为一张CBCT二维图像；

图2为用户选取牙齿范围示意图；

图3为经过平滑处理后的图像；

图4为图像的梯度算子图，颜色越浅(白)表示梯度算子值越大；

图5为初始状态的水平集函数图像；

图6三张小图分别为迭代次数为0,20和40次时的水平集函数图像；

图7为最终分割结果，白色轮廓即为分割得到的牙齿轮廓；

图8为当前层图像分割结果扩张的示意图；

图9为保存的牙冠至牙根其中六层的分割结果图像；

图10为重建后的三维网格模型。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细说明。

按照本发明方法的实施例及其具体实施过程如下：

CBCT图像中牙齿和牙槽骨的灰度值十分接近，同时不同层牙齿的灰度值差异较大，难以通过单一阈值进行分割，故设计此算法进行自动分割。图1为一张CBCT图像，该层图像上牙齿的轮廓较清晰，故选做初始层。首先用户通过选取四个点确定一颗牙齿的大致区域，如图2。为了避免噪点，对图像进行光滑处理，光滑处理后的图像如图3所示。

根据公式1计算得到图像的梯度算子如图4，其中一个像素点的颜色越深(黑)表示该像素点的梯度算子值越小，反正，颜色越浅(白)表示该像素点的梯度算子值越大。初始状态下的水平集函数由用户刚刚选定的范围确定，范围内水平集函数的值为-2，范围外水平集函数的值为2，将图像灰度值范围设置成-10到10，则此时水平集函数图像如图5所示。

对水平集函数进行迭代变换，其函数图像值如图6所示，左图为第一次迭代的结果，中图为第20次迭代的结果，右图为第40次迭代的结果，水平集函数越来越接近分割结果。根据水平集函数值的正负形，确定边界，将边界绘制于CBCT图像上，得到图7的结果，此时一层图像的分割完成。

将当前层的分割形状进行扩张操作，如图8所示，左图为当前层的分割形状，右图为扩张后的分割形状，根据右图形状构建水平集函数作为下一层图像的初始水平集函数，重复上一层图像的操作。每一层图像的分割结果均以二值图像进行保存，牙齿区域内的像素值置位255，区域外的像素值置位0，图9是牙齿不同部位的6张分割结果示意图，从左上到右下代表从牙冠到牙根方向。

将所有层CBCT图像的分割结果，即一个二维图像序列按序排列，组成一个三维图像，使用Marching Cubes算法进行重建，得到整颗牙齿的三维网格模型，绘制效果如图10所示。

由此本实施例实现了从CBCT图像抽取并重建整颗牙齿的三角网格模型过程，包含牙根信息，有利于后期的矫正实现。目前常用的手工方法需要医生等专业人士逐层手动分割牙齿，单颗牙耗时约10分钟，且需要逐层手工操作，而我们的方法可以在3到4分钟内完成分割，并且只需要进行一次交互，故可以大大提升效率。

Claims (4)

1.一种CBCT牙齿分割及建模方法，其特征在于：

1)采集CBCT二维图像序列，对图像序列中各层的图像进行平滑处理；

2)构建所需要分割牙齿的初始轮廓，构建水平集函数；

3)利用底层特征对水平集函数进行迭代计算直至收敛，获得分割轮廓结果；

4)重复上述步骤2)～3)对每层图像进行处理，第一次处理时选取CBCT二维图像序列中牙冠中点所在层作为初始层，每次重复步骤2)～3)处理时将当前层图像的分割轮廓结果进行扩张作为下一相邻层图像的初始轮廓，使得层之间传递迭代，获得所有层图像的分割轮廓结果；

5)由所有层的分割轮廓结果组成三维图像，重建三维模型；

所述步骤3)具体是包括：

3-1)根据图像分别沿X和Y方向的梯度值，X和Y方向是指图像的水平和竖直方向，采用以下公式计算梯度算子g：

其中，I表示灰度值矩阵，灰度值矩阵是由CBCT二维图像的所有像素值经归一化后组成的矩阵，G_X和G_Y表示对灰度值矩阵在X和Y方向上求导的结果；

3-2)利用梯度算子g采用以下公式计算图像的区域算子Area：

其中，矩阵的

符号是矩阵的逐元素乘法，表示两个矩阵同样位置的元素相乘，Phi表示图像的水平集函数，Phi_i,j表示水平集函数中第i列第j行的元素，i表示图像X方向像素点的序号，j表示图像Y方向像素点的序号，Del表示水平集函数的筛选函数，Del_i,j表示筛选函数Del中第i列第j行的元素，ε为预先设定的水平集函数阈值；

3-3)利用梯度算子g采用以下公式计算图像的边界算子Edge：

其中，

运算表示对矩阵的每个元素值计算开方，small为边界参数矩阵，Phi_X表示水平集函数在X方向的偏导，Phi_Y表示水平集函数在Y方向的偏导，v_x表示梯度算子在X方向的偏导，v_Y表示梯度算子在Y方向的偏导，N_x表示经过归一化的Phi_X函数，N_Y表示经过归一化的Phi_Y函数，curv表示函数N_x的曲率函数；

3-4)利用梯度算子g和边界算子Edge采用以下公式计算距离归一化算子Reg：

de(p)＝p',

de(S)＝S',

式中，lap表示水平集函数的拉普拉斯算子，S表示水平集函数的归一化矩阵，S'表示经过布尔筛选函数筛选后的水平集函数的归一化矩阵，S'_i,j表示S’矩阵的第i列第j行的元素，S_i,j表示S矩阵的第i列第j行的元素，A和B分别表示水平集函数的两个筛选矩阵，p表示双阱势函数，p’表示经过布尔筛选函数筛选后的p，p'_i,j表示p’矩阵的第i列第j行的元素，p_i,j表示p矩阵的第i列第j行的元素，de表示矩阵的布尔筛选函数，dp表示双阱势函数的归一化导函数；

3-5)如果当前层是第一层，则以所有项均为0的零矩阵作为形态学算子，直接进行步骤3-6)；

如果当前层不是第一层，则根据上一层的分割结果和当前的水平集函数计算形态学算子Shape：

Shape＝Phi-Phi₀

其中，Phi₀是上一层图像的分割结果，Phi是当前层图像的分割结果；

3-6)利用归一化算子、区域算子、边界算子和形态学算子，采用以下公式计算水平集函数进行迭代计算：

其中，Phi’表示迭代计算前的水平集函数，step表示水平集函数演化速度参数，α、β和γ分别表示区域算子的权重、边界算子的权重和形态学算子的权重；

迭代计算分为先后两个阶段：

第一阶段，取α、β和γ的值为1.5、5和0.25，迭代40次，使得水平集函数所表示的形态持续改变至收敛；

第二阶段，取α和β和γ的值为0、5和0.25，迭代10次，使得水平集函数所表示的形态不再改变，水平集函数的数值趋于稳定，完成水平集函数的演化；

3-7)然后对步骤3-6)最终获得的水平集函数进行处理，将其中正数的元素值置位0，将负数的元素值置位255，获得分割轮廓结果的二维位图图像。

2.根据权利要求1所述的一种CBCT牙齿分割及建模方法，其特征在于：

所述步骤2)具体是包括：

2-1)建立在牙齿轮廓附近的凸多边形作为初始轮廓，并使得凸多边形仅包围该颗牙齿，而不包围其他牙齿；

2-2)构造一个以图像长宽为维度的矩阵，矩阵的大小和图像大小相同；

2-2)对矩阵赋值，使得初始轮廓外的元素为正值，初始轮廓内的元素为负值，以矩阵作为水平集函数Phi。

3.根据权利要求1所述的一种CBCT牙齿分割及建模方法，其特征在于：

所述步骤4)中，将当前层图像的分割轮廓结果进行扩张作为下一相邻层图像的初始轮廓具体为：将当前层的分割轮廓向外扩张形成新轮廓，新轮廓内所对应的水平集函数中的各个元素值均置为负值，新的轮廓外所对应的水平集函数中的各个元素值均置为正值。

4.根据权利要求1所述的一种CBCT牙齿分割及建模方法，其特征在于：

所述步骤5)具体是包括：

5-1)将分割轮廓结果按照顺序依次叠加形成三维图像；

5-2)使用Marching Cubes算法，以三维图像作为输入，以255作为Marching Cubes算法中的阈值参数，计算得到三维网格模型。

Images