CN103576115A - 一种核磁共振图像梯度场变形校正方法 - Google Patents

Abstract

本发明公开了一种核磁共振图像梯度场变形校正方法，包括如下步骤：a)获取核磁共振图像信息，并将图像空间网格化；b)在图像坐标系下计算各网格点变形量和变形后的位置；c)结合网格点变形后的位置与重建后生成的图像，计算每个像素空间中的原始像素信号幅度值，回填进每个像素，得到校正后的图像。本发明提供的核磁共振图像梯度场变形校正方法，根据信号采集时的视野，图像矩阵信息将图像空间网格化，结合网格点变形后的位置与重建后生成的图像，计算原始采集时每个像素空间中的信号幅度值，回填进每个像素，得到校正后的图像，从而能最大限度的还原每一个采集像素内的信号，提高校正后图像的精度和信号均匀度。

Description

一种核磁共振图像梯度场变形校正方法

技术领域

本发明涉及一种图像变形后的校正方法，尤其涉及一种核磁共振图像梯度场变形校正方法。

背景技术

MR（核磁共振）成像过程中，因为梯度场的非线性，重建生成的图像会产生一定程度上的畸变，如图1和图2所示，因此需要对梯度场非线性造成的畸变图像进行校正。

专利申请号为201110192584.3，发明名称为一种核磁共振成像梯度场校正方法的中国专利，公开了一种采用近似的方法处理像素（以像素上的某一点来代表此像素），这样可以对畸变的图像进行部分校正，但是该校正方法用一阶导数对信号强度进行校正，是选取的一个像素区域的某一点的一阶导数值进行校正，用部分代替了整体，校正效果跟选取的参考点有关系，因此精度不高；同时该校正方法要分成位置校正和幅度校正两个步骤对一幅图像进行校正，增加了处理的困难。因此，有必要提供新的核磁共振图像梯度场变形校正方法。

发明内容

本发明所要解决的技术问题是提供一种核磁共振图像梯度场变形校正方法，能最大限度的还原每一个采集像素内的信号，提高校正后图像的精度和信号均匀度。

本发明为解决上述技术问题而采用的技术方案是提供一种核磁共振图像梯度场变形校正方法，包括如下步骤：a)获取核磁共振图像信息，并将图像空间网格化；b)在图像坐标系下计算各网格点变形量和变形后的位置；c)结合网格点变形后的位置与重建后生成的图像，计算每个像素空间中的原始像素信号幅度值，回填进每个像素，得到校正后的图像。

进一步地，所述步骤a）中图像信息包括视野信息和图像矩阵信息。

进一步地，所述步骤a)中图像空间网格点位置按下述式（1）、式（2）和式（3）计算：

RO(i,j,k)=-0.5*FoVro+i*FOVro/Xro    （1）

PE(i,j,k)=-0.5*FoVpe+j*FOVpe/Ype    （2）

SS(i,j,k)=-0.5*FoVss+k*FOVss/Zss    （3）

其中RO(i,j，k)表示第i行，第j列，第k层网格点在图像坐标系中RO方向位置，PE(i,j，k)表示第i行，第j列，第k层网格点在图像坐标系中PE方向位置，SS(i,j，k)表示第i行，第j列，第k层网格点在图像坐标系中SS方向位置，FOVro为频率编码方向的视野，FOVpe为相位编码方向的视野，FOVss为选层方向的视野，单位为毫米；Xro为频率编码方向的像素数，Ype为相位编码方向的像素数，Zss为选层方向的像素数，单位为pixel，i，j，k为正整数。进一步地，所述步骤b)中图像坐标系下，各网格点的变形量计算如下：b1)根据设备的磁场描述参数获取该设备在物理坐标系（x,y,z）下，各网格点x方向偏移量Δ_x，y方向偏移量Δ_y，z方向偏移量Δ_z；b2)根据物理坐标系和图像坐标系之间的转换矩阵M，计算图像坐标系中RO方向，PE方向和SS方向的网格点变形量。

进一步地，所述图像坐标系和物理坐标系之间的转换关系如下：

