一种心脏灌注磁共振图像的处理方法与流程

本发明涉及医学图像领域，尤其是涉及心脏灌注磁共振图像的后处理方法。

背景技术：

近年来，心血管疾病的发病率和死亡率正逐年增加。全世界每年几千万人死于心血管疾病。心肌缺血是引发心血管疾病发病和死亡的最重要原因。心脏灌注磁共振成像被认为是检查缺血性心脏病的首选无创型检查方法。根据心脏灌注磁共振图像的基本原理，放射科医生可以通过直接观察左心室心肌部分的信号强度随时间的变化情况进行心肌缺血的临床诊断。但是，在图像获取过程中，心脏运动、呼吸及病人位置移动在心脏灌注磁共振图像中产生的形变和运动伪影、快速时间采样、及相对较低的图像对比度，限制了直接进行视觉诊断的可靠性和效率。

计算机辅助量化分析方案有助于提高放射科医生基于四维心肌灌注磁共振图像对心肌缺血诊断的可靠性和效率。在研究计算机辅助量化分析方案中，时间序列图像的配准和左心室心肌相关的心内外膜分割是重要且具有挑战性的工作，其结果将直接影响左心室心肌部分信号强度随时间变化曲线的相关分析和可视化结果，进而影响心肌缺血诊断的准确性和效率。

应用计算机辅助量化分析方案首先要解决的是由心脏运动、呼吸及病人位置移动在心肌灌注磁共振图像中产生的形变。目前已有手动、半自动的方法对 心肌灌注磁共振图像进行形变校正。

经过形变校正后，需要对心肌相关的心内外膜进行分割。在四维心肌灌注磁共振图像成像过程中，为了能够捕捉到造影剂首次流入心脏各个腔室的时刻，就要保证较高的时间采样率，代价则是图像的空间分辨率、信噪比和容积效应，另外心肌周围复杂的组织结构(如：肺、气管、胃等)，这些都增加了心肌相关的心内外膜分割的难度。其中，手动和需用户辅助的半自动检测方法比较繁复、效率低、并且存在较大的观察者之间和观察者自身的差异。

技术实现要素：

本发明所要解决的技术问题是提供一种心脏灌注磁共振图像的全自动化量化分析方法。

本发明为解决上述技术问题而采用的技术方案是：一种心脏灌注磁共振图像的处理方法，其特征在于包括以下步骤：

获取包含左心室心肌的若干层切片的磁共振图像I_NT，其中N表示在同一心跳周期中切片所在层的序号，T表示不同心跳周期的序号，N、T均为大于或等于1的整数；

从起始层开始，分别分割N层切片的磁共振图像在T个心跳周期内的心肌内膜，获得每一层切片对应的磁共振图像I_MT的参考图像I_Mr；

以每层切片的磁共振图像I_Mr中的参考图像I_Mr为基准，完成本层切片在T个心跳周期内的磁共振图像的配准；以每层切片的磁共振图像的心肌内膜的分割结果为依据，完成相对应层心肌外膜的分割；

计算左心室腔区域和心肌区域的平均灰度-时间曲线，并平滑处理；以心肌内膜所界定的血池区域的质心点和右心室边界的中心点连线为基线，按照牛眼图划分方法，将心肌分块。

优选的，还包括在所述平均灰度-时间曲线上计算造影剂(药物)开始流入左心室的心跳周期(序号)，在所述平均灰度-时间曲线上计算造影剂(药物)完全流出左心室的心跳周期(序号)。

优选的，所述心肌外膜分割包括下步骤：

计算每层切片在T个心跳周期内的磁共振图像I_MT的最大强度投影图I_MIP(M)，

以当前层图像已分割的心肌内膜所界定的血池区域质心为中心，将当前层参考图像I_Mr和最大投影图I_MIP(M)分别转换为极坐标图像P_Mr和P_MIP(M)；

在极坐标图像P_Mr中，寻找右心室区域和禁区，在极坐标图像P_MIP(M)中应用动态规划完成心肌外膜的分割，并将分割结果变换至直角坐标系中；

优选的，所述右心室区域和禁区通过以下方式设定：

