CN114757956B - 一种种植体自动提取以及种植体中轴线自动确定方法 - Google Patents

Abstract

本发明提供一种种植体自动提取以及种植体中轴线自动确定方法，包括：获取术后患者种植体位置的CT图像；对滤波处理后的CT图像进行阈值分割，得到二值图；对三维斑块图中的各个三维斑块进行初步过滤，得到初步过滤的三维斑块图；识别三维斑块是否为种植体三维斑块：分别求解种植体斑块的边界盒与该种植体中轴线方程的两个交点，即为提取到的种植体下端点和上端点。本发明提供一种种植体自动提取以及种植体中轴线自动计算的方法，与传统的种植体提取以及种植体中轴线确定方法相比，具有实现简单，计算精度高，实时性强，通过简单的计算就可以分割出来种植体，并对种植体的中轴线进行提取，进而可以快速准确地评判出本次种植手术的种植精度。

Description

一种种植体自动提取以及种植体中轴线自动确定方法

技术领域

本发明属于三维图像分割技术领域，具体涉及一种种植体自动提取以及种植体中轴线自动确定方法。

背景技术

种植牙手术后，需要对本次种植手术后的种植精度和效果进行评判，也就是需要将种植后的种植体的实际位置和术前的规划位置进行对比评判，具体包括实际种植体的上端点和术前规划的种植体上端点、以及实际种植体的下端点和术前规划的种植体下端点进行位置误差的评判，实际种植体的中轴线方向和术前规划的种植体中轴线方向进行姿态误差的评判。因此，需要对术后患者牙骨内的种植体进行分割提取。

目前做法是在二维图像上直接进行种植体的轴向切面分割，此种方法对于种植体倾斜比较厉害的情况是无效的；或者直接使用人工智能中的三维深度学习网络，比如3D-uNet，3D-vNet等对种植体进行三维的分割提取，此种方法训练时需要大量的人工标注好的样本，费时费力，训练时间过长，而且需要配备高性能的计算设备，准确率也不能完全达到100％，并且无法对分割出来的种植体轴线进行计算和提取。

发明内容

针对现有技术存在的缺陷，本发明提供一种种植体自动提取以及种植体中轴线自动确定方法，可有效解决上述问题。

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

本发明提供一种种植体自动提取以及种植体中轴线自动确定方法，包括以下步骤：

步骤S1，获取术后患者种植体位置的CT图像；对所述CT图像进行归一化处理，并对归一化处理后的CT图像进行盒式滤波，得到滤波处理后的CT图像；

步骤S2，对所述滤波处理后的CT图像进行阈值分割，得到经过阈值分割后的二值图；

步骤S3，对于步骤S2得到的二值图，沿着CT扫描方向的轴向转化为三维斑块图；其中，所述三维斑块图中，包含多个三维斑块；

对所述三维斑块图中的各个三维斑块进行初步过滤，滤除明显不属于种植体的三维斑块，得到初步过滤的三维斑块图；

步骤S4，对于步骤S3得到的初步过滤的三维斑块图，共包含m个三维斑块，对于每个三维斑块，表示为三维斑块P_i，均采用以下方法，识别三维斑块P_i是否为种植体三维斑块：

步骤S4.1，求解三维斑块P_i相对于斑块质心O_i的三维转动惯量，得到三维转动惯量矩阵J_i；

步骤S4.1.1，三维斑块P_i共包含n个体素点，任意体素点的三维坐标表示为 (x_j，y_j，z_j)；斑块质心O_i的三维坐标为(x_i，y_i，z_i)；

通过以下公式，分别计算I_xx、I_yy、I_zz、I_xy、I_yx、I_xz、I_zx、I_yz和I_zy：

其中：

三维斑块P_i对应一个等效惯性椭球；该等效惯性椭球具有三个主轴，三个主轴按长度从短到长的方向，分别为x轴，y轴和z轴；

I_xx，I_yy，I_zz分别为绕着x轴，y轴和z轴旋转的主惯性量；

I_xy代表x轴和y轴的耦合惯性积；

I_xz代表x轴和z轴的耦合惯性积；

I_yx代表y轴和x轴的耦合惯性积；

I_yz代表y轴和z轴的耦合惯性积；

I_zx代表z轴和x轴的耦合惯性积；

I_zy代表z轴和y轴的耦合惯性积；

步骤S4.1.2，由此建立三维转动惯量矩阵J_i：

