CN104318563A - 一种基于医学图像的器官骨架提取方法 - Google Patents

Abstract

本发明提供一种基于医学图像的器官骨架提取方法，在利用DICOM格式医学图像的基础上，以提取一个完整的器官骨架为目标，围绕基于结构张量的局部结构提取、基于K-means模糊聚类的点云模型分割、基于点云的骨架提取方法展开研究，首先从DICOM格式医学图像中读取三维数据，并根据三维数据的局部结构特征对三维数据进行特征提取得到点云，然后对得到的点云通过K-means聚类进行模型分割，最终对得到的分割点云进行骨架提取。

Description

一种基于医学图像的器官骨架提取方法

技术领域

本发明涉及一种基于医学图像的器官骨架提取方法。

背景技术

随着计算机在各行业的广泛应用，人们开始不满足于仅能显示二维的图像，更希望计算机能表现出真实感强三维世界。计算机图形学、虚拟现实等技术的发展，三维数据成为继图像、音频、视频以来又一重要的多媒体数据形式。三维数据就是使用各种采集仪采集得到的，它记录了有限体表面在离散点上的各种物理参量。它包括最基本信息是物体的各离散点三维坐标，其它的可以包括物体表面颜色、透明度、纹理特征等。由于三维模型的信息量很大，这也使得在许多应用中出现占用内存空间过大、运行计算负载过重，所以需要一种“紧凑的”方式来尽可能完整、全面地表示描述三维模型的结构特征信息。其中最常用一种简化表示方式就是使用一维曲线，一般称为骨架。利用物体的骨架来描述对象是一种既能强调物体的结构特征，也能提高内存使用效率与数据压缩的好方法。

为了找到“紧凑的”一维曲线(骨架)，本发明提出了一种基于医学图像的器官骨架提取方法，以提取完整的器官骨架结构为目标，围绕三维数据局部结构、模糊分类、骨架提取展开研究。首先通过对医学文件读取得到CT的原始扫描数据，并对三维数据计算其二阶黑塞矩阵的特征值，选取特征值满足一定条件的三维数据，然后结合K-means模糊聚类，对得到的三维数据进行分类，最后通过L1-中值定义和最小化能量泛函，并加入正则化得到一个完整的器官骨架结构。

发明内容

本发明解决的技术问题是：提供了一种基于医学图像的器官骨架提取方法，且可以直接对原始扫描数据进行骨架提取，适用度广。

本发明采用的技术方案为：一种基于医学图像的器官骨架提取方法，主要包括以下三个步骤：

步骤(1)、读取三维医学图像数据：通过对DICOM格式的医学CT数据读取，得到医学器官的三维数据，作为整个算法原始输入数据。

步骤(2)、提取三维数据的局部特征结构：通过对得到三维数据固有的结构特征进行分析来提取局部结构；利用三维数据的黑塞矩阵来刻画局部几何结构，并通过计黑塞矩阵的特征值进行组合来提取合适的局部几何结构，最大特征值所对应的特征向量表征了从一种材质过度到另一种材质时最为迅速的变化方向，而最小特征值所对应的方向则指明了某种材质在图像内部流动的方向；因此可以对三维图像内部的材质属性变化进行定量刻画；得到的三维离散数据(也就是点云模型)作为整个骨架提取的输入数据；

步骤(3)、点云模型分割：基于步骤(1)中的数据对数据进一步进行优化，利用K-means对步骤(1)中数据进行分类，得到一个基于灰度值的分类结果数据，以此作为骨架提取的最终输入数据；

步骤(4)、骨架结构提取：步骤(3)中得到的分类结果数据作为骨架提取的输入数据，基于步骤(3)中的分类结果数据，利用L1-中值理论，并定义一个能量泛函，通过求解能量最小的中心点来代表骨架分支；通过结合正则化，得到一个一维连续的骨架结构；

