CN106934821A - 一种基于icp算法和b样条的锥形束ct和ct图像配准方法 - Google Patents

Abstract

本发明涉及一种基于ICP算法和B样条的锥形束CT和CT图像配准方法，包括：获取待配准图像和参考图像；对所述待配准图像和参考图像进行图像分割，将所述待配准图像中的目标物分割出来并提取出待配准图像点集数据，将所述参考图像中的目标物分割出来并提取参考图像点集数据；对所述待配准图像点集数据和参考图像点集数据进行ICP仿射配准，根据所述ICP配准结果获取仿射变换阵；根据所述仿射变换阵对所述待配准图像进行仿射配准；对所述仿射配准结果和所述参考图像进行B样条弹性配准。本发明融合了ICP算法快速的优点和B样条弹性配准精确的优点，可以在保证配准精度的前提下，提高配准效率，从而实现快速、准确的图像配准。

Description

一种基于ICP算法和B样条的锥形束CT和CT图像配准方法

技术领域

本发明涉及医学图像、计算机视觉、图像处理等领域，具体为一种基于ICP算法和B样条的锥形束CT和CT图像配准方法。

背景技术

锥形束CT具有成像快、非侵入性、低辐射等优点，并且它能够实时的显示病人软组织和骨骼信息，因此锥形束CT成像技术可用于实时手术导航中。然而在骨科手术导航中，由于CT成像对骨骼的清晰度较高，所以术中锥形束CT图像常常要与术前CT图像进行图像配准以进行精准治疗。但是骨科患者的病变部位可能产生变形，从而导致锥形束CT和CT进行配准时的精度成为难点，并且配准的速度也是手术导航中弹性配准的技术重点。

目前锥形束CT和CT图像的配准大致分为两种：基于特征点的配准和基于体素的配准。这两种配准方法的优势有所不同，基于特征的配准方法只使用局部特征作为配准要素，所以其配准的速度较快，但由于舍弃一部分图像信息，因此其精确度有所欠缺。基于体素的配准方法由于使用了图像的全部像素点信息，所以其配准准确度较高。但是计算量庞大导致基于体素的配准的速度较慢。近年来，将两种方法结合的混合算法相继出现，但仍未应用到锥形束CT和CT图像配准中。

发明内容

针对以上问题，本发明提供了一种基于ICP算法和B样条的锥形束CT和CT图像配准方法，这种方法结合基于点集的配准(ICP算法)和基于体素的配准(B样条配准)，将基于点集的配准结果作为基于体素配准的初值可以快速进行配准，同时基于体素的配准也使本发明的方法具有很高的精度，因此本发明可以快速准确的配准图像，满足其在实时骨科手术导航中的应用。

为实现上述目的，本发明提供如下技术方案：一种基于ICP和B样条的锥形束CT和CT图像配准方法，步骤如下：

为实现上述目的，本发明提供如下技术方案：一种基于ICP和B样条的锥形束CT和CT图像配准方法，步骤如下：

步骤1：获取待配准图像(CT图像)I₁和参考图像(锥形束CT图像)I₂；

步骤2：对所述待配准图像和参考图像进行图像分割，将所述待配准图像中的目标物分割出来并提取其坐标点，然后对坐标点进行采样生成点集数据P＝{p₁,p₂,…,p_t}，其中p₁,p₂,…,p_t表示目标物的坐标点，t表示点集P中坐标点个数，将所述参考图像中的目标物分割出来并提取其坐标点，然后对坐标点进行采样生成点集数据Q＝{q₁,q₂,…,q_t},其中q₁,q₂,…,q_t表示目标物的坐标点，点集Q中坐标点个数同点集P中的个数一样；

步骤3：对所述待配准图像点集数据P和参考图像点集数据Q进行基于点集的配准即ICP仿射配准，其步骤主要为：首先对所述仿射变换阵M_affine进行初始化；利用仿射变换阵M_affine对所述待配准点集P进行变换得到点集P，即P，＝P*M_affine；然后使用公式计算所述参考图像点集Q和P，之间的均方误差，其中p_i和q_i分别表示第i个坐标点；使用LBFGS算法不断优化M_affine对均方误差求最小化，得到最优解M_affine；

步骤4：根据所述仿射变换阵M_affine对所述待配准图像进行仿射配准得到ICP配准图像I_icp；

