CN102283651B - 一种核磁共振成像梯度场校正方法 - Google Patents

Abstract

本发明公开了一种核磁共振成像梯度场校正方法，只需要厂家提供的梯度场谐波参数，即可对由于梯度非线性度引起的图像畸变，包括像素位置和幅度值都进行校正，适用性强，不增加额外成本。

Description

一种核磁共振成像梯度场校正方法

技术领域

本发明属于核磁共振成像领域，具体涉及一种核磁共振成像梯度场校正方法。

背景技术

在核磁共振成像系统中，梯度线圈产生的三个方向的梯度场对于成像空间起到线性编码的作用。理想的梯度线圈产生的梯度场是完全线性的，但是考虑到梯度线圈的其它性能参数，包括梯度效率，上升时间和梯度切换率，梯度线圈功率损耗等，在梯度线圈的设计中往往需要对各个参数性能进行综合的优化，最终形成一个折中方案，其中，梯度的线性度往往会做出一定的牺牲，尤其是目前为了设计更加开放的新型磁体，线圈设计要求长度短，内径大，对梯度线性度的牺牲更大一些。除开在设计过程中预计内的梯度非线性度，在实际线圈加工生产过程的各个环节中还会引入一系列的公差，以及梯度线圈在强大的主磁场中通电后受到巨大的磁力矩造成的扭曲变形，都会影响最终的梯度场线性度。因此，实际核磁共振成像系统中三维梯度线圈产生的并不是理想的线性梯度场，而具有一定的非线性度的非理想梯度场。

图像的畸变是指所得到的图像信息与实际真实图像的不一致，主要包括两个方面：像素位置的空间平移和像素幅度值的不准确。在核磁共振成像中，引起图像畸变的原因有很多，包括主磁场的不均匀度，化学位移，不同介质的磁化率效应，梯度线圈涡流效应，以及梯度场的非线性。由于核磁共振成像的图像重建算法是基于梯度场是理想线性的假设，非线性梯度会导致图像的几何畸变，这种畸变具有以下几个特点：图像的畸变同时表现在像素位置的平移和像素幅度值的不准确；空间相同位置的图像畸变仅与梯度线圈产生的场的非线性度相关，与具体的成像序列无关，即对于同一台核磁共振成像系统，任意方法序列成出的图像，其由于梯度非线性度引起的空间几何畸变都是一样的；一般在成像区域的边缘位置梯度非线性度更大，所引起的图像畸变也更加显著。

对于梯度非线性引起的图像畸变的校正一般分为两种方法：根据梯度场分布计算得到畸变信息对图像进行校正；或者用特制的水模测量图像畸变从而计算得到校正数据。Andrew Janke等人提出用球谐函数的方法表示梯度场非线性度(Andrew Janke et al.Use of Spherical Harmonic DeconvolutionMethods to Compensate for Nonlinear Gradient Effects on MRI Images.Magnetic Resonance in Medicine，2004，52：115-122)，在此基础上得到图像畸变的校正方法，但是该方法中像素位置的校正需要进行迭代计算，而且只有对位置的校正，没有对像素幅度值的校正。Simon Doran等人提出一种特制的水模对图像畸变进行标定(Simon Doran et al.A complete distortioncorrection for MR images：I.Gradient warp correction.Phys.Med.Biol.，2005，50：1343-1361)，根据测到的图像畸变数据得到像素位置和幅度值的校正公式，但是该方法自制的水模增加了成本，而且需要用CT扫描得到参考图像来校正，还存在图像配准的问题，另外，水模的大小限制了该方法运用的范围，在水模以外区域的图像畸变不能得到校正。

发明内容

本发明目的是：提供一种用于核磁共振成像梯度场几何畸变的校正方法，只需要厂家提供的梯度场谐波参数，即可对由于梯度非线性度引起的图像畸变，包括像素位置和幅度值都进行校正，适用性强，不增加额外成本。

本发明的技术方案是：一种核磁共振成像梯度场校正方法，其特征在于，包括以下步骤：

(1)由核磁共振成像得到未校正的带有畸变的图像数据；

(2)从系统读出成像参数：包括频率编码、相位编码和选层，三个方向的视野以及空间分辨率，成像采集的旋转矩阵，成像中心偏离磁体中心的位置，计算得到每个像素点对应于以磁体中心为原点的绝对坐标系的坐标值；

