CUDA加速分形火焰绘制

摘 要： 提出了一种基于CUDA的并行分形火焰绘制算法，该算法利用了GPU的单指令多线程的特点，将常用于迭代函数系统（IFS）的传统的混沌游戏（Chaos Game）算法作了并行化修改，并在图形输出时利用了CUDA与OpenGL互操作加速分形火焰绘制。实验证明，该并行方法比CPU上运行的普通算法快了15倍左右，能够实时绘制分形火焰图形。在上述基本算法的基础上，又进一步研究了消除分支分歧的改进算法，改进算法的运行时间具有相对于变换函数数量的恒定性，多数情况下比基本算法性能更优越。
关键词： 分形火焰；迭代函数系统；统一计算设备架构；图形处理器

　分形通常被定义为“一个粗糙或零碎的几何形状，可以分成数个部分，且每一部分都（至少近似地）是整体缩小后的形状”，即具有自相似的性质。分形一词于1975年由曼德博创造出，来自拉丁文“frāctus”，有“零碎”、“破裂”之意。分形有几种类型，可以分别依据表现出的精确自相似性、半自相似性和统计自相似性来定义。虽然分形是一个数学构造，但也可以在自然界中找到。分形在艺术作品、医学、土力学、地震学和技术分析中都有应用。
　分形火焰是迭代函数系统IFS（Iterated Function System）分形的一种更复杂的变种形式，与普通的IFS相比，分形火焰能产生变化多端、引人入胜的图形，但与此同时，其计算复杂度很高。
　图形处理器（GPU）原本是处理计算机图形的专用设备，近十年来，由于高清晰度复杂图形实时处理的需求，GPU发展成为高并行度、多线程、多核的处理器。目前，主流GPU的运算能力已超过主流通用CPU，从发展趋势上来看，将来差距会越拉越大。GPU卓越的性能对开发GPGPU（使用GPU进行通用计算）非常具有吸引力。统一计算设备架构CUDA（Compute Unified Device Architecture）是NVIDIA公司伴随着统一渲染架构而推出的一种通用的GPU编程模型，可以将GPU视为一个并行数据计算的设备，对所进行的计算进行分配和管理[1]。在CUDA的架构中，通用计算不再像过去的GPGPU那样必须将计算映射到图形API中，开发者无需学习复杂的显示芯片的指令或是特殊的结构，CUDA编程语言只是对标准的C语言作了少量扩展，因此开发门槛大大降低了。目前有很多基于CUDA的GPU计算的研究成果，有不少计算问题特别适合这种计算模式。关于CUDA的相关应用可参阅参考文献[2]和参考文献[3]。
　CUDA加速的分形绘制的研究并不多，其中NVIDIA的SDK提供了基于CUDA的加速Mandelbrot、Julia集生成的实例[4]。其算法的基本过程是平面图的每一个像素对应CUDA并行计算中的一个线程，该并行方法能比CPU上的普通方法加速几十倍。参考文献[5]描述了基于CUDA的并行分形IFS点递归绘制算法。
本文的研究是利用GPU的计算能力加速分形火焰绘制，使之能更加实用。
1 背景及相关知识
1.1 经典的迭代函数系统
　数学中，迭代函数系统是构造分形的一种方法，由此产生的结构是自相似的。IFS分形可以是任何维度[6]，但通常在二维空间中计算和绘制。分形由一些自身的拷贝联合构成，每个拷贝根据一个函数来作变换。函数通常是收缩的，意思是它们会使图像点更靠近，使形状更小。因此，IFS分形的形状由一些可能重叠的自身拷贝组成，其中每个拷贝也由自身的一些拷贝组成，如此无限下去。这也是自相似分形特性的来源。

1.2 分形火焰简介
　分形火焰是1992年由DRAVES S提出的[7]，它可以说是分形迭代函数的一种形式，但有以下3方面不同：
　（1）迭代时使用非线性函数而不是仿射变换；
　（2）色调映射显示是对数密度而不是线性的；
　（3）颜色是根据结构（即通过采取的递归路径）决定的，而不是采用单色的或根据点的密度决定。
色调映射和着色旨在显示尽可能多的分形的详细信息，这样通常会产生更美观、更绚丽的图像。
分形火焰的每个函数具有如下的形式：

