蒲云洁，王学渊

（西南科技大学 信息工程学院，绵阳 621000)

0 引言

随着社会医疗技术的发展，使用计算机影像辅助医生进行外科手术已经逐渐成为趋势，而这其中最重要的一环就是图像的配准。医学图像以模态划分可以分成多种：CT，MRI，PET，X光，内窥镜图像等，以图像维度划分可以分为二维图像和三维图像。二维图像较三维图像缺乏很多有用的空间信息，因此将这两种模态的图像进行信息配准融合，从而定位患者手术过程中的靶区，显得尤为重要。本文研究的内容就是CT断层图像和X光图像之间的配准，这是一种多模态配准模式，也是一种2D/3D配准模式。目前常用的2D/3D配准方法主要分为以下3种：

1）基于图像灰度的配准算法[1]，该方法一般是先用DRR算法模拟X光图像的生成，从患者术前采集的三维体数据中得出多幅DRR图像来和术中采集的X光照片进行相似性判定配准，此算法精度高，但耗费的时间较长，适用于刚体之间的配准。

2）基于图像特征的配准算法[2]，该方法依赖在病人身上提前做好的标记点或者人体内固有的标记点进行定位，计算出相应的空间变换，进而计算出待配准图像中对应标记点的最短特征距离，从而达到配准的目的，此方法的配准速度虽然很快，但是精度与标记点的选取有很大的关联，鲁棒性也比较差。

3）基于以上两类方法结合的配准算法，即既需要计算出待配准图像对应标记点的空间变换和最短特征距离，也需要利用相似性原理进行配准，配准速度快精度也好，但此方法捕捉范围较小，鲁棒性较差。

本文根据以上三种配准方法提出一种基于CUDA加速的多模态配准方法，该方法基于图像灰度，在DRR图像生成上采用光线投射算法，先进行刚体变换粗配准，再采用梯度方向（Gradient orientation，GO）测度进行相似性测度，进而精配准，在参数优化上采用自适应协方差矩阵进化算法（CMA-ES），并且在DRR图像生成和相似性测度这两个耗时最大的过程采用支持CUDA架构的GPU处理器进行加速，最大化的减少配准时间。

1 DRR图像生成

DRR（Digitally Reconstructed Radiograph）数字影像重建技术是2D/3D医学图像配准的关键技术，也是最耗时最影响配准精度的步骤，它最终目的是为了从三维体数据中获取一系列模拟的二维X射线图像。获取DRR图像的算法主要有：抛雪球法（Splatting）、光线投射法（Ray-Casting）、剪切变形法（Shear-Warp）等。本文使用的是光线投射算法。

1.1 光线投射算法

光线投射算法是模拟X射线穿透三维体数据元（本文用的数据是CT体数据），经过衰减和吸收后投影到DRR平面，就形成了一张DRR图像，形成过程如图1所示。

图1 光线投影示意图

光线投射算法生成DRR图像具体过程如下：当虚拟X光源发出射线以一定的步长穿透3D体数据时，沿着其行进方向采集CT值，再把得到的CT值转化成相应的衰减系数μj，累加此方向上的衰减系数，从而计算得到射线到达DRR平面上的电子强度I，最后形成DRR图像。

DRR图像生成过程中，若入射光线电子强度Ii，它与穿过3D体数据后到达DRR平面上的电子强度Io关系如式(1)所示：

式(1)中μ(x)是对应组织的衰减系数，0和d分别是光线经过3D数据的穿入和穿出点。然而由于人体的组织不是均匀的，因而μ(x)是不连续的，所以式(1)可以改写成：

式(2)中μj是光线穿过不同组织的衰减系数，xj是穿过这种组织的有效长度。而衰减系数μ与CT值之间又有如下的转换关系：

式(3)中μw是纯水的衰减系数，通常是定值，HU是组织的CT值。最后将CT值再转化成灰度值就可以完成DRR图像的生成。

1.2 DRR图像和X射线图像之间的刚体变换

对于膝关节图像，无论是DRR图像还是X射线图像，这两者之间的配准可以看成是刚体之间的配准，则两者之间的变换就是刚体变换。刚体变换进一步可以分解为旋转变换和平移变换，从而可以使用六参数矢量来表示，其中θx，θy，θz分别表示物体围绕x，y，z轴的旋转变换，tx，ty，tz表示物体沿着x，y，z轴的平移变换。进一步将旋转变换和平移变换写成矩阵的形式：

则刚体变换矩阵为：

2 相似性测度

图像配准的相似度测度是指通过几何空间变换如刚体变换之后的两幅图像相似程度，常用的相似性测度方法有：互信息[3]（Mutual Information，MI）、绝对误差和[4]（Sum of Absolute Differences，SAD）、归一化互相关[5]（Normalized Cross Correlation，NCC），梯度相关[6]（Gradient Correlation，GC）、模式强度[7]（Pattern Intensity，PI）等。