对极坐标图像P_Mr进行模糊C均值聚类；

若当前层为起始层，在最亮的类中，选取除血池区域外最大的区域为右心室区域；若当前层不是起始层，则在最亮的类中，选取与前一层右心室区域重叠最大的区域为当前层的右心室区域；

将血池区域膨胀三个像素后的区域和最亮的类中除去血池和右心室外区域及其以下像素均设为禁区。

优选的，所述聚类数目分别与当前层在进行心肌内膜分割中寻找血池区域时采用的聚类数目相同。

优选的，按照牛眼图划分方法，所述心脏的心尖层与中间层及心底层处的心肌被划分为若干块。

优选的，所述心脏的心尖层处的心肌被划分的块数小于中间层及心底层处的心肌被划分的块数。

优选的，根据心肌的每一块计算平均灰度-时间曲线，平滑并基线校正， 根据曲线计算峰值、达峰时间和最大上升梯度；绘制参数牛眼图。

优选的，对配准后图像中心肌的每一个像素点计算灰度-时间曲线，平滑并基线校正，根据曲线计算峰值、达峰时间和最大上升梯度，在参考图像中用颜色图映射这些参数。

优选的，将心肌平分成心肌内圈和心肌外圈，对心肌内圈和外圈分别分块。

本发明对比现有技术有如下的有益效果：本发明专利可自动确定参考相位图像、感兴趣区域和左心室心肌相关的心内膜，结合采用刚体和非刚体配准方法实现四维心肌灌注磁共振图像在时序上的全局和局部配准，并采用动态规划技术实现左心室心肌相关的心外膜的分割，进而实现左心室心肌部分信号强度随时间变化曲线的可视化及全自动量化分析。

【附图说明】

图1为本发明实施例中获取若干层切片在若干个心跳周期内的心脏磁共振图像的示意图；

图2为本发明实施例中心脏灌注磁共振图像分割过程所对应图像；

图3为本发明实施例中心脏灌注磁共振图像的配准前、配准后的图像；

图4为在参考图像中进行分块后得到结果图；

图5为左心室腔区域和心肌区域的平均灰度-时间曲线；

图6为本发明实施例中心脏灌注磁共振图像的处理方法流程图。

【具体实施方式】

下面结合附图和实施例对本发明作进一步的描述。

一种心脏灌注磁共振图像的处理方法，其包括以下主要步骤：

获取包含左心室心肌的若干层切片的磁共振图像I_NT，其中N表示在同一心 跳周期中切片所在层的序号，T表示不同心跳周期的序号，N、T均为大于或等于1的整数；

从起始层开始，分别分割N层切片的磁共振图像在T个心跳周期内的心肌内膜，获得每一层切片对应的磁共振图像I_MT的参考图像I_Mr，其中M为大于或等于1，小于或等于N的整数；

以每层切片的磁共振图像I_Mr中的参考图像I_Mr为基准，完成本层切片在T个心跳周期内的磁共振图像的配准；以每层切片的磁共振图像的心肌内膜的分割结果为依据，完成相对应层心肌外膜的分割；

计算左心室腔区域和心肌区域的平均灰度-时间曲线，并平滑处理；以心肌内膜所界定的血池区域的质心点和右心室边界的中心点连线为基线，将心肌分块。

本发明主要关注心脏灌注磁共振图像的处理方法，即心脏灌注磁共振图像全自动量化分析方法，以下为主要步骤的具体说明。

【分割心肌内膜】

一种心脏灌注磁共振图像的心肌内膜的分割方法，其包括以下步骤：

获取包含左心室心肌的若干层切片的磁共振图像I_NT，其中N表示在同一心跳周期中切片所在层的序号(数)，T表示不同心跳周期的序号(数)，N、T均为大于或等于1的整数；对N层切片中的每一相应的层(M层)的切片的磁共振图像I_MT进行以下处理：

选定起始层(当M＝1时)切片在T个心跳周期内的磁共振图像I_1T，分割磁共振图像I_1T中的心肌内膜，确定第r个心跳周期所对应的磁共振图像为参考图像I_1r以及感兴趣区域ROI₁；其中，r为小于T的整数；

