CN103020944A - 一种磁共振图像的二维几何校正的方法 - Google Patents

Abstract

本发明公开了一种磁共振图像的二维几何校正的方法，包括：形成水模，并记录水模的网格点的数量与位置信息，作为初始控制点的数量与位置信息；利用水模在磁共振成像的不同区域采集多幅水模图像，提取每幅水模图像中的控制点的数量与位置信息；水模图像的控制点对应于水模的网格点；根据磁共振成像的不同区域组合多幅水模图像，计算网格点的位置信息与控制点的位置信息之间的形变量，并根据形变量估算形变场；根据形变场，修正磁共振图像。本发明具有可以使用较小尺寸的水模进行较大尺寸的成像区域进行无参数类校正的优势，改进了磁共振图像校正的标定过程。

Description

一种磁共振图像的二维几何校正的方法

技术领域

本发明涉及磁共振成像领域，尤其涉及一种磁共振图像的二维几何校正的方法。

背景技术

几何校正是一种把在几何形状上发生畸变的图像进行校正的一种数字图像处理技术。几何畸变在电视屏幕边缘、相机镜头边缘部分和磁共振图像中广泛存在，本发明特指磁共振图像的几何畸变。磁共振图像的几何畸变主要由于磁共振成像仪器硬件上产生梯度场非线性导致的，在越远离磁体中心的区域畸变越严重。磁共振图像的几何校正技术对医学诊断具有很大的帮助。

解决磁共振图像的几何畸变问题的通常方法是利用水模进行成像实验，采用合适的插值算法估算控制点的形变场，进而利用数字图像处理的几何变换进行几何校正。其中比较有代表性的就是王德民提出的王德明的算法，这种算法根据水模控制点的位置信息可以预知的特点，利用控制点真实位置为确定图像中控制点的位置初步定位，在这个位置附近精确寻找控制点在图像中的位置。但是，当图像的视野(FOV)较大时，图像形变也较大，控制点可能从所在的行变形到另外一行，这时上述算法将不能有效地提取相应的控制点，有可能会将另外一行的控制点作为当前行的控制点，这是第一个缺陷，在图6当中有呈现，图中对于形变量小的控制点A、B用该方法能够进行初步定位，对于形变量大的点C、E、G或D、F、H，进行初步定位时，得到的点是相邻行的点或者相邻列的点；这种几何校正方法采用单纯的插值算法不能估算控制点覆盖范围以外的形变场，这是第二个缺陷。

发明内容

本发明针对上述现有技术中无法有效提取控制点和无法估算控制点覆盖范围以外的形变场等缺陷，提出了一种磁共振图像的二维的几何校正的方法，包括：

步骤一：设计水模，采集实验数据并记录所述水模的网格点的数量与位置信息，作为初始控制点的数量与位置信息；

步骤二：利用所述水模在磁共振成像的不同区域采集多幅水模图像，提取每幅所述水模图像中的控制点的数量与位置信息，所述控制点的位置信息是所述水模中的控制点在图像形变后的位置；所述水模图像的控制点对应于所述水模的网格点；

控制点的提取方式主要结合局部最大值迭代搜索法和排序法，通过局部最大值对控制点进行定位，用迭代地方式修改参数以确保搜索到所有的控制点，并且将搜索到的控制点进行排序，将视觉上的控制点的位置和真实的位置对应起来，这是本发明的第一个创新点，在具体的实施方式中有具体的描述。

步骤三：根据所述磁共振成像的不同区域组合所述多幅水模图像，计算所述水模网格点即所述控制点的位置信息及形变量，并根据所述形变量估算形变场；

步骤四：根据所述形变场，修正存在几何形变的磁共振图像。

其中，所述步骤二中提取水模图像中控制点的数量及位置信息的步骤包括：

步骤A1：采用十字形模板对所述水模图像进行卷积处理，用于突出网格交叉点中心的信号，便于后续步骤确定准确的控制点位置；

步骤A2：设置阈值滤除所述水模图像的背景，将小于所述阈值的背景的灰度值设置为0；

步骤A3：遍历所述水模图像内一定尺寸的矩形窗范围内的亮度最大值的像素点为所述控制点，记录所述控制点的数量与位置信息；