3 分形火焰并行算法详细描述
　在分形火焰直方图生成阶段，使用串行混沌游戏算法的过程是：选择一个点开始，并根据概率挑选一个函数将该点代入来计算，得到下一个点，然后用新点继续进行下一次迭代，如此不断迭代下去。
　本文提出的并行算法是传统混沌游戏算法的扩展：算法开始时选择n个而不是一个起始点，参考文献[9]说明了选择n个起始点的方法开始迭代和选择1个起始点开始迭代得到的结果一致。并行的实现是通过将一个线程分配给某一个起始点，n个线程同时分别独立执行混沌游戏。通过这种方式，并行算法提高了速度。
　当然，上述只是并行算法的主要思想，因为分形火焰比较复杂，具体的算法还有很多细节问题需要考虑。实现主要包括3个在GPU上并行执行的内核函数：warm_up、iterate_batch和output_for_rending。与常见的串行混沌游戏一样，在warm_up函数中的每个线程选择一个起始点，迭代十几次，但只保持最后的结果，放弃前面迭代得到的值。warm_up函数计较简单，没必要描述实现细节。Iterate_batch函数接收warm_up函数生成的点并迭代数十或数百遍，然后计算得到的每个点的直方图。Iterate_batch函数的核心思想就是n个点同时迭代的并行混沌游戏算法，也是生成分形火焰图形的核心。其伪代码如下：
//输入：ractalInfo存储该fractal flame信息的结构体；iterPosStateBuffer是warm_up后存储点位置的数组；iterColorStateBuffer是warm_up后存储颜色的数组；randBuffer用户指定的每个函数选择的概率
//输出：accumBuffer存储累计直方图信息，该信息在绘制阶段会用到
Iterate_batch（）
{lid=threadIdx.x；
gid=（blockIdx.x*blockDim.x+threadIdx.x）；
pos=iterPosStateBuffer[gid]；
color=iterColorStateBuffer[gid]；
for（iter=0；iter<ITERCOUNT；iter++）
{fIndex=chooseRandomBranch（randBuffer+lid，fractalInfo）；//根据randBuffer选择一个函数
iterate（&pos，&color，fIndex，fractalInfo）；
//迭代生成新的点和颜色
screenPos=getPos（fractalInfo，pos）；
//将计算得到的点的位置转换为屏幕上点的位置
calhistogram（color，screenPos，accumBuffer）；
//计算每个屏幕点的统计直方图和颜色，保存到accumBuffer
}}
Output_for_rending函数主要完成图形绘制功能，它接收从iterate_batch传递来的直方图信息，作色调映射和伽玛校正，并将每个像素的最终显示颜色放在outputBuffer。其伪代码如下：
//输入：ractalInfo存储了该fractal flame信息的结构体；
accumBuffer是从Iterate_batch函数传递过来的，保存了直方图信息
//输出：outputBuffer保存了最终数据，用于通过OpenGL绘制最终图形
output_for_rending（ractalInfo，accumBuffer）
{x=（blockIdx.x*blockDim.x+threadIdx.x）；
y=（blockIdx.y*blockDim.y+threadIdx.y）；
pix=tonemap（fractalInfo，accumBuffer，x，y）；
//根据直方图信息，使用色调映射和伽玛校正产生实际
//显示的颜色
setoutput（outputBuffer，x，y，pix）；
//将每个像素的显示颜色保存到outputBuffer
}
　根据实验观察，图像绘制阶段大约消耗总时间的70%，为了加速绘制过程，需要利用CUDA和Open GL互操作[10]。互操作的基本方法是将Open GL缓冲区映射到CUDA的内存空间。CUDA用于计算和数据生成，Open GL用来绘制像素或顶点，因为它们共享同一内存空间，无需在CPU内存和GPU内存间移动数据，所以速度非常快。调用output_for_rending函数之前，需要做一些工作使CUDA和Open GL相关联，并设置outputBuffer与CUDA和Open GL共享。
4 对基本并行算法的改进
　分形火焰的基本并行实现是每个线程对应一个数据点，各线程随机选择一个变换函数计算，用计算得到的数据点替换原数据点，这个过程重复m次。该并行算法需要每个线程根据随机数选择变换函数，因而会导致CUDA严重的分支分歧[1]。变换函数越多，越有可能绘制出更有趣的图像，因此，应该尽可能允许有更多的变换函数。然而用到的变换函数越多，分支分歧问题就越严重。
　本文提出一种改进的算法，该算法通过预分类能移除随机选择函数，可以消除分支分歧。具体方法是通过随机数据访问来替代随机选定函数，即每个线程分配一个固定的函数，在每次迭代中随机选择一个数据点。这种选择是通过数据和线程索引之间的一个随机双射函数映射实现的。通过这种方法，指令能够以最佳方式静态地分配给线程，预先计算好固定的置换，它们不依赖于动态数据。每个线程使用分配的函数，通过置换索引间接地访问数据数组，然后将该函数的计算结果写回。为了避免写访问时的条件竞争，要么结果数据写回读取的位置，要么需要第二个数组来存储数据点，每次迭代挑选一个新的置换。
　上述的过程需要通过创建多个包含数据点索引及随机数的数组来进行置换，随机数可由Mersenne Twister方法生成[11]，数组根据随机数排序。所有的变换函数是在一个大的switch语句中实现的。用一个结构描述变化索引及缩放因子，该结构存储在常量内存中。
5 实验结果
　上述并行分形火焰算法是在NVIDIA GeForce GTX9800+上实现的，该GPU有128个频率为1.836 GHz的流处理器，分为16个SM，896 MB显存，装配在CPU为Intel Core 2 Duo E7400 2.8 GHz，内存为1 GB的计算机上。为了作对比，同时实现了CPU上运行的串行版本。实验结果表明，本文提出的基于CUDA的基本并行算法在大多数情况下比CPU上的串行算法快了约15倍，能做到分形火焰实时绘制。
　实验表明，消除分支分歧的算法与基本并行算法相比，基本并行算法的运行时间与变换函数的数目呈现线性增长关系，而消除分支分歧的改进算法运行时间恒定，多数情况下性能都比基本算法高。在变换函数数目为20个时，消除了分支分歧的算法比基本算法快了约两倍。
　本文提出了一种并行的分形火焰算法，该方法利用了GPU的计算能力及CUDA和Open GL互操作方式，速度快到足够实时高帧速率绘制。该算法使得在生成分形火焰图形时，能实时更改参数并观察绘制效果，可以有效地帮助加深对迭代函数系统的认识。本文方法对其他图像实时绘制的应用也有很高的参考价值。
参考文献
[1] NVIDIA Corporation. NVIDIA CUDA programming guide version 3.2[EB/OL]. （2011-03）. http：//developer、nvidia、com/cuda、
[2] Hwu Wenmei， RODRIGUES C， RYOO S， et al. Compute unified device architecture application suitability[J]. Computing in Science and Engineering， 2009，11（3）： 16-26.
[3] Hwu Wenmei. GPU computing gems[M]. New York： Morgan Kaufmann，2011.
[4] GRANGER M. Mandelbrot CUDA SDK code samples[EB/OL].http：//developer、nvidia、com/cuda、2011-03.
[5] KWANIEWSKI B. CUDA based parallel version of point recursive rendering algorithm of IFS attractor[C]. 12th International PhD Workshop OWD， 2010： 23-26.
[6] WHITELAW M. Metacreation： art and artificial life[M]. MIT Press， 2004.
[7] DRAVES S， RECKASE E. The fractal flame algorithm[EB/OL].http：//flam3、com/flame draves.pdf. 2008-11.
[8] NIKIEL S. Iterated function systems for real-time image synthesis[M]. London： Springer， 2007.
[9] DUTIL N. Construction of fractal objects with iterated function systems[EB/OL].（2000-10）.http：//www。cs。mcgill。ca/~ndutil/project.pdf .
[10] 刘进锋，郭雷.CUDA和OpenGL互操作的实现及分析[J].微型机与应用，2011（23）：40-43.
[11] MATSUMOTO M， NISHIMURA T. Mersenne twister： a 623-dimensionally equidistributed uniform pseudorandom number generator[J]. ACM Transactions on Modeling and Computer Simulation， 1998，8（1）：3-30.