以上的所有操作都基于交互的模式，用户通过调整参数，提取合适的局部特征结构，而且对于数据中含有噪声和孤立点，甚至对于数据部分缺失都可以提取一个完整的骨架借结构。

以上所有操作都是用Qt和OpenGL来实现的，具有很好的交互性，用户可以通过选定合适参数对提取不同形状，具有很高的适用性。

本发明的原理在于：

(1)通过微分几何理论来对图像特征进行分析，源自图像微分分析的结构张量可紧凑的刻画图像局部材质变化的方向以及材质间的过度变化关系，所以可以利用黑塞矩阵来表达图像特征，通过选取黑塞矩阵特征值满足一定条件的点，可以实现对三维数据的局部几何结构进行提取。

(2)为了解决复杂多分支模型情况下对骨架提取的影响，本发明利用K-means聚类方法，对从(1)中获取的数据进行聚类，该方法可以对点云模型进行有意义的分割。

(3)为了对点云模型进行骨架提取，本发明依据L1-中值理论，定义一个能量泛函，并求解使能量泛函在局部领域中最小值的点作为骨架节点，通过加入正则化，最终得到一个完整的一维连续的骨架结构。

本发明与现有技术相比的优点在于：

1、本发明可以直接使用原始CT扫描数据，对于复杂的医学CT扫描数据可以直接应用而不用进行噪声等的预处理。本发明对于含有噪声、孤立点，甚至是数据部分缺失都能提取完整的骨架结构。

2、良好的交互性，本发明基于Qt来进行交互框架的开发，用户可以通过手动调整参数来提取图像的不同局部结构，适用度高。

附图说明

图1为基于医学图像的器官骨架提取方法流程图；

图2为正方体盒子点云模型图；

图3为动脉血管的点云模型图；

图4k＝3正方体盒子聚类结果图；

图5k＝4正方体盒子聚类结果图；

图6k＝3动脉血管聚类结果图；

图7k＝4动脉血管聚类结果图；

图8正方体盒子点云模型的骨架提取结果图；

图9动脉血管点云模型的骨架提取结果图。

具体实施方式

图1给出了基于医学图像的器官骨架提取方法的总体处理流程图，下面结合其他附图及具体实施方式进一步说明本发明。

本发明提供一种基于医学图像的器官骨架提取方法，主要步骤介绍如下

1、提取三维数据局部特征结构方法

本发明采用微分几何理论来对图像特征进行分析，利用图像微分分析的结构张量可紧凑的刻画图像局部材质变化的方向以及材质间的过渡变化关系，由图像的二阶导数构成的具有对称性质的黑塞矩阵H能够很好地表示曲面法线的变化情况，进而可以用图像的二阶导数构成的黑塞矩阵来刻画三维数据的局部几何结构。

本发明重点在于利用黑塞矩阵的不同组合来提取合适的局部几何结构，从某种意义上说，最大特征值所对应的特征向量表征了从一种材质过度到另一种材质时最为迅速的变化方向，而最小特征值所对应的方向则指明了某种材质在图像内部流动的方向。并且，最大和最小特征值恰恰蕴含了图像的主曲率信息。因此，从各向异性热传递的角度来说，同时考虑黑塞矩阵的特征值和特征向量无疑可以克服各向同性滤波器的一些局限性。同时，提取三维图像的局部结构并对三维图像内部的材质属性变化进行定量刻画也正是我们的目的所在。

通过选择黑塞矩阵H特征值满足不同条件的三维数据来提取不同局部几何结构，通过将黑塞矩阵的特征值进行不同形式的组合运算，分析了它们与图像局部结构(平面、块和噪声)之间的关系：

R_sheet＝(|λ₃|+|λ₂|)/|λ₁|

R_blob＝(2|λ₁|-|λ₂|-|λ₃|)/2|λ₁|

R_tube＝(|λ₁|-|λ₂|-|λ₃|)/|λ₁|

