心脏动态MRI图像分割的时空多尺度网络

徐佳陈，肖志勇

江南大学人工智能与计算机学院， 无锡 214122

0 引 言

心脏动态磁共振成像技术在临床诊断中具有重要的意义。根据其分割结果，医生可以有效地获取心肌质量和厚度、射血分数以及心室体积等诊断指标。然而准确地对其进行分割是一项具有挑战性的任务，因为在成像的过程中由于不均匀的磁场强度，容易产生伪影，造成边界模糊。不同受试者的器官差异很大，尤其是右心室的形状大小多变，容易产生体积效应。此外动态磁共振成像技术本身就具有短轴切片数量少，切片厚度过厚等特点，导致图像的短轴分辨率过低、信息量稀疏等问题。现阶段的心脏分割方法依旧以手动分割方式为主。但这种手动分割的方法耗时费力且对医生的专业素质要求严格。因此，如何自动对心脏磁共振图像(cardiac magnetic resonance imaging，CMRI)进行精确的分割已经成为一个重要的医学图像问题。

深度学习在医学图像分割领域大放异彩。与此同时，提出了许多基于深度学习的心脏分割算法。其中Chang等人(2018)提出了利用YOLO(you only look once (Redmon等，2016))网络提取感兴趣区域，再输入到U型神经网络的方法，从而减少无关信息对分割效果的影响。Zotti 等人(2019)将心脏结构的先验信息嵌入网络，并为心脏分割提出了新的损失函数。Liu等人(2020)提出了将残差模块(He等，2016)与双向长短记忆模块相结合的方法以减少不符合物理特征的分割异常点的出现。Hasan和Linte(2020)提出群卷积与权重剪枝技术相融合的U型神经网络，增加网络深度的同时减少了参数量。高强等人(2020)提出了一种基于组归一化与最近邻插值的分割算法，能快速准确地提取特征信息。刘畅等人(2021)为了在分割心肌内膜与外模时保证其空间位置关系，提出了先用胶囊网络提取出相对位置信息和大小等信息，再输入到全连接网络进行预测的方法。然而现存的方法通常存在以下几个问题：1)针对动态心脏MRI图像本身短轴(z轴)切片厚度过厚、信息稀疏的特点，现有方法为避免3D卷积带来的边界效应(Baumgartner等，2018)往往会使用2维网络来进行分割，使用单帧的3维图像在z轴的切片作为网络的输入。这样会丢失图像在3维层面的上下文信息，产生目标器官之间的类内不确定、类间不清晰的问题。2)由于现有数据集仅有动态MRI图像在收缩期末期和舒张期末期两个时间点的标签值(ground truth)，现有网络通常也仅将这两个时间点的图像作为分割对象，这忽略了动态MRI图像在时间序列上的信息。3)由于心脏结构中细小器官的形状、大小多变的特性，现有方法往往会产生一些不符合目标器官相对物理特性的异常点(左心肌包含了右心室的情况)。

为解决以上3个问题，提出了一种新的多输入、多分支和多任务的伪3D分割网络ST UNet(spatio-temporal UNet)。该网络的主要创新点如下：1)为了充分利用心脏MRI影像的时间信息解决目标器官的类内不确定性问题，以连续时间内的多帧输入代替传统的单帧输入，让网络学习更为丰富的时域特征。2)利用目标帧学习可变形偏移量，并将其融合到时域特征，使其具有更为精准的边缘信息。3)为了减小细小器官的分割错误，更好地融合时域特征与边缘信息，提出可变全局注意力连接来代替传统的长连接，为网络提供更加广泛的全局上下文信息。4)针对由2维卷积造成的上下文信息丢失所引起的目标器官边界模糊、类间不清晰的问题，提出根据方向场特征的学习来进行特征向量的矫正。

1 本文方法

1.1 网络结构阐述

鉴于心脏动态MRI图像的切片过厚、层间信息稀疏的特点，提出了基于U型神经网络的伪3D分割网络(ST UNet)的结构如图1(a)所示。该网络的基本骨架来源于Ronneberger等人(2015)提出的UNet结构。该结构由编码层、解码层与跳跃连接组成。其中编码层又称为下采样层，其作用为提取高级语义特征，并利用最大池化层减小特征图尺寸来扩大卷积的感受野。解码层又称为上采样路径，其作用是将高级语义特征图逐步还原成原始图片的大小，实现端到端的分割效果。跳跃连接则是将编码层学习到的高级特征图扩充到解码路径，从而避免特征提取过程中信息丢失的问题，使得解码层还原出来的结果更为精准。UNet解码层的还原结果可表示为

