一种心脏图像序列的回顾性脱机呼吸门控方法
【技术领域】
[0001] 本发明涉及一种对在自由呼吸条件下采集的心脏图像序列进行回顾性脱机呼吸 门控，得到去除呼吸运动伪影的图像序列的方法，属于医学成像技术领域。
【背景技术】
[0002] 目前，成像技术在心脏疾病的临床诊断和治疗中起着至关重要的作用。在对心脏 及冠状动脉血管进行成像时，心脏有节律地收缩和舒张导致图像中产生心脏运动伪影，影 响疾病的准确诊断和介入治疗的效果。可以通过心电（ECG)门控技术抑制心脏运动伪影， 即在每个心动周期的同一时相（一般是R波到来时）采集图像。
[0003] 呼吸运动会导致心脏沿上下、左右和前后方向的移动，主运动为上下方向的平移 运动。在呼气末和吸气末之间心脏容易发生旋转和非刚性运动（如仿射形变或非线性形 变）。呼吸运动导致心脏图像中产生呼吸运动伪影，使得图像的分辨率和信噪比降低，在图 像指导的介入治疗过程中则会出现静态指导信息和运动结构的位置不一致的现象。
[0004] 目前临床常用的抑制心脏图像中呼吸运动伪影的方法是控制呼吸法（即屏气法） 和在线门控法。控制呼吸法是在患者憋气的过程中获取图像，使得图像采集不受呼吸运动 的影响。由于该方法限制了图像获取或介入治疗的时间，因此在很多情况下都不适用。比 如，当病人的呼吸控制能力较差，进而严重影响图像的时空分辨率时；或者病人仅能够保持 较短的屏气时间，而采集图像所需时间较长时，该方法均不适用。在线门控法包括呼吸门控 和导航门控，呼吸门控是在采用外置呼吸探测垫或腹带检测呼吸运动的基础上，控制图像 采集设备在特定时间内采集图像，其缺点是操作比较复杂；导航门控不需要任何外接部件， 但需要预先对横膈膜进行正常呼吸下的矢状位和冠状位的扫描，寻找横膈膜置顶位置作为 呼吸末期的导航窗位置，该技术操作复杂，扫描周期较长，适用范围较小。
【发明内容】
[0005] 本发明的目的在于针对现有技术之弊端，提供一种心脏图像序列的回顾性脱机呼 吸门控方法，所述方法操作简单、可在自由呼吸（即不需要患者憋气）条件下进行图像的采 集而又避免产生运动伪影。
[0006] 本发明所述问题是以下述技术方案解决的：
[0007] -种心脏图像序列的回顾性脱机呼吸门控方法，所述方法首先采用流形学习 (ManifoldLearning)方法中的拉普拉斯特征映射（LaplacianEigenmaps，LE)对存储ECG 门控心脏图像序列数据的矩阵进行降维处理，得到嵌入在高维观测数据点集中的低维坐标 矩阵；然后计算低维坐标矩阵中相邻特征向量之间的欧氏距离，检测欧氏距离的局部极大 值，并将其作为门控帧的选取位置，从而得到去除了呼吸运动伪影的门控图像序列。
[0008] 上述心脏图像序列的回顾性脱机呼吸门控方法，具体处理步骤如下：
[0009] a.利用拉普拉斯特征映射的方法对描述心脏图像序列的矩阵进行降维处理：
[0010] ①用一个二维矩阵表示图像序列中所有像素的灰度值，具体步骤如下：
[0011] 假设图像序列包括N帧图像，每帧图像由D=MXM个像素组成，定义一个NXD维 的矩阵X，按照由上至下、由左至右的顺序将每帧图像的像素灰度值排成一行存储到矩阵X 的对应行中；
[0012] ②通过计算矩阵X中各行之间的欧氏距离，建立邻域关系图，并据此来重构数据 流形的局部特征，具体步骤如下：
[0013] 对于矩阵X中的每Rx1GRd，计算\与、之间的欧氏距离d(x^xj):
[0014]d(Xi，Xj) =IlXi-XjIl2，
[0015] 其中，」_ = 1，2，...，~;||*||2表示2范数，将计算出的欧氏距离按照由小到大的 顺序排列，取出前k个所对应的向量（不包括X1本身）作为最近邻点，记为X1,，其 中j= 1，2, . . .，k，以X1为顶点，欧氏距离d(xpXl])为边，建立邻域关系图；
[0016] ③构建NXN维的权值矩阵W= {Wlj|i，j= 1，2，...，N}:
[0017] 如果^与\在邻域关系图中互为邻近点，则将它们之间的权值赋值为W1,= 1 ;否 贝1J，将权值赋值为Wlj= 0 ;
[0018] ④计算拉普拉斯矩阵的特征向量和特征值，完成特征映射，具体步骤如下：
[0019] 高维观测数据点集为X= {XlGRD，i=l，2，...，N}，低维坐标矩阵为Y= {yiGRd， i= 1，2, . . .，N}，其中#是d维实数空间，且d〈〈D，低维嵌入的目标函数为：
[0021] 计算使上式最小的矩阵Y，具体方法如下：
[0022] 设L为拉普拉斯矩阵：
[0023] L=B-W,
[0024] 其中，B=IbljIi,j= 1，2, ...，N}是一个NXN维的稀疏对角阵，其对角线上的 元素为：
[0026] 利用下式求解矩阵L的特征值A和特征向量f:
[0027] Lf=入f
[0028] 计算矩阵L的前d个非零的最小特征值，其所对应的特征向量就是低维坐标矩阵 Y；
[0029] b.通过检测呼吸运动信号的局部极值提取门控帧：
[0030] ①计算矩阵Y中相邻特征向量之间的欧氏距离A(m):
[0031] A(m) = [Y(l，m)-Y(l，m-l)]2+[Y(2,m)-Y(2,m-1)]2+. ? ? [Y(d，m)-Y(d，m-1)]2,其 *m= 2,3，...，N，A(l)=0;
[0032] ②检测A(m)的局部极大值，将其作为门控帧的选取位置，完成门控帧的选取。
[0033] 上述心脏图像序列的回顾性脱机呼吸门控方法，根据患者的呼吸周期长度Cb(单 位：s)和心动周期长度(；（单位：s)，计算ECG门控图像序列中一个呼吸周期内包含的图像 帧数n:
[0035] 如果A(m)的两个局部极大值所对应的帧间隔小于n/2,则选取两个局部极大值中 的较大者作为门控帧的提取位置。
[0036] 上述心脏图像序列的回顾性脱机呼吸门控方法，所述低维坐标矩阵Y={yiGRd， i= 1，2，...，N}中，实数空间Rd的维数d= 2。
[0037] 本发明不需要对各帧图像进行分割和特征提取等预处理，而是直接对由图像中所 有像素的灰度值组成的矩阵进行分析，提取高维图像数据空间的低维流形，得到心脏图像 序列中隐含的呼吸运动信息。本方法只需求解稀疏矩阵的特征值，而且无需操作者的手动 参与，具有计算复杂度低、自动化程度高、简单易行、应用成本低等优点。此外，本方法只利 用了局部距离信息，因而门控结果对噪声不敏感。
【附图说明】
[0038] 下面结合附图对本发明作进一步说明。
[0039] 图1是本发明方法的流程图；
[0040] 图2是邻域关系图G(V，E)的不意图；
[0041] 图3是LE二维嵌入后计算得到的A(m)的时域波形。
[0042] 文中所用符号为：N、原始ECG门控心脏图像序列的总帧数；M、每帧图像中像素的 总行数和总列数；D、每帧图像中的总像素数；X={XlGRD，i= 1，2，. . .，N}、由原始图像序 列中所有像素点的灰度值组成的NXD维的图像数据矩阵；Xl、Xj、矩阵X中的第i行和j行 数据；RD、D维实数空间；d(Xl，X])、\与X,之间的欧氏距离；d、低维流形的维数；Rd、d维实 数空间；k、样本点的最近邻点选取的个数；Xl]、X1的第j个最近邻点；G(V，E)、邻域关系图； V、邻域关系图中的结点集，其元素称为结点；V&V、集合V的无序积；E、邻域关系图中无序积 V&V的多重子集，其元素称为无向边（简称边）；ek、在邻域关系图的图形表示中代表无向边 (Vl， Vj) ;W、权值矩阵；Wlj、权值矩阵W中第i行、第j列的元素；L、拉普拉斯矩阵；A、矩阵 L的特征值；f、矩阵L的特征向量；Y=IyiGRd，i= 1，2, . . .，N}、d维嵌入结果的坐标矩 阵；yi、矩阵Y中的第i行和j行数据；B、通过权值矩阵W计算得到的NXN维的稀疏对 角阵；by矩阵B对角线上的元素；Cb、呼吸周期的长度（单位：s);(；、心动周期的长度（单 位：s) ;n、一个呼吸周期中包含的图像帧数;A(m)、相邻特征向量之间的欧氏距离。
【具体实施方式】
[0043] 下面结合附图1详细说明本发明方法的步骤：
[0044] 步骤1、利用拉普拉斯特征映射的方法对描述心脏图像序列的矩阵进行降维处 理：
[0045] 拉普拉斯特征映射算法（BelkinM，NiyogiP.Laplacianeigenmapsand spectraltechniquesforembeddingandclustering.NeuralInformationProcessing Systems. 2002,14:585-591.)是一种局部流形学习算法，通过保持高维数据点之间的近邻 关系来发现低维流形结构。具体方法如下：