步骤S4.2，对三维转动惯量矩阵J_i进行特征值分解，得到三维转动惯量矩阵 J_i的特征值和特征向量；

其中，特征值为三个，分别表征等效惯性椭球三个主轴的长度；

特征向量V形式为：V＝[v₁，v₂，v₃]；

其中：v₁表示等效惯性椭球最短的主轴的方向向量，称为第一惯性椭球轴；

v₂表示等效惯性椭球第二短主轴的方向向量，称为第二惯性椭球轴；

v₃表示等效惯性椭球长轴的方向向量，称为第三惯性椭球轴；

步骤S4.3，根据步骤S4.2的三个特征值，采用下式，分别建立第一椭球偏心率E₁和第二椭球偏心率E₂：

其中：ρ₁，ρ₂，ρ₃代表由大到小排列的三个特征值；

步骤S4.4，将步骤S4.2得到的最大特征值对应的方向向量作为种植体中轴线的方向向量C_i，并且该方向向量C_i经过斑块质心O_i；

根据种植体中轴线的方向向量C_i和斑块质心O_i，建立种植体中轴线方程W_i；

步骤S4.5，对三维斑块P_i进行CT轴向扫描，得到N个切片；每个切片为二值图，将每个切片表示为S_k，其中，k＝1，2，...，N；

求解切片S_k的质心O_k；

计算质心O_k到种植体中轴线的距离L_i；

计算N个切片的质心到种植体中轴线的距离的平均值L^；

采用下式，计算得到三维斑块P_i切片质心到种植体中轴线距离的方差V：

步骤S4.6，根据种植体细长以及轴对称的几何特点，分别对第一椭球偏心率 E₁、第二椭球偏心率E₂和方差V建立上下限范围，构造形成带通特征滤波器F；

采用带通特征滤波器F对三维斑块P_i进行识别，如果三维斑块P_i在第一椭球偏心率E₁、第二椭球偏心率E₂和方差V的上下限范围内，则三维斑块P_i为种植体斑块；否则，三维斑块P_i不是种植体斑块，滤除；

步骤S5，经过步骤S4，从m个三维斑块中识别得到种植体斑块；确定种植体斑块的边界盒；

分别求解种植体斑块的边界盒与该种植体中轴线方程的两个交点，表示为： A₁和A₂；A₁和A₂即为提取到的CT图像中种植体的下端点和上端点；A₁和A₂的连线为种植体的中轴线线段。

优选的，步骤S3具体为：

根据种植体的几何尺寸，得到种植体所形成的种植体斑块体素数目的上下限，根据所述种植体斑块体素数目的上下限，对所述三维斑块图中的各个三维斑块进行初步过滤，滤除明显不属于种植体的三维斑块，得到初步过滤的三维斑块图。

本发明提供的一种种植体自动提取以及种植体中轴线自动确定方法具有以下优点：

本发明提供一种种植体自动提取以及种植体中轴线自动计算的方法，与传统的种植体提取以及种植体中轴线确定方法相比，具有实现简单，计算精度高，实时性强，通过简单的计算就可以分割出来种植体，并对种植体的中轴线进行提取，进而可以快速准确地评判出本次种植手术的种植精度。

附图说明

图1为使用阈值分割后的三维种植体图像；

图2为根据连通域中体素的数量分割后的图像；

图3为根据三维转动惯量矩阵得到的椭球偏心率和对称度分割后的种植体图像；

图4为惯性矩阵特征向量和边界盒交点；

图5为分割出来的种植体轴线。

具体实施方式

为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

本发明提供一种种植体自动提取以及种植体中轴线自动计算的方法，与传统的种植体提取以及种植体中轴线确定方法相比，具有实现简单，计算精度高，实时性强，通过简单的计算就可以分割出来种植体，并对种植体的中轴线进行提取，进而可以快速准确地评判出本次种植手术的种植精度。

本发明提供一种种植体自动提取以及种植体中轴线自动确定的方法，包括以下步骤：

首先对于种植牙手术后的患者的CT图像进行预处理，包括归一化以及使用盒式滤波的方法；

接着采用阈值分割的方法对术后患者CT图像进行粗分割1，对于粗分割1得到的图像根据26联通域方法，建立连通域，根据种植体的几何尺寸得到连通域中所含像素数量的上下限，根据这一上下限对患者CT图像进行粗分割2，从而完成对患者CT图像的粗分割。

