CN101810480A - 基于缺失数据重构的磁共振图像截断伪影消除方法 - Google Patents

Abstract

一种图像处理技术领域的基于缺失数据重构的磁共振图像截断伪影消除方法，包括以下步骤：对含截断伪影的磁共振图像进行检测，得到低频K数据；对低频K数据进行频率编码方向的傅里叶反变换，得到相位编码方向上一维傅里叶变换的截断频谱G(k_x，y)；对G(k_x，y)的每一行数据进行该行图像的奇异点提取和奇异值提取，并进行重构得到高频部分K数据；对高频部分K数据F_y(k)进行傅里叶反变换，得到去伪影磁共振图像的行信息f_y(x)，将f_y(x)按行存放组合成为去伪影后的磁共振图像g(x，y)。本发明克服补零法成像中所存在的伪影问题，确保图像的高信噪比，有效降低图像误差，精确显示原磁共振图像，为医学核磁共振检测提供了高质量的可靠图像信息，有助于医学成像检测技术的发展和普及。

Description

基于缺失数据重构的磁共振图像截断伪影消除方法

技术领域

本发明涉及的是一种图像处理技术领域的方法，具体是一种基于缺失数据重构的磁共振图像截断伪影消除方法。

背景技术

随着现代医学技术的不断发展，磁共振成像(Magnetic Resonance Imaging，MRI)技术已经成为医学成像检测领域中不可或缺的手段。磁共振信号空间(原始数据空间)称为K空间，即为傅里叶变换空间，K空间采样信号经过傅里叶反变换后再取模，即得磁共振图像。磁共振成像临床应用中，常常用仅扫描K空间的低频数据进行成像，以节省缩短扫描时间。对于仅用低频K空间数据成像的临床方法是补零法成像，即用将欠缺的高频K空间用零填补，然后傅里叶反变换成像。这类图像的缺点是图像含截断伪影，也称吉布斯伪影(Gibbs artifacts)，它容易出现于图像像素值变化大的地方，这是因为这些地方原本高频成份丰富，低频部分K数据采集的方法把高频分量截掉了，因而产生较大的截断伪影。

截断频谱在数学上等价于在频域乘上一个矩形窗。截断频谱补零法的图像就等价于完整谱图像与一个矩形窗的傅里叶反变换函数的卷积结果。因而消除截断伪迹方法可以是对矩形进行修正为汉明窗(Hamming windows，汉明窗)，以有效消除陡峭频谱边缘引起的截断伪影。但是这种方法进一步损伤频谱数据，引入更大的误差，使图像的边缘特性进一步受到损伤，使图像变得更加模糊。到目前为止，还没有一种保证图像不失真前提下，并能适合临床应用的消除截断伪影方法，在临床只使用补零法得到的有伪影图像。如何从有伪影医学图像中，消除截断伪影是一项有实际应用价值的课题。

图像截断伪影是因频谱数据采集过程中为求扫描速度，只扫低频相位的部分K数据，没有采集高频相位K数据所致。因而从所给图像中找回图像的高频分量，从而恢复完整频谱的图像，消除截断伪影是最佳的除截断伪影方法。目前部分K数据成像方法主要有：补零法、迭代正则化方法、稀疏表达方法等等。补零法是将缺失的K数据用零补充，然后直接用傅里叶反变换得到图像。优点是简单可靠，缺点是图像有伪影，迭代正则化法是法在一定约束条件下，优化目标函数获得最优解，如最大熵方法。稀疏变换法也是通过优化方法获图像稀疏表达，然后再由稀疏表达函数来重构图像。后两种方法的成像质量都依赖于目标函数和优化策略的选择，与实际图像有相当误差。

经对现有文献检索发现，未发现基于缺失数据重构高频K数据的磁共振图像截断伪影消除方法的公开报道。

发明内容

