基于图像轮廓的图像畸变校正方法和装置

摘要

本发明实施例公开了一种图像畸变校正方法和装置，该方法包括步骤：提取源图像的图像轮廓；根据图像轮廓获取畸变信息；以及采用畸变信息对源图像进行校正，以与目标图像配准。采用本发明方法和装置，避免了磁共振扩散加权图像中对比度、b值和方向变化等对校正的影响，并且不需要额外的参考扫描，实现了方便、快速、准确的图像畸变校正。

说明书

技术领域

本发明涉及扩散加权成像领域，具体涉及对扩散加权成像形成的图像进行畸变校正的方法和装置。

背景技术

扩散成像是一种非常具有临床价值的功能磁共振成像技术，它可以探测水分子或可以在磁共振成像设备中实现扩散成像的其它分子在人体组织中扩散的生物学特性，从而检测许多常规医学影像手段不能发现的病变(非专利文献1)。扩散成像通常使用回波平面成像(EPI)序列来实现，EPI序列成像速度非常快，但是容易受到各种伪影尤其是图像畸变的影响(非专利文献2和3)，EPI序列中，由于主磁场的非均匀性等引起的额外磁场ΔB在频率编码方向和相位编码方向上引起的图像畸变可以分别表达为(非专利文献4和5)：频率编码方向上的畸变通常很小，可以忽略，一般只考虑在相位编码方向上引起的畸变Δy。

扩散加权回波平面成像(DW-EPI)序列是通过在EPI序列中加入扩散梯度来实现扩散加权的。增加的扩散梯度会在成像空间中引起额外的涡流场，这种额外的涡流场是图像畸变产生的主要原因。令表示涡流场在主磁场B₀方向上的分量，可以展开为(非专利文献3)：

$>B_e(\overrightarrow{r},t)=b_0\left(t\right)+\overrightarrow{r}·\overrightarrow{g}\left(t\right)+...---\left(1\right)$>

其中第一项b₀(t)与空间位置无关的，通常称为零阶项，第二项为线性项，或者称为一阶项，分量为gx，gy，gz，分别表示在x，y，z三个方向上涡流场的梯度。更高阶的项通常比较小，可以忽略。此外，线性项中的z分量虽然会造成选层相位回聚的不完善，从而使得图像整体信号强度下降，降低图像的对比度，但是不会造成图像畸变。因此，对于图像畸变的校正，通常只考虑零阶项和一阶项中的x，y分量，从而上式可以重写为：

B_e(x，y)＝b₀+x·g_x+y·g_y    (2)

代入图像畸变Δy的表达式中，可得到由涡流场引起的图像畸变的表达式：

$>\mathrm{\Delta y}(x,y)$>

$>=B_e(x,y)t_{\mathrm{esp}}/\overline{G_y}\tau$>

$>=(b_0+x·g_x+y·g_y)·t_{\mathrm{esp}}/\overline{G_y}\tau ---\left(3\right)$>

从上式可以看出，零阶涡流场会引起图像整体的位移，沿着频率编码方向的一阶项会使图像产生剪切形变，沿着相位编码方向的一阶项使图像产生伸缩形变。

扩散加权成像中需要在成像序列中加入扩散梯度，作为测量扩散系数等参数的手段。b值是与所加扩散梯度脉冲形状和强度相关的集总参数，通常b值越大，梯度脉冲的宽度和高度也大对最后形成的图像对比度越好；然而，大b值需要幅度大、持续时间长的扩散梯度脉冲，因此会引起更大的涡流场，导致更为严重的图像畸变。扩散加权成像通常需要采集b值不同和方向不同的多个图像，由此，当在脉冲序列中使用与b值相关的不同扩散梯度参数值时，会引起不同的涡流场分布，由此导致的图像畸变程度和类型不完全相同。因此，对这些图像畸变的校正需要逐个图像分别进行，才能得到清晰的扩散加权成像(DWI)图像，使得表观扩散系数图计算(ADC mapping)更加准确。

抑制涡流场的方法有使用自屏蔽的梯度线圈和梯度脉冲的预加强(非专利文献6-8)。自屏蔽线圈包括一个主线圈和一个屏蔽线圈，屏蔽线圈在成像时产生的磁场的方向与主线圈相反，经过特殊设计，使得这一对线圈在成像区域内产生期望的梯度磁场，而在成像区域外的磁场尽可能接近零，这就从根本上减小了涡流场的分布和大小。而预加强技术是通过改变梯度波形的形状，来抵消涡流场的影响。通过使用这种方法，可以减小涡流场的影响，但是并不能完全消除。而扩散加权成像对涡流场非常敏感，为了提高成像的质量，必须对得到的图像进行畸变校正。

