CN103729875B - 心脏磁共振图像的左心室三维轮廓重建方法和系统 - Google Patents

Abstract

本发明涉及一种心脏磁共振图像的左心室三维轮廓重建方法和系统。所述方法包括：建立心脏磁共振图像的混合高斯模型的步骤；初始化活动轮廓模型的步骤；确定左心室内外表面轮廓的步骤；以及重建左心室三维轮廓的步骤。上述心脏磁共振图像的左心室三维轮廓重建方法和系统，采用混合高斯模型将磁共振图像划分成多个区域，再采用移动正方形法初始化活动轮廓模型，采用活动轮廓模型的能量最小化方程求解得到左心室内外表面轮廓，再由左心室内外表面轮廓重建三维轮廓，因采用混合高斯模型及活动轮廓模型求取的左心室内外表面轮廓较为准确，进而重建的三维轮廓的准确度也较高。

Description

心脏磁共振图像的左心室三维轮廓重建方法和系统

技术领域

本发明涉及图像处理领域，特别是涉及一种心脏磁共振图像的左心室三维轮廓重建方法和系统。

背景技术

随着生活水平的提高，心血管疾病已成为人类死亡的主要原因之一。尽早确诊疾病是提高治愈和减少发病危险的一个重要因素。为了了解心脏的内部情况，采用磁共振成像技术对心脏扫描成像。磁共振图像可观察到一个心跳周期内左心室的运动情况。

传统的心脏磁共振图像一般为二维图像，对该二维图像进行分割可得到心脏的轮廓，然而该心脏轮廓分割准确度低，且无法准确重建左心室三维轮廓。

发明内容

基于此，有必要提供一种能提高准确度的心脏磁共振图像的左心室三维轮廓重建方法。

此外，还有必要提供一种能提高准确度的心脏磁共振图像的左心室三维轮廓重建系统。

一种心脏磁共振图像的左心室三维轮廓重建方法，包括：

建立心脏磁共振图像的混合高斯模型的步骤，对获取的心脏磁共振图像建立混合高斯模型，采用预设数量高斯分布将所述心脏磁共振图像划分成相应的灰度区域图像；

初始化活动轮廓模型的步骤，采用移动正方形法对所述灰度区域图像的边缘进行处理得到初始化的活动轮廓模型；

确定左心室内外表面轮廓的步骤，根据所述活动轮廓模型的能量最小化方程分别求解得到所述左心室内外表面轮廓；以及

重建左心室三维轮廓的步骤，获取一个心跳周期内多幅空间位置相邻的心脏磁共振图像的横轴图像，获取每幅横轴图像的左心室内外表面轮廓，将多个横轴图形的左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓。

一种心脏磁共振图像的左心室三维轮廓重建系统，包括：

模型建立模块，用于对获取的心脏磁共振图像建立混合高斯模型，采用预设数量高斯分布将所述心脏磁共振图像划分成相应的灰度区域图像；

初始化活动轮廓模块，用于采用移动正方形法对所述灰度区域图像的边缘进行处理得到初始化的活动轮廓模型；

左心室内外轮廓获取模块，用于根据所述活动轮廓模型的能量最小化方程分别求解得到所述左心室内外表面轮廓；以及

三维轮廓重建模块，用于获取一个心跳周期内多幅空间位置相邻的心脏磁共振图像的横轴图像，获取每幅横轴图像的左心室内外表面轮廓，将多个横轴图形的左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓。

上述心脏磁共振图像的左心室三维轮廓重建方法和系统，采用混合高斯模型将磁共振图像划分成多个区域，再采用移动正方形法初始化活动轮廓模型，采用活动轮廓模型的能量最小化方程求解得到左心室内外表面轮廓，再由左心室内外表面轮廓重建三维轮廓，因采用混合高斯模型及活动轮廓模型求取的左心室内外表面轮廓较为准确，进而重建的三维轮廓的准确度也较高。

此外，采用重建的左心室三维轮廓测量出泵血量和泵血率参数，测量方便且准确。

附图说明

图1为一个实施例中心脏磁共振图像的左心室三维轮廓重建方法的流程图图；

图2为建立心脏磁共振图像的混合高斯模型的流程图；

图3为四种加入短线段的情况；

图4为移动正方形法初始化的心脏轮廓示意图；

图5A为采用Canny边界检测方法检测的边界线条；

图5B为边界提供给活动轮廓模型的受力情况；

图6A为活动轮廓模型分割的左心室中部磁共振图像；

图6B为左心室尖端磁共振图像；

图7为左心室内外表面轮廓三维网格重建的结果；

图8为点云法测量体积的示意图；

图9为该测量心室的泵血量和泵血率的流程图；

图10为左心室三维轮廓划分的等距的小方格；

图11A为穿过点的直线与网格的交点的示意图；

图11B为另一穿过点的直线与网格的交点的示意图；

图11C为另一穿过点的直线与网格的交点的示意图；

图12为一个实施例中心脏磁共振图像的左心室三维轮廓重建系统的结构框图；

图13为图12中模型建立模块的内部结构框图；

图14为另一个实施例中心脏磁共振图像的左心室三维轮廓重建系统的结构框图；

图15为图14中参数获取模块的内部结构框图。

具体实施方式

下面结合具体的实施例及附图对一种心脏磁共振图像的左心室三维轮廓重建方法和系统的技术方案进行详细的描述，以使其更加清楚。

磁共振图像是在电流激活时偏离原来的方向，电流撤去后氢原子逐渐回到原来的方向，放出一个同频率的信号，在梯度磁场中氢原子转速的不同选取不同位置拍摄得到的图像。

