CN105069785A - 一种心脏图像序列的回顾性脱机呼吸门控方法 - Google Patents

Abstract

一种心脏图像序列的回顾性脱机呼吸门控方法，所述方法首先采用流形学习方法中的拉普拉斯特征映射对存储ECG门控心脏图像序列数据的矩阵进行降维处理，得到嵌入在高维观测数据点集中的低维坐标矩阵；然后计算低维坐标矩阵中相邻特征向量之间的欧氏距离，检测欧氏距离的局部极大值，并将其作为门控帧的选取位置，从而得到去除了呼吸运动伪影的门控图像序列。本发明直接对由图像中所有像素的灰度值组成的矩阵进行分析，得到心脏图像序列中的呼吸运动信息。本方法只需求解稀疏矩阵的特征值，而且无需操作者的手动参与，具有计算复杂度低、自动化程度高、应用成本低等优点。此外，本方法只利用了局部距离信息，因而门控结果对噪声不敏感。

Description

一种心脏图像序列的回顾性脱机呼吸门控方法

技术领域

本发明涉及一种对在自由呼吸条件下采集的心脏图像序列进行回顾性脱机呼吸门控，得到去除呼吸运动伪影的图像序列的方法，属于医学成像技术领域。

背景技术

目前，成像技术在心脏疾病的临床诊断和治疗中起着至关重要的作用。在对心脏及冠状动脉血管进行成像时，心脏有节律地收缩和舒张导致图像中产生心脏运动伪影，影响疾病的准确诊断和介入治疗的效果。可以通过心电(ECG)门控技术抑制心脏运动伪影，即在每个心动周期的同一时相(一般是R波到来时)采集图像。

呼吸运动会导致心脏沿上下、左右和前后方向的移动，主运动为上下方向的平移运动。在呼气末和吸气末之间心脏容易发生旋转和非刚性运动(如仿射形变或非线性形变)。呼吸运动导致心脏图像中产生呼吸运动伪影，使得图像的分辨率和信噪比降低，在图像指导的介入治疗过程中则会出现静态指导信息和运动结构的位置不一致的现象。

目前临床常用的抑制心脏图像中呼吸运动伪影的方法是控制呼吸法(即屏气法)和在线门控法。控制呼吸法是在患者憋气的过程中获取图像，使得图像采集不受呼吸运动的影响。由于该方法限制了图像获取或介入治疗的时间，因此在很多情况下都不适用。比如，当病人的呼吸控制能力较差，进而严重影响图像的时空分辨率时；或者病人仅能够保持较短的屏气时间，而采集图像所需时间较长时，该方法均不适用。在线门控法包括呼吸门控和导航门控，呼吸门控是在采用外置呼吸探测垫或腹带检测呼吸运动的基础上，控制图像采集设备在特定时间内采集图像，其缺点是操作比较复杂；导航门控不需要任何外接部件，但需要预先对横膈膜进行正常呼吸下的矢状位和冠状位的扫描，寻找横膈膜置顶位置作为呼吸末期的导航窗位置，该技术操作复杂，扫描周期较长，适用范围较小。

发明内容

本发明的目的在于针对现有技术之弊端，提供一种心脏图像序列的回顾性脱机呼吸门控方法，所述方法操作简单、可在自由呼吸(即不需要患者憋气)条件下进行图像的采集而又避免产生运动伪影。

本发明所述问题是以下述技术方案解决的：

一种心脏图像序列的回顾性脱机呼吸门控方法，所述方法首先采用流形学习(ManifoldLearning)方法中的拉普拉斯特征映射(LaplacianEigenmaps，LE)对存储ECG门控心脏图像序列数据的矩阵进行降维处理，得到嵌入在高维观测数据点集中的低维坐标矩阵；然后计算低维坐标矩阵中相邻特征向量之间的欧氏距离，检测欧氏距离的局部极大值，并将其作为门控帧的选取位置，从而得到去除了呼吸运动伪影的门控图像序列。

