CN108051766A - 基于奇异谱分析的speed快速磁共振成像方法 - Google Patents

基于奇异谱分析的speed快速磁共振成像方法 Download PDF

Info

Publication number
CN108051766A
CN108051766A CN201711337327.8A CN201711337327A CN108051766A CN 108051766 A CN108051766 A CN 108051766A CN 201711337327 A CN201711337327 A CN 201711337327A CN 108051766 A CN108051766 A CN 108051766A
Authority
CN
China
Prior art keywords
ghost
data
image
low
space
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201711337327.8A
Other languages
English (en)
Inventor
金朝阳
向清三
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hangzhou Dianzi University
Original Assignee
Hangzhou Dianzi University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hangzhou Dianzi University filed Critical Hangzhou Dianzi University
Priority to CN201711337327.8A priority Critical patent/CN108051766A/zh
Publication of CN108051766A publication Critical patent/CN108051766A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/4818—MR characterised by data acquisition along a specific k-space trajectory or by the temporal order of k-space coverage, e.g. centric or segmented coverage of k-space
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/5602—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution by filtering or weighting based on different relaxation times within the sample, e.g. T1 weighting using an inversion pulse
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/561—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution by reduction of the scanning time, i.e. fast acquiring systems, e.g. using echo-planar pulse sequences

