CN104134210A - 基于组合相似性测度的2d-3d医学图像并行配准方法 - Google Patents

Abstract

本发明公开了一种基于组合相似性测度的2D-3D医学图像并行配准方法。该方法首先使用CUDA并行计算模型完成DRR图像的快速生成过程，并组合差值绝对值和SAD与模式强度PI作为新的相似性测度在GPU上进行并行计算，最后将组合相似性测度值传递到CPU上采用基于细菌趋化行为的果蝇优化算法进行优化来寻找最优配准参数。通过实验对本方法性能进行验证表明：由于本发明方法在GPU中实现DRR快速生成及混合相似性测度的计算，有效地提高了本发明方法的执行速度，同时与单一相似性测度相比，本发明采用混合相似性测度提高了配准结果的精确性。

Description

基于组合相似性测度的2D-3D医学图像并行配准方法

技术领域

本发明涉及医学图像快速配准，尤其涉及的是一种基于组合相似性测度的2D-3D医学图像并行配准方法。

背景技术

图像引导放疗(Image-GuidedRadiotherapy,IGRT)是继三维适形放疗(3DCRT)和调强放疗(IMRT)之后的新的放疗技术。在图像引导放射治疗中，为了确定分次治疗的实际位置，物理师利用电子射野影像装置等获取X线射野影像。X线图像反映了治疗期间的实际位置，而治疗前CT体数据代表治疗的理想位置，通过X线射野图像和CT体数据进行2D-3D配准，计算出治疗床的偏移量，最后根据偏移量调整治疗床位置，从而保证放疗过程中摆位的精确性。

在2D-3D图像配准中，需要对三维CT体数据进行重建生成DRR，由于DRR图像的生成过程非常耗时，且在2D-3D图像配准算法中需要从不同角度生成多张DRR图像，故DRR的生成过程直接影响了整个2D-3D配准算法的效率，另外相似性测度的计算对整个配准过程的运行效率和结果精确性也有重要的影响，因此如何提高2D-3D配准算法的运行效率和精确性来满足IGRT需要成为了很多学者研究的热点。GraemeP.Penney等人对六种应用到2D-3D医学图像配准中的基于灰度的相似性测度进行了比较，即使当两幅待配准图像中软组织出现差异的情况下，模式强度和梯度差分两种相似性测度仍旧能取得精确的配准结果，并具有很好的鲁棒性；W.Birkfellner等人提出了一种基于splatrendering的DRR图像生成方法并应用到2D-3D医学图像配准中，在PC机上针对大小为30MB左右的体数据，绘制能用于2D-3D配准的轻度模糊的DRR图像的时间为100ms；Daniel B.Russakoff等人提出了一种基于衰减域的DRR快速生成方法并将其应用到2D-3D配准中，其2D-3D配准方法比使用基于光线投射生成DRR图像的配准算法在执行效率和精确性方面都有一定的优势；OsamaDorgham等人使用八叉树结构来压缩体数据，减少了投射光线和体数据内部空间的交点的数量，从而提高了光线投射生成DRR图像的速度，并将其应用到2D-3D配准中也能达到有效的配准精度，能够满足校正病人摆位误差的精度要求；以上配准算法在采用不同方法改进了DRR图像的生成过程，提高了配准算法的效率，但不能满足图像引导放疗中的实时性要求；kubiasA等人实现了一种基于GPU的2D-3D医学图像配准算法，该算法将DRR生成和相似性测度计算都在GPU上完成，与传统的基于CPU的配准算法相比大大地提高了配准速度，但由于受流水线的限制，GPU的加速能力仍旧没有得到充分的发挥；OsamaM.Dorgham等人在支持CUDA的GPU实现了DRR图像的并行加速生成，同时对光线投射过程进行压缩采样等处理来进一步提高DRR生成速度，最后将该DRR加速生成方法应用到了2D-3D医学图像配准中，并且取得了临床放疗可接受的精度，但是该方法并没有将相似性测度计算过程也搬到GPU上运行。

医学图像配准是将一幅图像的像素空间位置进行变换使其与另一幅图像像素的空间位置对齐的过程，是医学图像融合的前提和基础，在病情诊断、图像引导放疗及手术实施中具有重要的应用价值。在图像引导放射治疗应用中2D-3D医学图像配准是将放射治疗前制定放疗计划用的CT体数据集与放射治疗中代表病人当前治疗摆位的X线图像进行配准，得到相应的配准参数来纠正放疗实施过程中的摆位误差。整个配准过程包括三个主要的步骤：DRR图像的生成、DRR图像与X线图像之间的相似度计算和相似性测度函数寻优。