$\left(\begin{array}{c}\mathrm{ro}\\ \mathrm{pe}\\ \mathrm{ss}\end{array}\right)=M^{\′}\left(\begin{array}{c}x\\ y\\ z\end{array}\right),$ M为转换矩阵M的逆矩阵。

进一步地，所述步骤b2)中图像坐标系中网格点RO方向变形量Δ_ro，PE方向变形量Δ_pe和SS方向的变形量Δ_ss计算如下：

$\left(\begin{array}{c}{\mathrm{\Δ}}_{\mathrm{ro}}\\ {\mathrm{\Δ}}_{\mathrm{pe}}\\ {\mathrm{\Δ}}_{\mathrm{SS}}\end{array}\right)=M^{\′}\left(\begin{array}{c}{\mathrm{\Δ}}_x\\ {\mathrm{\Δ}}_y\\ {\mathrm{\Δ}}_z\end{array}\right).$

进一步地，所述步骤c)中每个像素空间中的原始像素信号幅度值按如下公式计算：其中i的取值为1-N，N为变形后的网格覆盖的重建后像素个数，w_i为第i个像素的权值，V_i为第i个像素的像素幅值。

进一步地，对于二维图像，所述像素的权值为校正后网格覆盖变形后网格的面积与重建后网格面积的比值。

进一步地，每个像素的权值为该像素点在x、y方向上权值的积与网格面积的比值。

进一步地，对于三维图像，所述像素的权值为校正后立方体网格覆盖变形后立方体网格的体积与重建后立方体网格体积的比值。

进一步地，每个像素的权值为该像素点在x、y、z方向上权值的积与立方体网格体积的比值。

进一步地，所述步骤c)中每个像素空间中的原始像素信号幅度值的计算方法如下：将未变形的网格划分为m*P个，取每个未变形网格中心点的像素幅值代表该网格的像素幅值，计算网格中心点变形后落在重建后的那一个像素里面，校正后的像素信号幅度值取重建后像素信号值的1/m，然后将来自同一个未变形网格的m个对应像素值之和作为校正后像素的像素值，P为未变形图像的像素个数，m为自然数。

本发明对比现有技术有如下的有益效果：本发明提供的核磁共振图像梯度场变形校正方法，根据信号采集时的视野，图像矩阵信息将图像空间网格化，利用磁场描述系数文件计算图像坐标系下各网格点变形量和变形后的位置，结合网格点变形后的位置与重建后生成的图像，计算原始采集时每个像素空间中的信号幅度值，回填进每个像素，得到校正后的图像，从而能最大限度的还原每一个采集像素内的信号，提高校正后图像的精度和信号均匀度。

附图说明

图1为本发明MR在理想梯度场下采集图像网格示意图；

图2为本发明由MR梯度场非线性导致的网格变形示意图；

图3为本发明MR图像梯度场变形校正方法流程图；

图4为本发明MR梯度场变形后的网格与重建后的像素网格重叠示意图。

图5为本发明MR梯度场变形后网格覆盖的重建后像素个数为4时像素信号幅度值计算示意图。

具体实施方式

下面结合附图和实施例对本发明作进一步的描述。

图3为本发明MR图像梯度场变形校正方法流程图。

请参见图3，本发明提供的核磁共振图像梯度场变形校正方法包括如下步骤：

步骤S301，获取核磁共振图像信息，包括视野信息和图像矩阵信息，将图像空间网格化；图像空间网格点位置按下述式（1）、式（2）和式（3）计算：

RO(i,j,k)=-0.5*FoVro+i*FOVro/Xro    （1）

PE(i,j,k)=-0.5*FoVpe+j*FOVpe/Ype    （2）

SS(i,j,k)=-0.5*FoVss+k*FOVss/Zss    （3）

