一种2D-3D影像配准算法的制作方法

本发明涉及手术导航技术领域，尤其涉及一种3dct影像和2d透视影像的配准算法。

背景技术：

3dct影像和2d透视影像配准算法计算3dct影像和2d透视影像之间的位姿变换值，在影像导航置入手术中被用于将术前信息（基于ct的手术计划数据）整合到术中x射线投影影像中，以提高到达目标位置的精度。当前的算法实现通常使用3dct影像生成2ddrr（数字重建放射）影像，计算drr影像和透视影像之间的相似度，找到ct的某个位姿使生成的drr影像和透视影像具有最佳的相似度。由于目标变量有6个自由度（其中3个表示旋转，另外3个表示平移），传统的基于模板匹配的相似度评价的算法或基于随机采样搜索的算法并不能得到精准的结果，且非常耗时，影响算法在实际场景的应用。

技术实现要素：

发明目的：本发明针对上述不足，提出了一种快速、准确的2d-3d影像配准算法。

技术方案：

一种2d-3d影像配准算法，包括步骤：

（1）在术前扫描患者患处得到术前3d影像；在术中分别对患者相对应患处进行正侧位透视扫描得到患者患处的正侧位2d透视影像；

（2）构建术中拍摄正侧位透视影像设备的空间模型，并根据步骤（1）计算得到术前3d影像在空间模型中的初始位姿；

（3）采用步骤（2）构建的空间模型生成术前3d影像的正侧位drr影像，并将其与其对应的2d透视影像进行模板匹配，得到术前3d影像的粗匹配变换位姿；

（4）使用梯度下降优化算法进行精准优化匹配，得到术前3d影像的最终变换位姿。

所述步骤（3）中，位姿包括三维欧拉角和三维平移向量；具体计算得到术前3d影像的粗匹配变换位姿如下：

通过模板匹配计算生成的drr影像和其对应的2d透视影像的归一化相关匹配值，得到最佳的匹配位置，据此计算得到术前3d影像粗匹配后的平移向量；

对于初始欧拉角，分别以生成的drr影像的每个分量角度以预设的偏差取值，并分别计算其与对应的2d透视影像之间的归一化相关匹配值，得到最高匹配值对应的欧拉角；

综上计算得到术前3d影像的粗匹配变换位姿。

所述步骤（4）中：

以术前3d影像的粗匹配变换位姿作为待优化变量；

分别计算正侧位drr影像的梯度图及经过粗匹配变换位姿变换后的3d影像生成的对应drr影像的梯度图，并以二者之间的相似度值作为评价函数；

使用梯度下降优化算法计算术前3d影像的最终变换位姿。

使用scharr算子分别计算所述drr影像的梯度图和术中2d透视影像的梯度图。

所述drr影像的梯度图和2d透视影像的梯度图之间的相似度值metric计算如下：

正位drr影像和正位2d透视影像的单椎骨梯度相似度值计算方法为：

其中，表示正位2d透视影像的梯度图，表示正位drr影像的梯度图，是正位透视影像的方差值；侧位相似度的计算方法类似，最后计算总的相似度值。

使用梯度下降优化算法计算术前3d影像的最终变换位姿具体为：

（41）计算当前值t的梯度，梯度由每个分量t[i]的导数组成；

（42）设定偏差δ，则t[i]的导数为：

（43）沿梯度方向优化当前值t，计算优化后的变量t[i]'为：

其中为梯度因子，gradientmagnitude为梯度dt的模长，steplength为优化步长；

（44）计算震荡系数shakefactor：

dpt[i]为上次迭代计算得到的pt每一分量的导数，初始值设为0；若shakefactor<0，则优化步长steplength调整为之前的一半，并重新确定偏差，重复步骤（42）；否则转至步骤（45）；

（45）判断优化步长steplength是否满足终止条件，若满足则停止迭代，输出最终的t；若不满足，则重新确定偏差δ，返回步骤（42）。

所述步骤（45）中的终止条件为steplength<0.03。

有益效果：本发明的算法首先采用模板匹配在一定范围内进行粗匹配，然后使用梯度下降优化算法得到最优结果，可以在保证配准速度的同时，得到非常精准的配准结果。

附图说明

图1为本发明的2d-3d影像配准算法框架流程图。

图2为术前ct影像中模拟x光机生成drr影像的原理图。

图3为在图像x和y方向上进行模板匹配计算的归一化相关匹配值。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明。

本发明主要应用在人体椎骨的ct透视影像配准中，图1为本发明的流程图，如图1所示，本发明的2d-3d影像配准算法包括步骤：