分别选定第二层至第N层切片在T个心跳周期内的磁共振图像I_MT，设定第 M层(当前层)切片在T个心跳周期内的磁共振图像为I_MT，以第M-1层(其前一层)切片在T个心跳周期内的磁共振图像中选择出的参考图像I_(M-1)r的心跳周期的序号r，感兴趣区域ROI_M-1为参考依据，分割第M层切片磁共振图像I_MT中的心肌内膜，获得第M层切片I_MT磁共振图像中的参考图像I_Mr；其中，M的取值为大于1且小于等于N的整数；

以每层切片的磁共振图像中的参考图像为基准，完成本层切片在T个心跳周期内的磁共振图像的配准。

其中，所述磁共振图像I_NT可通过以下方式获得，首先，向被扫描对象注射示踪剂(或其他药物)，在对比剂到达心脏之前、之中、之后的时间，用磁共振成像设备获得心脏不同部位的切片(slice)的磁共振图像。本实施例中，沿心脏的长轴方向(大致为上下方向)，获得6个不同层面的切片图像，即N的取值为6，但是，可根据具体的需求，对N取不同的值，例如N的值也可以为4、8、10等。另外，需要在不同个连续的心跳周期T或心脏时相(phase)获得每一层面的切片磁共振图形，T的取值范围为40-60，假定每个心跳周期的时间为Δt，则完成磁共振扫描所需的总时间大约为Δt*T。在每个心跳周期内，需要完成一次N个层面的切片的磁共振图像，因此，一共需要有N*T幅磁共振图像。

其中，起始层切片位于左心室中间层以下，特别的所述起始层切片优选为心底层，此处由心肌内膜所界定(包围)的血池(blood pool)最大。

其中，分割起始层(定义M＝1时的层为起始层)切片在不同心跳周期的磁共振图像I_1T中的心肌内膜，包括以下步骤：

a)从心跳周期总数的1/4(T/4)处开始到2/3(2*T/3)处结束，以每间隔2个心跳周期所对应的磁共振图像I_1T被选定为候选图像，并以候选图像中心确 定初始感兴趣区域ROI₁’(参图2a)，在初始感兴趣区域ROI₁’中做模糊C均值聚类；优选的，所述初始感兴趣区域ROI₀的形状为正方形，所述初始感兴趣区域ROI₀边长为111像素，设定聚类数目为2。

b)在每一个候选图像对应的初始感兴趣区域ROI₁’的聚类二值图中选取圆度最大的区域，此区域被定义为候选图像的心肌内膜(参图2c)所界定的血池区域；圆度的定义为(周长*周长)/(4*PI*面积)；

c)根据所有血池区域的质心点位置、平均灰度和圆度，选择出起始层切片的磁共振图像的血池区域(参图2b)；其中，所述质心点位置越远，平均灰度越大且圆度越大的区域为起始层血池区域的概率越大，且以该血池区域所在的图像作为参考图像；

d)对血池区域求凸包并平滑，其中心点和长轴可用于设置最终的感兴趣区域ROI₁的位置和大小。

进一步的，分割第二层至第N层切片在T个心跳周期内的磁共振图像I_MT中的心肌内膜，包括以下步骤：

a)从获取第M层切片(当前层)的磁共振图像所需的心跳周期总数的1/4处(T/4)开始到第M-1(上一层)层切片磁共振图像中的参考图像对应的心跳周期总数(序数)加3处结束，每间隔2个心跳周期所对应的磁共振图像被选定为第M层切片的候选图像，第M层切片的候选图像分别以第M-1层切片所对应的参考图像的心肌内膜所界定的血池的中心和长轴加若干个像素作为第M层切片所对应的候选图像的心肌内膜所界定的血池的初始感兴趣区域ROI_M’的中心和边长，在初始感兴趣区域中做模糊C均值聚类，目前设定聚类数目为2类；所述像素数为5个至20个，优选为10个。

