基于CUDA架构的FDK算法的研究

基于CUDA架构的FDK算法的研究
骆岩红;王建华
【摘要】该文研究了一种利用GPU并行架构的CUDA来完成FDK三维图像重建算法的加速.分析了FDK三维图像重建算法的可并行性特点,设计了适合CUDA的并行方法,分别在算法加权、滤波和反投影3个阶段,实现FDK的加速.经过实验验证,该文提出的方法与算法与CPU单独实现图像重建相比,不仅获得了150倍以上的加速效果,并且两种方式完成的重建图像,质量接近,平均误差小于10-4.由此可得出结论,利用GPU的三维锥束CT图像重建能够得到较满意的结果.
【期刊名称】《工业仪表与自动化装置》
【年(卷),期】2015(000)006
【总页数】5页(P3-7)
【关键词】三维锥束CT;FDK算法;图像处理器;CUDA
【作者】骆岩红;王建华
【作者单位】西北民族大学电气工程学院,兰州730030;西北民族大学电气工程学院,兰州730030
【正文语种】中文
【中图分类】TP391
0 引言
锥束CT技术的迅速发展与重建算法的发展是分不开的，由于FDK算法实现方法
简单，所以得到了广泛的应用。

但是利用该算法完成图像的重建是十分耗时的，这样就极大地限制了它的实际应用及商业化生产。

因此，如何提高CT图像重建速度是一个亟待解决的问题，它也吸引了很多学者致力于这方面的研究。

加速FDK图像的重建，可以从两个方面实现:一方面可以从算法本身的改进着手，另一方面可以从基于硬件实现加速。

对此，国内外学者已提出了多种FDK重建加速方法，例如几何参数法、对称法、采用DSP的加速方法以及采用FPGA的加速方法等［1－3］。

随着图像处理器GPU硬件的加速发展，如何利用该硬件的可编程和并行性来实现FDK的复杂计算，已成为当前热点研究之一。

GPU是一种大规模并行计算的多核处理器体系结构，它拥有高度的并行性、强大的计算能力及超大的存储带宽。

很多学者利用它在图像重建的加速方面做了研究。

例如Cabral等人利用GPU纹理映射研究加速锥束 CT图像重建算法［3］，Fang Xu和Klaus Mueller等人也使用纹理映射对该方法进行了研究［1，4］，戴智晟等人提出了 T － FDK 重建算法的GPU快速实现方法［5］，马东平等人提出结合使用GPU多重纹理加速三维锥束CT的FDK图像重建［6］。

2006年，NVIDIA公司推出了CUDA，即基于并行编程模型和指令集的一种GPU 的通用计算架构。

CUDA利用GPU的并行计算能够比CPU更高效地解决许多复杂任务，而且使得它比传统的GPU具有更好的可编程性。

针对CUDA架构，Holger Scherl等人提出CUDA实现FDK重建算法的加速方法［7］，但是其算法复杂度较高，而且与CUDA架构的结合程度不高;重庆大学的王珏等人利用CUDA实现了FDK算法的加速［8］，取得了很好的效果。

中国民航大学孙毅刚等人在该方法的基础上，提出基于角度相关性的优化方法［9］，虽然这种方案提升了算法的性能，但是重建图像质量不高。

分析以上各位提出FDK重建算法的加速实现，发现他们在FDK滤波实现部分几乎采用了相同的方法，即直接利用CUFFT库中的FFT函数对投影数据进行滤波。

由于CUFFT内部结构是封闭性的，
所以只能被动的使用它，而且CUFFT没有专门优化数据的传输，使其性能受限。

该文分析了FDK三维重建算法可并行的特点，研究运用GPU技术加速FDK算法。

在FDK算法的加权处理、滤波和反投影3个部分，分别设计了适合CUDA的并行计算方法。

1 GPU并行计算理论和CUDA软件架构
1.1 GPU并行计算理论
早期计算机处理图像和视频等方面的运算时，都是通过CPU来完成的。

但是进入
3D时代后，人们需要处理实时的、高清晰度的3D图形，如果只靠CPU的运算能力，很难达到对庞大的3D图形数据的处理，在这样情况下，需要研制一种能够处理海量数据和计算量的芯片，以便满足快速处理大量图形计算的需求。

由此，出现了一种专门负责图形计算的处理器GPU(Graphics Processing Unit)。

