CN102663387A - 基于牙科全景图的皮质骨宽度自动计算方法 - Google Patents

Abstract

本发明提供一种基于牙科全景图片的皮质骨宽度CW计算方法，包括如下步骤：1)基于OTSU和DLS的颏孔识别：a）图像预处理，使图像的颏孔和轮廓变得明显；b）采用OTSU方法对图像的阈值分割；c）形态学处理，让阈值分割后的图像边缘更加平滑；d）利用Canny算子检测图像边缘；e）最小二乘法椭圆拟合；f）结束；2)基于图像聚类和DLS的皮质骨宽度计算；a）提取ROI，提取左ROI时是以左颏孔为右上角，以颏孔到全景图下边缘为边长；而右ROI是以颏孔为左上角，以颏孔到全景图下边缘为边长；b)ROI预处理；c)图像聚类；d)形态学处理；e)最小二乘法曲线拟合；f计算皮质骨宽度，即计算颏孔点到上下边缘的距离差；g)结束。本发明在基于牙科全景图片下，实现了颏孔的自动识别和基于ROI的CW计算，在CW计算方法中提高了CW计算的准确性和效率，为医学图像处理提供了方便。

Description

基于牙科全景图的皮质骨宽度自动计算方法

技术领域

本发明涉及一种图像处理领域，是一种基于牙科全景图的皮质骨宽度自动计算方法

背景技术

基于牙科全景图的皮质骨宽度计算方法是在骨质疏松检测的基础上提出来的。

骨质疏松是一种妇女和老年人经常面临的一种疾病，由于其危害极大，而且不易察觉，因此对其预测显得尤为重要。在近几年的研究中多种研究成果显示骨质疏松疾病与皮质骨宽度值有着很密切的关系，但是皮质骨宽度的测量如果每次都是手工测量不仅费时而且可能会不准确，因此利用计算机对于皮质骨宽度CW的自动测量就显得很有必要。

在对于CW计算中，技术难题主要就是图像的分割。前人对图像分割中有了一些比较成熟的技术，包括了基于阈值的分割、基于边缘检测的图像分割、基于区域特性的图像分割、基于模式识别的图像分割等等。其中基于阈值的分割中OTSU(基于类间方差最大法)有着不错的口碑。基于边缘检测中有很多的算子比如说Canny，Sobel，Laplacian等。基于区域特性的图像分割可分为两种区域生长和分裂合并，分水岭算法算是里面比较新的一种分割方法。基于模式识别的分割技术中监督算法有神经网络，支持向量机等。无监督算法又叫做聚类算法。之前也有人做过CW计算方法的研究，但是对于颏孔的手动标识和在对上下边缘进行检测时需要用户给出有效的像素点使得方法存在着不足，使得用户的体验显得很不足，用户希望是能够直接得到CW的计算结果，越方便越好。

发明内容

本发明所要解决的技术问题是自动计算基于牙科全景图的皮质骨宽度(CW)值，主要关注其自动性、准确性及其效率。

为解决上述问题，本发明基于牙科全景图的皮质骨宽度(CW)自动计算方法包括如下步骤：

1)基于OTSU和和DLS(Direct Least Square)的颏孔识别：

a图像预处理，使图像的颏孔和轮廓变得明显；

b采用OTSU方法对图像的阈值分割；

c形态学处理，让阈值分割后的图像边缘更加平滑；

d利用Canny算子检测图像边缘；

e最小二乘法椭圆拟合，将颏孔的外观定义成椭圆，然后对其边缘进行拟合；

f结束；

2)基于聚类和DLS的皮质骨宽度计算：

a颏孔识别提取ROI(Region of Interest)，提取左ROI时是以左颏孔为右上角，以颏孔到全景图下边缘为边长；而右ROI是以颏孔为左上角，以颏孔到全景图下边缘为边长；

b ROI预处理；

c图像聚类，将图像分成两类，将各个灰度值看成一类，然后将离得较近的、相邻的类合并直至变成两类为止；

d形态学处理；

e最小二乘法曲线拟合，将图像模拟成为多项式曲线；

f计算皮质骨宽度，即从颏孔处向两条边缘线做垂线，两个距离之差即为皮质骨宽度；

g结束。