如图1所示，为一个实施例中心脏磁共振图像的左心室三维轮廓重建方法的流程图。该心脏磁共振图像的左心室三维轮廓重建方法，包括：

步骤S102，建立心脏磁共振图像的混合高斯模型的步骤，对获取的心脏磁共振图像建立混合高斯模型，采用预设数量高斯分布将该心脏磁共振图像划分成相应的灰度区域图像。

在一个实施例中，如图2所示，该建立心脏磁共振图像的混合高斯模型的步骤包括：

步骤S202，将心脏磁共振图像分成预设数量个高斯分布的灰度区域图像，建立混合高斯模型，计算混合高斯模型的先验概率。

因心脏磁共振图像的左心室内壁边界比较清晰可分，心脏磁共振图像按灰度大体上可分为三个区域，灰度最亮的区域为肺部的气管、心室内部的血液和脂肪，灰度比较暗的区域为肺部气泡区域，灰度处于中间的区域为心肌。心脏磁共振图像的灰度可采用三个高斯分布表示，用期望最大化可求解混合高斯模型的参数，该混合高斯模型的参数包括均值及方差。

本实施例中，将心脏磁共振图像分成三个区域，然后计算每个区域的先验概率。该先验概率可通过三种方法得到：第一种，对整个心脏磁共振图像的高斯拟合结果加随机扰动假设得到方差；第二种，采用K均值方法找到各区域中心，通过高斯拟合得到方差；第三种，通过经验抽样选取点作为区域中心，通过高斯拟合得到方差。该先验概率为采用三种方式根据以往数据分析得到的概率。

步骤S204，计算心脏磁共振图像中每个点属于某个区域的概率。

具体的，心脏磁共振图像的灰度为一个统计分布，灰度为x的点在区域内出现的概率为：

$p\left(x\right|i)=-\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_i}\exp (-\frac{|x-{\mathrm{\μ}}_i{|}^2}{2{\mathrm{\σ}}_i^2})---\left(1\right)$

其中，p(x|i)为灰度为x的点在i区域的概率，μ_i为在i区域的均值，σ_i ²为在i区域的方差，σ_i为标准差。

步骤S206，重新计算每个区域中各点的灰度均值和方差，然后返回步骤S204。

步骤S208，根据每个点属于某个区域的概率采用贝叶斯公式计算后验概率，并将该后验概率作为下次计算先验概率的先验假设。

具体的，计算每个区域中各点的灰度均值公式为式(2)：

μ_i＝∑_xxp(x|i) (2)

计算每个区域中各点的方差为式(3)：

${\mathrm{\σ}}_i^2=E\left(\right|x-{\mathrm{\μ}}_i{|}^2)={\mathrm{\Σ}}_x|x-{\mathrm{\μ}}_i{|}^{2}p\left(x\right|i)---\left(3\right)$

采用贝叶斯公式计算后验概率，计算公式为式(4)：

$p\left(i\right|x)=\frac{p\left(x\right|i)P\left(i\right)}{{\mathrm{\Σ}}_i^{M}p\left(x\right|i)P\left(i\right)}---\left(4\right)$

其中，p(i)＝∑_xp(i|x)，M为总区域数，i表示第i个区域。

将后验概率作为下一轮先验概率的先验假设。

循环计算直到高斯模型的均值及方差变化量趋近于零，将各点分入后验概率最大值的区域，得到更新后的心脏磁共振图像的预设数量个灰度区域图像。

具体的，通过循环计算直到高斯模型的均值及方差变化量趋近于零，将各点分入到后验概率最大值的区域。如此，可得到更新后的心脏磁共振图像的三个灰度区域图像，即将所有点分别标记对应的区域标识，如属于第一区域的点采用1标记，属于第二区域的点采用2标记，属于第三区域的点采用3标记。

步骤S104，初始化活动轮廓模型的步骤，采用移动正方形法对该灰度区域图像的边缘进行处理得到初始化的活动轮廓模型。

在一个实施例中，步骤S104包括：将该灰度区域图像以预设移动正方形大小的网格划分，在该网格中插入短线段，将插入的短线段连接成曲线，取最长的圆形曲线初始化活动轮廓模型。

具体的，将左心室心肌内部用斜线表示，左心室心肌外部用空白表示的图像以网格划分，根据每个网格四个顶点的值，加入顶点在网格正方形边中点的短线段。该加入短线段的方法有多种，如图3所示，为四种加入短线段的情况。将插入的短线段连接成长度不同的曲线，取最长的圆形曲线初始化活动轮廓模型。如图4所示为移动正方形法初始化的心脏轮廓示意图。图5A为采用Canny边界检测方法检测的边界线条，图5B为边界提供给活动轮廓模型的受力情况。

步骤S106，确定左心室内外表面轮廓的步骤，根据该活动轮廓模型的能量最小化方程分别求解得到该左心室内外表面轮廓。

在一个实施例中，步骤S106包括：获取并根据活动轮廓模型的内力函数、图像力函数和交互作用力函数通过微积分建立活动轮廓模型的能量最小化方程，沿能量梯度方向移动，直到图像力函数与内力函数的值平衡，得到左心室内外表面轮廓。

具体的，活动轮廓模型的能量最小化方程如式(5)所示：

$E_{snake}={\∫}_0^1\left(E_{\mathrm{int}ernal}\right(v\left(s\right))+E_{image}\left(v\right(s\left)\right)+E_{constraint}\left(v\left(s\right)\right))ds---\left(5\right)$

