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

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

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

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