杨 玲 李 宇 王中科 饶妮妮

1（电子科技大学生命科学与技术学院，成都 6 10054）2（成都信息工程学院电子工程学院，成都 6 10225）3（成都信息工程学院网络工程学院，成都 6 10225）

引言

牙齿的三维建模在牙种植领域中有着十分重要的意义。目前，基于特征点的三维重建算法由于其数据量少而精和可快速重建目标物的结构特征，具有较强的实用性［1］。目标物的特征点从几何方面可分为点、线、面三种，其中点特征是最常采用的一种，如灰度的局部最大值点、边缘点、角点和线交叉点等［2－8］。本研究中提及的牙齿特征点拟指在牙齿CT各断层图像中表示牙齿形状轮廓的边缘特征点。曲线或曲面演化是解决图像分割和目标检测的一种有效方法。曲线演化是曲线在自身内力与外力（代表图像信息）的共同作用下，变化收敛到图像的边缘或其它特征处的过程。因此，曲面演化常在医学上用于特征点的提取。曲线演化的方式有多种，如Sethian等提出的Level Set算法［9］，Kass等提出的Snake方法［10］，Chan等提出的基于Mumford－Shah理论的CV模型［11］，以及近期提出的窄带法和无重复初始化速度函数等方法［12－13］。上述方法为处理不同拓扑结构物体的变形提供了方便，且抗噪性能较强，然而这些曲线演化方法仍然存在许多的不足之处。例如，基于水平集的曲线演化依赖于图像梯度的局部信息，在梯度不明显处或对比度不大的区域容易产生曲线泄漏［14］；而CV模型没有利用图像局部的梯度信息，当初始水平集位于平滑区域或者极度凹陷区域进行分割时，水平集曲线的收敛速度会非常缓慢［15］。尽管在近期提出的算法中，对速度函数、参加运算的数据量进行了改善，但运算效率仍是该类算法需要重点优化的问题。

本研究结合曲线演化算法抗噪能力强的优势和牙齿CT序列图像数据的特点，改进离散曲线演化（discrete curve evolution，DCE）模型，并基于改进的DCE提出一种提取牙齿边缘特征点的新方法。新方法通过控制不同层图像的特征点数，能快速、准确地提取牙齿的三维特征点，既保证了特征点提取的信息量足以重建出整颗牙齿，又能降低数据的冗余量，简化了演化过程。最后，采用真实的磨牙CT图像数据验证新方法的可行性和计算效率。

1 离散曲线演化算法原理

DCE算法是Latecki等为提取目标物骨架而提出的一种快速获得复杂多边形拓扑结构的算法［16］。由于图像中物体的轮廓常受噪声影响而形变，DCE实际上是轮廓预处理过程，包括两方面的含义：一是除去边界噪声；二是保留视觉感官上的重要部分。它对复杂多边形进行离散曲线演化，既能降低噪声的影响，又能排除不重要的形状特征，从而简化形状。文献［17，18］证明了这个结论。

DCE算法在获取边缘信息时，首先采用形态学算法对图像进行去噪、二值化，提取出感兴趣的连通域，然后采用梯度算子和轮廓跟踪法获取连通域的封闭轮廓线，以此作为DCE演化的初始特征点，然后从最大凸（或凹）弧线的起点开始演化，获取图像边缘上的突变点，并在每一步演化中去除相关性最小的顶点，最后保留该图形的拓扑结构。其中初始轮廓点的选取对曲线演化结果有着重要的影响，算法中采用了多种去噪方法去除图像的伪边缘，降低了初始轮廓点对噪声的敏感性，以便获取单一的闭合轮廓曲线，确保曲线演化结果的可靠性。

由于任意数字曲线可以看作一个含有大量顶点的多边形，所以能将多边形的演化用数字曲线来代替。设C为实平面上的一闭合多边形，S1…，Sm为构成C的连续线段，P1…Pm，为C的顶点。其中，为Si和Si+1的公共顶点。为得到闭合曲线，可约定Pm+i=Pi。在曲线演化的每一步中，找出相关度最小的顶点Pi，用一段连接Pi+1和Pi－1的线段代替原来的两段相邻线段Si和Si+1。这样每一步将去掉多边形的一个顶点Pi。

由此可见，离散曲线演化的基本思想是：在演化的每一步，对连续线段S1、S2用单一线段S1∪S2连接端点替代，如图1所示。图中是在曲线演化过程中将相关性最小的Pi点去掉。