步骤A4：若所述控制点的数量大于所述网格点的数量，则提高阈值并重新执行步骤A2；若所述控制点的数量小于所述网格点的数量，则降低阈值并重新执行步骤A2；若所述控制点的数量等于所述网格点的数量，则完成提取控制点的数量与位置信息。

其中，所述步骤三中，对控制点范围以内的各所述控制点的形变量采用双线性插值的方法估算每一像素点的形变量，形成形变场；对于所述控制点范围以外的各点采用最小二乘法进行拟合估算每一像素点的形变量，形成形变场。

其中，所述步骤四中修正磁共振图像的步骤是利用所述存在几何形变的磁共振图像计算最终图像(即校正后的图像)上的每个像素点强度的过程，具体包括：

步骤C1：对于所述最终图像上的一个像素点，依据所述形变场计算所述其对应于所述存在几何形变的磁共振图像中的位置；

步骤C2：将所述最终磁共振图像上像素点的灰度值赋值为所述存在几何形变的磁共振图像中形变位置的像素点的灰度值；

步骤C3：重复执行所述步骤C1至步骤C2，直至所述最终图像上所有像素点都被处理后终止。

其中，当所述像素点真实位置超过所述磁共振图像时，所述像素点的灰度值赋值为0。

其中，当所述像素点真实位置在图像范围内，但同时又超出了形变场描述的范围，则在所述像素点的邻近区域选取至少两个形变场数据点，采用最小二乘法线性拟合得到所述像素点的形变位置。这是本发明的第二个创新点。

其中，所述最小二乘法的公式如下表示：

f(x，y)＝bx+y+a

$a=\frac{\mathrm{\Σ}{x}_i^2\mathrm{\Σ}{y}_i-\mathrm{\Σ}{x}_i\mathrm{\Σ}{x}_i{y}_i}{\mathrm{n\Σ}{x}_i^2-{\left(\mathrm{\Σ}{x}_i\right)}^2}$

$b=\frac{-\mathrm{\Σ}{x}_i\mathrm{\Σ}{y}_i+n\mathrm{\Σ}{x}_i{y}_i}{\mathrm{n\Σ}{x}_i^2-{\left(\mathrm{\Σ}{x}_i\right)}^2}$

式中，x表示形变位置的横坐标，y表示形变位置的纵坐标，xi，yi表示拟合的数值。

本发明在控制点的提取上采用局部最大值的迭代搜索方法，克服了常规方法难以提取形变量较大的控制点的缺陷；本发明中对于控制点范围以外的区域采用最小二乘法线性拟合的方法估算形变场，有效克服了常规方法中区域条件限制的不足，取得较好的校正效果。

附图说明

图1是本发明二维几何校正方法的流程图。

图2是水模数据采集过程中采用的定位板的示意图。

图3是左下方水模图像的示意图。

图4是水模图像经过十字形模板卷积后的结果图。

图5是水模图像经过阈值过滤后的结果图。

图6是提取的控制点的示意图。

图7是本实施例中处理水模图像的流程图。

图8是控制点在原水模图像中相对显示的示意图。

图9是在控制点范围内的点采用磁共振图像自动的二维几何校正的方法校正之后的图像。

图10是在控制点范围内的点采用磁共振图像自动的二维几何校正的方法校正，控制点范围以外的点采用控制点边界的形变量等值外延的校正结果图。

图11是采用本发明磁共振图像的二维几何校正的方法校正之后的结果图。

图12是从西门子3.0T仪器当中采集得到的校正的结果图。

图13是未经几何校正的真实MRI人体腹部原图。

图14是经过根据本发明的实施例的几何校正方法进行几何校正后的真实MRI人体腹部图。

图15是图12中从西门子仪器当中采集得到的校正的结果对应的校正之后的图像。

具体实施方式

结合以下具体实施例和附图，对本发明作进一步的详细说明。实施本发明的过程、条件、实验方法等，除以下专门提及的内容之外，均为本领域的普遍知识和公知常识，本发明没有特别限制内容。

本发明的磁共振图像的二维几何校正的方法，如图1所示，包括：

