基于特征追踪的磁共振衍生心肌应变评估

摘要

测量心肌组织应变等参数的准确且实用的方法具有重要的临床价值，因为已有研究表明，与常用指标射血分数（EF）相比，应变是检测收缩功能障碍更敏感且更早期的指标。

摘要

目的：准确且实用地测量心肌组织应变等参数在临床上具有重要价值，因为已有研究表明，与常用指标射血分数（EF）相比，应变是检测收缩功能障碍更敏感且更早期的指标。目前用于心脏磁共振（CMR）的技术耗时较长，难以在临床实践中推广应用。特征追踪（feature tracking）是一种可实现医学影像定量分析更高自动化和更强稳健性的技术，与现有方法相比，其分析耗时更短。

方法：通过单相位的自动或手动输入进行初始化，系统由此开始追踪代表解剖结构的各个图像模式随时间的位移。该方法的特点在于，图像无需在前期进行任何处理，例如心脏磁共振（CMR）图像的标记。

结果：该方法非常适合用于追踪肌肉组织，从而实现对心肌以及血流的定量分析。

结论：这种新方法提供了一种稳健且节省时间的流程，可在常规心脏磁共振（CMR）成像序列上对心肌组织和血液进行位移、速度和形变参数的定量分析，因此可应用于临床实践。

方案

1. 引言

边界自动检测是图像分析中的一个基本问题。在心脏成像中，若能实现左心室成像中内膜边界的自动检测，便可客观地测量心室容积及心肌形变（应变）。这一目标已通过超声心动图中的斑点追踪技术得以实现。然而，开发可靠的心脏磁共振（CMR）图像自动边界检测方法仍是一项具有挑战性的任务，目前尚未获得普遍可靠的解决方案。事实上，在临床实践中，边界通常由操作者手动勾画，或由软件检测心肌与心腔之间的界面^1,2。本文介绍一种不同的方法，其中边界并非“检测”得出，而是“追踪”获得，即从一个可靠且已存在的瞬时轮廓开始，随时间连续追踪其变化。该初始轮廓通常——但不一定——由经验丰富的操作者在单帧图像上手动绘制。构成该初始可靠轮廓的各个点，通过在后续帧中搜索其邻域内具有相似特征的对应点而实现时间上的追踪。所追踪的特征可以是心腔-组织边界，也可以是沿组织分布的不同解剖结构。这些特征通过在两帧图像间的两个感兴趣区域中采用最大似然法进行识别。

局部帧间位移等效于评估局部速度（位移与时间间隔的比值）。某一点的速度自动评估是通过比较连续两帧中该点周围图像数据的位移来确定的。这类方法在多个不同研究领域中已有多种不同形式的应用。它们属于高级图像分析中称为“光流法”（Optical Flow）的一般范畴^3,4。当利用此类速度追踪生理运动时，在超声成像中通常称之为“斑点追踪”（Speckle Tracking）^{5, 6}；但这些方法同样适用于其他任何图像模态，例如在心脏磁共振（CMR）中，这些方法被称为“特征追踪”或“边界追踪”。

2. 材料与方法

特征追踪方法

在某一任意帧上手动描记二维心脏磁共振电影（CMR cine）的心内膜或心外膜边界（见图1）。也可以描记心肌中层特征。该边界随后被定义为由N个点组成的序列，每个点由其坐标对（x_i,y_i）表示，其中i=1...N。边界追踪通过逐点追踪实现，该追踪基于多尺度的分层算法，并结合一维追踪技术（可保证更高的准确性）与二维追踪技术（对于准确检测空间上延展的二维特征是必要的）。

为了首先捕捉边界的大尺度几何位移，跟踪在垂直于边界的方向上进行，在该方向上心腔-组织边界最为清晰可辨。沿此方向的跟踪采用跨壁切线法进行，具体如下（见图3）。绘制一条穿过心壁、经过该点且与之垂直的直线。沿此跨壁线获取的像素被排列成列，每一列对应图像序列中的一个帧。通过这种方式，跨壁切线的演变过程可在一幅二维图像中同时呈现，其中一个轴表示沿直线的距离，另一个轴表示时间（见图2）。这种表示方法类似于超声心动图中所称的M型图像，在心脏磁共振（CMR）中则更接近于“定位像”（scout）功能。为了提高分析质量，在图像质量较差、信噪比较低的情况下，可采用一条厚度为5像素的线来构建时空图像。随后，边界跟踪即在该时空图像上进行。