然而，如上所述，扩散加权成像通常需要采集b值不同的多个图像，造成图像对比度差异较大；此外，不同的b值和方法造成不同的涡流场分布，导致的图像畸变程度和类型不完全相同，不得不逐个图像分别进行图像畸变校正，这会带来额外的参考扫描，增加扫描时间。

需要一种方便、快速、误差小的图像畸变校正方法。

非专利参考文献1：包尚联，现代医学影像物理学，北京：北京大学医学出版社，2004

非专利参考文献2：俎栋林，核磁共振成像学，北京：高等教育出版社，2004

非专利参考文献3：Bernstein M，King K，Zhou X，Fong W.Handbook of MRI pulse sequences.San Diego：Elsevier Academic Press；2004.

非专利参考文献4：Jezzard P，Balaban RS.Correction for geometricdistortion in echo-planar images from B-0field variations.MagneticResonance in Medicine 1995；34(1)：65-73.

非专利参考文献5：Zeng HR，Constable RT.Image distortioncorrection in EPI：Comparison of field mapping with point spreadfunction mapping.Magnetic Resonance in Medicine 2002；48(1)：137-146.

非专利参考文献6：Mansfield P，Chapman B.Active magneticscreening of coils for static and time-dependent magnetic-field generationin NMR imaging.Journal of Physics E-Scientific  Instruments1986；19(7)：540-545.

非专利参考文献7：Vanvaals JJ，Bergman AH.Optimization ofeddy-current compensation.Journal of Magnetic Resonance1990；90(1)：52-70.

非专利参考文献8：Jehenson P，Westphal M，Schuff N.Analyticalmethod for the compensation of eddy-current effects induced by pulsedmagnetic-field gradients in NMR systems.Journal of Magnetic Resonance1990；90(2)：264-278.

发明内容

为解决上述问题，本发明提供了一种图像畸变校正方法和装置，能够基于图像轮廓来得到用于校正的畸变信息，避免了图像对比度、b值和方向变化等对校正的影响，并且不需要额外的参考扫描，实现了方便、快速、准确的图像畸变校正。

根据本发明一方面，提供了一种图像畸变校正方法，包括如下步骤：

a.提取源图像的图像轮廓；

b.根据图像轮廓获取畸变信息；

c.采用畸变信息对源图像进行校正，以与目标图像配准。

根据本发明实施例，步骤b包括：将源图像的图像轮廓与目标图像的图像轮廓进行比较，计算位移系数和伸缩系数，作为畸变信息。

根据本发明实施例，该方法在步骤a之前还包括：

d.对源图像进行直方图均衡化处理；

其中，在步骤a中，提取处理后的源图像的图像轮廓。

根据本发明实施例，在步骤a中，使用水平集方法来提取源图像的图像轮廓。

根据本发明实施例，在步骤b中，针对图像轮廓的所有像素行，逐行地计算位移系数和伸缩系数，以得到整个源图像的像素偏移图，作为畸变信息。

根据本发明实施例，在步骤b中，针对图像轮廓的一部分像素行，逐行地计算位移系数和伸缩系数；

其中，在步骤b与步骤c之间所述方法还包括：将所述一部分像素行的位移系数和伸缩系数进行拟合并外推到其他像素行，以得到整个源图像的像素偏移图，作为畸变信息。

根据本发明实施例，源图像是扩散加权成像形成的图像或图像序列。

根据本发明另一方面，提供了一种图像畸变校正装置，包括：

提取单元，提取源图像的图像轮廓；

获取单元，根据图像轮廓获取畸变信息；

校正单元，采用畸变信息对源图像进行校正，以与目标图像配准。

根据本发明实施例，获取单元将源图像的图像轮廓与目标图像的图像轮廓进行比较，计算位移系数和伸缩系数，作为畸变信息。

根据本发明实施例，该装置还包括：

均衡化单元，对源图像进行直方图均衡化处理；

其中，提取单元提取均衡化单元处理后的源图像的图像轮廓。

根据本发明实施例，针对图像轮廓的所有像素行，获取单元逐行地计算位移系数和伸缩系数，以得到整个源图像的像素偏移图，作为畸变信息。

根据本发明实施例，针对图像轮廓的一部分像素行，获取单元逐行地计算位移系数和伸缩系数；

其中，所述装置还包括：外推单元，将获取单元计算的所述一部分像素行的位移系数和伸缩系数进行拟合并外推到其他像素行，以得到整个源图像的像素偏移图，作为畸变信息。