其次将粗分割后得到的二值图像转换到三维视角上来，对于得到的三维二值图像求解三维转动惯量，并对于三维转动惯量进行特征值分解，根据特征向量的方向建立等效的惯性椭球以及椭球的三个主轴长度，从而进一步建立用于图像分割的两个特征即：第一椭球偏心率，第二椭球偏心率。

紧接着根据椭球的最大主轴方向建立椭球中心线方程，将得到的前述粗分割后的三维图像转换到二维图像层中，并对于每一层求解质心，从而求解出来每一层质心到椭球中心线的距离，根据距离建立出来前述距离的方差，从而建立用于图像分割的第三个特征。

根据种植体细长的几何结构，针对上述三个特征，选择合适的三个特征值上下限，从而可以将种植体分割出来。

对于分割出来的种植体建立包络它的最小边界盒，求解上述椭球中心线方程和该边界盒的交点，就可以得到椭球种植体的中轴线两个端点，从而完成对种植牙手术后患者CT图像中种植体的提取和种植体中轴线的确定。

本发明求解过程简单，求解精度高，实时性强，是一种简便快捷的术后患者CT图像种植体的提取方法。

本发明提供一种种植体自动提取以及种植体中轴线自动确定方法，具体包括以下步骤：

步骤S1，获取术后患者种植体位置的CT图像；对所述CT图像进行归一化处理，并对归一化处理后的CT图像进行盒式滤波，得到滤波处理后的CT图像；

步骤S2，对所述滤波处理后的CT图像进行阈值分割，得到经过阈值分割后的二值图；如图1所示，在图1中，除图中标记的两个种植体外，其余斑块分区为头骨中骷骼的干扰。

步骤S3，对于步骤S2得到的二值图，沿着CT扫描方向的轴向转化为三维斑块图；其中，所述三维斑块图中，包含多个三维斑块；

其中，沿着CT扫描方向的轴向转化为三维斑块图，具体为：沿着CT扫描方向的轴向根据26连通域的方法求解三维连通域，得到三维斑块图；如图2所示，为根据连通域中体素的数量分割后的图像；在图2中，除标记的两个种植体外，其余的均为头骨骨骼干扰。

对所述三维斑块图中的各个三维斑块进行初步过滤，滤除明显不属于种植体的三维斑块，得到初步过滤的三维斑块图；

步骤S3具体为：

根据种植体的几何尺寸，得到种植体所形成的种植体斑块体素数目的上下限，根据所述种植体斑块体素数目的上下限，对所述三维斑块图中的各个三维斑块进行初步过滤，滤除明显不属于种植体的三维斑块，完成对三维斑块图的粗过滤，得到初步过滤的三维斑块图。

例如，种植体斑块体素数目的上下限为[h₁，h₂]，如果三维斑块的体素数目在 [h₁，h₂]范围，则初步认为该三维斑块为种植体三维斑块，保留；否则，滤除。

步骤S4，对于步骤S3得到的初步过滤的三维斑块图，共包含m个三维斑块，对于每个三维斑块，表示为三维斑块P_i，均采用以下方法，识别三维斑块P_i是否为种植体三维斑块：

步骤S4.1，求解三维斑块P_i相对于斑块质心O_i的三维转动惯量，得到三维转动惯量矩阵J_i；

步骤S4.1.1，三维斑块P_i共包含n个体素点，任意体素点的三维坐标表示为 (x_j，y_j，z_j)；斑块质心O_i的三维坐标为(x_i，y_i，z_i)；

通过以下公式，分别计算I_xx、I_yy、I_zz、I_xy、I_yx、I_xz、I_zx、I_yz和I_zy：

其中：

三维斑块P_i对应一个等效惯性椭球；该等效惯性椭球具有三个主轴，三个主轴按长度从短到长的方向，分别为x轴，y轴和z轴；

I_xx，I_yy，I_zz分别为绕着x轴，y轴和z轴旋转的主惯性量；

I_xy代表x轴和y轴的耦合惯性积；

I_xz代表x轴和z轴的耦合惯性积；

I_yx代表y轴和x轴的耦合惯性积；

I_yz代表y轴和z轴的耦合惯性积；

I_zx代表z轴和x轴的耦合惯性积；

I_zy代表z轴和y轴的耦合惯性积；

步骤S4.1.2，由此建立3×3的三维转动惯量矩阵J_i：

步骤S4.2，对三维转动惯量矩阵J_i进行特征值分解，得到三维转动惯量矩阵 J_i的特征值和特征向量；