步骤一：形成水模，并记录水模的网格点的数量与位置信息，作为初始控制点的数量与位置信息，这里选择网格点作为控制点。

步骤二：利用水模根据磁共振成像的不同区域采集多幅水模图像，提取每幅水模图像中的控制点的数量与位置信息；水模图像的控制点对应于水模的网格点。

步骤三：根据磁共振成像的不同区域组合多幅水模图像，计算网格点的位置信息与控制点的位置信息之间的形变量，并根据形变量估算形变场。

步骤四：根据形变场，修正磁共振图像。

本实施例中，设计的水模的内尺寸为190mm*190mm*190mm，外尺寸为210mm*210mm*210mm，水溶液填充水模形成18行18列，选择网格点为控制点，经过将水模图像拼接在一起之后，网格点的个数为324个。

本实施例中，优选地以主磁场的方向为纵深方向，人体左右方向为x轴方向，人体前后方向为y轴方向，如图2所示。实验分为四次进行，水模摆放以此摆放在仪器的床面与磁体中心定位线的交叉点所在的位置O的左上、左下、右上、右下四个位置，并确保能无间隙地拼接在一起。四次水模实验均在西门子3.0T磁共振扫描仪上采用相同的成像协议完成。从四次实验数据中分别选择y＝0层的四幅水模图像，图3显示的是位于左下方位的水模图像。四幅水模图像的324*4个控制点覆盖了以磁体中心为中心的y＝0平面上380mm*380mm的范围。

随后通过以下步骤对取得的水模图像进行处理，如图7所示，包括：

步骤A1：采用3行9列和9行3列的十字形模板的归一化的均值模板对水模图像进行卷积，加重网格点上的灰度值，由图3卷积后的图像如图4所示。

步骤A2：设置阈值T滤除水模图像的背景，将小于阈值T的背景的灰度值设为0，滤除后的图像如图5所示。

步骤A3：遍历水模图像上在一定尺寸的矩形窗范围内的亮度最大值的像素点为控制点，如图6所示，记录控制点的数量与位置信息。矩形范围为窗口W*Wmm的矩形范围，本实施例中控制点之间的间隔在理论上为10mm*10mm，因此设置W＝9。从而每个矩形范围内只会存在一个控制点，求取该矩形范围内亮度值最高的像素点为控制点。对控制点进行从上到下，从左到右排序，保存控制点的位置。

步骤A4：若控制点的数量大于网格点的数量，则提高阈值T并重新执行步骤A2；若控制点的数量小于网格点的数量，则降低阈值并重新执行步骤A2；若控制点的数量等于，则完成提取控制点的数量与位置信息。

图8显示的是控制点在原水模图像中相对显示的示意图，其中黑色方框点为控制点所在的位置。

本发明中估算形变场的步骤包括步骤B：对控制点覆盖范围以内的各点结合形变量采用双线型插值的方法估算形变场；对于控制点范围以外的各点采用最小二乘法进行拟合估算形变场。其中，控制点范围是最外围的控制点所围成的区域。

本发明中步骤四中修正磁共振图像的步骤是利用存在几何形变的磁共振图像计算最终图像(即校正后的图像)上的每个像素点强度的过程，具体包括：

步骤C1：对于最终图像上的一个像素点，依据形变场计算其对应于存在几何形变的磁共振图像中的位置；

步骤C2：将最终磁共振图像上像素点的灰度值赋值为存在几何形变的磁共振图像中形变位置的像素点的灰度值；

步骤C3：重复执行步骤C1至步骤C2，直至最终图像上所有像素点都被处理后终止。

其中，当像素点形变位置超过磁共振图像时，像素点的灰度值赋值为0。

其中，当像素点的真实位置超出了形变场的范围，则在像素点的邻近区域选取至少两个形变场的数据点。本实施例中选用10个形变场的数据点，采用最小二乘法线性拟合得到像素点的形变位置。图10显示的是在控制点范围内的点采用磁共振图像自动的二维几何校正的算法校正，控制点范围以外的点采用控制点边界的形变量等值外延的校正结果图。

