一种磁共振成像系统梯度场球谐系数的获取方法

摘要

本发明涉及一种磁共振成像系统梯度场球谐系数的获取方法，其步骤包括：(a)建立水模坐标系，以最接近图像中心的标志点Pc(u0，v0)为原点Op；以MRI图像中u轴方向为其c轴方向；以图像中v轴方向为其r轴方向；以与c、r轴构成右手坐标系的方向为s轴方向；(b)求标志点的MRI坐标，对标志点Pi，其MRI坐标可通过水模坐标系求得，其分为以下步骤：(i)不考虑水模在自身平面内的旋转和坐标原点的平移；(ii)考虑水模自身平面内以Op为圆心，以设置扫描的层slice方向为轴旋转θ角度；(iii)考虑水模坐标原点的平移；(c)如已知某点的磁场强度参数γ、θ、φ和图像坐标，则可得到关于球谐系数a、b的线性方程组，通过解线性方程组可以求出系数a、b。

说明书

技术领域

本发明涉及一种系数的获取方法，特别涉及一种磁共振成像系统梯度场球谐系数的获取方法。

背景技术

在核共振系统中，梯度场的非线性往往会引起图像的几何变形。不论是在不同的磁共振成像设备上还是相同磁共振成像设备而不同视野(FOV)内，图像的几何变形是不同的。因此，针对梯度场的非线性问题的图像校正技术对提高磁共振图像导航手术的精度以及多磁共振之间图像研究的可行性判断在实际应用中都具有相当的必要性。磁共振图像梯度变形校正的基本过程是：先通过某种方法求出磁共振成像设备空间中有限个控制点的位置偏移量，接下来计算待校正图像的每一个像素在磁共振成像设备空间内的坐标，通过插值方法根据每一个像素点的坐标计算该像素空间的偏移量；根据磁共振成像设备坐标系与图像坐标系的转换关系，将该像素空间偏移量换算为图像坐标系内的偏移量，之后将该偏移量补偿到该像素的图像坐标中，即完成了磁共振图像梯度变形校正过程。

计算磁共振成像空间内控制点的偏移量的方法可以分为两种：第一种是使用三维立体水模，三维立体水模提供了三维空间内的控制点，通过对立体水模的磁共振图像内控制点的分析，可以得到有限空间控制点上的图像变形量，根据磁共振成像设备坐标系与图像坐标系的转换关系，将该像素的图像坐标系内的偏移量换算为空间偏移量；另一种方法是应用磁场的描述函数-球谐函数来计算成像空间内有限个控制点的偏移量。由于第一种方法依赖三维水模，实施成本远高于第二种方法，而且需要对三维水模多次成像和分析才能获得控制点的位置偏移量，所以实施的效率很低。在应用第二种方法，即球谐函数进行计算之前，首先要给出球谐函数内的相关校正参数。哈佛大学的研究人员Jorge Jovicich、SilvesterCzanner在他们的论文，Reliability in Multi-Site Structural MRI Studies：Effects of Gradient Non-linearity Correction on Phantom and Human Data，Neuron Image(in press)(多点结构性磁共振成像可靠性研究：梯度非线性校正在水模和人体图像中的效果，《神经图像》中文译文)中提供了计算相关校正参数的方法。该方法是根据梯度线圈的设计参数计算用于图像校正的校正参数，但梯度线圈的设计参数只是图像变形的主要因素之一，此外还有其他因素没有被考虑进去。所以，实际的校正参数与设计参数之间仍存在一定差异。

发明内容

针对上述问题，本发明的目的是提供一种针对特定的梯度线圈在图像空间内取有限的点精确计算校正参数，从而获取真实的磁场梯度参数的磁共振成像系统梯度场球谐系数的获取方法。

为实现上述目的，本发明采取以下技术方案：一种磁共振成像系统梯度场球谐系数的获取方法，是根据从实际图像空间内取有限的点，反求磁共振校正参数即球谐系数a和b，其包括以下步骤：

(a)建立水模坐标系，MRI图像原点为V(x₀，y₀，z₀)，以水模图像内最接近MRI图像中心V(x_c，y_c，z_c)的标志点P_c(u₀，v₀)为水模坐标系原点O_p；以MRI图像中u轴方向为其c轴方向；以图像中v轴方向为其r轴方向；以与c、r轴构成右手坐标系的方向为s轴方向；任意标志点P_i的水模坐标系坐标为(c_i，r_i，s_i)；根据O_p的图像坐标计算其MRI坐标，得到

