一种基于心脏磁共振电影短轴图像的右心室多图谱分割方法与流程

本发明属于磁共振图像处理领域，具体涉及一种基于心脏磁共振电影短轴图像的右心室多图谱分割方法。

背景技术：

心脏磁共振电影短轴成像除了具有无电离辐射损伤、多方位成像、软组织对比度高等优点外，还能提供时空分辨率高的动态电影图像。心脏磁共振短轴电影图像一方面能够精确的显示心脏的结构，另一方面能够用于心脏分割和量化心脏功能参数。

右心室因其心肌壁薄，边界不清，心肌结构变化大，与周围组织对比度低等特点，成为心脏分割的一大难点，精确分割右心室更是全球心脏分割领域的一大挑战。

传统医学图像分割算法对心脏右心室的分割均存在不足：主动轮廓模型收敛条件难以确定、水平集算法的初始轮廓非常敏感。传统多图谱算法融合阶段多使用加权融合策略或staple融合策略进行处理，在对脑部海马体或肝脏分割是能获得较好的结果，但对右心室分割时不能达到良好的效果。

技术实现要素：

为解决右心室分割的技术问题，本发明的目的在于提供一种基于心脏磁共振电影短轴图像的右心室多图谱分割方法。

本发明是通过以下技术方案实现的：

一种基于心脏磁共振电影短轴图像的右心室多图谱分割方法，包括以下步骤：

s1、利用磁共振成像系统采集被试人一定数量的心脏原始磁共振图像，对所述原始磁共振图像进行预处理，提取感兴趣区域；

s2、在所述原始磁共振图像中选取固定数量的右心室的图谱图像加入图谱集，专家手动分割所述图谱图像获得专家手动分割结果；

s3、将所述图谱图像与目标图像进行配准，所述配准采用基于归一化互信息的b样条变换，将变换参数应用于与所述图谱图像对应的所述标记图像上，获得述配准后的图谱图像，得到右心室粗分割结果；

s4、对所述粗分割结果采用collate融合，先对完整数据进行对数似然估计，然后再利用最大期望算法迭代求解直至收敛，并进行修正处理获得右心室最终分割结果；

s5、重复步骤s1～s4，计算右心室容积和不同生理指标，选取相关参数作为客观的评价指标对实验结果进行评价和分析。

进一步的，所述提取感兴趣区域具体方法为：利用霍夫圆变换确定左心室中心位置，通过左右心室的在短轴图像上的位置关系确定右心室位置并截取同时包含左右心室的部分图像作为感兴趣区域。

进一步的，所述对数似然估计采用公式：

θ＝argmaxθlnf(d，t，c|θ)

其中，θ表示所述粗分割结果对应的性能水平估计，f(d，t，c|θ)表示右心室完整数据对应的随机向量的概率质量分布函数，d是由所述粗分割结果和所述配准后的图谱图像的像素构成的矩阵，t表示所述像素的未知的真实分割结果，c表示每个所述像素的一致性或混乱性。

进一步的，所述最大期望算法迭代求解采用公式：