式(5)中，E_snake为活动轮廓模型的能量；v(s)为活动轮廓上等距离点的坐标值；E_internal为活动轮廓内力函数，该内力受到曲线收缩和曲线平滑的限制；E_image为活动轮廓的图像力函数，该图像力函数是指活动轮廓模型移动到的位置的图像距图像边缘的图像灰度梯度力，将图像力函数作为活动轮廓模型的外力；E_constraint为交互作用力函数。图像力函数与内力函数的值平衡是指在外力和内力的作用下活动轮廓逐渐停止移动。

求解该活动轮廓的能量最小化方程，使得能量最小化，即轮廓的一阶导数、二阶导数和图像力的和最小。将能量最小化方程通过欧拉公式转化为式(6)：

$\begin{array}{c}{\mathrm{\αx}}_{ss}+{\mathrm{\βx}}_{ssss}+\∂E_{imag}/\∂x=0\\ {\mathrm{\αy}}_{ss}+{\mathrm{\βy}}_{ssss}+\∂E_{imag}/\∂y=0\end{array}---\left(6\right)$

式(6)中，α和β为多项式系数，x_ss和y_ss分别为x和y方向上的二次差分，x_ssss和y_ssss分别为x和y方向上的四次差分，是图像灰度梯度力。

对式(6)进行循环求解，将导数用邻近点的差分表示。活动轮廓点上的系数方程是一个由活动轮廓的一阶导数和二阶导数组成的矩阵在活动轮廓移动前计算得到的，在活动轮廓移动的过程中保持不变。如式(7)：

$\begin{array}{c}{x}_t={(A+\mathrm{\γ}I)}^{-1}({\mathrm{\γx}}_{t-1}-f_x(x_{t-1},y_{t-1}))\\ y_t={(A+\mathrm{\γ}I)}^{-1}({\mathrm{\γy}}_{t-1}-f_y(x_{t-1},y_{t-1}))\end{array}---\left(7\right)$

式(7)中，(f_x(x_t-1,y_t-1)，f_y(x_t-1,y_t-1))为(t-1)轮循环时的图像梯度向量，(x_t,y_t)为第t轮循环的计算活动轮廓模型每点的坐标值，(x_t-1,y_t-1)第(t-1)轮循环的计算活动轮廓模型每点的坐标值，A、γ、I均为系数。

图6A为活动轮廓模型分割的左心室中部磁共振图像，图6A中两圆形曲线为内外表面轮廓曲线；图6B为左心室尖端磁共振图像，图6B中两圆形曲线为内外表面轮廓曲线。

步骤S108，重建左心室三维轮廓的步骤，获取一个心跳周期内多幅空间位置相邻的心脏磁共振图像的横轴图像，获取每幅横轴图像的左心室内外表面轮廓，将多个横轴图像的左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓。

因内外轮廓之间是心肌部分，内外壁都会在运动过程中随着心脏的收缩而收缩。获取一个心跳周期内多幅空间位置相邻的磁共振图像的横轴图像，将右心室的横轴图像上的尖角作为重建三维网格的标识点，每个横轴图像上获取的左心室内外表面轮廓都以该标识点对齐，相邻横轴图像上的轮廓用三角面片连接起来，然后按照时间序列将多个横轴图像的左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓，如图7所示为左心室内外表面轮廓三维网格重建的结果。例如，每幅横轴图像的轮廓上取50个点，且每幅横轴图像的内或外表面轮廓由与左右心室间隔上右心室月牙形心后侧尖角标识点为内或外表面轮廓上的第一点组成三角网格。

上述心脏磁共振图像的左心室三维轮廓重建方法，采用混合高斯模型将磁共振图像划分成多个区域，再采用移动正方形法初始化活动轮廓模型，采用活动轮廓模型的能量最小化方程求解得到左心室内外表面轮廓，再由左心室内外表面轮廓重建三维轮廓，因采用混合高斯模型及活动轮廓模型求取的左心室内外表面轮廓较为准确，进而重建的三维轮廓的准确度也较高。

在另一个实施例中心脏磁共振图像的左心室三维轮廓重建方法，还包括：测量心室的泵血量和泵血率的步骤，根据该左心室三维轮廓采用点云法测量该左心室扩张末期的容积和左心室收缩末期的容积，再根据该左心室扩张末期的容积和左心室收缩末期的容积计算心室的泵血量和泵血率。

具体的，心室的泵血量为心脏扩张末期的血液容积与心室收缩末期的血液容积的差值。心室的泵血率为心脏扩张末期的血液容积与心室收缩末期的血液容积的差值与心脏扩张末期的血液容积的比值。

如图8所示为点云法测量体积的示意图。图8中，计算球体体积时，在原来球的周边立方体内等距取点阵，通过球内点阵和立方体内点阵中点的个数的比值，再乘以立方体体积可推算出球体的体积。

如图9所示，该测量心室的泵血量和泵血率的步骤包括：

步骤S902，左心室三维轮廓网格化及取点步骤，将该左心室三维轮廓以闭合网格表示，在该闭合网格内等距离取点，建立左心室三维轮廓网格。

将左心室三维轮廓以闭合网格表示。闭合网格内部的容积采用点云法进行计算。在该闭合网格内等距离取点，可将闭合网格分成多个相同的空间。如图10所示，为左心室三维轮廓划分的等距的小方格。

步骤S904，选取最小边界长方体的步骤，选取包含心室扩张末期的三维轮廓网格及与坐标轴平行的第一最小边界长方体和包含心室收缩末期的三维轮廓网格及与坐标轴平行的第二最小边界长方体。

步骤S906，统计点数量的步骤，获取第一最小边界长方体内点的数量、第二最小边界长方体内点的数量、心室扩张末期内点的数量和心室收缩末期内点的数量。