在2D-3D图像配准过程中，将实际拍摄的X线图像作为参考图像，计算机生成的DRR图像作为浮动图像，通过改变对CT体数据的虚拟透视位置及角度，获取相应的DRR图像，并逐个与真实X线图像进行相似性测度的比较直到出现最优解，此时CT数据的六自由度空间参数即为所求。2D-3D医学图像配准的数学表示如下：

$T^{*}=\underset{T}{\arg \max }S(X,T\left(\mathrm{VolumeData}\right))---\left(1\right)$

式中S是相似性测度函数，T是包含了六个自由度的空间变换矩阵，X是放疗实施过程中代表病人当前位置的X线图像，VolumeData为三维CT体数据，T(VolumeData)即为对三维CT体数据进行矩阵T所示的刚性变换后得到的DRR图像，整个配准过程就是寻求最佳空间变换T使得相似性测度函数S达到最大的过程。2D-3D配准是刚体坐标变换，研究对象的位置和方向变换采用六参数矢量来描述即六自由度。其中，T_x,T_y,T_z分别表示沿着X、Y、Z轴方向的平移量，R_x,R_y,R_z分别表示围绕X、Y、Z轴的旋转量。

发明内容

本发明针对图像引导放疗对2D-3D配准算法的实时性和高精度的要求，且由于DRR生成过程的运行速度严重制约了2D-3D医学图像配准算法的效率，首先使用CUDA并行计算模型设计实现了2D-3D医学图像配准算法的DRR并行快速生成过程，并在此基础上，由于SAD计算相对简单且模式强度PI能有效提高配准精度，本发明提出了一种基于SAD和PI的混合相似性测度函数并将该混合相似性测度函数也在CUDA上并行实现，最后在CPU上使用全局寻优能力强且计算方便的细菌趋化果蝇优化算法来进行整个配准过程的全局寻优，从而在保证算法运行效率的前提下提高了配准结果的精确性。

本发明的技术方案如下：

一种基于组合相似性测度的2D-3D医学图像并行配准方法，采用基于SAD和PI相结合的相似性测度M_new作为2D-3D医学图像配准的目标函数；其表达如式(1)所示：

$M_{\mathrm{new}}=\frac{1}{2}(-\mathrm{SAD}+\mathrm{PI})---\left(1\right)$

在式(1)中，平方差和SAD的计算方法如式(2)所示：

$\mathrm{SAD}=\frac{1}{N\×M}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}|R(i,j)-F_T(i,j)|---\left(2\right)$

在式(2)中,F_T(i,j)是进行空间变换后的浮动图像，N表示图像的宽度，M表示图像的高度。在式(1)中，模式强度PI的计算方法如式(3)所示：

$P_{r,\mathrm{\σ}}=\underset{i,j}{\mathrm{\Σ}}\underset{d^2\≤{\mathrm{rr}}^2}{\mathrm{\Σ}}\frac{{\mathrm{\σ}}^2}{{\mathrm{\σ}}^2+{(I_d(i,j)-I_d(v,w))}^2}---\left(3\right)$

其中：d²＝(i-v)²+(j-w)²   (4)

式(3)中，常数σ用来消除噪声的干扰，通常取σ＝10，半径r通常取3、4或5。式(3)中I_d为将X线图像I_X与DRR图像sI_DDR相减得到差值图像，如式(5)所示：

I_d＝I_r-sI_f   (5)

所述的基于组合相似性测度的2D-3D医学图像并行配准方法，包括以下步骤：

Step1:获取X线图像作为配准过程的参考图像，获取医学图像CT序列作为配准过程的浮动图像；

Step2:对X线图像进行去噪预处理，将医学图像CT序列图像构造生成三维体数据；

Step3:将处理后的X线图像数据和三维体数据装载到GPU中；

Step4:确定初始变换参数P₀＝(T_x,T_y,T_z,R_x,R_y,R_z)和CT值与衰减系数对应关系的查找表，以及SAD、SFD参数；

Step5:根据当前变换参数P在GPU中由多线程并行执行内核过程完成对三维体数据进行光线投射产生DRR图像；

Step6:在GPU中根据式(5)并行计算DRR图像与X线图像的混合相似性测度，并将计算得到的混合相似性测度值传递给CPU；

Step7:在CPU中利用基于细菌趋化行为的果蝇优化算法对相似性测度函数进行寻优；

Step8:如果当前相似性测度函数已达到最优，则输出变换参数P作为最终的配准结果；否则根据优化策略来更新当前的变换参数P，转至Step5执行。

所述Step7中的果蝇优化算法优化流程如下：

1)初始化果蝇种群规模SizeofPop、最大迭代次数MaxGen、初始迭代次数j＝0及所有果蝇个体的位置P_i,0＝(T_x,T_y,T_z,R_x,R_y,R_z)等参数；