(3)根据厂家给出的梯度线圈的谐波场分解系数，计算得到梯度场的空间非线性度函数；

(4)根据成像参数和梯度场的空间非线性度函数计算出梯度畸变校正公式；

(5)将梯度非线性函数带入梯度畸变校正公式，对畸变图像进行校正。

进一步的，对于用二维序列采集的图像，可在每个层面内运用二维畸变校正公式对图像进行校正；对于用三维序列采集的图像，需用三维畸变校正公式对三维图像进行校正。

进一步的，所述二维畸变校正公式为：

$\mathrm{\ρ}(x,y)=\frac{{\∂f}_x(x,y)}{\∂x}\frac{{\∂f}_y(x,y)}{\∂y}s(f_x(x,y),f_y(x,y)).$

进一步的，所述三维畸变校正公式为：

$\mathrm{\ρ}(x,y,z)=\frac{{\∂f}_x(x,y,z)}{\∂x}\frac{{\∂f}_y(x,y,z)}{\∂y}\frac{{\∂f}_z(x,y,z)}{\∂z}s(f_x(x,y,z),f_y(x,y,z),f_z(x,y,z)).$

本发明的优点是：

本发明可对核磁共振成像中由于梯度非线性度引起的图像畸变进行校正，其中对于像素位置和幅度值都作出修正，对于同一台磁体梯度线圈可用同样的校正参数，任意成像序列和参数扫描获得的图像都可以用本方法校正，适用性强，校正公式只需要梯度线圈生产厂家提供的梯度场谐波分解参数即可计算，无须水模等额外工具和其它测量工作，不增加成本。

附图说明

下面结合附图及实施例对本发明作进一步描述：

图1为本发明的核磁共振成像梯度场校正方法流程图。

具体实施方式

实施例：一种核磁共振成像梯度场校正方法，如图1所示，对于任意成像扫描序列和参数获得的带有由于梯度非线性度引起的几何畸变的图像，用扫描得到的图像参数结合梯度非线性函数计算校正公式，对像素的位置和像素幅度值都进行修正，得到校正后的图像，其具体流程如下：

1)核磁共振成像中用任意扫描序列以及参数得到k空间数据，直接用反傅里叶变换图像重建，即得到未校正的带有畸变的图像数据。图像数据是一个N_r×N_p×N_s规模的三维复数数据D，三个维度的个数分别为频率编码、相位编码和选层方向的采样数，一般情况只考虑其模像M的校正，即取M＝|D|，M中的元素都为非负实数。

2)核磁共振成像系统进行序列扫描之后，不仅可以读出图像数据，还可以读出成像的参数：包括频率编码、相位编码和选层，三个方向的FOV：FOV_r，FOV_p，FOV_s，以及空间分辨率： $\mathrm{dr}=\frac{{\mathrm{FOV}}_r}{N_r},$ $\mathrm{dp}=\frac{{\mathrm{FOV}}_p}{N_p},$ $\mathrm{ds}=\frac{{\mathrm{FOV}}_s}{N_s},$ 成像采集的旋转矩阵RotMat，成像中心偏离磁体中心的位置Xoffset，Yoffset，Zoffset，可以由每个像素点在成像相对坐标位置P′(x′，y′，z′)计算其对应于以磁体中心为原点的绝对坐标系的坐标值P(x，y，z)：

$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\mathrm{RotMat}\×\left(\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right)+\left(\begin{array}{c}\mathrm{Xoffset}\\ \mathrm{Yoffset}\\ \mathrm{Zoffset}\end{array}\right).$

3)在核磁共振磁体成像区域的磁场可以用谐波分解为如下形式，其中谐波系数a_nml以为μT单位：

$B_z(r,\mathrm{\θ},\mathrm{\φ})=\underset{l=\mathrm{1,2}}{\underset{m=0,n}{\underset{n=0,n_{\max }}{\mathrm{\Σ}}}}{P}_{\mathrm{nm}}\left(\cos \mathrm{\θ}\right)\×{\left(\frac{r}{r_{\mathrm{ref}}}\right)}^n\×[a_{\mathrm{nml}}\cos \left(\mathrm{m\φ}\right)+a_{\mathrm{nm}2}\sin \left(\mathrm{m\φ}\right)],$

其中P_nm(cosθ)表示连带勒让德函数，定义如下：