在第二步中，为校正边界的二维位移，对每个点独立地在 MxM 移动窗口上执行标准的二维追踪（基于光流法），该窗口始终以先前估计的边界点为中心。二维追踪分两步进行，其中第一次估计结果的一半用于在第二次追踪过程中对移动窗口进行中心定位。随后，在另外两次迭代中，窗口大小从 32 逐步减小至 16。

为了提高沿边界运动的追踪精度，从而更准确地估计旋转和扭转，采用沿时空图像进行一维追踪的方法，这些时空图像是通过沿与弯曲边界“平行”的较厚切片构建的（见图3）。在每个点上，独立地选取沿移动边界、以边界运动点为中心的像素，并将其按列排列，每列对应图像序列中的一个帧。为了提高分析质量，并更好地捕捉边界处的特征，该线向组织内部延伸5个像素（心内膜下）。随后，沿时空图像采用上述相同的方法进行边界追踪。为确保追踪边界的空间一致性，在每一步对相邻点计算出的位移应用三点中值滤波器和三点高斯滤波器（权重为0.25、0.5、0.25）。

沿二维时空图像进行追踪

本部分介绍一种从某一时刻的已知位置出发，在二维图像（类似M模式）中沿一个方向追踪边界的步骤。

x 定义为水平方向，y 定义为垂直方向。各列为 x_i，其中 i = 1...M，M 为图像中的列数。追踪通过确定一组离散的实数序列 y_i = y(x_i) 实现，该序列从已知点 y_k 开始，对应于列 x_k。

通过计算位于 x_k 处的整列与位于 x_k+1 处的整列之间的互相关，来估计从已知点 y_k 到点 y_k+1 的位移。互相关函数会出现一个最大值，该最大值的位置给出了使两列之间相似性最大的所需垂直位移值，因此通过将该位移量加到 y_k 上来估计 y_k+1。此过程在所有相邻列对之间重复进行，最终得到整个边界 y_i（i=1...M）的估计值。此处使用快速傅里叶变换算法计算互相关，以减少计算时间。

首个估计值 y_i 将通过迭代进一步优化。为实现此目标，通过在前一个估计值 y_i 上下各取若干点来提取图像的一个子集，并生成一幅新的图像，其中心对应于序列 y_i，用于校正跟踪。该优化过程将重复进行，直至不再发现校正量为止。

随后，通过最终的蛇行算法[5]，在时空图像中追踪经过固定点 y_k 的图像亮度水平，从而获得更优且更自然的结果。整个过程利用时间周期性，以确保结果具有周期性，并避免漂移效应。

特征追踪的技术局限性

边界追踪技术与任何散斑追踪方法类似，其原理是基于对连续帧之间像素亮度变化的量化。该原理决定了速度的下限：必须能够观察到某一帧中占据一个像素的散斑，在下一帧中移动到相邻像素。因此，该下限表示为

公式 1
其中 Δx 为像素尺寸，Δt 为两帧之间的时间间隔。系数 k 取决于追踪算法的质量及其评估亚像素动态变化的能力。该下限意味着：远高于此限值的速度可被非常精确地估算；而当速度接近或低于此限值时，估算的准确性将下降。

这一限制还意味着，一方面提高采集帧率（减小 Δt）有助于更准确地评估较大的速度及其快速变化（例如在等容阶段期间）。另一方面，帧率的提高（Δt 的减小）会增大该限制，导致对较低速度的评估精度下降，除非同时相应提高空间分辨率（减小 Δx）。

模体图像制备

已制备一系列人工计算机生成的环状结构，以在简单且完全可控的条件下测试图像分析流程。为此，制备了理想左心室短轴投影的模型，具体步骤如下。

心内膜和心外膜边界分别由半径为 R₀(t) 和 R₁(t) 的两个同心圆表示。图像通过在黑色背景上将代表两边界之间组织的环形区域设置为均匀的灰色来制备。随后应用8×8的顶帽线性滤波器，以避免出现非物理性的不连续现象。

心外膜运动以[mm]为单位，表示为 R₀(t)=10+5cos(2πt/T)，其中 T 为心动周期，取为 T=1s。理论心内膜运动学在边界上恒定，仅依赖于时间，速度仅为径向，由 V0(t)= dR₀/dt=-π sin(2πt/T) 给出，单位为[cm/s]。相对于边界在零时刻长度计算的应变百分比为 St₀(t)=100x(R₀(t)-R₀(0))/R₀(0)=100(cos(2πt/T) -1)/3，应变率根据式(1)可得为 SR₀(t)=10 V0/R₀，单位为[s^-1]。假设心外膜运动遵循恒定厚度，即 R₁(t)= R₀(t)+5mm，或保持静止，即 R₁(t)= R₀(0)+5mm。