上述心脏图像序列的回顾性脱机呼吸门控方法，具体处理步骤如下：

a.利用拉普拉斯特征映射的方法对描述心脏图像序列的矩阵进行降维处理：

①用一个二维矩阵表示图像序列中所有像素的灰度值，具体步骤如下：

假设图像序列包括N帧图像，每帧图像由D＝M×M个像素组成，定义一个N×D维的矩阵X，按照由上至下、由左至右的顺序将每帧图像的像素灰度值排成一行存储到矩阵X的对应行中；

②通过计算矩阵X中各行之间的欧氏距离，建立邻域关系图，并据此来重构数据流形的局部特征，具体步骤如下：

对于矩阵X中的每行x_i∈R^D，计算x_i与x_j之间的欧氏距离d(x_i，x_j)：

d(x_i，x_j)＝‖x_i-x_j‖²，

其中，j＝1，2，...，N；‖·‖²表示2范数，将计算出的欧氏距离按照由小到大的顺序排列，取出前k个所对应的向量(不包括x_i本身)作为x_i的k-最近邻点，记为x_ij，其中j＝1，2，...，k，以x_i为顶点，欧氏距离d(x_i，x_ij)为边，建立邻域关系图；

③构建N×N维的权值矩阵W＝{w_ij|i，j＝1，2，...，N}：

如果x_i与x_j在邻域关系图中互为邻近点，则将它们之间的权值赋值为w_ij＝1；否则，将权值赋值为w_ij＝0；

④计算拉普拉斯矩阵的特征向量和特征值，完成特征映射，具体步骤如下：

高维观测数据点集为X＝{x_i∈R^D，i＝1，2，...，N}，低维坐标矩阵为Y＝{y_i∈R^d，i＝1，2，...，N}，其中R^d是d维实数空间，且d<<D，低维嵌入的目标函数为：

$min\underset{i,j}{\&Sigma;}\left|\right|y_i-y_j|{|}^2\&CenterDot;w_{ij},$

计算使上式最小的矩阵Y，具体方法如下：

设L为拉普拉斯矩阵：

L＝B－W，

其中，B＝{b_ij|i，j＝1，2，...，N}是一个N×N维的稀疏对角阵，其对角线上的元素为：

$b_{ii}=\underset{j=1}{\overset{N}{\&Sigma;}}{w}_{ij},$

利用下式求解矩阵L的特征值λ和特征向量f：

Lf＝λf

计算矩阵L的前d个非零的最小特征值，其所对应的特征向量就是低维坐标矩阵Y；

b.通过检测呼吸运动信号的局部极值提取门控帧：

①计算矩阵Y中相邻特征向量之间的欧氏距离A(m)：

A(m)＝[Y(1，m)-Y(1，m-1)]²+[Y(2，m)-Y(2，m-1)]²+...[Y(d，m)-Y(d，m-1)]²，其中m＝2，3，...，N，A(1)＝0；

②检测A(m)的局部极大值，将其作为门控帧的选取位置，完成门控帧的选取。

上述心脏图像序列的回顾性脱机呼吸门控方法，根据患者的呼吸周期长度C_b(单位：s)和心动周期长度C_r(单位：s)，计算ECG门控图像序列中一个呼吸周期内包含的图像帧数n：

$n=\frac{C_b}{C_r},$

如果A(m)的两个局部极大值所对应的帧间隔小于n/2，则选取两个局部极大值中的较大者作为门控帧的提取位置。

上述心脏图像序列的回顾性脱机呼吸门控方法，所述低维坐标矩阵Y＝{y_i∈R^d，i＝1，2，...，N}中，实数空间R^d的维数d＝2。

