摘  要： 提出了一種基于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.