其中，特征值为三个，分别表征等效惯性椭球三个主轴的长度；

特征向量V形式为：V＝[v₁，v₂，v₃]；

其中：v₁表示等效惯性椭球最短的主轴的方向向量，称为第一惯性椭球轴；

v₂表示等效惯性椭球第二短主轴的方向向量，称为第二惯性椭球轴；

v₃表示等效惯性椭球长轴的方向向量，称为第三惯性椭球轴；

步骤S4.3，根据步骤S4.2的三个特征值，采用下式，分别建立第一椭球偏心率E₁和第二椭球偏心率E₂：

其中：ρ₁，ρ₂，ρ₃代表由大到小排列的三个特征值；

步骤S4.4，将步骤S4.2得到的最大特征值对应的方向向量作为种植体中轴线的方向向量C_i，并且该方向向量C_i经过斑块质心O_i；

根据种植体中轴线的方向向量C_i和斑块质心O_i，建立种植体中轴线方程W_i；

步骤S4.5，对三维斑块P_i进行CT轴向扫描，得到N个切片；每个切片为二值图，将每个切片表示为S_k，其中，k＝1，2，...，N；

求解切片S_k的质心O_k；

计算质心O_k到种植体中轴线的距离L_i；

计算N个切片的质心到种植体中轴线的距离的平均值L^；

采用下式，计算得到三维斑块P_i切片质心到种植体中轴线距离的方差V：

步骤S4.6，根据种植体细长以及轴对称的几何特点，分别对第一椭球偏心率 E₁、第二椭球偏心率E₂和方差V建立上下限范围，构造形成带通特征滤波器F；

采用带通特征滤波器F对三维斑块P_i进行识别，如果三维斑块P_i在第一椭球偏心率E₁、第二椭球偏心率E₂和方差V的上下限范围内，则三维斑块P_i为种植体斑块；否则，三维斑块P_i不是种植体斑块，滤除；如图3所示，为根据三维转动惯量矩阵得到的椭球偏心率和对称度分割后的种植体图像；

步骤S5，经过步骤S4，从m个三维斑块中识别得到种植体斑块；确定种植体斑块的边界盒；

分别求解种植体斑块的边界盒与该种植体中轴线方程的两个交点，表示为： A₁和A₂；A₁和A₂即为提取到的CT图像中种植体的下端点和上端点；A₁和A₂的连线为种植体的中轴线线段。如图4和5所示。图4为惯性矩阵特征向量和边界盒交点；图5为分割出来的种植体轴线。

本发明首先将对患者的CT图像进行预处理包括归一化以及阈值分割，并根据26连通域的方法建立三维斑块图，对于三维斑块图根据体素的数量进行粗分割，然后对于粗分割后的三维斑块图建立三维转动惯量矩阵；随后对三维转动惯量矩阵进行特征值分解，根据特征值和特征向量建立等效的惯性椭球，随后对于该惯性椭球进行特征描述，也就是：第一椭球偏心率和第二椭球偏心率；接着根据质心到惯性椭球的距离建立方差特征，根据第一椭球偏心率、第二椭球偏心率他方差这三个特征进行三维种植体分割；可以看到整个过程计算简单，不需要任何的迭代步骤，而且整个分割提取种植体的过程不需要购置昂贵的计算性能高的设备，由于分割种植体和计算种植体轴线都需要种植体等效惯性椭球的长轴，因此该方法是将种植体的分割和种植体轴线的计算融合在一起，因此所提出的方法具有很强的计算实时性；另外该算法引入等效三维转动惯量和特征值分解相结合的方法，使得该算法具有很高的种植体分割和轴线计算精度，可以消除一些离群点以及噪声的干扰。因此，所提出的算法具有很强的鲁棒性和很高的精度。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (2)

1.一种种植体自动提取以及种植体中轴线自动确定方法，其特征在于，包括以下步骤：

步骤S1，获取术后患者种植体位置的CT图像；对所述CT图像进行归一化处理，并对归一化处理后的CT图像进行盒式滤波，得到滤波处理后的CT图像；

步骤S2，对所述滤波处理后的CT图像进行阈值分割，得到经过阈值分割后的二值图；

步骤S3，对于步骤S2得到的二值图，沿着CT扫描方向的轴向转化为三维斑块图；其中，所述三维斑块图中，包含多个三维斑块；

