基于GPGPU的锥束螺旋CT DBPF重建算法的实现与加速方法研究

沈 乐 邢宇翔,2

1 (清华大学工程物理系 北京 100084)

2 (清华大学粒子技术与辐射成像教育部重点实验室 北京 100084)

锥束螺旋X射线CT近年来受到广泛关注，其重建算法是当前CT领域的研究热点。目前，解析重建方法有近似重建方法和精确重建方法，后者是大家追求的目标[1,2]。2004年，ZOU 等[3]提出了微分反投影滤波重建方法(DBPF)，能在最少数据条件下进行精确重建，很多相关工作(包括轨道的推广和局部重建等)随之展开。但其具体实现面临计算复杂度和离散化等问题，至今，商用系统仍未使用此类方法。因此，寻找一种好的算法实现策略以及如何对算法进行加速使其应用于实际非常关键。

利用显卡硬件的并行能力对CT重建算法加速是目前加速效率最高的一种方式。由于重建算法具有数据相互独立性及可高度并行处理特点，在GPGPU上设计并行重建算法，能获得很好的性能提升。对于平行束、扇形束的FBP算法，以及锥束FDK算法，已获得良好结果[4–6]。对于PI线BPF重建算法，文献[7,8]用 Cg语言处理重建的反投影过程，重建速度提升了100倍。

我们对DBPF方法的实施技巧和硬件加速进行了研究，给出了一种新的实现方法和加速策略，以在最少不均匀采样基础上提供可控分辨率重建，并最大程度地利用 GPU的并行加速性能全面提高容积重建速度。

1 DBPF算法和PI线离散方法

螺旋锥束CT的扫描坐标系和参数的定义如图1所示，射线源和探测器按螺旋轨道u运v动，射线源的位置在直角坐标系{x,y,z}中表示为Rsinλ,hλ/2π)，其中，R为射线源旋转的半径，λ为旋转角度，h为螺距，即射线源旋转一周相对于置物台的轴向移动距离。定义探测器旋转坐标系{u,v,w}，其基向量为：

图1 螺旋锥束CT坐标系定义Fig.1 Definition of coordinates for the helical conebeam CT.

2004年芝加哥大学潘晓川教授[3]提出的基于PI线的反投影滤波算法的步骤如下：

(1) 对投影数据求导：

(2) 沿PI线进行加权反投影：由求导后的数据G{λ,u,v}在 PI线所在角度范围λ1–λ2进行积分：

(3) 对加权反投影结果进行Hilbert变换，得到PI线上的重建结果：

这里，附加项xA,xB为PI线与重建区域的交点，P{u,v,λ}为PI线所在位置的射线投影，我们取P{u,v,λ} = [P{u,v,λ1}+P{u,v,λ2}]/2。

(4) 将PI线重建结果按三线性插值方式重新采样成为均匀网格的三维重建图像。在BPF重建中，常规的 PI线选取为从一个光源角度采样点到其它多个光源角度采样点形成的一个扇形曲面形式，这不利于建立重建图像分辨率与PI线密度和PI线上采样密度的关系。如图2，使一组PI线(图2中实线)在俯视图中(图2c)互相平行且距离为DG，采样点在XY平面上投影的距离为DP。我们把这样一组PI线构成的曲面称为一个PI面，图2的虚线所示构成另一个PI面。在选择轴方向上相邻的PI面时，令它们的PI线在XY平面上的投影成一角度Dθ，由重建需要的层厚确定。取两个PI面的中心点(0,0,z)作为基准点(也可选其它点作为基准点，如图2b所示，效果相同)，则Dθ由Z轴采样间隔DZ和螺距的线性关系确定。如此得到的PI线采样点在重建区域内具有全局一致性，可方便地控制最后的重建分辨率。如DG=DP=DZ，可得到各向同性的具有均匀分辨率效果的 PI线重建采样点。用式(4)得到的重建结果是沿PI线的数据，须按照规则网格重新采样为三维体数据。由于我们的PI线选取方式得到全局基本均匀的重建点，可方便地使用三线性插值方式重采样得到重建图像。

731 Comparison of different separation protocols for exosomes derived from urine of patients with aortic dissection

图2 PI 线选取和采样点 (a)全景图，(b)侧视图，(c)俯视图Fig.2 PI-lines and sampling point selection. (a) Perspective view, (b) Top view, (c) Side view.

式中，(xΠ,yΠ,zΠ)表示 PI线上像素的坐标，ax,ay,az为该像素与均匀网格像素之间距离的三个分量，D为均匀网格的层内离散化步距，Ds为均匀网格的层厚，Neighbor(i,j,k)表示(i,j,k)邻域(图 3)。

图3 三线性插值重采样的示意图Fig.3 Illustration of tri-linear interpolation resampling.

2 CUDA加速优化策略