GPU相对于CPU而言，它拥有更多核可以并行计算，在GPU内部，80%晶体管用作运算，剩余的拥有流控和数据缓存，从而决定了GPU具有更高密度的计算能力;而且GPU在内存宽度上(超过CPU带宽的10倍以上)也更适合大规模高性能的并行运算;现代GPU都是用户可编程的，用户可以通过编程来确定GPU的功能。

CUDA是NVIDIA公司为满足自家的GPU编写的一套编译器及相关的库函数。

它作为一种新型的软硬件架构平台，可以直接调用GPU的计算资源，而且增强了GPU可编程性，也大幅提高了其计算性能。

使用CUDA编程模型，可以在不了解GPU硬件架构的基础上，专注算法的编程，从而提高了开发效率。

1.2 CUDA软件架构介绍
CUDA架构的平台，由于可以直接调用GPU的计算资源，因此成为高性能的软件开发环境。

在CUDA框架中，CPU称为主机(Host)，GPU用来完成高度并行线程化的处理任务，故被作为主机的协处理器或设备(Device)。

CUDA通过允许编程人员定义内核(Kernel)，即扩展了C的C函数，当内核被调用时，将会被N个
CUDA线程执行N次，达到高度并行的效果。

在CUDA执行时，Host里面的Kernel按照线程网格(Grid)形式在显卡硬件(GPU)上执行，其中每个线程网格可包含多个线程块(Block)，每个线程块又可以包含多个线程(Thread)。

在网格中的线程块和线程可以是一维或多维的。

运行在GPU上的CUDA并行计算核函数Kernel，被编译成其设备上的指令集，通过一批线程并行执行。

图1是一个CUDA线程层次结构图。

图1 CUDA线程层次结构
在CUDA中，当指定网格和线程块数的维数后，启动CUDA Kernel，被并行执行的线程指令的次数为其中的线程块的维数与网格维数的乘积，与此同时，CUDA 网格中每个线程都具有一个索引号，即用(thread Idx，block Idx)来表示，以便用户通过该索引号来完成GPU的寻址操作。

2 FDK算法
2.1 FDK加速实现的必要性
FDK算法是一种近似三维重建算法，由于它计算速度快，数学形式简单，所以目前它仍然是公认的最有效的三维图像重建算法。

相比其他三维重建算法，FDK算法计算速度快，然而其时间复杂度高达O(N4)。

N是重建三维图像的尺寸，若重建分辨率为512×512×512的三维图像，重建算法的内部循环次数将近700亿次，仅利用PC机上的CPU来计算是非常耗时的。

人们在这方面做了很多工作，进行优化和改进算法，并从软、硬件方面提出了多种加速策略。

2.2 FDK算法并行性分析
在CT图像重建中，FDK算法一直是应用的主流。

FDK算法实际上是在扇束FBP 算法上做修改得到的。

FDK算法重建具体过程为:
1)对投影数据加权
物体上待重建点位f(x，y，z)的射线源在探测器平面的坐标为(Y，Z)，pβ(Y，Z)表示物体点的投影数据，R为射线源S绕中心轴Z的旋转半径。

修正的权系数的几
何意义是任一射线与中心射线之间夹角的余弦。

2)对二维投影数据，进行滤波，其滤波函数为3)沿射线方向进行反投影重建，获
得密度图像f(x，y，z)。

式中为加权因子，β为射线源与x轴正方向的夹角(即旋转角度)，dβ为在中心平
面上的扇束的投影角增量。

根据上述3个步骤分析，可以发现各个旋转角度下的计算互不影响，相互独立;另
外在反投影步骤中，物体构成的体素间的独立性，它们的计算也相互独立，由此可以采用并行的思想，结合CUDA技术的GPU硬件架构完成算法的加速实现。

3 FDK算法CUDA技术的加速方法
3.1 基于CUDA技术的加速方法
GPU拥有多个流处理器，GPU执行指令的方式是单指令多数据的方式。

该文的加速方法要求首先在CPU端准备好投影数据，然后输入到GPU设备存储器中，待GPU计算完毕以后，再传回到CPU，最后将重建的三维数据保存。

CUDA的编程和优化过程中，需要决定算法中任务的串行和并行的执行部分，由于在FDK算法
3个步骤中实现并行加速策略是不相同的，所以对于这3个步骤按照串行方式进行，它们分别通过3个不同CUDA核进行计算。

基于CUDA的FDK三维图像重建流
程如图2所示。