图1 网络结构图Fig.1 Network structure diagram((a)structure diagram of deformable global attention connection；(b)overall structure of the algorithm)

xi=c([αi,up(xi+1)])

(1)

式中，x表示解码层的还原结果，α表示编码层的高级语义特征，i表示当前所在的层数，[·]表示通道连接操作，up(·)表示上采样操作，c(·)表示解码操作。

ST UNet保留了UNet中编码块、解码块与跳跃连接层的基本结构框架。但为充分利用心脏动态MRI图像的时空信息，ST UNet对这3个结构做了不同改进，形成了一种多输入、多路径和多任务的新型分割网络。首先在编码阶段，ST UNet提出了多路径、多输入的时间信息聚合编码块。将目标帧和关键序列分别输入到不同的路径。规范化卷积对关键序列提取丰富的时域特征，可变形卷积(Dai等，2017)提取目标帧中更为精准的边缘信息。其次将特征通过可变形全局注意力连接扩充到解码路径。

最后在网络解码过程中，规定ST UNet在预测原始分割图的同时还需额外预测分割图所对应的空间方向场，再利用预测到的空间方向场对原始的分割图进行特征矫正，得到改进后的分割结果。

1.2 时间增强编码块与可变形全局连接

现有的大多数方法都忽略了心脏动态MRI图像的时间信息，仅以收缩期末期和舒张期末期的图像作为单独的输入，从而导致分割性能下降。因此ST UNet选择将需要分割的图像作为目标帧，以包含目标帧的连续时间序列一同作为网络的输入，该连续序列称为关键序列。

心脏动态MRI图像的维度表示为[x,y,z,t]，其中x-y为短轴平面，z为短轴，t为时间维度。考虑到MRI心脏图像在z轴上存在较大的层间间隙，选择在时间维度t上连续的x-y平面的图像组成关键序列。鉴于数据中心脏MRI舒张期末期的时间维度t均为0，选择从目标帧开始的后r帧组成关键序列。即假设ct0∈RH×W表示在t0时刻的目标帧，则关键序列为{ct0,…,ct0+r}。其中H，W分别代表了特征向量的高和宽。

关键序列的编码块由两个3×3的卷积核组成，并在每个卷积核后面添加批量归一化(batch-normalization)和ReLU激活函数。时间信息增强后的特征可以表示为

(2)

式中,Fkey(·)表示增强后的特征向量，p0表示特征向量中的任意位置，pn∈R表示卷积核中的任意位置，w(pn)表示卷积核在位置pm的权重值，而Ct(p0+pn)表示图像Ct在(p0+pn)位置的值。

此外，由于心脏随着时间变化的跳动，目标器官的边缘信息也会随着时间的变化产生变化，从而造成关键序列中不同时间点的边界信息冲突，因此在对关键序列的卷积过程中可能会产生边缘语义信息的冲突问题。为解决此问题，本文在编码过程中增加了一条针对目标帧Ct0的可变形编码路径，利用可变形卷积学习偏移量δ来增强卷积对于目标器官形状的适应能力，更精准地提取目标器官的信息。

规范卷积中卷积核对于输入的特征图的采样位置都是固定的，如3×3大小的卷积核的采样位置固定为R={(-1,-1),(-1,0),…,(0,1),(1,1)}，其输出的特征图可以直接表示为

(3)

式中，x(p0+pn)表示图像x在(p0+pn)位置的值。

固定的采样方式限制了卷积核对不同位置、不同尺寸或变形物体的特征提取能力,可变形卷积通过对卷积核中每一个采样点增加一个偏移变量δ，改变卷积核采样点的位置。

由图2可以看出，相较于常规的卷积，可变形卷积更加适应几何形变，其原理是让网络在训练过程中额外学习每一个卷积核的位置偏移矩阵δ，使得卷积的采样点发生偏移，实现在当前位置附近的随意采样，提高网络的灵活性。可变形卷积中输出的特征图表示为

图2 可变形卷积与规范卷积Fig.2 Deformable convolution and standard convolution((a)standard convolution；(b)deformable convolution)

(4)

式中，Δpn∈δ表示在pn位置的权重的位置偏移量。

再将由关键序列生成的特征与目标帧生成的特征进行融合，生成当前层编码块的特征图，最后将融合后的特征图作为下一层关键序列编码块的输入，该特征图表示为

F(p0)=Fkey(p0)+Ftar(p0)

(5)