$P_n^m\left(x\right)=\frac{{(1-x^2)}^{m/2}}{2^{n}n!}\×\frac{d^{m+n}}{{\mathrm{dx}}^{m+n}}{(x^2-1)}^n.$

梯度线圈的生产厂家提供在一定参考半径r_ref内梯度场的谐波分解系数a_nml，根据上式可以计算得到梯度场的空间非线性度函数，以x梯度为例，非线性度函数f_x(x)定义为实际梯度场G_r，x(x)与理想梯度强度G_x之比：

$f_x\left(x\right)=\frac{G_{r,x}\left(x\right)}{G_x}$

$=\frac{\underset{m=1,n}{\underset{n=1,n_{\max }}{\mathrm{\Σ}}}{P}_{\mathrm{nm}}\left(\cos \mathrm{\θ}\right)\×{\left(\frac{r}{r_{\mathrm{ref}}}\right)}^n\×a_{\mathrm{nml}}\cos \left(\mathrm{m\φ}\right)}{P_{\mathrm{1,1}}\left(\cos \mathrm{\θ}\right)\×\frac{r}{r_{\mathrm{ref}}}\×a_{\mathrm{1,1,1}}\cos \mathrm{\φ}}x,$

其中 $\left\{\begin{array}{c}r=\sqrt{x^2+y^2+z^2}\\ \mathrm{\θ}={\cos }^{-1}\left(\frac{z}{r}\right)\\ \mathrm{\φ}={\tan }^{-1}\left(\frac{y}{x}\right)\end{array}\right..$

4)由梯度非线性函数可以计算几何畸变校正公式，其推导过程如下：

先考虑一维情况，已知x梯度非线性函数为实际信号分布即真实的图像为ρ(x)，因此可以写出时域信号：

$s\left(t\right)=\∫\mathrm{dx\ρ}\left(x\right)e^{-\mathrm{i\γ}\underset{0}{\overset{t}{\∫}}{G}_x\left(t^{\′}\right)f_x\left(x\right)d{t}^{\′}},$

由k空间的定义可以将时域信号转化为k空间信号：

$s\left(k\right)=\∫\mathrm{dx\ρ}\left(x\right)e^{-2\mathrm{\πik}{f}_x\left(x\right)},$

将此k空间数据直接进行反傅里叶变换图像重建，得到的是带有几何畸变的未校正图像：

$s\left(x\right)=F^{-1}\left\{s\left(k\right)\right\}$

$=\∫\mathrm{dk}{e}^{2\mathrm{\πikx}}\∫{\mathrm{dx}}^{\′}\mathrm{\ρ}\left(x^{\′}\right)e^{-2\mathrm{\πik}{f}_x\left(x^{\′}\right)}$

$=\∫{\mathrm{dx}}^{\′}\mathrm{\ρ}\left(x^{\′}\right)\∫\mathrm{dk}{e}^{2\mathrm{\πik}[x-f_x\left(x^{\′}\right)]},$

$=\∫{\mathrm{dx}}^{\′}\mathrm{\ρ}\left(x^{\′}\right)\mathrm{\δ}[f_x\left(x^{\′}\right)-x]$

在此做一次变量代换，