本发明不需要对各帧图像进行分割和特征提取等预处理，而是直接对由图像中所有像素的灰度值组成的矩阵进行分析，提取高维图像数据空间的低维流形，得到心脏图像序列中隐含的呼吸运动信息。本方法只需求解稀疏矩阵的特征值，而且无需操作者的手动参与，具有计算复杂度低、自动化程度高、简单易行、应用成本低等优点。此外，本方法只利用了局部距离信息，因而门控结果对噪声不敏感。

附图说明

下面结合附图对本发明作进一步说明。

图1是本发明方法的流程图；

图2是邻域关系图G(V，E)的示意图；

图3是LE二维嵌入后计算得到的A(m)的时域波形。

文中所用符号为：N、原始ECG门控心脏图像序列的总帧数；M、每帧图像中像素的总行数和总列数；D、每帧图像中的总像素数；X＝{x_i∈R^D，i＝1，2，...，N}、由原始图像序列中所有像素点的灰度值组成的N×D维的图像数据矩阵；x_i、x_j、矩阵X中的第i行和j行数据；R^D、D维实数空间；d(x_i，x_j)、x_i与x_j之间的欧氏距离；d、低维流形的维数；R^d、d维实数空间；k、样本点的最近邻点选取的个数；x_ij、x_i的第j个最近邻点；G(V，E)、邻域关系图；V、邻域关系图中的结点集，其元素称为结点；V&V、集合V的无序积；E、邻域关系图中无序积V&V的多重子集，其元素称为无向边(简称边)；e_k、在邻域关系图的图形表示中代表无向边(v_i，v_j)；W、权值矩阵；w_ij、权值矩阵W中第i行、第j列的元素；L、拉普拉斯矩阵；λ、矩阵L的特征值；f、矩阵L的特征向量；Y＝{y_i∈R^d，i＝1，2，...，N}、d维嵌入结果的坐标矩阵；y_i、y_j、矩阵Y中的第i行和j行数据；B、通过权值矩阵W计算得到的N×N维的稀疏对角阵；b_ii、矩阵B对角线上的元素；C_b、呼吸周期的长度(单位：s)；C_r、心动周期的长度(单位：s)；n、一个呼吸周期中包含的图像帧数；A(m)、相邻特征向量之间的欧氏距离。

具体实施方式

下面结合附图1详细说明本发明方法的步骤：

步骤1、利用拉普拉斯特征映射的方法对描述心脏图像序列的矩阵进行降维处理：

拉普拉斯特征映射算法(BelkinM，NiyogiP.Laplacianeigenmapsandspectraltechniquesforembeddingandclustering.NeuralInformationProcessingSystems.2002，14:585-591.)是一种局部流形学习算法，通过保持高维数据点之间的近邻关系来发现低维流形结构。具体方法如下：

首先，用一个二维矩阵表示图像序列中所有像素的灰度值，具体步骤如下：

假设图像序列包括N帧图像，每帧图像由D＝M×M个像素组成，定义一个N×D维的矩阵X，按照由上至下、由左至右的顺序将每帧图像的像素灰度值排成一行存储到矩阵X的对应行中。

其次，通过计算矩阵X中各行之间的欧氏距离，建立邻域关系图，并据此来重构数据流形的局部特征，具体步骤如下：

对于矩阵X中的每行x_i∈R^D，其中i＝1，...，N，R^D是D维实数空间，x_i是一个1行、D列的向量，计算x_i与x_j之间的欧氏距离d(x_i，x_j)：

d(x_i，x_j)＝‖x_i-x_j‖²(1)其中，j＝1，2，...，N；‖·‖²表示2范数。将计算出的欧氏距离按照由小到大的顺序排列，取出前k个所对应的向量(不包括x_i本身)作为x_i的k-最近邻点，记为x_ij，其中j＝1，2，...，k。以x_i为顶点，欧氏距离d(x_i，x_ij)为边，建立邻域关系图，它是一种无向图。下面举例说明无向图G(V，E)，V是结点集，E是无序积V&V＝{(x，y)|x∈V，y∈V}的多重子集，其元素称为无向边(简称边)。如附图2所示，设V＝{v₁，v₂，...，v₅}，E＝{(v₁，v₁)，(v₁，v₂)，(v₂，v₃)，(v₂，v₅)，(v₁，v₅)，(v₄，v₅)}，当无向图的集合定义转化成图形表示之后，常用e_k表示无向边(v_i，v_j)，附图2中的e₁表示(v₁，v₁)，e₂表示(v₁，v₂)，e₃表示(v₁，v₅)，e₄表示(v₂，v₅)，e₅表示(v₂，v₃)，e₆表示(v₄，v₅)。

