说明书建立釉牙本质界的统计平均模型的方法
技术领域
本发明涉及医学图像处理领域，具体地说涉及一种用于建立釉牙本质界的统计平均模型的方法。
背景技术
天然牙颜色和半透明渐变的美学效果是牙釉质和牙本质以釉牙本质界为交界面叠加产生的。模拟牙齿结构和光学特点研制一种能被CAD（Computer Aided Design，计算机辅助设计）/CAM（Computer Aided Manufacturing，计算机辅助制造）加工的光学颜色梯度渐变氧化锆陶瓷，是解决后牙纯氧化锆全瓷冠不能兼顾美观和强度问题的有效途径。而了解釉牙本质界完整曲面形态是对其进行结构仿真的基础。关于釉牙本质界形态的研究，文献[1]和[2]通过物理组织切片法研究牙齿某个轴面上釉牙本质界的形态，其中获得的釉牙本质界是线而不是面。文献[3]通过Micro-CT（Computed Tomography，计算机断层造影），文献[4]-[6]通过酸脱矿选择性去除釉质的方法获得完整的釉牙本质界曲面。但是因为Micro-CT射线剂量大，而酸脱矿的方法是破坏性操作，所以这两种方法只能用于离体牙，而不能用于活体天然牙，因此很难获得建立釉牙本质界数据库所需的样本量。文献[7]使用锥形束CT扫描并结合图像分割技术，可以在不破坏牙齿的条件下获得釉牙本质界表面形态。
文献[8]尝试在全瓷瓷块中构建类似天然牙冠的三维叠层结构，其中该瓷块结构是：内部瓷模拟牙本质光学特点，外部瓷模拟牙釉质光学特点，二者的结合界面是圆锥形，对天然牙冠的半透明和颜色渐变效果有一定程度的模拟。但是其釉牙本质界采用圆锥形是粗糙的。
统计形状模型是根据不同时间同一个体或同一时间不同个体间确定感兴趣解剖结构表面上对应的标志点集建立的模型。统计形状模型构建的前提是获得解剖结构间对应标志点集。文献[9]描述了一种球面调和函数（Spherical Harmonic Description，SPHARM），其利用参数空间的归一化确定结构的对应关系。但是采用SPHARM进行统计建模与分析需要输入网格是零亏格表面(genus-zero surfaces)，即需要模型具有球形拓扑结构。而釉牙本质界表面是上皮和外间充质两种不同矿化组织的交界面，其外形为贝壳状的不封闭表面。
发明内容
因此，需要实现一种能够适用于非封闭表面的统计平均模型建立方法。
根据本发明的一个方面，提供一种用于建立釉牙本质界的统计平均模型的方法，包括步骤：获得牙齿的CT图像数据；对所述CT图像数据进行分割，以获得釉牙本质界的表面；使用基于曲率的聚类算法对所获得的表面进行分割，以分割出釉牙本质界的底面；通过球面调和分析对分割出底面的釉牙本质界的表面进行球面参数化；以及将该牙齿的不同样本对齐，以获得该统计平均模型。
根据本发明的实施方式，采用水平集算法对所述CT图像数据进行分割。
根据本发明的另一实施方式，所述聚类算法是K-Means算法。
根据本发明的另一实施方式，采用优化的面积和长度失真控制算法来进行所述球面参数化。
根据本发明的实施方式，通过添加权重来优化所述面积和长度失真控制算法。
根据本发明的实施方式，还包括在完成所述球面参数化之后对参数化后的表面进行球谐傅立叶展开。
根据本发明的实施方式，所述对齐在球坐标系下进行。
根据本发明的实施方式，所述对齐采用使用最近点迭代进行SPHARM配准(SHREC)算法来执行。
使用本发明的平均模型建立方法，可以在非封闭表面上使用球谐函数进行物体表面的描述，完成结构间标志点集的对应关系建立。从而如果使用锥形束CT扫描离体牙样本，可以在不破坏牙齿的条件下获得描述釉牙本质界表面形态的统计平均模型。
附图说明
图1示出根据本发明方法的一个实施例的流程图。
图2示出根据本发明方法的一个实施例分割获得的釉牙本质界的表面的示意图。
图3示出根据本发明方法的一个实施例分割出底面的釉牙本质界的表面的示意图。
图4示出根据本发明方法的一个实施例的球面参数化与未通过本发明方法的球面参数化的比较。
图5示出根据本发明方法的一个实施例求取的平均模型的示意图。
具体实施方式
图1示出根据本发明方法的一个实施例的流程图。
在获得牙齿的CT（Computed Tomography，计算机断层造影）图像数据之后，首先在步骤S100中对所获得CT图像数据进行图像分割，去除表层釉质，重建釉牙本质界表面的三维形态。在一个实施例中，作为分割算法例如采用专业人员已知的基于稀疏场 (Sparse Field)的水平集 (Level Set)算法。
在分割CT图像时，感兴趣区域立方体对牙本质的截取会使本应呈开放形态的釉牙本质界表面封闭。但是如果将底面一起用于求釉牙本质界的统计平均模型，将大大降低平均模型建立的准确性。因此在本发明的该实施例中，需要分割出底面和可能存在的牙髓腔，从而在统计建模过程中降低底面和可能存在的牙髓腔对求取统计平均模型的影响。
因此在重建了釉牙本质界表面的三维形态之后，在步骤S200中使用基于曲率的聚类算法来分割出釉牙本质界的底面。
在本发明的一个实施例中，用三角网格模型来模拟釉牙本质界表面。该三角网格模型 可由一对线性集合表示：
，
其中是顶点集合，为顶点个数；是三角面片集合，为三角面片个数。
为了实施基于曲率的聚类算法，需要计算出该三角网格模型的顶点处的法向矢量，然后基于这些法向矢量计算三角网格顶点处的曲率。
顶点处的法向矢量的计算如下：
假定三角网格模型的顶点为，除顶点以外的顶点组成的集合记为。如果顶点，则是的相邻点。中的顶点个数记为。包含的三角面片集合记为。若三角面片，则与是相关的，中的三角面片个数记为。
为了反映三角面片形状对顶点处法向量的影响，采用等式（1）计算顶点的法向量：
（1）
其中，是在顶点处的内角，Ak为三角面片的面积，为三角面片的法向量，其由等式（2）确定：
（2）
其中，表示由顶点指向顶点的边矢量。
计算出顶点处的法向矢量之后，接着如下计算所述三角形网格顶点处的曲率：
首先求出顶点处的切平面，然后将该顶点区域内的邻点投影到所述切平面上。顶点处的法向量为，以为法向量且通过点的平面即为顶点处的切平面，记为。顶点区域中的顶点在该切平面上的投影记为。
点在切平面上沿方向的单位切矢为：
（3）
其中表示两个向量的点积。
通过该单位切矢来计算对称矩阵：
（4）
其中，权值，是点处的单位切矢，是点处沿的方向曲率。
对称矩阵的两个特征值分别为和。文献[10]描述了通过这两个特征值来获得点处的主曲率的等式：
（5）
由此获得三角网格的每一个顶点的主曲率值和，及其对应的曲率方向和。
在获得顶点的曲率之后，可以基于曲率对顶点进行聚类分割。在本发明的一个实施例中，使用K-Means算法对三角网格上的顶点根据主曲率值进行基于距离(即曲率空间欧氏距离)的聚类。首先随机选取任意两个对象作为初始聚类的中心，作为初始簇。在每次迭代中对由主曲率向量构成的数据集进行聚类。特别是，通过对数据集中剩余的每个对象，根据其与各个簇中心的距离将每个对象重新赋给最近的簇。当处理完所有数据对象后，一次迭代结束，获得新的聚类中心。在算法结束时，三角网格上的每一个顶点都被分类归属于一个簇。通过该基于曲率的聚类算法，可以从釉牙本质界的表面中分割出底面来。图3示出通过该聚类算法将底面从釉牙本质界表面中分割出的示意图。颜色最深的部分为分割出的底面，其由三角形网格组成。
在将底面剔除之后，在步骤S300中通过球面调和分析将剔除底面之后的釉牙本质界的表面进行球面参数化。
在该过程中，首先使用优化的CALD（Control of Area and Length Distortions，面积和长度失真控制）算法完成面片模型的球面参数化。文献[9]对CALD算法进行了描述。但在本发明的方法中，对该CALD算法进行了优化。
对于整个三角网格模型，对其进行参数化过程中的面积变形代价值为：
（6）
其中，表示笛卡尔坐标空间中原始模型的三角形集合。为一个可逆映射，将三角形映射为参数空间中的网格，。A(﹒)表示三角形或参数空间中网格的面积。
在此，表示权重，其由等式（7）定义：
（7）
其中，为底面三角形构成的集合；，为原始模型中除底面以外的三角形集合。
通过加入权重项，降低了底面区域在参数化过程中对结果的影响。图4中的（a）示出未加入权重项的球面参数化结果，而图4（b）示出加入权重项后的球面参数化结果。从图4中可以看出，对于非底面的区域，球面映射使得三角形网格连续均匀地映射到单位球上；而对于底面区域，由于保证了非底面区域映射的三角形面积失真最小，此区域映射的面积失真会变大。
在完成球面参数化后，还对该参数化的结果进行球谐傅里叶展开以获得对应的球面调和函数展开形式。球面调和函数是定义在球面坐标系的一组函数，构成球面上的一组标准正交基。将球面调和描述子参数化，使得任何一个单连通的曲面拓扑结构都等价于球面。球面上的点可以通过球面调和函数的线性组合来表示：
（8）
其中，可通过最小二乘估算求出。是对原始系数的估计，可以重构出原函数如下：
其中为用户指定的最大自由度。越大，还原出的就越精确。
通过等式（8）对三角网格的顶点进行了球谐重构，从而实现了对该模型的参数化过程。
最后在步骤S400，将同一牙齿个体的不同样本进行对齐，以求取该个体的统计平均模型。图5示出前磨牙釉牙本质界的统计平均模型。
在本发明的实施例中，采用由文献[11]已知的SHREC（SPHARM Registration with Iterative Closest Point，利用迭代最近点的SPHARM配准）算法进行样本模型空间和参数空间的对齐(Alignment)，从而建立不同样本模型间的对应关系(Correspondence)。在一个示例中，为了建立第一前磨牙釉牙本质界的统计平均模型，可以收集离体牙样本45个。使用该SHREC算法可以实现不同个体间同一结构对应标志点集的自动定位。
通过本发明的方法，可以对非封闭表面建立统计平均模型。此外，本方法可以对其它对象的锥形束CT图像获得相应的三维统计平均模型。
本发明的方法可以通过计算机软件、计算机硬件或它们的组合来实施。
在此公开的实施例仅仅以示例方式描述了用于对CT图像数据进行融合绘制的方法。本领域技术人员可以理解，在不脱离所附权利要求所限定的精神和范围的情况下，可以对所公开的实施例进行各种变化和修改。
参考文献
[1] Anthony J. Olejniczak , Lawrence B, etal. Quantification of dentine shape in anthropoid primates[J]. Ann Anat ,2004,186: 479-485.
[2] Smith T.M, Olejniczak A.J, Reid D.J, etal. Modern human molar enamel thickness and enamel-dentine junction shape[J]. Archives of Oral Biology , 2006, 51: 974-995.
[3] Rodrigo Borges Fonseca, Francisco Haiter-Neto, Alfredo J. Fernandes-Neto,etal. Radiodensity of enamel and dentin of human,bovine and swine teeth[J]. Archives of Oral Biology,2004,49:919-922.
[4] Robert S, Corruceini. The Dentinoenamel Junction in Primates[J]. International Journal of Primatology, 1987, 8(2): 99-114.
[5] Bertram S, Kraus. Morphologic relationships between enamel and dentin surfaces of lower first molar teeth[J].J.D.Res,1952,31(2):248-256.
[6] Panaghlotis Bazos, Pascal Magne. Bio-emulation:biomimetically emulating nature utilizing a histo-anatomic approach;structural analysis[J]. The European journal of esthetic dentistry, 2011, 6(1): 8-19.
[7] 王晓静. 锥形束CT初建釉牙本质界表面三维形态的可行性研究[D]. 北京: 北京大学医学部, 2012.
[8] Kurbad A. Three-dimensionally layered ceramic blocks[J]. Int J Comput Dent, 2010,13(4): 351-65.
[9] Shen L, Makedon FS. Spherical mapping for processing of 3-D closed surfaces. Image and Vision Computing, 24(7): 743-761, 2006.
[10] G. Taubin. Estimating the tensor of curvature of a surface from a polyhedral approximation[A]. In: WEL Grimson ed. Proceedings of the Fifth International Conference on Computer Version[C], Los Alamitos: IEEE Computer Society Press, 1995: 902-907.
[11] Shen, L., Huang, H., Makdeon, F., Saykin, A.J.: Efficient registration of 3D spharm surface. In: 4th Canadian Conf. on Computer & Robot Vision. (2007) 81-88。