$\left\{\begin{array}{c}X=f_x\left(x^{\′}\right)\\ x^{\′}{=f}_x^{-1}\left(X\right)\end{array}\right.,$

得到：

$s\left(x\right)=\∫d{f}_x^{-1}\left(X\right)\mathrm{\ρ}\left(f_x^{-1}\left(X\right)\right)\mathrm{\δ}(X-x)$

$=\∫\mathrm{dX}\frac{{\mathrm{df}}_x^{-1}\left(X\right)}{\mathrm{dX}}\mathrm{\ρ}\left(f_x^{-1}\left(X\right)\right)\mathrm{\δ}(X-x),$

$=\frac{{\mathrm{df}}_x^{-1}\left(x\right)}{\mathrm{dx}}\mathrm{\ρ}\left(f_x^{-1}\left(x\right)\right)$

也即：

$\mathrm{\ρ}\left(f_x^{-1}\left(x\right)\right)=\frac{\mathrm{dx}}{d{f}_x^{-1}\left(x\right)}s\left(x\right),$

再做一次变量代换，

$\left\{\begin{array}{c}x=f_x\left(X^{\′}\right)\\ X^{\′}=f_x^{-1}\left(x\right)\end{array}\right.,$

可以得到：

$\mathrm{\ρ}\left(X^{\′}\right)=\frac{{\mathrm{df}}_x\left(X^{\′}\right)}{d{X}^{\′}}s\left(f_x\left(X^{\′}\right)\right),$

即：

$\mathrm{\ρ}\left(x\right)=\frac{{\mathrm{df}}_x\left(x\right)}{\mathrm{dx}}s\left(f_x\left(x\right)\right).$

此式表示真实图像ρ(x)的像素位置可以用畸变的图像s(x)在梯度非线性函数f_x(x)的插值变换下得到s(f_x(x))，并且真实图像的信号幅度值需要用梯度非线性函数的一阶微分

来做修正。

由一维校正公式可以推广到二维校正：

$\mathrm{\ρ}(x,y)=\frac{{\∂f}_x(x,y)}{\∂x}\frac{{\∂f}_y(x,y)}{\∂y}s(f_x(x,y),f_y(x,y)),$

以及三维校正：

$\mathrm{\ρ}(x,y,z)=\frac{{\∂f}_x(x,y,z)}{\∂x}\frac{{\∂f}_y(x,y,z)}{\∂y}\frac{{\∂f}_z(x,y,z)}{\∂z}s(f_x(x,y,z),f_y(x,y,z),f_z(x,y,z)).$

对于三维图像数据，用厂家提供的谐波分解系数可以分别求出三个梯度线圈的非线性度函数f_x(x，y，z)，f_y(x，y，z)，f_z(x，y，z)。

将梯度非线性函数带入梯度畸变校正公式，根据此公式在图像域上对于畸变图像进行校正，校正具体分为两部分内容：

a)根据校正公式的位移部分，采用线性插值的方法将像素位置平移，使得像素点由畸变图像中的(f_x(x，y，z)，f_y(x，y，z)，f_z(x，y，z))位置，回归到真实图像的(x，y，z)位置；

b)根据校正公式的图像幅度权重部分： $\frac{{\∂f}_x(x,y,z)}{\∂x}\frac{{\∂f}_y(x,y,z)}{\∂y}\frac{{\∂f}_z(x,y,z)}{\∂z},$ 与像素的幅度相乘，得到真实图像的幅度强度值。

对于用二维序列采集的图像，可在每个层面内运用二维畸变校正公式对图像进行校正；对于用三维序列采集的图像，需用三维畸变校正公式对三维图像进行校正。

Claims (1)

1.一种核磁共振成像梯度场校正方法，其特征在于，包括以下步骤：

（1）由核磁共振成像得到未校正的带有畸变的图像数据；

（2）从系统读出成像参数：包括频率编码、相位编码和选层，三个方向的视野以及空间分辨率，成像采集的旋转矩阵，成像中心偏离磁体中心的位置，计算得到每个像素点对应于以磁体中心为原点的绝对坐标系的坐标值；

（3）根据厂家给出的梯度线圈的谐波场分解系数，计算得到梯度场的空间非线性度函数；

（4）根据成像参数和梯度场的空间非线性度函数计算出梯度畸变校正公式；

（5）将所述梯度场的空间非线性度函数代入梯度畸变校正公式，对畸变图像进行校正；

（6）其中，所述空间非线性度函数定义为实际梯度场与理想梯度强度之比，所述梯度畸变校正公式的一维矫正公式为

二维校正公式为 $\mathrm{\ρ}(x,y)=\frac{\∂f_x(x,y)}{\∂x}\frac{\∂f_y(x,y)}{\∂y}s(f_x(x,y),f_y(x,y)),$ 三维校正公式为 $\mathrm{\ρ}(x,y,z)=\frac{\∂f_x(x,y,z)}{\∂x}\frac{{\∂f}_y(x,y,z)}{\∂y}\frac{{\∂f}_z(x,y,z)}{\∂z}s(f_x(x,y,z),f_y(x,y,z),f_z(x,y,z)),$ 其中，f_x(x)表示一维情况下的x轴方向的非线性度函数，f_x(x,y)和f_y(x,y)分别表示二维情况下的x、y轴方向的非线性度函数，f_x(x,y,z)、f_y(x,y,z)和f_z(x,y,z)分别表示三维情况下的x、y、z轴方向的非线性度函数，ρ(x)、ρ(x,y)和ρ(x,y,z)分别表示一维、二维、三维校正的真实图像，函数s表示畸变的图像。

Images