2)由于果蝇个体主要利用嗅觉来搜索食物，针对每个果蝇个体随机赋予搜索食物的方向和距离，并采用式(7)在果蝇群体位置的基础上添加随机距离P来生成果蝇群体，其中P_i,j表示第i个果蝇在第j次迭代中的位置参数；

P_i,j＝P_i,j+RandomDist   (7)

3)当前迭代次数j是否小于最大迭代次数MaxGen，如果是则执行4)继续迭代，如果否则终止优化过程，转至10)；

4)将每个果蝇个体的位置参数P_i,j＝(T_x,T_y,T_z,R_x,R_y,R_z)传递到GPU中生成DRR图像，生成后并在GPU中根据式(10)计算X线图像和DRR图像的混合相似性测度函数值M_i,j，将计算得到的相似性测度函数值M_i,j传回给CPU，其中M_i,j表示第i个果蝇个体在第j次迭代中的相似性测度值；

5)根据式(8)对所有果蝇个体位置的对应的相似性测度值M_i,j求出其最小值和最大值，从而针对所有果蝇个体分别找出当前迭代中的混合相似性测度最优的果蝇和混合相似性测度最差的果蝇；

$\left\{\begin{array}{c}\left[\mathrm{bestMeasure\; bestIndex}\right]=\min \left(M_{i,j}\right)\\ \left[\mathrm{worstMeasure\; worstIndex}\right]=\max \left(M_{i,j}\right)\end{array}\right.---\left(8\right)$

6)根据式(9)记录并保留当前最佳相似性测度值bestMeasure与其P_b坐标，当前最差混合相似性测度值worstMeasure与其对应的P_w坐标值；

$\left\{\begin{array}{c}\mathrm{MeasureBest}=\mathrm{bestMeasure}\\ P_b=P_{\mathrm{bestIndex},j}\\ \mathrm{MeasureWorst}=\mathrm{worstMeasure}\\ P_w=P_{\mathrm{worstIndex},j}\end{array}\right.---\left(9\right)$

7)根据式(10)来计算果蝇群体的平均相似性测度值M_avg，然后根据平均相似性测度值M_avg利用式(11)来计算相似性测度方差σ²，并根据相似性测度方差σ²来判断果蝇群体朝最优个体的方向移动还是朝远离最差个体的方向移动；

$M_{\mathrm{avg}}=\frac{1}{\mathrm{SizeOfPop}}\underset{i=1}{\overset{\mathrm{SizeOfPop}}{\mathrm{\Σ}}}{M}_{i,j}---\left(10\right)$

${\mathrm{\σ}}^2=\frac{1}{\mathrm{SizeOfPop}}\underset{i=1}{\overset{\mathrm{SizeOfPop}}{\mathrm{\Σ}}}{(M_{i,j}-M_{\mathrm{avg}})}^2---\left(11\right)$

8)根据混合相似性测度方差来判断果蝇群体执行吸引操作还是逃离操作，从而来丰富种群多样性。若σ²≠0，则按式(12)更新所有果蝇个体位置P_i,j；否则按式(13)更新所有果蝇个体位置P_i,j；

P_i,j＝P_i,j+P_b   (12)

P_i,j＝P_i,j-P_w   (13)

9)当前迭代次数自增1，即j＝j+1，转至3)执行；

10)输出混合相似性测度函数的全局最优值及其对应的果蝇个体的位置，该位置参数即为2D-3D配准算法所求的最佳配准变换参数。

该算法首先使用CUDA并行计算模型完成DRR图像的快速生成过程，并组合差值绝对值和SAD与模式强度PI作为新的相似性测度在GPU上进行并行计算，最后将组合相似性测度值传递到CPU上采用基于细菌趋化行为的果蝇优化算法进行优化来寻找最优配准参数。通过实验对算法性能进行验证表明：由于本发明算法在GPU中实现DRR快速生成及混合相似性测度的计算，有效地提高了本发明算法的执行速度，同时与单一相似性测度相比，本发明采用混合相似性测度提高了配准结果的精确性。

附图说明

图1是基于CUDA实现DRR生成算法的主要流程；

图2是本发明2D-3D配准算法的主要流程。

具体实施方式

以下结合具体实施例，对本发明进行详细说明。

实施例1

CUDA是一种不需要借助图形学API就能使用类C语言进行通用计算的开发环境和软件体系结构,它的出现有效地发挥了GPU的并行计算性能，针对本质具有并行特征的算法能够有效地提高算法的执行效率。在基于光线投射的DRR生成算法中，每束光线穿透体数据来模拟X线入射人体后的衰减过程相互独立，具有良好的并行性，因此该过程可通过设计在GPU上执行的内核函数来实现。