经过以上步骤得到的结果图如图11所示。图11显示的是采用本发明磁共振图像的二维几何校正的方法校正之后的结果图。与图12从西门子仪器当中采集得到的校正的结果图相比，在偏离磁体中心较远的地方具有更好的校正效果，几何形变本身会引起灰度不均，需要对图像进行灰度校正以及图像加强等操作，当前的几何校正算法已经将灰度校正引入在几何校正当中，本发明的主旨在于几何校正算法的改进，对于灰度校正并为未进行详细探讨。

实施例：

本实施例将水模摆在如图2所示的定位板的左上位置，水模的外轮廊与十字标位线贴合，采用的扫描序列为tr_mpr_tra_p2_iso，体素为1.0mm*1.0mm*1.0mm，FOV为：256mm*256*mm，FOV中心偏离磁体中心R101.7mm，P18.8mm，H101.0mm，采集的层数为144，得到数据集Data1。

将水模摆在如图2所示的定位板的右上位置，水模的外轮廊与十字标位线贴合，采用的扫描序列为tr_mpr_tra_p2_iso，体素为1.0mm*1.0mm*1.0mm，FOV为：256mm*256*mm，FOV中心偏离磁体中心L120.7mm，P18.8mm，H109.0mm，采集的层数为144，得到数据集Data2。

将水模摆在如图2所示的定位板的右下位置，水模的外轮廊与十字标位线贴合，采用的扫描序列为tr_mpr_tra_p2_iso，体素为1.0mm*1.0mm*1.0mm，FOV为：256mm*256*mm，FOV中心偏离磁体中心L120.7mm，P18.8mm，F118.0mm，采集的层数为144，得到数据集Data3。

将水模摆在如图2所示的定位板图的左下位置，水模的外轮廊与十字标位线贴合，采用的扫描序列为tr_mpr_tra_p2_iso，体素为1.0mm*1.0mm*1.0mm，FOV为：256mm*256*mm，FOV中心偏离磁体中心R123.7mm，A12.3mm，F123.7mm，采集的层数为144，得到数据集Data4。

根据图像自带的参数信息Slice location(层位置)选择Data1当中y＝0所在的数据层，按照本发明步骤二提取324个控制点，根据参数pixel spacing(像素间距)和Imageorientation(图像方位)获得324个控制点形变之后的位置P2。以磁体中心为原点，控制点的准确位置P1是可以测量得到。左上部分324个控制点的形变量可以根据公式D＝P1-P2计算得到。

根据图像自带的参数信息Slice location选择Data2当中y＝0所在的数据层，按照本发明步骤二提取324个控制点，根据参数pixel spacing和Image orientation获得324个控制点形变之后的位置P2。以磁体中心为原点，控制点的准确位置P1是可以测量得到。右上部分324个控制点的形变量可以根据公式D＝P1-P2计算得到。

根据图像自带的参数信息Slice location(层位置)选择Data3当中y＝0所在的数据层，按照本发明步骤二提取324个控制点，根据参数pixel spacing(像素间距)和Imageorientation(图像方位)获得324个控制点形变之后的位置P2。以磁体中心为原点，控制点的准确位置P1是可以测量得到。右下部分324个控制点的形变量可以根据公式D＝P1-P2计算得到。

根据图像自带的参数信息Slice location选择Data4当中y＝0所在的数据层，按照本发明步骤二提取324控制点，根据参数pixel spacing和Image orientation获得324个控制点形变之后的位置P2。以磁体中心为原点，控制点的准确位置P1是可以测量得到。左下324个控制点的形变量可以根据公式D＝P1-P2计算得到。

将以上获得的四个水模图像的形变量存放在一个矩阵地图(map)当中，存放的具体方式为：假设矩阵中两个数据之间为1mm两个同一幅图像中的两个相邻控制点之间的间隔为9个数据点，相邻两幅图像如左上和右上，左上水模图提取的最右边的控制点距离右上水模图提取的最左边控制点的间隔为40mm，间隔39个数据点。控制点之间的数据点的值利用邻近的数据点采用双线性插值的方式进行填充。公式如下所示

$I(x,y)=\frac{I\left(Q_{11}\right)}{(x_2-x_1)(y_2-y_1)}(x_2-x)(y_2-y)$