图1 DCE演化示意图Fig.1 The evolution scheme of DCE

每个顶点的相关度参数K的计算式为

式中，β（S1，S2）是线段 S1和 S2的拐角，以逆时针为正；根据曲线的凹凸特性，β（S1，S2）的计算方法有所差异。

图2所示为曲线凹凸对拐角的影响，其中b点为当前顶点，a、c分别为b点前后→相邻的两个顶点，S1表示b和a两点构成的线段—ba、S2表示 b 和 c 两点构成的线段→—bc，β表示从线段 S 逆时针旋转至线1段S2的夹角。图（a）是凸曲线的拐角，小于90°，则β（S1，S2）的计算公式为

式中，l为线段的长度函数，计算公式为

式中，（ΔxS1，ΔyS2），（ΔxS2，ΔyS1）分别表示线段 S1和S2两端点x和y的坐标差值。

图（b）是凹曲线的拐角，大于 9 0°，则 β（S1，S2）的计算公式为

根据式（1）计算的K值越高，则表明弧S1∪S2对形状的贡献越大。DCE算法将K作为曲线演化的控制参数。如对采样后的闭合多边形C，若P1，…，Pm为C的顶点，可定义M（p）为其顶点个数，也即特征点个数，设定其阈值为T，使DCE算法演化时，通过K值的计算，依次去掉K较小的点，使M（p）=T，从而结束曲线的演化，产生相应图像的骨架。

2 DCE算法改进

DCE算法可以保留图形的拓扑结构。可知，采用DCE提取的图像特征点个数越多，表达的图像信息就越准确。然而，由于DCE算法停止演化是由设定的控制点数来确定，因此，对序列图像而言，这会造成不同结构具有相同的停止条件，也会造成数据点冗余或不足，最终导致三维重建达不到满意的效果。

图3（a）为一磨牙的第330层CT图像（512像素×512像素）。用DCE算法对其右上方连通域进行特征提取，结果如图3中（b）～（e）所示，其中（b）和（d）将拐角定义为凸点提取，分别提取了10和20个特征点；（c）和（e）将拐角定义为凹点提取，分别提取了10和20个特征点。图3结果表明，DCE算法会因拐角设置不同而导致提取的特征点存在差异，同时也表明特征点个数增加到某种程度时，拐角定义为凹或凸对特征点的提取影响不大。

为了克服DCE算法的特征点表达信息不够的问题，现有两种解决方案：一是定义两种拐角，对一幅图像先进行凹(凸)点提取，再进行凸(凹)点提取；二是将每幅图像特征点的数目提高到一定程度。然而，上述两种方案会使程序的执行效率降低，甚至产生特征点的冗余问题。

本研究提出改变DCE的停止准则M(p)=T，使之成为可控点的操作方式；并提出曲线特征因子量，使DCE特征点数量的选取具有自适应性和合理性。

设闭合曲线Y=F(X)在某一区间I上连续，曲线上共有N个拐点，第j个拐点为(x(j)，F(x(j)))(j=1，2，…，N)，则将 N相对于构成该曲线所有点的个数即曲线长度进行归一化处理，用来表示曲线的复杂程度，并称之为曲线特征因子量(characteristic factor of curve，CFOC)。即

式中，L为曲线上数据点总数。

如图4所示是一牙齿CT图像(分辨率为512像素×512像素)的13至483断层的曲线特征因子量。横坐标表示不同断层图像，纵坐标表示特征因子量。图4中(a)是13至483断层的曲线特征因子量，每一层图像对应的多个离散点分别表示多个连通域的特征因子量。取其各层曲线特征因子量中最大值对应的点，得到13至483断层的最大曲线特征因子量，如图4(b)所示。为自适应地确定不同断层图像的特征点数Z，可将特征因子量CFOC分成若干等级，不同等级提取不同数目的特征点。

在保留其拓扑前提下，改进算法对不同等级的特征因子量设定特征点数Z的原则是：若各断层曲线相差较大，即CFOC曲线随CT层数变化陡峭，则各层的Z值可设定为跨度较大的值；反之，当各断层曲线相差较小，即CFOC曲线随CT层数变化比较平缓，各层的Z值可取跨度较小的值。

分别将CFOC分成了5个和8个等级，并按式(7)设定了各断层的特征点数Z与最大曲线特征因子量等级n的关系：

