CN102609978A - 基于cuda架构的gpu加速锥束ct图像重建的方法 - Google Patents

Abstract

本发明涉及一种基于CUDA架构的GPU加速锥束CT图像重建的方法；具体为：在CUDA编程模型中，采用GPU进行锥束CT图像重建算法的运算，锥束CT图像重建采用FDK算法，构建CPU和GPU协作的CPU-GPU异构计算模式，将算法划分为M个子任务，每个子任务的结构包括两个部分：数据和在数据上施加的运算，根据各个子任务的特点将其调度到不同的硬件上执行，子任务分为串行子任务和并行子任务，其中，CPU负责子任务的调度和串行子任务的执行，GPU负责并行子任务的执行，CPU与GPU之间通过PCI-E总线进行通信；本发明提供了一种计算速度快的基于CUDA架构的GPU加速锥束CT图像重建的方法。

Description

基于CUDA架构的GPU加速锥束CT图像重建的方法

(一)、技术领域：本发明涉及一种CT图像重建的方法，特别是涉及一种基于CUDA架构的GPU加速锥束CT图像重建的方法。

(二)、背景技术：锥束CT图像重建加速方法简介：

锥束CT具有扫描速度快、空间分辨率高和射线利用率高的优点，在实际应用中已经越来越广泛的被使用。锥束CT三维重建算法计算量大、耗时长，比如最常用的FDK算法在没有任何加速的情况下重建512³规模的体数据大约需要90分钟，很难满足实际应用需求，因此如何提高重建速度是锥束CT亟待解决的问题。根据实现方式上是否脱离了单个CPU来看，目前三维重建算法的加速可以分为软件加速和硬件加速。软件加速通过对算法本身的改进来实现，比如利用参数表法减少每一次反投影过程中的计算量、利用对称性并结合递归技术降低计算复杂度，或是通过对处理器的更好利用来达到加速目的，但是，软件加速受限于硬件，加速效果有限。硬件加速主要是通过多CPU并行机、FPGA(Field Programmable Logic Array)、Cell处理器以及GPU，在这些硬件当中，基于PC的GPU具有运算快捷、价格低廉和更新换代快等优点，并且锥束CT的数据扫描过程与计算机图形学的投影过程本质上是一致的，因此GPU非常适合用于锥束CT图像重建的加速。

FDK算法简介：

FDK算法是Feldkamp、Davis和Kress三人在1984年针对锥束几何圆形扫描轨迹提出的一种近似三维图像重建算法。FDK算法简洁、高效、在小锥角下具有较好的重建效果，自提出以来就得到了广泛的使用。迄今为止，FDK算法仍是在完整数据条件下国际上公认的最有效的三维滤波反投影重建算法。

FDK算法主要包括3个步骤：加权、滤波和反投影。具体实现如下：

I、对投影数据做加权运算：

$P_{\mathrm{\θ}}^{\′}(u,v)=\frac{d}{\sqrt{d^2+u^2+v^2}}{P}_{\mathrm{\θ}}(u,v)---\left(1\right)$

其中，d为光源到旋转中心的距离，(u，v)为探测器平面坐标，P_θ(u，v)表示第θ角度的投影数据，P′_θ(u，v)表示对投影数据加权后的结果。

II、对加权后的投影数据做滤波运算：

P_θ ^*(u，v)＝P′_θ(u，v)*h(u)    (2)

其中，函数h(u)表示|ω|的傅里叶逆变换，P_θ ^*(u，v)表示对加权后的投影数据进行逐行滤波后的结果。

III、反投影运算：

$f(x,y,z)={\∫}_0^{2\mathrm{\π}}\frac{d^2}{{(d-x\sin \mathrm{\θ}+y\cos \mathrm{\θ})}^2}{P}_{\mathrm{\θ}}^{*}\left(u(x,y,\mathrm{\θ})v(x,y,z,\mathrm{\θ})\right)\mathrm{d\θ}$

$u(x,y,\mathrm{\θ})=\frac{d\×(x\cos \mathrm{\θ}+y\sin \mathrm{\θ})}{d-x\sin \mathrm{\θ}+y\cos \mathrm{\θ}}$

$v(x,y,z,\mathrm{\θ})=\frac{d\×z}{d-x\sin \mathrm{\θ}+y\cos \mathrm{\θ}}---\left(3\right)$

其中，f(x，y，z)表示被重建点，u(x，y，θ)和v(x，y，z，θ)表示点(x，y，z)投影到探测器上的坐标。