再次，构建N×N维的权值矩阵W＝{w_ij|i，j＝1，2，...，N}，如果x_i与x_j在邻域关系图中互为邻近点，则将它们之间的权值赋值为w_ij＝1；否则，将权值赋值为w_ij＝0。

最后，计算拉普拉斯矩阵的特征向量和特征值完成特征映射，具体步骤如下：

高维观测数据点集为X＝{x_i∈R^D，i＝1，2，...，N}，低维坐标矩阵为Y＝{y_i∈R^d，i＝1，2，...，N}，其中R^d是d维实数空间，且d《D。低维嵌入的目标函数为：

$min\underset{i,j}{\&Sigma;}\left|\right|y_i-y_j|{|}^2\&CenterDot;w_{ij}---\left(2\right)$

计算使式(2)最小的矩阵Y，具体方法如下：

设L为拉普拉斯矩阵：

L＝B-W(3)

其中，B＝{b_ij|i，j＝1，2，...，N}是一个N×N维的稀疏对角阵，是利用权值矩阵W计算得到的，其对角线上的元素为：

$b_{ii}=\underset{j=1}{\overset{N}{\&Sigma;}}{w}_{ij}---\left(4\right)$

矩阵B中除对角线外的其余元素均为0。求解矩阵L的特征值和特征向量：

Lf＝λf(5)

其中，λ和f分别是矩阵L的特征值和特征向量。计算矩阵L的前d个非零的最小特征值，其所对应的特征向量就是对高维数据的低维嵌入结果Y。

步骤2通过检测呼吸运动信号的局部极值提取门控帧：

首先，根据患者的呼吸周期长度C_b(单位：s)和心动周期长度C_r(单位：s)，计算ECG门控图像序列中一个呼吸周期内包含的图像帧数n。由于ECG门控图像采集方式是在每个心动周期中采集一帧图像，所以：

$n=\frac{C_b}{C_r}---\left(6\right)$

然后，计算矩阵Y中相邻特征向量之间的欧氏距离A(m)。对于高维空间d维嵌入的计算公式为：

A(m)＝[Y(1，m)-Y(1，m-1)]²+[Y(2，m)-Y(2，m-1)]²+...[Y(d，m)-Y(d，m-1)]²(7)其中，m＝2，3，...，N，A(1)＝0。实验证明，当d＝2时，A(m)能够较好地反映出呼吸运动的周期变化规律。因此本发明方法取d＝2，则A(m)的计算公式如下：

A(m)＝[Y(1，m)-Y(1，m-1)]²+[Y(2，m)-Y(2，m-1)]²(8)

其中m＝2，3，...，N，A(1)＝0。

最后，检测A(m)的局部极大值，它表示前一个呼吸周期的末尾和后一个呼吸周期的开始，将其作为门控帧的选取位置。如果两个局部极大值所对应的帧间隔小于n/2(如附图3中的局部峰值1和局部峰值2)，则选取两个局部极大值中的较大者作为门控帧的提取位置，完成门控帧的选取。

Claims (4)

1.一种心脏图像序列的回顾性脱机呼吸门控方法，其特征是，所述方法首先采用流形学习方法中的拉普拉斯特征映射对存储ECG门控心脏图像序列数据的矩阵进行降维处理，得到嵌入在高维观测数据点集中的低维坐标矩阵；然后计算低维坐标矩阵中相邻特征向量之间的欧氏距离，检测欧氏距离的局部极大值，并将其作为门控帧的选取位置，从而得到去除了呼吸运动伪影的门控图像序列。