其中RO(i,j，k)表示第i行，第j列，第k层网格点在图像坐标系中RO方向位置，PE(i,j，k)表示第i行，第j列，第k层网格点在图像坐标系中PE方向位置，SS(i,j，k)表示第i行，第j列，第k层网格点在图像坐标系中SS方向位置，FOVro为频率编码方向的视野，FOVpe为相位编码方向的视野，FOVss为选层方向的视野，单位为毫米；Xro为频率编码方向的像素数，Ype为相位编码方向的像素数，Zss为选择方向的像素数，单位为pixel，i，j，k为正整数。视野信息(FOVro,FOVpe，FOVss)，图像矩阵信息（Xro,Ype，Zss）由前端得到,即在图像重建组件前，通过界面设置得到的视野信息和图像矩阵信息。

步骤S302，计算图像坐标系下空间各网格点变形量和变形后的位置；空间各网格点变形量计算如下：

i）根据设备的磁场描述参数获取该设备在物理坐标（x,y,z）下点x方向偏移量Δ_x，y方向偏移量Δ_y，z方向偏移量Δ_z；

磁场设计部门会给出一套磁场描述参数，由此套参数可以得到设备在物理坐标下（x,y,z）三个方向二阶及二阶以上场相对于一阶场随空间位置变化的分布状况，如下式

Δ_x=f(x,y,z)    （4）

Δ_y=g(x,y,z)    （5）

Δ_z=k(x,y,z)    （6）

公式4表示在物理坐标系下，空间位置为(x,y,z)的点x方向偏移量为Δ_x；

公式5表示在物理坐标系下，空间位置为(x,y,z)的点y方向偏移量为Δ_y；

公式6表示在物理坐标系下，空间位置为(x,y,z)的点z方向偏移量为Δ_z。

ii)根据物理坐标系和图像坐标系之间的转换矩阵M，计算图像坐标系中RO方向，PE方向和SS方向的网格点变形量。

由信号采集时的相关信息可以得到图像平面在物理坐标系中的位置信息，以图像中心作为坐标原点，频率编码方向和相位编码方向及选层方向作为坐标系的方向可以创建出一个图像坐标系，由图像在物理坐标系中的位置可以确定出一个图像坐

标系到物理坐标系之间的转换矩阵M，如下面式（7）所示

利用上述仿射变换可以将公式4，公式5，公式6转换成图像坐标系下的转换，如下式：

Δ_x=f(ro,pe,ss)    (8)

Δ_y＝g(ro,pe,ss)   (9)

Δ_z=k(ro,pe,ss)    (10)

由矩阵M可以得到其逆矩阵M'，M'矩阵具有如下性质

$\left(\begin{array}{c}\mathrm{ro}\\ \mathrm{pe}\\ \mathrm{ss}\end{array}\right)=M^{\′}\left(\begin{array}{c}x\\ y\\ z\end{array}\right)---\left(11\right)$

利用公式11可以将公式（8），公式（9），公式（10）中对物理坐标x,y,z方向二阶及二阶以上磁场对一阶场随空间位置变化的分布情况转换成对图像坐标ro,pe,ss方向的描述，如下式（12）：

$\left(\begin{array}{c}{\mathrm{\Δ}}_{\mathrm{ro}}\\ {\mathrm{\Δ}}_{\mathrm{pe}}\\ {\mathrm{\Δ}}_{\mathrm{SS}}\end{array}\right)=M^{\′}\left(\begin{array}{c}{\mathrm{\Δ}}_x\\ {\mathrm{\Δ}}_y\\ {\mathrm{\Δ}}_z\end{array}\right)---\left(12\right)$

利用上述仿射变换可以将公式（8），公式（9），公式（10）变换成：

Δ_ro＝f′(ro,pe,ss)    （13）

Δ_pe＝g′(ro,pe,ss)    （14）

Δ_ss=k′(ro,pe,ss)     （15）

如果只做二维图像上的变形校正，公式（13），公式（14），公式（15）可以简化为：

Δ_ro＝f′(ro,pe)    （16）

Δ_pe＝g′(ro,pe)    （17）

将公式（1），公式（2）计算得到的图像坐标系下各网格点代入公式（16），公式（17）可以得到ro,pe方向各网格点的变形量。

