一种基于空间中心轴的髁突CBCT影像数据增强方法

一种基于空间中心轴的髁突cbct影像数据增强方法
技术领域
1.本发明属于医疗图像处理领域，尤其涉及了一种基于空间中心轴的髁突cbct影像数据增强方法。


背景技术：

2.颞下颌关节紊乱(tmd)是口腔颌面部最常见的疾病，是发病机制尚未完全明了，颞下颌关节(tmj)不能发挥正常功能的一组疾病的总称。且发病率高，在常见的口腔疾病中排在第四位。颞下颌关节炎(tmj-oa)是tmd的一种亚型，可导致严重的关节疼痛、功能障碍、牙错咬合和与健康相关的生活质量下降。
3.髁突又称髁状突或关节突，是下颌骨主要生长中心之一。髁突的结构可以分为髁突头与髁突颈两部，髁突的上端膨大为髁突头，髁突下部缩小部分称为髁突颈。其中沿下颌骨方向通过髁突中心的一条假想线，称为髁突的中心轴；髁突头上最长的直径称为髁突头的长轴，最短直径称为髁突头的短轴。在tmj-oa患者中发现其髁突形态学发生了改变，包括骨皮质的吸收变薄，骨小梁的增多增粗，骨小梁间囊样变形成等。所以可以通过分析髁突结构的变化，对tmj-oa进行诊断。此外，锥束计算机断层扫描(cbct)是一种成熟的头颈部区域的成像方式，由于它能三维(轴位、冠状位和矢状位)显示正常组织结构和病变组织，避免了二维图像上影像重叠等缺点，因此它被推荐为诊断tmj-oa的最可靠的方法之一，可以通过影像观察到髁突结构的完整性。
4.近年来，深度学习技术在计算机视觉应用中取得了很多突破，包括医学图像的分割。这一成功促使研究人员计划采用深度学习模型对cbct影像中的髁突进行识别并进行tmj-oa的诊断。然而，与其他计算机视觉任务相比，医学图像采集和标注的成本很高，由于cbct影像是三维的，在标注时需要医生逐片追踪髁突轮廓进行标注，比二维图像的标注工作更繁琐耗时。因此需要通过处理原始的髁突cbct影像以扩大可用的样本数量，从而减少医生标注的时间，提高工作效率。
5.然而，现阶段利用深度学习研究髁突分割的研究很少，所以通常采用通用的数据增强方式对髁突cbct影像进行处理。在现有通用的数据增强方法中，如：旋转、翻转、裁剪等几何变换，并不完全适用于髁突的cbct影像数据，因为区别于二维图像，在cbct这样的三维影像上，要考虑髁突在cbct影像中的空间位置关系。因此借助几何变换的数据增强方法得到的扩充后的数据训练模型，可能会导致模型效果比扩充前训练的模型效果还要差。此外，不同人之间髁突所在的相对空间位置并不相同，若采用固定的处理参数，可能会产生影像内没有髁突的感兴趣区域这一类数据，而这一类数据并不能提升髁突分割模型的分割效果。由于髁突cbct影像中髁突本身感兴趣区域比较小，而采用加椒盐噪声，高斯噪声等像素变换的数据增强方式，会对髁突分割任务产生干扰。所以目前适用于髁突cbct影像的数据增强方式是存在缺陷的。


技术实现要素：