随着半导体技术及芯片制造工艺的不断突破，图形处理器的流式核心数量、浮点运算能力不断提升。近年来快速发展的通用图形处理器(GPGPU)除图形图像的显示处理外，还可完成很多原本由CPU完成的计算任务。目前的高端图形处理器(如美国NVIDIA公司的 GF100核心和美国 AMD公司的RV870核心)的浮点运算能力已超过2 TFLOPS(1012浮点运算/s)，已大大超越多核CPU(如美国Intel公司芯片)的计算能力(表 1)。图形处理器本身系并行处理器，对于可并行算法，通用图形处理器的执行速度远高于CPU。

CUDA(Compute Unified Device Architecture)[9]是美国NVIDIA公司提出的并行架构，可利用图形处理器的并行处理器能力大幅度提升计算性能。CUDA包含一系列硬件架构、一套指令集和编程模型，针对C/C++/Fortran语言的应用程序编程接口和编译器。开发人员利用相应软件工具，可在CUDA架构的GPU上开发并行程序应用及加速算法。

我们利用CUDA架构，在GPU上完成：(1) 对投影数据求导；(2) 沿 PI线加权反投影；(3) 有限Hilbert变换；(4) 将PI线重建结果按图像网格重采样。具体步骤如下：

(1) 求导：将投影数据传入显卡的纹理显存中，将探测器面阵所得的投影数据映射到一个二维线程组中，每一个线程对应计算某一点的全导数。

(2) 反投影：在一个PI面内的一组PI线，不仅起点终点角度不同、长短不同，PI线上的采样点数目也不同，因此无法将这些PI线上的采样点直接映射为规则的二维线程组。若采用最长PI线的采样点数作为一个映射维度，则其他PI线映射的线程又不能充分利用，对于显卡并行模型是不经济的。CUDA并行编程中应尽可能使显卡上的线程达到负载平衡。为此，我们沿Z方向选取长度相同(图4)，即起始和终止角度差Dλ=λ2−λ1相等的一组 PI线并发执行反投影操作。不同的PI线仅λ1、λ2的参数不同，长度及采样点数量都相同。实验证明此法能有效匹配负载平衡，提高流处理器的利用率，大大加速反投影运算。

(3) 有限Hilbert变换：将一组PI线反投影结果映射到二维线程组中，每一线程计算一个采样点的卷积结果。计算前将反投影数据和卷积核绑定到一维纹理中，以利用纹理缓存命中提高数据读取速度。

(4) 重采样：将一组PI线的上述结果映射到二维线程组中，按图3方式将结果分配到采样点周围的8个邻域像素上，并累加8个像素所占权重值。当重建全部完成后，使用累加的权重值对结果归一化处理，得到最终图像。

图4 反投影的并行方式示意图Fig.4 Illustration of paralleled back-projection.

为提高算法效率，我们提出了优化手段，包括打包存储重建所需常量参数，应用支持线性插值方式的三维纹理存储器保存投影数据等，取得了较好效果。

首先，在加权反投影及Hilbert变换中可能用到的重建参数超过20个。若全部作为函数调用参数来传递，则kernel的寄存器占用数会增加27%，流处理器占用率明显减少。对于计算能力<1.3的硬件，可能会因硬件资源不足导致kernel无法启动。为减少寄存器占用，我们将重建所需的参数打包存于数组中，预先输入显卡的常数存储器或共享存储器中以供调用，有效解决了每个kernel硬件资源占用过多的问题。

其次，CUDA的纹理存储器是一种只读存储器，具有缓存机制，同时支持硬件线性插值。直接访问三维浮点坐标的数据即可得按三线性插值的结果，无需编程实现。缓存机制还可提高内存的实际访问带宽，降低延迟。

最后，在并行重采样时，若PI线上的采样点间隔较小，相邻的采样点可能具有公共的邻域像素。此时若并行进行重采样，会引起不同线程对同一内存地址进行写操作。为避免线程间的访问冲突，在反投影滤波结束后，将PI线重建数据传回内存，由CPU进行串行重采样。这一步可在子线程上进行，而CPU主线程则控制GPU进行下一轮反投影滤波运算，如此可充分利用GPU计算时CPU的空闲时间。只要在数据传输前设置同步，就可保证访问是安全的。这种CPU+GPU同时计算的方法，在重建数据较小、PI线采样步距较大时可有效隐藏重采样时间。

对于大数据量处理，同步导致的空闲时间明显增加。为此，我们设计了第二种GPU上执行的方法：分批并行重采样，在加速的同时避免写操作的冲突。在二维线程组映射时，把 PI线上的采样点分批处理。如图5所示，圆点标记的PI线采样点映射为一个批次的线程组，完成后把方形标记PI线采样点映射为另一批次的线程组进行计算，如此继续直至处理完所有的PI线上采样点。选取每一个批次内处理的PI线采样点的合适间距，可保证采样点不会有公共的被重建邻域像素，避免了线程访问冲突。

图5 并行重采样方式示意图Fig.5 Illustration of paralleled resampling.

3 结果与分析