$+\frac{I\left(Q_{21}\right)}{(x_2-x_1)(y_2-y_1)}(x-x_1)(y_2-y)$

$+\frac{I\left(Q_{12}\right)}{(x_2-x_1)(y_2-y_1)}(x_2-x)(y-y_1)$

$+\frac{I\left(Q_{22}\right)}{\left(x_2{-x}_1\right)(y_2-y_1)}(x-x_1)(y-y_1)$

其中，Q₁₁(x₁，y₁)、Q₂₁(x₂，y₁)、Q₁₂(x₁，y₂)、Q₂₂(x₂，y₂)为位置处于矩形框的四个角点，I(x，y)指处于(x，y)位置的灰度值，x、y为矩形内部点的横纵坐标。例如I(Q₁₁(1，2))＝1、I(Q₂₁(2，2))＝2、I(Q₁₂(1，1))＝3、I(Q₂₂(2，1))＝4，利于公式1.1计算可以得到

$I\left(\mathrm{1.5,1.5}\right)=\frac{1}{(2-1)(2-1)}(2-1.5)(2-1.5)$

$+\frac{2}{(2-1)(2-1)}(1.5-1)(2-1.5)$

$+\frac{3}{(2-1)(2-1)}(2-1)(1.5-1)$

$+\frac{4}{(2-1)(2-1)}(1.5-1)(1.5-1)$

$=2.5$

根据一幅未经几何校正MRI原图自带的参数信息Pixel Spacing和Image Orientation获得图像某一点的位置信息，图像如图13所示。此时获得的位置信息作为真实位置，根据位置信息从形变场中获得相应的形变量，计算得到这一点形变之后的位置，取形变之后这一点的灰度值作为当前真实位置的灰度值，遍历该图像所有的点，即可得到校正之后的图像。如果形变之后的位置超出原图，真实位置的该点的灰度值用0填充。当计算得到的真实位置超出map的边界时，取邻近的十个点作为参考点，采用最小二乘法线性拟合得到该真实位置形变之后的位置。采用最小二乘法拟合公式为

f(x，y)＝bx+y+a

$a=\frac{\mathrm{\Σ}{x}_i^2\mathrm{\Σ}{y}_i-\mathrm{\Σ}{x}_i\mathrm{\Σ}{x}_i{y}_i}{\mathrm{n\Σ}{x}_i^2-{\left(\mathrm{\Σ}{x}_i\right)}^2}$

$b=\frac{-\mathrm{\Σ}{x}_i\mathrm{\Σ}{y}_i+n\mathrm{\Σ}{x}_i{y}_i}{\mathrm{n\Σ}{x}_i^2-{\left(\mathrm{\Σ}{x}_i\right)}^2}$

其中xi，yi表示拟合的数值，拟合得到a、b的数值以后，根据公式可以根据坐标计算形变场的数值。参考点的选取方法包括：当真实位置位于map中心位置的上下方时，取真实位置所在的一列的十个点作为参考点；当真实位置位于map中心位置的左右方时，取真实位置所在的一行的十个点作为参考点；当真实位置位于map的中心位置左上、右上、右下、左下方时，取以真实位置和最近的矩阵对角点所在的直线上邻近的十个点作为参考点。当选取的参考点不在网格点上时，采用灰度级补插的方式计算得到。

经过校正后的图片如图14，与图15在西门子仪器上采集得到的图12对应的经几何校正的真实MRI人体腹部原图相比，在偏离磁体中心较远的地方具有更好的校正效果，几何形变本身会引起灰度不均，需要对图像进行灰度校正以及图像加强等操作，当前的几何校正算法已经将灰度校正引入在几何校正当中，本发明的主旨在于几何校正算法的改进，对于灰度校正并未进行详细探讨。

本发明的保护内容不局限于以上实施例。在不背离发明构思的精神和范围下，本领域技术人员能够想到的变化和优点都被包括在本发明中，并且以所附的权利要求书为保护范围。

Claims (8)

1.一种磁共振图像的二维几何校正的方法，其特征在于，包括：

