CN107843862B

CN107843862B - 一种propeller技术中参考位置图像的非迭代生成方法

Info

Publication number: CN107843862B
Application number: CN201610835427.2A
Authority: CN
Inventors: 朱高杰; 罗海; 查乐平
Original assignee: ALLTECH MEDICAL SYSTEMS LLC
Current assignee: ALLTECH MEDICAL SYSTEMS LLC
Priority date: 2016-09-20
Filing date: 2016-09-20
Publication date: 2021-04-23
Anticipated expiration: 2036-09-20
Also published as: CN107843862A

Abstract

本发明公开一种PROPELLER技术中参考位置图像的非迭代生成方法，包括以下步骤：采集数据、处理数据、转换数据、构建矩阵、奇异值分解、合成数据和图像重建，本发明不依赖于迭代过程，因此具有较高的计算效率；同时，这种方法根据N个输入的图像生成一个新的参考位置图像，通过计算得到的自适应加权因子算法能够将误差和运动的影响最小化，因此能够在临床多变的扫描条件下得到健壮、稳定的参考位置图像。

Description

一种PROPELLER技术中参考位置图像的非迭代生成方法

技术领域

本发明涉及一种PROPELLER技术中参考位置图像的计算方法，尤其涉及一种PROPELLER技术中参考位置图像的非迭代生成方法。

背景技术

磁共振成像技术是利用氢质子的核磁共振现象进行成像的一种技术。人体内包含单数质子的原子核，例如广泛存在的氢原子核，其质子具有自旋运动。带电原子核的自旋运动，在物理上类似于单独的小磁体，而且在没有外部条件影响下这些小磁体的方向性分布是随机的。当人体置于外部磁场中时，这些小磁体将按照外部磁场的磁力线重新排列具体为在平行于或反平行于外在磁场磁力线的两个方向排列，将上述平行于外在磁场磁力线的方向称为正纵向轴，将上述反平行于外在磁场磁力线的方向称为负纵向轴，原子核只具有纵向磁化分量，该纵向磁化分量既具有方向又具有幅度。

用特定频率的射频(RF，Radio Frequency)脉冲激发处于外在磁场中的原子核，使这些原子核的自旋轴偏离正纵向轴或负纵向轴，产生共振，这就是磁共振现象。上述被激发原子核的自旋轴偏离正纵向轴或负纵向轴之后，原子核具有了横向磁化分量。

停止发射射频脉冲后，被激发的原子核发射回波信号，将吸收的能量逐步以电磁波的形式释放出来，其相位和能级都恢复到激发前的状态，将原子核发射的回波信号经过空间编码等进一步处理即可重建图像。

在磁共振扫描成像过程中，由于检测时间比较长，病人往往会有自主的或无意识的运动。这种运动会造成图像变模糊，更严重的是产生伪影，影响医生对病人的诊断。消除运动伪影是磁共振成像领域研究的重点和技术难题之一。

1999年J.G.Pipe提出了PROPELLER(Periodically Rotated OverlappingParallel Lines with Enhanced Reconstruction)技术。该技术基于K空间非等间隔的数据采集，减少数据采集的时间；利用K空间中心数据带重叠区域的过采样数据估计运动形式并矫正运动伪影；运用网格化算法将非笛卡尔数据转化为笛卡尔数据，再进行傅立叶逆变换最终生成重建图像。这种方法对于刚性运动伪影的消除效果非常显著，已经在头颅磁共振成像中获得了成功应用。

具体的来说，PROPELLER包含以下几个主要步骤：数据采集、相位矫正、旋转矫正、平移矫正、相关性加权和图像重建。其中，最为关键的是对运动的矫正，包括旋转矫正和平移矫正。对运动矫正之前，要估计出每个K空间采集带相对于某个参考位置图像的运动量，包括旋转量和平移量。参考位置图像确定的精确性、健壮性、是否包含伪影等都会对接下来运动量的估计产生重要影响。

为了保证参考位置图像选取的质量，1999年J.G.Pipe采用了取平均的方法生成参考位置图像。有N个同中心的K空间数据带，它们在中心部分是重叠的。将这N个重叠的数据分别网格化分布在笛卡尔坐标系上，然后取平均值。这样，新生成的K空间数据代表将要用于运动量估计的参考位置图像。这种方法比较简单直接，但是缺点也很明显，如果N个数据带之间发生了运动或者某个数据带内部包含有运动伪影，那么采用简单平均方法找到的参考图像几乎必然含有运动伪影。