从图 4(b)可以看出，当取［0.0、0.2］，［0.2、0.4］，［0.4、0.6］，［0.6、0.8］，［0.8、1.0］5 个等级区间时，曲线变化相对平坦的区域归为同一区间，而不同区间变换相对陡峭，此时不同区间可采用跨距较大的特征点数 Z，即为 3、5、9、17和 33个。而当分成8个均匀等级区间时，就会降低各区间之间的陡峭性，使Z值跨距减少，冗余点增多。

图4 牙齿CT序列切片的断层曲线特征因子量。(a)各断层曲线特征因子量；(b)各断层最大曲线特征因子量Fig．4 CFOC of Dental CT slices．(a)CFOC of each slice；(b)The max CFOC of each slice

该方案既能保证牙齿图像特征点提取的有效性，又能减少数据的冗余度。

3 实验与结果分析

图5为一磨牙的第29、119、229和329层 CT图像。采用改进的DCE算法对图5中4层CT图像进行特征点提取，结果如图6所示。图6中(a)～(c)分别是从图5的第29层CT图像的3个连通域中提取的特征点，且从每个连通域中提取的特征点数均为5；图6(d)是从图5的第119层CT图中提取的特征点，数目为17个；图6(e)是从图5的第229层CT图中提取的特征点，数目为17个；图6中(f)～(h)分别是从图5的第329层图像的3个连通域中提取的特征点，每个连通域提取了9个特征点。从图6的结果可以看出，改进的DCE算法能够从磨牙CT中准确有效的提取出特征点。

接着可利用提取的特征点，重建牙齿的三维曲面，具体方法是：对每一层牙齿的特征点采用不同比率的样条插值得到各断层曲线数据，再将每层曲线数据扩展到三维空间，用样条曲线插值生成对应的曲面数据。图7是在VTK（Visualization ToolKit）开发平台下，利用图6的特征点，采用样条曲线重建出的牙齿表面图。可以看出，牙齿整个表面重建完整，牙冠和牙根明显，与原牙齿曲面结构吻合，证实了所提出算法的可行性和有效性。

图5 一磨牙的4层CT图像（512像素×512像素）。（a）第29层；（b）第119层；（c）第229层；（d）第329层Fig.5 Four CT slices of a molar（512 pixel×512 pixel）.（a）The 29thslice；（b）The 119thslice；（c）The 229thslice；（d）The 329thslice

4 算法效率分析与比较

表1给出了采用DCE和改进的DCE算法提取图5磨牙CT图像特征点的计算效率情况。DCE算法没有经过特征点个数判定，直接在每一层牙齿CT图的各个连通域提取固定的12个特征点，所需总时间为7 min 35 s，每一层平均所需时间为0.97 s。改进的DCE算法能够自适应地确定每一层牙齿CT图的各个连通域需要提取的特征点数目，降低了数据存储的冗余量，同时提高对目标特征表达的精确性。其中特征点数目确定方法所需的总时间为2 min 56 s和3 min 11 s，每一层平均所需时间分别为0.37 s和0.41 s。可见改进算法在每层上提取特征点所花时间约为原算法的50％。此外从提取的总特征量而言，改进算法提取的特征量约为原算法的80％左右，因此改进DCE算法的计算效率明显优于DCE算法。

图6 对图5中的4层CT图像进行特征点提取的结果。（a）对第29层图像的连通域1提取5点；（b）对第29层图像的连通域2提取5点；（c）对第29层图像的连通域3提取5点；（d）对第119层图像的连通域提取17点；（e）对第229层图像的连通域提取17点；（f）对第329层图像的连通域1提取9点；（g）对第329层图像的连通域2提取9点；（h）对第329层图像的连通域3提取9点Fig.6 The extracted feature points of Fig.5’s four CTslices.（a）5pointsofthefirst connected area of the 29thslice；（b）5 points of the second connected area of the 29thslice；（c）5 points of the third connected area of the 29th slice；（d）17 points of the connected area of the 119thslice；（e）17 points of the connected area of the 229thslice；（f）9 points of the first connected area of the 329thslice；（g）9 points of the second connected area of the 329thslice；（h）9 points of the third connected area of the 329th slice

图7 利用改进DCE算法提取的特征点重建的牙齿表面图Fig.7 The reconstructed tooth surface with the extracted feature points by improved DCE algorithm

表1 DCE和改进DCE算法提取牙齿13～483层CT图特征点效率的对比Tab.1 The contrast of efficiency off eature points extracted by DCE and improved DCE for 13rd～483rd CT slices

