专利名称：磁共振图像刚性平移运动伪影的自动逆向修正方法
技术领域：
本发明涉及一种磁共振(MR)图像运动伪影的修正方法，尤其是涉及一种磁共振(MR)图像刚性平移运动伪影的自动逆向修正方法。
背景技术：
由于MR具有无放射线伤害、分辨率很高等优点，故在临床诊断中的应用越来越广泛。但是MR成像时间较长，病人的自主性运动(如呼吸、心跳等)和非自主性运动会在图像中产生运动伪影，严重损害图像质量，从而影响对疾病的诊断。对于运动伪影的消除，基本分为两大类，即前处理与后处理。前处理技术中，主要使用快速扫描、门控扫描、欠扫描、过扫描、半扫描或可校正的重建模式，这些模式基本上也是成功的。但可校正重建去伪影的模式需要对各类运动建立起数学模型，而此类特定数学模型及其他扫描重建模式的事先选择参数是非常困难，因为患者的运动方式是随机的，不可预测的。至于后处理去伪影法，目前对此问题的研究，大都是在刚性运动的范畴内。现有方法需要人工参与，没有实现自动修正。事实上，任何实际的复杂刚性运动都是由两种简单的运动平移和旋转复合而成，因此对平移和旋转运动造成的运动伪影的研究和消除技术为进一步研究实际的复杂运动产生的运动伪影奠定了基础，具有重要的意义。
MR刚性平移运动伪影消除的有关理论和技术可参阅下面三篇描述论文。
1、德里.M，严.H，罗斯费尔德.D“一种二维平移运动伪影修正的改进方法”，IEEE Trans.Med.Ima.，vol.10，No.4，pp.548-553，Dec.19912、瑞查.宙洛斐.A等“一种MRI图像平面平移运动伪影修正的改进方法”，IEEE Trans.Med.Ima.，Vol.14，No.3，PP471-479，Sept.19953、瑞查.宙洛斐.A等“刚性运动引起的MR图像平面伪影的消除方法”，IEEE Tans.Med.Ima.，Vol.15，No.6，PP768-784，Dec.1996国际上对刚性平移运动伪影的消除研究较少，上面三篇论文提出的对刚性运动伪影的消除方法需要在图像域手工标识出图像空间约束区域，不能完全自动修正。国内目前尚无这方面的研究。