对所述三维斑块图中的各个三维斑块进行初步过滤，滤除明显不属于种植体的三维斑块，得到初步过滤的三维斑块图；

步骤S4，对于步骤S3得到的初步过滤的三维斑块图，共包含m个三维斑块，对于每个三维斑块，表示为三维斑块P_i，均采用以下方法，识别三维斑块P_i是否为种植体三维斑块：

步骤S4.1，求解三维斑块P_i相对于斑块质心O_i的三维转动惯量，得到三维转动惯量矩阵J_i；

步骤S4.1.1，三维斑块P_i共包含n个体素点，任意体素点的三维坐标表示为(x_j，y_j，z_j)；斑块质心O_i的三维坐标为(x_i，y_i，z_i)；

通过以下公式，分别计算I_xx、I_yy、I_zz、I_xy、I_yx、I_xz、I_zx、I_yz和I_zy：

其中：

三维斑块P_i对应一个等效惯性椭球；该等效惯性椭球具有三个主轴，三个主轴按长度从短到长的方向，分别为x轴，y轴和z轴；

I_xx，I_yy，I_zz分别为绕着x轴，y轴和z轴旋转的主惯性量；

I_xy代表x轴和y轴的耦合惯性积；

I_xz代表x轴和z轴的耦合惯性积；

I_yx代表y轴和x轴的耦合惯性积；

I_yz代表y轴和z轴的耦合惯性积；

I_zx代表z轴和x轴的耦合惯性积；

I_zy代表z轴和y轴的耦合惯性积；

步骤S4.1.2，由此建立三维转动惯量矩阵J_i：

步骤S4.2，对三维转动惯量矩阵J_i进行特征值分解，得到三维转动惯量矩阵J_i的特征值和特征向量；

其中，特征值为三个，分别表征等效惯性椭球三个主轴的长度；

特征向量V形式为：V＝[v₁，v₂，v₃]；

其中：v₁表示等效惯性椭球最短的主轴的方向向量，称为第一惯性椭球轴；

v₂表示等效惯性椭球第二短主轴的方向向量，称为第二惯性椭球轴；

v₃表示等效惯性椭球长轴的方向向量，称为第三惯性椭球轴；

步骤S4.3，根据步骤S4.2的三个特征值，采用下式，分别建立第一椭球偏心率E₁和第二椭球偏心率E₂：

其中：ρ₁，ρ₂，ρ₃代表由大到小排列的三个特征值；

步骤S4.4，将步骤S4.2得到的最大特征值对应的方向向量作为种植体中轴线的方向向量C_i，并且该方向向量C_i经过斑块质心O_i；

根据种植体中轴线的方向向量C_i和斑块质心O_i，建立种植体中轴线方程W_i；

步骤S4.5，对三维斑块P_i进行CT轴向扫描，得到N个切片；每个切片为二值图，将每个切片表示为S_k，其中，k＝1，2，...，N；

求解切片S_k的质心O_k；

计算质心O_k到种植体中轴线的距离L_i；

计算N个切片的质心到种植体中轴线的距离的平均值L^；

采用下式，计算得到三维斑块P_i切片质心到种植体中轴线距离的方差V：

步骤S4.6，根据种植体细长以及轴对称的几何特点，分别对第一椭球偏心率E₁、第二椭球偏心率E₂和方差V建立上下限范围，构造形成带通特征滤波器F；

采用带通特征滤波器F对三维斑块P_i进行识别，如果三维斑块P_i在第一椭球偏心率E₁、第二椭球偏心率E₂和方差V的上下限范围内，则三维斑块P_i为种植体斑块；否则，三维斑块P_i不是种植体斑块，滤除；

步骤S5，经过步骤S4，从m个三维斑块中识别得到种植体斑块；确定种植体斑块的边界盒；

分别求解种植体斑块的边界盒与该种植体中轴线方程的两个交点，表示为：A₁和A₂；A₁和A₂即为提取到的CT图像中种植体的下端点和上端点；A₁和A₂的连线为种植体的中轴线线段。

2.根据权利要求1所述的一种种植体自动提取以及种植体中轴线自动确定方法，其特征在于，步骤S3具体为：

根据种植体的几何尺寸，得到种植体所形成的种植体斑块体素数目的上下限，根据所述种植体斑块体素数目的上下限，对所述三维斑块图中的各个三维斑块进行初步过滤，滤除明显不属于种植体的三维斑块，得到初步过滤的三维斑块图。