从表1中也可看出，当选择不同的特征因子量等级划分方法和特征点数，改进算法效率会有所差异，故为发挥改进算法的最优效率，还需进一步优化特征因子量的等级划分和特征点数的选择方法。

5 结论

本研究利用改进后的离散曲线演化算法提取牙齿CT图像的边缘特征点，并在VTK平台上用样条插值方法对牙齿特征点进行曲面重建，验证了算法的有效性。所提出的特征点提取算法的优势主要体现在：一是由于采用基于曲线演化模型，所以算法的抗噪能力强，特征点的定位准确；二是算法在提取特征点之前自适应地计算所要提取的特征点数量。由于不同信息量的连通域所需的特征点数不同，因此，本算法既保证了获取的特征点信息能足以表达出整幅牙齿信息，也降低了特征点的冗余量。三是本算法相比传统曲线演化模型，效率较高，对多幅图像的处理具有优势。实验结果表明，该算法能够准确地提取出牙齿图像边界的特征点，且重建出的牙齿三维曲面结构完整清晰，其性能优于相关传统算法。

所提出算法还存在计算速度较慢的问题，今后的工作将致力于算法执行效率的提高。初步思路是进一步优化各断层特征点数与特征因子量等级的函数关系，并对同一层包含多个连通域的情况，可分别根据每个连通域的自身的曲线特征因子量选取不同的特征点数，以克服目前的统一处理模式，从而进一步降低特征点的冗余度，提高算法的执行效率。

［1］王付新，黄毓瑜，孟偲.三维重建中特征点提取算法的研究与实现［J］.工程图学学报，2007，3：91－96.

［2］谢萍，马小勇，张宪民.一种快速的复杂多边形匹配算法［J］.计算机工程，2003，29（16）：176－181.

［3］彭文，童若锋，钱归平.使用特征点与曲线配准医学图像［J］.计算机辅助设计与图形学学报，2007，19（9）：1126－1131.

［4］Baboulaz L，Dragotti PL.Exact feature extraction using finite rate of innovation principles with an application to image superresolution［J］.IEEE Transactions on Image Processing，2009，18（2）：281－298.

［5］邵巍，朱圣英，陈灵芝.一种快速多尺度特征点匹配算法［J］.中国图像图形学报，2009，14（12）：2572－2576.

［6］雷琳，王壮，粟毅.基于多尺度Gabor滤波器组的不变特征点提取新方法［J］.电子学报，2009，37（10）：2314－2319.

［7］李竹林，王文东，赵宗涛.改进的边缘特征点提取算法［J］.计算机工程与应用，2009，45（2）：63－65.

［8］王建文，李青.基于点和边缘相结合特征提取的图像配准算法［J］.计算机工程与设计，2009，30（4）：928－930.

［9］Osher S，Sethian J.Fronts propagation with curvature dependent speed：Algorithms based on Hamilton-Jacobi formulations［J］.Journal of Comp.Phys.1988，79：12－49.

［10］Kass M，Witkin A，Terzopoulos D.Snakes：Active contour models［J］.International Journal of Computer Vision，1988，1（4）：321－331.

［11］Chan TF，Vese LA.Active Contours Without Edges［J］.IEEE Transactions on Image Processing，2001，10（2）：266－276.

［12］Wang Hongyuan，Zhou Zeming.Study on the Narrow Band Level Set Method［J］.Journal of Communication and Computer，2006，3（9）：37－40.

［13］Li Chunming，Xu Chenyang，Gui Changfeng，et al.Level Set Evolution Without Re-initialization：A New Variational Formulation［C］//Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition.San Diego：IEEE，2005，1：430－436.

［14］温铁祥，杨丰.一种新的曲线演化混合模型图像分割算法［J］.电路与系统学报，2007，12（4）：48－52.

［15］Vese LA，Chan TF.A multiphase level set framework for image segmentation using the Mumford and Shah model［J］.International Journal of Computer Vision，2002，50（3）：271－293.

［16］Latecki LJ，Lakamper R.Shape similarity measure based on correspondence of visual parts［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2000，22（10）：1185－1190.

［17］Bai Xiang，Latecki LJ，Liu Wenyu.Skeleton Pruning by Contour Partitioning with Discrete Curve Evolution［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2007，29（3）：449－462.

［18］Krinidis S，Chatzis V.A skeleton family generator via physicsbased deformable models［J］.IEEE Transactions on Image Processing，2009，18：1－11.