梯度方向[8]测度已经被应用于特征匹配[9]，物体检测[10]和图像配准[11]，但是传统的梯度方向测度可能会在具有低梯度幅度的区域中受到噪声的影响，从而影响配准的成功率与精度。本文提出一种改进方式：预先对待配准的两幅图像分别设定两个阈值，在计算梯度方向测度时，只计算这两幅图像中大于设定阈值的像素，这两个阈值可以设定为对应图像所有梯度的中值，这样就消除了50%的包含低梯度幅度的像素，减少像素遍历数，从而降低了算法复杂度和减小了配准中噪声的影响。梯度方向测度函数如式(9)、式(10)所示：

式中，N是图像总像素，NL是人为规定的像素总数下界，以补偿像素数低于NL的损失，∇u是梯度算子，I1，I2是待配准的两幅图像，t1，t2是设定的梯度阈值，等于对应图像梯度的中值，θi是∇uI1(i)和∇uI2(i)之间的夹角。最后计算的GO_Value值越接近1配准效果越好。

3 参数优化

参数优化是医学图像配准的重要过程，其实质就是相似性测度目标函数进行不断迭代优化来寻找最佳的变换参数，以达到配准成功的目的。常用的优化算法[3]有模拟退火算法，Powell算法，梯度下降算法，遗传算法等。本文采用的是自适应协方差矩阵进化算法（CMA-ES）。

CMA-ES算法是Nikolaus Hansen等人[12]在2006年提出的一种模拟生物进化过程的算法。该算法的特点是假设不论基因发生何种变化，产生的结果（性状）总遵循这零均值，某一方差的高斯分布。CMA-ES算法实现简单，对目标函数要求宽松，适用于解决高维度复杂的优化问题。在实际操作中，CMA-ES算法需要先初始化参数，设置种群的相关参数：种群中心点（期望）、步长、进化路径、子代数量、协方差矩阵、最大迭代次数以及一些自适应参数。算法步骤[13]如下：

1）种群采样操作：是指在目标函数的解空间中随机选择一个解，并以此为中心生产正态分布种群。以式(10)进行采样并产生新解。

式中k=1，2，...λ，～表示式子两边服从相同分布，表示第k+1代的第n个后代，m(k)表示第k代搜索分布的均值，δ(k)表示第k代的步长，N(0，C(k))表示均值为0，协方差矩阵为C(k)的正态分布。第k代的协方差矩阵为C(k)，C(0)=In*n，I是单位矩阵。协方差矩阵是CMA-ES算法进化策略的关键。

2）选择和重组操作：是指在生成的总群中选择最优子群。当前新的最优子群通过加权重组得到新的期望：

3）更新操作：从一代代的种群中估计协方差矩阵并且更新步长。估计第代协方差矩阵的公式如下所示：

其中，c1≤1是协方差矩阵进化路径Pc的学习率，T表示矩阵转置。hδ表示Heaviside函数，它的作用是控制进化路径Pc的过大增长。通常用平滑指数来构造进化路径：

步长δ更新函数：

式中，e为自然底数，cδ表示步长学习率，dδ为步长更新的阻尼系数。E||N(0，I)||表示进化路径归一化后的期望模长，I为单位矩阵。

CMA-ES算法通过以上三个步骤的循环迭代，直到找到全局最优解。

4 基于CUDA加速的多模态图像配准

传统的图像配准方法通常是基于软件层面的，耗时非常大。因此，采取基于支持CUDA架构的GPU处理器进行硬件加速，弥补传统配准方法在耗时上面的不足，大大提高配准的速度。本文的配准算法流程如图2所示。

图2 本文算法流程图

5 结语

本文实验数据来源为佛山市中医院提供的4名患者膝关节CT断层体数据，单层图像大小为512×512，配套有同一患者的X射线图像。硬件配置为Intel(R) Core(TM) i7-8700K CPU，主机频率为3.7GHz，显卡是NVIDIA公司的GeForce GTX1060 6GB，软件配置为Win10下的Visual Studio 2015+Python 3.6+CUDA 9.0。实验以裁剪后的感兴趣区域X射线图像为参考图像，CT体数据生成的DRR图像为浮动图像，在Python环境下调用VS2015和CUDA生成的动态链接库完成配准，结果如图3～图6所示：

患者1：

图3 患者1配准图

患者2：

图4 患者2配准图

患者3：

图5 患者3配准图

患者4：

图6 患者4配准图

表1数据是每例患者数据均配准50次之后取的平均值，从4组数据的配准结果图来说，也可以看出本文算法对于即使拥有噪声（如患者3的X射线图像）或者有其他干扰物的情况下（如患者4）仍然能够取得很好的配准结果，这证明了本文算法具有高配准率和鲁棒性，利用CUDA进行加速处理后配准速度也有了质的飞跃，为临床应用提供了时效性。

表1 GPU与CPU配准对比