Ashish A.Tamhane等人在2009年提出采用单个数据带作为参考位置图像的方法[Motion correction in Periodically-Rotated Overlapping Parallel Lines WithEnhanced Reconstruction(PROPELLER)and Truboprop MRI.Mag Reson Med 62:174-182,2009]。在这种方法中，首先计算出N个数据带中心区域的平均值；其次，以该平均的数据带为基准，分别计算同N个数据带的相关系数，选择相关系数最大的数据带作为参考位置图像的数据。这种方法试图通过求出单个可靠的数据带，避免平均方法带来的伪影。但是，一方面，N个数据带中有可能每一个都包含了运动伪影，那么通过上述方法无法保证高质量的参考位置图像；另一方面，上述方法计算的相关系数都是基于平均参考图像的，如果平均参考图像本身含有伪影，那么即使相关系数最大，仍然不能保证被选择的参考位置图像最好。

在2001年J.G.Pipe还提出了一种基于迭代的方法来选择参考位置图像[ImprovedIn-plane motion correction for PROPELLER MRI.ISMRM 2001]。这种方法的思路是，在第一次迭代过程中，选择某个图像作为参考位置图像，然后其他图像都相对于该图像做运动矫正；在第二次迭代中，将上述经过矫正的图像取平均，生成的图像作为参考位置图像。迭代方法的好处是，通过第一次迭代中的运动矫正，提高了参与平均的第二次迭代的N个图像的质量，一定程度上减少了引入的误差和运动伪影。但是，仍然无法保证最终平均后的参考位置图像的质量，主要原因是一方面，仍然使用了平均的方法，那么如果N个输入图像中含有伪影，则平均后这种伪影仍然无法消除；其次，第一次迭代开始的时候，需要选择某个图像作为参考基准，实际上难以保证这个随意选出来的参考图像是否包含运动伪影，如果包含了伪影，那么接下来的运动矫正就会发生错误。

为了提高参考位置图像在双极形式的运动中的精确性，Zhe Liu等人在2014年提出采用一种基于群组分类的方法[Improved motioin correction in PROPELLER bygrouped Blades as Reference]。这种方法的基本思想是，首先对于所有的N个图像，分别计算其两两的相关系数，总共为C_n ²个相关系数。其次，基于上述相关系数的大小，通过分类算法将上述所有图像分成2～3组并且保证2*N/3的图像都被包含在上述组内。最后，包含最多图像的数据组取平均，得到的图像被选作参考位置图像。这种方法能够避免N个图像中某些图像含有伪影对最终参考位置图像的影响。但是，一方面，这种方法需要较大的计算量，分组的过程是迭代的；第二，在分组的过程中很多边界条件的设置具有随意性，难以保证临床中复杂情况下最终参考位置图像的高质量。例如，分组算法中需要根据相关系数将某些近似的组融合成为一个组，但是如果阈值设置太高，容易导致分组数目大于3；如果阈值设置较低，则容易导致分组数目太小而变为1，失去了该算法本身的含义。而阈值的设置又依赖于输入图像本身的质量，不是一个可精确计算的数值，因此，较难以控制其稳定性和精确性。

在2014年J.G.Pipe进一步改进了2001年提出的迭代方法[Revised MotionEstimation Algorithm for PROPELLER MRI.Mag Reso Med.72:430-437,2014]。这种方法的迭代过程是，首先，选择N个图像的平均值作为参考位置图像；其次，分别计算N个图像和上述平均参考位置图像的相关系数，选择相关系数最大作为参考位置图像。并且，将其他N-1个图像相对于该参考位置图像做运动矫正。之后，将上述相关系数作为权重因子，合成上述N个图像生成新的参考位置图像；最后，反复迭代上述步骤，直到中止条件满足，选出参考位置图像。这种迭代的方法本身需要较大的计算量，而且实际的步骤中还需要不断的对图像进行运动矫正，导致计算量更是成倍的提高了计算负担。

发明内容

本发明旨在提供一种PROPELLER技术中参考位置图像的计算方法，该方法能够不依赖于迭代的过程，具有较小的计算量；同时，能够规避N个输入图像中包含的误差和伪影，在临床复杂情况下得到健壮、稳定的参考位置图像。

为达到上述目的，本发明采用的技术方案如下：

