摘  要： 提出了一種基于CUDA的并行分形火焰繪制算法,，該算法利用了GPU的單指令多線程的特點(diǎn),，將常用于迭代函數(shù)系統(tǒng)（IFS）的傳統(tǒng)的混沌游戲（Chaos Game）算法作了并行化修改，并在圖形輸出時(shí)利用了CUDA與OpenGL互操作加速分形火焰繪制,。實(shí)驗(yàn)證明,，該并行方法比CPU上運(yùn)行的普通算法快了15倍左右，能夠?qū)崟r(shí)繪制分形火焰圖形,。在上述基本算法的基礎(chǔ)上,，又進(jìn)一步研究了消除分支分歧的改進(jìn)算法，改進(jìn)算法的運(yùn)行時(shí)間具有相對(duì)于變換函數(shù)數(shù)量的恒定性,，多數(shù)情況下比基本算法性能更優(yōu)越,。
關(guān)鍵詞： 分形火焰；迭代函數(shù)系統(tǒng),；統(tǒng)一計(jì)算設(shè)備架構(gòu),；圖形處理器

　分形通常被定義為“一個(gè)粗糙或零碎的幾何形狀，可以分成數(shù)個(gè)部分,，且每一部分都（至少近似地）是整體縮小后的形狀”,，即具有自相似的性質(zhì)。分形一詞于1975年由曼德博創(chuàng)造出，來自拉丁文“frāctus”,，有“零碎”,、“破裂”之意。分形有幾種類型,，可以分別依據(jù)表現(xiàn)出的精確自相似性,、半自相似性和統(tǒng)計(jì)自相似性來定義。雖然分形是一個(gè)數(shù)學(xué)構(gòu)造,，但也可以在自然界中找到。分形在藝術(shù)作品,、醫(yī)學(xué),、土力學(xué)、地震學(xué)和技術(shù)分析中都有應(yīng)用,。
　分形火焰是迭代函數(shù)系統(tǒng)IFS（Iterated Function System）分形的一種更復(fù)雜的變種形式,，與普通的IFS相比，分形火焰能產(chǎn)生變化多端,、引人入勝的圖形,，但與此同時(shí)，其計(jì)算復(fù)雜度很高,。
　圖形處理器（GPU）原本是處理計(jì)算機(jī)圖形的專用設(shè)備,，近十年來，由于高清晰度復(fù)雜圖形實(shí)時(shí)處理的需求,，GPU發(fā)展成為高并行度,、多線程、多核的處理器,。目前,，主流GPU的運(yùn)算能力已超過主流通用CPU，從發(fā)展趨勢上來看,，將來差距會(huì)越拉越大,。GPU卓越的性能對(duì)開發(fā)GPGPU（使用GPU進(jìn)行通用計(jì)算）非常具有吸引力。統(tǒng)一計(jì)算設(shè)備架構(gòu)CUDA（Compute Unified Device Architecture）是NVIDIA公司伴隨著統(tǒng)一渲染架構(gòu)而推出的一種通用的GPU編程模型,，可以將GPU視為一個(gè)并行數(shù)據(jù)計(jì)算的設(shè)備,，對(duì)所進(jìn)行的計(jì)算進(jìn)行分配和管理[1]。在CUDA的架構(gòu)中,，通用計(jì)算不再像過去的GPGPU那樣必須將計(jì)算映射到圖形API中,，開發(fā)者無需學(xué)習(xí)復(fù)雜的顯示芯片的指令或是特殊的結(jié)構(gòu)，CUDA編程語言只是對(duì)標(biāo)準(zhǔn)的C語言作了少量擴(kuò)展,，因此開發(fā)門檻大大降低了,。目前有很多基于CUDA的GPU計(jì)算的研究成果，有不少計(jì)算問題特別適合這種計(jì)算模式。關(guān)于CUDA的相關(guān)應(yīng)用可參閱參考文獻(xiàn)[2]和參考文獻(xiàn)[3],。
　CUDA加速的分形繪制的研究并不多,，其中NVIDIA的SDK提供了基于CUDA的加速M(fèi)andelbrot、Julia集生成的實(shí)例[4],。其算法的基本過程是平面圖的每一個(gè)像素對(duì)應(yīng)CUDA并行計(jì)算中的一個(gè)線程,，該并行方法能比CPU上的普通方法加速幾十倍。參考文獻(xiàn)[5]描述了基于CUDA的并行分形IFS點(diǎn)遞歸繪制算法,。
本文的研究是利用GPU的計(jì)算能力加速分形火焰繪制,，使之能更加實(shí)用。
1 背景及相關(guān)知識(shí)
1.1 經(jīng)典的迭代函數(shù)系統(tǒng)
　數(shù)學(xué)中,，迭代函數(shù)系統(tǒng)是構(gòu)造分形的一種方法,，由此產(chǎn)生的結(jié)構(gòu)是自相似的。IFS分形可以是任何維度[6],，但通常在二維空間中計(jì)算和繪制,。分形由一些自身的拷貝聯(lián)合構(gòu)成，每個(gè)拷貝根據(jù)一個(gè)函數(shù)來作變換,。函數(shù)通常是收縮的,，意思是它們會(huì)使圖像點(diǎn)更靠近，使形狀更小,。因此,，IFS分形的形狀由一些可能重疊的自身拷貝組成，其中每個(gè)拷貝也由自身的一些拷貝組成,，如此無限下去,。這也是自相似分形特性的來源。

