一种利用Spark-GPU加速CT图像重建的设计

熊威， 曾有灵， 李喆

(南方医科大学 生物医学工程学院， 广东 广州 510515)

CT(computed tomography, CT)图像重建速度是衡量CT系统的重要指标.使用高性能计算硬件是加速图像重建的重要手段，目前主要有使用CPU集群系统重建法、使用FPGA重建技术、使用GPU加速重建法以及使用细胞宽带引擎(cell broadband engine,CBE)重建法[1].GPU拥有数以千计的计算核心，能同时创建成千上万的线程，非常适合处理图像重建中巨大的计算量.许多学者对GPU加速CT图像重建进行了研究，Yan等[2]使用GPU加速FDK算法，并通过循环渲染到纹理、结合Z轴对称和多个渲染目标技术降低了切片几何映射和投影的计算成本，提升了重建速度.Lu等[3]使用GPU加速同步代数迭代算法，采用射线驱动投影和体素驱动反投影技术使重建速度和常用的滤波反投影算法的速度相当.

统一设备计算框架(compute unified device architecture,CUDA)是英伟达(NVIDIA)公司推出的支持NVIDIA GPU的通用并行计算架构，随着GPU性能的提高和CUDA的成熟，GPU集群并行计算成为研究热点，Fan等[4]率先开发出可扩展的GPU集群，用于高性能科学计算和大规模仿真.消息传递接口MPI(message passing interface,MPI)是目前广泛使用的并行计算开发环境.基于MPI-GPU[5-6]方案已经成为设计GPU高性能计算集群的主流选择.

大数据时代，传统的并行计算框架已经无法满足快速、高效率的要求，Spark[7]正是在这种背景下诞生的产物.与CPU相比，GPU有众多线程，在进行数据密集型计算时具有速度优势.Spark-GPU能够同时使用CPU和GPU两种计算资源，其性能高于以CPU为计算核心的Spark，因此利用GPU加速Spark应用具有良好的应用前景.

1 研究现状

当前两大通用GPU计算框架分别是OpenCL(open computing language)和CUDA，实现Spark-GPU有两种方案：Spark+OpenCL和Spark+CUDA.第一种方案具有代表性的是Segal等[8]提出的SparkCL框架，该框架包括Aparapi和Spark两部分，其中由Aparapi执行Java编写的核函数并与OpenCL进行关联.与OpenCL相比，CUDA抽象程度高，易于使用，因此Spark-GPU多采用Spark+CUDA的方案.

在JVM(Java virtual machine)环境下，Spark通过JNI(Java native interface)与其他语言编写的GPU核函数进行交互从而实现了对GPU的调用，如文献[9-10].PyCUDA和Numba是CUDA支持Python所依赖的两个主要的工具.有较多的研究者使用PyCUDA作为在Spark与GPU结合的主要工具，如文献[11-12].

基于GPU加速的Spark框架能适合不同的计算任务，Li等[9]提出了针对机器学习的CPU-GPU异构Spark计算平台—HeteroSpark，Yuan等[10]利用GPU加速机器学习和SQL查询，其中机器学习的速度提升了16.13倍，SQL查询的速度提升了4.83倍.周情涛等[13]探讨了Spark-GPU的实现，并阐述了利用PyCUDA实现K-means算法的流程.Ohno等[14]利用GPU 加速Spark中部分行动(action)和转换(transform)操作.Serrano[11]在Spark框架中利用GPU加速重建高分辨率三维图像(像素分辨率为2 048×2 048×2 048)，大幅缩短了重建时间.

利用GPU加速CT图像重建是提高速度的重要方法，在高分辨率图像重建[15]和批量图像重建上都有良好的应用.本文在此研究基础上将Spark-GPU用于批量CT图像重建，通过加速滤波反投影(fltered backprojection,FBP)算法和同时代数迭代重建技术(simultaneous algebraic reconstruction technique,SART)算法，使图像的重建速度得到显著提升.

2 Spark-GPU并行计算框架

2.1 Spark框架

Spark是加州大学伯克利分校的AMP实验室(UC Berkeley AMP lab)所开源的通用并行计算框架，在大数据计算领域有着广泛的应用，如网络流量分析，用户行为预测，广告推荐等业务上.Spark的通用性和易用性使其在生物医学领域也有着良好的应用前景，如人类基因组测序云平台[16]、基因组学大数据分析[17]、医学图像并行处理[18]、批量医学图像重建[19]等.