在该统计点数量的步骤之前，还包括：判断点是否在某一目标内部的步骤，获取穿过该点的直线，该直线的两端在该目标的最小边界长方体上，获取该直线在该点的两端与该目标的交点的数量，若交点的数量均为奇数，则该点在该目标内部，若交点的数量均为偶数，则该点在该目标外部，其中该目标为心室扩张末期或心室收缩末期。

直线与目标的交点即直线与扩张末期网格或心室收缩末期网格等网格的交点。直线与网格的交点可采用直线与网格面片的交点计算。可将网格的最小边界长方体等距划分为多个小立方体，建立两张对应表格，一张是每个小立方体对应的面片(该面片可为三角形或其他多边形)，一张是每个面片所在的小立方体。

求取每个面片所在的平面和直线的交点和判断每个面片所在平面和直线的交点是否在面片范围内过程为：

给定直线的两个端点为(x₁，y₁，z₁)和(x₂，y₂，z₂)，直线在三维空间的参数方程为：

$\begin{array}{c}x=x_1+t\×(x_2-x_1)\\ y=y_1+t\×(y_2-y_1)\\ z=z_1+t\×(z_2-z_1)\end{array}---\left(8\right)$

式(8)中，t∈[0,1]。

利用三个点可确定一个平面，平面方程可由平面上的一个顶点和平面的法向量确定。先由三角形的三个顶点V₁、V₂、V₃形成两个向量。法向量的方程是两个向量的叉乘，公式如式(9)：

V_n＝(V₂-V₁)×(V₃-V₁) (9)

平面方程是求出与平面内一点形成在平面内向量的所有点，与法向量垂直的向量都在平面内，平面方程为式(10)：

(V₂-V₁)·V_n＝0 (10)

将式(8)带入式(10)可求出交点。

检测一个交点是否在面片内可判断三角形和交点形成的三个夹角∠V₁VV₂、∠V₂VV₃和∠V₃VV₁之和是否等于360度。当∠V₁VV₂+∠V₂VV₃+∠V₃VV＝360°，交点在该面片内部，只有在面片内的点才是直线与网格的交点。

由于直线的两个端点在网格的最小边界长方体上，也就是在网格的外部，所以直线与网格若不相切，则总是有偶数个交点。故当直线在该点的两端与该网格的交点的数量为奇数，表示该点在网格内，交点的数量为偶数，表示该点在网格外。

如图11A、11B和11C所示，为穿过点的直线与网格的交点的示意图，三幅图中穿过点为112，交点为114。图11A和图11B中穿过点112的直线与网格的交点114在该点112两端的数量为一个，故该点112在网格内，图11C中穿过点112的直线与网格的交点114在该点112两端的数量为两个，故该点112在网格外。

步骤S908，计算容积的步骤，计算第一最小边界长方体的体积和第二最小边界长方体的体积，根据该心室扩张末期内点的数量及该第一最小边界长方体的体积和该第一最小边界长方体内点的数量计算得到心室扩张末期的容积，根据该心室收缩末期内点的数量及该第二最小边界长方体的体积和该第二最小边界长方体内点的数量计算得到该心室收缩末期的容积。

具体的，第一最小边界长方体和第二最小边界长方体的体积可根据对应的长宽高乘积进行求取。

心室扩张末期的容积等于该心室扩张末期内点的数量乘以该第一最小边界长方体的体积，再除以该第一最小边界长方体内点的数量。

同理，心室收缩末期的容积等于心室收缩末期内点的数量乘以该第二最小边界长方体的体积，再除以该第二最小边界长方体内点的数量。

步骤S910，计算泵血量和泵血率的步骤，根据该心室扩张末期的容积和心室收缩末期的容积计算得到该心室的泵血量和泵血率。

上述心脏磁共振图像的左心室三维轮廓重建方法，采用重建的左心室三维轮廓测量出泵血量和泵血率参数，测量方便且准确。

如图12所示，为一个实施例中心脏磁共振图像的左心室三维轮廓重建系统，包括模型建立模块120、初始化活动轮廓模块140、左心室内外轮廓获取模块160和三维轮廓重建模块180。其中：

模型建立模块120用于对获取的心脏磁共振图像建立混合高斯模型，采用预设数量高斯分布将该心脏磁共振图像划分成相应的灰度区域图像。

如图13所示，模型建立模块120包括先验概率获取子模块122、所属概率获取子模块124、更新子模块126、后验概率获取子模块128和循环计算子模块129。其中：

先验概率获取子模块122用于将心脏磁共振图像分成预设数量个高斯分布的灰度区域图像，建立混合高斯模型，计算混合高斯模型的先验概率。

因心脏磁共振图像的左心室内壁边界比较清晰可分，心脏磁共振图像按灰度大体上可分为三个区域，灰度最亮的区域为肺部的气管、心室内部的血液和脂肪，灰度比较暗的区域为肺部气泡区域，灰度处于中间的区域为心肌。心脏磁共振图像的灰度可采用三个高斯分布表示，用期望最大化可求解混合高斯模型的参数，该混合高斯模型的参数包括均值及方差。

本实施例中，将心脏磁共振图像分成三个区域，然后计算每个区域的先验概率。该先验概率可通过三种方法得到：第一种，对整个心脏磁共振图像的高斯拟合结果加随机扰动假设得到方差；第二种，采用K均值方法找到各区域中心，通过高斯拟合得到方差；第三种，通过经验抽样选取点作为区域中心，通过高斯拟合得到方差。该先验概率为采用三种方式根据以往数据分析得到的概率。

所属概率获取子模块124用于计算心脏磁共振图像中每个点属于某个区域的概率。

具体的，心脏磁共振图像的灰度为一个统计分布，灰度为x的点在区域内出现的概率为：