步骤5：将获取的ICP配准图像I_icp作为基于体素的配准初值，即对所述仿射配准结果I_icp和所述参考图像I₂进行B样条配准得到最终结果I_final。

优选的，所述步骤2中的对所述待配准图像和参考图像进行图像分割的步骤，包括：

对所述待配准图像和参考图像的像素灰度进行归一化为[0,255]；

对所述归一化处理后的待配准图像和参考图像进行阈值分割；

对所述阈值分割后的待配准图像和参考图像进行像素坐标点提取和采样并生成点集数据。

优选的，对所述待配准图像和参考图像的像素灰度进行归一化时采用的公式为：

其中g和g，分别表示归一化后灰度值和原始图像CT值，A和B分别用公式A＝wc-ww/2，B＝wc+ww/2，其中ww表示图像上CT值的窗口范围，wc表示图像上CT值窗口中心。

优选的，对所述归一化处理后的待配准图像和参考图像进行阈值分割的目标物为骨骼区域，其中骨骼的CT值范围为[100,1000],使用归一化公式计算其灰度值。

优选的，所述阈值分割后的待配准图像和参考图像进行像素坐标点提取和采样并生成点集数据时首先提取出256灰度级中骨骼范围的坐标点，然后对提取出的坐标点进行1/50等间距采样得到点集图像。

优选的，所述步骤5中，B样条配准采用多分辨率B样条弹性配准方法，步骤为：使用公式对输入图像进行高斯平滑处理，根据图像大小进行降采样实施多分辨率配准，先将图像在低分辨进行B样条粗配准。首先初始化控制点，然后对两幅图像的相似度测量，通过对相似度值的判断修改控制点，将修改后的控制点代入B样条变换函数对待配准图像进行变换。使用优化方法不断修改控制点对相似度值进行最小化求解，从而得到最终配准结果I_final。

优选的，所述多分辨率配准的分辨率层数分为1到3层，每层中XYZ轴最小分辨率参数分别设为128×128×10,256×256×20,512×512×40。

优选的，所述步骤5中的B样条配准所用变换模型为B样条函数，用公式求取像素点在(x，y，z)位置处的形变场，其中 其中表示向下取整，φ_i,j,k表示x，y，x轴控制点间距分别为δ_i,δ_j,δ_k，大小为n_x×n_y×n_z的网格中序号为ijk的控制点的值，其中i，j，k为了求像素点在(x,y,z)处的位移值T(x,y,z)仅与相邻的4³个控制点有关；B_l,B_m,B_n均为三次B样条函数基函数，B_l表示B样条函数中第l个基函数，其计算公式为B₀(u)＝(1-u)³/6，B₁(u)＝(3u³-6u²+4)/6，B₂(u)＝(-3u³+3u²+3u+1)/6，B₃(u)＝u³/6,其中0≤u≤1B_m和B_n的计算同B_l，这些函数起到权函数的作用，他们根据控制点到(x,y,z)的距离大小来加权每个控制点对T(x,y,z)的贡献。

优选的，所述多分辨率配准的每层分辨率的B样条网格间距分别设为：32×32×16,16×16×8,8×8×4。

优选的，所述相似度测量和优化方法采用灰度均方差作为相似度测量，优化方法选用LBFGS算法。

本发明与现有技术相比的优点在于：本发明基于ICP和B样条的锥形束CT和CT图像配准方法，采用结合基于点的配准和基于体素的配准，在对锥形束CT和CT图像配准时，相对于传统方法具有速度快的同时保持配准精度高的特点，可以很好的应用于基于锥形束CT成像技术的骨科手术导航系统中。

附图说明

图1为基于ICP和B样条的配准方法流程图；

图2为对图像分割和点云提取采样后的点集；其中a是从锥形束CT图像中获取的点集，b是从CT图像中获取的点集；

图3为对点集图像进行ICP配准的效果展示；其中a为ICP配准前的点集图像，b为ICP配准后的点集图像；

图4为本发明配准效果展示图；其中a为初始锥形束CT图像，b为初始CT图像，c为ICP配准结果，d为最终配准结果；

图5为仿真图像生成原理展示图；其中a为生成标准形变场过程示意图，b为生成目标形变场过程示意图；

图6为配准结果的差值对比；其中a为配准前初始锥形束CT和CT图像的差值图，b为ICP配准结果与初始锥形束CT图像的差值图，c为最终配准结果与初始锥形束CT图像的差值图；

图7为配准误差的直方图显示。

具体实施方式