图2 基于CUDA的FDK算法实现框图
在算法实现中，每次数据的传送都需一定时间开销，所以为了减少时间的消耗，把
许多小的数据传送合并为一次大的数据传送，如果在显存容量允许下，应尽可能将每次多个角度的投影数据从CPU输入GPU。

在这过程中，有效建立内存的访问模式显得尤为重要，尽量使用多种存储器相结合的方式，来有效组织内存访问。

3.2 加权过程的GPU加速
根据实现原理分析，在每个投影角度下的各个像素点的加权计算是相互不影响的，它们仅与计算的像素点坐标有关。

所以为了提高并行性，算法中采用按照像素点划分线程块。

若每个角度下的投影数据的分辨率为N×N，而GPU拥有n个多处理器，把K个角度下的投影数据输入到GPU全局存储器中，在软件架构中，设置每个GPU多处理器中包含a个线程块，每个线程块中又包含b个线程，这样就可以计算出每个线程可以完成K·N·N/(a·b)像素点的加权处理。

该文在算法加权步骤实
现的实验中，采用的数据分辨率为512×512，GTX260多处理器有24个，共设
有30个投影角度，每个多处理器包含512个线程块，每个线程块中拥有512个
线程。

3.3 滤波过程GPU加速
滤波过程可通过投影数据与滤波器(算法中采用的S－L滤波器)的卷积完成，而它
可以借助傅里叶变换来实现。

为了便于计算，该过程首先将投影和滤波器变换到频域，即使两者成为相乘的关系，在这过程中可调用CUDA的并行FFT核函数完成
频域的变换，滤波处理完毕后，用CUDA并行的FFT逆变换，再将其变换到时域。

该算法在完成滤波加速的实现中，实际上是采用了GPU完成FFT并行计算，通过加速FFT，最终完成滤波的加速。

由于在FFT实现中，采用蝶形运行方式，每级蝶形具有一定的相互独立性，所以
可以采用并行实现方法。

并行思想的提出，减少了FFT中实现的循环层数，从而
大大提高了其实现的速度。

FFT的并行实现的方法流程如图3所示。

图3 FFT的并行执行框图
3.4 反投影过程GPU加速
FDK算法滤波步骤完成之后，就进入到反投影的计算过程中。

实践表明，该步骤在FDK算法的实现中是最为耗时的。

因此对该步骤的合理设计，能够很好完成FDK的加速实现效果。

文中在反投影计算过程中，将物体体数据进行三维划分，形成线程块，线程块再次细分若干个线程，每个线程执行一个体素的反投影计算，每次可以并行完成计算K 幅投影数据。

这样做可以有效地降低算法的复杂度，而且提高了算法的并行性。

图4 体数据的三维划分方法
具体步骤描述如下:
1)设置核内块和线程数量，每行为一个线程块，块内拥有多个线程，计算一个或多个蝶形运算;
2)变换前要对数据进行重排，然后数据输入GPU全局存储器，计算过程为原位计算，不需另外开辟存储器空间;
3)循环进行M级蝶形运算，即M次调用CUDA的Kernel函数;
4)循环期间不必每次都要进行Host和Device端的数据传输，只需在循环开始前将数据从Host输入Device，循环结束后，由Device输出到Host。

在FFT并行计算处理中，为了有效地访问内存，考虑到共享存储器位于芯片上，它的存储速度要比全局存储器快很多，所以将输入的投影数从全局存储器载入到共享存储器中，然后在共享存储器内处理数据，最后将处理后的结果存储到全局存储器。

在算法实现中对于大分辨率的三维数据，由于GPU开辟的线程数量有限，所以无法通过一次反投影计算完成体数据的所有体素计算，但可以借助循环设置，通过反复调用Kernel函数来完成。

为了提高存储器的有效优化，在执行循环过程之前，可以将投影数据全部输入GPU，并绑定在其纹理存储器中，这样可以避免反复输
入循环所需数据。

当循环结束后，需要将体数据传输给CPU，进行累加运算。

整
个过程的实现可以提高CPU和GPU的有效合作。