本发明的目的是克服上述现有技术中的缺点，提供一种基于缺失数据重构的磁共振图像截断伪影消除方法。本发明从含截断伪影的图像中，提取低频K数据，然后从这些低频数据用一维奇异函数分析方法重构高频K数据，最后重构出完整频谱数据的图像，得到既保精度又能消除截断伪影的高精度图像。

本发明是通过以下技术方案实现的，包括以下步骤：

第一步，对含截断伪影的磁共振图像进行检测，得到低频K数据。

所述的检测，包括以下步骤：

(1)对含截断伪影的磁共振图像进行傅里叶变换得到二维K数据，即：

其中：

是尺寸为N×N含截断伪影的磁共振图像，

为二维K数据，k_x为磁共振扫描的相位编码，k_y为频率编码，为二维傅里叶变换算子；

(2)

在k_x方向大小突然变化的点就是截断频率c，0＜c＜N/2，从而得到截断频率的矩形函数R_c(k_x，k_y)，即： 

(3)将二维K数据

和截断频率的矩阵函数R_c(k_x，k_y)相乘，得到低频K数据。

第二步，对低频K数据进行频率编码k_y方向的傅里叶反变换，得到相位编码k_x方向上一维傅里叶变换的截断频谱G(k_x，y)。

所述的k_y方向的傅里叶反变换，是指：

$G(k_x,y)=\underset{k_y=-N/2}{\overset{N/2-1}{\mathrm{\Σ}}}\tilde{G}(k_x,k_y)R_c(k_x,k_y)e^{-j\frac{2\mathrm{\πy}{k}_y}{N}},y=\mathrm{1,2},...,N,k_x=-N/2,-N/2+1,...,N/2-1$

其中：G(k_x，y)是k_x方向上一维傅里叶变换的截断频谱，

是二维K数据，R_c(k_x，k_y)是截断频率的矩形函数，k_x为磁共振扫描的相位编码，k_y为频率编码，N是含截断伪影的磁共振图像的尺寸。

第三步，对G(k_x，y)的每一行数据进行该行图像的奇异点提取和奇异值提取，并进行重构得到高频部分K数据。

所述的奇异点提取，具体是：对G(k_x，y)的每一行数据进行补零法重构，得到该行的近似图像信号，并利用层析反卷积法检测得到每行近似图像信号中的奇异点集。

所述的补零法重构，是：

其中：

是第y行的近似图像信号，x表示该行的横坐标信息，k_x为磁共振扫描的相位编码。

所述的层析反卷积法，包括以下步骤：

1)估计近似图像信号

的差分信号，具体是：

$\widetilde{d_y^1}\left(x\right)={\tilde{f}}_y\left(x\right)-{\tilde{f}}_y(x-1),$

其中：

是第y行、横坐标为x的近似图像信号的差分信号，

是第y行、横坐标为x的近似图像信号，

是第y行、横坐标为(x-1)的近似图像信号，2≤x≤N。

2)找到

处的横坐标值x₁，将x₁加入到该行近似图像信号中的奇异点集中，根据

将

更新为

其中：

N是含截断伪影的磁共振图像的尺寸，c是截断频率。

3)当

时，T是系统预设的噪声阈值，则结束处理；

否则，按照2)的方法对进行处理，直至找到

处的横坐标值x_n，将x_n加入到该行近似图像信号中的奇异点集中，n是奇异点集中奇异点的序列号，1≤n≤N，当

时，则结束处理。

所述的奇异值估计是用伪逆矩阵法求解奇异值方程，得到最小误差解，获得每个奇异点对应的奇异值，其中的奇异值方程具体是：

$G(k_x,y)(1-e^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x})=\underset{l=\mathrm{1,2},...q}{\mathrm{\&Sigma;}}{a}_l{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_l},\left|k_x\right|<c$

其中：q是该行数据的奇异点总个数，x_l是第l个奇异点，a_l是奇异点x_l对应的奇异值，N是含截断伪影的磁共振图像的尺寸，c是截断频率，k_x为磁共振扫描的相位编码。