步骤S303，结合网格点变形后的位置与重建后生成的图像，计算每个像素空间中的原始像素信号幅度值，回填进每个像素，得到校正后的图像。

每个像素空间中的原始信号幅度值按如下公式计算：

$\mathrm{value}=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{w}_i{V}_i---\left(18\right)$

其中i的取值为1-N，N为变形后的网格覆盖的重建后像素个数，w_i为第i个像素的权值，V_i为第i个像素的像素幅值。

对于二维图像，其变形后网格覆盖重建后网格中，各个像素点所取的权值，为校正后网格覆盖变形后网格的面积与重建后网格面积的比值；对于三维图像，变形后立方体网格覆盖重建后立方体网格中，各个像素点所取的权值，为校正后立方体网格覆盖变形后立方体网格的体积与重建后立方体网格体积的比值。

结合图4、图5所示，进一步说明对于二维图像，原始像素幅值的具体计算方法。图4左上角变形网格与重建后网格，在重叠区域内含有四个像素，因此，校正后左上角像素的信号值就应该为它所覆盖的四个像素信号值加权之和。加权的权值由变形网格在重建后的像素中所占面积比重决定。若变形后一个像素如图5中斜框所示，它的信号值由A，B，C，D四个像素值按比例合成，若校正后的每个像素可看作边长为a的正方形，则原始像素信号幅值计算如下：

value=c*b*A/a²+d*e*B/a²+f*j*C/a²+k*l*D/a²

其中b、c、d、e、f、j、k、l分别为所对应像素在ro、pe方向上的权值，而ro、pe方向上权值的乘积与一个像素面积的比值：c*b/a²,d*e/a²,f*j/a²,k*l/a²分别为变形后的像素信号幅度值占A，B，C，D像素值的权重，即所对应像素所占权值。以上虽然是针对二维的情况做的描述，但是同样可以推广到三维的情况。只是在对像素幅值进行分割的时候不只考虑ro,pe方向，还要考虑ss方向上的权值。三个方向上的权值之积与一个像素体积的比值就是拆分时分得的此像素幅值的比重，即所对应像素所占权值。

此外，本发明也可以对未变形的网格进行比像素网格更密集的划分,如未变形像素个数为P，则将未变形的网格划分为m*P个，取每个未变形网格中心点代表该网格，计算网格中心点变形后落在重建后的那一个像素里面，校正后的像素信号幅度值取重建后像素信号值的1/m，然后将来自同一个未变形网格的m个对应像素值之和作为校正后像素的像素值，m为自然数。

综上，本发明提供的核磁共振图像梯度变形校正方法，根据信号采集时的视野，图像矩阵信息将图像空间网格化，利用磁场描述系数文件计算空间各网格点变形量和变形后的位置，结合网格点变形后的位置与重建后生成的图像，计算原始采集时每个像素空间中的信号幅度值，回填进每个像素，得到校正后的图像，从而能最大限度的还原每一个采集像素内的信号，提高校正后图像的精度和信号均匀度。

虽然本发明已以较佳实施例揭示如上，然其并非用以限定本发明，任何本领域技术人员，在不脱离本发明的精神和范围内，当可作些许的修改和完善，因此本发明的保护范围当以权利要求书所界定的为准。

Claims (12)

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

a)获取核磁共振图像信息，并将图像空间网格化；

b)在图像坐标系下计算各网格点变形量和变形后的位置；

c)结合网格点变形后的位置与重建后生成的图像，计算每个像素空间中的原始像素信号幅度值，回填进每个像素，得到校正后的图像。

2.如权利要求1所述的核磁共振图像梯度场变形校正方法，其特征在于，所述步骤a)中图像信息包括视野信息和图像矩阵信息。

3.如权利要求2所述的核磁共振图像梯度场变形校正方法，其特征在于，所述步骤a)中图像空间网格点位置按下述式（1）、式（2）和式（3）计算：

RO(i,j,k)=-0.5*FoVro+i*FOVro/Xro    （1）

PE(i,j,k)=-0.5*FoVpe+j*FOVpe/Ype    （2）

