CN102163327A - 一种医用的心脏ct图像分割方法 - Google Patents

Abstract

本发明公开了一种医用的心脏CT图像分割方法，包括以下步骤：读取DSCT采集的心脏断层图像，将N×N矩阵形式表示的待分割图像转换为以N²个元素表示的一维CT数值数组；对所述一维CT数值数组进行k聚类分析运算，计算图割能量函数所需的数据项和平滑项；构造所述待分割图像所表达的有向图，并根据数据项和平滑项建立能量函数，使用最大流/最小割算法最小化所述能量函数，求取所述有向图的最小割作为所述待分割图像的分割结果。本方法实现了心脏断层CT图像重点区域的精确分割，为重建心脏解剖结构，给与临床直观的参考数据和评价手段提供良好的基础。

Description

一种医用的心脏CT图像分割方法

技术领域

本发明涉及一种医学图像处理方法，具体来说涉及一种医用的心脏CT图像分割方法。

背景技术

随着人类生活水平的提高和预期寿命的延长，心血管疾病成为人类的头号死因，因此心血管疾病的早期诊断可有效降低病死率。了解心脏结构的影像学表现及其功能数据是正确诊断心脏疾病的重要前提，CT技术的发展，明显提高了时间分辨力，减少了心脏搏动伪影，在显示心脏细微结构方面显示出良好的应用潜力。

临床医学图像具有极其繁杂的多样性和复杂性，由于各类成像设备的固有特性，导致图像存在一定的噪声、场偏移效应和体效应，导致感兴趣目标物体的部分边缘存在局部不清晰，使得医学图像特别是CT图像分割一直是图像分割领域的难题。医学图像分割常用有阈值分割、基于边缘分割、基于区域分割和基于形态学分割等几种方法，但这些方法在分割高噪声、低对比度的心脏图像时难以达到理想的效果。基于广义梯度矢量流Snake模型的心脏MR图像分割方法，实现了MR左心室图像的精确分割。图割(Graph-Cuts)方法广泛运用于物理、化学、控制论、网络理论、计算机科学等许多领域，2001年BoyKov Y等，提出了最大流/最小割的快速计算方法，实现基于能量最小化进行目标提取的框架。这一方法具有快速、鲁棒、全局最优、抗噪强、扩展性好的特点，不仅可以用于图像分割，也可为其它视觉问题的处理提供统一的思路。k-均值聚类算法由Mac Queen于1967年提出，是聚类分析中一种基本的划分方法，方法简单、收敛速度快、并具有很强的局部搜索力，能有效地处理大数据集。

以下内容为现有CT图像分割方法的相关内容：

CT值：

CT值代表X线穿过组织被吸收后的衰减值。为了定量衡量组织对于X光的吸收率，Hounsfield定义了一个新的标度“CT值”。为了表示对他的敬意，后人将CT值的单位定为“HU”(Hounsfield Unit)。

CT值的计算：某物质的CT值等于该物质的衰减系数与水的衰减系数之差，再与水的衰减系数之比后乘以1000。即某物质CT值＝1000×(u-u水)/u水，其单位名称为HU，CT值不是一个绝对值，而是一个相对值。不同组织的CT值各异，各自在一定范围内波动。其中骨骼的CT值最高，为1000HU，软组织的CT值为20～70HU，水的CT值为0(±10)HU，脂肪的CT值为-50～-100HU以下，空气的CT值为-1000HU.各向异性扩散滤波算法：

Perona和Malik提出了基于偏微分方法的各向异性扩散方程，即在非线性尺度空间，根据图像不同方向上的梯度确定扩散系数。这一算法在滤除噪声和保留细节方面有了很大的提升，P-M方程表达式为：

$\frac{\∂u}{\∂t}=\mathrm{div}[c\left(\right||\▿u|\left|\right)\·\▿u]---\left(1\right)$

其中，div为散度算子，u为图像像素的灰度值，表示梯度，t是物理热扩散时间，在程序设计中采用迭代步长表示，c是扩散系数，是图像梯度

的非负单调递减函数。P-M方程的离散化形式如下：

$u_q^{t+1}=u_q^t+\mathrm{\λ}\underset{p\∈N_q}{\mathrm{\Σ}}c\left({\▿u}_{q,p}\right){\▿u}_{q,p}---\left(2\right)$