融合后的特征图不仅包含了心脏动态MRI图像中丰富的时间信息，也保留了目标器官在t0时刻(收缩期末期ed或舒张期末期es)的精准边缘信息。

考虑到目标器官之间的相对位置固定，而实验中会产生一些不符合该相对位置的异常点，因此本文受GC block(global context block)(Cao等，2019)的启发设计了可变形全局注意力连接，为解码层提供一个全局的注意力矩阵，让网络在解码还原的过程中获取图像全局上下文信息，从而使得还原的结果更为准确。

可变形全局注意力连接的结构如图1(b)所示。该结构主要分为上下文建模部分、特征转化部分以及可变形卷积层部分。

上下文建模部分利用1×1的卷积和sigmoid激活函数获取特征向量在x-y平面的注意力权值，然后将该权值与图像相乘得到全局的上下文特征。特征转化部分则通过两个1×1的卷积wv1与wv2捕获图像通道间的依赖，并利用层归一化来降低优化难度。该过程可以表示为

(6)

(7)

式中，Np为图像的位置数量，i∈[1,Np]，xi表示第i个位置的输入，zi是对应的输出。αj是采用了嵌入式高斯表达式后的注意力权值，LN表示层归一化操作。

最后利用可变形卷积层捕获心脏的几何信息，将其与zi以及原始特征进行融合，并将融合后的结果传给解码层，进行解码还原。

与原始网络仅仅将编码层的结果传给解码层相比，可变形全局连接后的特征是全局特征的加权和。在解码还原的过程中解码层可以根据权值进行有选择的聚合，提高类内紧凑型和语义一致性。

1.3 空间方向场

受Cheng等人(2020)的启发，为了进一步增强网络预测的结果的类间差异性与类内相似性，提出在预测原始分割图的同时预测对应的空间方向场。

空间方向场是由真实标签图所生成，其生成规则为对每一个前景像素点p，找到目标器官边界最近的一点b。然后对向量bp以b到p的距离为|bp|进行归一化，最终像素点p的方向场可以表示为

(8)

由式(8)可以看出，空间方向场为每一个像素提供了从边界指向中心区域的唯一路径，揭示了目标器官上像素与像素的方向关系。

用一个简单的方向场模块来学习空间方向场。该模块由一个1×1的卷积层组成，模块的输入为UNet最后一个解码层还原出来的64通道的特征，输出为双通道的方向场特征。其中一个通道反映像素点x轴之间的方向关系，另一个通道反映y轴之间的方向关系。

图3 真实标签图与其对应的空间方向场Fig.3 Ground truth and spatial directional field((a)ground truth；(b)spatial directional field)

根据空间方向场的指导，本文利用中心区域的特征对原始分割图进行矫正。如图4所示，最终特征图Ffinal是由原始特征图F0和预测的方向场结合得到的，矫正方式为根据方向场提供的位置信息对原始特征图进行双线性插值。换句话说，让原始特征图F0上的像素点p0根据学习到的方向场进行修改,其整个优化过程为

图4 特征矫正过程图Fig.4 The process of feature correction

∀p∈Ω,F(p)=F(px+DF(p)x,py+DF(p)y)

(9)

式中，Ω表示所有的前景像素点集，px表示像素点p的x轴坐标，py表示像素点p的y轴坐标，F(p)表示经过修正后像素点p的值。由式(9)可以看出，像素点p的矫正后的特征值是根据方向场位置偏移后的像素值。

1.4 损失函数

本文的损失函数涉及分割图的损失函数Lseg以及空间方向场的损失函数LDF。

分割图的损失函数为交叉熵损失函数与Dice损失函数的联合函数。其中交叉熵是将预测结果与真实标签图的每一个像素进行对比，然后对所有的对比结果求平均值。而Dice损失函数是评价预测图与真实标签值的重叠程度，重叠程度越高，分割效果越好。联合损失函数可以表示为

Lseg=0.5×Lcross+Ldice

(10)

(11)

(12)

式中，Lcross表示交叉熵损失函数，Ldice表示Dice损失函数。N表示样品的总数，yi表示第i个样本的预测结果，mi表示第i个样本的真实标签值。

对于空间方向场的学习，以L2范式距离和角度之间的距离作为方向场的损失函数，即

(13)

(14)

式中，|Ci|表示当前标签的所有像素的个数，N表示所有标签的个数。

2 实验结果与讨论

2.1 实验数据与预处理

