一种口腔软硬组织CT序列与三维网格模型配准方法技术领域
本发明涉及图像处理、计算机视觉和增强现实领域，具体地说是一种口腔软硬组织
CT序列与三维网格模型配准方法。
背景技术
目前，一些科研工作者在相关技术上开展了研究。在图像分割方面，2005年，康涅
狄格州的Chun-MingLi提出了一种基于变分法的改进的水平集(Level-set)方法[5]，此方
法无需对用户标定的初始边界的重新初始化。相对于传统的水平集方法，此方法有以下
好处：能够采用较大的时间戳数值，以提高算法的收敛速度，减少计算时间；能够采用
一般化的标定区域边界作为初始边界，而不受其必须为距离函数的限制；能够使用有限
差分法有效地实现。实验结果表明，该方法即使对相对模糊的物体边界也有较好的分割
效果。但是该方法对用户初始手动标定的粗略边界的容错率较低，容易陷入无法收敛到
正确边界的情况。2010年，康涅狄格州大学的Chun-MingLi采用距离一致化的方法[6]
对[5]进行了改进。在改进的水平集方法中，水平集函数的一致性在其演化过程中保持，
并且由最小化能量函数的距离一致化项的梯度流动方向驱使值为零的水平集边界向正
确的物体轮廓演进。距离一致化项被定义为一个势函数，由此导致的水平集演进的方向
具有独特的FAB(ForwardandBackWard)融合效果，保证了值为零的水平集边界的正确
性。
2010年，韩国庆熙大学的OksamChae基于水平集方法，提出了一种针对牙齿CT
序列的轮廓分割方法[7]。该方法对牙齿的不同部分(牙冠、牙根)分别采用单水平集分
割和多水平集分割的方法，克服了普通水平方法对于轮廓边界过于靠近的不同物体分割
效果不佳的缺点。同时，该方法利用了牙齿CT切片中轮廓边界附近的前后背景灰度作
为先验信息，在传统的水平集方法的能量函数中增添了区域灰度一致性项，以实现能量
函数的最小化时驱动值为零的水平集边界向正确的牙齿边界演进。该方法相对于传统的
水平集方法在牙齿CT序列的分割上有一定改进，但计算速度过慢。
在点集对准方面，迭代最近点(IterativeClosestPoints)算法是一种用于最小化两
个点集之间的差异的配准算法。ICP算法通过迭代计算两个点集之间的匹配关系以及最
小化距离来不断优化两个点集之间的相似程度，从而实现点集配准的效果。
2010年德国卡斯陆工学院提出一种基于专家知识的口腔种植手术中增强现实位置解析
方法。该方法通过获取手术器械的局部上下文信息构建了一个上下文敏感的虚拟现实系统，
通过跟踪人与手术器械的位置识别手术进行的阶段，从而有效的可视化虚实融合的结果。该
方法的示意图如图1所示。
2011年ISMAR会议上，德国慕尼黑工业大学介绍了一个使用增强现实和虚拟现实技术
的手术辅助系统，并详细分析了该系统在100次临床手术操作中的应用情况。该系统的组成
如图2所示。
2006年，瑞典林雪平大学的P.Ljung提出了一种基于交互式的大规模、高解析度的CT
扫描数据的虚拟解剖方法。该方法采用了世界领先的体绘制方法在GPU上实现了上G级的
模型数据的渲染，并利用了采用多解析度渲染技术以及层次细节选取方法的基于变换函数的
数据降维方法，如图3所示。
发明内容
根据上述实际需求和关键问题，本发明的目的在于：提出一种口腔软硬组织CT序
列与三维网格模型配准方法，能够利用头部的CT序列以及人牙齿的三维网格模型，将
三维网格模型配准到CT序列的局部坐标系内。
本发明采用的技术方案是：利用用户编辑标定的牙齿粗略轮廓，采用基于水平集的
图像分割方法，对于输入CT序列进行像素分割，获取CT序列中牙齿轮廓的坐标；在
CT序列序列的局部坐标空间内，通过ICP算法对输入牙齿三维网格模型的点集坐标与
上述步骤得到的牙齿轮廓点集坐标进行配准，得到牙齿三维网格模型在CT序列序列局
部坐标空间内的粗略位置；采用基于拟牛顿法的能量函数优化方法，利用CT序列内牙
齿部分的灰度信息对牙齿三维网格模型进行迭代优化，以获取三维网格模型在CT序列
坐标内的最佳位置。
用户编辑标定牙齿的粗略轮廓，主要是对于CT序列中的若干张图片，通过用户的
直观认识选取确定牙齿轮廓的多边形顶点。采用基于水平集方法的图像分割方法，在用
户标定的每张图片的平面坐标内分割得到牙齿的精细轮廓。
CT序列点集配准，主要是利用CT序列的像素间隔信息以及上述图像分割步骤得到
的牙齿轮廓的精细轮廓得到牙齿轮廓上的点在CT坐标空间内的三维坐标，与牙齿的三
维网格模型进行点集配准，以使牙齿的三维网格模型变换到在CT坐标空间内大致正确
的位置。
三维网格模型的精细配准，主要是利用CT序列的灰度信息，构建出CT序列软硬
组织的灰度梯度矩阵。根据口腔软硬组织的密度差异以及CT序列扫描的特点，在硬组
织牙冠的表面附近其对应的灰度值有较大变化率，即其对应的灰度梯度值为其局部最大
值。通过三维网格模型各顶点所在的CT序列坐标系内对应的CT像素坐标灰度梯度最
大这一约束，采用拟牛顿法对三维网格模型的位置不断迭代优化，从而最终得到三维网
格模型在CT序列坐标系内的精确坐标。
本发明与现有的技术相比，其有益的效果是：1，本发明完全利用了头部的CT序列
数据以及对应的咬模模型模型数据，无需任何其他数据；2，本发明充分利用了头部CT
序列的灰度信息，在利用迭代最近点对算法进行配准的基础上提高了配准精度；3，采
用拟牛顿方法进行迭代优化，能够快速计算出最优的配准位置；
附图说明
图1德国卡斯陆工学院提出的基于专家知识的口腔种植手术中增强现实位置解析方
法。
图2德国慕尼黑工业大学提出的一个使用增强现实和虚拟现实技术的手术辅助系统。
图3瑞典林雪平大学提出的CT扫描数据的虚拟解剖方法。
图4本发明而标注过程显示图。
图5本发明的配准结果截面显示图。
图6本发明的结构流程图；
图7本发明配准算法的流程图。
具体实施方式
下面结合附图对本发明作详细说明。
参阅本发明提出的能量函数定义。口腔软硬组织CT序列与三维网格模型的配准能
量函数是一种基于CT序列的密度分布以及三维网格模型在空间中的旋转、平移变化，
用以描述CT序列与三维网格模型的配准程度的一种能量函数。
设有待配准的CT序列A与三维网格数据B，其中A在三维空间的原点坐标为
{x0,y0,z0}，三个方向的像素间距为{sx0,sy0,sz0}，则A中位于体坐标{x,y,z}将映射到三
维空间内的坐标{x0+sx0x,y0+sy0y,z0+sz0z}。类似，对于B中的任意顶点，假设其在三
维空间内的坐标为{px,py,pz}，则其映射到CT序列内的体坐标为
{ px - x 0 sx 0 , py - y 0 sy 0 , pz - z 0 sz 0 } . ]]>
对于CT序列，设其密度分布函数为D(x,y,z)，其中x,y,z为CT序列的体坐标，可
求得其密度梯度的分布函数依据实际的CT序列的密度分布情况，密度较小
的部分对应病人的软组织、空气等，密度较大的部分对应病人的牙齿硬组织。则CT序
列在病人软硬组织的交界处的密度将会有较大幅度的变化，即对于硬组织轮廓上的点
{x,y,z}，的模取局部极大值。
设有旋转平移变换{α,β,γ,t1,t2,t3}，其中α,β,γ分别代表在x,y,z三个方向的旋转角
度，t1,t2,t3分别代表在x,y,z三个方向的平移量，则对应的旋转矩阵为R＝Rx*Ry*Rz，
其中 R x = 1 0 0 0 cos α - sin α 0 sin α cos α , R y = cos β 0 - sin β 0 1 0 sin β 0 cos β , R z = cos γ - sin γ 0 sin γ cos γ 0 0 0 1 , ]]>可求得
R = cos β cos γ - cos β sin γ - sin β - cos γ sin α sin β + cos α cos γ cos α cos γ + sin α sin β sin γ - cos β sin α cos α cos γ sin β + sin α sin γ cos γ sin α - cos α sin β sin γ cos α cos β . ]]>三维网格模
型上的点{x,y,z}通过该旋转平移变换得到的新点坐标为
{ x ′ , y ′ , z ′ } = t r a n s p o s e ( R * x y z + t 1 t 2 t 3 ) , ]]>计算新点映射到CT序列内的体坐标
则其对应的CT序列密度梯度幅值为
得到配准能量函数的定义此函数反映了在旋转平移变量
{α,β,γ,t1,t2,t3}下三维网格模型与CT序列的配准程度，E的值越低则配准结果越好。因
此，通过迭代优化E的参数{α,β,γ,t1,t2,t3}的值，能够得到{α,β,γ,t1,t2,t3}的最优解，也
就是配准的最佳值。
参阅图7本发明的配准方法流程图。依据能量函数的定义可
得能量函数关于{α,β,γ,t1,t2,t3}的表达式
| ∂ D ( t 1 + cos β cos γ ( - x 0 + vx ′ ) sx 0 - cos β sin γ ( - y 0 + vy ′ ) sy 0 - sin β ( - z 0 + vz ′ ) sz 0 , t 2 + ( - cos γ sin α sin β + cos α sin γ ) ( - x 0 + vx ′ ) sx 0 + ( cos α cos γ + sin α sin β sin γ ) ( - y 0 + vy ′ ) sy 0 - cos β sin α ( - z 0 + vz ′ ) sz 0 , ]]>
t 3 + ( cos α cos γ sin β + sin α sin γ ) ( - x 0 + vx ′ ) sx 0 + ( cos γ sin α - cos α sin β sin γ ) ( - y 0 + vy ′ ) sy 0 + cos α cos β ( - z 0 + vz ′ ) sz 0 ) | ]]>
。设当前的搜索位置为xn＝{αn,βn,γn,t1n,t2n,t3n},可通过E的定义求得xn的各分量的偏导数
以及E关于xn的雅可比行列式Jn，同时
也可求得E关于xn的海森矩阵Hn。线搜索的更新公式为xn+1＝xn-andn，其中dn为线搜
索的方向，且an为搜索的步长。
采用BFGS拟牛顿法进行线搜索时，采用迭代的方法以近似值近似替代
的递推公式为其中yn+1＝xn+1-xn，sn+1＝Jn+1-Jn。
至此，对xn的优化的线搜索方向已经确定，剩下的问题是确定线搜索的步长。根据Strong
Wolfe-Powell准则可以确定合适的搜索步长，使得步长既不太大以至于结果偏离最优值，
也不太小以至于陷入局部极值以及降低搜索效率。
参阅图6本发明的步骤图。
步骤一，显示了需要对咬模模型进行的预处理操作。实际应用中，采用手动模型分
割的方法将咬模模型分割得到不带底座的模型，从而从原始数据中分离得到了排除了非
人体组织结构的噪声数据的干扰，提高了配准的精确度。
步骤二，用户对CT序列各层的牙齿轮廓进行标定。通过读入用户逐点选择的CT
序列内牙齿轮廓的多边形顶点，绘制出其边界。
步骤三，采用水平集方法，利用上一步用户对CT序列的粗略标定结果，对CT序
列进行精确分割，以获取各张CT序列中牙齿的精细轮廓。根据前人的方法，水平集函
数应当保持为有向距离函数，即在指定的尺度空间(Ω,d)内，给定函数f以及空间上的
一个曲面Ωc，满足以下表达式： f ( x ) = d ( x , Ω c ) x ∈ Ω - d ( x , Ω ) x ∈ Ω c , ]]>其中， d ( x , Ω ) = i n f y ∈ Ω d ( x , y ) , ]]>
d(x,Ω)为x到曲面Ω的距离。
为了保证f为有向距离函数，定义来衡量f与真实的有向
距离函数的差异程度，其中为拉普拉斯算子。针对图像分割问题，定义如下边界检测
函数：其中Gσ代表一个以σ为标准差的高斯核，I代表待处理的CT
图像。然后定义能量函数
E ( f ) = μ ∫ Ω 1 2 ( | ▿ f | - 1 ) 2 d x d y + λ ∫ Ω g σ ( f ) | ▿ f | d x d y + ν ∫ Ω g H ( - f ) d x d y , ]]>其中σ为狄拉克函
数，H为Heaviside函数，μ，λ，ν为各项的参数。通过对此能量函数求f的梯度，并
取作为f更新的方向，不断迭代更新直至f收敛，所求得的函数f的零水平线
即为分割得到的边界。结合CT序列的像素间距信息，得到牙齿轮廓在CT序列坐标空
间内的三维坐标。
步骤四，采用迭代最近点算法，对CT序列各层的牙齿轮廓进行精细分割。迭代最
近点算法，即ICP，是IterativeClosestPoint的缩写，用于点集与点集的配准。利用已知
的牙齿三维网格模型的点集坐标，采用迭代最近点算法将三维网格模型对准到CT序列
坐标空间内。由于上一步通过水平集分割方法得到的CT序列各层的牙齿轮廓所包含的
点集与三维网格模型的点集坐标并不相同，这一步通常只能得到三维网格模型的粗略位
置。
步骤五，利用CT序列的灰度信息，构建三维网格模型点集对应于CT序列梯度值
之和的能量方程，以上一步通过迭代最近点算法得到的结果R0,t0作为配准的初始结果，
采用拟牛顿方法对R,t的值进行迭代优化，得到最终的配准位置。
以上所述仅为本发明的一些基本说明，依据本发明的技术方案所做的任何等效变
换，均应属于本发明的保护范围。