由于锥形束CT和CT图像取自不同成像装置并且拍摄时间不同，因此术前规划系统中的待配准图像即CT图像和术中图像即锥形束CT图像并不相同，例如两者成像视角不同，每层切片的对应关系可能有出入，图像会有偏移和形变，需要进行配准矫正。因此本发明采用如下的基于ICP和B样条的配准方法。

如图1所示，本发明一种基于ICP和B样条的配准方法流程图，包括以下步骤：

步骤S1,获取待配准图像(CT图像)I₁和参考图像(锥形束CT图像)I₂。

步骤S2：对所述待配准图像和参考图像进行图像分割，将所述待配准图像中的目标物分割出来并提取其坐标点，然后对坐标点进行采样生成点集数据P＝{p₁,p₂,…,p_t}，其中p₁,p₂,…,p_t表示目标物的坐标点，t表示点集P中坐标点个数，将所述参考图像中的目标物分割出来并提取其坐标点，然后对坐标点进行采样生成点集数据Q＝{q₁,q₂,…,q_t},其中q₁,q₂,…,q_t表示目标物的坐标点，点集Q中坐标点个数同点集P中的个数一样；

在本实例中，步骤S2中的对待配准图像和参考图像进行图像分割的步骤，包括：

(1)对所述待配准图像和参考图像的像素灰度进行归一化为[0,255]；

(2)对所述归一化处理后的待配准图像和参考图像进行阈值分割；

(3)对所述阈值分割后的待配准图像和参考图像进行像素坐标点提取和采样并生成点集数据。

由于图像是取自不同时间，不同设备或者不同视角的同一场景的图像，因此CT图像和锥形束CT图像的分辨率会有所不同，所以将图像的像素灰度值进行归一化处理。归一化处理所使用的公式为：

其中g和g，分别表示归一化后灰度值和原始图像CT值，A和B分别用公式A＝wc-ww/2，B＝wc+ww/2，其中ww表示图像上CT值的窗口范围，wc表示图像上CT值窗口中心。

在本实例中，阈值分割的目标可以是：肺部、骨骼。其中肺部的CT值范围为[-950,200],骨骼的CT值范围为[100,1000]。优选地，对骨骼进行分割。

在本实例中，提出256灰度级中相应器官范围的坐标点，然后对坐标点进行采样，采样方法的采用如下方式之一：随机采样、等间隔采样。优选地，采用1/50等间隔采样。

如图2所示，为对待配准图像和参考图像进行阈值分割后提取坐标点并进行采样后的实验展示图，其中a是从锥形束CT图像中获取的点集，b是从CT图像中获取的点集。

S3：对所述待配准图像点集数据P和参考图像点集数据Q进行基于点集的配准即ICP仿射配准，其步骤主要为：首先对所述仿射变换阵M_affine进行初始化；利用仿射变换阵M_affine对所述待配准点集P进行变换得到点集P，即P，＝P*M_affine；然后使用公式计算所述参考图像点集Q和P，之间的均方误差，其中p_i和q_i分别表示第i个坐标点；使用LBFGS算法不断优化M_affine对均方误差求最小化，得到最优解M_affine。

在本实例中，ICP配准是基于点的配准，其中所使用的参数优化方法为LBFGS算法，其搜索步长为0.1mm、梯度收敛公差为0.01，最大迭代次数为100次

在本发明实例中，ICP配准的变换模型采用如下方式之一：刚性变换(旋转、平移)，仿射变换(旋转、平移、缩放、剪切)。优选地，采用仿射变换。

如图3所示，为对待配准图像点集和参考图像点集进行ICP配准前(图3中的a)和配准后(图3中的b)的实验对比图。

步骤S4：根据仿射变换阵M_affine对待配准图像进行仿射配准得到ICP配准图像I_icp；

在本实例中，对待配准图像进行仿射配准时的插值方法可采用如下方式之一：线性插值、三次插值。优选地，选用线性插值。

步骤S5：将获取的ICP配准图像I_icp作为基于体素的配准初值，即对仿射配准结果I_icp和参考图像I₂进行多分辨率B样条弹性配准得到最终结果I_final。

将ICP配准图像和参考图像作为B样条弹性配准的初始图像以提高配准精度和配准速度。同时B样条配准前，对ICP配准图像和参考图像进行高斯滤波处理用以提高配准速度。