所述的重构，是指：

$F_y\left(k_x\right)=(a_1{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_1}+a_2{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_2}+...+a_q{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_q})/(1-e^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x}),\left|k_x\right|>c$

其中：F_y(k_x)是高频部分K数据，x_l是第l个奇异点，a_l是奇异点x_l对应的奇异值，N是含截断伪影的磁共振图像的尺寸，c是截断频率，q是该行数据的奇异点总个数，1≤l≤q。

第四步，对高频部分K数据F_y(k)进行傅里叶反变换，得到去伪影磁共振图像的行信息f_y(x)，将f_y(x)按行存放组合成为去伪影后的磁共振图像g(x，y)。

与现有技术相比，本发明的有益效果是：采用了该发明的基于缺失数据重构的磁共振图像截断伪影消除方法，由于对缺失部分的频谱数据进行了重构，因而能够克服补零法成像中所存在的伪影问题，确保了图像的高信噪比，有效降低图像误差，精确显示原磁共振图像，为医学核磁共振检测提供了高质量的可靠图像信息，有助于医学成像检测技术的进一步发展和应用普及。

附图说明

图1是实施例方法和现有技术方法的噪声级别-误差峰值信噪比的比较示意图。

图2是实施例去除伪影前后的比较示意图；

其中：(a)是参照图像，(b)是去伪影前的图像，(c)是用汉明窗法去伪影后得到的图像，(d)是用本实施例方法去伪影后得到的图像。

图3是图2对应的频谱图；

其中：(a)是图2(a)对应的频谱图，(b)是图2(b)对应的频谱图，(c)是图2(c)对应的频谱图，(d)是图2(d)对应的频谱图。

图4是实施例方法和现有技术方法的截断频率-误差峰值信噪比的比较示意图。

图5是实施例实际图像去除伪影前后的比较示意图；

其中：(a)是去伪影前的图像，(b)是采用汉明窗法去伪影得到的图像，(c)是采用本实施例方法去伪影得到的图像。

具体实施方式

以下结合附图对本发明的方法进一步描述：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例

本实施例包括以下步骤：

第一步，对含截断伪影的磁共振图像进行检测，得到低频K数据。

所述的检测，包括以下步骤：

(1)对含截断伪影的磁共振图像进行傅里叶变换得到二维K数据，即：

其中：

是尺寸为N×N含截断伪影的磁共振图像，为二维K数据，k_x为磁共振扫描的相位编码，k_y为频率编码，为二维傅里叶变换算子；

(2)在k_x方向大小突然变化的点就是截断频率c，0＜c＜N/2，从而得到截断频率的矩形函数R_c(k_x，k_y)，即：

(3)将二维K数据

和截断频率的矩阵函数R_c(k_x，k_y)相乘，得到低频K数据。

第二步，对低频K数据进行频率编码k_y方向的傅里叶反变换，得到相位编码k_x方向上一维傅里叶变换的截断频谱G(k_x，y)。

所述的k_y方向的傅里叶反变换，是指：

$G(k_x,y)=\underset{k_y=-N/2}{\overset{N/2-1}{\mathrm{\&Sigma;}}}\tilde{G}(k_x,k_y)R_c(k_x,k_y)e^{-j\frac{2\mathrm{\&pi;y}{k}_y}{N}},y=\mathrm{1,2},...,N,k_x=-N/2,-N/2+1,...,N/2-1$

其中：G(k_x，y)是k_x方向上一维傅里叶变换的截断频谱，

是二维K数据，R_c(k_x，k_y)是截断频率的矩形函数，k_x为磁共振扫描的相位编码，k_y为频率编码，N是含截断伪影的磁共振图像的尺寸。

第三步，对G(k_x，y)的每一行数据进行该行图像的奇异点提取和奇异值提取，并进行重构得到高频部分K数据。