根据上面的分析，具体FDK实现的伪代码如下:输入:Projection P1，P2，…，Pk，Projection size M，N(M=N)
Volume size X，Y，Z输出:Volume f CPU主函数
{Get projection data for i=0 to K/Ido
Transfer projection P1～PI GPU memory P ←Weighting_Kernel //加权P
←Filtering_Kernel //滤波
Bind Filtered projection P as a texture for j=0 to Z/z do
F ←Backprojection_Kernel //反投影end for
end for GPU memory→CPU memory}其中CUDA函数
{Weighting_Kernel(P，M，N)
4 实验结果分析
实验环境:Inter Core 2 Duo(2.53 GHz)CPU，内存为2 GB，Nvidia GTX260显
存为896MB;OS系统是Windows XP SP3;编程环境 Visual Studio2010，安装并配置CUDA工具包。

实验采用360幅512×512的32位浮点型仿真投影数据，三维重建图像大小为512×512×512，也是32位浮点型。

分别使用CPU与GPU完成CT三维重建过程，对比它们的运行时间，结果如表1所示。

表1 CPU和GPU重建时间对比FDK算法过程重建时间/s CPU－FDK GPU－FDK 加速比加权10.5 0.78 13滤波 260.3 3.54 73反投影 4093.6 23.8 172总过程4364.4 28.12 155
从表1可以看出，算法的加权部分提高了13倍，滤波部分提高了73倍，而在实
现最耗时的反投影中提高了172倍，根据计算得到整个算法的加速提高了155倍，而且三维图像完成重建用时28 s多。

从实验结果来看，算法并行得到了很好的实现。

5 总结
该文研究了一种利用GPU通用计算架构(CUDA)的FDK三维图像重建加速算法。

在实现FDK算法的加权、滤波和反投影3个阶段中都采用了并行计算思想。

实验
结果表明，并行图像重建加速方法与单独使用CPU计算方法相比获得了较高的加
速效果，而且还对它们进行了误差分析，最大误差10－3，平均误差保持在 10－4，可以看出两者重建的图像质量相差甚小。

所以文中提出的三维图像加速重建方法获得了较好的结果。

另外，对于CUDA的推出，使得GPU具有更好的可编程性，适合开发人员快速掌握编程方法，缩短程序开发周期。

相信随着CUDA架构逐步成熟和GPU性能的提高，利用GPU的锥束CT三维图像重建将会更快，更能满足实时准确的图像重建
要求。

参考文献:
［1］ Fang Xu，Klaus Mueller.Ultra － fast 3D filtered backprojection on commodity graphics hardware［J］.IEEE International Symposium on Biomedical Imaging:Macro to Nanno，Arlington United states，
2004(1):571 －574.
［2］陈雪松.CT图像重建可扩展多DSP并行计算系统结构［J］.清华大学学报:自然科学版，2004，44(3):44 －47.
［3］邹永宁，谭辉，黄亮.CT图像重建加速的几种方法［J］.计算机系统应用，2009，18(4):167 －170.
［4］ B Cabral，N Cam，JForan.Accelerated volume rendering and
tomographic reconstruction using texturemapping hardware
［J］.Proc.Symp.On Volume Visualization，1994(10):91－98.
［5］ Klaus Mueller，Fang Xu.Practical considerations for GPU －accelerated CT［J］.IEEE International Symposium on Biomedical Imaging:Macro to Nano，2006(7):1184－1187.
［6］戴智晟，陈志强.用通用显卡加速三维锥束T－FDK重建算法［J］.清华大学学报:自然科学版，2006，46(9):1589－1592.
［7］马东平，曾理.GPU多重纹理加速三维CT图像重建［J］.计算机工程与应用，2008，44(7):227 －230.
［8］ Holger Scherl，Benjamin Keck，Markus Kowarschik，et al.Fast GPU
－Based CTReconstruction using the Common Unified Device Architecture(CUDA)［J］.Nuclear Science Symposium Conference Record，2007，26(280):4464－4466.
［9］王珏，曹思远，邹永宁.利用CUDA技术实现锥束CT图像快速重建［J］.核电子学与探测技术，2010，30(3):315－320.
［10］孙毅刚，孙修宇，张红颖.基于现代GPU的实时锥束重建算法研究［J］.核电子学与探测技术，2010，30(9):1260－1265.
［11］肖江，胡柯良，邓元勇.基于CUDA的矩阵乘法和FFT性能测试［J］.计算机工程，2009，35(10):7 －10.
［12］ R De Beer，D Van Ormondt，FCesare Di，etal.Accelerating batched 1D－FFT with a CUDA－capable computer［C］.Imaging Systems and Techniques(IST)，IEEE International Conference，2010:446 －451.
［13］ Kauker D，Sanftmann H，Frey S.Memory Saving Discrete Fourier Transform on GPUs［J］.Computer and Information Technology(CIT)，
2010(6):1152 －1157.。