1.2 分形火焰簡介
　分形火焰是1992年由DRAVES S提出的[7],，它可以說是分形迭代函數(shù)的一種形式,，但有以下3方面不同：
　（1）迭代時(shí)使用非線性函數(shù)而不是仿射變換,；
?。?）色調(diào)映射顯示是對(duì)數(shù)密度而不是線性的；
?。?）顏色是根據(jù)結(jié)構(gòu)（即通過采取的遞歸路徑）決定的,，而不是采用單色的或根據(jù)點(diǎn)的密度決定。
色調(diào)映射和著色旨在顯示盡可能多的分形的詳細(xì)信息,，這樣通常會(huì)產(chǎn)生更美觀,、更絢麗的圖像。
分形火焰的每個(gè)函數(shù)具有如下的形式：

3 分形火焰并行算法詳細(xì)描述
　在分形火焰直方圖生成階段,，使用串行混沌游戲算法的過程是：選擇一個(gè)點(diǎn)開始,，并根據(jù)概率挑選一個(gè)函數(shù)將該點(diǎn)代入來計(jì)算,，得到下一個(gè)點(diǎn)，然后用新點(diǎn)繼續(xù)進(jìn)行下一次迭代,，如此不斷迭代下去,。
　本文提出的并行算法是傳統(tǒng)混沌游戲算法的擴(kuò)展：算法開始時(shí)選擇n個(gè)而不是一個(gè)起始點(diǎn)，參考文獻(xiàn)[9]說明了選擇n個(gè)起始點(diǎn)的方法開始迭代和選擇1個(gè)起始點(diǎn)開始迭代得到的結(jié)果一致,。并行的實(shí)現(xiàn)是通過將一個(gè)線程分配給某一個(gè)起始點(diǎn),，n個(gè)線程同時(shí)分別獨(dú)立執(zhí)行混沌游戲。通過這種方式,，并行算法提高了速度,。
　當(dāng)然，上述只是并行算法的主要思想,，因?yàn)榉中位鹧姹容^復(fù)雜,，具體的算法還有很多細(xì)節(jié)問題需要考慮。實(shí)現(xiàn)主要包括3個(gè)在GPU上并行執(zhí)行的內(nèi)核函數(shù)：warm_up,、iterate_batch和output_for_rending。與常見的串行混沌游戲一樣,，在warm_up函數(shù)中的每個(gè)線程選擇一個(gè)起始點(diǎn),，迭代十幾次，但只保持最后的結(jié)果,，放棄前面迭代得到的值,。warm_up函數(shù)計(jì)較簡單，沒必要描述實(shí)現(xiàn)細(xì)節(jié),。Iterate_batch函數(shù)接收warm_up函數(shù)生成的點(diǎn)并迭代數(shù)十或數(shù)百遍,，然后計(jì)算得到的每個(gè)點(diǎn)的直方圖。Iterate_batch函數(shù)的核心思想就是n個(gè)點(diǎn)同時(shí)迭代的并行混沌游戲算法,，也是生成分形火焰圖形的核心,。其偽代碼如下：
//輸入：ractalInfo存儲(chǔ)該fractal flame信息的結(jié)構(gòu)體；iterPosStateBuffer是warm_up后存儲(chǔ)點(diǎn)位置的數(shù)組,；iterColorStateBuffer是warm_up后存儲(chǔ)顏色的數(shù)組,；randBuffer用戶指定的每個(gè)函數(shù)選擇的概率
//輸出：accumBuffer存儲(chǔ)累計(jì)直方圖信息，該信息在繪制階段會(huì)用到
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）,；//根據(jù)randBuffer選擇一個(gè)函數(shù)
iterate（&pos,，&color，fIndex,，fractalInfo）,；
//迭代生成新的點(diǎn)和顏色
screenPos=getPos（fractalInfo，pos）,；
//將計(jì)算得到的點(diǎn)的位置轉(zhuǎn)換為屏幕上點(diǎn)的位置
calhistogram（color,，screenPos，accumBuffer）,；
//計(jì)算每個(gè)屏幕點(diǎn)的統(tǒng)計(jì)直方圖和顏色,，保存到accumBuffer
}}
Output_for_rending函數(shù)主要完成圖形繪制功能，它接收從iterate_batch傳遞來的直方圖信息,，作色調(diào)映射和伽瑪校正,，并將每個(gè)像素的最終顯示顏色放在outputBuffer。其偽代碼如下：
//輸入：ractalInfo存儲(chǔ)了該fractal flame信息的結(jié)構(gòu)體,；
accumBuffer是從Iterate_batch函數(shù)傳遞過來的,，保存了直方圖信息
//輸出：outputBuffer保存了最終數(shù)據(jù)，用于通過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）,；
//根據(jù)直方圖信息,，使用色調(diào)映射和伽瑪校正產(chǎn)生實(shí)際
//顯示的顏色
setoutput（outputBuffer，x,，y,，pix）；
//將每個(gè)像素的顯示顏色保存到outputBuffer
}
　根據(jù)實(shí)驗(yàn)觀察,，圖像繪制階段大約消耗總時(shí)間的70%,，為了加速繪制過程，需要利用CUDA和Open GL互操作[10],?；ゲ僮鞯幕痉椒ㄊ菍pen GL緩沖區(qū)映射到CUDA的內(nèi)存空間。CUDA用于計(jì)算和數(shù)據(jù)生成,，Open GL用來繪制像素或頂點(diǎn),，因?yàn)樗鼈児蚕硗粌?nèi)存空間，無需在CPU內(nèi)存和GPU內(nèi)存間移動(dòng)數(shù)據(jù),，所以速度非?？?。調(diào)用output_for_rending函數(shù)之前，需要做一些工作使CUDA和Open GL相關(guān)聯(lián),，并設(shè)置outputBuffer與CUDA和Open GL共享,。
4 對(duì)基本并行算法的改進(jìn)
　分形火焰的基本并行實(shí)現(xiàn)是每個(gè)線程對(duì)應(yīng)一個(gè)數(shù)據(jù)點(diǎn)，各線程隨機(jī)選擇一個(gè)變換函數(shù)計(jì)算,，用計(jì)算得到的數(shù)據(jù)點(diǎn)替換原數(shù)據(jù)點(diǎn),，這個(gè)過程重復(fù)m次。該并行算法需要每個(gè)線程根據(jù)隨機(jī)數(shù)選擇變換函數(shù),，因而會(huì)導(dǎo)致CUDA嚴(yán)重的分支分歧[1],。變換函數(shù)越多，越有可能繪制出更有趣的圖像,，因此,，應(yīng)該盡可能允許有更多的變換函數(shù)。然而用到的變換函數(shù)越多,，分支分歧問題就越嚴(yán)重,。
　本文提出一種改進(jìn)的算法，該算法通過預(yù)分類能移除隨機(jī)選擇函數(shù),，可以消除分支分歧,。具體方法是通過隨機(jī)數(shù)據(jù)訪問來替代隨機(jī)選定函數(shù)，即每個(gè)線程分配一個(gè)固定的函數(shù),，在每次迭代中隨機(jī)選擇一個(gè)數(shù)據(jù)點(diǎn)。這種選擇是通過數(shù)據(jù)和線程索引之間的一個(gè)隨機(jī)雙射函數(shù)映射實(shí)現(xiàn)的,。通過這種方法,，指令能夠以最佳方式靜態(tài)地分配給線程，預(yù)先計(jì)算好固定的置換,，它們不依賴于動(dòng)態(tài)數(shù)據(jù),。每個(gè)線程使用分配的函數(shù)，通過置換索引間接地訪問數(shù)據(jù)數(shù)組,，然后將該函數(shù)的計(jì)算結(jié)果寫回,。為了避免寫訪問時(shí)的條件競爭，要么結(jié)果數(shù)據(jù)寫回讀取的位置,，要么需要第二個(gè)數(shù)組來存儲(chǔ)數(shù)據(jù)點(diǎn),，每次迭代挑選一個(gè)新的置換。
　上述的過程需要通過創(chuàng)建多個(gè)包含數(shù)據(jù)點(diǎn)索引及隨機(jī)數(shù)的數(shù)組來進(jìn)行置換,，隨機(jī)數(shù)可由Mersenne Twister方法生成[11],，數(shù)組根據(jù)隨機(jī)數(shù)排序。所有的變換函數(shù)是在一個(gè)大的switch語句中實(shí)現(xiàn)的,。用一個(gè)結(jié)構(gòu)描述變化索引及縮放因子,，該結(jié)構(gòu)存儲(chǔ)在常量內(nèi)存中,。
5 實(shí)驗(yàn)結(jié)果
　上述并行分形火焰算法是在NVIDIA GeForce GTX9800+上實(shí)現(xiàn)的，該GPU有128個(gè)頻率為1.836 GHz的流處理器,，分為16個(gè)SM,，896 MB顯存，裝配在CPU為Intel Core 2 Duo E7400 2.8 GHz,，內(nèi)存為1 GB的計(jì)算機(jī)上,。為了作對(duì)比，同時(shí)實(shí)現(xiàn)了CPU上運(yùn)行的串行版本,。實(shí)驗(yàn)結(jié)果表明,，本文提出的基于CUDA的基本并行算法在大多數(shù)情況下比CPU上的串行算法快了約15倍，能做到分形火焰實(shí)時(shí)繪制,。
　實(shí)驗(yàn)表明,，消除分支分歧的算法與基本并行算法相比，基本并行算法的運(yùn)行時(shí)間與變換函數(shù)的數(shù)目呈現(xiàn)線性增長關(guān)系,，而消除分支分歧的改進(jìn)算法運(yùn)行時(shí)間恒定,，多數(shù)情況下性能都比基本算法高。在變換函數(shù)數(shù)目為20個(gè)時(shí),，消除了分支分歧的算法比基本算法快了約兩倍,。
　本文提出了一種并行的分形火焰算法，該方法利用了GPU的計(jì)算能力及CUDA和Open GL互操作方式,，速度快到足夠?qū)崟r(shí)高幀速率繪制,。該算法使得在生成分形火焰圖形時(shí)，能實(shí)時(shí)更改參數(shù)并觀察繪制效果,，可以有效地幫助加深對(duì)迭代函數(shù)系統(tǒng)的認(rèn)識(shí),。本文方法對(duì)其他圖像實(shí)時(shí)繪制的應(yīng)用也有很高的參考價(jià)值。
參考文獻(xiàn)
[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] 劉進(jìn)鋒,，郭雷.CUDA和OpenGL互操作的實(shí)現(xiàn)及分析[J].微型機(jī)與應(yīng)用，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.