6.为了克服已有技术的不足，本发明提供了一种提供一种基于空间中心轴的髁突cbct影像数据增强方法，能够在原始数据量较少的情况下通过找到髁突的中心轴对髁突cbct影像进行几何变换，从而生成可以作为可用样本的数据，以使髁突的cbct影像训练样本集得到数据增强，进而提升利用数据增强后髁突的cbct影像训练样本集进行髁突分割的精度以及模型的泛化能力。
7.本发明解决其技术问题所采用的技术方案是：
8.一种基于空间中心轴的髁突cbct影像数据增强方法，所述方法包括以下步骤：
9.s1、将作为训练集的髁突cbct影像以及人工标注结果d
train1
输入至中心轴计算模块中，得到d
train1
中所有髁突的中心轴，组成一个中心轴集l
center
；
10.s2、将训练集d
train1
和中心轴集l
center
输入至旋转模块中，得到处理后的数据集d
train2
；
11.s3、将d
train2
输入至长轴计算模块中，得到d
train2
中所有的长轴，组成一个长轴集l
long
；
12.s4、将训练集d
train2
和长轴集l
long
输入至旋转模块中，得到最终用于训练髁突分割模型的数据集d
train3
。
13.进一步，所述步骤s1中，将d
train1
输入至中心轴计算模块中，得到d
train1
中所有髁突的中心轴，组成一个中心轴集l
center
，过程如下：
14.步骤1.1、将大小为d1×
h1×
w1的人工标注结果m1(m1∈d
train1
)沿着z轴方向得到d1张大小为d1×
h1×
w1的人工标注结果的切片，得到切片集s1；
15.步骤1.2、将切片集s1中每张切片进行图像灰度化和二值化处理，得到d1张二值化后的图像，组成新的切片集s2；
16.步骤1.3、根据切片集s2，借助cv2库中的轮廓检测算法，得到所有切片中人工标注区域的轮廓的点坐标集p1；
17.步骤1.4、根据每张人工标注结果的切片s
1，i
(s
1，i
∈s2，i∈[0，d
1-1])中的轮廓点坐标集p
1，i
(p
1，i
∈p1，i∈[0，d
1-1])计算相应的几何中心，而轮廓的几何中心坐标可以根据切片s
1，i
的零阶矩与一阶矩计算得到；
[0018]
切片s
1，i
的(p q)阶矩的计算公式如下：
[0019][0020]
其中p
1i_j
(x，y)表示切片s
1i
的轮廓点坐标集p
1，i
中的一个点坐标(p
1，i_j
∈p
1，i
，i∈[0，d
1-1]，j＝0，1，2，...)，x、y分别表示轮廓点坐标集p
1，i
中所有x轴的坐标值和y轴坐标值；
[0021]
切片s
1i
轮廓的几何中心p
c，i
(c
x
，cy)的坐标计算公式如下：
[0022]cx
＝m
10
/m
00
[0023]cy
＝m
01
/m
00
[0024]
其中m
00
是二值化后的切片s
1，i
的零阶矩，表示切片中所有白色区域即像素值为1的和，m
10
是二值化后的切片s
1，i
的一阶矩，表示切片中所有白色区域像素的x坐标的和，m
01
是二值化后的切片s
1，i
的一阶矩，表示切片中所有白色区域像素的y坐标的和；
[0025]
步骤1.5、将切片集s2中所有的切片均进行步骤1.4的操作，得到每张切片中人工标注结果的轮廓几何中心，最后得到一个几何中心集p2(p
c，i
∈p2，i∈[0，d
1-1])；
[0026]
步骤1.6、根据得到的几何中心集p2拟合一条空间直线，其中空间直线标准方程的表达式如下：
[0027][0028]
将公式(1)进行转换，得到新的表达式如下：
[0029][0030]
其中m，n，k，x0，y0，z0均为常数，所以可以将公式(2)进行简化，得到新的表达式如下：
[0031][0032]
其中k1，k2，b1，b2表示常数；
[0033]
可以利用最小二乘法拟合原理，得到常数k1，k2，b1，b2的最佳值，其中最小二乘法拟合原理如下：拟合公式(3)的空间直线表达式，并使得残差平方和最小，其中残差平方和的公式如下：
[0034]
q1＝σ(x
i-(k1·
zi b1))2[0035]
q2＝σ(y
i-(k2·
zi b2))2[0036]
其中zi表示切片s
1，i
在人工标注结果m1中所对应的z轴上的坐标值，xi、yi分别表示切片s
1，i
几何中心p
c，i
的x轴和y轴上所对应的坐标值；
[0037]
当q1、q2最小，即残差平方和最小时，即可得到常数k1，k2，b1，b2的最佳值，因为q1、q2是关于k1，k2，b1，b2的函数，当q1、q2的导数为0时，可取到相应的最小值，所以，分别对k1，k2，b1，b2求偏导可以得到如下公式：
[0038][0039]
从而得到k1，k2，b1，b2的表达式如下：
[0040][0041][0042]
遍历几何中心集p2中所有几何中心，带入k1，k2，b1，b2的表达式中，求得常数k1，k2，
b1，b2的最佳值，最后即可得到拟合的公式(3)的空间直线表达式，拟合得到的空间直线即为m1的中心轴l
center，i
；
[0043]
步骤1.7、将训练集d
train1
中的所有人工标注结果进行步骤1.1～步骤1.6的操作，最后得到一个中心轴数据集l
center
(l
center，i
∈l
center
，i∈[0，d
1-1])。
[0044]
再进一步，所述步骤s2中，将训练集d
train1
和中心轴集l
center
输入至旋转模块中，得到处理后的数据d
train2
，过程如下：
[0045]
步骤2.1、根据中心轴计算模块得到的中心轴l
center，i
(l
center，i
∈l
center
，i＝0，1，2，...)在大小为d1×
h1×
w1的人工标注结果m1中的空间位置，计算中心轴l
center，i
与z轴的夹角，根据夹角，将髁突cbct影像i1以及人工标注结果m1(i1∈d
train1
，m1∈d
train1
)进行旋转，使得髁突的中心轴l
center，i
旋转至z轴，并进行重采样操作，最后得到大小仍为d1×
h1×
w1的新的髁突cbct影像i2以及人工标注结果m2；
[0046]
步骤2.2、将训练集d
train1
中的所有髁突cbct影像以及人工标注结果进行步骤2.1的操作，最后得到处理后的数据d
train2
(i2∈d
train2
，m2∈d
train2
)。
[0047]
更进一步，所述步骤s3中，将d
train2
输入至长轴计算模块中，得到d
train2
中所有的长轴，组成一个长轴集l
long
，过程如下：
[0048]
步骤3.1、给定一个阈值t1，将人工标注结果m2(m2∈d
train2
)的z轴坐标值小于t1的部分舍弃，从而得到仅含髁突头的感兴趣区域i
hroi
；
[0049]
步骤3.2、将髁突头的感兴趣区域i
hroi
沿着z轴方向得到d1张大小为d1×
h1×
w1的髁突头的感兴趣区域的切片，得到切片集s3；
[0050]
步骤3.3、将切片集s3中每张切片进行图像灰度化和二值化处理，得到d1张二值化后的图像，组成新的切片集s4；
[0051]
步骤3.4、根据切片集s4，借助cv2库中的轮廓检测算法，得到所有切片中髁突头的感兴趣区域的轮廓的点坐标集p3；
[0052]
步骤3.5、根据每张髁突头的感兴趣区域的切片s
2，i
(s
2，i
∈s4，i∈[0，d
1-1])中的轮廓点坐标集p
2，i
(p
2，i
∈p3，i∈[0，d
1-1])，记录轮廓中线段最长的两点的坐标，即记录髁突长轴的坐标，最后得到d1个线段组成的线段集l；
[0053]
步骤3.6、根据得到的线段集l拟合一条平面直线，其中平面直线方程的一般表达式如下：
[0054]
y＝k3·
x b3(4)
[0055]
根据最小二乘法原理，得到常数k3，b3的最佳值，其中最小二乘法拟合原理如下：拟合公式(4)的平面直线表达式，并使得残差平方和最小，其中残差平方和的公式如下：
[0056]
q3＝σ(k3·
xi b
3-yi)2[0057]
其中xi、yi分别表示切片s
2，i
上髁突长轴的x轴和y轴的坐标值，
[0058]
当q3最小，即残差平方和最小时，即得到常数k3，b3的最佳值，因为q3是关于k3，b3的函数，当q3的导数为0时，取到相应的最小值，所以，分别对k3，b3求偏导可以得到如下公式：
[0059]
[0060]
从而得到k3，b3的表达式如下：
[0061][0062]
其中表示的算术平均值，计算公式如下：
[0063][0064]
遍历线段集l中所有长轴，带入k3，b3的表达式中，求得常数k3，b3的最佳值，最后即可得到拟合的公式(4)的平面直线表达式，拟合得到的平面直线即为m2的髁突头的长轴；
[0065]
步骤3.7、将训练集d
train2
中的所有人工标注结果进行步骤3.1～
[0066]
步骤3.6的操作，最后得到d
train2
中髁突头的长轴数据集l
long
。
[0067]
所述步骤s4中，将训练集d
train2
和长轴集l
long
输入至旋转模块中，得到最终用于训练髁突分割模型的数据集d
train3
，过程如下：
[0068]
步骤4.1、根据长轴计算模块得到的髁突头长轴li(li∈l
long
，i＝0，1，2，...)在大小为d1×
h1×
w1的人工标注结果m2中的空间位置，计算长轴与x轴的夹角，根据夹角，将髁突cbct影像i2以及人工标注结果m2(i2∈d
train1
，m2∈d
train1
)进行旋转，使得髁突的长轴li旋转至x轴，并进行重采样操作，最后得到大小仍为d1×
h1×
w1的新的髁突cbct影像i3以及人工标注结果m3；
[0069]
步骤4.2、将训练集d
train2
中的所有髁突cbct影像以及人工标注结果进行步骤4.1的操作，最后得到处理后的数据d
train3
(i3∈d
train3
，m3∈d
train3
)。
[0070]
本技术提出的一种基于空间中心轴的髁突cbct影像数据增强方法，首先通过数据划分，得到用于髁突分割模型训练的数据。然后通过数据预处理，排除其他软组织的干扰。最后通过中心轴计算模块和图像旋转处理模块进行数据增强，最终得到数据增强后髁突的cbct影像训练样本集。
[0071]
本发明的有益效果主要表现在：提升利用数据增强后髁突的cbct影像训练样本集进行髁突分割的精度以及模型的泛化能力。
附图说明
[0072]
图1是本技术基于空间中心轴的髁突cbct影像数据增强方法流程图；
[0073]
图2是本技术基于空间中心轴的髁突cbct影像数据增强方法相关模块示意图；
[0074]
图3是本技术的髁突原始数据与进行数据增强后的数据在相同空间位置上的拟样例对比图，其中，(a)为髁突原始数据，(b)为经过本技术提出的方法得到的影像。
具体实施方式
[0075]
下面结合附图对本发明作进一步描述。
[0076]
参照图1～图3，一种基于空间中心轴的髁突cbct影像数据增强方法，包括以下步骤：
[0077]
步骤s1、将作为训练集的髁突cbct影像以及人工标注结果d
train1
输入至中心轴计算模块中，得到d
train1
中所有髁突的中心轴，组成一个中心轴集l
center
。
[0078]
本技术所述训练集d
train1
是由原始髁突cbct图像以及人工标注结果按照一定比例
进行划分，包括设置训练集、验证集以及测试集中数据数量的比例得到的，同时对d
train1
中的所有数据进行了预处理，包括图像裁剪和图像归一化处理的操作，如图2所示，所述将d
train1
输入至中心轴计算模块中，得到d
train1
中所有髁突的中心轴，组成一个中心轴集l
center
，过程如下：
[0079]
步骤1.1、将大小为d1×
h1×
w1的人工标注结果m1(m1∈d
train1
)沿着z轴方向得到d1张大小为d1×
h1×
w1的人工标注结果的切片，得到切片集s1；
[0080]
步骤1.2、将切片集s1中每张切片进行图像灰度化和二值化处理，得到d1张二值化后的图像，组成新的切片集s2；
[0081]
步骤1.3、根据切片集s2，借助cv2库中的轮廓检测算法，得到所有切片中人工标注区域的轮廓的点坐标集p1；
[0082]
步骤1.4、根据每张人工标注结果的切片s
1，i
(s
1，i
∈s2，i∈[0，d
1-1])中的轮廓点坐标集p
1，i
(p
1，i
∈p1，i∈[0，d
1-1])计算相应的几何中心，而轮廓的几何中心坐标可以根据切片s
1，i
的零阶矩与一阶矩计算得到；
[0083]
切片s
1，i
的(p q)阶矩的计算公式如下：
[0084][0085]
其中p
1i_j
(x，y)表示切片s
1i
的轮廓点坐标集p
1，i
中的一个点坐标(p
1，i_j
∈p
1，i
，i∈[0，d
1-1]，j＝0，1，2，...)，x、y分别表示轮廓点坐标集p
1，i
中所有x轴的坐标值和y轴坐标值；
[0086]
切片s
1i
轮廓的几何中心p
c，i
(c
x
，cy)的坐标计算公式如下：
[0087]cx
＝m
10
/m
00
[0088]cy
＝m
01
/m
00
[0089]
其中m
00
是二值化后的切片s
1，i
的零阶矩，表示切片中所有白色区域即像素值为1的和，m
10
是二值化后的切片s
1，i
的一阶矩，表示切片中所有白色区域像素的x坐标的和。m
01
是二值化后的切片s
1，i
的一阶矩，表示切片中所有白色区域像素的y坐标的和；
[0090]
步骤1.5、将切片集s2中所有的切片均进行步骤1.4的操作，得到每张切片中人工标注结果的轮廓几何中心，最后得到一个几何中心集p2(p
c，i
∈p2，i∈[0，d
1-1])；
[0091]
步骤1.6、根据得到的几何中心集p2拟合一条空间直线，其中空间直线标准方程的表达式如下：
[0092][0093]
将公式(1)进行转换，得到新的表达式如下：
[0094][0095]
其中m，n，k，x0，y0，z0均为常数，所以可以将公式(2)进行简化，得到新的表达式如下：
[0096]
[0097]
其中k1，k2，b1，b2表示常数；
[0098]
可以利用最小二乘法拟合原理，得到常数k1，k2，b1，b2的最佳值。其中最小二乘法拟合原理如下：拟合公式(3)的空间直线表达式，并使得残差平方和最小，其中残差平方和的公式如下：
[0099]
q1＝σ(x
i-(k1·
zi b1))2[0100]
q2＝σ(y
i-(k2·
zi b2))2[0101]
其中zi表示切片s
1，i
在人工标注结果m1中所对应的z轴上的坐标值。xi、yi分别表示切片s
1，i
几何中心p
c，i
的x轴和y轴上所对应的坐标值；
[0102]
当q1、q2最小，即残差平方和最小时，即可得到常数k1，k2，b1，b2的最佳值，因为q1、q2是关于k1，k2，b1，b2的函数，当q1、q2的导数为0时，可取到相应的最小值，所以，分别对k1，k2，b1，b2求偏导可以得到如下公式：
[0103][0104]
从而得到k1，k2，b1，b2的表达式如下：
[0105][0106][0107]
遍历几何中心集p2中所有几何中心，带入k1，k2，b1，b2的表达式中，求得常数k1，k2，b1，b2的最佳值，最后即可得到拟合的公式(3)的空间直线表达式，拟合得到的空间直线即为m1的中心轴l
center，i
；
[0108]
步骤1.7、将训练集d
train1
中的所有人工标注结果进行步骤1.1～
[0109]
步骤1.6的操作，最后得到一个中心轴数据集l
center
(l
center，i
∈l
center
，i∈[0，d
1-1])；
[0110]
步骤s2、将训练集d
train1
和中心轴集l
center
输入至旋转模块中，得到处理后的数据集d
train2
。
[0111]
如图2所示，本技术所述将训练集d
train1
和中心轴集l
center
输入至旋转模块，过程如下：
[0112]
步骤2.1、根据中心轴计算模块得到的中心轴l
center，i
(l
center，i
∈l
center
，i＝0，1，2，...)在大小为d1×
h1×
w1的人工标注结果m1中的空间位置，计算中心轴l
center，i
与z轴的夹角。根据夹角，将髁突cbct影像i1以及人工标注结果m1(i1∈d
train1
，m1∈d
train1
)进行旋转，使得髁突的中心轴l
center，i
旋转至z轴，并进行重采样操作，最后得到大小仍为d1×
h1×
w1的新的髁突cbct影像i2以及人工标注结果m2；
[0113]
步骤2.2、将训练集d
train1
中的所有髁突cbct影像以及人工标注结果进行步骤2.1的操作，最后得到处理后的数据d
train2
(i2∈d
train2
，m2∈d
train2
)；
[0114]
步骤s3、将d
train2
输入至长轴计算模块中，得到d
train2
中所有的长轴，组成一个长轴集l
long
。
[0115]
如图2所示，本技术所述长轴计算模块，过程如下：
[0116]
步骤3.1、给定一个阈值t1，将人工标注结果m2(m2∈d
train2
)的z轴坐标值小于t1的部分舍弃，从而得到仅含髁突头的感兴趣区域i
hroi
；
[0117]
步骤3.2、将髁突头的感兴趣区域i
hroi
沿着z轴方向得到d1张大小为d1×
h1×
w1的髁突头的感兴趣区域的切片，得到切片集s3；
[0118]
步骤3.3、将切片集s3中每张切片进行图像灰度化和二值化处理，得到d1张二值化后的图像，组成新的切片集s4；
[0119]
步骤3.4、根据切片集s4，借助cv2库中的轮廓检测算法，得到所有切片中髁突头的感兴趣区域的轮廓的点坐标集p3；
[0120]
步骤3.5、根据每张髁突头的感兴趣区域的切片s
2，i
(s
2，i
∈s4，i∈[0，d
1-1])中的轮廓点坐标集p
2，i
(p
2，i
∈p3，i∈[0，d
1-1])，记录轮廓中线段最长的两点的坐标，即记录髁突长轴的坐标，最后得到d1个线段组成的线段集l；
[0121]
步骤3.6、根据得到的线段集l拟合一条平面直线，其中平面直线方程的一般表达式如下：
[0122]
y＝k3·
x b3(4)
[0123]
根据最小二乘法原理，得到常数k3，b3的最佳值，其中最小二乘法拟合原理如下：拟合公式(4)的平面直线表达式，并使得残差平方和最小，其中残差平方和的公式如下：
[0124]
q3＝σ(k3·
xi b
3-yi)2[0125]
其中xi、yi分别表示切片s
2，i
上髁突长轴的x轴和y轴的坐标值；
[0126]
当q3最小，即残差平方和最小时，即可得到常数k3，b3的最佳值，因为q3是关于k3，b3的函数，当q3的导数为0时，可取到相应的最小值，所以，分别对k3，b3求偏导可以得到如下公式：
[0127][0128]
从而得到k3，b3的表达式如下：
[0129][0130]
其中表示的算术平均值，计算公式如下：
[0131][0132]
遍历线段集l中所有长轴，带入k3，b3的表达式中，求得常数k3，b3的最佳值，最后即可得到拟合的公式(4)的平面直线表达式，拟合得到的平面直线即为m2的髁突头的长轴；
[0133]
步骤3.7、将训练集d
train2
中的所有人工标注结果进行步骤3.1～
[0134]
步骤3.6的操作，最后得到d
train2
中髁突头的长轴数据集l
long
；
[0135]
步骤s4、将训练集d
train2
和长轴集l
long
输入至旋转模块中，得到最终用于训练髁突分割模型的数据集d
train3
。
[0136]
如图2所示，本技术所述训练集d
train2
和长轴集l
long
输入至旋转模块，过程如下：
[0137]
步骤4.1、根据长轴计算模块得到的髁突头长轴li(li∈l
long
，i＝0，1，2，...)在大小为d1×
h1×
w1的人工标注结果m2中的空间位置，计算长轴与x轴的夹角，根据夹角，将髁突cbct影像i2以及人工标注结果m2(i2∈d
train1
，m2∈d
train1
)进行旋转，使得髁突的长轴li旋转至x轴，并进行重采样操作，最后得到大小仍为d1×
h1×
w1的新的髁突cbct影像i3以及人工标注结果m3；
[0138]
步骤4.2、将训练集d
train2
中的所有髁突cbct影像以及人工标注结果进行步骤4.1的操作，最后得到处理后的数据d
train3
(i3∈d
train3
，m3∈d
train3
)。
[0139]
需要说明的是，本技术的髁突原始数据与进行数据增强后的数据在相同空间位置上，即x，y，z轴坐标相同情况下的拟样例对比图如图3所示，其中(a)为髁突原始数据，(b)为经过本技术提出的方法得到的影像。如(b)所示，在空间位置相同的情况下，(b)中各个轴面的髁突信息更为丰富，从而有利于提高髁突分割模型的精度。
[0140]
需要说明的是，本技术中d为图片的深度，h为图片的高度，w为图片的宽度，字母的下标表示序号，用以区别不同特征图的维度。
[0141]
本技术通过中心轴计算模块，计算得到了髁突的中心轴。通过图像旋转处理模块，结合髁突的中心轴，找到了髁突的最佳轴面。把处理后的数据扩充到训练样本集中，从而提升利用数据增强后髁突的cbct影像训练样本集进行髁突分割的精度以及模型的泛化能力。
[0142]
本说明书的实施例所述的内容仅仅是对发明构思的实现形式的列举，仅作说明用途。本发明的保护范围不应当被视为仅限于本实施例所陈述的具体形式，本发明的保护范围也及于本领域的普通技术人员根据本发明构思所能想到的等同技术手段。