q(θ|θ⁽⁰⁾)≡e[lnf(d，t，c|θ)|d，θ⁽⁰⁾＝∑tf(d，t，c|θ)f(d，t，c|θ⁽⁰⁾)

其中，q(θ|θ⁽⁰⁾)表示θ＝θ⁽⁰⁾时所述完整数据对数似然函数的条件期望，θ(0)表示所述粗分割结果对应的性能水平估计的初始值。

进一步的，所述修正处理方法为：对右心室最顶部易出现分割错误的的1～2层进行结果修正，中部图像结果中心坐标在顶层初始分割结果内时，进行形状约束，有效防止泄露；中部图像结果中点坐标不在顶层图像初始分割结果内时，进行手动分割。

进一步的，所述生理指标包括舒张末期容积edv、收缩末期容积esv、射血分数ef和每搏输出量sv。

进一步的，所述相关参数包括精确度、相关性和一致性。

更进一步的，所述精确度指标包括豪斯多夫距离hd，采用公式：

hd(a，b)＝max(maxa∈a(minb∈bd(a，b))，maxb∈b(mina∈ad(a，b)))

其中，a表示所述专家手动分割结果，b表示所述最终分割结果，d(a,b)为欧几里得距离，其中a是所述图谱图像中的点，b是所述目标图像中的点，hd(a，b)表示所述图谱图像和所述目标图像的最大不对称差异，用于反映右心室所述专家手动分割结果与本发明的所述最终分割结果的距离差异。

更进一步的，所述精确度指标包括矩阵相似度dm，采用公式：

其中，aa表示所述最终分割结果的体积大小，am表示所述专家手动分割结果的体积大小，a∩＝aa∩am表示重叠区域的体积大小，dm(aa，am)表示所述最终分割结果与所述专家手动分割结果的重合度。

相比现有技术,本发明具有如下有益效果：

本发明提供一种基于心脏磁共振电影短轴图像的右心室多图谱分割方法，采用多图谱框架下collate融合算法对一系列右心室粗分割结果进行融合。本发明在staple算法的基础上添加了右心室图像中像素一致性概率，较staple算法具有更高的鲁棒性，并能够提高融合的准确度和精确性，可以精确分割心脏右心室短轴图像。

附图说明

图1为右心室舒张末期短轴底部图像分割结果图；

图2为右心室舒张末期短轴中部图像分割结果图；

图3为右心室舒张末期短轴顶部图像分割结果图；

图4为右心室收缩末期短轴底部图像分割结果图；

图5为右心室收缩末期短轴中部图像分割结果图；

图6为右心室收缩末期短轴顶部图像分割结果图；

其中，a表示原始图像中的感兴趣区域，包含完整的左、右心室；b表示专家手动分割结果，即金标准；c表示本发明方法得到的分割结果；d表示本发明分割结果在感兴趣区域中的标记；

图7为右心室最终分割结果与专家手动分割结果的距离差异对比图；

图8为右心室最终分割结果与专家手动分割结果的相似度对比图。

具体实施方式

下面对本发明的实施例作详细说明，本实施例以本发明的技术方案为依据开展，给出了详细的实施方式和具体的操作过程。

本发明实施例用不同时相不同部位的多幅心脏磁共振短轴电影图像进行右心室分割，得到最终分割结果的具体实施过程。其中，用于右心室分割的心脏磁共振电影短轴图像数据来源于磁共振系统，经ssfp序列获得。实验数据中，男性7例，女性3例，年龄涵盖14到75岁。具体成像参数：图像大小256×256，层厚6-8mm，层间距2-4mm，每组数据包含6-10层，每层20-28个时相，包含多个心动周期。

本发明实施例包括以下步骤：

s1、利用磁共振成像系统采集被试人一定数量的心脏原始磁共振图像，对图像进行预处理。首先利用霍夫圆变换确定左心室中心位置，然后通过左右心室的在短轴图像上的位置关系确定右心室位置并截取同时包含左右心室的部分图像作为感兴趣区域。

s2、在所述原始磁共振图像中选取结构清晰，受周围组织影响较小的右心室短轴图像加入图谱集，专家手动分割图谱图像获得右心室分割结果，并生成相应的标记图像l1-ln。

s3、每张图谱图像具有对应的已分割好的标记图像，配准过程采用基于归一化互信息的b样条变换，将图谱图像分别配准到目标图像后，记录对应的变换参数q1-qn，将变换参数应用于与图谱图像对应的标记图像上，得到右心室粗分割结果l′1-l′n。

s4、对所述粗分割结果采用collate融合，先对完整数据进行对数似然估计，然后再利用最大期望算法迭代求解直至收敛，并进行修正处理获得右心室最终分割结果。

定义右心室配准后的图谱图像中有n个像素点，对应每个像素点有r个粗分割结果，像素用i表示，粗分割结果用j表示。一系列标签l表示某种粗分割结果将一系列可能的标准分配给所有的像素n。d是一个n*r的矩阵，用于r个粗分割结果对所有像素n的描述，其中dij∈{0，1，...l-1}。长度为n的向量t，代表所有像素n的未知的真实分割结果，其中ti∈{0，1，...l-1}。

在staple算法的基础上，定义一个具有n个元素的向量c，表示右心室图像中每个像素的一致性或混乱性，共有f个一致可能性等级。该向量中所有元素ci∈{0，1，...f-1}。其中ci＝0表示像素i具有混乱性，ci＝1表示像素i具有一致性。一致性和混乱性是从相反的方向描述同一个现象，像素i一致性减少，混乱性就增加。

定义θ表示r个右心室粗分割结果的特性，每一个像素θj是一个l*l的混合矩阵，矩阵中的每一个元素量化了右心室粗分割j中某一像素是真实分割结果的可能概率。作为参考，“金标准”是一个单位矩阵。将右心室完整数据定义为(d，t，c)，完整数据的概率群分布定义为f(d，t，c|θ)。

对完整数据的对数似然估计表示成公式：

θ＝argmaxθlnf(d，t，c|θ)

其中，θ＝[θ1，θ2，...，θr]表示与r个右心室粗分割结果对应的性能水平估计；f(d，t，c|θ)表示右心室完整数据对应的随机向量的概率质量分布函数。

e步计算右心室完整数据的对数似然函数的条件期望表示成公式：

q(θ|θ⁽⁰⁾)≡e[lnf(d，t，c|θ)|d，θ⁽⁰⁾]＝σtf(d，t，c|θ)f(d，t，c|θ⁽⁰⁾)

其中，q(θ|θ⁽⁰⁾)表示θ＝θ⁽⁰⁾时右心室完整数据对数似然函数的条件期望。θ⁽⁰⁾表示参数θ的初始值，即右心室粗分割结果对应的性能水平估计的初始值。

m步要求在θ的参数空间中求解q(θ|θ⁽⁰⁾)的最大值，即对于所有θ⁽⁰⁾，

存在θ⁽¹⁾满足公式：

q(θ|θ⁽¹⁾)≥q(θ|θ⁽⁰⁾)

其中，q(θ|θ⁽¹⁾)和q(θ|θ⁽⁰⁾)分别表示θ＝θ⁽¹⁾和θ＝θ⁽⁰⁾时右心室完整数据对数似然函数的条件期望。θ⁽¹⁾的估计值是通过q(θ|θ⁽⁰⁾)的最大化过程确定的。

本发明collate融合策略主要利用最大期望算法，分为e步和m步。e步计算完整数据的对数似然函数的条件期望，m步求解函数最大化时θ的值。每次迭代过程均由e步和m步交替重复进行，直至收敛为止。迭代收敛过程控制在20次左右，融合效果较好且用时较短，经实验确定截止参数ε＝1*10^-7。通过collate融合算法对右心室粗分割结果进行融合得到右心室最终分割结果。

对于结构复杂且不易分割的右心室最顶部1-2层图像，需要进行结果修正。以其相邻的中部图像结果作为模板，以模板中心点坐标是否在顶层初始分割结果内分为两种情况进行处理：中部图像结果中心坐标在顶层初始分割结果内时，进行形状约束，有效防止泄露；中部图像结果中点坐标不在顶层图像初始分割结果内时，进行手动分割，由于顶层数据面积小，层数少，手动分割在提高精度的同时对处理时间影响不大。

s5、重复以上步骤s1～s4，首先获得右心室某一时相内所有层面扫描图像的分割结果，由面积长度法计算每层的心室面积，再计算每层的心室面积与层厚的乘积，通过累加得到容积，如舒张末期容积edv和收缩末期容积esv。通过计算所有时相的容积可以绘制右心室充盈曲线。通过对右心室舒张末期和收缩末期图像的分割，可以获得每搏输出量和射血分数等生理信息。通过客观的评价指标：精确度、相关性、一致性等对实验结果进行评价和分析，计算得到精确度指标：豪斯多夫距离和矩阵相似度。

如图7为右心室舒张末期和收缩末期两状态下的最终分割结果与专家手动分割结果的距离差异对比图。横轴代表时间点，纵轴代表最终分割结果与专家手动分割结果的豪斯多夫距离hd。

如图8为右心室舒张末期和收缩末期两状态下的最终分割结果与专家手动分割结果的相似度对比图。横轴代表时间点，纵轴代表最终分割结果与专家手动分割结果的相似度。dm(aa，am)表示右心室最终分割结果与专家手动分割结果的相似度，变化范围从0(完全不匹配)到1(完全匹配)，重叠率越高，结果越精确。

以上实施例为本申请的优选实施例，本领域的普通技术人员还可以在此基础上进行各种变换或改进，在不脱离本申请总的构思的前提下，这些变换或改进都应当属于本申请要求保护的范围之内。