$R_{\mathrm{noise}}=\sqrt{{\mathrm{\λ}}_1^2+{\mathrm{\λ}}_2^2+{\mathrm{\λ}}_3^2}$

其中λ₁、λ₂、λ₃是黑塞矩阵H的特征值，且|λ₁|≥|λ₂|≥|λ₃|。上述四个公式分别代表平面、快、管状、噪声四中结构，当它们所代表的结构出现时，上述4种特征值的比率值会相应的变为0。据此，我们可以通过将特征向量、特征值和上述描述局部几何结构的比率值进行组合，进而我可以通过几何结构比率值进行判定局部几何结构。(如图2,3所示的为提取正方体盒子和动脉血管的管状点云模型)。

2、基于K-means点云模型分割

本文采用K-means聚类方法，将1中得到的点云模型进行分割，对于复杂的多分支模型,在骨架提取之前，为了减少各分支间的相互干扰，需对点云模型进行有意义的分割。

本发明采用K-means流程如下：

输入：聚类个数k，以及包含n个三维数据的体素数据。

输出：满足方差最小标准的k个聚类。

(1)首先选取聚类个数K，在我们实现中可以通过交互给出K的值，这样就可以合理的选取k值。

(2)在n个输入三维数据中选出K个点的体素作为初始的中心点。

(3)对n个输入三维数据每个数据都和(2)中K个初始中心点做比较，例如每个数点与初始中心点i最小，就把这个三维数据点归为i类。

(4)所有点都比较完后，对每类从新计算其中心点。将每个类中点的灰度值相加在除以每个类中点的个数作为新的中心点。

(5)将所得的中心点与旧中心点比较，如果小于一个阀值，就认为算法结束，否则经行(3)。(如图3,4为k＝3和k＝4正方体盒子点云模型的聚类结果，不同颜色代表不同的类，如图5,6为k＝3和k＝4动脉血管点云模型的聚类结果)。

3、基于点云的骨架提取

本发明根据L1-中值理论，定义一个能量泛函，通过求解使能量泛函满足局部最小值的中心点来代表点云骨架，并加入正则化理论，最终得到一个完整的一维连续的骨架曲线。

a)L1-中值点集的计算

本发明就是利用L1-中值从原始输入点Q＝{q_j}_j∈J来找出一组局部中心点集合，局部中值点集合可以定义为：

$\underset{X}{\arg \min }\underset{i\∈I}{\mathrm{\Σ}}\underset{j\∈J}{\mathrm{\Σ}}\left|\right|x_i-q_j\left|\right|\mathrm{\θ}\left(\right||x_i-q_j|\left|\right)+R\left(X\right)---\left(1\right)$

公式(1)可以看作局部版本的L1-中值，其中X＝{x_i}_i∈I为要寻找的局部中值点的集合。与局部相关的是函数是一个快速衰减的平滑函数，与局部相关的变量h代表局部邻域的半径。R(X)用来控制局部中心点的分布，后面会详细介绍。公式(1)就是通过L1-中值求出一组局部中心点集合到原始输入点集合的距离和最小。

b)正则化

正则化理论首先由统计学家提出，试图用模型来拟合严重欠约束解空间的数据。正则化是指在线性代数理论中，不适定问题通常是由一组线性代数方程定义的，而且这组方程组通常来源于有着很大的条件数的不适定问题。由于单独应用L1-中值，会造成中心点会聚积到一起从而产生一个稀疏的分布，所以需要一个“惩罚项”来控制中心点的分布。所以R(X)的公式定义为：

$R\left(X\right)=\underset{i\∈I}{\mathrm{\Σ}}{\mathrm{\γ}}_i\underset{i^{\′}\∈I/i}{\mathrm{\Σ}}\frac{\mathrm{\θ}\left(\right||x_i-q_j|\left|\right)}{{\mathrm{\δ}}_i\left|\right|x_i-x_{i^{\′}}\left|\right|}$