b)在每一个第M层切片的候选图像的聚类二值图中选取与第M-1层切片的 参考图像的血池重合度最大、且长轴小于第M-1层切片的参考图像的血池长轴的1.1倍的区域被定义为相对应的候选图像血池区域；如果没找到相对应的候选图像的血池，聚类数目加1，继续做模糊聚类，直到出现合适的候选图像血池区域；

c)根据所有第M层切片的候选图像的血池区域的质心点位置、平均灰度和圆度，即质心点位置越远，平均灰度越大且圆度越大的区域为血池区域的概率越大，选择出第M层切片的候选图像的血池区域，且该候选图像为第M层切片的参考图像I_Mr。

d)对第M层切片的候选图像血池区域求凸包并平滑，其中心点和长轴可用于设置最终感兴趣区域的位置和大小，具体为以血池区域质心为中心点，血池长轴加上2*心肌厚度为感兴趣区域的边长，可选的，心肌厚度范围为6-20mm。

【时序图像配准】

进一步的，以每层切片的磁共振图像中的参考图像I_Mr为基准，完成本层切片在T个心跳周期内的磁共振图像的配准。其中M为大于或等于1，小于或等于N的整数。

对于每层切片的磁共振图像的配准包括以本层中的参考图像为基准，采用刚体配准和非刚体配准，当刚体配准没有取得较好效果时，进一步采用非刚体配准。

所述刚体配准包括：从第M层(当前层)切片的参考图像I_Mr向两端各心跳周期的磁共振图像配准；即以第t个心跳周期所对应的参考图像I_Mr分别向心跳周期序号数递减和心跳周期序号数递增所对应的同一层切片的磁共振图像配准。

具体的，将第M层切片在第t个心跳周期的磁共振图像I_Mt进行形变，使得 形变图像I_Mt’与参考图像I_Mr及第t-1个(前一个)心跳周期配准后的磁共振图像I_M(t-1)‘的相似度之和达到最大，此时t＞r；或者使得形变图像I_Mt’与参考图像I_Mr及第t+1个(后一个)心跳周期配准后的磁共振图像I_M(t+1)’的相似度之和达到最大，此时t＜r。

其中，相似度准则为sobel梯度的幅值和角度信息，采用的优化方法为下降单纯型法。

其中，刚体配准分两步进行形变，针对第t个心跳周期的磁共振图像I_Mt，第一步为大尺度(具体为13像素)形变(平移+缩放)，第二步为小尺度(具体为3像素)形变(平移+缩放+旋转)。

进一步的，包括判断第t个心跳周期的磁共振图像I_Mt经过刚体配准后得到的磁共振图像I_Mt’是否需要非刚体配准。

进一步的，经过刚体配准后获得的第t个心跳周期的磁共振图像I_Mt’与参考图像I_Nr的相似度变量(s_ref)，与其第t-1个(前一个)或第t+1个(后一个)心跳周期形变后的磁共振图像I_N(t-1)’的相似度变量(s_pre)，若两变量(s_ref、s_pre)的一阶导的平均值＞0.03；则该经过刚体配准后得到的磁共振图像I_Mt’还需要进行非刚体配准。

从第M层(当前层)切片的参考图像I_Mr向两端同一层切片在各心跳周期的磁共振图像配准。

利用demons非刚体配准方法对当前图像进行形变，与其相应的伪真实图像进行配准。相似度准则为sobel梯度的幅值和角度信息。

若当前图像为参考图像I_Mr前、后一个心跳周期的图像I_M(r-1)、I_M(r+1)，则当前图像对应的伪真实图像P_M(r-1)、P_M(r+1)为参考图像I_Mr，否则，当前磁共振图像I_Mt的伪真实图像P_Mt为第t-1个(前一个)心跳周期的磁共振图像的形变图 像I_M(t-1)’(*0.3)与第t-1个(前一个)心跳周期的伪真实图像P_M(t-1)(*0.7)的加权平均，此时t＞r。或者，当前磁共振图像I_Mt的伪真实图像为第t+1个(后一个)心跳周期的磁共振图像的形变图像I_M(t+1)’(*0.3)与第t+1个(后一个)心跳周期的伪真实图像P_M(t+1)(*0.7)的加权平均，此时t＜r。

【分割心肌外膜】