其中q为离散化图像中像素的位置，u_q为图像中像素q处的值，

t是迭代的次数，λ为决定扩散速率的常数，N_q是像素q的邻域。

k聚类算法

聚类算法属于无监督学习，不依赖预定义的类和类标号，在处理海量数据的时候，代价较小。K聚类算法通过定量计算闵可夫斯基距离，确定两个可比较元素间的相异度，闵可夫斯基距离E的定义如下：

$E=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\left|\right|C_j-X_i\left|\right|}^p---\left(3\right)$

其中k为聚类中心点的个数，L为总的像素点个数。常用作度量标量相异度的有欧氏距离、曼哈顿距离等，其中欧氏距离利曼哈顿距离可认为是闵可夫斯基距离在p＝2和p＝1下的特例。

K聚类算法通过设定预分类数量k；然后从n个数据对象任意选择k个对象作为初始聚类中心；剩下其它数据对象，则根据它们与这些聚类中心的相似度，分别分配给与其最相似的聚类中心，如下所示则X_p属于θ_j(m)类：

X_p∈θ_j(m)if||X_p-C_j(m)||＜||X_p-C_i(m)||for i＝1，2，...，k，i≠j；(4)

然后再计算每个所获新聚类的聚类中心(聚类中心对应于该聚类中所有对象的均值)；不断重复这一过程直到标准测度函数开始收敛，各中心点不再变化为止。k个聚类具有各聚类本身尽可能的紧凑，而个同聚类之间尽可能分开的特点。

图割(Graph-Cuts)算法

图论算法中，一个图G由顶点集合V和边的集合E构成。流网络是一个有向图，其中给出两个特殊的顶点，分别称为源点和汇点。网络中的每条有向边(p，q)∈E均有一非负容量c(p，q)≥0。一个典型的流网络图见图1，其中s，t分别为源点和汇点，标记出的边为双向边，其中斜线左边的数字表示实际边上的流，右边的数字表示边的最大容量。

流出s或流进t的总流量被称为网络流量的值，记为|f|。网络G中的一个割(S，T)将V划分为不相连的两部分S和T＝V-S，使得s∈S，t∈T。割的容量定义为

$c(S,T)=\underset{p\∈S}{\mathrm{\Σ}}\underset{q\∈T}{\mathrm{\Σ}}c(p,q)---\left(5\right)$

最大流/最小割定理证明一个流网络的最大流量值等于最小割的容量，有max(|f|)＝min(c(S，T))。除s，t外，用网络中的每个节点来代表图像的像素点。像素节点之间的连接称为n-连接，像素节点与s和t的连接称为t-连接。n-连接代表了图像的邻域系统，而t-连接代表了像素与标记的关系。定义s为bck标记，t为obj标记，λV_ij为任意两相邻位置i和j之间n-连接的容量，D_i(W_i＝obj)为节点i和t之间t-连接的容量，D_i(W_i＝bck)为节点i和s之间t-连接的容量。则该网络的最小割(S，T)定义了最小化能量的一个分割。根据最大流/最小割的等价性定理，可以采用求解网络最大流的算法来获得最小割。

CT图像可以转化为用G(V，E)表示的有向图网络，其中V表示一系列顶点，对应于CT图像中的像素点，E表示一系列边，对应于像素之间的连接，(p，q)表示从顶点p到顶点q，c(p，q)表示连接顶点p到顶点q边的权重。一股地，图像像素对应于网络的顶点，像素特征之间的差异或相似性对应边上的容量，通过以上方法将图像映射为网络，使得网络割的容量对应视觉问题的能量函数。

假设上述有向图G的顶点的标号用f来表示，其中f∈{1，2，...，q}，q为总的标号数，则整个图G各顶点的状态可表示为F＝{f₁，f₂，...，f_N}，N为总的顶点数。

为计算各像素的标号，首先构造一个关于标号状态F的能量函数，根据数据约束和光滑约束，基于Potts模型的能量函数为：