由于生成DRR的输入数据是医学CT图像序列，其分辨率通常为512*512，且进行2D-3D配准过程中需要使用的DRR的分辨率也为512*512，因此可从光源发出512*512条射线，每条射线对应DRR图像的一个像素，将一条射线穿透体数据进行采样并累加计算CT值的过程封装在内核函数中由一个线程完成，最终所有射线穿过体数据进行光线衰减产生二维的DRR图像。为了充分发挥GPU的并行计算能力，可将整个DRR图像屏幕作为一个网格，并将屏幕分为若干块区域，屏幕上每个像素点对应一个线程的执行结果，这样可保证整个并行计算过程充分利用CUDA的三层体系结构。因此根据上述分析确定线程分配方案：网格大小为32*32，线程块大小为16*16，每条光线在X和Y方向上都有32*16＝512条射线，每条光线路径上的采样点计算和CT值累加计算过程由一个线程独立执行内核函数完成。在完成线程分配和数据存储之后，根据每条光线的投射采样和CT值衰减累加计算的过程来设计内核函数，在CUDA中多线程来执行内核函数实现光线投射过程的并行计算。基于CUDA实现DRR生成算法的主要流程如图1所示。

相似性测度函数是医学图像配准过程中非常重要的一步，相似性测度的选择直接影响到配准结果的精确性，因此寻找一种良好的相似性测度函数对整个2D-3D配准非常重要。在2D-3D医学图像配准算法中，常用的相似性测度有差值平方和SSD、差值绝对值和SAD、互信息MI、归一化互信息NMI、模式强度PI、梯度相关系数GC及梯度差分GD等。由于SAD只需要两幅图像对应像素求差，不需要进行导数运算，计算十分简单，在同模态图像配准中使用比较广泛，且在2D-3D配准中的两幅待配准图像X线图像和DRR图像的成像原理相似，可以看做同模态图像，故可采用SAD作为相似性测度函数；另外模式强度PI考虑了两幅待配准图像的像素灰度的同时还考虑了图像的空间信息，配准精度高且计算相对方便，因此为了获取更好的配准精度，结合待配准图像的空间信息和像素强度信息，在本发明算法中采用了基于SAD和PI相结合的相似性测度M_new作为本发明2D-3D医学图像配准的目标函数。其数学表达如式(2)所示：

$M_{\mathrm{new}}=\frac{1}{2}(-\mathrm{SAD}+\mathrm{PI})---\left(2\right)$

在式(2)中，平方差和SAD的计算方法如式(3)所示：

$\mathrm{SAD}=\frac{1}{N\×M}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}|R(i,j)-F_T(i,j)|---\left(3\right)$

在式(3)中,F_T(i,j)是进行空间变换后的浮动图像，N表示图像的宽度，M表示图像的高度。在式(2)中，模式强度PI的计算方法如式(4)所示：

$P_{r,\mathrm{\σ}}=\underset{i,j}{\mathrm{\Σ}}\underset{d^2\≤{\mathrm{rr}}^2}{\mathrm{\Σ}}\frac{{\mathrm{\σ}}^2}{{\mathrm{\σ}}^2+{(I_d(i,j)-I_d(v,w))}^2}---\left(4\right)$

其中：d²＝(i-v)²+(j-w)²   (5)

式(4)中，常数σ用来消除噪声的干扰，通常取σ＝10，半径r通常取3、4或5。式(4)中I_d为将X线图像I_X与DRR图像sI_DDR相减得到差值图像，如式(6)所示：

I_d＝I_r-sI_f   (6)

从式(4)可以看出：当该模式趋近于零时，模式强度测度的值趋近于一个最大值；当该模式增大时，模式强度测度趋于零。且模式强度采用函数1/(1+x²)的形式进行构造，根据该函数的渐进性质可知，当两幅图像的灰度差别很大时对该测度影响很小，从而当参考图像和浮动图像之间具有较大灰度差别时模式强度测度也能体现出良好的适应性。

新的相似性测度M_new在保持SAD计算相对简单的前提下引入了模式强度PI，可以有效地提高配准结果的精确性。本发明算法在对CT序列图像进行变换并生成DRR图像后，根据X线图像与DRR图像在GPU上通过执行内核函数来计算混合相似性测度M_new，将计算得到的混合相似性测度值传到CPU上，供优化算法完成配准过程的全局优化。

为了提高2D-3D医学图像配准算法执行效率的同时兼顾配准结果的精度，本发明提出了一种基于CUDA并行计算模型和组合相似性测度的2D-3D医学图像快速配准算法。本发明算法首先在支持CUDA并行计算模型的GPU上完成DRR图像的生成过程；并鉴于SAD计算简单的特点与模式强度PI能有效提高配准精度的优势，将两者进行组合作为新的相似性测度，并在GPU上完成参考图像X图像和DRR图像的组合相似性测度的计算；最后将混合相似性测度值传递到CPU上采用基于细菌趋化行为的果蝇优化算法来寻找最佳配准变换参数。