$p\left(x\right|i)=-\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_i}\exp (-\frac{|x-{\mathrm{\μ}}_i{|}^2}{2{\mathrm{\σ}}_i^2})---\left(1\right)$

其中，p(x|i)为灰度为x的点在i区域的概率，μ_i为在i区域的均值，σ_i ²为在i区域的方差，σ_i为标准差。

更新子模块126用于重新计算每个高斯分布的均值和方差。由所属概率获取子模块124重新计算每个点属于某个区域的概率。

后验概率获取子模块128用于根据每个点属于某个区域的概率采用贝叶斯概率计算后验概率，并将该后验概率作为下次计算先验概率的先验假设。

具体的，计算每个区域中各点的灰度均值公式为式(2)：

μ_i＝∑_xxp(x|i) (2)

计算每个区域中各点的方差为式(3)：

${\mathrm{\σ}}_i^2=E\left(\right|x-{\mathrm{\μ}}_i{|}^2)={\mathrm{\Σ}}_x|x-{\mathrm{\μ}}_i{|}^{2}p\left(x\right|i)---\left(3\right)$

采用贝叶斯公式计算后验概率，计算公式为式(4)：

$p\left(i\right|x)=\frac{p\left(x\right|i)P\left(i\right)}{{\mathrm{\Σ}}_i^{M}p\left(x\right|i)P\left(i\right)}---\left(4\right)$

其中，p(i)＝∑_xp(i|x)，M为总区域数，i表示第i个区域。

将后验概率作为下一轮先验概率的先验假设。

循环计算子模块129用于循环计算直到高斯模型的均值及方差变化量趋近于零，将各点分入后验概率最大值的区域，得到更新后的心脏磁共振图像的三个灰度区域图像。

具体的，通过循环计算直到高斯模型的均值及方差变化量趋近于零，将各点分入到后验概率最大值的区域。如此，可得到更新后的心脏磁共振图像的三个灰度区域图像，即将所有点分别标记对应的区域标识，如属于第一区域的点采用1标记，属于第二区域的点采用2标记，属于第三区域的点采用3标记。

初始化活动轮廓模块140用于采用移动正方形法对该灰度区域图像的边缘进行处理得到初始化的活动轮廓模型。

具体的，初始化活动轮廓模块140还用于将该灰度区域图像以预设移动正方形大小的网格划分，在该网格中插入短线段，将插入的短线段连接成曲线，取最长的圆形曲线初始化活动轮廓模型。

具体的，将左心室心肌内部用斜线表示，左心室心肌外部用空白表示的图像以网格划分，根据每个网格四个顶点的值，加入顶点在网格正方形边中点的短线段。该加入短线段的方法有多种，如图3所示，为四种加入短线段的情况。将插入的短线段连接成长度不同的曲线，取最长的圆形曲线初始化活动轮廓模型。如图4所示为移动正方形法初始化的心脏轮廓示意图。图5A为采用Canny边界检测方法检测的边界线条，图5B为边界提供给活动轮廓模型的受力情况。

左心室内外轮廓获取模块160用于根据该活动轮廓模型的能量最小化方程分别求解得到该左心室内外表面轮廓。

在一个实施例中，步骤S106包括：获取并根据活动轮廓模型的内力函数、图像力函数和交互作用力函数通过微积分建立活动轮廓模型的能量最小化方程，沿能量梯度方向移动，直到图像力函数与内力函数的值平衡，得到左心室内外表面轮廓。

具体的，活动轮廓模型的能量最小化方程如式(5)所示：

$E_{snake}={\∫}_0^1\left(E_{\mathrm{int}ernal}\right(v\left(s\right))+E_{image}\left(v\right(s\left)\right)+E_{constraint}\left(v\left(s\right)\right))ds---\left(5\right)$

式(5)中，E_snake为活动轮廓模型的能量；v(s)为活动轮廓上等距离点的坐标值；E_internal为活动轮廓内力函数，该内力受到曲线收缩和曲线平滑的限制；E_image为活动轮廓的图像力函数，该图像力函数是指活动轮廓模型移动到的位置的图像距图像边缘的图像灰度梯度力，将图像力函数作为活动轮廓模型的外力；E_constraint为交互作用力函数。图像力函数与内力函数的值平衡是指在外力和内力的作用下活动轮廓逐渐停止移动。

求解该活动轮廓的能量最小化方程，使得能量最小化，即轮廓的一阶导数、二阶导数和图像力的和最小。将能量最小化方程通过欧拉公式转化为式(6)：

$\begin{array}{c}{\mathrm{\αx}}_{ss}+{\mathrm{\βx}}_{ssss}+\∂E_{imag}/\∂x=0\\ {\mathrm{\αy}}_{ss}+{\mathrm{\βy}}_{ssss}+\∂E_{imag}/\∂y=0\end{array}---\left(6\right)$

式(6)中，α和β为多项式系数，x_ss和y_ss分别为x和y方向上的二次差分，x_ssss和y_ssss分别为x和y方向上的四次差分，是图像灰度梯度力。

对式(6)进行循环求解，将导数用邻近点的差分表示。活动轮廓点上的系数方程是一个由活动轮廓的一阶导数和二阶导数组成的矩阵在活动轮廓移动前计算得到的，在活动轮廓移动的过程中保持不变。如式(7)：

$\begin{array}{c}{x}_t={(A+\mathrm{\γ}I)}^{-1}({\mathrm{\γx}}_{t-1}-f_x(x_{t-1},y_{t-1}))\\ y_t={(A+\mathrm{\γ}I)}^{-1}({\mathrm{\γy}}_{t-1}-f_y(x_{t-1},y_{t-1}))\end{array}---\left(7\right)$