每幅图像均为边长48 mm的正方形，以组织环为中心，分辨率为N×N。示例图像见图4的图版a和b；应变和应变率的时间曲线见图4的图版c和d。通过改变分辨率N、帧率FR以及心外膜运动类型来制备这些环状图像。
心内膜追踪法通过在第一帧图像上选取一定数量的点来应用于此类图像 N_p 沿圆形心内膜均匀分布的点。

3. 代表性结果

模体研究

此处分析了图像分析方法在计算机生成的模拟图像上的应用。通过均方根百分比差异来计算总体误差的度量。心内膜应变的均方根误差、平均误差和最大误差定义为

公式 2
其中 St₀(t) 为精确值，St(t) 为图像分析计算得到的值，求和范围覆盖所有帧 N_F=FRxT。该定义同样适用于半径、速度和应变率。追踪结果几乎不受心内膜位置的影响，不同点之间的差异远低于 1%。

表I总结了15个模型的结果，这些模型具有不同的空间分辨率、帧率以及心外膜边界的运动类型；同时还展示了用于追踪心内膜边界点数变化的影响。

对于积分量（半径和应变），所有情况下的误差都非常小；而对于与前者导数相关的微分量（速度和应变率），误差则稍大一些。这是预料之中的，因为求导运算会放大误差。当分辨率降低时，结果质量会下降；事实上，精度与像素尺寸相关，像素尺寸在某种意义上代表了从一帧到下一帧可识别的最小位移。只要帧率足够，时间分辨率对结果的影响并不显著；当帧率非常高时，结果不再改善，因为逐帧之间的位移变得小于像素尺寸。这表明，若不同时提高空间分辨率，单纯增加帧率几乎无法带来任何效益。

然而，此处所考虑的简单正弦运动并不需要极高的时间分辨率。同样，仅使用8个点就足以追踪简单的圆形心内膜轮廓。心内膜的结果几乎不受心外膜所经历的运动类型的影响。我们还验证了结果不受所采用图像滤波方法的显著影响。

图4展示了在两个时间点的模拟图像上叠加计算得到的心内膜边界（图a和图b），以直观呈现结果。对于案例#1和分辨率较低的案例#8，应变和应变率分别在图c和图d中给出。在案例#1中（方块符号），应变和应变率与理论值高度一致，平均误差分别为0.6%和3%。在案例#8中，图像分辨率降低一半，一致性略低，应变和应变率的误差分别为0.9%和4.5%。

临床验证 ¹。

我们对191名杜氏肌营养不良症（DMD）患者和42名年龄匹配的健康对照者进行了比较，根据年龄和心脏功能障碍严重程度分组，采用谐波相位成像（HARP）和特征追踪（FT）技术评估左心室中段整层心肌的周向应变（ε_cc）。对匹配的标记图像（tagged）和稳态自由进动（SSFP）图像切片进行了回顾性离线分析。在全部研究人群（n=233）中，FT测得的平均ε_cc（-13.3 ± 3.8%）与HARP测得的平均ε_cc（-13.6 ± 3.4%）高度相关，Pearson相关系数为0.899。DMD患者组中，HARP测得的平均ε_cc（-12.52 ± 2.69%）与FT测得的结果（-12.16 ± 3.12%）差异无统计学意义（p=NS）。同样，对照组中HARP测得的平均ε_cc（-18.85 ± 1.86%）与FT测得的结果（-18.81 ± 1.83%）差异亦无统计学意义（p=NS）。我们得出结论：在心脏功能障碍程度广泛的大量DMD患者人群中，基于FT的ε_cc评估结果与基于标记图像的ε_cc具有高度相关性。