$V_{o_p}(x,y,z)=v_c*u_0*p*v_r*v_0*p+V(x_0,y_0,z_0)---\left(1\right)$

其中p为相邻像素之间的距离。v_c，v_r，v_s分别为c、r和s轴方向的单位向量。

(b)求标志点的MRI坐标，对标志点P_i(c_i，r_i，s_i)，其MRI坐标可通过水模坐标系分为以下步骤：

(i)不考虑水模在自身平面内的旋转和坐标原点的平移，标志点P_i(c_i，r_i，s_i)在MRI坐标系中的坐标为

$V_{p_i}(x,y,z)=c_i*d*v_c+r_i*d*v_r+s_i*d*v_s---\left(2\right)$

(ii)考虑水模自身平面内以O_p为圆心，以设置扫描的层方向为轴旋转θ角度，标志点P_i在MRI坐标系中的坐标需要再乘以旋转矩阵T_rot

$V_{p_i}=(x^{\prime },y^{\prime },z^{\prime },1)=V_{p_i}(x,y,z,1)*T_{\mathrm{rot}}---\left(3\right)$

$T_{\mathrm{rot}}=\left(\begin{array}{cccc}t*x*x+c& t*x*y+s*z& t*x*z-s*y& 0\\ t*x*y-s*z& t*y*y+c& t*y*z+s*x& 0\\ t*x*z+s*y& t*y*y+c& t*z*z+c& 0\\ 0& 0& 0& 1\end{array}\right)$

c＝cos(θ)，s＝sin(θ)，t＝1-c

其中x，y，z是层方向向量的坐标；

(iii)考虑水模坐标原点的平移，标志点P_i在MRI坐标系中的坐标为

$V_{p_i}=(x^{\prime \prime },y^{\prime \prime },z^{\prime \prime })=V_{p_i}(x^{\prime },y^{\prime }{z}^{\prime })+V_{o_p}(x,y,z)---\left(4\right)$

其中，V_pi(x″，y″，z″)即为标志点P_i在MRI坐标系中的实际坐标。磁场强度可以表达为：

B_r(n，m)(r，θ，_)＝rⁿ[a_v(n，m)cos(m_)+b_v(n，m)sin(m_)]×P_(n，m)(cosθ)    (5)

磁场梯度函数：

$G_{v\left(r\right)}\equiv \frac{{\partial \mathrm{\&Bgr;}}_{\mathrm{zv}\left(r\right)}}{\partial v}\equiv \frac{\partial {\mathrm{\&Bgr;}}_{\mathrm{zv}\left(r\right)}^L}{\partial v}+\frac{{\partial \mathrm{\&Bgr;}}_{\mathrm{zv}\left(r\right)}^N}{\partial v}\equiv G_v^L+G_{v\left(r\right)}^L---\left(6\right)$

由公式(5)(6)(7)司推出：

V_x＝∑rⁿ[a_x(n，m)cos(m_)+b_x(n，m)sin(m_)]×P_(n，m)(cosθ)/a_(1，1)

V_y＝∑rⁿ[a_y(n，m)cos(m_)+b_y(n，m)sin(m_)]×P_(n，m)(cosθ)/b_(1，1)    (8)

V_z＝∑rⁿ[a_z(n，m)cos(m_)+b_z(n，m)sin(m_)]×P_(n，m)(cosθ)/a_z(1，0)

其中a_v(n，m)、b_v(n，m)是常数，且a_v(n，m)、b_v(n，m)是v方向n阶m级展开项的系数，P_(n，m)(cosθ)为勒让德多项式；公式(8)左端是标志点P_i在MRI坐标系下的坐标值；

(c)如上所述，如已知某点的γ、θ、_和MRI图像坐标，则可得到关于球谐系数a、b的线性方程组，通过解线性方程组可以求出系数a、b。

所述步骤(c)中γ、θ、_是标志点在MRI坐标系中的坐标。可以由V_pi(x″，y″，z″)通过以下换算得到

所述步骤(a)中相邻标志点为等距设置，其间隔d取定值。

所述步骤(b)中，设P₁(u₁，v₁)为水模上c方向上与P_c(u₀，v₀)相邻的标志点；根据P₁(u₁，v₁)和P_c(u₀，v₀)的图像计算θ角，则θ角的计算公式

θ＝tg^-1((v₁-v₀)/(u₁-u₀))。