本发明中，步骤1)-a中所述的图像预处理是先将24位图变成8位的，然后对其进行灰度拉伸，灰度拉伸后进行图像的锐化。

本发明中，步骤1)-b中进行OTSU图像分割方法如下：

计算图像各个灰度值的概率；

用T将灰度值分成两组C₀＝{0～t}和C₁＝{T+1～m}，计算组间方差；

将T进行变化(0～m)，分别计算阈值选择函数δ²(T)＝w₀(μ₀-μ)²+w₁(μ₁-μ)²，求出最大值，就选择其为阈值点进行返回。

本发明中，步骤1)-c中进行形态学处理是指进行形态学腐蚀膨胀处理.

本发明中，步骤2)-b中ROI图像预处理对图像进行灰度化、锐化、去除前景和背景操作。

本发明中，步骤2)-c中图像聚类的步骤如下：

令C_k为灰度图中升序排第k类的一个集群，T_k是C_k类中的最大的灰度值，因此C_k包含了从[T_k-1+1，T_k]的灰度值，这里我们定义T₀＝-1，接下来是一个归并的过程，如下：

1)假设直方图包含了k个不同的非空的灰度值，在开始的阶段，每个类中只包含一个灰度值；

2)下面的过程将会持续k-t次，来将k类变成t类：

2.1)每个相邻类之间的距离将会被计算，这个距离就表明了相邻类的相异程度；

2.2)找出相邻距离最短的路径，并且将这两个类合并成一个类。由于类会越变越小，C和T将会被重新赋值；

3)最后将会获得t个类。

本发明中，每个相邻类之间的距离的计算方法如下：令h(z)表示目标图像的直方图，z＝0，1，...L-1，z表示的是灰度等级，而L为最大可以表示的最大等级个数，定义p(z)＝h(z)/N，N表示的是图像中像素的个数，p(z)表示的是灰度等级为z的像素的概率，如下所示：

$P\left(C_k\right)=\underset{z=T_{k-1}+1}{\overset{T_k}{\mathrm{\Σ}}}p\left(z\right),\underset{k=1}{\overset{K}{\mathrm{\Σ}}}P\left(C_k\right)=1.$

而Ck1和Ck2之间的距离定义如下：

$\mathrm{Dist}(C_{k1},C_{k2})={\mathrm{\δ}}_1^2(C_{k1}\∪C_{k2}){\mathrm{\δ}}_A^2(C_{k1}\∪C_{k2}).$

而 ${\mathrm{\δ}}_1(C_{k1}\∪C_{k2})=\frac{P\left(C_{k1}\right)}{P\left(C_{k1}\right)+P\left(C_{k2}\right)}{[m\left(C_{k1}\right)-M(C_{k1}\∪C_{k2})]}^2$

$+\frac{P\left(C_{k2}\right)}{P\left(C_{k1}\right)+P\left(C_{k2}\right)}{[m\left(C_{k2}\right)-M(C_{k1}\∪C_{k2})]}^2+$

$\frac{P\left(C_{k1}\right)P\left(C_{k2}\right)}{{(P\left(C_{k1}\right)+P\left(C_{k2}\right))}^2}{[m\left(C_{k1}\right)-m\left(C_{k2}\right)]}^2$

在公式里m(C_k)表示的是C_k的平均值如下表示：

$m\left(C_k\right)=\frac{1}{P\left(C_k\right)}\underset{z=T_{k-1}+1}{\overset{T_k}{\mathrm{\Σ}}}\mathrm{zp}\left(z\right).$

而M(C_k1∪C_k2)则表示是类C_k1与C_k2的平均值，如下表示：

$M(C_{k1}\∪C_{k2})=\frac{P\left(C_{k1}\right)m\left(C_{k1}\right)+P\left(C_{k2}\right)m\left(C_{k2}\right)}{P\left(C_{k1}\right)+P\left(C_{k2}\right)}$

内部类方差

表示如下所示：

${\mathrm{\δ}}_A^2(C_{k1}\∪C_{k2})=\frac{1}{P\left(C_{k1}\right)+P\left(C_{k2}\right)}*\underset{z=T_{k1-1}+1}{\overset{T_{k2}}{\mathrm{\Σ}}}\left[{(z-M(C_{k1}\∪C_{k2}))}^{2}p\left(z\right)\right].$