发明内容
本发明的目的在于提供一种磁共振图像刚性平移运动伪影的自动逆向修正方法，使得因刚性平移而损害的数据经修正后能重建出不含运动伪影的高质量的图像。
本发明的目的可通过以下的技术措施来实现，依次包括以下步骤(1)获取含有刚性平移运动伪影的磁共振图像的原始频域矢量矩阵数据，即K空间数据；(2)以频率编码方向为Kx方向，相位编码方向为Ky方向，把相位编码梯度Gy为零时定为Ky等于零确定K空间坐标；运用傅立叶变换的周期特性，将所获各原始频域矢量矩阵转换成在K空间坐标中以相位编码方向Ky等于零为对称轴的上下共扼对称的矩阵数据；并获得呈矢量共扼对称的矩阵中每个频域矢量的相位，即K空间坐标中每个K空间线的原相位；(3)以K空间坐标中Ky等于零为参照K空间线，沿坐标轴向正方向的第一个K空间线开始，逐个修正每一个K空间线已偏移的相位；(4)除设定为必须要逆向估计线间运动距离的第一个K空间线之外，对要修正的K空间线，首先确定该K空间线是否需要进行逆向估计线间运动距离dx和dy，所述的线间运动距离dx和dy为两个K空间线线间的运动距离，即所需修正的K空间线相对于前一已修正好的K空间线的运动距离；(5)若需要逆向估计线间运动距离dx和dy的K空间线时，采用逆向修正得出所求的线间运动距离dx和dy，然后根据此线间运动距离dx和dy加上前一已被修正好的K空间线的实际运动距离x_和y_得出该被修正K空间线的实际运动距离x和y，据此实际运动距离x和y修正该K空间线；(6)若不需要逆向估计线间运动距离dx和dy的K空间线时，可直接把前一已修正好的K空间线的实际运动距离x_和y_加上该被修正K空间线前两个已修正好K空间线的实际运动距离之差作为该被修正K空间线的实际运动距离x和y，然后据此实际运动距离x和y修正该K空间线；(7)重复步骤4、5、6中的处理过程，修正全部沿坐标轴向正方向的K空间线，对Ky等于零为参照的向负方向的各个K空间线的修正，以Ky等于零为对称轴取正方向对称行的共扼矢量即可；(8)将所有修正好的K空间线的频域矢量运用傅立叶变换的周期特性，变换成跟初始矩阵一致的非共扼对称矩阵，然后通过傅立叶反变换变换到图像域，并将矢量转化为向量，得到修正好的图像。
本发明中判断K空间线是否需要进行逆向估计线间运动距离dx和dy的具体步骤如下(1)先将该未修正该K空间线的矩阵运用傅立叶变换的周期特性，变换成跟初始各矩阵一致的呈非共扼对称的矩阵，并将呈非共扼对称的矩阵通过傅立叶反变换变换到图像域，最后将矢量转换成向量而得图像，获得该图像的熵；(2)再在前一已修正好的K空间线的实际运动距离x_和y_的基础上加上该被修正K空间线前两个已修正好K空间线的实际运动距离之差，得到该被修正K空间线相对于Ky等于零的实际运动距离x和y，据此修正好该K空间线，所述的实际运动距离为相对于Ky等于零的运动距离；将修正后所得矩阵运用傅立叶变换的周期特性，变换成跟初始各矩阵一致的呈非共扼对称的矩阵，并将呈非共扼对称的矩阵通过傅立叶反变换变换到图像域，最后将矢量转换成向量而得图像，获得该图像的熵；(3)将上述步骤2中获得的熵值与步骤1中获得的未修正该K空间线时的熵值比较，如果熵值增大，则此K空间线需要进行逆向估计线间运动距离dx和dy。
本发明步骤5中所述逆向修正得出所求的线间运动距离dx和dy具体步骤如下(1)在所需修正的K空间线的前一已被修正好的K空间线的实际运动距离x和y的基础上加上估计的线间运动距离dx′和dy′，得到估计的所需修正的K空间线的实际运动距离x′和y′，对所需修正的K空间线用每次估计的实际运动距离x′和y′来进行修正；
(2)将上述每次估计修正后所得矩阵再次运用傅立叶变换的周期特性，变换成跟初始各矩阵一致的呈非共扼对称的矩阵，并将呈非共扼对称的矩阵通过傅立叶反变换变换到图像域，最后将矢量转换成向量而得图像，求出所得图像的熵；在所得所有的熵中找出最小值，该最小熵值所对应的估计的线间运动距离dx′和dy′就是所求的线间运动距离dx和dy。
为了使得估计的线间运动距离dx′和dy′估计得更加精确，本发明所述得出线间运动距离dx和dy的过程中，在修正需要逆向估计线间运动距离dx和dy的K空间线的同时，沿坐标轴向正方向对紧接该K空间线的余下未修正的K空间线中依次选取一部分也分别同时进行估计修正，估计修正时每个选取的K空间线的估计的实际运动距离分别为前一K空间线的估计的实际运动距离加上与估计的线间运动距离dx′和dy′，对估计修正每个选取的的K空间线用各自每次估计的实际运动距离x′和y′来进行估计修正。这样，在每次逆向迭代修正的过程中，就相当于有多个K空间线同时被估计修正了，以便使熵值之间的差值增大，免受积累误差和计算误差的影响，从而使得最小熵值所对应的估计的线间运动距离dx′和dy′更加接近所求的线间运动距离dx和dy的实际值。
本发明中所述根据所得实际运动距离修正K空间线的具体过程为该K空间线数据的幅值保持不变，相位为在步骤2中所得的该K空间线原相位的基础上加上一个值，该值为X方向的实际运动距离x和所在K空间坐标Kx的乘积与Y方向的实际运动距离y和所在K空间坐标Ky的乘积的和，再乘上一个常量C；其中常量C为2π/N，N为图像大小。
本发明中所述根据所得估计的实际运动距离修正K空间线的具体过程为该K空间线数据的幅值保持不变，相位为在步骤2中所得的该K空间线原相位的基础上加上一个值，该值为X方向估计的实际运动距离x′和所在K空间坐标Kx的乘积与Y方向估计的实际运动距离y′和所在K空间坐标Ky的乘积的和，再乘上一个常量C；其中常量C为2π/N，N为图像大小。
本发明中所述获得熵值的具体过程为得出每个象素点灰度值Bj与Bmax的比值，其中Bmax为图像所有像素的灰度值平方和开根而得，然后用每个象素算得的比值乘以该比值的以e为底的对数值，将图像中每个象素所得的乘积相加，对所得的和取负值就得到该图像的熵。
本发明方法与现有技术相比具有以下的优点1、由于本发明方法能自动搜索出运动轨迹，因此不需知道运动的先验信息；即不需知道运动的形式和过程 2、由于本发明方法是对任意方向进行逆向的估计，因此不管是周期的还是任意的都可估计出来；3、由于本发明方法是对整个图像区域进行自动修正，不用人工给出感兴趣区(ROI)或人工勾画出近似模板图像。