Spark运行架构包括运行在主节点(master)上的集群资源管理器(cluster manager)、运行作业任务的工作节点或从节点(worker node、slave node)、每个应用的任务控制节点(driver)和每个工作节点上负责具体任务的执行进程(executor).其中，集群资源管理器可以是Spark自带的资源管理器，也可以是YARN(yet another resource negotiator)或Mesos等资源管理框架.Spark定义了两大类方法：转换(transform)和行动(action).转换操作是惰性的，即通过构造有向无环图记录弹性分布式数据集(resilient distributed data set,RDD)间的依赖关系，在行动操作时才会进行RDD变换.首先由SparkContext构建有向无环图(directed acyclic graph,DAG)，作业调度模块(DAG scheduler)将DAG图分解成阶段(stage)，任务调度模块(TaskScheduler)则将stage分解为task.任务控制节点向集群管理器申请资源，启动Executor并向其发送应用程序代码和文件，由Executor执行任务并将最终结果返回给任务控制节点.

2.2 Spark-GPU架构

Spark-GPU总体架构如图 1所示，在Spark集群的基础上，每个节点均结合了CPU和GPU两种计算单元.当一个节点内有多个GPU时，为了充分发挥多GPU的性能，需要对GPU进行调度以选择最合适的GPU进行计算.但在本架构中，一个节点只配置一个GPU故无须调度.主节点将任务分发给各工作节点，各节点从本机硬盘读取数据，在节点内部由CPU负责复杂的串行逻辑计算，GPU负责数据密集的并行计算，通过两者结合能充分提升集群性能.

图1 Spark-GPU并行架构

Spark-GPU程序可看作单机程序和Spark程序的结合，其中Spark程序负责读取数据和保存数据，而单机程序是进行计算的部分.本文设计的程序分为3个部分：①主程序使用thunder读取数据并创建RDD， RDD中的每个元素为一个完整的投影数据.在创建RDD时可以设置n个元素为1个task，默认task数等于元素数.当n较大时，执行一个task会消耗更多内存.主程序还负责将各节点返回的结果保存到本地. ②图像重建程序负责接收RDD中的每个元素并进行重建，通过网络将重建结果返回给Master. ③图像重建程序将部分数据传到GPU显存，调用GPU加速接口执行计算任务.主节点负责第一部分并不参与计算，但仍需要配置显卡，从节点负责执行①和②.

2.3 开源工具Numba

早期CUDA支持C、C++和Fortran语言，随着Python在科学计算领域的突出表现，NVIDIA公司先后推出PyCUDA[20]和Numba[21]两款工具，使Python成为CUDA支持的第4种语言.PyCUDA是支持RTCC(GPU run-time code generation)技术的Python开源工具.其原理是将C语言编写的核函数源码当作字符串传给Source Module的构造函数，通过动态编译生成机器代码.PyCUDA采用了编译器缓存机制，在重复调用核函数时能减少编译次数，节省编译时间.Numba是基于底层虚拟机(low level virtual machine，LLVM)的即时编译工具.Numba提供和CUDA C类似的编程模型，称为Numba CUDA模型，该模型包含了CUDA C的大部分功能，并且完全使用Python语言，支持在核函数内部使用Cmath、 Math函数库，因此Numba能充分发挥Python的灵活性和扩展性.Spark并不支持直接调用GPU，通过导入Numba工具库可以直接用Python编写核函数，从而在Spark中使用GPU进行计算.

3 利用GPU加速CT图像重建分析

3.1 FBP 算法及并行化

CT重建算法按照其重建方式可分为两类：解析重建算法和迭代重建算法.FBP是经典的解析重建算法.根据中心切片定理，二维图像可表达为：

t=xcosθ+ysinθ

(1)

f(x,y)为x-y平面上定义的密度函数，P(ω,θ)是密度函数的二维傅里叶变换与ωy轴夹角为θ的切片,|ω|是斜坡滤波器.

FBP算法的步骤可分为滤波和反投影两步.滤波的时间复杂度为O(N2log2N),反投影的时间复杂度为O(N3)，可见反投影计算量大、耗时长.GPU滤波一般采用cuFFT库，时间复杂度为O(Nlog2N).对短序列进行快速傅里叶变换(fast Fourier transform,FFT)FFT变换，GPU的速度要慢于CPU的速度.