本发明中，步骤2)-f)中图像CW值计算是从颏孔处对已经拟合出来的两条抛物线做垂线，两者距离差就是CW值。

本发明的有益效果是：本发明的基于牙科全景图片的CW计算方法在保持计算精度的准确性的前提下，在自动识别方面做出了改进，对于类算法也作出了改进，提高了计算的准确性和好的用户体验。

附图说明

图1是本发明一实施例的基于牙科全景图的皮质骨宽度自动计算方法的流程图。

图2为图1中颏孔识别的流程图。

图3为图1中CW计算的流程图。

图4为图3中图像聚类步骤的具体流程图。

具体实施方式

为了更了解本发明的技术内容，特举具体实施例并配合所附图式说明如下。

图1是本发明一实施例的基于牙科全景图的皮质骨宽度自动计算方法的流程图。基于牙科全景图的皮质骨宽度自动计算方法包括两个阶段：基于OTSU和DLS的颏孔识别、基于聚类和DLS的皮质骨宽度计算。

图2为本发明基于牙科全景图的皮质骨宽度(CW)自动计算方法的颏孔识别的流程图。如图2所示，在颏孔识别流程图中：

1)图像预处理中先将24位图变成8位的，然后对其进行灰度拉伸，灰度拉伸主要是针对一些颏孔不是很清晰的图片，拉伸灰度后会使得颏孔更加明显。图像锐化也是为了增强对比度，使得图像的轮廓更加明显，便于后面的阈值分割。

2)图像的阈值分割是采用OTSU方法，对于计算出来的阈值T进行如下操作转换成二值图像。

$g(x,y)=\left\{\begin{array}{cc}0& f(x,y)\≤T\\ 255& f(x,y)>T\end{array}\right.$

3)形态学的处理主要是针对图像的噪声，让阈值分割后的图像边缘更加平滑

4)Canny算子是通过寻找图像梯度的局部极大值。它用两个阈值分别来检测强边缘和弱边缘，而且仅仅当弱边缘和强边缘相连时，弱边缘才会包含在输出中。因此，Canny边缘算子不容易受到噪声的干扰，这也是本发明选择Canny算子的原因。

5)在进行椭圆拟合时，所用的数据便是边缘检测所检测的边缘像素。由于所看见的颏孔是近似椭圆的形状，因此将颏孔的外观定义成椭圆，然后对其边缘进行拟合。其实拟合的结果可能不止一个，因为除了颏孔处还会有其他区域会有一块比其他地方灰度值低，造成边缘，我们可以通过颏孔的特殊性将其他点进行排除，颏孔的位置是出于第四颗与第五颗牙齿下方，颏孔的大小虽不固定，但是有个范围)。

图3是图1中基于图像聚类和DLS的皮质骨宽度计算(CW)的流程图。如图3所示，在CW计算流程图中：

1.提取左ROI时是以左颏孔为右上角，以颏孔到全景图下边缘为边长；而右ROI是以颏孔为左上角，以颏孔到全景图下边缘为边长。

2.ROI预处理中要对图像进行灰度化，锐化和去除前景和背景等操作。图像的锐化是为了突出边缘，增强对比度。而去除前景和背景则是为了减少客观原因对最后产生的误差。

3.图像聚类是将图像分成两类。主要思想是将各个灰度值先看成一类，然后将离得较近的相邻的类合并直至变成两类为止。其中聚类的具体操作过程如图4所示，后文会对图4进行具体描述。

4.这里的曲线模拟是将图像模拟成为多项式曲线，如

形式。

5.图像CW值计算是从颏孔处对已经拟合出来的两条抛物线做垂线，两者距离差就是CW值。假设二次曲线的方程组为y＝a+bx+cx²，而颏孔的坐标为(x₀，y₀)，则可以用下面方程组求出垂线与曲线的交点(切线与垂线垂直，点在曲线上)。