所述的奇异点提取，具体是：对G(k_x，y)的每一行数据进行补零法重构，得到该行的近似图像信号，并利用层析反卷积法检测得到每行近似图像信号中的奇异点集。

所述的补零法重构，是：

其中：

是第y行的近似图像信号，x表示该行的横坐标信息，k_x为磁共振扫描的相位编码。

所述的层析反卷积法，包括以下步骤：

1)估计近似图像信号

的差分信号，具体是：

$\widetilde{d_y^1}\left(x\right)={\tilde{f}}_y\left(x\right)-{\tilde{f}}_y(x-1),$

其中：

是第y行、横坐标为x的近似图像信号的差分信号，

是第y行、横坐标为x的近似图像信号，

是第y行、横坐标为(x-1)的近似图像信号，2≤x≤N。

2)找到

处的横坐标值x₁，将x₁加入到该行近似图像信号中的奇异点集中，根据

将更新为

其中：

N是含截断伪影的磁共振图像的尺寸，c是截断频率。

3)当

时，T是系统预设的噪声阈值，则结束处理；

否则，按照2)的方法对

进行处理，直至找到

处的横坐标值x_n，将x_n加入到该行近似图像信号中的奇异点集中，n是奇异点集中奇异点的序列号，1≤n≤N，当

时，则结束处理。

所述的奇异值估计是用伪逆矩阵法求解奇异值方程，得到最小误差解，获得每个奇异点对应的奇异值，其中的奇异值方程具体是：

$F_y\left(k_x\right)(1-e^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x})=\underset{l=\mathrm{1,2},...q}{\mathrm{\&Sigma;}}{a}_l{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_l},\left|k_x\right|<c$

其中：q是该行数据的奇异点总个数，x_l是第l个奇异点，a_l是奇异点x_l对应的奇异值，N是含截断伪影的磁共振图像的尺寸，c是截断频率，k_x为磁共振扫描的相位编码，当|k_x|＜c，F_y(k_x)＝G(k_x，y)。因此写成矩阵形式为：WA＝D，

其中：

A＝[a₁，a₂，...，a_q]^T， $D={[F_y(-c)(1-e^{i\frac{2\mathrm{\&pi;}}{N}c}),F_y(-c+1)(1-e^{i\frac{2\mathrm{\&pi;}}{N}(c-1)}),...,F_y\left(c\right)(1-e^{-i\frac{2\mathrm{\&pi;}}{N}c})]}^T$

利用伪逆矩阵法求解上述方程组，可以得到方程的一个伪逆解：

A＝W⁺D

其中：W⁺＝(W^*W)^-1W^*为W的伪逆矩阵，W^*为W的共轭转置，A是该行数据的奇异点集。

所述的重构，是指：

$F_y\left(k_x\right)=(a_1{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_1}+a_2{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_2}+...+a_q{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_q})/(1-e^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x}),\left|k_x\right|>c$

其中：F_y(k_x)是高频部分K数据，x_l是第l个奇异点，a_l是奇异点x_l对应的奇异值，N是含截断伪影的磁共振图像的尺寸，c是截断频率，q是该行数据的奇异点总个数，1≤l≤q。

第四步，对高频部分K数据F_y(k)进行傅里叶反变换，得到去伪影磁共振图像的行信息f_y(x)，即：

将f_y(x)按行存放组合成为去伪影后的磁共振图像g(x，y)，即：

将http://www.bic.mni.mcgill.ca/brainweb/上的噪声分别是0-9等级的PD卷中的第100片取出作为10幅参照图像，对每一幅参照图像进行用截断频率c＝N/4(N为256)的低频部分频谱用补零法成像，获得9幅为去伪影前的图像。分别用汉明窗法及本实施例方法消除伪影，去伪影前后的图像误差峰值信噪比如图1所示，其中：横坐标为噪声级别，纵坐标表示误差峰值信噪比。从图1可以看出，用汉明窗法去伪影的图像在各级噪声下都不及去伪影前的图像，而本实施例方法去伪影的图像误差峰值信噪比比去伪影前有明显提高，对噪声小的图像本实施例方法的误差峰值信噪比提升更为显著。本实施例方法去伪影后图像误差峰值信噪比大为提高，而汉明窗法的去伪影后图像的峰值信噪比不升反而下降，其原因是汉明窗法为了消除伪影进一步加大频谱误差，这是牺牲图像的分辨率换取无截断伪影的效果。

