分形/計算機圖形技術/2D/最佳化

 "Your current attempt runs in about O(N^2). So even if you had enough memory, it won't finish in your lifetime." Alexander Yee

   "premature optimization is the root of all evil" - Donald Knuth in paper "Structured Programming With Go To Statements"

           "good algorithms beat optimized code" Bruce Dawson ( Fractal extreme )

軟體最佳化

• 程式碼剖析
• 基準測試

如何詢問最佳化問題 ?

• "如果你已經實際剖析了程式碼，有具體的程式碼片段，以便每個人都可以執行相同的程式碼來檢視其效能，並且你在某個地方釋出了這個庫，可以檢視（GitHub、Bitbucket 等），那麼我不一定認為在程式碼審查中進行討論有問題。包括你選擇的分析器中的結果，並識別瓶頸，將有助於保持主題的集中性。
• 如果你剛開始編寫程式碼，但已經剖析了一小段程式碼，並展示了異常的效能，那麼在 Stack Overflow 上提問是可接受的。包括你選擇的分析器中的結果，並識別瓶頸，將有助於保持主題的集中性。
• 如果你只是讓別人幫你最佳化程式碼，那就不要在任何地方提問。得到一個程式碼轉儲並被要求找到效能瓶頸，這只是我在專業工作中才會做的事，而不是在我的空閒時間。" (Makoto)

另請參閱 : 程式碼審查

內部迴圈

• 維基百科中的內部迴圈

數值計算的最佳化

數字型別

• double 比 float 快
• "在我的測試（使用擾動）中，在 Fraktaler-3 中，x87 長雙精度浮點數在 AMD 2700X CPU 上比雙精度浮點數慢大約 4.2 倍。AMD RX 580 GPU 上的 floatexp（float + int32）比 AMD 2700X CPU 上的 x87 長雙精度浮點數快大約 2.3 倍，因此根據硬體的不同，x87 長雙精度浮點數可能沒有必要..."fractalforums.org: from-java-to-cplusplus

距離

• 距離近似
• 科比 1/4 定理

分解

• 四叉樹分解

"I did a web-based quad-tree Mandelbrot once.  I labelled the quadrants (eg a,b,c,d) the tile files were in directories (something like a.png b.png a/a.png a/b.png a/c/a/b.png), along with an sqlite database that said which tiles  were boring (all same colour), along with their period (0 for 100% exterior, 1 for 100% inside the main cardioid, 2 for 100% inside the main circle, etc).  The client was really simple, just a web page with tile paths; the server  would do a redirect to the relevant boring tile 0.png 1.png 2.png etc if any parent of the tile was boring, otherwise it'd try to serve the tile.png.  Meanwhile there was another process that just generated the tiles and populated the database.  I can't remember if I made 404 tile not found errors bump the priority of the tile for the renderer, but maybe.

Eventually i abandoned the project because recalculating from scratch is typically faster and more flexible for shallow zooms, and for deep zooms the storage is prohibitive." Claude Heiland-Allen 

對稱性

映象對稱

引數平面和曼德布洛特集合 在平面中關於 x 軸對稱 :

colour(x,y) = colour(x,-y)

這裡有一個 C 程式碼示例，展示瞭如何將對稱影像（其軸位於其高度的中間）分成 3 部分 

• 對稱軸
• 2 個對稱部分 : 軸上和軸下

// fill array using mirror symmetry of image 
// uses global var :  ...
int FillArray(unsigned char data[] )
{
   
 unsigned char Color; // gray from 0 to 255

// draw axis of symmetry
iy = iyAxisOfSymmetry; 
for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy));

// draw symmetric parts : above and below axis 
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
  for(ix=ixMin; ix<=ixMax; ++ix)

  { // above axis compute color and save it to the array
    iy = iyAxisOfSymmetry + iyAbove;
    Color = GiveColor(ix, iy);
    PlotPoint(ix, iy, Color ); 
    // below the axis only copy Color the same as above without computing it 
    PlotPoint(ix, iyAxisOfSymmetry - iyAbove , Color ); 
   }   
 return 0;
}

另請參閱用於一般情況的 Pascal 程式碼 (Delphi)

旋轉對稱

動力學平面和朱利亞集合 具有旋轉對稱性。它可以用來加速計算 

// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char data[] )
{
   
 unsigned char Color; // gray from 0 to 255

printf("axis of symmetry \n"); 
iy = iyAxisOfSymmetry;

for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy));

// above and below axis 
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
  {printf(" %d from %d\n", iyAbove, iyAboveMax); //info 
  for(ix=ixMin; ix<=ixMax; ++ix)

  { // above axis compute color and save it to the array
    iy = iyAxisOfSymmetry + iyAbove;
    Color = GiveColor(ix, iy);
    PlotPoint(ix, iy, Color ); 
    // below the axis only copy Color the same as above without computing it 
    PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color ); 
   } 
}  
 return 0;
}

複數

不顯式定義複數或不使用複數型別的計算通常更快。

平行計算

在曼德布洛特和朱利亞集合中，每個點/畫素都是獨立計算的，因此可以輕鬆地將其分解成更小的任務並同時執行。

平行計算型別 : 維基百科描述

• 多核（計算）
  • Parallella 和 OpenCl
  • 英特爾 MIC
  • MPI
  • OpenMP
• 通用圖形處理器計算 (GPGPU)
  • OpenCl
  • CUDA
  • WebGl