SS(i,j,k)=-0.5*FoVss+k*FOVss/Zss    （3）

其中RO(i,j，k)表示第i行，第j列，第k层网格点在图像坐标系中RO方向位置，PE(i,j，k)表示第i行，第j列，第k层网格点在图像坐标系中PE方向位置，SS(i,j，k)表示第i行，第j列，第k层网格点在图像坐标系中SS方向位置，FOVro为频率编码方向的视野，FOVpe为相位编码方向的视野，FOVss为选层方向的视野，单位为毫米；Xro为频率编码方向的像素数，Ype为相位编码方向的像素数，Zss为选层方向的像素数，单位为pixel，i，j，k为正整数。

4.如权利要求2所述的核磁共振图像梯度场变形校正方法，其特征在于，所述步骤b)中图像坐标系下各网格点的变形量计算如下：

b1)根据设备的磁场描述参数获取该设备在物理坐标系（x,y,z）下，各网格点x方向偏移量Δ_x，y方向偏移量Δ_y，z方向偏移量Δ_z；

b2)根据物理坐标系和图像坐标系之间的转换矩阵M，计算图像坐标系中RO方向，PE方向和SS方向的网格点变形量。

5.如权利要求4所述的核磁共振图像梯度场变形校正方法，其特征在于，所述图像坐标系和物理坐标系之间的转换关系如下：

$\left(\begin{array}{c}\mathrm{ro}\\ \mathrm{pe}\\ \mathrm{ss}\end{array}\right)=M^{\′}\left(\begin{array}{c}x\\ y\\ z\end{array}\right),$ M′为转换矩阵M的逆矩阵。

6.如权利要求4所述的核磁共振图像梯度场变形校正方法，其特征在于，所述步骤b2)中图像坐标系中网格点RO方向变形量Δ_ro，PE方向变形量Δ_pe和SS方向的变形量Δ_SS计算如下：

$\left(\begin{array}{c}{\mathrm{\Δ}}_{\mathrm{ro}}\\ {\mathrm{\Δ}}_{\mathrm{pe}}\\ {\mathrm{\Δ}}_{\mathrm{SS}}\end{array}\right)=M^{\′}\left(\begin{array}{c}{\mathrm{\Δ}}_x\\ {\mathrm{\Δ}}_y\\ {\mathrm{\Δ}}_z\end{array}\right).$

7.如权利要求2所述的核磁共振图像梯度场变形校正方法，其特征在于，所述步骤c)中每个像素空间中的原始像素信号幅度值按如下公式计算：

其中i的取值为1-N，N为变形后的网格覆盖的重建后像素个数，w_i为第i个像素的权值，V_i为第i个像素的像素幅值。

8.如权利要求7所述的核磁共振图像梯度场变形校正方法，其特征在于，对于二维图像，所述像素的权值为校正后网格覆盖变形后网格的面积与重建后网格面积的比值。

9.如权利要求8所述的核磁共振图像梯度场变形校正方法，其特征在于，每个像素的权值为该像素点在x、y方向上权值的积与网格面积的比值。

10.如权利要求7所述的核磁共振图像梯度场变形校正方法，其特征在于，对于三维图像，所述像素的权值为校正后立方体网格覆盖变形后立方体网格的体积与重建后立方体网格体积的比值。

11.如权利要求10所述的核磁共振图像梯度场变形校正方法，其特征在于，每个像素的权值为该像素点在x、y、z方向上权值的积与立方体网格体积的比值。

12.如权利要求2所述的核磁共振图像梯度场变形校正方法，其特征在于，所述步骤c)中每个像素空间中的原始像素信号幅度值的计算方法如下：将未变形的网格划分为m*P个，取每个未变形网格中心点的像素幅值代表该网格的像素幅值，计算网格中心点变形后落在重建后的那一个像素里面，校正后的像素信号幅度值取重建后像素信号值的1/m，然后将来自同一个未变形网格的m个对应像素值之和作为校正后像素的像素值，P为未变形图像的像素个数，m为自然数。

Images