使用ACDC(Automated Cardiac Diagnosis Challenge)(Bernard等，2018)自动心脏分割挑战数据集，该数据集是法国Dijon大学医院根据当地委员会制定的规定进行采集的。采集设备为两台磁场强度不同的MRI扫描仪：1.5 T磁场与3.0 T磁场。该数据集包含150个志愿者的心脏核磁共振图像，150名志愿者又分为5个亚组，分别为：1)健康的志愿者;2)既往有心肌梗塞的志愿者;3)扩张型心肌病的志愿者;4)肥厚型心肌病的志愿者;5)右心室异常的志愿者。每个志愿者的数据包含心脏整个动态电影MRI图像，以及图像中舒张期末期以及收缩期末期的真实标签值。该数据集的真实标签值是由两名资深医生手动绘制的。动态MRI图像的空间分辨率为1.37～1.68 mm2/像素之间，并且动态图像由28～40幅图像组成，完全或部分覆盖了心跳周期。该数据集将100名志愿者划分为训练集，50名志愿者划分为测试集，训练集与测试集各亚组的分布比例相同并且并未公布测试集的真实标签值。本文所有实验结果都来自于ACDC的官方在线测试平台(https://acdc.creatis.insa-lyon.fr/)。

对ACDC数据集进行如下预处理：1)本文中关键序列参数r取值为1，因此从短轴动态心脏MRI视频中获取收缩期末期tes与舒张期末期ted以及其各自后1帧的时间点ted+1与tes+1时间点的3维图像。2)将3维图像沿z轴进行2维切片，单独将ted或tes时间点所得的切片作为目标帧，再将目标帧对应的ted+1或tes+1时间点所得的切片与目标帧组成关键序列。3)为增加训练集的数据量，提升模型的泛化能力与鲁棒性，本文同时对目标帧与关键序列做相同的数据增强，数据增强的方式有随机仿射变化与随机旋转。4)由于数据集中图像的尺寸各不相同，将所有的切片都剪切或者填补成256×256像素。5)对图像进行归一化操作。

2.2 实验环境与参数设置

实验硬件环境为Intel Core i7处理器与NVIDIA GTX1080ti显卡。编码语言为python，深度学习框架选用的是pytorch 1.4.0版本。权重的初始化方式为pytorch的默认初始化，优化器选用了Adam优化器。学习率α初始化为0.001，权重衰减率为1E-5。网络的通道数从上到下分别为64，128，256，512，1 024。

2.3 评判指标

本文选用ACDC在线测试平台给出的两种评判标准来定量评估ST UNet的性能。

这两个评判标准分别为Dice相似性系数(Dice similarity coefficient, DSC)与豪斯多夫距离(Hausdorff distance, HD)，其中Dice相似性系数DSC表示预测结果与真实标签图的重叠程度，其值在[0,1]之间。而HD系数解释了预测结果与真实标签图的最大不匹配程度。DSC的数学表达式为

(15)

HD的数学表达式为

(16)

HD=max{dpy,dyp}

(17)

式中，P与Y分别表示预测结果与真实标签图，p、y为预测结果与真实标签图上的点。

2.4 消融实验

为探究本文算法的可行性，做了4组消融实验来说明空间方向场、时间增强编码块和可变形全局连接对于实验结果的影响，实验结果如表1和表2所示，LV表示左心室，RV表示右心室，MYO表示左心肌。图5为ST UNet模型在测试集上的DSC箱型图。图6是4种方法的分割效果对比图，分别选取了收缩期末期与舒张期末期4个不同患者进行对比。在分割图中，红色、蓝色和绿色分别表示右心室、左心肌和左心肌。

图5 ST UNet 测试集DSC箱型图Fig.5 Test set DSC box diagram of ST UNet

表1 消融实验DSC的对比结果Table 1 DSC contrast results of ablation experiment

表2 消融实验HD的对比结果Table 2 HD contrast results of ablation experiment

造成类内不确定、类间不清晰的主要原因是心脏动态MRI本身短轴切片数量少，分辨率过低，上下文信息量稀疏。稀疏的信息量造成以往的网络难以适应目标器官复杂的形状和成像时不同的磁场强度。而相邻时间的切片之间的形变较小，尤其是背景与目标器官中心位置的信息并未因心脏跳动产生较大的变化，该信息依旧有效，可用于分割网络。时间增强编码块中的关键序列编码路径正是利用同一切片在相邻时间点的时间信息来弥补稀疏的3维上下文信息，使得网络可以学习到的信息量倍增，增强了网络的适应能力，以此来解决类内不确定的问题。