式(7)中，(f_x(x_t-1,y_t-1)，f_y(x_t-1,y_t-1))为(t-1)轮循环时的图像梯度向量，(x_t,y_t)为第t轮循环的计算活动轮廓模型每点的坐标值，(x_t-1,y_t-1)第(t-1)轮循环的计算活动轮廓模型每点的坐标值，A、γ、I均为系数。

图6A为活动轮廓模型分割的左心室中部磁共振图像，图6A中两圆形曲线为内外表面轮廓曲线；图6B为左心室尖端磁共振图像，图6B中两圆形曲线为内外表面轮廓曲线。

三维轮廓重建模块180用于获取一个心跳周期内多幅相邻的心脏磁共振图像的横轴图像，获取每幅横轴图像的左心室内外表面轮廓，按照时间序列将多个左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓。

因内外轮廓之间是心肌部分，内外壁都会在运动过程中随着心脏的收缩而收缩。获取一个心跳周期的磁共振图像的横轴图像，将右心室的横轴图像上的尖角作为重建三维网格的标识点，每个横轴图像上获取的左心室内外表面轮廓都以该标识点对齐，相邻横轴图像上的轮廓用三角面片连接起来，然后按照时间序列将多个左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓，如图7所示为左心室内外表面轮廓三维网格重建的结果。例如，每幅横轴图像的轮廓上取50个点，且每幅横轴图像的内或外表面轮廓由与左右心室间隔上右心室月牙形心后侧尖角标识点为内或外表面轮廓上的第一点组成三角网格。

上述心脏磁共振图像的左心室三维轮廓重建系统，采用混合高斯模型将磁共振图像划分成多个区域，再采用移动正方形法初始化活动轮廓模型，采用活动轮廓模型的能量最小化方程求解得到左心室内外表面轮廓，再由左心室内外表面轮廓重建三维轮廓，因采用混合高斯模型及活动轮廓模型求取的左心室内外表面轮廓较为准确，进而重建的三维轮廓的准确度也较高。

如图14所示，为另一个实施例中心脏磁共振图像的左心室三维轮廓重建系统，包括模型建立模块120、初始化活动轮廓模块140、左心室内外轮廓获取模块160和三维轮廓重建模块180，还包括参数获取模块190。其中：

参数获取模块190用于根据该左心室三维轮廓采用点云法测量该左心室扩张末期的容积和左心室收缩末期的容积，再根据该左心室扩张末期的容积和左心室收缩末期的容积计算心室的泵血量和泵血率。

具体的，心室的泵血量为心脏扩张末期的血液容积与心室收缩末期的血液容积的差值。心室的泵血率为心脏扩张末期的血液容积与心室收缩末期的血液容积的差值与心脏扩张末期的血液容积的比值。

如图15所示，该参数获取模块190包括网格化子模块192、选取子模块194、统计子模块196、容积计算子模块198和参数计算子模块199。其中：

网格化子模块192用于将该左心室三维轮廓以闭合网格表示，在该闭合网格内等距离取点，建立左心室三维轮廓网格。将左心室三维轮廓以闭合网格表示。闭合网格内部的体积采用点云法进行计算。在该闭合网格内等距离取点，可将闭合网格分成多个相同的空间。

选取子模块194用于选取包含心室扩张末期的三维轮廓网格及与坐标轴平行的第一最小边界长方体和包含心室收缩末期的三维轮廓网格及与坐标轴平行的第二最小边界长方体。

统计子模块196用于获取第一最小边界长方体内点的数量、第二最小边界长方体内点的数量、心室扩张末期内点的数量和心室收缩末期内点的数量。

容积计算子模块198用于计算第一最小边界长方体的体积和第二最小边界长方体的体积，根据该心室扩张末期内点的数量及该第一最小边界长方体的体积和该第一最小边界长方体内点的数量计算得到心室扩张末期的容积，根据该心室收缩末期内点的数量及该第二最小边界长方体的体积和该第二最小边界长方体内点的数量计算得到该心室收缩末期的容积。

参数计算子模块199用于据该心室扩张末期的容积和心室收缩末期的容积计算得到该心室的泵血量和泵血率。

该参数获取模块190还包括判断子模块。该判断子模块，用于获取穿过该点的直线，该直线的两端在该目标的最小边界长方体上，获取该直线在该点的两端与该目标的交点的数量，若交点的数量均为奇数，则该点在该目标内部，若交点的数量均为偶数，则该点在该目标外部，其中该目标为心室扩张末期或心室收缩末期。

直线与目标的交点即直线与扩张末期网格或心室收缩末期网格等网格的交点。直线与网格的交点可采用直线与网格面片的交点计算。可将网格的最小边界长方体等距划分为多个小立方体，建立两张对应表格，一张是每个小立方体对应的面片(该面片可为三角形或其他多边形)，一张是每个面片所在的小立方体。

求取每个面片所在的平面和直线的交点和判断每个面片所在平面和直线的交点是否在面片范围内过程为：

给定直线的两个端点为(x₁，y₁，z₁)和(x₂，y₂，z₂)，直线在三维空间的参数方程为：

$\begin{array}{c}x=x_1+t\×(x_2-x_1)\\ y=y_1+t\×(y_2-y_1)\\ z=z_1+t\×(z_2-z_1)\end{array}---\left(8\right)$

式(8)中，t∈[0,1]。

利用三个点可确定一个平面，平面方程可由平面上的一个顶点和平面的法向量确定。先由三角形的三个顶点V₁、V₂、V₃形成两个向量。法向量的方程是两个向量的叉乘，公式如式(9)：

V_n＝(V₂-V₁)×(V₃-V₁) (9)

平面方程是求出与平面内一点形成在平面内向量的所有点，与法向量垂直的向量都在平面内，平面方程为式(10)：