本发明2D-3D医学图像配准算法的主要流程如下：

Step1:获取X线图像作为配准过程的参考图像，获取医学图像CT序列作为配准过程的浮动图像；

Step2:对X线图像进行去噪预处理，将医学图像CT序列图像构造生成三维体数据；

Step3:将处理后的X线图像数据和三维体数据装载到GPU中；

Step4:确定初始变换参数P₀＝(T_x,T_y,T_z,R_x,R_y,R_z)和CT值与衰减系数对应关系的查找表，以及SAD、SFD参数；

Step5:根据当前变换参数P在GPU中由多线程并行执行内核过程完成对三维体数据进行光线投射产生DRR图像；

Step6:在GPU中根据式(6)并行计算DRR图像与X线图像的混合相似性测度，并将计算得到的混合相似性测度值传递给CPU；

Step7:在CPU中按照实施例2的优化流程利用基于细菌趋化行为的果蝇优化算法对相似性测度函数进行寻优；

Step8:如果当前相似性测度函数已达到最优，则输出变换参数P作为最终的配准结果；否则根据优化策略来更新当前的变换参数P，转至Step5执行。

本发明2D-3D配准算法的主要流程如图2所示。

实施例2

医学图像配准本质上是将相似性测度函数作为目标函数寻求最优解的过程。由于很多相似性测度函数都具有多个局部极值，因此在医学图像配准过程中使用一种具有很强全局寻优能力的智能优化算法，对寻找全局最优值和提高配准结果精度具有重要的意义。

果蝇优化算法(FruitFlyAlgorithm,FOA)是由中国台湾学者潘文超教授提出的一种模拟果蝇觅食行为的全局优化算法，该算法具有参数少、运行时间短、容易理解和实现等特点，被广泛应用于科学和工程很多领域。但是针对医学图像配准中的具有多极值特征的相似性测度函数进行优化时，FOA与遗传算法GA、粒子群算法PSO一样容易陷入局部最优，从而降低了医学图像配准的精确性。韩俊英等将细菌趋化过程中的吸引与排斥行为引入到FOA中，提出了基于细菌趋化的果蝇优化算法(BCFOA)，有效地解决了容易陷入局部最优的问题。

本发明在2D-3D配准算法的优化过程中采用基于细菌趋化的果蝇优化算法BCFOA来实现混合相似性测度函数寻优，主要优化流程如下：

1)初始化果蝇种群规模SizeofPop、最大迭代次数MaxGen、初始迭代次数j＝0及所有果蝇个体的位置P_i,0＝(T_x,T_y,T_z,R_x,R_y,R_z)等参数；

2)由于果蝇个体主要利用嗅觉来搜索食物，针对每个果蝇个体随机赋予搜索食物的方向和距离，并采用式(7)在果蝇群体位置的基础上添加随机距离来生成果蝇群体，其中P_i,j表示第i个果蝇在第j次迭代中的位置参数；

P_i,j＝P_i,j+RandomDist   (7)

3)当前迭代次数j是否小于最大迭代次数MaxGen，如果是则执行4)继续迭代，如果否则终止优化过程，转至10)；

4)将每个果蝇个体的位置参数P_i,j＝(T_x,T_y,T_z,R_x,R_y,R_z)传递到GPU中生成DRR图像，生成后并在GPU中根据式(10)计算X线图像和DRR图像的混合相似性测度函数值M_i,j，将计算得到的相似性测度函数值M_i,j传回给CPU，其中M_i,j表示第i个果蝇个体在第j次迭代中的相似性测度值；

5)根据式(8)对所有果蝇个体位置的对应的相似性测度值M_i,j求出其最小值和最大值，从而针对所有果蝇个体分别找出当前迭代中的混合相似性测度最优的果蝇和混合相似性测度最差的果蝇；

$\left\{\begin{array}{c}\left[\mathrm{bestMeasure\; bestIndex}\right]=\min \left(M_{i,j}\right)\\ \left[\mathrm{worstMeasure\; worstIndex}\right]=\max \left(M_{i,j}\right)\end{array}\right.---\left(8\right)$

6)根据式(9)记录并保留当前最佳相似性测度值bestMeasure与其P_b坐标，当前最差混合相似性测度值worstMeasure与其对应的P_w坐标值；

