CN1327809C

CN1327809C - 在t1加权磁共振成像中用propeller采样方式消除运动伪影的方法

Info

Publication number: CN1327809C
Application number: CNB2005100337726A
Authority: CN
Inventors: 冯衍秋; 陈武凡
Original assignee: Southern Medical University
Current assignee: Southern Medical University
Priority date: 2005-03-28
Filing date: 2005-03-28
Publication date: 2007-07-25
Anticipated expiration: 2025-03-28
Also published as: CN1686050A

Abstract

本发明公开了一种在T1加权磁共振成像中用PROPELLER采样方式消除运动伪影的方法，包括步骤如下：(1)基于PROPELLER采集方式对T1加权成像数据进行采集；(2)在图像域进行相位矫正；(3)由中心重叠采样区域的K空间数据，根据相关测度，进行基于频域相关的运动参数估计；(4)以步骤3中结果出发，基于图像域最大化互信息量算法，进行运动参数的迭代优化搜索；(5)由步骤4中结果，进行运动补偿；(6)通过非笛卡尔数据网格化重建算法，重建出不受运动伪影干扰的图像。本发明可以提高运动估计与补偿的精度与稳健性，相应地更好地消除T1加权图像的运动伪影干扰。

Description

在T1加权磁共振成像中用PROPELLER采样方式消除运动伪影的方法

技术领域

本发明涉及一种图象处理方法，特别涉及在T1加权磁共振成像中用PROPELLER采样方式消除运动伪影的方法，用于医学磁共振成像处理。

背景技术

磁共振成像(MRI)是当前医学影像学的重要检查手段之一，由于其数据采集的时间比较长，病人的运动常常导致图像中出现伪影，图像质量严重恶化。据文献报道，在头部磁共振检查中，为了避免运动的影响，大约14％的受检者需要服用镇静剂或通过麻醉手段，才能得到满意的图像。因此，有效地矫正运动伪影在临床上存在巨大的需求，也一直是医学磁共振成像的研究热点同时也是技术难题之一。

James G.Pipe于1999年提出PROPELLER采样方式(图1)，对于运动伪影有很好的消除效果。其基本思想是：将K空间采样分成几个部分顺序进行，每个部分叫做一个K空间条，在采集每个K空间条的时候，由于采集很快完成，可以近似认为受检者是静止的；而运动仅仅发生在采集K空间条之间，而运动的方向与幅度可以由K空间条之间重叠采样区域的数据来估计得到；对运动加以有效补偿后，就可以重建出不受运动伪影干扰的优质图像。目前PROPELLER方式基本上都是在基于快速成像序列做T2加权磁共振成像。但对T1加权成像，PROPELLER方式的采集及其重建算法中，基于K空间域相关的运动估计的精度与稳健性较低，造成运动伪影消除的效果不理想。

发明内容

本发明的目的在于提出一种在T1加权磁共振成像中用PROPELLER采样方式消除运动伪影的方法，可以提高运动估计与补偿的精度与稳健性，相应地更好地消除T1加权图像的运动伪影干扰。

本发明目的可通过以下技术措施来实现：包括步骤如下：

1、基于PROPELLER采集方式对T1加权成像数据进行采集；

2、通过在图像域进行相位矫正，以保证数据的最大值出现在每个K空间条的中心；

3、由中心重叠采样区域的K空间数据，根据相关测度，进行基于频域相关的运动参数估计；

4、以步骤3中结果出发，基于图像域最大化互信息量算法，进行运动参数的迭代优化搜索；

5、由步骤4中结果，进行运动补偿；

6、通过非笛卡尔数据网格化重建算法，重建出不受运动伪影干扰的图像。

本发明步骤4中基于图像域最大化互信息量的运动参数的迭代优化搜索具体步骤为：

a)将PROPELLER采集到的每个K空间带的数据，补零后进行逆傅立叶变换，取模后得到一组模糊的临时图像Si；

b)对每幅图像Si进行模糊增强；

c)以互信息量为测度，POWELL多维优化算法，以步骤3得到的初始参数出发，迭代搜索使互信息量最大的变换参数x_trans，y_trans，Theta，即运动参数。

本发明中模糊增强的具体过程为：首先通过正弦函数变换将图像灰度变换到模糊隶属度空间(0，1)，模糊隶属度f＝sin(S*PI/(2*Max))，Max为所取图像中最大灰度值，然后选择一个阈值参数β，如果模糊隶属度f取值大于或等于β，增强后模糊隶属度f取值参数α*f²，如果模糊隶属度f取值小于β，增强后f’取值(1-α)*(1-f)²，其中α＝1/(β²+(1-β)²)；最后将增强后模糊隶属度f’变换回图像灰度域即可。