$E\left(F\right)=E_{\mathrm{data}}\left(F\right)+E_{\mathrm{smooth}}\left(F\right)=\underset{p\∈P}{\mathrm{\Σ}}{D}_p\left(f_p\right)+\underset{\{p,q\∈N\}}{\mathrm{\Σ}}{D}_{p,q}(f_p,f_q)---\left(6\right)$

其中，E_data(F)为数据项，用于衡量F和所观察到的数据的不一致性，E_smooth(F)为平滑项，用于衡量F非分片光滑的程度，N表示相邻像素对构成的集合，D_p(f_p)为像素p为标号f_p的概率，D_p，q(f_p，f_q)表示像素p和q之间的相似度。

基于图割(Graph Cuts)的图像分割的基本概念就是把图像的分割描述为基于能量函数的最小化问题，如图2所示，将图像中的每一个像素设为图中的一个结点，再各设一个虚拟的源点S、汇点t。

为获得最优分割，需要最小化全局能量式。以上形式的能量最小化问题，转化为图割问题求解。图割(Graph-Cuts)中源点{s}和汇点{t}，对应于分割所用的目标像素点和背景像素点。通过CT图像到有向图G的转换，图的分割就是关于所有顶点针对目标和背景像素点的划分，作为图像分割就是要找到关于G的最小代价分割。

根据临床经验，心脏CT成像中，心肌组织CT值一股为-240左右。由于心脏结构的复杂性，以及周围软组织与心肌CT值相似，分割难度较高。图割(Graph Cuts)在仅包含强度信息的医学图像分割中很容易在待分割对象的弱边界、闭塞处分割失败，尤其是图像中不连续的部份有可能分割成相同的物体，同时算法难于区分强度信息接近的心脏周围组织。而公式(6)的能量函数要求目标特征不同于背景特征，才能取得较好的分割效果。

发明内容

本发明的目的在于提供一种医用的心脏CT图像分割方法，本方法可对西门子双源CT(Dual source CT，DSCT)心脏图像进行分割，具有较好的分割结果，且具有计算简单，速度较快，鲁棒性好的特点，实现了心脏断层CT图像重点区域的精确分割，为重建心脏解剖结构，给与临床直观的参考数据和评价手段提供良好的基础。

本发明的目的可通过以下的技术措施来实现：

一种医用的心脏CT图像分割方法，包括以下步骤：

(1)读取DSCT采集的心脏断层图像，作为待分割图像；

(2)将N×N矩阵形式表示的待分割图像转换为以N²个元素表示的一维CT数值数组；

(3)对所述一维CT数值数组进行k聚类分析运算，并得到所述待分割图像中的每个像素对应的聚类号；

(4)根据步骤(3)的k聚类分析运算结果，计算图割能量函数所需的数据项和平滑项；

(5)构造所述待分割图像所表达的有向图，并根据步骤(4)中的数据项和平滑项建立能量函数，使用最大流/最小割算法最小化所述能量函数，求取所述有向图的最小割作为所述待分割图像的分割结果。

所述步骤(3)中的k聚类分析运算过程为：对待分割图像的一维CT数值数组中的每个数值进行K均值聚类，并根据聚类结果为每个像素标上所属聚类标号，从而得到CT图像中各像素聚类标号，聚类标号为1～k；所述K值取为4。

上述K均值聚类为Matlab的kmeans函数：其中kmeans函数中的像素差异度采用平方欧几里德距离计算；逐个将待分割图像中各像素的CT数值按照平方欧几里德距离结果分配给k个聚类区域中的一个；每迭代计算完一次，kmean函数重新计算聚类中心，直到聚类中心不再变化或达到迭代次数为止，则算法收敛，得到各类的聚类中心和CT图像各像素点所属聚类的标号；对于每个聚类依次赋予标号1至k，并将待分割图像的每个像素标记为对应的聚类号，则待分割图像被映射到一个标记图上，此标记图即为下一步图割计算所需的有向标号图；所述标记图中的每个像素在1至k之间取值，该值等于该像素对应的聚类号。

所述步骤(4)中的数据项和平滑项的计算过程为：

(a)数据项的计算采取对具有相同聚类标号的图像CT值进行协方差矩阵计算；得到从统计意义上反映图像CT值相关性程度的信息，反映的是待分割图像中每个像素归属于某一聚类区域的概率。