FDK算法的算法复杂度为n×o(N³)，n为投影数，N为重建规模的维长，当N较大时，FDK算法的计算量迅速增加。在FDK算法中反投影部分占据了绝大部分时间，由公式(3)可知，反投影是基于每个体素的，不同体素间的运算互不相关，因此反投影运算具有很好的并行性。

GPU开发方式简介：

1.可编程图形流水线

早期的GPU只能通过汇编语言进行开发，开发效率低是其一大缺点。直到2002年，Nvidia公司与微软密切协作，开发了图形硬件的第一个高级语言Cg(Cfor Graphics)。Cg为开发人员提供了一整套编程平台，并且完全与OpenGL或Direct3D结合在一起，跨平台性大大提升，是图形处理器程序开发的一个里程碑。可编程图形流水线的开发语言还有微软的高级着色语言(High Level ShadingLanguage，HLSL)和OpenGL架构评估委员会(Architecture Review Board，ARB)开发的OpenGL着色语言(OpenGL Shading Language，GLSL)。

尽管这些高级语言将GPU带入了通用计算的时代，但它们还无法离开图形API(如OpenGL)，数据只能以纹理存储，遵从固定的图形流水线(如图1所示)，因此对于非图形学专业的人来说还有一定难度。

2.CUDA编程模型

针对可编程图像流水线的缺点，Nvidia公司于2007年推出了CUDA(Compute Unified Device Architecture)编程模型，降低了开发门槛，促进了GPU通用计算的发展。CUDA编程模型具有类C的编程风格，开发人员无需重新学习新语法和掌握晦涩的图形API，只需要对并行计算有一定了解就能够获得性能上的较多提升。然而，CUDA在带来自由与灵活的同时，也要求开发人员必须掌握其并行机制与GPU架构方面的知识才能开发出高性能的并行程序，主要包括GPU的硬件架构、CUDA软件体系、存储器模型及线程模型，下面将对其进行详细介绍。

(1)支持CUDA的GPU硬件架构

支持CUDA的GPU较以往的GPU有两点不同：一是它采用了统一处理架构，即不再区分顶点着色器和片段着色器；二是引入了片内共享存储器，支持线程间通信和随机存储。例如，Tesla系列GPU的核心GT200，其硬件由流处理器阵列(Scalable Streaming Processor Array，SPA)和存储器系统组成。SPA的结构可以分为两层：第一层是线程处理器群(Thread Processing Cluster，TPC)；第二层是流多处理器(Streaming Multiprocessor，SM)。多个SM组成1个TPC，其典型结构如图2所示。其中，每个SM拥有多个SP，单个SP有独立的寄存器和指令指针，相当于CPU中的一条流水线。

(2)CUDA存储器模型

在CUDA编程模型中，CPU与GPU各自拥有相互独立的存储器空间，通过PCI-E总线实现数据的交互。

支持CUDA编程模型的GPU的存储器主要有寄存器(register)、共享存储器(shared memory)、常数存储器(constant memory)、纹理存储器(texturememory)、本地存储器(local memory)和全局存储器(global memory)，如图3所示。其中寄存器和共享存储器位于GPU芯片内部，访问速度最快；常数存储器和纹理存储器虽然位于显存中，但由于它们可以被缓存，访问速度比直接读取存储器快；本地存储器和全局存储器位于显存中，且没有缓存，所以访问速度最慢，一般存在几百个时钟周期的延迟。

(3)CUDA软件体系

CUDA的软件体系由CUDA库函数、CUDA运行时API和CUDA驱动API三部分构成，如图4所示。CUDA的Nvcc编译器在编译程序时首先将源文件中的主机端和设备端代码分离开来，设备端代码由Nvcc编译成ptx代码或二进制代码，主机端代码由其它合适的高性能编译器进行编译。

目前CUDA函数库主要有CUFFT、CUBLAS和CUDPP三个函数库。CUFFT(CUDA Fast Fourier Transform)函数库是一个利用GPU进行傅立叶变换的函数库，CUBLAS(CUDA Basis Linear Algebra Subprograms)函数库是一个基本的矩阵和向量的运算库，CUDPP(CUDA Data Parallel Primitives)函数库提供了很多基本的常用的并行操作函数，如排序、搜索等。使用这些函数库中的函数，就可以方便地解决实际中的计算问题。

(4)CUDA线程模型