图2(a)是一参照图像，图2(b)是去伪影前的图像，图2(c)是用汉明窗法去伪影后得到的图像，图2(d)是用本实施例方法去伪影后得到的图像。通过比较图2(b)和图2(c)可以发现汉明窗法可以去掉伪影，但图像变得有更加模糊。而图2(d)不但伪影可以去掉，而且图像也变得更加清晰的，比去伪影前更加接近图2(a)。图2(a)所对应的频谱图如图3(a)所示，图2(b)所对应的频谱图如图3(b)所示，图2(c)所对应的频谱图如图3(c)所示，图2(d)所对应的频谱图如图3(d)所示。比较图3(c)和图3(b)就可以很明显看到汉明窗法去伪影后图像高频频谱不如伪影前图像的频谱，汉明窗法比去伪影前的图像丢失更多的频谱分量，而本实施例方法图像的频谱近乎和参照图像的一样。

将http://www.bic.mni.mcgill.ca/brainweb/上噪声等级为4的T1加权的第100横切面图像作为参考图像，并使用截断频率c从40到100(对于256X256的图像最大截止频率为128)范围以5为间隔构造含截断伪影图像，分别用汉明窗法和本实施例方法进行去伪影，得到的去伪影前后图像的误差峰值信噪比如图4所示，其中：横坐标表示截断频率变化，纵坐标表示误差峰值信噪比。由图4可知：在各种截断频率下，本实施例方法去伪影后图像的误差峰值信噪比较去伪影前有显著提高，而汉明窗法去伪影后图像的误差峰值信噪比比去伪影前下降。随着截断频率上升，去伪前图像的误差峰值信噪比提高，本实施例方法所能提高的误差峰值信噪比变小，这是因为截止频率越高，去伪影前图像的伪影变得越小。故本实施例方法在各种截止频率下都有很好的去伪影效果。

图5(a)所示是c＝64、像素为256x256的实际医学头部图像去伪影前的图像，图5(b)是采用汉明窗法去伪影后得到的图像，图5(c)是采用本实施例方法去伪影后得到的图像。比较这三幅图像可知：本实施例方法在消除伪影的基础上，图像清晰度有所提高，而汉明窗法则使去伪影后图像变得模糊。

Claims (8)

1.一种基于缺失数据重构的磁共振图像截断伪影消除方法，其特征在于，包括以下步骤：

第一步，对含截断伪影的磁共振图像进行检测，得到低频K数据；

第二步，对低频K数据进行频率编码k_y方向的傅里叶反变换，得到相位编码k_x方向上一维傅里叶变换的截断频谱G(x_k，y)；

第三步，对G(k_x，y)的每一行数据进行该行图像的奇异点提取和奇异值提取，并进行重构得到高频部分K数据F_y(k)；

第四步，对高频部分K数据F_y(k)进行傅里叶反变换，得到去伪影磁共振图像的行信息f_y(x)，将f_y(x)按行存放组合成为去伪影后的磁共振图像g(x，y)。

2.根据权利要求1所述的基于缺失数据重构的磁共振图像截断伪影消除方法，其特征是，第一步中所述的检测，包括以下步骤：

(1)对含截断伪影的磁共振图像进行傅里叶变换得到二维K数据，即：

k_x，k_y＝-N/2，-N/2+1，...，N/2-1，

其中：

是尺寸为N×N含截断伪影的磁共振图像，为二维K数据，k_x为磁共振扫描的相位编码，k_y为频率编码，