本发明由于采取以上技术方案，其具有以下优点：1、本发明根据实际图像反求磁场梯度参数，克服了梯度线圈设计参数与真实的磁场梯度参数之间存在差异的缺点。2、本发明在反求磁场梯度参数的过程中，同时考虑了水模摆放位置对算法精确度的影响，提高了算法的准确率。3、本发明采用梯度场球谐参数的搜索算法，大大简化了算法实施的步骤，提高了算法效率。本发明方法针对在磁共振成像应用中由梯度场非线性导致的图像变形校正，对提高磁共振图像导航手术的精度和多磁共振之间图像研究的可行性分析具有重要的意义。

附图说明

图1是本发明水模结构示意图

图2是图1的侧视剖视示意图

图3是本发明扫描得到的水模图像以及基于此建立的水模坐标系示意图

具体实施方式

下面结合附图和实施例，对本发明进行详细的描述。

本发明所提供的磁共振成像系统梯度场球谐系数的获取方法，是根据实际图像空间内取有限的点，反求磁共振校正参数也即球谐系数，出于计算量和校正精确度的综合考虑，仅进行5阶以下的参数计算。

如图1、2所示，梯度校正水模是用于梯度校正的工具，它既用于搜索系统参数，也用于测量校正后的误差。梯度校正水模为一个方形的盒体1，在盒体1内设置按照等间距方阵的方式排列成方阵的圆柱2，圆柱2内充满硫酸铜溶液，可以在磁共振成像设备中成像。

将水模放置在有效的成像空间中，保证磁场中心在水模内，设置扫描平面使扫描得到的水模图像，如图3所示。在水模图像内最接近MRI图像中心V(x_c，y_c，z_c)的标志点P_c，其水模图像坐标为(u₀，v₀)。MRI图像原点为V(x₀，y₀，z₀)，此处V(x，y，z)中的x，y，z均为MRI坐标系下坐标。建立水模坐标系时，以标志点P_c为坐标系原点O_p，以MRI图像中u轴方向为c轴方向，以图像中v轴方向为r轴方向，以与c、r轴构成右手坐标系的方向为s轴方向。由于按照所述水模摆放方式以及扫描的要求，水模原点O_p距离图像中心V(x_c，y_c，z_c)很近，因此可以忽略由于梯度非线性引起的图像变形对O_p位置的影响，那么可以直接根据O_p的图像坐标计算其MRI坐标如下所示

$V_{o_p}(x,y,z)=v_c*u_0*p+v_r*v_0*p+V(x_0,y_0,z_0)---\left(1\right)$

其中p为相邻像素之间的距离，v_c，v_r，v_s分别为c、r和s轴方向的单位向量。

而对于任意一个标志点P_i，由于其位置的梯度非线性引起的图像变形比较显著，无法再使用图像坐标计算其MRI坐标，所以P_i的MRI坐标需通过水模坐标系计算。标志点P_i的水模坐标系坐标为(c_i，r_i，s_i)，相邻标志点的间隔为d，计算可分为以下步骤：

(i)不考虑水模在自身平面内的旋转和坐标原点的平移，标志点P_i在MRI坐标系中的坐标为

$V_{p_i}(x,y,z)=c_i*d*v_c+r_i*d*v_r+s_i*d*v_s---\left(2\right)$

(ii)考虑水模自身平面内以O_p为圆心，以设置扫描的层slice方向为轴旋转θ角度，标志点P_i在MRI坐标系中的坐标需要再乘以旋转矩阵T_rot，

$V_{p_i}=(x^{\prime },y^{\prime },z^{\prime },1)=V_{p_i}(x,y,z,1)*T_{\mathrm{ros}}---\left(3\right)$

$\left(T_{\mathrm{rot}}=\left(\begin{array}{cccc}t*x*x+c& t*x*y+s*z& t*x*z-s*y& 0\\ t*x*y+s*z& t*y*y+c& t*y*z+s*x& 0\\ t*x*z+s*y& t*y*z-s*x& t*z*z+c& 0\\ 0& 0& 0& 1\end{array}\right)\right)$

c＝cos(θ)，s＝sin(θ)，t＝1-c

其中x，y，z是slice方向向量的坐标；设P₁(u₁，v₁)为水模上c方向上与P_c(u₀，v₀)相邻的标志点；根据P₁(u₁，v₁)和P_c(u₀，v₀)的图像计算θ角，则θ角的计算公式

θ＝tg^-1((v₁-v₀)/(u₁-u₀))；

(iii)考虑水模坐标原点的平移，标志点P_i在MRI坐标系中的坐标为