(V₂-V₁)·V_n＝0 (10)

将式(8)带入式(10)可求出交点。

检测一个交点是否在面片内可判断三角形和交点形成的三个夹角∠V₁VV₂、∠V₂VV₃和∠V₃VV₁之和是否等于360度。当∠V₁VV₂+∠V₂VV₃+∠V₃VV＝360°，交点在该面片内部，只有在面片内的点才是直线与网格的交点。

由于直线的两个端点在网格的最小边界长方体上，也就是在网格的外部，所以直线与网格若不相切，则总是有偶数个交点。故当直线在该点的两端与该网格的交点的数量为奇数，表示该点在网格内，交点的数量为偶数，表示该点在网格外。

如图11A、11B和11C所示，为穿过点的直线与网格的交点的示意图，三幅图中穿过点为112，交点为114。图11A和图11B中穿过点112的直线与网格的交点114在该点112两端的数量为一个，故该点112在网格内，图11C中穿过点112的直线与网格的交点114在该点112两端的数量为两个，故该点112在网格外。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (12)

1.一种心脏磁共振图像的左心室三维轮廓重建方法，包括：

建立心脏磁共振图像的混合高斯模型的步骤，对获取的心脏磁共振图像建立混合高斯模型，采用预设数量高斯分布将所述心脏磁共振图像划分成相应的灰度区域图像；

初始化活动轮廓模型的步骤，采用移动正方形法对所述灰度区域图像的边缘进行处理得到初始化的活动轮廓模型；

确定左心室内外表面轮廓的步骤，根据所述活动轮廓模型的能量最小化方程分别求解得到所述左心室内外表面轮廓，包括获取并根据活动轮廓模型的内力函数、图像力函数和交互作用力函数通过微积分建立活动轮廓模型的能量最小化方程，沿能量梯度方向移动，直到图像力函数与内力函数的值平衡，得到左心室内外表面轮廓，图像力函数是指活动轮廓模型移动到的位置的图像距图像边缘的图像灰度梯度力；以及

重建左心室三维轮廓的步骤，获取一个心跳周期内多幅空间位置相邻的心脏磁共振图像的横轴图像，获取每幅横轴图像的左心室内外表面轮廓，将多个横轴图像的左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓，包括：获取一个心跳周期内多幅空间位置相邻的磁共振图像的横轴图像，将右心室的横轴图像上的尖角作为重建三维网格的标识点，每个横轴图像上获取的左心室内外表面轮廓都以该标识点对齐，相邻横轴图像上的轮廓用三角面片连接起来，然后按照时间序列将多个横轴图像的左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓，每幅横轴图像的内或外表面轮廓由与左右心室间隔上右心室月牙形心后侧尖角标识点为内或外表面轮廓上的第一点组成三角网格。

2.根据权利要求1所述的心脏磁共振图像的左心室三维轮廓重建方法，其特征在于，所述方法在重建左心室三维轮廓的步骤之后还包括：

测量心室的泵血量和泵血率的步骤，根据所述左心室三维轮廓采用点云法测量左心室扩张末期的容积和左心室收缩末期的容积，再根据所述左心室扩张末期的容积和左心室收缩末期的容积计算心室的泵血量和泵血率。

3.根据权利要求2所述的心脏磁共振图像的左心室三维轮廓重建方法，其特征在于，所述测量心室的泵血量和泵血率的步骤包括：

左心室三维轮廓网格化及取点步骤，将所述左心室三维轮廓以闭合网格表示，在所述闭合网格内等距离取点，建立左心室三维轮廓网格；

选取最小边界长方体的步骤，选取包含心室扩张末期的三维轮廓网格及与坐标轴平行的第一最小边界长方体和包含心室收缩末期的三维轮廓网格及与坐标轴平行的第二最小边界长方体；

统计点数量的步骤，获取第一最小边界长方体内点的数量、第二最小边界长方体内点的数量、心室扩张末期内点的数量和心室收缩末期内点的数量；

计算容积的步骤，计算第一最小边界长方体的体积和第二最小边界长方体的体积，根据所述心室扩张末期内点的数量及所述第一最小边界长方体的体积和所述第一最小边界长方体内点的数量计算得到心室扩张末期的容积，根据所述心室收缩末期内点的数量及所述第二最小边界长方体的容积和所述第二最小边界长方体内点的数量计算得到所述心室收缩末期的容积；以及

计算泵血量和泵血率的步骤，据所述心室扩张末期的容积和心室收缩末期的容积计算得到所述心室的泵血量和泵血率。

4.根据权利要求3所述的心脏磁共振图像的左心室三维轮廓重建方法，其特征在于，在所述统计点数量的步骤之前，还包括：

判断点是否在某一目标内部的步骤，获取穿过所述点的直线，所述直线的两端在所述目标的最小边界长方体上，获取所述直线在所述点的两端与所述目标的交点的数量，若交点的数量均为奇数，则所述点在所述目标内部，若交点的数量均为偶数，则所述点在所述目标外部，其中所述目标为心室扩张末期或心室收缩末期。

5.根据权利要求1所述的心脏磁共振图像的左心室三维轮廓重建方法，其特征在于，所述建立心脏磁共振图像的混合高斯模型的步骤包括：

将心脏磁共振图像分成预设数量个高斯分布的灰度区域图像，建立混合高斯模型，计算混合高斯模型的先验概率；

计算心脏磁共振图像中每个点属于某个区域的概率；

重新计算每个高斯分布的均值和方差，再重新计算每个点属于某个区域的概率；

根据每个点属于某个区域的概率采用贝叶斯概率计算后验概率，并将所述后验概率作为下次计算先验概率的先验假设；以及