向量化

C

• gcc 向量擴充套件
• 英特爾® 高階向量擴充套件

c++

• Eigen 中的 Mandelbrot - 顯式向量化

R

這張圖片是使用向量化程式碼製作的。它的生成時間非常短。

向量化

以下是上面圖片的 R 程式碼，包含註釋：

## based on the R code by Jarek Tuszynski
## http://commons.wikimedia.org/wiki/File:Mandelbrot_Creation_Animation_%28800x600%29.gif

iMax=20; # number of frames and iterations

## world coordinate ( float )
cxmin = -2.2;
cxmax = 1.0;
cymax = 1.2;
cymin = -1.2;

## screen coordinate ( integer)

dx=800; dy=600 ;           # define grid size

## sequences of  values = rows and columns
cxseq = seq(cxmin, cxmax, length.out=dx);
cyseq = seq(cymin, cymax, length.out=dy);
## sequences of  values = rows * columns
cxrseq = rep(cxseq, each = dy);
cyrseq = rep(cyseq,  dx);
C = complex( real=cxrseq, imag = cyrseq); # complex vector
C = matrix(C,dy,dx)       # convert from vector to matrix
# Now C is a matrix of c points coordinates (cx,cy)

# allocate memory for all the frames
F = array(0, dim = c(dy,dx,iMax)) # dy*dx*iMax array filled with zeros

Z = 0                     # initialize Z to zero
for (i in 1:iMax)         # perform iMax iterations
{
  Z = Z^2+C               # iteration; uses n point to compute n+1 point
  
  F[,,i]  = exp(-abs(Z))              # an omitted index is used to represent an entire column or row
     # store magnitude of the complex number, scale it using an exponential function to emphasize fine detail, 
}

# color and save F array to the file 
library(caTools)          # load library with write.gif function
# mapped to color palette jetColors . Dark red is a very low number, dark blue is a very high number
jetColors = colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
write.gif(F, "Mandelbrot.gif", col=jetColors, delay=100, transparent=0)

多核

使用 OpenMP 和對稱性的圖片和原始碼

"我總是建立與我擁有的核心數量相同的執行緒。超過這個數量，你的系統會花費太多時間進行任務切換" - Duncan C

使用 OpenMP

// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char data[] )
{
   
 unsigned char Color; // gray from 0 to 255

printf("axis of symmetry \n"); 
iy = iyAxisOfSymmetry; 
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy));

/*
The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
The use of ‘private(variable, variable2)’ specifies that these variables should have a separate instance in each thread.
*/

#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)

// above and below axis 
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
  {printf(" %d from %d\n", iyAbove, iyAboveMax); //info 
  for(ix=ixMin; ix<=ixMax; ++ix)

  { // above axis compute color and save it to the array
    iy = iyAxisOfSymmetry + iyAbove;
    Color = GiveColor(ix, iy);
    PlotPoint(ix, iy, Color ); 
    // below the axis only copy Color the same as above without computing it 
    PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color ); 
   } 
}  
 return 0;
}

單元

"索尼 PlayStation 3 是消費市場上最便宜的平行計算機之一。"

GPGPU

• CUDA

  • 適用於 CUDA 的 Yellow Dog Linux

  • par4all
• OpenCl
• GLSL - GLSL 編輯工具
  • Fragmentarium : 用於使用基於 GPU（GLSL）的畫素圖形構建的 C++、OpenGL/GLSL 和 Qt 4 的遊樂場。
  • 著色器編輯器 使用 kickjs
  • GLSL 沙盒 由 Ricardo Cabello Miguel 使用 three.js
  • jsfiddle 和 示例
  • WebGl : 透過 html 使用 GLSL。它是一個基於 OpenGL ES 2.0 API 的 3D 圖形 API，開發人員可以使用它來建立完整的 3D 網路應用程式。
    • Inigo Quilez、Pol Jeremias 的 shadertoy : 一個基於瀏覽器的片段著色器編輯器，用於建立 2D 效果。
      • shadertoy 如何運作？
    • Chrome 實驗
    • WebGL 遊樂場 - 新的不工作！！
    • 示例
      • brainjam - webgl 螺旋
      • Evgeny Demidov 的 webgl 示例
      • 示例 由 Hvidtfeldt Christensen

Canvas : 支援瀏覽器的列表

效果

在 2 核 CPU 上的效果（每個核心 1 個執行緒）

方法	時間 [分鐘]	相對速度
CPU & 無最佳化	80	1
CPU & 對稱性	39	2
CPU & 對稱性 & OpenMP	24	3
GPU	? (待辦)

GPU 應該是

• 比 CPU 快 62 倍
• 快 760 倍

擾動演算法 與級數逼近（在 e100 附近的縮放級別快 100 倍！）

使用 WinXP Intel 2.9 GHz CPU（使用 1 個 CPU）和 GTX 480 GPU，我使用 1000x1000 圖和 1000 次最大迭代得到以下結果：

型別	時間	描述
gpu	0.07 秒	gpu 是 GPU 上的純 CUDA 解決方案
gpuarray	3.4 秒	gpuarray 在 GPU 上使用 Python 中類似 numpy 的 CUDA 包裝器
numpy	43.4 秒	numpy 是 CPU 上的純 Numpy（基於 C）解決方案
python（序列）	1605.6 秒	python 是 CPU 上的純 Python 解決方案，使用 numpy 陣列