本发明的图像畸变校正方法和装置基于图像轮廓来获得畸变信息，用于校正，不需要图像有相似的对比度，从而能够用于b值较大的情况；此外，不需要序列的改动和额外的参考扫描，不会增加扫描时间。校正后的图像可以用于计算各向同性的DWI图像和ADC图，消除了这些图像中的伪影。

附图说明

下面的附图表明了本发明的实施方式。这些附图和实施方式以非限制性、非穷举性的方式提供了本发明的一些实施例。

图1示出了根据本发明实施例的图像畸变校正方法的流程图；

图2示出了目标图像和源图像之间的位移和缩放关系；

图3示出了根据本发明实施例的图像畸变校正方法中图像轮廓提取的计算机模拟图；

图4是本发明实施例的图像畸变校正装置的示意结构图；

图5示出了采用本发明实施例的图像畸变校正方法的校正结果的示意图；以及

图6示出了经过本发明实施例的图像畸变校正后的ADC图。

具体实施方式

下面参照附图，详细说明本发明实施例的图像畸变校正方法和装置。在接下来的说明中，一些具体的细节，例如实施例中的具体元件、步骤、具体参数等，都用于对本发明的实施例提供更好的理解。本发明不限于这些具体的元件、步骤或参数。此外，省略了对现有技术中公知的算法、手段等的具体描述。本技术领域的技术人员可以理解，即使在缺少一些细节或者其他元件、步骤、参数等结合的情况下，本发明的实施例也可以被实现。

为了清楚、简要地描述本发明实施例，以下主要以磁共振成像(MRI)中扩散加权的回波平面成像(DW-EPI)技术为例进行描述。然而，本发明实施例不限于此，而是适用于各种扩散加权成像形成的带有畸变或伪影的图像的校正，包括其他类型的扩散加权成像(DWI)、表观扩散系数计算(ADC mapping)、扩散张量成像(DTI)、纤维束跟踪等。

图1示出了根据本发明实施例的图像畸变校正方法的流程图，图2示出了目标图像I₀(x，y)和源图像I₁(x，y)之间的位移和缩放关系。使用图像畸变校正方法100对来自DW-EPI序列的源图像进行校正，以消除源图像中由于扩散加权或扩散梯度引起的畸变，使之与目标图像配准。图像畸变校正方法100包括：步骤102，提取源图像的图像轮廓；步骤104，根据图像轮廓获取畸变信息；以及步骤106，采用畸变信息对源图像进行校正，以与目标图像配准。

下面结合图2对本发明实施例的图像畸变校正方法100进行具体描述。图2中，a指示了目标图像I₀(x，y)，为b值为零的、无畸变的图像，b指示了源图像I₁(x，y)，为b值不为零的、带有畸变的图像，其中，x表示沿着频率编码方向的像素坐标，y表示沿着相位编码方向的像素坐标。在校正方法100的步骤102中，采用水平集(1evel set)方法提取源图像I₁(x，y)的轮廓。对于水平集方法的具体描述，可以参见非专利文献9(Chan TF，Vese LA.Active contours without edges.IEEETransactions on Image Processing，2001；10(2)：266-277.)。本领域技术人员可以理解，本发明实施例不限于水平集轮廓提取方法，而是可以采用本领域技术公知的其他任何适合的轮廓提取方法来提取图像轮廓。

在步骤102中，目标图像I₀(x，y)是预先设定的，并且目标图像I₀(x，y)的轮廓可以是预先提取的。备选地，在步骤102中，也可以与源图像I₁(x，y)的轮廓提取一起，来提取目标图像I₀(x，y)的轮廓。

以图2所示图像为例，本领域技术人员可以清楚看到，目标图像I₀(x，y)和源图像I₁(x，y)的轮廓分别是圆形和椭圆形的，源图像I₁(x，y)的轮廓相对于目标图像I₀(x，y)的轮廓而言发生了位移，剪切和伸缩，即，如背景技术部分所述，扩散加权中涡流场引起了整体位移，剪切和伸缩。如背景技术部分所述，频率编码方向上的畸变通常很小，可以忽略，一般只考虑在相位编码方向上引起的畸变。此外，对得到的图像轮廓进行逐行处理，由此消除了剪切畸变的影响，只需考虑每个像素行的位移和伸缩两种形变，把I₁(x)进行适量的伸缩和平移，就可以与I₀(x)对齐：

I₀(x)＝S(x)·I₁(x)+T(x)    (4)

其中，S(x)，T(x)分别表示第x行的伸缩系数和位移系数。

基于以上分析，在本实施例中，在步骤104，将源图像的图像轮廓与目标图像的图像轮廓进行比较，计算位移系数和伸缩系数，作为畸变信息，来用于校正。