在CUDA模型下，运行在GPU上的程序称为Kernel(内核)函数，Kernel函数以Grid(线程网格)的形式组织，每个Grid由若干Block(线程块)组成，每个Block又包含若干Thread(线程)。在程序执行时，Grid被加载到SPA上，Block则被分发到各个SM，Thread以Warp(线程束，包含32个线程)为单位运行。Block可以是1维、2维或3维，所有的Block又以1维或2维的形式组成Grid，一个Block内的线程可彼此协作，通过共享存储器来共享数据，并同步其执行来协同存储器访问。通过Block的粗粒度分解和Thread的细粒度并行执行就可以完成任务的并行计算，达到加速的目的。

(三)、发明内容：

本发明要解决的技术问题是：克服现有技术的缺陷，提供一种计算速度快的基于CUDA架构的GPU加速锥束CT图像重建的方法。

本发明的技术方案：

一种基于CUDA架构的GPU加速锥束CT图像重建的方法，在CUDA编程模型中，采用GPU进行锥束CT图像重建算法的运算，锥束CT图像重建算法采用FDK算法，构建CPU和GPU协作的CPU-GPU异构计算模式，CPU被称为主机，GPU被称为设备；在CPU-GPU异构计算模式中，利用GPU进行锥束CT图像重建算法的运算时，将算法划分为M个子任务，M为大于等于1的自然数，每个子任务的结构包括两个部分：数据和在数据上施加的运算，根据各个子任务的特点将其调度到不同的硬件上执行，子任务分为串行子任务和并行子任务，其中，CPU负责子任务的调度和串行子任务的执行，GPU负责并行子任务的执行，CPU与GPU之间通过PCI-E总线进行通信，如图5所示。

尽管GPU在解决计算密集型问题上比CPU具有较大优势，然而不可能仅仅依靠GPU来完成计算任务，很多串行指令和控制指令仍然需要CPU来完成。因此，基于CUDA编程模型的GPU通用计算实际上是充分发挥CPU与GPU在解决特定问题中各自的优势，构建CPU和GPU协作的CPU-GPU异构计算模式来更好地解决计算问题。

CPU-GPU异构计算模式有两个特点：一是CPU与GPU的执行是异步的，即CPU将任务分配给GPU以后，它将继续执行下面的任务而不会等待GPU端的运算完成，直到执行同步指令。这种模式虽然需要人为地对CPU与GPU的运算进行同步，但却提高了运算效率。二是通信开销量较大。在执行任务的过程中，CPU将数据通过PCI-E总线发送给GPU，GPU在完成运算之后再将结果回传到内存中，这模式将增加大数据量应用场合的执行时间。

FDK算法中含有加权步骤和滤波步骤，采用下列方式对加权步骤和滤波步骤进行加速：

加权步骤和滤波步骤所用时间占整个FDK算法消耗时间的比例较少，图6为数据规模为256³和512³时加权步骤与滤波步骤分别占总重建时间的比例，并且随着重建规模的增大，所占比例还会继续降低。加权步骤的计算量小，对整个FDK算法运行时间的影响很小，因此，加权步骤选择在CPU或GPU中完成，在本发明中，加权步骤选择在GPU中完成。

利用CUDA编程模型中自带的CUFFT库对投影数据做频域滤波，滤波时还利用了FFT函数能够同时处理一批一维傅里叶变换的特性，并且根据复数FFT变换的特点，一次对投影数据的两行做滤波。

目前，对于FDK算法中反投影部分的CUDA加速策略主要有：

1、在线程分配时考虑反投参数的Z轴相关性。每个线程各自完成对应Z轴上一系列点的重建任务，这种线程分配方式可以减少GPU的重复计算量；

2、纹理存储器的使用。纹理存储器能够被缓存，因此使用纹理存储器可以大大提高访问速度，并且纹理存储器带有优化过的寻址系统，支持方便高效的二维插值操作，可用来完成反投影中的双线性插值运算；

3、对全局存储器采用合并访问优化。在满足合并访问条件的情况下，只需要一次传输便可以处理来自一个Half-Warp的16个线程对全局内存的访问请求；

4、对三角函数计算的优化。在GPU中完成一次三角函数运算需要32个时钟周期，而一次乘加只需要4个时钟周期，因此在计算三角函数时采用递归技术，将不同角度下的三角函数运算转换为一次乘加运算，便可在一定程度上减少GPU的运算时间。

由于GPU架构本身在高性能计算方面的优势，在应用上述优化策略后，便可获得相比CPU上百倍的加速比。

本发明在上述并行策略的基础上，做了更加细致的优化，使得优化后的加速效果获得更大幅度的提升。相比上述优化策略，在本发明中，对FDK算法中反投影部分的加速有以下列三个策略：