$\left\{\begin{array}{c}y=a+\mathrm{bx}+{\mathrm{cx}}^2\\ \frac{y-y_0}{x-x_0}(2\mathrm{cx}+b)=-1\end{array}\right.$

根据交点值与颏孔值确定了这个直线，直线与两个抛物线之间的距离差就是即为CW值。

图4是在CW计算中图像聚类步骤的具体流程图。如图4所示，图像聚类的过程如下：

1)开始我们先初始化将图像分成k个类(有k个不同的灰度值)，即类的个数就是不同灰度值的个数，每个灰度值就是一个类。初始化C，T，C_k为灰度图中升序排第k类的一个集群，T_k是C_k类中的最大的灰度值，因此C_k包含了从[T_k-1+1，T_k]的灰度值。初始化Num＝0，Num表示的是循环的次数。转2)

2)判断Num是否小于k-2，若小于转3)；若不小于转7)

3)计算相邻两个类之间的距离，这个距离就表明了相邻类的相异程度，这个距离的计算将会在下面定义，两者的距离越近表示越相似。转4)

4)重复步骤3)k-Num次，然后跳5)

5)找出k-Num个距离中的最小值为第M个类，转6)

6)将第M和第M+1个类合并，由于类会越变越小，C和T将会被重新赋值，Num++转2)

7)返回T1值，转8)

8)结束

每个相邻类之间的距离的计算方法如下：我们令h(z)表示目标图像的直方图，z＝0，1，...L-1，z表示的是灰度等级，而L为最大可以表示的最大等级个数。我们定义p(z)＝h(z)/N(这里的h(z)是灰度值为z的像素的个数)，这里N表示的是图像中像素的个数，p(z)表示的是灰度等级为z的像素的概率，如下所示：

$P\left(C_k\right)=\underset{z=T_{k-1}+1}{\overset{T_k}{\mathrm{\Σ}}}p\left(z\right),\underset{k=1}{\overset{K}{\mathrm{\Σ}}}P\left(C_k\right)=1.$

(1)

而C_k1和C_k2之间的距离定义如下：

$\mathrm{Dist}(C_{k1},C_{k2})={\mathrm{\δ}}_1^2(C_{k1}\∪C_{k2}){\mathrm{\δ}}_A^2(C_{k1}\∪C_{k2}).$

(2)

而 ${\mathrm{\δ}}_1(C_{k1}\∪C_{k2})=\frac{P\left(C_{k1}\right)}{P\left(C_{k1}\right)+P\left(C_{k2}\right)}{[m\left(C_{k1}\right)-M(C_{k1}\∪C_{k2})]}^2$

$+\frac{P\left(C_{k2}\right)}{P\left(C_{k1}\right)+P\left(C_{k2}\right)}{[m\left(C_{k2}\right)-M(C_{k1}\∪C_{k2})]}^2+$

$\frac{P\left(C_{k1}\right)P\left(C_{k2}\right)}{{(P\left(C_{k1}\right)+P\left(C_{k2}\right))}^2}{[m\left(C_{k1}\right)-m\left(C_{k2}\right)]}^2---\left(3\right)$

在公式里m(C_k)表示的是C_k的平均值如下表示：

$m\left(C_k\right)=\frac{1}{P\left(C_k\right)}\underset{z=T_{k-1}+1}{\overset{T_k}{\mathrm{\Σ}}}\mathrm{zp}\left(z\right).$

(4)

而M(C_k1∪C_k2)则表示是类C_k1与C_k2的平均值，如下表示：

$M(C_{k1}\∪C_{k2})=\frac{P\left(C_{k1}\right)m\left(C_{k1}\right)+P\left(C_{k2}\right)m\left(C_{k2}\right)}{P\left(C_{k1}\right)+P\left(C_{k2}\right)}$

(5)

内部类方差

表示如下所示：

${\mathrm{\δ}}_A^2(C_{k1}\∪C_{k2})=\frac{1}{P\left(C_{k1}\right)+P\left(C_{k2}\right)}*\underset{z=T_{k1-1}+1}{\overset{T_{k2}}{\mathrm{\Σ}}}\left[{(z-M(C_{k1}\∪C_{k2}))}^{2}p\left(z\right)\right].$

(6)

根据如上的6个公式可以算出相邻两个类之间的距离。

本发明的基于牙科全景图片的CW计算方法在保持计算精度的准确性的前提下，在自动识别方面做出了改进，对于类算法也作出了改进，提高了计算的准确性和好的用户体验。

虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因此，本发明的保护范围当视权利要求书所界定者为准。

Claims (7)

1.一种基于牙科全景图片的皮质骨宽度计算方法，其特征在于，包括如下步骤：

1)基于OTSU和DLS的颏孔识别：

a图像预处理，使图像的颏孔和轮廓变得明显；

b采用OTSU方法对图像的阈值分割；

c形态学处理，让阈值分割后的图像边缘更加平滑；

d利用Canny算子检测图像边缘；

e最小二乘法椭圆拟合，将颏孔的外观定义成椭圆，然后对其边缘进行拟合；

f结束；

2)基于图像聚类和DLS的皮质骨宽度计算：

a颏孔识别提取ROI，提取左ROI时是以左颏孔为右上角，以颏孔到全景图下边缘为边长；而右ROI是以颏孔为左上角，以颏孔到全景图下边缘为边长；

b ROI预处理；

c图像聚类，将图像分成两类，将各个灰度值看成一类，然后将离得较近的、相邻的类合并直至变成两类为止；

d形态学处理；

e最小二乘法曲线拟合，将图像模拟成为多项式曲线；

f计算皮质骨宽度，即从颏孔处向两条边缘线做垂线，两个距离之差即为皮质骨宽度； 

g结束。

2.根据权利要求1所述的基于牙科全景图片的皮质骨宽度计算方法，其特征在于，其中步骤1)-a中所说的图像预处理是先将24位图变成8位的，然后对其进行灰度拉伸，灰度拉伸后进行图像的锐化。

3.根据权利要求1所述的基于牙科全景图片的皮质骨宽度计算方法，其特征在于，其中步骤1)-b中进行OTSU图像分割方法如下

计算图像各个灰度值的概率；

用T将灰度值分成两组C₀＝{0～T}和C₁＝{T+1～m}，计算组间方差；

将T进行变化(0～m)，分别计算阈值选择函数δ²(T)＝w₀(μ₀-μ)²+w₁(μ₁-μ)²，求出最大值，就选择其为阈值点进行返回。

4.根据权利要求1所述的基于牙科全景图片的皮质骨宽度计算方法，其特征在于，其中步骤1)-c中进行形态学处理是指进行形态学腐蚀膨胀处理。.

5.根据权利要求1所述的基于牙科全景图片的皮质骨宽度计算方法，其特征在于，其中步骤2)-b中ROI图像预处理对图像进行灰度化、锐化、去除前景和背景操作。

6.根据权利要求1所述的基于牙科全景图片的皮质骨宽度计算方法，其特征在于，其中步骤2)-c中图像聚类的步骤如下：令C_k为灰度图中升序排第k类的一个集群，T_k是C_k类中的最大的灰度值，因此C_k包含了从[T_k-1+1，T_k]的灰度值，定义T₀＝-1，

1)假设直方图包含了k个不同的非空的灰度值，在开始的阶段，每个类中只包含一个灰度值； 

2)下面的过程将会持续k-t次，来将k类变成t类：

2.1)每个相邻类之间的距离将会被计算，这个距离就表明了相邻类的相异程度；

2.2)找出相邻距离最短的路径，并且将这两个类合并成一个类。由于类会越变越小，C和T将会被重新赋值；

3)最后将会获得t个类。

7.根据权利要求6所述的基于牙科全景图片的皮质骨宽度计算方法，其特征在于，其中步骤2.1)中每个相邻类之间的距离的计算方法如下：

令h(z)表示目标图像的直方图，z＝0，1，...L-1，z表示的是灰度等级，而L为最大可以表示的最大等级个数，定义p(z)＝h(z)/N，N表示的是图像中像素的个数，p(z)表示的是灰度等级为z的像素的概率，如下所示：

而Ck1和Ck2之间的距离定义如下：

而

在公式里m(C_k)表示的是C_k的平均值如下表示： 

而M(C_k1∪C_k2)则表示是类C_k1与C_k2的平均值，如下表示：

内部类方差 

表示如下所示：

。

8.根据权利要求1所述的基于牙科全景图片的皮质骨宽度计算方法，其特征在于，其中步骤2)-f)中图像CW值计算是从颏孔处对已经拟合出来的两条抛物线做垂线，两者距离差就是CW值。 

Images