（1）通过c臂机或ct机在术前扫描患者患处（本实施例为脊椎）得到术前3d影像；对术前3d影像进行自动分割得到椎骨图像，进而得到椎骨的中心点位置（即椎骨中心点在影像坐标系中的坐标）和包围盒（即椎骨的外接长方体）；

（2）通过c臂机或ct机在术中分别对患者患处进行正侧位透视扫描，得到患者患处的正侧位术中2d透视影像，并使用透视变换矩阵得到单个椎骨影像，在术中2d透视影像中的椎骨中心点加上标记；

（3）构建术中拍摄正侧位透视影像的c臂机或ct机的空间模型（以x射线发射管中心作为空间模型光心，c臂机成像平面作为空间模型成像平面），根据构建的空间模型中成像平面位姿及步骤（2）标记的术中2d透视影像的椎骨中心点的坐标可以计算出术中2d透视影像中的椎骨中心点在构建的空间模型中的位置坐标；

（4）根据步骤（1）得到的椎骨中心点在影像坐标系中的坐标，可以计算出3d影像在c臂机或ct机空间模型中的初始位姿t0[r0,t0]，其中r0是表示旋转的3维欧拉角，t0是表示平移的3维平移向量；

（5）采用步骤（3）构建的c臂机或ct机空间模型生成术前3d影像正侧位2ddrr影像，生成方法如图2所示：

从x射线源向2ddrr成像平面每个像素的方向采样，由以下公式计算最终的像素值ixy：

式中i是术前3d影像的采样值，作为衰减系数。生成算法需要大量计算，并且是可以并行计算的，所以使用cuda加速计算；

（6）利用模板匹配进行粗配准：

通过模板匹配计算生成的2ddrr影像和其对应的术中2d透视影像的归一化相关匹配值，结果如图3所示，最亮的像素点对应最佳的匹配位置，这样就可以得到术前3d影像粗匹配后的平移向t1；对于初始欧拉角r1，分别以生成的2ddrr影像的每个分量角度以预设的偏差取值，并分别计算其与对应的术中2d透视影像之间的归一化相关匹配值，得到最高匹配值对应的欧拉角r1，最终得到粗匹配后的术前3d影像的变换位姿t1[r1,t1]；

（7）粗匹配可以将平移向量的误差控制在0.5毫米以内（透视影像一个像素对应的长度），下面需要使用梯度下降优化算法进行精准优化匹配，具体优化步骤为：

（71）将粗匹配后的变换位姿t1[r1,t1]传递给梯度下降优化算法的初始值t，t就是优化算法的待优化变量；

（72）由当前的变换位姿t得到变换后的3d影像，并采用步骤（5）生成drr影像，并使用透视变换矩阵得到单个椎骨影像（脊柱是非刚体，所以算法进行单椎骨配准），使用scharr算子分别计算该drr影像的梯度图和透视影像的梯度图，再计算两幅梯度图的相似度值metric；

以正位透视为例，正位drr影像和正位2d透视影像的单椎骨梯度图相似度值计算方法为：

其中，表示正位2d透视影像的梯度图，表示正位drr影像的梯度图，是正位透视影像的方差值。该方法计算每个对应像素的相似度，然后将所有像素的相似度值累加。显然，透视图和drr图越相似，越接近0，像素相似度值越接近1，的值越接近梯度的像素数，它们是正相关的；侧位相似度的计算方法类似，最后计算总的相似度值。

本算法框架在计算drr、梯度图和相似度值都使用cuda完成，这样不仅可以利用gpu加速计算，还可以减少数据在主机内存和设备内存（gpu显存）间传递所花费的时间，最终只需要将相似度值从显存中拷贝至主机内存。

（73）计算当前6自由度t的梯度，梯度由每个分量t[i]的导数组成；具体为：定义偏差δ，则t[i]的导数为：

（74）调整优化步长steplength，具体为计算震荡系数shakefactor：

dpt[i]为上次迭代计算得到的t（定义为pt）每一分量的导数，初始值设为0。如果shakefactor<0，则steplength调整为之前的一半，并重新确定偏差，重复步骤（73）；否则转至步骤（75）；

（75）沿梯度方向优化当前值t，计算优化后的变量t[i]'为：

其中为梯度因子，gradientmagnitude为梯度dt的模长，steplength为优化步长；

（76）判断优化步长steplength是否满足终止条件（即steplength<0.03），若满足则停止迭代，输出最终的变换位姿t；若不满足，则重新确定偏差δ，返回步骤（73）。

以上详细描述了本发明的优选实施方式，但是本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换（如数量、形状、位置等），这些等同变换均属于本发明的保护范围。