循环计算直到高斯模型的均值及方差变化量趋近于零，将各点分入后验概率最大值的区域，得到更新后的心脏磁共振图像的预设数量个灰度区域图像。

6.根据权利要求1所述的心脏磁共振图像的左心室三维轮廓重建方法，其特征在于，所述初始化活动轮廓模型的步骤包括：

将所述灰度区域图像以预设移动正方形大小的网格划分，在所述网格中插入短线段，将插入的短线段连接成曲线，取最长的圆形曲线初始化活动轮廓模型。

7.一种心脏磁共振图像的左心室三维轮廓重建系统，其特征在于，包括：

模型建立模块，用于对获取的心脏磁共振图像建立混合高斯模型，采用预设数量高斯分布将所述心脏磁共振图像划分成相应的灰度区域图像；

初始化活动轮廓模块，用于采用移动正方形法对所述灰度区域图像的边缘进行处理得到初始化的活动轮廓模型；

左心室内外轮廓获取模块，用于根据所述活动轮廓模型的能量最小化方程分别求解得到所述左心室内外表面轮廓，包括：获取并根据活动轮廓模型的内力函数、图像力函数和交互作用力函数通过微积分建立活动轮廓模型的能量最小化方程，沿能量梯度方向移动，直到图像力函数与内力函数的值平衡，得到左心室内外表面轮廓，图像力函数是指活动轮廓模型移动到的位置的图像距图像边缘的图像灰度梯度力；以及

三维轮廓重建模块，用于获取一个心跳周期内多幅空间位置相邻的心脏磁共振图像的横轴图像，获取每幅横轴图像的左心室内外表面轮廓，将多个横轴图形的左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓，包括：获取一个心跳周期内多幅空间位置相邻的磁共振图像的横轴图像，将右心室的横轴图像上的尖角作为重建三维网格的标识点，每个横轴图像上获取的左心室内外表面轮廓都以该标识点对齐，相邻横轴图像上的轮廓用三角面片连接起来，然后按照时间序列将多个横轴图像的左心室内外表面轮廓连接起来，形成心脏磁共振图像的左心室三维轮廓，每幅横轴图像的内或外表面轮廓由与左右心室间隔上右心室月牙形心后侧尖角标识点为内或外表面轮廓上的第一点组成三角网格。

8.根据权利要求7所述的心脏磁共振图像的左心室三维轮廓重建系统，其特征在于，所述系统还包括：

参数获取模块，用于根据所述左心室三维轮廓采用点云法测量左心室扩张末期的容积和左心室收缩末期的容积，再根据所述左心室扩张末期的容积和左心室收缩末期的容积计算心室的泵血量和泵血率。

9.根据权利要求8所述的心脏磁共振图像的左心室三维轮廓重建系统，其特征在于，所述参数获取模块包括：

网格化子模块，用于将所述左心室三维轮廓以闭合网格表示，在所述闭合网格内等距离取点，建立左心室三维轮廓网格；

选取子模块，用于选取包含心室扩张末期的三维轮廓网格及与坐标轴平行的第一最小边界长方体和包含心室收缩末期的三维轮廓网格及与坐标轴平行的第二最小边界长方体；

统计子模块，用于获取第一最小边界长方体内点的数量、第二最小边界长方体内点的数量、心室扩张末期内点的数量和心室收缩末期内点的数量；

容积计算子模块，用于计算第一最小边界长方体的容积和第二最小边界长方体的容积，根据所述心室扩张末期内点的数量及所述第一最小边界长方体的容积和所述第一最小边界长方体内点的数量计算得到心室扩张末期的容积，根据所述心室收缩末期内点的数量及所述第二最小边界长方体的容积和所述第二最小边界长方体内点的数量计算得到所述心室收缩末期的容积；以及

参数计算子模块，用于据所述心室扩张末期的容积和心室收缩末期的容积计算得到所述心室的泵血量和泵血率。

10.根据权利要求9所述的心脏磁共振图像的左心室三维轮廓重建系统，其特征在于，所述参数获取模块还包括：

判断子模块，用于判断点是否在某一目标内部，获取穿过所述点的直线，所述直线的两端在所述目标的最小边界长方体上，获取所述直线在所述点的两端与所述目标的交点的数量，若交点的数量均为奇数，则所述点在所述目标内部，若交点的数量均为偶数，则所述点在所述目标外部，其中所述目标为心室扩张末期或心室收缩末期。

11.根据权利要求7所述的心脏磁共振图像的左心室三维轮廓重建系统，其特征在于，所述模型建立模块包括：

先验概率获取子模块，用于将心脏磁共振图像分成预设数量个高斯分布的灰度区域图像，建立混合高斯模型，计算混合高斯模型的先验概率；

所属概率获取子模块，用于计算心脏磁共振图像中每个点属于某个区域的概率；

更新子模块，用于重新计算每个高斯分布的均值和方差，再由所述所属概率获取子模块重新计算每个点属于某个区域的概率；

后验概率获取子模块，用于根据每个点属于某个区域的概率采用贝叶斯概率计算后验概率，并将所述后验概率作为下次计算先验概率的先验假设；以及

循环计算子模块，用于循环计算直到高斯模型的均值及方差变化量趋近于零，将各点分入后验概率最大值的区域，得到更新后的心脏磁共振图像的预设数量个灰度区域图像。

12.根据权利要求7所述的心脏磁共振图像的左心室三维轮廓重建系统，其特征在于，所述初始化活动轮廓模块还用于将所述灰度区域图像以预设移动正方形大小的网格划分，在所述网格中插入短线段，将插入的短线段连接成曲线，取最长的圆形曲线初始化活动轮廓模型。