本发明所提基于频域相关与图像域最大化互信息量相结合的运动估计算法，可以更充分地利用现有数据信息，运动估计的精度获得很大提高；采用频域相关为图像域最大互信息量算法提供初始参数估计，先粗略地收敛到全局最优值附近，既可以避免局部极值的影响，又有利于减小搜索步长，从而加速算法运行；而且，模糊增强可以压缩图像中过亮与过暗区域，有助于进一步提高配准的精度与稳健性。

本发明方法与现有技术的实验比较及结果如下：

首先应用本发明方法对静止对象的PROPELLER数据进行了重建实验。图4左边给出了Pipe教授所述算法的重建结果，中间是给定实际运动参数时的重建结果图像，最右边图像是本发明算法的最终结果。图4中第一行图像为头部T1加权数据的重建实验结果，原有算法重建出的图像很模糊(图4a)，特别是在大脑纵裂处等处，而本发明的重建图像(图4c)更为清晰，非常接近给定真实参数时的结果(图4b)；图4中第二行图像是体模T1加权数据的重建实验结果，原有算法的重建图像(图4d)受到了严重的星条状伪影污染，而本发明算法的重建图像(图4f)大大改善了这一情况，更加接近给定真实参数时的重建图像(图4e)。图5a、图5b为两种算法旋转参数误差比较曲线，可以明显看出本发明所提算法，其旋转参数估计的精度与稳健性都大大优于原有的频域相关法。

为了检验本发明算法对运动伪影的消除效果，让志愿者躺在检查床上，并让他在检查过程中进行2～3次摇头运动，进行了T1加权的PROPELLER MRI数据采集。重建图像结果比较见图6，左边为常规采集重建后的图像，很明显受到了严重的运动伪影污染；中间为现有Pipe所提频域算法的重建图像，运动伪影受到了一定的抑制，但伪影消除的并不彻底；右边为本发明所提算法的重建结果，可以看到所成图像更清晰，对运动伪影的消除效果更彻底。图7给出了两个算法估计出来的旋转运动参数，可以看到，在大部分地方，两种算法参数估计结果比较接近，这就为本发明利用频域相关做图像域互信息量搜索的初始参数估计提供了依据。

附图说明

图1 PROPELLER MRI的K空间采样轨迹图；(加粗的线表示采样得到的一个K空间条，然后围绕K空间的中心每旋转一个角度再采集下一个K空间条，直至完成整个K空间的采样)

图2相位矫正流程图；

图3基于图像域最大化互信息量的运动估计算法流程图；

图4受检查物体静止时最终重建图像对照；((a)与(d)Pipe所提算法重建结果(b)与(e)给定真实参数重建结果(c)与(f)本发明算法重建结果)

图5a为头部T1加权数据的旋转参数估计误差曲线图；

图5b为体模T1加权数据的旋转参数估计误差曲线图；

图6志愿者运动时T1加权PROPELLER采集数据的重建结果；((a)T1加权常规采集数据重建后结果，含严重的运动伪影；(b)现有Pipe教授所提算法重建结果；(c)本发明算法的重建结果)

图7志愿者运动时旋转参数估计曲线图。

具体实施方式

本发明实施有六个步骤，具体如下：

步骤1数据采集。按照图1设计采样方式，K空间带数与带内相位编码线数可以灵活设定，本发明的实验中设定为：每10度采集一个K空间条，共采集18个，每个K空间条采集24行相位编码线，每条相位编码线采样256个数据，信号叠加平均次数为1。将采样数据存入一3维数组。

步骤2相位矫正。由于梯度系统的不理想与涡流的存在，使得每个K空间带的最大值并不是出现在其中心位置，表现在图像上就是一个缓慢变化的相位变化，因此可以通过在图像域进行相位矫正，以保证数据的最大值出现在每个K空间条的中心。通过高通滤波器滤除图像中缓慢变化的相位差，以将每个K空间带的最大值调整到中心位置，算法的基本流程见图2。以一个K空间条的数据R(kx，ky)为例，具体如下：

a)对R(kx，ky)补零后进行二维FFT变换，得到图像，设为M1(x，y)；

b)构造金字塔函数：函数在一个256×24的矩阵上取值，函数在中心位置(129，13)处取值1，边缘取值为0，从中心到边缘呈线性递减，记为H(kx，ky)；

c)将R(kx，ky)与H(kx，ky)逐点相乘后进行二维FFT变换，得到图像，记为M2(x，y)；

d)在M1(x，y)数据中减去M2(x，y)相位，得修正后的图像M(x，y)；

e)对M(x，y)进行二维FFT逆变换，得到相位矫正后的K空间条数据。

步骤3运动参数的初始估计。应用频域相关算法，具体为：