从起始层开始分割所有层参考相位图像的心肌外膜。

(1)计算当前层切片(即对应的每一层切片，从起始层切片开始)在T个心跳周期内(时间方向)的磁共振图像I_MT的最大强度投影图I_MIP(M)(参图2d)；

(2)以当前层图像已分割的血池区域质心为中心，将当前层参考图像I_Mr和最大投影图I_MIP(M)分别转换为极坐标图像P_Mr和P_MIP(M)(参图2e)；

(3)在极坐标图像P_Mr中，寻找右心室区域和禁区。

a)对极坐标图像P_Mr进行模糊C均值聚类，聚类数目分别与当前层(第M层)在进行心肌内膜分割中寻找血池区域时采用的聚类数目相同；

b)如果当前层(M＝1时)为起始层，在最亮的类中，选取除血池区域外最大的区域为右心室区域；若当前层不是起始层，则在最亮的类中，选取与前一层右心室区域重叠最大的区域为当前层的右心室区域。

c)将血池区域膨胀三个像素后的区域和最亮的类中除去血池和右心室外区域及其以下像素均设为禁区。

(4)在最大投影图的极坐标图像P_MIP(M)中，在禁区的限制下，采用动态规划寻找心肌外膜(参图2f)最优路径；

在动态规划中，射线扫描方法顺序获得的射线(二维转换图像的每一列)被认为是阶段，射线上的点(二维转换图像每一列上的点)被认为是阶段上的候选点，从第一阶段到最后阶段(二维转换图像的第一列到最后一列)具有最小累 积局部能量的路径被认为是最优路径，即心肌外膜。内部能量和外部能量共同决定了局部能量，其中内部能量决定了最优路径的平滑性，而外部能量决定了最优路径位于梯度大的位置。在具体实施过程中：内部能量为相邻候选边界点的跳动距离。外部能量为极坐标图像P_MIP(M)在列方向上的梯度值，右心室所在列的梯度计算方向与其它列的梯度的计算方向相反，禁区内的像素点外部能量设为无穷大。

(5)将最优路径逆转换到直角坐标系中并进行平滑处理。

【灌注分析】

对每层切片在T个心跳周期内的磁共振图像I_MT：

(1)计算左心室腔区域(参图4血池区域中心附件小范围的圆形区域)和心肌区域的平均灰度-时间曲线(参图5)，并平滑；在左心室腔区域的平均灰度-时间曲线上计算起始相位造影剂(药物)开始流入左心室的心跳周期(第一个上升拐点)，在心肌区域的平均灰度-时间曲线上计算造影剂(药物)最后流出左心室的心跳周期(心肌灰度峰值所在心跳周期的序数加10个心跳周期)；

(2)以血池区域的质心点和右心室边界的中心点连线为基线，按照牛眼图划分方法，将心肌分块；心尖层分4块，中间层和心底层分6块；并在参考图像中画出辐条显示分块结果(参图4)；

(3)基于分块结果可以进行分块分析。对心肌的每一块计算平均灰度-时间曲线，平滑并基线校正，根据曲线计算峰值、达峰时间和最大上升梯度；绘制参数牛眼图。

(4)局部像素点分析。对心肌的每一个像素点计算灰度-时间曲线，平滑并基线校正，根据曲线计算峰值、达峰时间和最大上升梯度，在参考图像中用颜色图映射这些参数。

本发明专利提出方法步骤中涉及到的参数根据实际医学图像特点可以进行任意设置；

本发明专利中刚体配准和非刚体配准采用的相似性测度可以进行变化，如应用互信息；

本发明专利中动态规划技术相关的局部能量定义可以根据图像的灰度、梯度和形状等；

本发明专利中心肌部分信号强度随时间变化曲线量化分析相关参数不限于已提供的，如可将心肌平分成内外两圈，心肌内圈和心肌外圈。可对心肌内圈和外圈分别分块分析。

虽然本发明已以较佳实施例揭示如上，然其并非用以限定本发明，任何本领域技术人员，在不脱离本发明的精神和范围内，当可作些许的修改和完善，因此本发明的保护范围当以权利要求书所界定的为准。