Landscapes

  • Physics & Mathematics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Condensed Matter Physics & Semiconductors (AREA)
  • General Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Engineering & Computer Science (AREA)
  • Signal Processing (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于奇异谱分析的SPEED快速磁共振成像方法,本发明包括k空间数据采集、填零重建、奇异谱分析重建、鬼影定位、基于双层鬼影模型的SPEED图像重建，采用本发明方法，通过对k空间中心数据进行奇异谱分析重建来抑制k空间中心数据直接填零重建带来的吉布斯环状伪影，基于不同分辨率图像都具有相位成片连贯性的特点，通过合成高质量的定位图像来提高鬼影定位的精度，进而提高了SPEED成像质量；通过鬼影定位方法，可将常规SPEED快速成像的数据采集组数减少三分之一，进一步提高了SPEED成像速度。

Description

基于奇异谱分析的SPEED快速磁共振成像方法
技术领域
本发明属于磁共振的图像成像领域，涉及一种基于奇异谱分析的SPEED快速磁共振成像方法。
背景技术
SPEED(Skipped Phase Encoding and Edge Deghosting)成像技术是一种通过在相位编码(Phase Encoding,PE)方向上减少数据点来缩短采集时间的快速磁共振成像(Magnetic Resonance Imaging,MRI)方法(QS Xiang,Accelerating MRI by skippedphase encoding and edge deghosting(SPEED)，Magnetic Resonance in Medicine,53:1112-1117,2005)。SPEED通过在k空间进行简单规则的欠采样来采集数据，基于差分图像的稀疏性和双层鬼影模型，可通过解析法来重建图像，其解析求解过程非常快速。SPEED成像方法不但易于实现，而且也易于和现有的采集方式结合，是一种很有应用潜力的MRI快速成像方法。
目前已申请的关于SPEED快速成像方面的MRI专利有：基于小波域稀疏表示的SPEED快速磁共振成像方法(授权号：ZL 2013102071971.1)，提出基于小波域的数据稀疏特性来提高SPEED快速成像方法的成像质量。基于离散余弦变换的SPEED快速磁共振成像方法(授权号：ZL201310719667.2)，提出基于离散余弦变换来提高SPEED快速成像方法的成像质量。基于k空间中心鬼影定位的SPEED磁共振成像方法(申请号：201610920503.X)，提出一种通过利用k空间中心数据来进行鬼影定位，进一步减小SPEED数据采集组数，加快SPEED成像速度的方法。
目前还未能查询到任何基于奇异谱分析的SPEED快速磁共振成像方法的授权发明专利或申请。
国内外已发表的关于SPEED成像方面的文章有：2016年，金朝阳等人提出了基于离散余弦变换和离散小波变换的SPEED快速成像方法(Jin Z,Ye H,Du YP,XiangQS.Improving image quality for skipped phase encoding and edge deghosting(SPEED)by exploiting several sparsifying transforms.Magnetic Resonance inMedicine.75:2031-2045，2016)，该方法利用离散余弦变换和离散小波变换对数据进行稀疏表示，相比于常规SPEED方法采用基于离散差分变换的数据稀疏表示方法，该方法获得了更好的成像质量。2013年，金朝阳和向清三提出了通用G-SPEED(General-SPEED)采样方法(Jin Z,Xiang QS.Accelerated MRI by SPEED with generalized samplingschemes.Magnetic Resonance in Medicine.70:1674-1681，2013)，G-SPEED方法突破了传统SPEED方法的采样间隔周期N必须是质数(例如：N＝5、7、11)的限制，通过秩判据的方式，使得N不但可为质数，也可为合数(例如：N＝2、4、6、8、9)。2009年，常征等人提出EMA-SPEED(Efficient Multiple Acquisition by SPEED)算法(Chang Z,Xiang QS,Ji J,and YinFF.Efficient multiple acquisitions by skipped phase encoding and edgedeghosting(SPEED)using shared spatial information.Magnetic Resonance inMedicine.61:229-233,2009)，通过共享多个采集间的相似空间信息进一步缩短了SPEED的数据采集时间，从而可获得比单次采集更高的加速比。2007年，基于MRA数据本身就非常稀疏的特性，常征和向清三将SPEED的双层模型简化到单层模型(Chang Z and XiangQS.Simplified skipped phase encoding and edge deghosting(SPEED)for imagingsparse objects with applications to MRA.Med Phys.34:3173-3182,2007)，提出了S-SPEED(Simplified-SPEED)算法，该算法适用于数据本身就非常稀疏的场合，例如暗背景亮信号的MRA应用(Chang Z,Xiang QS,Shen H and Yin FF.Accelerating non-contrast-enhanced MR angiography with inflow inversion recovery imaging by skippedphase encoding and edge deghosting(SPEED).Journal of Magnetic ResonanceImaging.31:757-765,2010)。2006年，常征和向清三将SPEED算法与并行成像技术进一步结合，提出了SPEED-ACE成像法(Chang Z and Xiang QS.Highly accelerated MRI byskipped phase encoding and edge deghosting with array coil enhancement(SPEED-ACE).Med Phys.33:3758-3766,2006)，通过采用多个采集线圈来共同采集k-空间欠采样数据，从而提高成像速度。
以上发表的关于SPEED快速成像方面的文章或已授权的发明专利，是基于三组欠采样数据来定位鬼影或基于两组欠采样数据和低分辨率图像来定位鬼影，再进行图像重建的方法，还未公开过任何基于奇异谱分析的SPEED快速磁共振成像方法。
发明内容
本发明针对现有SPEED技术的不足，通过对k空间中心数据进行奇异谱分析重建来抑制k空间中心数据直接填零重建成低分辨率图像中的吉布斯环状伪影(Gibbs RingingArtifacts)，基于不同分辨率图像都具有相位成片连贯性的特点，本发明通过合成高质量的定位图像来提高鬼影定位的精度，进而提高了SPEED成像质量；通过鬼影定位方法，可将常规SPEED快速成像的数据采集组数减少三分之一，进一步提高了SPEED成像速度；本发明主要包括五个步骤：k空间数据采集、填零重建、奇异谱分析重建、鬼影定位、基于双层鬼影模型的SPEED图像重建。
步骤1：k空间数据采集
在k空间的相位编码方向，即PE方向，每隔N行采集一行数据，共采集两组，分别用S1和S2表示，用kx表示在k空间频率编码FE(Frequency Encoding)方向的位置，ky表示在PE方向的位置，则S1(kx,ky)和S2(kx,ky)分别表示S1和S2中的一个数据点的值，用d1,d2表示每组欠采样数据在PE方向上的偏移量，采样方式可用N(d1,d2)表示。根据图像大小，在PE方向的k空间中心区域采集16至64行数据，用Sc表示，Sc(kx,ky)表示Sc中的一个数据点。在FE方向的数据都为全采集。
步骤2：填零重建
填零重建包含3个步骤：两组欠采样数据的填零重建、差分变换、k空间中心数据填零重建成低分辨率图像。
步骤2-1：两组欠采样数据的填零重建
对于两组欠采样数据S1和S2，其对应k空间中没有进行数据采集的点用0表示，进行常规的填零傅立叶重建，重建后图像分别用I1和I2表示，用(x,y)表示图像空间中每个像素点的坐标值，则I1(x,y)和I2(x,y)表示I1和I2中的一个像素点的值。由于k空间中每隔N行采集一行数据，因此每组数据对应的填零傅立叶重建图像中有N层重叠的鬼影，每个像素点上最多有N层重叠的鬼影。例如，当N＝5时，I1和I2上分别有5层重叠的鬼影。
步骤2-2：差分变换
对步骤2-1得到的图像I1和I2分别进行差分变换，得到稀疏的边缘鬼影图像E1和E2。稀疏的鬼影图像E1和E2中的每个像素点E1(x,y)和E2(x,y)可分别表示为：
E1(x,y)＝I1(x,y)-I1(x,y-1) [1]
E2(x,y)＝I2(x,y)-I2(x,y-1) [2]
差分变换后，E1(x,y)和E2(x,y)上通常只有两层重叠的鬼影。
步骤2-3：k空间中心数据填零重建成低分辨率图像
采集到的k空间中心部分数据Sc，其对应k空间中没有进行数据采集的点用0表示，进行常规的填零傅立叶重建，形成一个低分辨率的重建图像Ilow，其对应的每一个像素点用Ilow(x,y)表示。
步骤3：奇异谱分析重建
奇异谱分析重建包含8个步骤：差分变换、确立采集矩阵、寻找奇异点、层析法更新差分图像、构建奇异谱函数、计算奇异度、重建出k空间缺失数据、基于离散傅立叶逆变换的图像重建。
步骤3-1：差分变换
对填零重建图像Ilow在y方向(对应于k空间中的ky方向)上进行差分变换，得到稀疏的差分图像ΔIlow，对应的每个像素点可表示为ΔIlow(x,y)：
ΔIlow(x,y)＝Ilow(x,y)-Ilow(x,y-1) [3]
阈值T设为差分图像ΔIlow中幅值最大值的十分之一。
步骤3-2：确立采集矩阵
采集矩阵Mask中的每一点Mask(kx,ky)表示为：
采集矩阵对应的频域点表示为：
i(x,y)＝IDFT(Mask(kx,ky)) [5]
步骤3-3：寻找奇异点
取ΔIlow中幅值最大的点为奇异点j(xj,yj)，加入奇异点集SP中。把i(x,y)矩阵的原点移动至奇异点j(xj,yj)处，更新采集矩阵的频域表示为i(x-xj,y-yj)。
步骤3-4：层析法更新差分图像
更新后的差分图像ΔIlow(x,y)new可表示为，
ΔIlow(x,y)new＝ΔIlow(x,y)-α·i(x-xj,y-yj) [6]
α＝ΔIlow(x,y)||Mask||/(M×N) [7]
其中M和N表示x和y方向的像素点数，Mask表示采集矩阵中非零点的总数。
如果ΔIlow(x,y)new≥T，则令ΔIlow(x,y)＝ΔIlow(x,y)new，转到步骤3-3继续寻找奇异点；如果ΔIlow(x,y)new＜T，则停止奇异点的寻找，奇异点集SP的奇异点的总个数为Q，继续执行步骤3-5。
步骤3-5：构建奇异谱函数
对奇异点集SP中的奇异点分别构建奇异谱函数Wj(kx,ky)，用公式表示为：
Wj(kx,ky)＝IDFT(δ(y-yj)u(x-xj)) [8]
其中，
IDFT表示离散傅立叶逆变换；
步骤3-6：计算奇异度
利用已采集到的k空间数据Sc通过求解方程[11]来得到各个奇异点对应的奇异度aj，
步骤3-7：重建出k空间缺失数据
用奇异谱函数和奇异度来估算出k空间中没有进行过数据采集处的数据Sr(kx,ky)：
步骤3-8：基于离散傅立叶逆变换的图像重建
对于k空间中进行过数据采集的点用实际采集到的数据Sc(kx,ky)，没有进行过采集的点用步骤3.7估算得到的Sr(kx,ky)，然后利用离散傅立叶逆变换进行图像重建，得到奇异谱重建图像ISFA，其对应的像素点可表示为，
步骤4：鬼影定位
鬼影定位包含4个步骤：合成定位图像、差分变换、建立重叠鬼影图和确立鬼影阶数。
步骤4-1：合成定位图像
将步骤2-3得到的低分辨率图像Ilow与步骤3-8得到的重建图像ISFA进行对应像素点之间的相乘操作，合成用于鬼影定位的定位图像Ic，对应的像素点可表示为：
Ic(x,y)＝Ilow(x,y)×ISFA(x,y) [14]
步骤4-2：差分变换
对Ic进行差分变换，得到稀疏的边缘鬼影图像Ec。
步骤4-3：建立重叠鬼影图
在相位编码方向对Ec分别进行长度为Ny×n/N的平移，其中Ny表示沿PE方向的数据矩阵的大小，n表示边缘鬼影的阶数(不同的阶表示鬼影位置不同)，n＝0,1,2,…,N-1。这n个稀疏的边缘鬼影相加后形成一个重叠的鬼影映射图Ec,n。
步骤4-4：确立鬼影阶数
在鬼影映射图Ec,n中，为每个像素点找出两个最强的鬼影，并记录下它们对应的鬼影阶数(n1,n2)。
步骤5：基于双层鬼影模型的SPEED图像重建
基于双层鬼影模型的SPEED图像重建包含4个步骤：双层鬼影模型求解、重叠鬼影的分离、多个分离鬼影图像的配准求和、逆滤波重建。
步骤5-1：双层鬼影模型求解
稀疏边缘鬼影图像E1和E2中，由于每个像素点上通常只有两层鬼影的重叠，因此可采用双层稀疏边缘鬼影模型来描述，表示为：
公式[15]中为相位因子，Gn1和Gn2分别为每个像素点上需要确定的不同阶的鬼影，n1和n2分别表示不同的鬼影阶数。定义为：
公式[16]中d表示每组欠采样数据在PE方向上的偏移量d1和d2，n为鬼影阶数。
在公式[15]中，由于E1、E2、d和N已知，步骤3-3得到了鬼影阶数(n1,n2)，即公式[15]中的两个方程，仅有两个未知数Gn1和Gn2，因此可直接解出公式[15]中的两个重叠的鬼影Gn1和Gn2。
步骤5-2：重叠鬼影的分离
对步骤5-1得到的Gn1和Gn2中的像素点，按不同的鬼影阶数n进行分类，产生N个分离的鬼影映射图Gn,其中n＝0,1,…,N-1。
步骤5-3：多个分离鬼影图像的配准求和
步骤5-2得到的N个鬼影映射图Gn，各自对应的鬼影位置不同，可通过移位和对齐来配准。配准后各鬼影图对应的像素点求和后，得到没有重叠鬼影的边缘映射图像E0。
步骤5-4：逆滤波重建
步骤5-3得到的边缘映射图像E0经离散傅立叶变换(DFT)到k空间，其对应k空间中实际进行数据采集的点的值用实际采集的数据替代，得到k空间数据R0。基于逆滤波公式[17]重建出最终的SPEED图像I0。
公式[17]中ky表示沿PE方向的k空间位置。
采用本发明方法，即基于奇异谱分析合成的定位图可以提高重叠鬼影定位的精度，从而在缩短SPEED数据采集时间的情况下，提高了SPEED重建图的质量。同时本发明具有以下特点：
(1)本发明基于奇异谱分析法恢复出k空间中心区域外的数据，避免了k空间中心区域直接填零带来的信号截断现象，因此相比k空间中心区域填零重建形成的低分辨率图像，基于奇异谱分析法重建的图像能有效抑制吉布斯环状伪影。
(2)本发明基于低分辨率图像和奇异谱分析法重建图像合成的定位图像，不但能抑制吉布斯环状伪影，还能有效减少高频噪声。高质量的定位图像能提高鬼影定位的精度，从而提高最终的SPEED成像质量。
(3)本发明提高了SPEED成像方法的数据采集速度。常规SPEED技术需采集三组k空间欠采样数据，而本发明通过合成的定位图像来计算鬼影阶数，只需采集两组k空间欠采样数据即可进行SPEED双层鬼影模型的求解，提高了SPEED成像方法的数据采集速度。
附图说明
图1是本发明数据采集方式的示意图；
图2基于奇异谱分析的SPEED快速磁共振成像方法示意图；
图3是采用本发明进行数据采集和图像重建实例的结果图。
具体实施方式
以下结合附图对本发明作进一步说明。
本发明主要包括五个步骤：k空间数据采集、填零重建、奇异谱分析重建、鬼影定位、基于双层鬼影模型的SPEED图像重建。
步骤1：k空间数据采集
如图1所示，在k空间的相位编码方向，即PE方向，每隔N行采集一行数据，共采集两组，分别用S1和S2表示，用kx表示在k空间频率编码FE(Frequency Encoding)方向的位置，ky表示在PE方向的位置，则S1(kx,ky)和S2(kx,ky)分别表示S1和S2中的一个数据点的值，用d1,d2表示每组欠采样数据在PE方向上的偏移量，采样方式可用N(d1,d2)表示。根据图像大小，在PE方向的k空间中心区域采集16至64行数据，用Sc表示，Sc(kx,ky)表示Sc中的一个数据点。在FE方向的数据都为全采集。
步骤2：填零重建
填零重建包含3个步骤：两组欠采样数据的填零重建、差分变换、k空间中心数据填零重建成低分辨率图像。
步骤2-1：两组欠采样数据的填零重建
对于两组欠采样数据S1和S2，其对应k空间中没有进行数据采集的点用0表示，进行常规的填零傅立叶重建，重建后图像分别用I1和I2表示，用(x,y)表示图像空间中每个像素点的坐标值，则I1(x,y)和I2(x,y)表示I1和I2中的一个像素点的值。由于k空间中每隔N行采集一行数据，因此每组数据对应的填零傅立叶重建图像中有N层重叠的鬼影，每个像素点上最多有N层重叠的鬼影。例如，当N＝5时，I1和I2上分别有5层重叠的鬼影，如图2所示。
步骤2-2：差分变换
如图2所示，对步骤2-1得到的图像I1和I2分别进行差分变换，得到稀疏的边缘鬼影图像E1和E2。稀疏的鬼影图像E1和E2中的每个像素点E1(x,y)和E2(x,y)可分别表示为：
E1(x,y)＝I1(x,y)-I1(x,y-1) [1]
E2(x,y)＝I2(x,y)-I2(x,y-1) [2]
差分变换后，E1(x,y)和E2(x,y)上通常只有两层重叠的鬼影。
步骤2-3：k空间中心数据填零重建成低分辨率图像
采集到的k空间中心部分数据Sc，其对应k空间中没有进行数据采集的点用0表示，进行常规的填零傅立叶重建，形成一个低分辨率的重建图像Ilow，其对应的每一个像素点用Ilow(x,y)表示，如图2所示。
步骤3：奇异谱分析重建
奇异谱分析重建包含8个步骤：差分变换、确立采集矩阵、寻找奇异点、层析法更新差分图像、构建奇异谱函数、计算奇异度、重建出k空间缺失数据、基于离散傅立叶逆变换的图像重建。
步骤3-1：差分变换
如图2所示，对填零重建图像Ilow在y方向(对应于k空间中的ky方向)上进行差分变换，得到稀疏的差分图像ΔIlow，对应的每个像素点可表示为ΔIlow(x,y)：
ΔIlow(x,y)＝Ilow(x,y)-Ilow(x,y-1) [3]
阈值T设为差分图像ΔIlow中幅值最大值的十分之一。
步骤3-2：确立采集矩阵
采集矩阵Mask中的每一点Mask(kx,ky)表示为：
采集矩阵对应的频域点表示为：
i(x,y)＝IDFT(Mask(kx,ky)) [5]
步骤3-3：寻找奇异点
取ΔIlow中幅值最大的点为奇异点j(xj,yj)，加入奇异点集SP中。把i(x,y)矩阵的原点移动至奇异点j(xj,yj)处，更新采集矩阵的频域表示为i(x-xj,y-yj)。
步骤3-4：层析法更新差分图像
更新后的差分图像ΔIlow(x,y)new可表示为，
ΔIlow(x,y)new＝ΔIlow(x,y)-α·i(x-xj,y-yj) [6]
α＝ΔIlow(x,y)||Mask||/(M×N) [7]
其中M和N表示x和y方向的像素点数，||Mask||表示采集矩阵中非零点的总数。
如果ΔIlow(x,y)new≥T，则令ΔIlow(x,y)＝ΔIlow(x,y)new，转到步骤3-3继续寻找奇异点；如果ΔIlow(x,y)new＜T，则停止奇异点的寻找，奇异点集SP的奇异点的总个数为Q，继续执行步骤3-5。
步骤3-5：构建奇异谱函数
对奇异点集SP中的奇异点分别构建奇异谱函数Wj(kx,ky)，用公式表示为：
Wj(kx,ky)＝IDFT(δ(y-yj)u(x-xj)) [8]
其中，
步骤3-6：计算奇异度
利用已采集到的k空间数据Sc通过求解方程[11]来得到各个奇异点对应的奇异度aj，
步骤3-7：重建出k空间缺失数据
用奇异谱函数和奇异度来估算出k空间中没有进行过数据采集处的数据Sr(kx,ky)：
步骤3-8：基于离散傅立叶逆变换的图像重建
如图2所示，对于k空间中进行过数据采集的点用实际采集到的数据Sc(kx,ky)，没有进行过采集的点用步骤3.7估算得到的Sr(kx,ky)，然后利用离散傅立叶逆变换进行图像重建，得到奇异谱重建图像ISFA，其对应的像素点可表示为，
步骤4：鬼影定位
鬼影定位包含4个步骤：合成定位图像、差分变换、建立重叠鬼影图和确立鬼影阶数。
步骤4-1：合成定位图像
如图2所示，将步骤2-3得到的低分辨率图像Ilow与步骤3-8得到的重建图像ISFA进行对应像素点之间的相乘操作，合成用于鬼影定位的定位图像Ic，对应的像素点可表示为：
Ic(x,y)＝Ilow(x,y)×ISFA(x,y) [14]
步骤4-2：差分变换
对Ic进行差分变换，得到稀疏的边缘鬼影图像Ec，如图2所示。
步骤4-3：建立重叠鬼影图
在相位编码方向对Ec分别进行长度为Ny×n/N的平移，其中Ny表示沿PE方向的数据矩阵的大小，n表示边缘鬼影的阶数(不同的阶表示鬼影位置不同)，n＝0,1,2,…,N-1。这n个稀疏的边缘鬼影相加后形成一个重叠的鬼影映射图Ec,n，如图2所示。
步骤4-4：确立鬼影阶数
在鬼影映射图Ec,n中，为每个像素点找出两个最强的鬼影，并记录下它们对应的鬼影阶数(n1,n2)，如图2所示。
步骤5：基于双层鬼影模型的SPEED图像重建
基于双层鬼影模型的SPEED图像重建包含4个步骤：双层鬼影模型求解、重叠鬼影的分离、多个分离鬼影图像的配准求和、逆滤波重建。
步骤5-1：双层鬼影模型求解
稀疏边缘鬼影图像E1和E2中，由于每个像素点上通常只有两层鬼影的重叠，因此可采用双层稀疏边缘鬼影模型来描述，表示为：
公式[15]中为相位因子，Gn1和Gn2分别为每个像素点上需要确定的不同阶的鬼影，n1和n2分别表示不同的鬼影阶数。定义为：
公式[16]中d表示每组欠采样数据在PE方向上的偏移量d1和d2，n为鬼影阶数，如图2所示。
在公式[15]中，由于E1、E2、d和N已知，步骤3-3得到了鬼影阶数(n1,n2)，即公式[15]中的两个方程，仅有两个未知数Gn1和Gn2，因此可直接解出公式[15]中的两个重叠的鬼影Gn1和Gn2。
步骤5-2：重叠鬼影的分离
对步骤5-1得到的Gn1和Gn2中的像素点，按不同的鬼影阶数n进行分类，产生N个分离的鬼影映射图Gn,其中n＝0,1,…,N-1。
步骤5-3：多个分离鬼影图像的配准求和
步骤5-2得到的N个鬼影映射图Gn，各自对应的鬼影位置不同，可通过移位和对齐来配准。配准后各鬼影图对应的像素点求和后，得到没有重叠鬼影的边缘映射图像E0，如图2所示。
步骤5-4：逆滤波重建
步骤5-3得到的边缘映射图像E0经离散傅立叶变换(DFT)到k空间，其对应k空间中实际进行数据采集的点的值用实际采集的数据替代，得到k空间数据R0。基于逆滤波公式[17]重建出最终的SPEED图像I0。
公式[17]中IDFT表示离散傅立叶逆变换，ky表示沿PE方向的k空间位置，如图2所示。
以下结合人体头部的MRI数据，对基于奇异谱分析的SPEED磁共振成像方法进行实例说明。假设要采集的MRI图像的矩阵大小为kx×ky＝256×256。首先以采样方式N(d1,d2)＝5(0,2)进行数据采集，在k空间的相位编码PE方向每隔N＝5行采集一行k空间数据，共采集两组，分别得到欠采样的k空间数据S1和S2。在信息量集中的k空间中心区域进行全采样，共采集32行相位编码数据。接下来，对两组欠采样数据S1和S2和k空间中心数据Sc分别进行常规的填零傅立叶重建，重建后图像分别为I1、I2和Ilow。然后对图像I1和I2分别进行差分变换，得到稀疏的鬼影图像E1和E2。然后进行奇异谱分析重建，先对Ilow进行差分变换，得到稀疏的ΔIlow图像，阈值T设为ΔIlow中幅值最大值的十分之一，通过确立采集矩阵、基于阈值T和层析法来寻找ΔIlow中的奇异点集SP，基于k空间中心的32行数据构建奇异谱函数并计算奇异度后，重建出Sc对应的k空间缺失数据，再基于离散傅立叶逆变换重建出奇异谱图像ISFA。将Ilow和ISFA图像的对应点相乘合成定位图像Ic图像，对进行差分变换，得到稀疏的边缘图Ec，在相位编码方向对Ec分别进行长度为256×n/5的平移，其中n＝0,1,2,3,4。这5个边缘鬼影相加后生成一个重叠的鬼影映射图Ec,n＝Ec,5。在鬼影映射图Ec,5中，为每个像素点找出两个最强的鬼影，并记录下它们对应的鬼影阶数(n1,n2)。基于鬼影阶数(n1,n2)和稀疏鬼影图E1和E2，可直接解出双层鬼影模型中的两个主要的重叠鬼影Gn1和Gn2。根据不同的鬼影阶数值，对Gn1和Gn2中的像素点分类，产生5个分离的鬼影映射图像Gn，n＝0,1,2,3,4。这5个鬼影映射图像经过移位、配准和求和后得到无重叠鬼影的边缘映射图E0。E0经过DFT变换到k空间，然后用所有实际采集到的k空间数据替换部分k空间数据，再经过逆滤波公式[17]重建出最终的SPEED图像I0。
如图3所示，图3(a)为人体头部的全采样参考图，图3(b)为Sc图像，图3(c)为稀疏鬼影图Ilow，图3(d)ΔIlow图像，可以看出有明显的吉布斯环状伪影，图3(e)基于低分辨率图像来定位重叠鬼影的SPEED方法的重建结果图，图3(f)为图3(e)与图3(a)相比较的误差映射图，图3(g)为奇异谱分析重建恢复的Sc的k空间分布，图3(h)为奇异点集，图3(i)为Ec，可以看出吉布斯环状伪影得到很好的抑制，图3(j)为基于SFA的SPEED重建结果，图3(k)为图3(j)与图3(a)相比较的误差映射图。从图3中可以看出，图3(j)的质量优于图3(e)，从误差映射图(f)和(k)的主观评价可以看出，本发明的方法误差比较小。从客观的相对均方差误差测量值(TRE)也可看出，3(j)的TRE值为8.22e-4，明显小于图3(e)的TRE值1.13e-3。可见本发明通过奇异谱分析重建抑制了低分辨率图中的吉布斯环状伪影，提高了定位图像的质量和鬼影的定位精度，通过定位图像来确定鬼影阶数，只需采集两组欠采样数据，提高了SPEED方法的数据采集速度，相比低分辨率图像作为定位图像的SPEED方法，获得了更高的成像质量。

Claims (1)

1.基于奇异谱分析的SPEED快速磁共振成像方法,其特征在于，该方法包括以下步骤：
步骤1：k空间数据采集
在k空间的相位编码方向，即PE方向，每隔N行采集一行数据，共采集两组，分别用S1和S2表示，用kx表示在k空间频率编码FE方向的位置，ky表示在PE方向的位置，则S1(kx,ky)和S2(kx,ky)分别表示S1和S2中的一个数据点的值，用d1,d2表示每组欠采样数据在PE方向上的偏移量，采样方式用N(d1,d2)表示；根据图像大小，在PE方向的k空间中心区域采集16至64行数据，用Sc表示，Sc(kx,ky)表示Sc中的一个数据点；在FE方向的数据都为全采集；
步骤2：填零重建
填零重建包含3个步骤：两组欠采样数据的填零重建、差分变换、k空间中心数据填零重建成低分辨率图像；
步骤2-1：两组欠采样数据的填零重建
对于两组欠采样数据S1和S2，其对应k空间中没有进行数据采集的点用0表示，进行常规的填零傅立叶重建，重建后图像分别用I1和I2表示，用(x,y)表示图像空间中每个像素点的坐标值，则I1(x,y)和I2(x,y)表示I1和I2中的一个像素点的值；由于k空间中每隔N行采集一行数据，因此每组数据对应的填零傅立叶重建图像中有N层重叠的鬼影，每个像素点上最多有N层重叠的鬼影；
步骤2-2：差分变换
对步骤2-1得到的图像I1和I2分别进行差分变换，得到稀疏的边缘鬼影图像E1和E2；稀疏的鬼影图像E1和E2中的每个像素点E1(x,y)和E2(x,y)分别表示为：
E1(x,y)＝I1(x,y)-I1(x,y-1) [1]
E2(x,y)＝I2(x,y)-I2(x,y-1) [2]
步骤2-3：k空间中心数据填零重建成低分辨率图像
采集到的k空间中心部分数据Sc，其对应k空间中没有进行数据采集的点用0表示，进行常规的填零傅立叶重建，形成一个低分辨率的重建图像Ilow，其对应的每一个像素点用Ilow(x,y)表示；
步骤3：奇异谱分析重建
奇异谱分析重建包含8个步骤：差分变换、确立采集矩阵、寻找奇异点、层析法更新差分图像、构建奇异谱函数、计算奇异度、重建出k空间缺失数据、基于离散傅立叶逆变换的图像重建；
步骤3-1：差分变换
对填零重建图像Ilow在y方向上进行差分变换，得到稀疏的差分图像ΔIlow，对应的每个像素点表示为ΔIlow(x,y)：
ΔIlow(x,y)＝Ilow(x,y)-Ilow(x,y-1) [3]
阈值T设为差分图像ΔIlow中幅值最大值的十分之一；
步骤3-2：确立采集矩阵
采集矩阵Mask中的每一点Mask(kx,ky)表示为：
采集矩阵对应的频域点表示为：
i(x,y)＝IDFT(Mask(kx,ky)) [5]
步骤3-3：寻找奇异点
取ΔIlow中幅值最大的点为奇异点j(xj,yj)，加入奇异点集SP中；把i(x,y)矩阵的原点移动至奇异点j(xj,yj)处，更新采集矩阵的频域表示为i(x-xj,y-yj)；
步骤3-4：层析法更新差分图像
更新后的差分图像ΔIlow(x,y)new表示为，
ΔIlow(x,y)new＝ΔIlow(x,y)-α·i(x-xj,y-yj) [6]
α＝ΔIlow(x,y)||Mask||/(M×N) [7]
其中M和N表示x和y方向的像素点数，||Mask||表示采集矩阵中非零点的总数；
如果ΔIlow(x,y)new≥T，则令ΔIlow(x,y)＝ΔIlow(x,y)new，转到步骤3-3继续寻找奇异点；如果ΔIlow(x,y)new＜T，则停止奇异点的寻找，奇异点集SP的奇异点的总个数为Q，继续执行步骤3-5；
步骤3-5：构建奇异谱函数
对奇异点集SP中的奇异点分别构建奇异谱函数Wj(kx,ky)，用公式表示为：
Wj(kx,ky)＝IDFT(δ(y-yj)u(x-xj)) [8]
其中，
IDFT表示离散傅立叶逆变换；
步骤3-6：计算奇异度
利用已采集到的k空间数据Sc通过求解方程[11]来得到各个奇异点对应的奇异度aj，
步骤3-7：重建出k空间缺失数据
用奇异谱函数和奇异度来估算出k空间中没有进行过数据采集处的数据Sr(kx,ky)：
步骤3-8：基于离散傅立叶逆变换的图像重建
对于k空间中进行过数据采集的点用实际采集到的数据Sc(kx,ky)，没有进行过采集的点用步骤3.7估算得到的Sr(kx,ky)，然后利用离散傅立叶逆变换进行图像重建，得到奇异谱重建图像ISFA，其对应的像素点表示为，
步骤4：鬼影定位
鬼影定位包含4个步骤：合成定位图像、差分变换、建立重叠鬼影图和确立鬼影阶数；
步骤4-1：合成定位图像
将步骤2-3得到的低分辨率图像Ilow与步骤3-8得到的重建图像ISFA进行对应像素点之间的相乘操作，合成用于鬼影定位的定位图像Ic，对应的像素点表示为：
Ic(x,y)＝Ilow(x,y)×ISFA(x,y) [14]
步骤4-2：差分变换
对Ic进行差分变换，得到稀疏的边缘鬼影图像Ec；
步骤4-3：建立重叠鬼影图
在相位编码方向对Ec分别进行长度为Ny×n/N的平移，其中Ny表示沿PE方向的数据矩阵的大小，n表示边缘鬼影的阶数，n＝0,1,2,…,N-1；这n个稀疏的边缘鬼影相加后形成一个重叠的鬼影映射图Ec,n；
步骤4-4：确立鬼影阶数
在鬼影映射图Ec,n中，为每个像素点找出两个最强的鬼影，并记录下它们对应的鬼影阶数(n1,n2)；
步骤5：基于双层鬼影模型的SPEED图像重建
基于双层鬼影模型的SPEED图像重建包含4个步骤：双层鬼影模型求解、重叠鬼影的分离、多个分离鬼影图像的配准求和、逆滤波重建；
步骤5-1：双层鬼影模型求解
稀疏边缘鬼影图像E1和E2中，采用双层稀疏边缘鬼影模型来描述，表示为：
公式[15]中为相位因子，Gn1和Gn2分别为每个像素点上需要确定的不同阶的鬼影，n1和n2分别表示不同的鬼影阶数；定义为：
公式[16]中d表示每组欠采样数据在PE方向上的偏移量d1和d2，n为鬼影阶数；
在公式[15]中，由于E1、E2、d和N已知，步骤3-3得到了鬼影阶数(n1,n2)，即公式[15]中的两个方程，仅有两个未知数Gn1和Gn2，因此直接解出公式[15]中的两个重叠的鬼影Gn1和Gn2；
步骤5-2：重叠鬼影的分离
对步骤5-1得到的Gn1和Gn2中的像素点，按不同的鬼影阶数n进行分类，产生N个分离的鬼影映射图Gn,其中n＝0,1,…,N-1；
步骤5-3：多个分离鬼影图像的配准求和
步骤5-2得到的N个鬼影映射图Gn，各自对应的鬼影位置不同，通过移位和对齐来配准；配准后各鬼影图对应的像素点求和后，得到没有重叠鬼影的边缘映射图像E0；
步骤5-4：逆滤波重建
步骤5-3得到的边缘映射图像E0经离散傅立叶变换到k空间，其对应k空间中实际进行数据采集的点的值用实际采集的数据替代，得到k空间数据R0；基于逆滤波公式[17]重建出最终的SPEED图像I0；
公式[17]中ky表示沿PE方向的k空间位置。
CN201711337327.8A 2017-12-14 2017-12-14 基于奇异谱分析的speed快速磁共振成像方法 Pending CN108051766A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711337327.8A CN108051766A (zh) 2017-12-14 2017-12-14 基于奇异谱分析的speed快速磁共振成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711337327.8A CN108051766A (zh) 2017-12-14 2017-12-14 基于奇异谱分析的speed快速磁共振成像方法

Publications (1)

Publication Number Publication Date
CN108051766A true CN108051766A (zh) 2018-05-18

Family

ID=62132897

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711337327.8A Pending CN108051766A (zh) 2017-12-14 2017-12-14 基于奇异谱分析的speed快速磁共振成像方法

Country Status (1)

Country Link
CN (1) CN108051766A (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101051075A (zh) * 2007-04-24 2007-10-10 骆建华 基于复奇异谱分析的磁共振部分k数据图像重建方法
CN103728581A (zh) * 2013-12-20 2014-04-16 杭州电子科技大学 基于离散余弦变换的speed快速磁共振成像方法
CN106526511A (zh) * 2016-10-21 2017-03-22 杭州电子科技大学 基于k空间中心鬼影定位的SPEED磁共振成像方法
CN107290699A (zh) * 2017-06-27 2017-10-24 杭州电子科技大学 基于奇异谱分析的快速动静脉数据同时成像的方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101051075A (zh) * 2007-04-24 2007-10-10 骆建华 基于复奇异谱分析的磁共振部分k数据图像重建方法
CN103728581A (zh) * 2013-12-20 2014-04-16 杭州电子科技大学 基于离散余弦变换的speed快速磁共振成像方法
CN106526511A (zh) * 2016-10-21 2017-03-22 杭州电子科技大学 基于k空间中心鬼影定位的SPEED磁共振成像方法
CN107290699A (zh) * 2017-06-27 2017-10-24 杭州电子科技大学 基于奇异谱分析的快速动静脉数据同时成像的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
JIANHUA LUO等: ""MRI Reconstruction From 2D Truncated k-Space"", 《JOURNAL OF MAGNETIC RESONANCE IMAGING》 *

Similar Documents

Publication Publication Date Title
US11079456B2 (en) Method of reconstructing magnetic resonance image data
CN111081354B (zh) 用于通过深度学习网络对医疗图像进行去噪的系统和方法
Liu et al. Highly undersampled magnetic resonance image reconstruction using two-level Bregman method with dictionary updating
CN106526511B (zh) 基于k空间中心鬼影定位的SPEED磁共振成像方法
US20140091793A1 (en) Method for diffusion magnetic resonance imaging
CN113096208A (zh) 基于双域交替卷积的神经网络磁共振图像的重建方法
CN112526423B (zh) 基于共轭和层间信息的并行磁共振成像算法
CN102590773B (zh) 磁共振成像的方法和系统
CN103584864A (zh) 一种磁共振成像方法和装置
US20190317172A1 (en) Method and apparatus for reconstructing magnetic resonance image data
CN105678822B (zh) 一种基于Split Bregman迭代的三正则磁共振图像重构方法
Falvo et al. A multimodal dense u-net for accelerating multiple sclerosis mri
CN107942271B (zh) 基于迭代的speed快速磁共振成像方法
Yan et al. SMIR: A Transformer-Based Model for MRI super-resolution reconstruction
Cha et al. Geometric approaches to increase the expressivity of deep neural networks for MR reconstruction
Lv et al. Parallel imaging with a combination of sensitivity encoding and generative adversarial networks
Gan et al. Deep image reconstruction using unregistered measurements without groundtruth
Xiao et al. SR-Net: a sequence offset fusion net and refine net for undersampled multislice MR image reconstruction
Virtue et al. Learning contrast synthesis from MR fingerprinting
CN111161370B (zh) 一种基于ai的人体多核dwi联合重建方法
CN103728581A (zh) 基于离散余弦变换的speed快速磁共振成像方法
CN108051766A (zh) 基于奇异谱分析的speed快速磁共振成像方法
CN112557981B (zh) 一种并行磁共振成像的改进算法
Cao et al. CS-GAN for high-quality diffusion tensor imaging
Cheng et al. deep MR parametric mapping with unsupervised multi-tasking framework

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20180518