模型编号	帧率	分辨率 N	心外膜运动	Np	εR	εRmax	εV	εVmax	εSR	εSRmax	εSt	εStmax
1	32	401	含心内膜	16	0.57	1.10	3.98	8.26	2.95	6.55	0.59	1.34
2	32	401	无运动	16	0.27	0.60	3.62	8.53	2.90	5.99	0.75	1.29
3	64	401	含心内膜	16	0.35	0.62	1.75	4.92	3.08	10.47	0.34	0.86
4	64	401	无运动	16	0.39	0.73	1.78	5.48	2.38	6.46	0.32	0.63
5	128	401	含心内膜	16	0.22	0.40	1.57	3.66	2.73	8.39	0.35	0.89
6	16	401	含心内膜	16	0.41	0.92	12.10	18.27	7.12	15.17	0.88	1.93
7	32	201	含心内膜	16	0.43	0.82	3.84	7.03	4.21	7.79	0.83	1.52
8	32	201	无运动	16	0.49	0.90	4.00	8.54	4.54	10.84	0.93	1.64
9	32	101	含心内膜	16	2.58	3.87	5.70	17.44	9.22	18.27	4.44	6.77
10	32	401	含心内膜	64	0.32	0.73	3.98	8.09	3.11	8.86	0.48	1.10
11	32	401	无运动	64	0.33	0.53	3.56	7.83	2.78	7.44	0.70	0.98
12	32	401	含心内膜	32	0.44	1.00	3.99	8.16	2.99	6.82	0.63	1.49
13	32	401	含心内膜	8	0.20	0.41	3.64	6.84	3.09	9.10	0.43	0.74
14	32	401	无运动	8	0.22	0.41	3.36	7.42	2.76	5.93	0.58	1.01
15*	32	401	含心内膜	16	1.35	2.42	5.54	9.68	7.02	16.58	2.62	4.46

表1. 心内膜边界追踪的模拟体分析：针对不同模拟体参数，计算主要参数的均方根误差和最大百分比误差 [%]。加粗标记的参数表示与模拟体 #1 的差异。考虑了帧率、分辨率和追踪点数量的影响。针对心外膜运动类型的两种极限情况，分析其影响：心外膜边界不移动（无运动）或与心内膜同步移动（无增厚）。最后一个模拟体 (*) 在构建时未对亮度在一个像素内突变的基本阶跃图像进行滤波处理。误差超过 10% 的数值以加粗标出。

图1. 左心室的cMR图像，左侧为长轴视图，右侧为短轴视图，图像上叠加了勾画出的心内膜边界。

图2. 图像序列沿心肌壁切面的空间-时间表示。心肌壁切面的起始点与图3中的起始点相同。起始点随时间的变化通过自动追踪获得并展示。

图3. 左心室长轴切面图像，显示透壁切面及沿曲边平行的切面。

图 4. 模体研究。两张图像（病例 #2）分别显示最大扩张（图 a）和最大收缩（图 b）时的计算心内膜边界点重叠情况。应变（图 c）和应变率（图 d）通过两个不同模体（病例 #1 和 #8）计算得出，并与有效值进行对比。

图 5. 正常患者中的整体周向应变（黑色曲线）与节段性应变（彩色曲线）示例（a）。左心室功能减低合并左束支传导阻滞患者的整体周向应变（黑色曲线）与节段性周向应变（彩色曲线）示例（b）。注意峰值周向应变时间的不同，提示存在左心室不同步

讨论

我们介绍一种新方法——“特征追踪”（feature tracking），用于测量由SSFP图像获得的心肌组织的形变和位移参数。该方法最近在评估杜氏肌营养不良症（Duchenne Muscular Dystrophy）患者群体的峰值整体圆周应变时，已与HARP技术进行了验证对比（1）。我们发现这两种技术之间具有极好的相关性。特征追踪技术操作相对简便，无需任何额外的成像序列，如标记（tagging）、SENSE或DENSE序列。除了经过临床验证的峰值整体圆周应变（Ecc）评估外，还可获得区域应变数据。特征追踪（FT）似乎是一种可靠的技术，能够估测左心室的纵向和径向应变，以及右心室的纵向应变。

为了准确定义该技术在图像时间分辨率和空间分辨率方面的局限性，还需要开展更多研究。由于基本原理相同，傅里叶变换（FT）可实现与其他模态（如超声心动图）之间的比较。

结论

特征追踪技术有望被整合到临床实践中，因为它无需耗时的分析过程和额外的采集步骤。FT可用于评估左心室和右心室的整体及节段性应变，以及左心室不同步性。目前正在进行更多的临床研究，以评估FT在临床心脏病学中的应用价值。

材料

本文使用的材料清单

姓名	公司	目录编号	评论
Diogenes® 特征追踪软件	TomTec Imaging Systems, Fulda, Germany