策略1、线程分配方式的优化；

策略2、通过使用常数存储器存储计算反投影参数时的中间量，减少GPU的重复计算量；

策略3、通过优化一个内核函数(Kernel)中反投影的张数，减少访问全局存储器的次数。

策略1具体为：

GPU的任务划分可以按照输入划分或者按输出划分；对于反投影算法来说，输入为投影数据，输出为重建体数据，GPU的两种任务分配方式实质上反映了两种不同的反投影实现方式：基于射线驱动和基于体素驱动；基于射线驱动的方式按投影数据进行任务划分，一个或几个线程完成一条射线上所有体素的反投，反投时首先计算出当前射线穿过的所有体素，然后将这些体素的值赋为当前射线投影值，由于不同射线对应的线程或同一射线对应的不同线程可能同时给一个体素赋值，即这种方式存在“写竞争”问题；基于体素驱动的方式按体数据划分任务，一个线程完成一个或几个体素的反投，反投时首先计算出当前体素的正投影位置，然后取该位置的投影值赋给当前体素；基于体素驱动的反投影方式不存在“写竞争”问题，不需要设计额外的归约步骤，因此这种任务分配方式更适合于GPU加速。

本发明采用基于体素的驱动方式，即按输出对GPU的任务进行划分；在线程分配时需要考虑一个很重要的原则，即SM的占用率不能太低，SM占用率指的是每个SM上的活动Warp数量与允许的最大活动Warp数量之比；这主要是因为GPU是通过线程间的切换来隐藏长延时操作(访问全局存储器、线程间同步等)的，当一个Warp中的线程进行长延时操作时，另外一个活动Warp中的线程会自动进入运算状态，这样便可以隐藏一部分延迟；但这并不代表SM占用率越高越好，SM占用率越高则每个线程占用的GPU资源越少，即每个线程完成的计算量越少，而每个SM上的最大活动Warp数量是一定的，因此就可能出现即使SM中活动的Warp数量达到了最大，由于Warp中每个线程计算量太少，从而使所有活动Warp中线程同时进入长延时操作，不能充分的隐藏延迟；因此，在GPU的流多处理器(SM)占用率和每个线程完成的计算量间选择一个平衡点，使GPU的性能发挥到最佳。

策略1中的线程分配方式为：当体数据的规模为N³时，则线程块(Block)的大小恒定为(16，16，2)，线程网格(Grid)的尺寸随N的改变而变化为(N/16，N/16，1)。

策略2具体为：

反投影需要计算体数据中任意一点(x，y，z)在投影角度为θ时投射到探测器上的点(u，v)，在这个运算过程中需要多次计算只与投影角度和几何关系有关的三角函数和其他中间变量；对于每个体素，这些相同的变量都被计算一次，假设体数据规模为N³，则这些值会被重复计算N³次，造成系统资源的极大浪费；针对这个问题，本发明将反投影运算中与体素(x，y，z)无关的变量进行分离和合并，并在反投影前将这些变量计算出来存于GPU的常数存储器中，反投影运算时，直接读取常数存储器中的变量参与计算。具体过程如下：

投影点(u，v)的计算为：

$w=\frac{((x-\mathrm{vCenter})\×\mathrm{pix}\×(-\sin \left(\mathrm{\θ}\right))+(y-\mathrm{vCenter})\×\mathrm{pix}\×\left(\cos \left(\mathrm{\θ}\right)\right)+s)}{s}$