a)只考虑每个K空间带内中心附近圆形重叠采样区域内的数据，选择卷积函数为二维Keiser-Bessel函数，通过卷积插值方法将其转化到相同的笛卡尔坐标系下，存入大小为24×24的矩阵，设为Ri，i为K空间带序数；

b)估计K空间带之间的旋转运动参数：以估计R₁与R₂之间旋转参数为例，将取两个矩阵中数据取模后，将R₂旋转某一角度θ后插值生产新矩阵R₂’，计算R₂’和R₁之间的相关系数，记为CORR，很明显CORR为θ的函数，而且CORR的最大值对应的θ值就是R₁与R₂之间的旋转运动参数；

c)估计K空间带之间的平移运动参数：以估计R₁与R₂之间平移参数为例，求矩阵R₂的复数共轭R₂ ^*，将R₂ ^*与R₁相乘后，进行傅立叶变换，得到R₁与R₂之间相对于平移x_trans与y_trans的相关能量CORR，然后分别找到CORR最大值对应的x_trans与y_trans值，就是要求的平移运动参数；

步骤4运动参数估计的进一步细化。由第3步所得结果出发，应用图像域的最大化互信息量算法(见图3)，具体为：

b)对每幅图像Si进行模糊增强，以图像S为例，首先通过正弦函数变换将图像灰度变换到模糊隶属度空间(0，1)，模糊隶属度f＝sin(S*PI/(2*Max))，Max为S中最大灰度值，然后选择一个阈值参数β，如果模糊隶属度f取值大于或等于β，增强后模糊隶属度f’取值参数α*f²，如果模糊隶属度f取值小于β，增强后f’取值(1-α)*(1-f)²，其中α＝1/(β²+(1-β)²)；最后将增强后模糊隶属度f’变换回图像灰度域即可；

c)估计K空间带之间的运动参数：以估计S1与S2之间旋转参数为例，对图像S2进行参数的空间变换S2’，然后计算S2与之间的S1互信息量测度，记为MI，则MI为变换参数(x_trans，y_trans，Theta)的函数；由POWELL多维优化算法，以步骤3得到的初始参数出发，迭代寻找使MI最大的变换参数(x_trans，y_trans，Theta)，也即运动参数。

步骤5运动补偿。由步骤4得到的运动估计结果，根据傅立叶变换的平移与旋转特性，进行运动补偿。具体为：

a)旋转运动补偿，由步骤4估计得到的旋转参数θ，只需将K空间数据相相反方向移动θ角度即可：

b)平移运动补偿，由步骤4估计得到的平移参数(x_trans，y_trans)，

只需改变K空间数据的相位大小即可实现补偿：

步骤6网格化重建。组合所有的K空间带，应用卷积插值的办法，同样选择卷积函数为二维Keiser-Bessel函数，将非笛卡尔坐标系下的数据插值到笛卡尔坐标系下表示；然后对网格化后的K空间数据，进行二维傅立叶逆变换，取模后得到最终重建图像。

Claims

1、一种在T1加权磁共振成像中用PROPELLER采样方式消除运动伪影的方法，其特征在于包括步骤如下：

(1)基于PROPELLER采集方式对T1加权成像数据进行采集；

(2)通过在图像域进行相位矫正，以保证数据的最大值出现在每个K空间条的中心；

(3)由中心重叠采样区域的K空间数据，根据相关测度，进行基于频域相关的运动参数估计；

(4)以步骤3中结果出发，基于图像域最大化互信息量算法，进行运动参数的迭代优化搜索；

(5)由步骤4中结果，进行运动补偿；

(6)通过非笛卡尔数据网格化重建算法，重建出不受运动伪影干扰的图像。

2、根据权利要求1所述的在T1加权磁共振成像中用PROPELLER采样方式消除运动伪影的方法，其特征在于步骤4中基于图像域最大化互信息量的运动参数的迭代优化搜索具体步骤为：

b)对每幅图像Si进行模糊增强；

c)以互信息量为测度，POWELL多维优化算法，以所述步骤(3)中得到的运动参数估计结果出发，迭代搜索使互信息量最大的变换参数x_trans，y_trans，Theta，即运动参数。

3、根据权利要求2所述的在T1加权磁共振成像中用PROPELLER采样方式消除运动伪影的方法，其特征在于模糊增强的具体过程为：首先通过正弦函数变换将图像灰度变换到模糊隶属度空间(0，1)，模糊隶属度f＝sin(S*PI/(2*Max))，Max为所取图像中最大灰度值，然后选择一个阈值参数β，如果模糊隶属度f取值大于或等于β，增强后模糊隶属度f’取值参数α*f²，如果模糊隶属度f取值小于β，增强后f’取值(1-α)*(1-f)²，其中α＝1/(β²+(1-β)²)；最后将增强后模糊隶属度f’变换回图像灰度域即可。