GPU可将所有角度的反投影一步完成，时间复杂度减少为O(1)，可见将反投影并行化是加速FBP重建的关键.主机将每个角度滤波后的投影数据传输到GPU全局内存，由GPU进行反投影并将各个角度的投影数据叠加，再将最终数据传输到主机内存.例如180个角度的投影数据，每个角度的数据长度为725，将投影数据进行FFT变换和反投影的时间对比如图 2.

图2 FFT和反投影时间对比

由图2可以看出在CPU端进行FFT的时间远小于在GPU端进行FFT的时间，由于GPU内核在首次启动时耗时更长，因此在CPU端进行滤波是最优的选择.反投影计算量大，因此GPU多线程更有优势，利用GPU进行反投影的时间约为CPU的1/3.通过分析可知在主机端进行滤波，在GPU端进行反投影是最佳的组合.

3.2 SART算法及并行化

SART算法是将原图像看作待求解的矩阵，并通过方程来求解.设原图像大小是像素，如式(2)所示：

WX=P

(2)

X即为待求的N个像素值(N=n2)，是一维列向量.P是所有射线投影数据向量，P=(p1,p2,…,pM),M是投影射线数量.W是M×N维矩阵，其元素wij表示第j个体素对第i个投影数据pi的贡献值.

同时代数迭代算法(SART)具体步骤如下：

1)对未知图像进行赋初值：

(3)

2)计算图像的投影值：

(4)

3)计算误差：

(5)

4)计算修正值：

(6)

5)对第j个像素值进行修正：

(7)

λ为松弛因子.

6)以上迭代结果作初始值，重复2)-5)步骤.

SART可以简化为4个步骤，即正投影、修正、反投影、更新，基于GPU加速的SART算法多采用体素驱动正投影，并行遍历以及射线驱动反投影[22].正投影和反投影总时间复杂度均为O(N3)，并行化后时间总的复杂度变为2O(N)，复杂度降低，速度提升大.FBP中所有角度的反投影是同时进行的，而SART的反投影是逐次进行，一次子迭代过程包含4个步骤，只对一个角度的数据进行反投影，算法的差异导致SART的速度远慢于FBP的速度.

根据两种图像重建算法共同的特点，本文设计了两种利用GPU加速的函数：正投影函数和反投影函数.将数据密集型计算任务转移到GPU，能够充分利用了GPU多线程，大幅提升程序运行速度.针对SART迭代计算的特点，使用GPU进行修正和更新，通过将部分数据常驻GPU内存，在调用核函数时能减少主机和GPU的数据传输次数，加快程序运行速度.实验表明使用GPU缓存能够使SART程序运行时间减少32%.

4 实验

4.1 实验环境配置

集群硬件配置如表 1所示.集群由一台PC作为Master节点，3台PC作为Slave节点，各节点通过百兆交换机连接.每个节点均运行Centos 7.0系统，系统所需环境为Spark-2.3，Numba-36.0，thunder-1.4.2，CUDA-9.1.

表1 集群硬件配置Table 1 Cluster hardware configuration

本文所用的CT图像数据为Shepp-Logan模型，图像像素分辨率分别为 256×256、512×512.投影角度0°～180°，共180个角度，每个角度投影下的射线数分别为363和725.SART算法松弛因子为0.02，迭代次数为2×180次.重建图像效果如图3.

图3 原图像、FBP重建图像和SART重建图像

4.2 实验结果

本实验一次性重建图像数量分别为30，60，90，300，600，900，1 200张，图像数量为3的整数倍，便于主节点分配.记录3个节点下重建上述数量图像的时间，详细数据如表 2所示，重建时间以秒为单位，数据保留小数点后两位.

表2 重建图像所需的时间Table 2 Time required to reconstruct the image

使用SART和FBP算法重建一张图像的时间随数量的变化如图4.可以看出，随着图像数量的增加，平均重建一张图像所需的时间逐渐减少，速度越来越快.这是因为从提交Spark应用开始，任务控制节点向集群管理器申请资源，任务调度，启动executor以及集群间的通信等都需要一定的时间，这些时间都是系统开销.在图像数量较少时，系统开销在总时间中占比大，随着任务执行时间变长，系统开销相对减少，重建速度逐渐提升并趋于稳定.FBP算法的时间最快为0.03 s，SART算法的时间最快为0.85 s.同样参数下，在单机上用CPU运行FBP算法的时间为1.31秒，运行SART算法的时间为8.34 s，经过Spark-GPU加速后，FBP算法的速度有近40倍的提升，SART算法的速度有近10倍的提升.由以上分析可知，随着数据增多，Spark-GPU加速效果明显，速度逐渐达到理想状态，适合大规模计算任务.