$\left\{\begin{array}{c}\mathrm{MeasureBest}=\mathrm{bestMeasure}\\ P_b=P_{\mathrm{bestIndex},j}\\ \mathrm{MeasureWorst}=\mathrm{worstMeasure}\\ P_w=P_{\mathrm{worstIndex},j}\end{array}\right.---\left(9\right)$

7)根据式(10)来计算果蝇群体的平均相似性测度值M_avg，然后根据平均相似性测度值M_avg利用式(11)来计算相似性测度方差σ²，并根据相似性测度方差σ²来判断果蝇群体朝最优个体的方向移动还是朝远离最差个体的方向移动；

$M_{\mathrm{avg}}=\frac{1}{\mathrm{SizeOfPop}}\underset{i=1}{\overset{\mathrm{SizeOfPop}}{\mathrm{\Σ}}}{M}_{i,j}---\left(10\right)$

8)根据混合相似性测度方差来判断果蝇群体执行吸引操作还是逃离操作，从而来丰富种群多样性。若σ²≠0，则按式(12)更新所有果蝇个体位置P_i,j；否则按式(13)更新所有果蝇个体位置P_i,j；

P_i,j＝P_i,j+P_b   (12)

P_i,j＝P_i,j-P_w   (13)

9)当前迭代次数自增1，即j＝j+1，转至3)执行；

10)输出混合相似性测度函数的全局最优值及其对应的果蝇个体的位置，该位置参数即为2D-3D配准算法所求的最佳配准变换参数。

实施例3

为了评估本发明2D-3D配准算法的执行效率和配准精度，本发明针对不同病人采集了多组临床放疗X线图像和制定放疗计划用的CT图像序列作为本实验的原始数据，在配准实验中X线图像作为参考图像，制定放疗计划用的CT图像序列作为浮动图像，通过对CT图像序列进行6自由度的刚性变换来生成DRR图像，将X线图像与DRR图像计算混合相似性测度并利用果蝇优化算法进行全局寻优，从而实现整个配准过程。由于本发明2D-3D医学图像配准算法中的DRR生成过程和混合相似性测度函数计算两大步骤需要在支持CUDA并行编程模型的GPU中计算，采用果蝇优化策略进行迭代寻优的过程在CPU上进行，故本实验在配置有Intercoreduoprocessor2.0*2GHzCPU、内存为2GB及GeForceGTX650Ti的显卡上进行，本实验环境的详细配置如表1所示。

表1 本实验的硬件环境配置

3.1配准过程执行效率

首先对本发明2D-3D配准算法的执行速度进行评价，并与基于CPU的2D-3D配准算法进行比较。本发明算法在支持CUDA的GPU上进行DRR快速生成以及X线和DRR图像的混合相似性测度计算，并在CPU上采用基于细菌趋化行为的果蝇优化算法进行整个配准过程优化，其中最大迭代次数MaxGen＝150，种群规模SizeofPop＝15，对果蝇群体的位置在区间[0,15]上进行随机初始化，迭代过程中果蝇搜寻最优位置的随机飞行方向与距离区间为[-1,1]；基于CPU的2D-3D配准算法在CPU上完成整个配准过程，包括DRR生成、相似性测度计算及优化过程。本实验中对8组实验数据分别使用了以上两种方法进行了配准，实验中选取X线图像的分辨率512*512，且在配准过程中在GPU上生成DRR图像的分辨率512*512，8组实验数据中的三维体数据的分辨率及使用两种配准算法的执行速度如表2所示。

表2 针对8组实验的执行效率比较

由表2可以看出，针对8组测试数据进行实验，与CPU下的2D-3D配准算法相比，本发明算法的计算速度提高了20倍左右。由于本发明2D-3D医学图像配准算法在支持CUDA并行计算模型的GPU上完成了DRR图像的快速生成过程及DRR图像与X线图像的混合相似性测度函数计算，充分发挥了GPU的并行计算性能，有效提高了2D-3D配准算法的运行效率。

3.2配准结果精度

为了有效测试本发明算法中提出的混合相似性测度对配准结果精确性的影响，针对第一组实验分别使用本发明算法及相似性测度分别为差值绝对值和SAD、互信息MI、模式强度PI的三种配准算法连续运行十次，然后分别针对六个自由度将四种算法的实际运行结果与理想位置做差来计算每次运行误差，最后分别根据式(14)、式(15)统计十次运行的平均误差和均方误差，平均误差的统计结果如表3所示,均方误差的统计结果如表4所示。在第一组实验中，理想的六自由度配准参数为(T_x,T_y,T_z,R_x,R_y,R_z)＝(3,0,5,0,10,0)。

$\stackrel{\‾}{X}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_i---\left(14\right)$

$\stackrel{\‾}{S}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{X})}^2}---\left(15\right)$

表3 四种不同相似性测度的平均误差统计结果