步骤一：设计水模，采集实验数据并记录所述水模的网格点的数量与位置信息，作为初始控制点的数量与位置信息；

步骤二：利用所述水模在磁共振成像的不同区域采集多幅水模图像，提取每幅所述水模图像中的控制点的数量与位置信息，所述控制点的位置信息是所述水模中的控制点在图像形变后的位置；所述水模图像的控制点对应于所述水模的网格点；

步骤三：根据所述磁共振成像的不同区域组合所述多幅水模图像，计算所述水模网格点即所述控制点的位置信息及形变量，并根据所述形变量估算形变场；

步骤四：根据所述形变场，修正存在几何形变的磁共振图像。

2.如权利要求1所述的磁共振图像的二维几何校正的方法，其特征在于，所述水模的内尺寸的各边长为190mm，外尺寸的各边为210mm，外壁的厚度为10mm，所述水模的网格点的数量为81个，所述水模的形状为正立方体，所述外壁为玻璃，所述水模的网格为网格状玻璃填充物，所述水模中的填充物为水或者溶液。

3.如权利要求1所述的磁共振图像的二维几何校正的方法，其特征在于，所述步骤二中提取控制点的数量与位置信息的步骤包括：

步骤A1：采用十字形模板对所述水模图像进行卷积处理；

步骤A2：设置阈值滤除所述水模图像的背景；

步骤A3：遍历所述水模图像内一定尺寸的矩形窗范围内的亮度最大值的像素点为所述控制点，记录所述控制点的数量与位置信息；

步骤A4：若所述控制点的数量大于所述网格点的数量，则提高所述阈值并重新执行所述步骤A2；若所述控制点的数量小于所述网格点的数量，则降低所述阈值并重新执行所述步骤A2；若所述控制点的数量等于所述网格点的数量，则完成所述提取控制点的数量与位置信息。

4.如权利要求1所述的磁共振图像的二维几何校正的方法，其特征在于，所述估算形变场的步骤包括：

步骤B：对控制点范围以内的各所述控制点的形变量采用双线型插值的方法估算每一像素点的形变量，形成形变场；对于所述控制点范围以外的各点采用最小二乘法进行拟合估算每一像素点的形变量，形成形变场。

5.如权利要求1所述的磁共振图像的二维几何校正的方法，其特征在于，所述步骤四中修正磁共振图像的步骤包括：

步骤C1：对于所述最终图像上的一个像素点，依据所述形变场计算所述其对应于所述存在几何形变的磁共振图像中的位置；

步骤C2：将所述最终磁共振图像上像素点的灰度值赋值为所述存在几何形变的磁共振图像中形变位置的像素点的灰度值；

步骤C3：重复执行所述步骤C1至步骤C2，直至所述最终图像上所有像素点都被处理后终止。

6.如权利要求5所述的磁共振图像的二维几何校正的方法，其特征在于，当所述像素点形变位置超过所述磁共振图像时，所述像素点的灰度值赋值为0。

7.如权利要求5所述的磁共振图像的二维几何校正的方法，其特征在于，当所述像素点真实位置在图像范围内，但同时又超出了形变场描述的范围，则在所述像素点的邻近区域选取至少两个形变场数据点，采用最小二乘法线性拟合得到所述像素点的形变位置。

8.如权利要求7所述的磁共振图像的二维几何校正的方法，其特征在于，所述最小二乘法的公式如下表示：

f(x，y)＝bx+y+a

$a=\frac{\mathrm{\Σ}{x}_i^2\mathrm{\Σ}{y}_i-\mathrm{\Σ}{x}_i\mathrm{\Σ}{x}_i{y}_i}{\mathrm{n\Σ}{x}_i^2-{\left(\mathrm{\Σ}{x}_i\right)}^2}$

$b=\frac{-\mathrm{\Σ}{x}_i\mathrm{\Σ}{y}_i+n\mathrm{\Σ}{x}_i{y}_i}{\mathrm{n\Σ}{x}_i^2-{\left(\mathrm{\Σ}{x}_i\right)}^2}$

式中，x表示形变位置的横坐标，y表示形变位置的纵坐标，xi，yi表示拟合的数值。

Images