仿真实验中，我们采用三维Shepp Logan 模型生成螺旋锥束投影数据，光源到旋转轴和探测器的距离分别为650 mm和1000 mm，每个2π采样1080个角度，探测器为 256×80，每个单元尺寸为 1.8 mm×1.8 mm，扫描螺距46.7 mm。算法使用Visual Studio 2005，CUDA Toolkit 3.0编程设计，程序在一台AMD Athlon II X4 620 (2.60 GHz)计算机上进行测试，具有DDR2 800 MHz 2 GB RAM和NVIDIA Tesla C1060 GPU。我们分别使用了纯CPU执行和GPU加速测试重建256×256×64的结果，D = 1 mm,采样间隔DG=DP=DZ=0.8D。有效重建区域覆盖52层(3.4×106像素)，在CPU上重建时间787 s，使用GPU加速后的重建时间为2.469 s，重建加速比达到318倍。图6为重建图像比较，显示CPU与GPU重建结果非常吻合。

图6 第28层重建结果比较(a) CPU重建结果, (b) CUDA加速后的结果, (c) 图像(a)(b)第92行的水平剖面线Fig.6 Reconstructions of Slice 28.(a) CPU result, (b) CUDA accelerated result, (c) Horizontal profiles at Row 92.

在以上的重建中，重建速度及质量与 Hilbert卷积核选取、附加项计算、采样间隔等参数密切相关。为避免混叠，卷积核长度通常为2Np−1，Np表示 PI线上的采样点数目。缩短卷积核长度可加快Hilbert变换的运算速度，但会引入伪影。固定最后的重建分辨率Δ条件下，重建过程中PI线间的距离DG，PI线上的采样间距ΔP，Z轴方向采样间隔DZ三个参数的取值会对重建结果及速度产生影响。在需要重建图像各向同性分辨率情况下，通常取(但不严格要求)DG=DP=DZ。显然，采样间隔越小，重建图像质量越好，重建时间越长。根据经验，一般可设采样间隔(DG, DP, DZ)～0.8D，即使 PI线的采样略细于最后的重建分辨率要求即可得到满意的重建图像和高的重建速度。图7显示不同采样参数组的重建效果。图8是不同参数情况下的重建时间变化情况，这里我们设置ΔG= ΔP。

图7 第20层不同采样参数的重建图像Fig.7 Reconstructions of Slice 20 with different parameters.

图8 不同PI线选取和采样网格下的重建时间比较(ΔG=ΔP的情况)Fig.8 Reconstruction time comparison using different sampling grids for PI line setting (ΔG = ΔP).

4 结论

本文从DBPF重建算法的原理和目前通用图形处理器与并行重建算法的技术入手，介绍了一种新的实现方法，以及使用CUDA平台下算法进行高效加速的策略和优化手段。充分利用了GPU的并行处理与CPU+GPU混合计算的计算能力，达到了很好的加速效果。

该算法仍有进一步改进的可能，重建过程中的计算精度由于使用内置快速数学函数会带来一定误差。2010年4月美国NVIDIA公司发布的最新图形处理核心 GF-100具有更高的单精度和双精度浮点运算能力，并且支持IEEE 754-2008单精度浮点数标准。在最新架构下，可尝试使用高精度的标准数

学函数，以获取更优的图像质量。此外，反投影过程是所有步骤中耗时最长的操作，有可能通过继续提高反投影过程中的纹理缓存命中率，减少分支判断而提高速度。我们将在后续工作中继续研究这些技术。

1 Tuy H, SIAM J. Appl Math, 1983, 43(3): 546–552

2 Katsevich A. Phys Med Biol, 2002, 47: 2583–2597

3 ZOU Yu, PAN Xiaochuan. Phys Med Biol, 2004, 49:941–959

4 Wang Y J, Hu H F, Xing Y X. Strategy for GPU acceleration of massive data cone beam CT reconstruction.The 10thinternational meeting on fully 3D image reconstruction in radiology and nuclear medicine, Beijing,2009. 57–60

5 Bi W Y, Xing Y X, Chen Z Q,et al. 2008 IEEE Nucl Sci Symp and Med Imag Conf (2008 NSS/MIC). 2009, 1–9:2605–2608

6 Xu F, Mueller K. Phys Med Biol, 2007, 52: 3405–3419

7 Zheng H, Yu Y, Kang Y,et al. Proc SPIE, 2010, 7622:76222G-1–76222G-12

8 Zheng H, Kang Y, Liu J,et al. Implementation of Helical Cone-Beam Back-Projection Filtered Reconstruction Algorithm on GPU. The 10thinternational meeting on fully 3D image reconstruction in radiology and nuclear medicine, Beijing, 2009. 45–48

9 NVIDIA Corporation. NVIDIA CUDA Best Practices Guide [EB/OL]. http://developer. download. nvidia.com/compute/cuda/3_0/toolkit/docs/NVIDIA_CUDA_BestPra cticesGuide.pdf, 2010