(b)平滑项采用高斯(gauss)和索贝尔(sobel)算子的二维卷积对CT图像进行平滑和边缘检测，计算得到平滑项的结果。

所述步骤(5)中构造待分割图像表达的有向图的过程：所述有向图中的每个顶点对应于所述待分割图像的一个像素点，所述的顶点中包括两个终点顶点以及其余的非终点顶点，其中两个终点分别为源点和汇点；所述有向图中的边包括一组n-连接和一组t-连接，所述的每一组n-连接用于连接一对非终点顶点，对应于有向图中相邻对应像素点的连接；所述每一组t-连接用于连接一个非终点顶点和一个终点；

所述步骤(5)中最小化步骤(4)中的能量函数的过程为：

(51)为所述待分割图像所对应的有向图中的每一个n-连接和t-连接均分配一个容量；其中，包含源点的t-连接的容量为R_p(″bkg″)，包含汇点的t-连接的容量为R_p(″obj″)，n-连接的容量为βB_p，q，其中p、q为该连接的两个顶点所对应的像素标号，具体为：

R_p(″obj″)＝-ln Pr(I_p|″obj″)

R_p(″bkg″)＝-ln Pr(I_p|″bkg″)

$B_{p,q}=\left\{\begin{array}{c}1,I_p\≤I_q\\ \exp (-\frac{{(I_p-I_q)}^2}{{2\mathrm{\σ}}^2})\frac{1}{\left|\right|p-q\left|\right|},I_p>I_q\end{array}\right.$

P_r(I_p|″obj″)表示像素属于某一聚类的概率，即像素p的CT值I_p属于这一聚类的概率，；P_r(I_p|″bkg″)为像素属于其它聚类的概率，即像素p的CT值I_p属于其它聚类的概率，计算每个聚类区域的协方差矩阵，统计其像素点CT值直方图，并将其转化成概率得到；||p-q||为像素p，q在带分割图像中的欧式距离；σ²为像素p，q的CT值相似程度的度量，其值为：

σ²＝E((I_p-I_q)²)

其中E(S)表示S的数学期望。

(52)根据步骤(51)得到的各个连接的容量运用最小割、最大流定理将所述有向图中的顶点划分为子集，每个子集包含一个终点；最终实现待分割图像的目标和背景的分割。

所述步骤(2)中的待分割图像I对应于512×512的矩阵形式，且所述矩阵所对应的一维CT数值数组包括262144个元素。

本方法的医用的心脏CT图像分割方法所述的步骤(1)进行之前还包括对待分割图像的去噪滤波和增强预处理。

所述去噪滤波预处理过程中采用各向异性扩散滤除同质区域的噪声，并采用中值滤波去除脉冲噪声；所述各向异性扩散采用图像在不同方向上梯度的单调函数作为扩散系数。

本发明对比现有技术，有如下优点：本发明方法结合Graph-cuts和k聚类算法的优点，通过这两种方法有针对性地对CT图像全局信息和局部信息进行提取，结合两种算法的优点，将复杂的问题分解，提高了图像分割的准确性和速度。

附图说明

图1是一个典型的流网络图：

图2是Graph Cuts算法构建有向图的示意；

图3是本发明CT图像分割方法的处理流程图；

图4(a)是一幅DSCT原图示意图；

图4(b)是对图4(a)所示原图进行各向异性滤波预处理后得到的图像；

图5(1)为另一幅DSCT待分割原图示意图；

图5(2)至图5(8)是采用本方法医用的心脏CT图像分割方法对图5(1)所示待分割图像分别采用聚类数k为2、3、4、5、6、7、8时进行分割处理后的效果图。

具体实施方式

根据临床经验，心脏CT成像中，心肌组织CT值一股为-240左右。由于心脏结构的复杂性，以及周围软组织与心肌CT值相似，分割难度较高。图割(GraphCuts)在仅包含强度信息的医学图像分割中很容易在待分割对象的弱边界、闭塞处分割失败，尤其是图像中不连续的部份有可能分割成相同的物体，同时算法难于区分强度信息接近的心脏周围组织。而能量函数(公式(6))要求目标特征不同于背景特征，才能取得较好的分割效果。

因此，本发明的医用的心脏CT图像分割方法中包括如下处理过程：

处理过程一：提供一种标记待分割CT图像中像素的方法：

为获得目标和背景的纹理统计模型，先对心脏CT图像的所有CT值进行K均值聚类。K个聚类中心代表典型的CT值。为了将目标和背景区分开，一个好的K值应尽可能地将目标和背景种子点分到不同的聚类中。一股将K值取为4，实验证明这个取值在保证速度的同时，均能取得较好的效果。

用K聚类算法得到每个像素点的初始标号。算法首先将图像CT值矩阵转换为CT值向量数组，大小为512*512，然后在向量数组中选择K个初始聚类中心，并确定迭代运算的次数；算法逐个将各个像素按照聚类准则分配给K个聚类中心中的一个；每聚类完一次，就重新计算聚类中心，直到聚类中心不再变化或达到迭代次数为止，则算法收敛，得到各类的聚类中心和CT图像各像素点所属聚类的标号。

对于每个聚类依次赋予标号1-K，并将图像I的每个像素标记为它的特征向量对应的聚类号，则I被映射到一个标记图上。标记图的每个像素在1-K之间取值，该值等于该像素对应的聚类号。

处理过程二：构造基于CT值的有向网络图：

(a)根据图2示例构造待分割的CT图像表达的有向图。有向图包含一系列顶点，每一个顶点对应于所表示CT图像的一个像素；两个终点(汇点t和源点S)，对应于分割所对应的目标种子点和背景种子点；一组n-连接，每一个连接一对顶点；一组t-连接，每个连接终点之一和图中的顶点之一；

(b)为每一个n-连接分配一个容量；为每一个t-连接分配一个容量，建立对应能量函数；

(c)运用最小割、最大流方法，将节点划分为子集，每一个子集包含一个终点。

(d)根据初始分割结果，添加新的目标或背景像素种子点，重新分派n-连接和t-连接容量。

(e)动态更新第(c)步的方法，计算改变的容量，重新划分节点，从而实现CT图像目标和背景的分割。

根据上述两个处理过程，具体对待分割图像进行分割的过程如下：

一、建立分割模型

针对CT图像的特点，心脏和周围组织具有比较接近的CT值。单元势能，即数据代价，为根据CT图像像素CT值确定的像素所属聚类簇的负对数概率。邻域势能，即通常所说的平滑项，目的是为了保留边缘，对边缘处的平滑进行惩罚，该势能度量邻域像素点CT值的差异，避免对CT值变化较大区域进行平滑。

本发明算法引入基于聚类结果的目标种子局部信息，能量函数中平滑项定义为：

E_p，q＝g(p，q)·|I_p-I_q|                    (7)

$g(p,q)=\exp (-\frac{{\left|\right|I\left(p\right)-I\left(q\right)\left|\right|}^2}{{2\mathrm{\σ}}^2})---\left(8\right)$

其中，I表示像素的强度，即CT值。

能量函数数据项定义为：

$\left\{\begin{array}{c}{E}_p(I_p=\mathrm{bkg}_{\′\′}^{\′\′})=-\log [P_r\left(I\left(p\right)\right|p\∈O)\×\exp (-\left(\frac{{d(p,O)}^2}{{\mathrm{\σ}}_a}\right))]\\ E_p(I_p=\mathrm{obj}_{\′\′}^{\′\′})=-\log \left[P_r\left(I\left(p\right)\right|p\∈B)\right]\end{array}\right.---\left(9\right)$

公式中，O和B分别对应在待分割图像I中标记出的目标和背景种子像素点集合，在本方法中目标和背景种子像素点分别对应为第一步聚类结果的聚类内部区域和聚类外部区域。根据所标记的目标和背景种子点，分别建立目标和背景的统计模型，其中σ_a＝10表示种子点的扩张和区域影响。函数d(p，O)＝min{d(p，q)|q∈O}表示像素p到目标种子点之间的距离，可采用欧几里德几何距离或者测地线距离进行计算。

数据项和平滑项的计算过程为：

(a)数据项的计算采取对具有相同聚类标号的图像CT值进行协方差矩阵计算；得到从统计意义上反映图像CT值相关性程度的信息，反映的是待分割图像中每个像素归属于某一聚类区域的概率。

(b)平滑项采用高斯(gauss)和索贝尔(sobel)算子的二维卷积对CT图像进行平滑和边缘检测，计算得到平滑项的结果。平滑项的计算主要是为了保留边缘。

二、模型的求解

对图像进行去噪和增强等预处理，为后续分割提供良好的图像基础，减少噪声等原因引起的分割误差。本方法中取二维图像上、下、左、右四邻域，以得到更好的效果。各向异性扩散采用图像在不同方向上梯度的单调函数作为扩散系数，在同质区域扩散系数较大，可有效滤除同质区域的噪声；在边缘区域，由于灰度值变化剧烈，梯度较大，扩散系数较小，可以较好地保留图像中的边缘信息。CT图像所含噪声主要为脉冲噪声和高斯噪声，为尽量保持边缘的清晰，采用中值滤波去除脉冲噪声，各向异性去噪法使用各向异性扩散的非线性滤波方法可以在抑制高斯噪声的同时保持细小结构。从图4(b)可看出图像中各区域边缘清晰，区域内部明显比原图4(a)平滑。

图3示出了本发明的CT图像分割方法的流程图，该流程也就是对上述建立的模型求解，具体过程如下：

(1)读取DSCT采集的心脏断层图像，作为待分割图像I；

(2)将512×512矩阵形式表示的待分割图像I转换为以262144个元素表示的一维CT数值数组；

(3)对一维CT数值数组进行k聚类分析运算，并得到待分割图像I中的每个像素对应的聚类号；k聚类分析运算过程即采用上述处理过程一，其中，本发明中根据具体试验的结果确定K值取为4；通过K均值聚类将CT图像中的所有像素(262144个像素)的CT数值划分为k个聚类，使得同一聚类中的像素相似度较高，不同聚类中的像素相似度较小。CT数值的K均值聚类主要用于CT图像的粗分割。

具体地，K均值聚类采用Matlab的kmeans函数：其中函数中的像素差异度采用平方欧几里德距离计算；逐个将待分割图像中各像素的CT数值按照平方欧几里德距离结果分配给k个聚类区域中的一个；每迭代计算完一次，由于聚类区域中分配了新的像素CT数值，就加入新的CT数值，kmean函数重新计算聚类中心，直到聚类中心不再变化或达到迭代次数为止，则算法收敛，得到各类的聚类中心和CT图像各像素点所属聚类的标号；对于每个聚类依次赋予标号1至k，并将待分割图像的每个像素标记为对应的聚类号，则待分割图像被映射到一个标记图上，用于作为后续图割运算的有向图；标记图中的每个像素在1至k之间取值，该值等于该像素对应的聚类号。

(4)根据步骤(3)的聚类分析运算结果，计算图割能量函数所需的数据项和平滑项；

(5)构造待分割图像I所表达的有向图，并根据步骤(4)中的数据项和平滑项建立能量函数，使用最大流/最小割算法最小化能量函数，求取图的最小割作为待分割图像I的分割结果。

步骤(5)中构造待分割图像表达的有向图的过程：根据前述图割算法有向图的构造方法，将上述的标记图转换为有向图，有向图中的每个顶点对应于待分割图像的一个像素点，顶点中包括两个终点顶点以及其余的非终点顶点，其中两个终点分别为源点和汇点；有向图中的边包括一组n-连接和一组t-连接，的每一组n-连接用于连接一对非终点顶点，对应于有向图中相邻对应像素点的连接；每一组t-连接用于连接一个非终点顶点和一个终点；

步骤(5)中最小化步骤(4)中的能量函数的过程为：

(51)为待分割图像所对应的有向图中的每一个n-连接和t-连接均分配一个容量；其中，包含源点的t-连接的容量为R_p(″bkg″)，包含汇点的t-连接的容量为R_p(″obj″)，n-连接的容量为βB_p，q，其中p、q为该连接的两个顶点所对应的像素标号，具体为：

R_p(″obj″)＝-ln Pr(I_p|″obj″)

R_p(″bkg″)＝-ln Pr(I_p|″bkg″)

$B_{p,q}=\left\{\begin{array}{c}1,I_p\≤I_q\\ \exp (-\frac{{(I_p-I_q)}^2}{{2\mathrm{\σ}}^2})\frac{1}{\left|\right|p-q\left|\right|},I_p>I_q\end{array}\right.$

P_r(I_p|″obj″)表示目标对象中像素p的CT值I_p出现的概率，即像素属于某一聚类的概率；P_r(I_p|″bkg″)为背景中像素p的CT值I_p出现的概率，即像素属于其他聚类的概率，它们的值都是在(3)聚类结果的基础上，计算每个聚类区域的协方差矩阵，统计其像素点CT值直方图，并将其转化成概率得到；||p-q||为像素p，q在带分割图像中的欧式距离；σ²为像素p，q的CT值相似程度的度量，其值为：

σ²＝E((I_p-I_q)²)                    (10)

其中E(S)表示S的数学期望。

(52)根据步骤(51)得到的各个连接的容量运用最小割、最大流定理将有向图中的顶点划分为子集，每个子集包含一个终点；最终实现待分割图像的目标和背景的分割。

下面采用本方法对图5(1)所示的DSCT原图进行分割处理，去噪和增强预处理后，如图5(2)至5(8)所示的为分别采用聚类数k为2、3、4、5、6、7、8时进行分割处理后的效果图，计算过程结果如下表1所示，随着聚类数的增加，由于与边界相关的边缘数量增加，平滑势能值显著增加，而数据势能与像素点CT值相关，因此变化不大，分割时计算时间随着总代价的增加，响应增加；当k取值4或5时，具有较好的分割结果。综合计算时间和分割质量，k值选定为4，即当k值取为4时得到最好的分割效果图。

区域数k

平滑代价

数据代价

计算时间(s)

2	0.0076	1.3092e+005	0.5798
				3	1.5529	1.2583e+005	0.9172
4	4.2714	1.1339e+005	1.3477
				5	4.7926	1.2318e+005	1.8809
6	18.4284	1.1816e+005	2.1190
				7	25.0507	1.1880e+005	2.5074
8	23.5622	1.2179e+005	3.0004

表1

本方法根据DSCT心脏图像高噪声、低对比度的特点，利用临床心室、心房的先验解剖结构知识，在阈值分割的基础上，为克服传统阈值对图像灰度分布不成明显双峰时分割不理想的缺点，将图像梯度算子与阈值分割结合起来，首先将DSCT心脏图像利用梯度算子计算得到梯度图像，然后在得到的梯度图像上运用k聚类算法进行分析获得各像素的初始标号，建立一个关于标号的能量函数，并构造相应网络，最后根据最大流/最小割算法最小化能量函数实现CT心脏图像精确分割。

本发明的实施方式不限于此，在本发明上述基本技术思想前提下，按照本领域的普通技术知识和惯用手段对本发明内容所做出其它多种形式的修改、替换或变更，均落在本发明权利保护范围之内。

Claims (10)

1.一种医用的心脏CT图像分割方法，其特征在于包括以下步骤：

(1)读取DSCT采集的心脏断层图像，作为待分割图像；

(2)将N×N矩阵形式表示的待分割图像转换为以N²个元素表示的一维CT数值数组；

(3)对所述一维CT数值数组进行k聚类分析运算，并得到所述待分割图像中的每个像素对应的聚类号；

(4)根据步骤(3)的k聚类分析运算结果，计算图割能量函数所需的数据项和平滑项；

(5)构造所述待分割图像所表达的有向图，并根据步骤(4)中的数据项和平滑项建立能量函数，使用最大流/最小割算法最小化所述能量函数，求取所述有向图的最小割作为所述待分割图像的分割结果。

2.根据权利要求1所述的医用的心脏CT图像分割方法，其特征在于：所述步骤(3)中的k聚类分析运算过程为：对待分割图像的一维CT数值数组中的每个数值进行K均值聚类，并根据聚类结果为每个像素标上所属聚类标号，从而得到CT图像中各像素聚类标号，聚类标号为1～k；所述K值取为4。

3.根据权利要求1所述的医用的心脏CT图像分割方法，其特征在于：所述K均值聚类采用Matlab的kmeans函数：其中kmeans函数中的像素差异度采用平方欧几里德距离计算；逐个将待分割图像中各像素的CT数值按照平方欧几里德距离结果分配给k个聚类区域中的一个；每迭代计算完一次，kmean函数重新计算聚类中心，直到聚类中心不再变化或达到迭代次数为止，则算法收敛，得到各类的聚类中心和CT图像各像素点所属聚类的标号；对于每个聚类依次赋予标号1至k，并将待分割图像的每个像素标记为对应的聚类号，则待分割图像被映射到一个标记图上，此标记图即为下一步图割计算所需的有向标号图；所述标记图中的每个像素在1至k之间取值，该值等于该像素对应的聚类号。

4.根据权利要求2或3所述的医用的心脏CT图像分割方法，其特征在于：所述步骤(4)中的数据项和平滑项的计算过程为：

(a)数据项的计算采取对具有相同聚类标号的图像CT值进行协方差矩阵计算；得到从统计意义上反映图像CT值相关性程度的信息，反映的是待分割图像中每个像素归属于某一聚类区域的概率。

(b)平滑项采用高斯和索贝尔算子的二维卷积对CT图像进行平滑和边缘检测，计算得到平滑项的结果。

5.根据权利要求4所述的医用的心脏CT图像分割方法，其特征在于：所述步骤(5)中构造待分割图像表达的有向图的过程：所述有向图中的每个顶点对应于所述待分割图像的一个像素点，所述的顶点中包括两个终点顶点以及其余的非终点顶点，其中两个终点分别为源点和汇点；所述有向图中的边包括一组n-连接和一组t-连接，所述的每一组n-连接用于连接一对非终点顶点，对应于有向图中相邻对应像素点的连接；所述每一组t-连接用于连接一个非终点顶点和一个终点；

6.根据权利要求5所述的医用的心脏CT图像分割方法，其特征在于：所述步骤(5)中最小化步骤(4)中的能量函数的过程为：

(51)为所述待分割图像所对应的有向图中的每一个n-连接和t-连接均分配一个容量；

(52)根据步骤(51)得到的各个连接的容量运用最小割、最大流定理将所述有向图中的顶点划分为子集，每个子集包含一个终点；最终实现待分割图像的目标和背景的分割。

7.根据权利要求6所述的医用的心脏CT图像分割方法，其特征在于：所述步骤(51)中的n-连接和t-连接的容量为：其中，包含源点的t-连接的容量为R_p(″bkg″)，包含汇点的t-连接的容量为R_p(″obj″)，n-连接的容量为βBp，q，其中p、q为该连接的两个顶点所对应的像素标号，具体为：

R_p(″obj″)＝-ln Pr(I_p|″obj″)

R_p(″bkg″)＝-ln Pr(I_p|″bkg″)

$B_{p,q}=\left\{\begin{array}{c}1,I_p\≤I_q\\ \exp (-\frac{{(I_p-I_q)}^2}{{2\mathrm{\σ}}^2})\frac{1}{\left|\right|p-q\left|\right|},I_p>I_q\end{array}\right.$

P_r(I_p|″obj″)表示像素属于某一聚类的概率，即像素p的CT值I_p属于这一聚类的概率，；P_r(I_p|″bkg″)为像素属于其它聚类的概率，即像素p的CT值I_p属于其它聚类的概率，计算每个聚类区域的协方差矩阵，统计其像素点CT值直方图，并将其转化成概率得到；||p-q||为像素p，q在带分割图像中的欧式距离；σ²为像素p，q的CT值相似程度的度量，其值为：

σ²＝E((I_p-I_q)²)

其中E(S)表示S的数学期望。

8.根据权利要求7所述的医用的心脏CT图像分割方法，其特征在于：所述步骤(2)中的待分割图像I对应于512×512的矩阵形式，且所述矩阵所对应的一维CT数值数组包括262144个元素。

9.根据权利要求1或8所述的医用的心脏CT图像分割方法，其特征在于：所述的步骤(1)进行之前还包括对待分割图像的去噪滤波和增强预处理。

10.根据权利要求9所述的医用的心脏CT图像分割方法，其特征在于：所述去噪滤波预处理过程中采用各向异性扩散滤除同质区域的噪声，并采用中值滤波去除脉冲噪声；所述各向异性扩散采用图像在不同方向上梯度的单调函数作为扩散系数。

Images