摘  要： 提出了一種基于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ā)展趨勢(shì)上來看，將來差距會(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 分形火焰簡(jiǎn)介
　分形火焰是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ì)較簡(jiǎn)單，沒必要描述實(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]。互操作的基本方法是將Open 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í)的條件競(jìng)爭(zhēng)，要么結(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ǔ)在常量?jī)?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)線性增長(zhǎng)關(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.