$V_{p_i}=(x^{\prime \prime },y^{\prime \prime },z^{\prime \prime })=V_{p_i}(x^{\prime },y^{\prime }{z}^{\prime })+V_{o_p}(x,y,z)---\left(4\right)$

V_pi，(x″，y″，z″)即为标志点P_i在MRI坐标系中的实际坐标。

由公式(5)(6)(7)即

磁场强度

B_r(n，m)(r，θ，_)＝rⁿ[a_v(n，m)cos(m_)+b_v(n，m)sin(m_)]×P_(n，m)(cosθ)(5)

磁场梯度函数

$G_{v\left(r\right)}\equiv \frac{{\partial \mathrm{\&Bgr;}}_{\mathrm{zv}\left(r\right)}}{\partial v}\equiv \frac{\partial {\mathrm{\&Bgr;}}_{\mathrm{zv}\left(r\right)}^L}{\partial v}+\frac{{\partial \mathrm{\&Bgr;}}_{\mathrm{zv}\left(r\right)}^N}{\partial v}\equiv G_v^L+G_{v\left(r\right)}^L---\left(6\right)$

由公式(5)(6)(7)可推出

V_x＝∑rⁿ[a_x(n，m)cos(m_)+b_z(n，m)sin(m_)]×P_(n，m)(cosθ)/a_x(1，1)

V_y＝∑rⁿ[a_y(n，m)cos(m_)+b_y(n，m)sin(m_)]×P_(n，m)(cosθ)/b_y(1，1)    (8)

V_z＝∑rⁿ[a_z(n，m)cos(m_)+b_z(n，m)sin(m_)]×P_(n，m)(cosθ)/a_z(1，0)

公式(8)的左端是某一成像点在MRI坐标系下变形后的坐标值，可以通过成像点的图像坐标计算得到。其中a_v(n，m)、b_v(n，m)是常数，a_v(n，m)、b_v(n，m)是v方向n阶m级展开项的系数，是磁场非线性梯度的固有特性。P_(n，m)(cosθ)为勒让德多项式。右端中的γ、θ、_是该点在MRI坐标系(球坐标系)中的坐标。可以由V_pi(x″，y″，z″)通过公式(9)换算得到。

如上所述，已知某点的γ、θ、_和图像坐标，则可得到关于a、b的线性方程组，通过解线性方程组可以求出系数a、b。接下来可以通过标志点P_i的图像坐标(u，v)计算该点变形后的MRI坐标V_pi(x，y，z)。

$V_{p_i}(x,y,z)=u*v_c*p+v*v_r*p+V(x_0,y_0,z_0)---\left(10\right)$

出于计算量和矫正精确度的综合考虑，仅进行5阶以下的参数计算。在球谐函数中，参数a_v(n，m)，b_v(n，m)在5阶以下仅有部分值非零，在计算中零项可以不予考虑，则非零项在表1中列出：

表1

   X    a    b  (n＝1，m＝1)    1    0  (n＝3，m＝1)    10^-4    0  (n＝5，m＝1)    10^-7    0   Y    a    b  (n＝1，m＝1)    0    1  (n＝3，m＝1)    0    10^-4  (n＝5，m＝1)    0    10^-7  Z    a    b  (n＝1，m＝0)    1    0  (n＝3，m＝0)    10^-4    0  (n＝5，m＝0)    10^-7    0

注：其中的a、b非零值均为估计值。

以计算a_x、b_x为例，公式(8)可以简化为：

         (11)

令拉格朗日函数公式(11)可以写成矩阵相乘的形式：

$V_x=\left(\begin{array}{ccc}\cos \left(\theta \right)*\gamma *P_1& \cos \left(\theta \right)*{\gamma }^3*P_3& \cos \left(\theta \right)*{\gamma }^5*P_5\end{array}\right)*\left(\begin{array}{c}{a}_{x\left(\mathrm{1,1}\right)}\\ a_{x(3,1)}\\ a_{x\left(\mathrm{5,1}\right)}\end{array}\right)---\left(12\right)$

以不同标志点的γ、θ、_，V_pi(x，y，z)代入公式(12)可以得到线性方程组，

通过解线性方程组就可以计算得到：$\left(\begin{array}{c}{a}_{x\left(\mathrm{1,1}\right)}\\ a_{x\left(\mathrm{3,1}\right)}\\ a_{x\left(\mathrm{5,1}\right)}\end{array}\right).$由于y、z方向的参数计算方法与此类似，所以在此不再赘述。