$u=\frac{((x-\mathrm{vCenter})\×\cos \left(\mathrm{\θ}\right)+(y-\mathrm{vCenter})\sin \left(\mathrm{\θ}\right)}{w}+\mathrm{pCenter}$

$v=\frac{z-\mathrm{vCenter}}{w}+\mathrm{pCenter}---\left(4\right)$

分离和合并变量后为：

w＝x×A[0]+y×A[1]+A[2]

$u=\frac{x\×A\left[3\right]+y\×A\left[4\right]+A\left[5\right]}{w}+\mathrm{pCenter}$

$v=\frac{z-\mathrm{vCenter}}{w}+\mathrm{pCenter}.---\left(5\right)$

其中，vCenter表示体数据中心，pCenter为投影数据中心，pix为体素大小，θ为投影角度。

由式(5)可知，分离和合并变量后每一个角度的反投影运算可以提取出6个与体素(x，y，z)无关的中间变量，假设投影数为360，则整个反投影共有2160个变量(单精度浮点型)需要在反投影前计算出来并存于GPU常数存储器中；常数存储器是GPU特有的只读存储空间，能够被缓存，并且来自同一half-Warp的线程访问常数存储器中的同一数据时，如果发生缓存命中，只需要一个周期就可以获得数据；一般来说，常数存储器空间较小(如Tesla C1060中只有64KB)，但完全能够满足本文方法的需要；反投影时，只需要读取常数存储器中的变量值参与常规的乘加运算便可得到最终的反投影参数；因此，这种加速策略不仅可以避免用GPU去计算开销极大的三角函数，而且避免了GPU的重复计算，在提升反投影运算效率方面具有较好效果。

策略3具体为：

重建体数据通常存放于GPU全局存储器内，全局存储器占据了显存绝大部分，可以用来存放大规模的数据，但全局存储器没有缓存加速，虽然合并访问方式可以极大地提高访问速度，但通常仍然会有几百个周期的访问延迟；研究表明，GPU用于高性能计算的瓶颈不是计算消耗而是访存消耗；因此，如何缩减访问全局存储器的时间是GPU加速的关键。

通常情况下，一个内核函数(Kernel)中只对1幅投影图像进行运算，对于360幅投影数据，整个反投影过程需要读/写全局存储器360×N³次；在本发明中，一个内核函数(Kernel)完成m个角度投影图像的反投影，m为大于等于2的自然数，每个内核函数(Kernel)需要计算m个角度投影图像的反投影参数，但只读/写全局存储器N³次，即可以将读/写全局存储器的次数变为原来的1/m。

在减少全局存储器读写次数的同时，算法会增加内核函数(Kernel)中每个线程的计算负担，若一味的增大m，势必会降低整个GPU中活动Block和活动Warp的数量，活动Block和活动Warp数量减少又会反过来影响GPU隐藏长延时操作(访问全局存储器)的优势；所以，寻找一个适中的m显得尤为重要；在本发明在实验平台上通过多次试验发现，当一个内核函数(Kernel)中同时完成的反投影图像为6幅时，即m为6时，双方达到一个平衡，加速效果最为理想；图7为一个内核函数(Kernel)同时完成不同角度数投影图像运算时全部反投影消耗的时间。

本发明的有益效果：

1、本发明以CPU和GPU协作的CPU-GPU异构计算模式为平台，利用CUDA为编程模型，实现了对锥束CT的FDK重建算法的并行加速；完全优化后的并行加速方法，重建256³规模的数据需要0.5秒，重建512³规模的数据需要2.5秒，与单CPU的重建速度相比有上千倍的提高。

2、本发明中，支持CUDA的GPU在架构上有了显著改进，可以更加有效地利用过去分布在顶点渲染器和像素渲染器的计算资源，使GPU更加适合于通用计算。

(四)、附图说明：

图1为可编程图形流水线操作流程示意图；

图2为支持CUDA的GPU的硬件架构示意图；

图3为CUDA存储器模型示意图；

图4为CUDA软件体系示意图；

图5为CPU-GPU异构计算模式结构示意图；

图6为不同重建规模下加权步骤与滤波步骤占总重建时间比例的示意图；

图7为一个内核函数同时完成不同反投影角度数时反投影的消耗时间示意图；

图8为重建体数据在Z＝0时的切片图(数据规模256³)；

图9为对应图8的第128和256行数据剖线图；

图10为重建体数据在Z＝0时的切片图(数据规模5123)；

图11为对应图10的第128和256行数据剖线图。

(五)、具体实施方式：

参见图5～图11，图中，基于CUDA架构的GPU加速锥束CT图像重建的方法为：在CUDA编程模型中，采用GPU进行锥束CT图像重建算法的运算，锥束CT图像重建算法采用FDK算法，构建CPU和GPU协作的CPU-GPU异构计算模式，CPU被称为主机，GPU被称为设备；在CPU-GPU异构计算模式中，利用GPU进行锥束CT图像重建算法的运算时，将算法划分为M个子任务，M为大于等于1的自然数，每个子任务的结构包括两个部分：数据和在数据上施加的运算，根据各个子任务的特点将其调度到不同的硬件上执行，子任务分为串行子任务和并行子任务，其中，CPU负责子任务的调度和串行子任务的执行，GPU负责并行子任务的执行，CPU与GPU之间通过PCI-E总线进行通信，如图5所示。

尽管GPU在解决计算密集型问题上比CPU具有较大优势，然而不可能仅仅依靠GPU来完成计算任务，很多串行指令和控制指令仍然需要CPU来完成。因此，基于CUDA编程模型的GPU通用计算实际上是充分发挥CPU与GPU在解决特定问题中各自的优势，构建CPU和GPU协作的CPU-GPU异构计算模式来更好地解决计算问题。

CPU-GPU异构计算模式有两个特点：一是CPU与GPU的执行是异步的，即CPU将任务分配给GPU以后，它将继续执行下面的任务而不会等待GPU端的运算完成，直到执行同步指令。这种模式虽然需要人为地对CPU与GPU的运算进行同步，但却提高了运算效率。二是通信开销量较大。在执行任务的过程中，CPU将数据通过PCI-E总线发送给GPU，GPU在完成运算之后再将结果回传到内存中，这模式将增加大数据量应用场合的执行时间。

FDK算法中含有加权步骤和滤波步骤，采用下列方式对加权步骤和滤波步骤进行加速：

加权步骤和滤波步骤所用时间占整个FDK算法消耗时间的比例较少，图6为数据规模为256³和512³时加权步骤与滤波步骤分别占总重建时间的比例，并且随着重建规模的增大，所占比例还会继续降低。加权步骤的计算量小，对整个FDK算法运行时间的影响很小，因此，加权步骤选择在GPU中完成。

利用CUDA编程模型中自带的CUFFT库对投影数据做频域滤波，滤波时还利用了FFT函数能够同时处理一批一维傅里叶变换的特性，并且根据复数FFT变换的特点，一次对投影数据的两行做滤波。

目前，对于FDK算法中反投影部分的CUDA加速策略主要有：

1、在线程分配时考虑反投参数的Z轴相关性。每个线程各自完成对应Z轴上一系列点的重建任务，这种线程分配方式可以减少GPU的重复计算量；

2、纹理存储器的使用。纹理存储器能够被缓存，因此使用纹理存储器可以大大提高访问速度，并且纹理存储器带有优化过的寻址系统，支持方便高效的二维插值操作，可用来完成反投影中的双线性插值运算；

3、对全局存储器采用合并访问优化。在满足合并访问条件的情况下，只需要一次传输便可以处理来自一个Half-Warp的16个线程对全局内存的访问请求；

4、对三角函数计算的优化。在GPU中完成一次三角函数运算需要32个时钟周期，而一次乘加只需要4个时钟周期，因此在计算三角函数时采用递归技术，将不同角度下的三角函数运算转换为一次乘加运算，便可在一定程度上减少GPU的运算时间。

由于GPU架构本身在高性能计算方面的优势，在应用上述优化策略后，便可获得相比CPU上百倍的加速比。

本发明在上述并行策略的基础上，做了更加细致的优化，使得优化后的加速效果获得更大幅度的提升。相比上述优化策略，在本发明中，对FDK算法中反投影部分的加速有以下列三个策略：

策略1、线程分配方式的优化；

策略2、通过使用常数存储器存储计算反投影参数时的中间量，减少GPU的重复计算量；

策略3、通过优化一个内核函数(Kernel)中反投影的张数，减少访问全局存储器的次数。

策略1具体为：

GPU的任务划分可以按照输入划分或者按输出划分；对于反投影算法来说，输入为投影数据，输出为重建体数据，GPU的两种任务分配方式实质上反映了两种不同的反投影实现方式：基于射线驱动和基于体素驱动；基于射线驱动的方式按投影数据进行任务划分，一个或几个线程完成一条射线上所有体素的反投，反投时首先计算出当前射线穿过的所有体素，然后将这些体素的值赋为当前射线投影值，由于不同射线对应的线程或同一射线对应的不同线程可能同时给一个体素赋值，即这种方式存在“写竞争”问题；基于体素驱动的方式按体数据划分任务，一个线程完成一个或几个体素的反投，反投时首先计算出当前体素的正投影位置，然后取该位置的投影值赋给当前体素；基于体素驱动的反投影方式不存在“写竞争”问题，不需要设计额外的归约步骤，因此这种任务分配方式更适合于GPU加速。

本发明采用基于体素的驱动方式，即按输出对GPU的任务进行划分；在线程分配时需要考虑一个很重要的原则，即SM的占用率不能太低，SM占用率指的是每个SM上的活动Warp数量与允许的最大活动Warp数量之比；这主要是因为GPU是通过线程间的切换来隐藏长延时操作(访问全局存储器、线程间同步等)的，当一个Warp中的线程进行长延时操作时，另外一个活动Warp中的线程会自动进入运算状态，这样便可以隐藏一部分延迟；但这并不代表SM占用率越高越好，SM占用率越高则每个线程占用的GPU资源越少，即每个线程完成的计算量越少，而每个SM上的最大活动Warp数量是一定的，因此就可能出现即使SM中活动的Warp数量达到了最大，由于Warp中每个线程计算量太少，从而使所有活动Warp中线程同时进入长延时操作，不能充分的隐藏延迟；因此，在GPU的流多处理器(SM)占用率和每个线程完成的计算量间选择一个平衡点，使GPU的性能发挥到最佳。

策略1中的线程分配方式为：当体数据的规模为N³时，则线程块(Block)的大小恒定为(16，16，2)，线程网格(Grid)的尺寸随N的改变而变化为(N/16，N/16，1)。

策略2具体为：

反投影需要计算体数据中任意一点(x，y，z)在投影角度为θ时投射到探测器上的点(u，v)，在这个运算过程中需要多次计算只与投影角度和几何关系有关的三角函数和其他中间变量；对于每个体素，这些相同的变量都被计算一次，假设体数据规模为N³，则这些值会被重复计算N³次，造成系统资源的极大浪费；针对这个问题，本发明将反投影运算中与体素(x，y，z)无关的变量进行分离和合并，并在反投影前将这些变量计算出来存于GPU的常数存储器中，反投影运算时，直接读取常数存储器中的变量参与计算。具体过程如下：

投影点(u，v)的计算为：

$w=\frac{((x-\mathrm{vCenter})\×\mathrm{pix}\×(-\sin \left(\mathrm{\θ}\right))+(y-\mathrm{vCenter})\×\mathrm{pix}\×\left(\cos \left(\mathrm{\θ}\right)\right)+s)}{s}$

$u=\frac{((x-\mathrm{vCenter})\×\cos \left(\mathrm{\θ}\right)+(y-\mathrm{vCenter})\sin \left(\mathrm{\θ}\right)}{w}+\mathrm{pCenter}$

$v=\frac{z-\mathrm{vCenter}}{w}+\mathrm{pCenter}---\left(4\right)$

分离和合并变量后为：

w＝x×A[0]+y×A[1]+A[2]

$u=\frac{x\×A\left[3\right]+y\×A\left[4\right]+A\left[5\right]}{w}+\mathrm{pCenter}$

$v=\frac{z-\mathrm{vCenter}}{w}+\mathrm{pCenter}.---\left(5\right)$

其中，vCenter表示体数据中心，pCenter为投影数据中心，pix为体素大小，θ为投影角度。

由式(5)可知，分离和合并变量后每一个角度的反投影运算可以提取出6个与体素(x，y，z)无关的中间变量，假设投影数为360，则整个反投影共有2160个变量(单精度浮点型)需要在反投影前计算出来并存于GPU常数存储器中；常数存储器是GPU特有的只读存储空间，能够被缓存，并且来自同一half-Warp的线程访问常数存储器中的同一数据时，如果发生缓存命中，只需要一个周期就可以获得数据；一般来说，常数存储器空间较小(如Tesla C1060中只有64KB)，但完全能够满足本文方法的需要；反投影时，只需要读取常数存储器中的变量值参与常规的乘加运算便可得到最终的反投影参数；因此，这种加速策略不仅可以避免用GPU去计算开销极大的三角函数，而且避免了GPU的重复计算，在提升反投影运算效率方面具有较好效果。

策略3具体为：

重建体数据通常存放于GPU全局存储器内，全局存储器占据了显存绝大部分，可以用来存放大规模的数据，但全局存储器没有缓存加速，虽然合并访问方式可以极大地提高访问速度，但通常仍然会有几百个周期的访问延迟；研究表明，GPU用于高性能计算的瓶颈不是计算消耗而是访存消耗；因此，如何缩减访问全局存储器的时间是GPU加速的关键。

通常情况下，一个内核函数(Kernel)中只对1幅投影图像进行运算，对于360幅投影数据，整个反投影过程需要读/写全局存储器360×N³次；在本发明中，一个内核函数(Kernel)完成m个角度投影图像的反投影，m为大于等于2的自然数，每个内核函数(Kernel)需要计算m个角度投影图像的反投影参数，但只读/写全局存储器N³次，即可以将读/写全局存储器的次数变为原来的1/m。

在减少全局存储器读写次数的同时，算法会增加内核函数(Kernel)中每个线程的计算负担，若一味的增大m，势必会降低整个GPU中活动Block和活动Warp的数量，活动Block和活动Warp数量减少又会反过来影响GPU隐藏长延时操作(访问全局存储器)的优势；所以，寻找一个适中的m显得尤为重要；在本发明在实验平台上通过多次试验发现，当一个内核函数(Kernel)中同时完成的反投影图像为6幅时，即m为6时，双方达到一个平衡，加速效果最为理想；图7为一个内核函数(Kernel)同时完成不同角度数投影图像运算时全部反投影消耗的时间。

根据本发明提出的方法，进行如下实验：

测试体模采用3D Shepp-Logan体模，重建规模为256³和512³，投影图像规模分别为360×256³和360×512³。测试平台为：2.27GHz Intel Xeon双核CPU，16GB内存，Tesla C1060GPU；开发环境为：Visual Studio 2008，cuda 3.0runtime API。重建结果如图8～图11所示，重建时间(10次统计求均值)如表1所示，本发明方法反投影时间与单CPU反投影时间比较如表2所示。

重建结果只与原始数据做了对比，并没有和CPU重建结果作对比，这是因为我们更关心的是重建算法相对原始值得表现，而不是不同平台上表现的比较。由图9和图11中剖线图可知，GPU加速后的重建结果与原始值比较接近，与FDK算法在小锥角情况下的重建表现一致，即本发明针对FDK算法的加速策略并没有损失算法原有的重建精度。

表1本发明重建时间(单位：S)

表2反投影时间对比

Claims (7)

1.一种基于CUDA架构的GPU加速锥束CT图像重建的方法，在CUDA编程模型中，采用GPU进行锥束CT图像重建算法的运算，锥束CT图像重建算法采用FDK算法，其特征是：构建CPU和GPU协作的CPU-GPU异构计算模式，CPU被称为主机，GPU被称为设备；在CPU-GPU异构计算模式中，利用GPU进行锥束CT图像重建算法的运算时，将算法划分为M个子任务，M为大于等于1的自然数，每个子任务的结构包括两个部分：数据和在数据上施加的运算，根据各个子任务的特点将其调度到不同的硬件上执行，子任务分为串行子任务和并行子任务，其中，CPU负责子任务的调度和串行子任务的执行，GPU负责并行子任务的执行，CPU与GPU之间通过PCI-E总线进行通信；

FDK算法中含有加权步骤和滤波步骤，采用下列方式对加权步骤和滤波步骤进行加速：

加权步骤选择在CPU或GPU中完成；

利用CUDA编程模型中自带的CUFFT库对投影数据做频域滤波，滤波时还利用了FFT函数能够同时处理一批一维傅里叶变换的特性，并且根据复数FFT变换的特点，一次对投影数据的两行做滤波；

对FDK算法中反投影部分的加速有以下列三个策略：

策略1、线程分配方式的优化；

策略2、通过使用常数存储器存储计算反投影参数时的中间量，减少GPU的重复计算量；

策略3、通过优化一个内核函数中反投影的张数，减少访问全局存储器的次数。

2.根据权利要求1所述的基于CUDA架构的GPU加速锥束CT图像重建的方法，其特征是：所述策略1具体为：

采用基于体素的驱动方式，即按输出对GPU的任务进行划分；在GPU的流多处理器占用率和每个线程完成的计算量间选择一个平衡点，使GPU的性能发挥到最佳。

3.根据权利要求2所述的基于CUDA架构的GPU加速锥束CT图像重建的方法，其特征是：所述策略1中的线程分配方式为：当体数据的规模为N3时，则线程块的大小恒定为(16，16，2)，线程网格的尺寸随N的改变而变化为(N/16，N/16，1)。

4.根据权利要求1所述的基于CUDA架构的GPU加速锥束CT图像重建的方法，其特征是：所述策略2具体为：

将反投影运算中与体素无关的变量进行分离和合并，并在反投影前将这些变量计算出来存于GPU的常数存储器中，反投影运算时，直接读取常数存储器中的变量参与计算。

5.根据权利要求1所述的基于CUDA架构的GPU加速锥束CT图像重建的方法，其特征是：所述策略3具体为：

一个内核函数完成m个角度投影图像的反投影，m为大于等于2的自然数，每个内核函数需要计算m个角度投影图像的反投影参数。

6.根据权利要求5所述的基于CUDA架构的GPU加速锥束CT图像重建的方法，其特征是：所述m为6。

7.根据权利要求1所述的基于CUDA架构的GPU加速锥束CT图像重建的方法，其特征是：所述加权步骤选择在GPU中完成。

Images