然而相邻时间切片之间，目标器官的边界信息由于心脏跳动发生了变化，不同时间点的边界信息是相互冲突的。针对这一问题，时间增强编码块中的目标帧编码路径利用可变形卷积提取出目标器官在ed时刻或es时刻精准的边界信息，并通过可变形全局连接将精确的边缘信息与丰富的时间信息融合后扩散到解码路径，使得网络解码还原过程中可以根据权值进行有效的聚合，还原出更为精准的边缘信息，从而解决目标器官与目标器官之间、目标器官与背景之间的类间不清晰的问题。

对比表1和表2中ST UNet的实验结果和没有添加时间增强编码块与可变形全局连接结构的实验结果可以看出，添加后的实验结果中3个目标器官的DICE值和HD距离都有了明显的改善。结合图6的分割结果图，可以发现，ST UNet的分割结果更加接近真实标签，心肌不连贯、心室误分割等情况也都得到了改善。这充分体现了提出的利用时间信息来弥补稀疏的3维上下文信息的特征增强方式是可行的。

同时从表1，表2与图6可以看出原始UNet分割容易出现边界模糊、将无关背景误分割成目标器官以及相邻器官分割错误等问题。通过加入额外的空间方向场训练与矫正后，边界模糊的现象有所改善，3个目标器官的分割精度也有所提升，但是依旧有边缘不平滑以及误分割的现象。在此基础上再添加多输入的时间分支，其分割精度有了进一步的提升，这样更加证明了利用关键序列来进行特征增强能提高网络学习特征的能力，改善网络的分割效果。虽然此时的分割效果相较于前两个模型，边界平滑，没有明显的误分割现象，但是由于关键序列中不同时间点的图像的目标器官的边缘冲突问题，此时的分割结果的边缘出现了欠分割的情况。而加入全局注意力连接后的模型，利用目标帧的全局信息弥补了关键序列边缘冲突的问题，从而使得分割结果达到了更好的效果，更加接近真实标签。

图6 不同模型的分割结果Fig.6 Segmentation results of different models((a)ground truth；(b)ST UNet；(c)UNet + time branch + spatial directional field；(d)UNet + spatial directional field；(e)UNet)

2.5 与其他算法的对比

不同算法对比结果如表3和表4所示。其中Liu等人(2020)是利用残差网络与双向记忆网络来提取相邻切片的3维空间特征。Painchaud等人(2019)是在对抗性编码器(aVAE)中添加了一个潜在的心脏形态。Sun等人(2020)使用二次形流来捕获心脏的形态特征。Simantiris和Tziritas(2020)提出了新扩展卷积网络和自定义的损失函数。

由表3和表4可以看出，本文算法的DSC与HD值都要优于其他方法，这是因为这些算法都仅仅以收缩期末期与舒张期末期的图像作为网络的训练对象，忽视了心脏视频在时间序列上的信息。而HD值的降低也从侧面表示利用空间方向场来修正边界的可行性。因此提出利用心脏动态MRI图像的时间信息来进行特征增强，再利用空间方向场提取出的空间信息来进行特征矫正的方法是可行的。截止2020年12月7日，与ACDC排行榜上的算法对比，ST UNet的3个目标器官的平均DSC与平均HD值分别达到了92.5%与8.90 mm的成绩，仅次于第1名。

表3 不同算法的DSC比对结果Table 3 DSC comparison results of different algorithms

表4 不同算法的HD对比结果Table 4 HD comparison results of different algorithms

3 结 论

针对心脏磁共振图像自动分割的类内不确定、类间模糊等问题，本文提出了利用时间信息来进行特征增强，再利用空间信息来进行特征矫正的ST UNet。ST UNet的主要特点如下：1)相较于现有算法，本文将需要分割的3维图像作为目标帧，将包含目标帧的连续时间序列作为关键序列一起输入网络。2)利用关键序列学习丰富的类内特征，再利用可变形路径学习目标帧上更为精准的边缘信息，通过编码路径与可变形全局注意力连接将两者充分融合。3)网络额外学习空间方向场信息，并利用学习到的空间方向场对原始的预测结果进行矫正。实验结果表明，相较于其他算法，ST UNet能够得到更为精准的分割结果，各器官的分割精度均取得了提升。

尽管本文方法取得了一定进展，但依旧存在如下问题：收缩期末期的分割效果普遍差于舒张期末期的分割效果，原因是在收缩期末期，目标器官排血后体积变小，尤其是心尖与心室基层的形变较大，从而导致目标器官与背景的标签不平衡。后续工作将针对这一问题展开。