在本实例中，B样条配准的主要步骤为：对输入图像进行平滑处理，根据图像大小进行降采样实施多分辨率配准，先将图像在低分辨进行B样条粗配准。首先初始化控制点，然后对两幅图像的相似度测量，通过对相似度值的判断修改控制点，将修改后的控制点代入B样条变换函数对待配准图像进行变换。使用优化方法不断修改控制点对相似度值进行最小化求解，从而得到最终配准结果I_final。

在本实例中，B样条弹性配准部分采用多分辨率配准，分辨率根据具体图像大小可分为1到3层，其中每层分辨率的最小分辨率参数分别为128×128×10,256×256×20,512×512×40。例如当图像大小等于或超过512×512×40(如512×512×40)时可以分为三层进行配准，相邻层数之间为倍数关系，(即第一层128×128×10，第二层256×256×20，第三层512×512×40)。首先进行低分辨率配准，然后得出的结果作为高分辨率的输入再进行高分辨率配准直至结束。

在本实例中，每层分辨率的B样条网格间距分别选为：32×32×16,16×16×8,8×8×4。

在本实例中，相似度测量方法可选以下方式之一：互信息、互相关、灰度均方差。优选地，选用灰度均方差作为相似度测量。参数优化选用LBFGS算法。

如图4所示为整个配准实验其中一个切片的结果展示，其中图4中的(a)为参考图像(锥形束CT)，图4中的(b)为待配准图像(CT图像)，图4中的(c)为ICP配准结果，图4中的(d)为最终配准结果。可以看出在经过配准后，待配准图像与参考图像的接近程度逐渐增加。

如表1示，为CT(512*512*53)和锥形束CT(410*410*42)图像配准的结果，其中灰度均方差值作为相似度测量。配准输出参数包括配准前参考图像和待配准图像的灰度均方差值，ICP配准结果和B样条配结果和参考图像的灰度均方差值对比以及配准所需时间。从配准前后的灰度均方差值对比可以看出灰度均方差值提升较多，说明配准后待配准图像和参考图像的一致性有很大提高，此配准方法对于CT和锥形束CT图像配准效果很好。

表1

	配准前	ICP配准结果	B样条配准结果	配准时间
					灰度均方差值	2595	1208	575	630

为对本发明的准确度进行测量，设计仿真实验。如图5为仿真图像生成的流程展示，图5中的(a)将锥形束CT和CT进行配准时产生的形变场作为标准值，图5中的(b)中用标准形变场的逆变换对锥形束CT进行形变生成仿真CT图像，然后对仿真CT图像和锥形束CT图像进行配准生成目标形变场作为目标值。图6为对仿真图像进行ICP配准和B样条配准后的图像差值显示，图6中的(a)为配准前的差值效果展示，图6中的(b)为ICP配准结果和锥形束CT之间像素差值效果展示，图6中的(c)为最终配准结果和锥形束CT之间像素差值效果展示。从图中可以看出，ICP配准后参考图像和待配准图像的空间位置大致重合，但仍有局部差异，在进行B样条配准之后其配准结果和参考图像基本重合。

将标准值和目标值相减求出形变误差，图7画出误差的直方图，从图中可以看出误差大多分布在较低值，通过求均值算出误差值为0.51mm。准确度较高，满足手术导航的精度要求。

表2为比较本实例中配准方法与传统基于体素配准方法的配准速度而设计的实验结果。该实验中传统的配准方法也是使用放射变化和B样条配准作为配准流程，同时相似度测量选用灰度均方差，优化方法选用LBFGS算法。在此基础上，两种方法的收敛条件设为相同。传统方法的配准时间和本发明的配准方法的配准时间如下表所示：

表2

	传统方法	本发明配准方法
			时间(s)CT/锥形束CT	1610	630

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于，所述方法包括：

步骤1：获取待配准图像，即CT图像I₁和参考图像，即锥形束CT图像I₂；

步骤2：对所述待配准图像和参考图像进行图像分割，将所述待配准图像中的目标物分割出来并提取其坐标点，然后对坐标点进行采样生成点集数据P＝{p₁,p₂,…,p_t}，其中p₁,p₂,…,p_t表示目标物的坐标点，t表示点集P中坐标点个数，将所述参考图像中的目标物分割出来并提取其坐标点，然后对坐标点进行采样生成点集数据Q＝{q₁,q₂,…,q_t},其中q₁,q₂,…,q_t表示目标物的坐标点，点集Q中坐标点个数同点集P中个数相同；