相似性测度	Tx	Ty	Tz	Rx	Ry	Rz
							差值绝对值和SAD	0.408	0.652	0.235	0.562	0.362	0.413
互信息MI	0.286	0.354	0.248	0.368	0.467	0.349
							模式强度PI	0.325	0.239	0.203	0.346	0.298	0.315
SAD+PI(本发明混合相似测度)	0.289	0.248	0.189	0.308	0.256	0.320

表4 四种不同相似性测度的均方误差统计结果

相似性测度	Tx	Ty	Tz	Rx	Ry	Rz
							差值绝对值和SAD	0.125	0.158	0.085	0.245	0.133	0.124
互信息MI	0.082	0.107	0.078	0.114	0.271	0.116
							模式强度PI	0.099	0.077	0.065	0.105	0.163	0.095
SAD+PI(本发明混合相似性测度)	0.091	0.081	0.059	0.087	0.124	0.097

从表3和表4中可以看出，对同一组参考图像和浮动图像分别采用四种不同的相似性测度进行十次配准，再对配准结果进行统计，本发明算法采用SAD与PI相结合的相似性测度，所得到的配准结果的平均误差及均方误差比其他三种单一相似性测度SAD、MI及PI都要小。由于本发明采用基于SAD和PI的混合相似性测度作为整个配准过程的目标寻优函数在GPU完成计算，充分发挥了图像中像素灰度信息和空间位置信息的作用，因此与差值绝对值和SAD、互信息MI、模式强度PI三种单一相似性测度函数比较，有效地提高了配准结果的精确性。

综上对本发明算法的执行效率和精确性进行验证可知，本发明算法在CUDA上计算DRR生成过程及混合相似性测度，在有效利用2D-3D配准算法的并行性的基础上充分发挥了CUDA的并行计算优势，显著地提高了配准速度，同时采用基于SAD和PI的混合相似性测度，并采用具有很强全局寻优能力果蝇优化算法进行寻优，在几乎不影响配准速度的同时有效地提高配准结果的精确性。

评估结论

本发明针对2D-3D医学图像配准算法在图像引导放射治疗中的运行效率和配准精度要求，提出了一种基于CUDA并行计算模型和组合相似性测度的2D-3D医学图像快速配准算法。该算法使用CUDA并行计算模型完成DRR图像的生成和混合相似性测度函数的计算，并将相似性测度值传递到CPU上采用具有细菌趋化行为的果蝇优化算法寻找最佳配准参数。通过实验对算法验证表明：本发明算法将DRR生成过程和混合相似性测度计算两步都在GPU中进行，同时使用了基于SAD和PI的混合相似性测度和果蝇优化算法，有效地保证了算法的执行速度的同时改善了配准结果的精确性。从实验结果可以看出，本发明算法已有效的提高了2D-3D医学配准过程的速度，但是还有进一步改善的空间。后续工作将对DRR生成过程尝试使用其他体绘制方法实现进一步优化和加速，以及对DRR生成过程中的采样方式进行优化处理，期望能够进一步提高算法的实时性。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (2)

1.一种基于组合相似性测度的2D-3D医学图像并行配准方法，其特征在于，

采用基于SAD和PI相结合的相似性测度M_new作为2D-3D医学图像配准的目标函数；其表达如式(1)所示：

$M_{\mathrm{new}}=\frac{1}{2}(-\mathrm{SAD}+\mathrm{PI})---\left(1\right)$

在式(1)中，平方差和SAD的计算方法如式(2)所示：

$\mathrm{SAD}=\frac{1}{N\×M}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{M-1}{\mathrm{\Σ}}}|R(i,j)-F_T(i,j)|---\left(2\right)$

在式(2)中,F_T(i,j)是进行空间变换后的浮动图像，N表示图像的宽度，M表示图像的高度；在式(1)中，模式强度PI的计算方法如式(3)所示：

$P_{r,\mathrm{\σ}}=\underset{i,j}{\mathrm{\Σ}}\underset{d^2\≤{\mathrm{rr}}^2}{\mathrm{\Σ}}\frac{{\mathrm{\σ}}^2}{{\mathrm{\σ}}^2+{(I_d(i,j)-I_d(v,w))}^2}---\left(3\right)$

其中：d²＝(i-v)²+(j-w)²   (4)

式(3)中，常数σ用来消除噪声的干扰，通常取σ＝10，半径r通常取3、4或5；式(3)中I_d为将X线图像I_X与DRR图像sI_DDR相减得到差值图像，如式(5)所示：

I_d＝I_r-sI_f   (5)

所述的基于组合相似性测度的2D-3D医学图像并行配准方法，包括以下步骤：

Step1:获取X线图像作为配准过程的参考图像，获取医学图像CT序列作为配准过程的浮动图像；