具体而言，参见图2，考虑目标和源图像中沿着相位编码方向的像素行，分别用I₀(x)和I₁(x)表示，图2a中的实线指示了I₀(x)，图2b中的虚线指示了I₁(x)，图2c指示了目标图像的像素行I₀(x)和源图像的像素行I₁(x)之间的畸变信息的比较，从图2c所示的像素行的左右两个端点之间的线段y₀y₀’与y₁y₁’的比较可以清楚看出，y₁y₁’相对于y₀y₀’发生了向右平移，并且长度发生了缩放形变，这在图2c中由“缩放”和“位移”示出。将y₀、y₀’与y₁、y₁’带入上式(4)，得到S(x)，T(x)为：

$>S\left(x\right)=\frac{(y_0^{\prime }-y_0)}{(y_1^{\prime }-y_1)}$>

$>T\left(x\right)=\frac{(y_0^{\prime }+y_0)-(y_1^{\prime }+y_1)}{2}---\left(5\right)$>

这样，逐行计算出伸缩和位移系数，作为畸变信息对图像的畸变进行校正。

根据本发明一个实施例，在步骤104，可以基于提取的图像轮廓，针对源图像I₁(x，y)图像轮廓的所有像素行，逐行地计算位移系数和伸缩系数，以得到整个源图像I₁(x，y)的像素偏移图，作为畸变信息。结合以上示例，用S(x)和T(x)如下计算出像素偏移图psm(x，y)：

psm(x，y)＝(y-y_center)·(S(x)-1)+T(x)    (6)

其中，y_center是第x行中心像素的坐标。得到像素偏移图psm(x，y)后，在步骤106，就可以使用插值的方法对畸变的源图像I₁(x，y)进行校正，例如使用插值方法。

根据本发明另一实施例，进一步考虑到图像中会存在除涡流场引起的畸变之外的其它畸变或噪声。例如，假设图像共有256行，其中的若干行是背景噪声，例如第1-20行和240-256行，那么获取畸变信息时最好只使用第21-239行，以减小其他畸变或噪声带来的误差。此外，图像轮廓提取也会引入新的噪声或误差。为了进一步减少逐行计算伸缩和位移系数过程中的误差，在只考虑零阶和一阶线性误差的条件下，可以用对S(x)、T(x)进行拟合，然后外推到图像中所有的行。据此，备选地，在步骤104，只针对源图像图像轮廓的一部分像素行，逐行地计算位移系数和伸缩系数。然后，在步骤104与步骤106之间，校正方法100还包括如下步骤：将该一部分像素行的位移系数和伸缩系数进行拟合并外推到其他像素行，以得到整个源图像的像素偏移图，作为畸变信息。这里，拟合和外推可以采用本领域公知的任何适当方法，考虑到S(x)与x的线性关系，优选的是线性拟合方法，例如最小二乘法。

用S′(x)，T′(x)分别表示拟合和外推之后的伸缩和位移系数，按照下式计算出像素偏移图psm(x，y)：

psm(x，y)＝(y-y_center)·(S′(x)-1)+T′(x)    (7)

其中，y_center是第x行中心像素的坐标。同样，得到像素偏移图psm(x，y)后，在步骤106，就可以对畸变的源图像I₁(x，y)进行校正，例如与源图像I₁(x，y)相加，来得到与目标图像I₀(x，y)配准的图像。

由于根据本发明实施例的图像畸变校正方法100是基于图像轮廓的，所以准确的图像轮廓提取是精确计算像素偏移图的前提。对于这一前提的重要性，结合图3所示示例进行说明。

在使用体模和健康志愿者对校正方法100进行测试时，利用水平集方法能够准确提取出图像轮廓。但是，在对有脑部肿瘤的病人进行扫描时，由于肿瘤在扩散加权图中表现为高信号，可能会影响到图像边界的提取，这时需要用到一些图像处理技术，比如直方图均衡化等，降低肿瘤与正常组织的对比度，以便正确找到图像边界。图3示出了根据本发明实施例的图像畸变校正方法100中图像轮廓提取的计算机模拟图，图3a是计算机模拟产生的Shepp-Logan模型，其中一个椭圆被设置为高信号，用来模拟扩散成像时的肿瘤区域，另外图像中加入了高斯分布的噪声；图3b是直接用水平集方法提取出的轮廓(如白色箭头所指)，此时得到的并非是想要的结果，因为轮廓仅仅包括肿瘤区域；图3c是对图3a进行直方图均衡化后的结果，图3d是对图3c的图像进行轮廓提取，得到了实际需要的图像轮廓，如图中白色箭头所示。