为二维傅里叶变换算子；

(2)

在k_x方向大小突然变化的点就是截断频率c，0＜c＜N/2，从而得到截断频率的矩形函数R_c(k_x，k_y)，即：

(3)将二维K数据和截断频率的矩阵函数R_c(k_x，k_y)相乘，得到低频K数据。

3.根据权利要求1所述的基于缺失数据重构的磁共振图像截断伪影消除方法，其特征是，第二步中所述的k_y方向的傅里叶反变换，是指：

y＝1，2，...，N，k_x＝-N/2，-N/2+1，...，N/2-1，

其中：G(k_x，y)是k_x方向上一维傅里叶变换的截断频谱，

是二维K数据，R_c(k_x，k_y)是截断频率的矩形函数，k_x为磁共振扫描的相位编码，k_y为频率编码，N是含截断伪影的磁共振图像的尺寸。

4.根据权利要求1所述的基于缺失数据重构的磁共振图像截断伪影消除方法，其特征是，第三步中所述的奇异点提取，具体是：对G(k_x，y)的每一行数据进行补零法重构，得到该行的近似图像信号，并利用层析反卷积法检测得到每行近似图像信号中的奇异点集。

5.根据权利要求4所述的基于缺失数据重构的磁共振图像截断伪影消除方法，其特征是，所述的补零法重构，是：

其中：

是第y行的近似图像信号，x表示该行的横坐标信息，k_x为磁共振扫描的相位编码。

6.根据权利要求4所述的基于缺失数据重构的磁共振图像截断伪影消除方法，其特征是，所述的层析反卷积法，包括以下步骤：

1)估计近似图像信号

的差分信号，具体是：

${\tilde{d}}_y^1\left(x\right)={\tilde{f}}_y\left(x\right)-{\tilde{f}}_y(x-1),$

其中：

是第y行、横坐标为x的近似图像信号的差分信号，

是第y行、横坐标为x的近似图像信号，

是第y行、横坐标为(x-1)的近似图像信号，2≤x≤N；

2)找到

处的横坐标值x₁，将x₁加入到该行近似图像信号中的奇异点集中，根据将

更新为

其中：

N是含截断伪影的磁共振图像的尺寸，c是截断频率；

3)当

时，T是系统预设的噪声阈值，则结束处理；

否则，按照2)的方法对

进行处理，直至找到处的横坐标值x_n，将x_n加入到该行近似图像信号中的奇异点集中，n是奇异点集中奇异点的序列号，1≤n≤N，当

时，则结束处理。

7.根据权利要求1所述的基于缺失数据重构的磁共振图像截断伪影消除方法，其特征是，第三步中所述的奇异值估计是用伪逆矩阵法求解奇异值方程，得到最小误差解，获得每个奇异点对应的奇异值，其中的奇异值方程具体是：

$G(k_x,y)(1-e^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x})=\underset{l=\mathrm{1,2},...q}{\mathrm{\&Sigma;}}{a}_l{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_l},\left|k_x\right|<c$

其中：q是该行数据的奇异点总个数，x_l是第l个奇异点，a_l是奇异点x_l对应的奇异值，N是含截断伪影的磁共振图像的尺寸，c是截断频率，k_x为磁共振扫描的相位编码。

8.根据权利要求1所述的基于缺失数据重构的磁共振图像截断伪影消除方法，其特征是，第三步中所述的重构，是指：

$F_y\left(k_x\right)=(a_1{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_1}+a_2{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_2}+...+a_q{e}^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x{x}_q})/(1-e^{-i\frac{2\mathrm{\&pi;}}{N}{k}_x}),\left|k_x\right|>c$

其中：F_y(k_x)是高频部分K数据，x_l是第l个奇异点，a_l是奇异点x_l对应的奇异值，N是含截断伪影的磁共振图像的尺寸，c是截断频率，q是该行数据的奇异点总个数，1≤l≤q。

Images