Step2:对X线图像进行去噪预处理，将医学图像CT序列图像构造生成三维体数据；

Step3:将处理后的X线图像数据和三维体数据装载到GPU中；

Step4:确定初始变换参数P₀＝(T_x,T_y,T_z,R_x,R_y,R_z)和CT值与衰减系数对应关系的查找表，以及SAD、SFD参数；

Step5:根据当前变换参数P在GPU中由多线程并行执行内核过程完成对三维体数据进行光线投射产生DRR图像；

Step6:在GPU中根据式(5)并行计算DRR图像与X线图像的混合相似性测度，并将计算得到的混合相似性测度值传递给CPU；

Step7:在CPU中利用基于细菌趋化行为的果蝇优化算法对相似性测度函数进行寻优；

Step8:如果当前相似性测度函数已达到最优，则输出变换参数P作为最终的配准结果；否则根据优化策略来更新当前的变换参数P，转至Step5执行。

2.根据权利要求1所述的方法，其特征在于，所述Step7中的果蝇优化算法优化流程如下：

1)初始化果蝇种群规模SizeofPop、最大迭代次数MaxGen、初始迭代次数j＝0及所有果蝇个体的位置P_i,0＝(T_x,T_y,T_z,R_x,R_y,R_z)等参数；

2)由于果蝇个体主要利用嗅觉来搜索食物，针对每个果蝇个体随机赋予搜索食物的方向和距离，并采用式(6)在果蝇群体位置的基础上添加随机距离来生成果蝇群体，其中P_i,j表示第i个果蝇在第j次迭代中的位置参数；

P_i,j＝P_i,j+RandomDist   (6)

3)当前迭代次数j是否小于最大迭代次数MaxGen，如果是则执行4)继续迭代，如果否则终止优化过程，转至10)；

4)将每个果蝇个体的位置参数P_i,j＝(T_x,T_y,T_z,R_x,R_y,R_z)传递到GPU中生成DRR图像，生成后并在GPU中根据式(1)计算X线图像和DRR图像的混合相似性测度函数值M_i,j，将计算得到的相似性测度函数值M_i,j传回给CPU，其中M_i,j表示第i个果蝇个体在第j次迭代中的相似性测度值；

5)根据式(7)对所有果蝇个体位置的对应的相似性测度值M_i,j求出其最小值和最大值，从而针对所有果蝇个体分别找出当前迭代中的混合相似性测度最优的果蝇和混合相似性测度最差的果蝇；

$\left\{\begin{array}{c}\left[\mathrm{bestMeasure\; bestIndex}\right]=\min \left(M_{i,j}\right)\\ \left[\mathrm{worstMeasure\; worstIndex}\right]=\max \left(M_{i,j}\right)\end{array}\right.---\left(7\right)$

6)根据式(8)记录并保留当前最佳相似性测度值bestMeasure与其P_b坐标，当前最差混合相似性测度值worstMeasure与其对应的P_w坐标值；

$\left\{\begin{array}{c}\mathrm{MeasureBest}=\mathrm{bestMeasure}\\ P_b=P_{\mathrm{bestIndex},j}\\ \mathrm{MeasureWorst}=\mathrm{worstMeasure}\\ P_w=P_{\mathrm{worstIndex},j}\end{array}\right.---\left(8\right)$

7)根据式(9)来计算果蝇群体的平均相似性测度值M_avg，然后根据平均相似性测度值M_avg利用式(10)来计算相似性测度方差σ²，并根据相似性测度方差σ²来判断果蝇群体朝最优个体的方向移动还是朝远离最差个体的方向移动；

$M_{\mathrm{avg}}=\frac{1}{\mathrm{SizeOfPop}}\underset{i=1}{\overset{\mathrm{SizeOfPop}}{\mathrm{\Σ}}}{M}_{i,j}---\left(9\right)$

${\mathrm{\σ}}^2=\frac{1}{\mathrm{SizeOfPop}}\underset{i=1}{\overset{\mathrm{SizeOfPop}}{\mathrm{\Σ}}}{(M_{i,j}-M_{\mathrm{avg}})}^2---\left(10\right)$

8)根据混合相似性测度方差来判断果蝇群体执行吸引操作还是逃离操作，从而来丰富种群多样性；若σ²≠0，则按式(11)更新所有果蝇个体位置P_i,j；否则按式(12)更新所有果蝇个体位置P_i,j；

P_i,j＝P_i,j+P_b   (11)

P_i,j＝P_i,j-P_w   (12)

9)当前迭代次数自增1，即j＝j+1，转至3)执行；

10)输出混合相似性测度函数的全局最优值及其对应的果蝇个体的位置，该位置参数即为2D-3D配准算法所求的最佳配准变换参数。