图1本发明的流程图；图2模板原始图像；图3模拟刚性平移运动后重建出的伪影图像；图4运用本发明所提方法消除运动伪影后重建出的图像。
具体实施例方式实施例一含有刚性平移运动伪影的磁共振图像的原始频域矢量矩阵数据可从磁共振成像系统中直接获得，也可将原始模板图像根据MR半傅立叶成像策略模拟刚性平移运动来获得。如图2，大小为256*256的原始模板图像根据MR半傅立叶成像策略模拟刚性平移运动得到图3的伪影图像的频域矢量矩阵数据，模拟过程中x和y方向的运动范围在[-10，10](单位像素)之间，对伪影图像的修正依次按如下步骤1、将所得频域矢量矩阵数据以频率编码方向为Kx方向，相位编码方向为Ky方向，把相位编码梯度Gy为零时定为Ky等于零确定K空间坐标；运用傅立叶变换的周期特性，将所获各原始频域矢量矩阵转换成在K空间坐标中以相位编码方向Ky等于零为对称轴的上下对称即的矩阵数据；并获得呈矢量共扼对称的矩阵中每个频域矢量的相位，即K空间坐标中每个K空间线的原相位；2、以K空间坐标中Ky等于零为参照K空间线，沿坐标轴向正方向的第一个K空间线开始，逐个修正每一个K空间线已偏移的相位；
3、修正沿坐标轴向正方向的第一个K空间线，该第一个K空间线必须要逆向估计线间运动距离dx和dy，采用逆向修正得出所求的线间运动距离dx和dy，具体过程如下在所需修正的第一个K空间线的前一已被修正好的K空间线的实际运动距离x和y的基础上加上估计的线间运动距离dx′和dy′，对应第一个所需修正的K空间线时前一已被修正好的K空间线的实际运动距离x_和y_为零，取dx′和dy′的估计值分别为[-3，-2，-1，0，1，2，3]，得到7*7共49个估计的所需修正的第一个K空间线的估计的实际运动距离x′和y′，对所需修正的第一个K空间线用每次估计的实际运动距离x′和y′进行估计修正，根据所得估计的实际运动距离修正K空间线的具体过程为该K空间线数据的幅值保持不变，相位为在步骤2中所得的该K空间线原相位的基础上加上一个值，该值为X方向估计的实际运动距离x′和所在K空间坐标Kx的乘积与Y方向估计的实际运动距离y′和所在K空间坐标Ky的乘积的和，再乘上一个常量C；其中常量C为2π/N，N为图像大小。
将每次估计修正后所得矩阵再次运用傅立叶变换的周期特性，变换成跟初始各矩阵一致的呈非共扼对称的矩阵，并将呈非共扼对称的矩阵通过傅立叶反变换变换到图像域，最后将矢量转换成向量而得图像，求出49次估计修正后对应所得的49个估计修正图像的熵。在所得所有的熵中找出最小值，它所对应的估计的线间运动距离dx′和dy′就是所求的线间运动距离dx和dy。然后根据此线间运动距离dx和dy与前一已被修正好的K空间线的实际运动距离x和y(零)，得出该被修正第一个K空间线的实际运动距离x和y，据此修正该第一个K空间线。
根据所得实际运动距离x和y修正K空间线的具体过程为该K空间线数据的幅值保持不变，相位为在步骤2中所得的该K空间线原相位的基础上加上一个值，该值为X方向的实际运动距离x和所在K空间坐标Kx的乘积与Y方向的实际运动距离y和所在K空间坐标Ky的乘积的和，再乘上一个常量C；其中常量C为2π/N，N为图像大小256。
4、对要修正的第二个K空间线，首先判断它是否需要进行估计线间运动距离dx和dy，这里dx和dy为两个K空间线线间的运动距离，即所需修正的K空间线相对于前一已修正好的K空间线的运动距离。
判断步骤如下先将该未修正的第二个K空间线的矩阵运用傅立叶变换的周期特性，变换成跟初始各矩阵一致的呈非共扼对称的矩阵，并将呈非共扼对称的矩阵通过傅立叶反变换变换到图像域，最后将矢量转换成向量而得图像，计算该图像的熵；再在前一已修正好的K空间线的实际运动距离x_和y_(此运动距离是相对于Ky等于零的运动距离)的基础上加上该被修正K空间线前两个已修正好K空间线的运动实际运动距离之差，得到该被修正K空间线相对于Ky等于零的实际运动距离x和y，据此修正好该第二个K空间线，将所得矩阵运用傅立叶变换的周期特性，变换成跟初始各矩阵一致的呈非共扼对称的矩阵，并将呈非共扼对称的矩阵通过傅立叶反变换变换到图像域，最后将矢量转换成向量而得图像，计算该图像的熵；上述获得熵值的具体过程为得出每个象素点灰度值Bj与Bmax的比值，其中Bmax为图像所有像素的灰度值平方和开根而得，然后用每个象素算得的比值乘以该比值的以e为底的对数值，将图像中每个象素所得的乘积相加，对所得的和取负值就得到该图像的熵。
将上述得到的两个修正前后熵值比较，如果熵值增大，则此第二K空间线需要进行逆向估计线间运动距离dx和dy，其修正过程重复上述步骤3即可。
5、如果熵值减小，修正不需要估计线间运动距离dx和dy的K空间线，可直接把前一已修正好的K空间线的实际运动距离x_和y_(即上述已修正好的第一K空间线的实际运动距离)加上该被修正K空间线前两个已修正好K空间线的运动实际运动距离之差(即上述已修正好的第一K空间线的实际运动距离与参照空间线的零运动距离之差)作为该被修正K空间线的实际运动距离x和y，然后据此修正该第二K空间线。
根据所得实际运动距离修正K空间线的具体过程为该K空间线数据的幅值保持不变，相位为在步骤2中所得的该K空间线原相位的基础上加上一个值，该值为X方向的实际运动距离x和所在K空间坐标Kx的乘积与Y方向的实际运动距离y和所在K空间坐标Ky的乘积的和，再乘上一个常量C；其中常量C为2π/N，N为图像大小256。
6、重复步骤4、5中的处理过程，修正全部沿坐标轴向正方向的K空间线，对Ky等于零为参照的向负方向的各个K空间线的修正，以Ky等于零为对称轴取正方向对称行的共扼矢量即可；7、将所有修正好的K空间线的频域矢量运用傅立叶变换的周期特性，变换成跟初始矩阵一致的非共扼对称矩阵，然后通过傅立叶反变换变换到图像域，并将矢量转化为向量，得到修正好的图像。
实施例二步骤1、2与实施例一相同；步骤3，为了使估计出的线间运动距离dx′和dy′估计得更加精确，根据运动一致性和连续性，在上述实施例一的步骤3中修正需要逆向估计线间运动距离dx和dy的第一K空间线的同时，对余下未修正的K空间线中紧接该第一个K空间线沿坐标轴向正方向依次选取即将被修正的七个K空间线也分别依次进行估计修正，选取余下未修正的K空间线的多少可任意，以满足图像的熵值变化不受积累误差和计算误差的影响为准。估计修正时每个选取的K空间线的估计的实际运动距离分别为前一K空间线的估计的实际运动距离加上估计的线间运动距离dx′和dy′，对估计修正每个选取的的K空间线用各自每次估计的实际运动距离x′和y′来进行修正。这样加上该被修正的第一个K空间线为8个，根据正负共扼对称，负方向也为8个，所以每次估计修正时同时对16个K空间线进行了估计修正。经上述每次估计修正后相应得到的49个矩阵在变换成图像后，所得的熵值之间的差值会变大，免受积累误差和计算误差的影响，从而使得最小熵值所对应的估计的线间运动距离dx′和dy′更加接近所求的线间运动距离dx和dy的实际值。
然后根据上述得出的更加接近实际值的线间运动距离dx和dy与前一已被修正好的K空间线的实际运动距离x和y(零)，得出该被修正第一个K空间线的实际运动距离x和y，据此修正该第一个K空间线。
余下步骤的处理过程与实施例一相同，得到修正好的图像(如图4)。
权利要求
1.一种磁共振图像刚性平移运动伪影的自动逆向修正方法，其特征在于依次包括以下步骤(1)获取含有刚性平移运动伪影的磁共振图像的原始频域矢量矩阵数据，即K空间数据；(2)以频率编码方向为Kx方向，相位编码方向为Ky方向，把相位编码梯度Gy为零时定为Ky等于零确定K空间坐标；运用傅立叶变换的周期特性，将所获各原始频域矢量矩阵转换成在K空间坐标中以相位编码方向Ky等于零为对称轴的上下共扼对称的矩阵数据；并获得呈矢量共扼对称的矩阵中每个频域矢量的相位，即K空间坐标中每个K空间线的原相位；(3)以K空间坐标中Ky等于零为参照K空间线，沿坐标轴向正方向的第一个K空间线开始，逐个修正每一个K空间线已偏移的相位；(4)除设定为必须要逆向估计线间运动距离的第一个K空间线之外，对要修正的K空间线，首先判断该K空间线是否需要进行逆向估计线间运动距离dx和dy，所述的线间运动距离dx和dy为两个K空间线线间的运动距离，即所需修正的K空间线相对于前一已修正好的K空间线的运动距离；(5)若需要逆向估计线间运动距离dx和dy的K空间线时，采用逆向修正得出所求的线间运动距离dx和dy，然后根据此线间运动距离dx和dy加上前一已被修正好的K空间线的实际运动距离x_和y_得出该被修正K空间线的实际运动距离x和y，据此实际运动距离x和y修正该K空间线；(6)若不需要逆向估计线间运动距离dx和dy的K空间线时，可直接把前一已修正好的K空间线的实际运动距离x_和y_加上该被修正K空间线前两个已修正好K空间线的实际运动距离之差作为该被修正K空间线的实际运动距离x和y，然后据此实际运动距离x和y修正该K空间线；(7)重复步骤4、5、6中的处理过程，修正全部沿坐标轴向正方向的K空间线，对Ky等于零为参照的向负方向的各个K空间线的修正，以Ky等于零为对称轴取正方向对称行的共扼矢量即可；(8)将所有修正好的K空间线的频域矢量运用傅立叶变换的周期特性，变换成跟初始矩阵一致的非共扼对称矩阵，然后通过傅立叶反变换变换到图像域，并将矢量转化为向量，得到修正好的图像。
2.根据权利要求1所述的磁共振图像刚性平移运动伪影的自动逆向修正方法，其特征在于在步骤4中判断K空间线是否需要进行逆向估计线间运动距离dx和dy的具体步骤如下(1)先将该未修正该K空间线的矩阵运用傅立叶变换的周期特性，变换成跟初始各矩阵一致的呈非共扼对称的矩阵，并将呈非共扼对称的矩阵通过傅立叶反变换变换到图像域，最后将矢量转换成向量而得图像，获得该图像的熵；(2)再在前一已修正好的K空间线的实际运动距离x_和y_的基础上加上该被修正K空间线前两个已修正好K空间线的实际运动距离之差，得到该被修正K空间线相对于Ky等于零的实际运动距离x和y，据此修正好该K空间线，所述的实际运动距离为相对于Ky等于零的运动距离；将修正后所得矩阵运用傅立叶变换的周期特性，变换成跟初始各矩阵一致的呈非共扼对称的矩阵，并将呈非共扼对称的矩阵通过傅立叶反变换变换到图像域，最后将矢量转换成向量而得图像，获得该图像的熵；(3)将上述步骤2中获得的熵值与步骤1中获得的未修正该K空间线时的熵值比较，如果熵值增大，则此K空间线需要进行逆向估计线间运动距离dx和dy。
3.根据权利要求1所述的磁共振图像刚性平移运动伪影的自动逆向修正方法，其特征在于步骤5中所述逆向修正得出所求的线间运动距离dx和dy具体步骤如下(1)在所需修正的K空间线的前一已被修正好的K空间线的实际运动距离x和y的基础上加上估计的线间运动距离dx′和dy′，得到估计的所需修正的K空间线的估计的实际运动距离x′和y′，对所需修正的K空间线用每次估计的实际运动距离x′和y′来进行修正；(2)将上述每次修正后所得矩阵再次运用傅立叶变换的周期特性，变换成跟初始各矩阵一致的呈非共扼对称的矩阵，并将这些呈非共扼对称的矩阵通过傅立叶反变换变换到图像域，最后将矢量转换成向量而得图像，求出所得图像的熵；在所得所有的熵中找出最小值，该最小熵值所对应的估计的线间运动距离dx′和dy′就是所求的线间运动距离dx和dy。
4.根据权利要求3所述的磁共振图像刚性平移运动伪影的自动逆向修正方法，其特征在于所述得出线间运动距离dx和dy的过程中，在修正需要逆向估计线间运动距离dx和dy的K空间线的同时，沿坐标轴向正方向对紧接该K空间线的余下未修正的K空间线中依次选取一部分也分别依次进行估计修正，估计修正时每个选取的K空间线的估计的实际运动距离分别为前一K空间线的估计的实际运动距离加上估计的线间运动距离dx′和dy′，对估计修正每个选取的的K空间线用各自每次估计的实际运动距离x′和y′来进行估计修正。
5.根据权利要求1或2所述的磁共振图像刚性平移运动伪影的自动逆向修正方法，其特征在于所述根据所得实际运动距离修正K空间线的具体过程为该K空间线数据的幅值保持不变，相位为在步骤2中所得的该K空间线原相位的基础上加上一个值，该值为X方向的实际运动距离x和所在K空间坐标Kx的乘积与Y方向的实际运动距离y和所在K空间坐标Ky的乘积的和，再乘上一个常量C；其中常量C为2π/N，N为图像大小。
6.根据权利要求3或4所述的磁共振图像刚性平移运动伪影的自动逆向修正方法，其特征在于所述根据所得估计的实际运动距离修正K空间线的具体过程为该K空间线数据的幅值保持不变，相位为在步骤2中所得的该K空间线原相位的基础上加上一个值，该值为X方向估计的实际运动距离x′和所在K空间坐标Kx的乘积与Y方向估计的实际运动距离y′和所在K空间坐标Ky的乘积的和，再乘上一个常量C；其中常量C为2π/N，N为图像大小。
7.根据权利要求2或3所述的磁共振图像刚性平移运动伪影的自动逆向修正方法，其特征在于所述获得熵值的具体过程为得出每个象素点灰度值Bj与Bmax的比值，其中Bmax为图像所有像素的灰度值平方和开根而得，然后用每个象素算得的比值乘以该比值的以e为底的对数值，将图像中每个象素所得的乘积相加，对所得的和取负值就得到该图像的熵。
全文摘要
本发明公开了一种磁共振图像刚性平移运动伪影的自动逆向修正方法，获取伪影图像的原始频域矢量矩阵数据，并确定K空间坐标；以Ky等于零为参照，从第一个K空间线开始，逐个修正每一个K空间线偏移的相位；对要修正的K空间线，确定是否需要进行逆向估计线间运动距离，若需要，采用逆向修正得出实际运动距离来修正该K空间线；若不需要，可直接得出其实际运动距离来修正该K空间线；最后将所有修正好的K空间线变换到图像域，得到修正好的图像。本发明使得因刚性平移而损害的数据经修正后能重建出不含运动伪影的高质量的图像。
文档编号G01R33/56GK1515912SQ03140278
公开日2004年7月28日 申请日期2003年8月27日 优先权日2003年8月27日
发明者侯正松, 江贵平, 陈武凡 申请人:中国人民解放军第一军医大学