其中{γ_i}_i∈I是中心点集的平衡常数，δ_i被称作方向度数，用来指导骨架的行程过程。现在已经确定了R(X)的函数形式，剩下就是对于公式(1)的求解。对于每一个中心点为它所在邻域内使公式(1)最小的点。也就是每一个中心点为公式(1)在其邻域内极小值，也就是公式(1)在每个中心点的梯度为0。求得满足上述条件的点为中心点，也就是用其来代表骨架(如图8、9分别为正方体盒子的骨架结构图和动脉血管骨架结构图)。

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (5)

1.一种基于医学图像的器官骨架提取方法，其特征在于包括以下三个步骤：

步骤(1)、读取三维医学图像数据：通过对DICOM格式的医学CT数据读取，得到医学器官的三维数据，作为整个算法原始输入数据。

步骤(2)、提取三维数据的局部特征结构：通过对得到三维数据固有的结构特征进行分析来提取局部结构；利用三维数据的黑塞矩阵来刻画局部几何结构，并通过计黑塞矩阵的特征值进行组合来提取合适的局部几何结构，最大特征值所对应的特征向量表征了从一种材质过度到另一种材质时最为迅速的变化方向，而最小特征值所对应的方向则指明了某种材质在图像内部流动的方向；因此可以对三维图像内部的材质属性变化进行定量刻画；得到的三维离散数据点云模型作为整个骨架提取的输入数据；

步骤(3)、点云模型分割：基于步骤(1)中的数据对数据进一步进行优化，利用K-means对步骤(1)中数据进行分类，得到一个基于灰度值的分类结果数据，以此作为骨架提取的最终输入数据；

步骤(4)、骨架结构提取：步骤(3)中得到的分类结果数据作为骨架提取的输入数据，基于步骤(3)中的分类结果数据，利用L1-中值理论，并定义一个能量泛函，通过求解能量最小的中心点来代表骨架分支；通过结合正则化，得到一个一维连续的骨架结构；

以上的所有操作都基于交互的模式，用户通过调整参数，提取合适的局部特征结构，而且对于数据中含有噪声和孤立点，甚至对于数据部分缺失都能提取一个完整的骨架借结构。

2.根据权利要求1所述的基于医学图像的器官骨架提取方法，其特征在于：步骤(2)中所述的提取三维数据的局部特征结构，该方法使用一组从CT上读取的原始扫描数据，利用三维数据的局部几何结构，通过求解每个三维数据的黑塞矩阵的特征值，将其提取局部特征结构转化为选取满足一定特征值条件的三维数据。

3.根据权利要求1所述的基于医学图像的器官骨架提取方法，其特征在于：步骤(3)中所述的点云模型分割方法，利用K-means聚类对步骤(2)中得到具有几何结构的数据进行分类，减少分支间在提取骨架时的互相干扰。

4.根据权利要求1所述的基于医学图像的器官骨架提取方法，其特征在于：步骤(4)中所述的骨架结构提取利用L1-中值理论，找出一组局部邻域最小值点作为骨架分支点、定义能量泛函来刻画局部邻域最小值，并通过最小化能量泛函求出骨架的分支点；利用L1-中值理论求局部邻域最小值点对噪声不敏感，且当点的维数大于二时，求得的局部领域最小值点唯一；由于单独应用L1-中值理论，会造成中心点聚积到一起从而产生一个稀疏的分布，所以需要用正则化理论控制中心点的分布，来实现骨架结构提取。

5.根据权利要求1所述的基于医学图像的器官骨架提取方法，其特征在于：步骤(2)(3)(4)中所述所有操作都是基于交互的模式；步骤(2)中可以选择点的黑塞矩阵特征值满足不同条件的三维离散点来实现提取感兴趣的局部结构；步骤(3)中用户可以选择类的个数，不同的类可以产生不同的效果；步骤(4)用户可以选择采样的数量等利于优化骨架结构的参数；所有操作和交互都是利用OpenGL和Qt技术进行处理。