一种PROPELLER技术中参考位置图像的非迭代生成方法，包括以下步骤：

步骤01：采集数据，采集N个数据带A，所述N个数据带A为旋转的同中心数据带，所述N个数据带A在K空间的中心区域重叠，N为正整数；

步骤02：处理数据，将步骤01采集到的N个数据带A先进行相位矫正，再进行网格化处理，处理后得到N个数据带B；

步骤03：转换数据，将N个数据带B在K空间中心重叠区域的二维数据转换为一维数据，所述一维数据为一列；每个数据带B含有D个数据，N个数据带B转换为N个D*1的列向量，转换后得到N个数据带M，D为正整数；

步骤04：构建矩阵，构建矫正矩阵C，其大小为D*N，即D行N列；

将步骤03中所述的N个D*1列向量形式的数据带M分别放置于矫正矩阵C的对应列，将第1个数据带M放置在矫正矩阵C的第1列，第i个数据带M放置在矫正矩阵C的第i列，i为正整数；

步骤05：奇异值分解，对矫正矩阵C进行奇异值分解，分解过程采用以下方程：

C＝UΣV^H (Ⅰ)

在式(Ⅰ)中，C为D*N大小的矫正矩阵，U为D*D大小的酉矩阵，V为N*N大小的酉矩阵，V^H为V的共轭转置，Σ为D*N大小的对角阵，对角阵中的每个元素均为矩阵C的奇异值；

步骤06：合成数据，将所述对角阵Σ中的奇异值元素按照降序排列，选择对应于最大奇异值的右奇异向量R中的对应列R_max，所述R_max为N*1大小的列向量；R_max为权重因子，使用R_max对C中的N列数据进行合成，采用以下方程：

D_cop＝C*R_max (Ⅱ)

在式(Ⅱ)中，C为D*N大小的矫正矩阵，R_max为N*1大小列向量形式的权重因子，D_cop为D*1大小形式的合成数据；

步骤07：图像重建，将步骤06合成的一维数据D_cop转换为二维数据，将所述一维数据D_cop转换后的二维数据执行傅立叶变换，得到的图像为基于奇异值分解生成的自适应参考位置图像。

优选的，步骤02中所述相位矫正采用金字塔窗函数滤波的方法。

优选的，步骤02中所述网格化处理为将非笛卡尔坐标系下的数据转换到笛卡尔坐标系下的数据，所述网格化处理采用Jackson网格化处理算法。

优选的，步骤03中所述二维数据转换为一维数据采用列优先原则或者行优先原则，所述N个数据带B均采用相同的转换原则。

本发明提出的用于PROPELLER技术中参考位置图像的计算方法，由于不依赖于迭代过程，对于采集到的N个K空间数据，构造矫正矩阵并对其进行奇异值分解，选择最大奇异值对应的奇异向量作为权重因子，合成N个输入的K空间数据，得到伪影及误差最小化的参考位置图像。因此具有较高的计算效率；同时，这种方法根据N个输入的图像生成一个新的参考位置图像，通过计算得到的自适应加权因子算法能够将误差和运动的影响最小化，因此能够在临床多变的扫描条件下得到健壮、稳定的参考位置图像。

附图说明

图1为本发明的方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。

如图1所示，本发明为解决上述技术问题而采用的技术方案是提供一种PROPELLER技术中的参考位置图像的非迭代生成方法，描述了该方法实现的主要过程，包括如下步骤：

步骤01：按照PROPELLER技术的采集方式，需要采集N个旋转的同中心数据带，这N个数据带在K空间的中心区域是重叠的。如果物体没有发生运动，在排除网格化处理带来的计算误差的情况下，这些重叠的数据应该是等价的。相反，如果在不同的数据带采集期间，被扫描图像发生了旋转或平移运动，那么这些重叠的数据就不等价了。不等价的关键在于两个对比的数据带之间有运动误差。

步骤02：对上述N个数据带进行相位矫正和网格化处理。PROPELLER数据采集中，旋转方向不断改变，实际扫描中沿着K空间编码方向上的梯度系统不均匀以及常见的涡流影响，会导致每个K空间数据带的中心与旋转中心不一致，也就是每个K空间数据带的中心并不总是出现在采集数据的中心上。这种K空间的偏移会导致图像空间有一个缓慢的相位变化。金字塔窗函数滤波的方法能较好的实现相位的矫正。