4.3 集群性能分析

通过计算相对加速比[23]可以判断集群的性能和效率，其定义式为：

相对加速比反映了当集群节点数目增加时集群性能提升的大小，在理想情况下，加速比与节点数量成正比.在重建900张分辨率为512×512的图像时，集群的加速比如图5所示，可以看出加速比与节点数量并不严格成正比，这是由于各节点的CPU性能不完全相同导致加速比降低.SART算法的运行时间远高于FBP算法，系统开销相对减少，所以SART算法的加速比大于FBP算法的加速比.使用性能相同的处理器、保证硬件配置一致，防止集群出现性能短板，合理地分配资源和设置任务大小等方法都可以提升加速比.

图5 集群加速比

以Slave03节点为例，图 6记录该节点运行的前200 s内CPU和内存利用率，可以看出在程序开始后CPU的使用率接近100%，内存使用率在30%上下浮动.在默认情况下，Spark给每个节点分配的任务数与该节点CPU的核数相同，这使得CPU性能被充分利用.例如某节点的CPU为四核，Spark可以发起4个线程同时执行4个Task.Spark默认的任务调度策略为FIFO(first in first out)，一个task所需的内存大小与该task的分区大小有关，当一个task完成后，结果会立即返回给主节点并释放内存，从而避免内存占用过高.

与从节点相比，由于主节点并不参与计算，而是接收各个节点传回的数据，因此对网络带宽的要求高，图 7记录了FBP和SART程序开始后的45 s内主节点网络下行速度，可以看出FBP程序基本将网络带宽全部占用，这是因为FBP算法的速度更快，短时间内大量的数据要传回主节点.SART的速度较慢，从节点有充分的时间传输数据，因此主节点的下载速度较低，只出现有规律的波动.网络带宽对Spark任务调度亦有较大影响，由于FIFO策略的特点，后续任务等待被执行直至从节点将结果传给主节点，当网络带宽较小时，任务调度时间会明显变长.在FBP程序中，任务调度时间和执行时间相当导致计算效率严重降低.通过提升带宽能缩短传输时间，或者将数据直接保存在从节点的磁盘中，后者能将总时间缩短50%.

主节点接收传回的数据，因此内存占用大于从节点.Spark对内存的依赖使得内存大小直接影响了能够处理数据的多少，为了能处理更多的数据，可以将数据集拆分，分批创建RDD，待前一个RDD任务完成后再执行下一个RDD任务.除了集群硬件性能外，Spark系统配置参数的设置也会对任务的运行性能有很大的影响[24]，如maxResultSize太小会导致driver没有足够空间容纳序列化结果，太大会导致driver内存溢出.

图6 slave03系统资源使用情况

图7 Master节点下行速度

5 讨论

大规模图像重建对计算机性能有了更高的要求，本文构建由4节点组成的Spark-GPU集群，实现CPU和GPU结合的并行计算，大幅提升了图像重建的速度，并可通过增加节点的方式获得更大的加速比以满足更高的任务要求.批量图像重建是将数据由大化小，分批重建，其核心思想是MapReduce，在超高分辨率图像和三维图像的加速重建中也有应用价值.本文使用Spark-CUDA-Numba-Thunder的结合模式是为整合Spark与GPU提供了新方向.

由于各节点通过百兆交换机相连，不能满足大数据量的传输要求，网络配置还需进一步升级，同时设计新的调度策略提升数据的传输效率，降低调度时间.GPU内存一般不及主机内存大且不可扩展，如何优化算法减少GPU内存使用，对充分发挥GPU性能具有重要意义.Spark还无法在应用间管理和调度GPU资源，在应用内对GPU的使用和单机程序并无区别，当一个节点内有多个GPU时，对GPU的管理调度十分重要.有文献提出将GPU作为单独的计算资源进行调度，无论在应用间和应用内都能自动选择最合适的GPU[10-11]，但需要额外设计GPU的调度程序.如何将GPU融合到Spark框架中，不需要外部工具包并在多语言平台实现统一是重要的研究方向.