2.根据权利要求1所述的一种心脏图像序列的回顾性脱机呼吸门控方法，其特征是，具体处理步骤如下：

a.利用拉普拉斯特征映射的方法对描述心脏图像序列的矩阵进行降维处理：

①用一个二维矩阵表示图像序列中所有像素的灰度值，具体步骤如下：

假设图像序列包括N帧图像，每帧图像由D＝M×M个像素组成，定义一个N×D维的矩阵X，按照由上至下、由左至右的顺序将每帧图像的像素灰度值排成一行存储到矩阵X的对应行中；

②通过计算矩阵X中各行之间的欧氏距离，建立邻域关系图，并据此来重构数据流形的局部特征，具体步骤如下：

对于矩阵X中的每行x_i∈R^D，计算x_i与x_j之间的欧氏距离d(x_i,x_j)：

d(x_i,x_j)＝||x_i-x_j||²；

其中，j＝1,2,...,N；||·||²表示2范数，将计算出的欧氏距离按照由小到大的顺序排列，取出前k个所对应的向量(不包括x_i本身)作为x_i的k-最近邻点，记为x_ij，其中j＝1,2,…,k，以x_i为顶点，欧氏距离d(x_i,x_ij)为边，建立邻域关系图；

③构建N×N维的权值矩阵W＝{w_ij|i,j＝1,2,...,N}：

如果x_i与x_j在邻域关系图中互为邻近点，则将它们之间的权值赋值为w_ij＝1；否则，将权值赋值为w_ij＝0；

④计算拉普拉斯矩阵的特征向量和特征值，完成特征映射，具体步骤如下：

高维观测数据点集为X＝{x_i∈R^D,i＝1,2,...,N}，低维坐标矩阵为Y＝{y_i∈R^d,i＝1,2,...,N}，其中R^d是d维实数空间，且d＜＜D，低维嵌入的目标函数为：

$min\underset{i,j}{\&Sigma;}\left|\right|y_i-y_j|{|}^2{w}_{ij},$

计算使上式最小的矩阵Y，具体方法如下：

设L为拉普拉斯矩阵：

L＝B-W，

其中，B＝{b_ij|i,j＝1,2,...,N}是一个N×N维的稀疏对角阵，其对角线上的元素为：

$b_{ii}=\underset{j=1}{\overset{N}{\&Sigma;}}{w}_{ij},$

利用下式求解矩阵L的特征值λ和特征向量f：

Lf＝λf，

计算矩阵L的前d个非零的最小特征值，其所对应的特征向量就是低维坐标矩阵Y；

b.通过检测呼吸运动信号的局部极值提取门控帧：

①计算矩阵Y中相邻特征向量之间的欧氏距离A(m)：

A(m)＝[Y(1,m)-Y(1,m-1)]²+[Y(2,m)-Y(2,m-1)]²+…[Y(d,m)-Y(d,m-1)]²，

其中m＝2,3,...,N，A(1)＝0；

②检测A(m)的局部极大值，将其作为门控帧的选取位置，完成门控帧的选取。

3.根据权利要求2所述的一种心脏图像序列的回顾性脱机呼吸门控方法，其特征是，根据患者的呼吸周期长度C_b(单位：s)和心动周期长度C_r(单位：s)，计算ECG门控图像序列中一个呼吸周期内包含的图像帧数n：

$n=\frac{C_b}{C_r},$

如果A(m)的两个局部极大值所对应的帧间隔小于n/2，则选取两个局部极大值中的较大者作为门控帧的提取位置。

4.根据权利要求3所述的一种心脏图像序列的回顾性脱机呼吸门控方法，其特征是，所述低维坐标矩阵Y＝{y_i∈R^d,i＝1,2,...,N}中，实数空间R^d的维数d＝2。