因此，在本发明实施例的图像畸变校正方法100中，为了准确地提取图像轮廓，在步骤102之前，对源图像和/或目标图像进行直方图均衡化处理，然后在步骤102中对处理后的源图像和/或目标图像进行轮廓提取。由此，可以进一步减小轮廓提取引起的误差，提高像素偏移和校正的精确性。

以上描述了根据本发明实施例的图像畸变校正方法100的主要或备选步骤。对于扩散加权成像形成的图像序列中的每一源图像，可以针对扫描中所有的b值和所有方向，进行上述处理，然后再基于处理后的图像计算扩散加权图像或ADC值，就可以很大程度上减小图像扩散梯度涡流引起的图像模糊，使得ADC值的计算更加准确。

下面参照图4描述本发明实施例的图像畸变校正装置400，如图4的结构示意图所示，图像畸变校正装置400包括：提取单元402，提取源图像的图像轮廓；获取单元404，根据图像轮廓获取畸变信息；以及校正单元406，采用畸变信息对源图像进行校正，以与目标图像配准。

根据本发明实施例，获取单元404将源图像的图像轮廓与目标图像的图像轮廓进行比较，计算位移系数和伸缩系数，作为畸变信息。

根据本发明实施例，该装置400还可以包括均衡化单元(未示出)，对源图像进行直方图均衡化处理。提取单元402提取均衡化单元处理后的源图像的图像轮廓。

根据本发明实施例，针对源图像图像轮廓的所有像素行，获取单元404逐行地计算位移系数和伸缩系数，以得到整个源图像的像素偏移图，作为畸变信息。

根据本发明实施例，针对源图像图像轮廓的一部分像素行，获取单元404逐行地计算位移系数和伸缩系数。在这种情况下，装置400还包括外推单元408，将获取单元404计算的一部分像素行的位移系数和伸缩系数进行拟合并外推到其他像素行，以得到整个源图像的像素偏移图，作为畸变信息。

下面参照图5和6的示例来说明本发明实施例的图像畸变校正方法100和装置400的校正效果。

图5示出了采用本发明实施例的图像畸变校正方法的校正结果的示例。在该实验示例中，在体实验的对象为健康志愿者，成像参数为：TE/TR＝105/3400msec；视野：230×230mm；成像矩阵：160×160；带宽：895Hz/Pixel；平均次数：3；扩散成像使用的b值分别为：0，500，1000。图5显示了对b＝500的图像和b＝1000的图像进行畸变校正的结果，其中白色虚线表示提取的轮廓线。图5a和图5d是b＝0的无畸变的目标图像，轮廓线与图像边界吻合地很好；图5b和图5c是未经过校正的b＝500和b＝1000的源图像，可以看出，在白色箭头所指的区域，轮廓线与图像边界有比较大的差异；图5e和图5f是采用本发明图像畸变校正方法校正过的b＝500和b＝1000的图像，提取的轮廓线与图像边界吻合地很好，也与b＝0的图像的轮廓线吻合地很好。可见，图5示例说明本发明图像畸变校正方法有效校正了涡流场造成的图像畸变。

图6示出了利用本发明实施例的图像畸变校正方法校正后的图像计算得到的ADC图。图6a是使用未经校正的扩散加权图计算得到的ADC图，白色箭头所指的部分有明显的错误，这是由于b＝500和b＝1000的扩散加权图的畸变造成的。图6b是使用本发明校正方法校正过的扩散加权图计算得到的ADC图，图像质量有了明显改善。图6示例显示了本发明实施例的图像畸变校正方法对ADC图的改善作用。

结合以上具体描述，本发明的图像畸变校正方法和装置基于图像轮廓来获得畸变信息，用于校正，不需要图像有相似的对比度，能够用于多种不同b值，特别是b值较大的情况；此外，不需要序列的改动和额外的参考扫描，不会增加扫描时间。校正后的图像可以用于计算各向同性的DWI图像和ADC图，消除了这些图像中的伪影。可见，本发明的图像畸变校正方法和装置能够方便、快速和准确地校正扩散梯度造成的畸变和伪影，适用于多种类型的扩散加权成像技术。

上述本发明的说明书和实施方式仅仅以示例性的方式对本发明实施例的图像畸变校正方法和装置进行了说明，并不用于限定本发明的范围。对于公开的实施例进行变化和修改都是可能的，其他可行的选择性实施例和对实施例中元件的等同变化可以被本技术领域的普通技术人员所了解。本发明所公开的实施例的其他变化和修改并不超出本发明的精神和保护范围。