步骤3：对所述待配准图像点集数据P和参考图像点集数据Q进行基于点集的配准即ICP仿射配准，其步骤主要为：首先对所述仿射变换阵M_affine进行初始化；利用仿射变换阵M_affine对所述待配准点集P进行变换得到点集P，即P，＝P*M_affine；然后使用公式计算所述参考图像点集Q和P，之间的均方误差，其中p_i和q_i分别表示第i个坐标点；使用LBFGS算法不断优化M_affine对均方误差求最小化，得到最优解M_affine；

步骤4：根据所述仿射变换阵M_affine对所述待配准图像进行仿射配准得到ICP配准图像I_icp；

步骤5：将获取的ICP配准图像I_icp作为基于体素的配准初值，即对所述仿射配准结果I_icp和所述参考图像I₂进行B样条配准得到最终结果I_final。

2.根据权利要求1所述的基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于：所述步骤2中的对所述待配准图像和参考图像进行图像分割的步骤，包括：

对所述待配准图像和参考图像的像素灰度进行归一化为[0,255]；

对所述归一化处理后的待配准图像和参考图像进行阈值分割；

对所述阈值分割后的待配准图像和参考图像进行像素坐标点提取和采样并生成点集数据。

3.根据权利要求2所述的基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于：所述的对所述待配准图像和参考图像的像素灰度进行归一化，采用如下公式：

$g=\left\{\begin{array}{cc}0& -1000<g^{,}<A\\ 255\&times;\frac{g^{,}-A}{B-A}& A<g^{,}<B\\ 255& B<g^{,}<1000\end{array}\right.$

其中g和g，分别表示归一化后灰度值和原始图像CT值，A和B分别用公式A＝wc-ww/2，B＝wc+ww/2，其中ww表示图像上CT值的窗口范围，wc表示图像上CT值窗口中心。

4.根据权利要求2所述的基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于：所述对所述归一化处理后的待配准图像和参考图像进行阈值分割的目标物为骨骼区域，其中骨骼的CT值范围为[100,1000],使用归一化公式计算其灰度值。

5.根据权利要求2所述的基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于：所述阈值分割后的待配准图像和参考图像进行像素坐标点提取和采样并生成点集数的方法为：首先提取出256灰度级中骨骼范围的坐标点，然后对提取出的坐标点进行1/50等间距采样得到点集图像。

6.根据权利要求1所述的基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于：所述步骤5中，B样条配准采用多分辨率B样条弹性配准方法，步骤为：使用公式：

对输入图像进行高斯平滑处理，根据图像大小进行降采样实施多分辨率配准，先将图像在低分辨进行B样条粗配准，首先初始化控制点，然后对两幅图像进行相似度测量，通过对相似度值的判断修改控制点，将修改后的控制点代入B样条变换函数对待配准图像进行变换,使用LBFGS优化方法不断修改控制点对相似度值进行最小化求解，从而得到最终配准结果I_final。

7.根据权利要求6所述基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于：所述B样条变换函数为：采用公式求取像素点在(x，y，z)位置处的形变场，其中：

其中表示向下取整，φ_i,j,k表示x，y，x轴控制点间距分别为δ_i,δ_j,δ_k，大小为n_x×n_y×n_z的网格中序号为ijk的控制点的值,其中i，j，k为了求像素点在(x,y,z)处的位移值T(x,y,z)仅与相邻的4³个控制点有关；B_l,B_m,B_n均为三次B样条函数基函数，B_l表示B样条函数中第l个基函数，其计算公式为B₀(u)＝(1-u)³/6，B₁(u)＝(3u³-6u²+4)/6，B₂(u)＝(-3u³+3u²+3u+1)/6，B₃(u)＝u³/6,其中0≤u≤1B_m和B_n的计算同B_l，起到权函数的作用，根据控制点到(x,y,z)的距离大小来加权每个控制点对T(x,y,z)的贡献。

8.根据权利要求6所述的基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于：所述多分辨率层数分为1到3层，每层XYZ轴最小分辨率参数分别设为128×128×10,256×256×20,512×512×40。

9.根据权利要求8所述的基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于：每层分辨率的B样条网格间距分别设为：32×32×16,16×16×8,8×8×4。

10.根据权利要求6所述的基于ICP算法和B样条的锥形束CT和CT图像配准方法，其特征在于：对两幅图像进行相似度测量时，采用灰度均方差作为相似度测量。