Nvidia从GeForce8系列开始，在显卡上推出统一计算设备框架技术，使GPU的通用计算（GPGPU）从图形硬件流水线和高级绘制语言中解放出来，开发人员无须掌握图形学编程方法即可在单任务多数据模式（SIMD）下完成高性能并行计算。在医学图像分析中，图像配准通常是一个耗时的过程，不利于临床应用，为了加速医学图像的2D-3D配准过程，研究了CUDA的设计思想和编程方式，提出了一种基于CUDA并行编程模型的加速配准新技术，在构建的虚拟X线摄像系统下，采用并行计算的方式快速生成高质量DRR图像，以对应像素的灰度值残差作为相似性测度，使用Powell优化方法寻找最优变换。实验结果表明，该技术既很好地保持了配准精度，同时又大大提高了配准速度，加速比达到了十几甚至几十倍。
2010,46(11) Computer Engineering and Applications计算机工程与应用 换量,提高配准速度;iwzλwo指临床ⅹ线图像(标准图像); GPU Kernel output中存储输出图像;cosh记录临床Ⅹ线图像与DRR图像 GPU/Device 间对应像素的灰度值残差;,h分别表示配准图像的宽和高; Grid density, bright, offset, scale保证生成的DRR图像的质量,配准 Block Block Block 过程中不需要变动 (1,0) 计算线程索引 每个线程在GPU上执行的流程如图5所示,由于每个线 Block Block Block 程处理的数据对象和线程之间是一一对应的,图5首先计算该 (0, (1,1) 视线方向的空间变换 线程的索引,接着进行视线方向的空间变换,并依据GPU自带 算DRR图像的像素灰度值 Block (1, 1) 的三线性插值方式确定穿过CT体数据的光线上的重采样点 的CT值,用来计算DRR图像上对应象素点的灰度值,最后计 Thread(0, 0)Thread(1, 0) 计算相似性测度 算ⅹ线图像与DRR图像间对应像素的灰度差值,输出残差供 Thread(0,1)Thread(1,1) 输出残差,记录当前图像 CPU优化使用,并记录由此时的变换矩阵生成的DRR图像,方 图4CUDA编程模型图5GPU上每个线程的处理流程图便CPU保存图像时使用。 设数据大小为mXxp,配准并行化实现步骤 1)CPU将三维数据场读入主机内存中,然后调用CUD3实验结果及讨论 函数将数据拷贝到设备内存中。 该文所用计算机配置为:内存1.0G,显卡是 Geforce9800 (2)根据显卡处理器的数目和生成图像的大小指定Gid GTX/9800GX+,其静态存储器大小为512M,该显卡拥有16 和 Block大小。在CUDA编程模式下,每个 Block块中的线程个多处理器,每个多处理器中含有8个流处理器。实验分为模 的总数最大不超过512,而GP发挥最佳性能的线程总数是拟配淮实验和临床ⅹ线图像配准实验两部分,实验编程环境 64~256个,结合医学图像的大小,为每个Bock分配256个 为 Windows XP下的Ⅴ isual studio o+NET205,实验中以两 线程。 幅图像的差值以及拼接结果作为评估配准效果好坏的依据。 (3)根据设定的优化初始值,调用 OpengL函数生成变换3.1模拟配准实验 矩阵T,并将T传入GPU上进行空间变换。 模拟实验部分,首先生成一幅DRR参考图像,记录此时的 (4)GPU上按分配好的线程进行并行操作,并计算相似性位置参数,将其作为评估配准精确性的金标准,然后对CT体 测度的值。 数据进行空间变换,通过优化得到最佳的空间位置,最终实现 (5)判断(4)中得到的值是否最优,若是,则停止搜索,输出 -3D的配准。其中CT体数据是16位的膝部裸数据,大小为 此时的最优位置,否则,变换位置生成新的变换矩阵T并将其256×256×256。空间变换的单位:平移为0.7032mm,旋转为度。 传入GPU,然后转(4) 32临床配准实验 (6)输出最佳变换后的DRR图像。 临床ⅹ线图像配准部分与模拟实验部分步骤相同,不同 用到的线程索引计算公式如下 处是采用临床X线图像作为参考图像,CT体数据是16位的 uinx=_umul24( blockade.rx, block Dim.x)+ hreadida.x(5)膝鄙裸数据,大小为512×512×360。空间变换的单位:平移为 uint y=_ml24( blockade.y, blockDim.y)+ threadedy(6)0.3516mm,旋转单位为度。 下面所示的是图像配准的核心函数 实验中采用 Powell)法,寻找最佳的配准参数Tx,Ty, d render<< midsIze, block size, Share>>( biaozhun, output,Tz,Rx,Ry,Rz,其中Tx,7y,7z表示沿相应轴的平移距离,Rx cost, w, h, density, bright, offset, scale) Ry,Rz表示绕相应轴的旋转角度。由于沿X轴的平移对成像影 其中, midsIze, block size分别表示设备中分配的Grid和 Block响比较小,对参数的搜索顺序为(Ty,T,Rx,Ry,Rz,Tx) 的个数; Share表示GPU上分配的共享存储器大小,目的是将 从表1的实验数据可以看出,一方面,基于CUDA的配准 差值结果存储到共享存储器,减少CPU与GPU之间的数据交保持很高的配准精度,配准精度达到了亚像素机,这一点也可 表1膝部模拟实验 第2组 第3组 第4组 第5组 理想位置(0,20,20,10,10,10)(0,10,10,0,0,0)(0,10,-10,10,10,10)(0,10,10,5,5,5)(0,10,9,8,7,6) 初始位置 (0,0,0,0,0,0) (0,0,0,0,0,0 (0,0,0,0,0,0) (0,0,0,0,0,0)(0,0,0,0,0,0) -0.034 0013 0.024 -0.012 0.001 19.998 9.980 9.998 9.980 9.999 Tz 19.990 9.998 -10.003 9.997 R 9.997 0.000 10.003 5.000 7.999 R 10.000 0.003 10.006 5.003 7.000 R 10.003 0.013 9.998 4.988 5.999 GPU配准时间/s 7.3 10.9 10.4 CPU配准时间/s 260.2 300.4 418.0 322.2 174.9 加速比 14.06 41.15 38.35 30.98 14.70 徐建,秦安,卜祥磊,等:基于CUDA的2D-3D配准技术的研究 2010,46(11)59 (a)参考图像(b)配准前的DRR图像(c)配准后的DR图像 (a)临床ⅹ线图像(b)模拟DRR图像(c)配准后的DRR图像 (模拟Ⅹ线图像) (d)配准前差值图像(e)配准后差值图像(f)配准后的拼接图像 (d)配准前的差值图像(e)配准后的差值图像(f)配准后的拼接图像 图7临床X线图像配准的效果图 图6模拟配准的效果图 注:为了显示清晰的拼接效果,人为设定两幅图像中物体重叠部分的 应用D]南京:东南大学,2007 灰度值相差20,拼接效果如图(f)所示。 [3] Khamene A, Bloch P Automatic registration of portal images and volumetric CT for patient positioning in radiation therapy[J]. Medi- 以从配准后的效果图看出来,如图6所示;另一方面,基于 cal Image Analysis, 2006, 10: 96-112 CUDA的配准速度是基于CPU的配准的十几倍到几十倍。其14 Sherouse g, Novins k, Chaneyetal B. Computatipn of digitally recon 实,随着数据规模的增大,计算量的増多,CUDA并行计算的优 structed radiographs for use in radiotherapy treatment design [J].Int 势将会发挥得更加彻底,从而加速效应也将会更加明显。从实 J Radiation Oncology Biol Phys,1990,18(3):651-658 验效果图可以看出,该技术不仅对模拟实验到达很好的配准效{5]GalⅷinJM, Sims O. The use of digitally reconstructed radiographs 果,对于临床实验,配准效果也是非常不错的,这主要是由于生 for three-dimensional treatment planning and CT-simulation[JIn 成的DRR图像的质量很接近于临床ⅹ线图像的缘故。临床Ⅹ ternational Journal of Radiation Oncology*Biology*Physics, 1995 线图像的配准效果如图7所示。从上面的分析可以得出:基于 31:935-942. CUDA的并行加速配准技术,在保持了配准精度的前提下,很6梁玮鲍旭东,罗立民基于互信息的2D-3D医学图像配准门生物 好地加速了配准的进程。 医学工程研究,2004,23(1):14-16 「η张薇,黄毓瑜,栾胜,等基于灰度的二维厂三维图像配准方法及其在 4结束语 骨科导航手术中的实现小中国医学影像技术,2007,23(7):1080- DRR图像的生成是2D-3D配准过程中的重要组成部分, 直接关系到配准结果的好坏。由于DRR图像的生成是基于光 [8 Rohlfing T, Russakoff D B, Denzler J, et al. Progressive attenuation 线投射算法,所以非常耗时,提出了一种基于CUDA加速的 fields: Fast 2D-3D image registration without precomputation [JI Medical Physics, 2005, 32(9): 2870-2880 医学图配准新技术,加速了DRR图像的生成,实现了2D-3D [9]刘磊基于GPU的医学图像三维重建及可视化技术研究D]广州 图像的自动快速配准,该技术对CT及ⅹ线图像(临床和模拟 南方医科大学,2008 进行了实验得到了不错的配准结果。此技术也适合其他模态10桂叶晨,冯前进刘磊等基于CUDA的双三次B样条缩放方法 图像的配准。从实验结果可以看出,该技术虽然能够大大提高 计算机工程与应用,2009,45(1):183-185 2D-3D配准的速度,但离临床的实时性要求还是有一定的差 [ll] Kubias A, Deinzer F, Feldmann T, et al. 2D/3D Image Registration 距。下一步可对体数据进行预处理以减少计算量,并进一步提 on the GPU[J]. Pattern Recognition and Image Analysis, 2008, 18 高DRR图像的质量,提高配准的精度。同时,寻找一种快速算 (3):381-389 法实现图像的粗配准,在此基础上进行精确寻优,进一步优化12] NVIDIA CUDA Computational Unified Device Architecture Pro 代码的执行,期望以该技术为依托的2D-3D配准能够满足临 gramming Guide Version 2.0[.[S.: NVIDIA Corporation, 2008 床的实时性要求。 [13] Mathews J H, Fink KK Numerical methods using Matlab[M]. 4th ed. S.I. Prentice Hall, 2004 434-439 参考文献: 41薛毅最优化原理与方法[M北京:北京工业大学出版社,2004 刘洋医学图像2D3D配准及可视化研究[D]广州:南方医科大学, 215-221 2007 ∏15]杨虎,马斌荣,任海萍,等基于互信息的人脑图像配准研究J中 [2]陆忞X线透视图像(2D)与CT体积图像(3D)配准方法研究及其 国医学物理学杂志,2001,18(2):69-73.