由于PROPELLER数据采集的过程中，采集方向不断旋转，是一种非笛卡尔采集。为了接下来基于傅立叶逆变换的重建方式，必须将非笛卡尔形式的数据转换到笛卡尔坐标系下。将上述非笛卡尔坐标系下数据转换到笛卡尔坐标系数据的过程就是网格化。PROPELLER方法中一般可采用Jackson网格化处理算法。

步骤03：针对N个数据带,其在K空间中心重叠区域的二维数据转换为一维数据。经过上述步骤的网格化处理，K空间中的数据都已经被转换到了笛卡尔坐标系中。本步骤是要将N个二维数据转化为一维数据。一维数据可以是一行或者一列，本发明以一列数据为例进行描述。

将二维数据转换为一列数据可以采用多种方式，例如可以采用列优先原则或者行优先原则。但无论采用何种方式，都要确保N个K空间中心数据带都采用相同的转换原则。假设每个K空间中心数据带含有D个数据，则转换后为D*1的列向量。

步骤04：构建矫正矩阵C，其大小为D*N，即D行N列。将上述N个列向量形式的K空间中心数据带分别放置于矫正矩阵的对应列。例如，将第1个数据带放置在矫正矩阵C的第一列，第i个数据带放置在矫正矩阵C的第i列。此时的矫正矩阵C包含了N个K空间中心数据带的所有信息。

步骤05：对矫正矩阵C进行奇异值分解。分解过程采用以下方程：

C＝UΣV^H (Ⅰ)

(Ⅰ)中，C为D*N大小的矫正矩阵，U为D*D大小的酉矩阵，V为N*N大小的酉矩阵，V^H为V的共轭转置，Σ为D*N大小的对角阵，对角阵中的每个元素均为矩阵C的奇异值；

步骤06：上述对角阵Σ中的元素是奇异值，一般可按照降序排列。选择对应于最大奇异值的右奇异向量R中的对应列R_max，其大小为N*1列向量。这个列向量代表了一种权重因子,利用该权重因子对C中的N列数据进行合成，采用以下方程：

D_cop＝C*R_max (Ⅱ)

(Ⅱ)中，C为大小为D*N的矫正矩阵，R_max为大小为N*1的权重因子，D_cop为大小为D*1的合成后的数据。物理上，该数据代表了对误差及运动伪影最小化优化后的参考位置图像。由于该方法不需要迭代就能够自主的排除伪影等非主要信号成分的干扰，而且也不需要人为的设置阈值进行剔除、分类等处理，因此称作自适应算法。

上述步骤中描述了基于奇异值分解的自适应参考位置图像的计算方法。将输入的多个相似数据，譬如矫正矩阵C中的每一列，它们包含有较多近似成分，但也包含了各自独特的运动特征，在数学上进行分析，分析的依据是计算后的奇异值及其奇异向量。

对应于较大奇异值的奇异向量，代表了输入矩阵包含的主要的成分，而较小的奇异值对应的奇异向量则共同代表了较为次要的成分。在参考位置图像的计算中，图像本身代表了主要成分，而误差以及运动伪影代表了次要成分。通过选择最大的奇异值对应的奇异向量，就能够排除不利效应的影响，得到稳定、健壮、非迭代的参考位置图像。

步骤07：将上述合成的一维数据转换为二维数据，再将二维数据执行傅立叶变换后，得到的图像即为采用基于奇异值分解生成的自适应参考位置图像。

基于奇异值分解的自适应参考位置图像，既可以用于计算用于旋转矫正的参考位置图像，也能够用于计算用于平移矫正的参考位置图像。

当然，本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims

1.一种PROPELLER技术中参考位置图像的非迭代生成方法，其特征在于包括以下步骤：

C＝UΣV^H (Ⅰ)

D_cop＝C*R_max (Ⅱ)

步骤07：图像重建，将步骤06合成的一维数据D_cop转换为二维数据，将所述一维数据D_cop转换后的二维数据执行傅立叶变换，得到的图像为基于奇异值分解生成的自适应参考位置图像；

步骤02中所述的相位矫正采用金字塔窗函数滤波的方法；步骤02中所述的网格化处理为将非笛卡尔坐标系下的数据转换到笛卡尔坐标系下的数据，所述网格化处理采用Jackson网格化处理算法；

步骤03中所述二维数据转换为一维数据采用列优先原则或者行优先原则，所述N个数据带B均采用相同的转换原则。