經典的分形演算法
轉載自：http://www.douban.com/note/230496472/
被譽為大自然的幾何學的分形（Fractal）理論，是現代數學的一個新分支，但其本質卻是一種新的世界觀和方法論。它與動力系統的混沌理論交叉結合，相輔相成。它承認世界的區域性可能在一定條件下，在某一方面（形態，結構，資訊，功能，時間，能量等）表現出與整體的相似性，它承認空間維數的變化既可以是離散的也可以是連續的，因而拓展了視野。
分形幾何的概念是美籍法國數學家曼德布羅（B.B.Mandelbrot）1975年首先提出的，但最早的工作可追朔到1875年，德國數學家維爾斯特拉斯（K.Weierestrass）構造了處處連續但處處不可微的函式，集合論創始人康託（G.Cantor，德國數學家）構造了有許多奇異性質的三分康託集。1890年，義大利數學家皮亞諾（G.Peano）構造了填充空間的曲線。1904年，瑞典數學家科赫（H.von
Koch）設計出類似雪花和島嶼邊緣的一類曲線。1915年，波蘭數學家謝爾賓斯基（W.Sierpinski）設計了象地毯和海綿一樣的幾何圖形。這些都是為解決分析與拓樸學中的問題而提出的反例，但它們正是分形幾何思想的源泉。1910年，德國數學家豪斯道夫（F.Hausdorff）開始了奇異集合性質與量的研究，提出分數維概念。1928年布利幹（G.Bouligand）將閔可夫斯基容度應用於非整數維，由此能將螺線作很好的分類。1932年龐特里亞金（L.S.Pontryagin）等引入盒維數。1934年，貝塞考維奇（A.S.Besicovitch）更深刻地提示了豪斯道夫測度的性質和奇異集的分數維，他在豪斯道夫測度及其幾何的研究領域中作出了主要貢獻，從而產生了豪斯道夫-貝塞考維奇維數概念。以後，這一領域的研究工作沒有引起更多人的注意，先驅們的工作只是作為分析與拓撲學教科書中的反例而流傳開來。
真正令大眾瞭解分形是從計算機的普及肇始，而一開始，分形圖的計算機繪製也只是停留在二維平面，但這也足以使人們心馳神往。近來，一個分形體愛好者丹尼爾•懷特（英國一鋼琴教師）提出一個大膽的方法，創造出令人稱奇的3D分形影像，並將它們命名為芒德球（mandelbulb）。
 |
 |
在芒德球極其繁瑣的外表下，這個集合實際上是由一種非常基礎的演算法得出的。那是一種利用複數的演算法。就曼德布羅集而言，它是直接由最簡單的乘方運算得出的——對複數進行乘方。但問題在於無法在三維空間恰當地擴充套件數的概念。與複數和平面點之間的關係不同，19世紀的數學家們曾證明，立體空間中的點是無法用適宜傳統加法和乘法運算的代數工具來表示的。既然無法定義數字計算，自然也就無法勾畫曼德布羅集的三維形象。解決方案之一是在四維空間中進行計算，然後將結果投射到三維空間中。四維空間中的每個點都可與 “四元數”（quaternion）匹配，對它們可以進行傳統算術操作。儘管四維空間無法用肉眼看到，但利用四元數便能輕而易舉地列出與曼德布羅集相對應的演算法，之後去掉一個分量，就能使結果顯示成三維效果。但這個方案也令人失望，得到的畫面比二維影象好不了多少。
為了避開這個難題，丹尼爾•懷特兩年前冒出一個古怪的想法。徹底擺脫數學的羈絆，他在三維空間的點與點之間憑空構建出一種“偽分形”。儘管其處理手段算不上中規中矩的乘法，但至少將與曼德布羅集相對應的演算法擴充套件到了三維空間中所有的點。丹尼爾•懷特對幾百萬個點進行了計算，之後又追加了光影和紋理以體現立體效果，終於，在他的螢幕上呈現出第一個芒德球，形狀與嚴格的曼德布羅集十分近似。遺憾的是，這一結果沒能滿足他的期望：“圖形令人驚歎，但我期望的是更精緻的細節。”
嘗試並未就此止步。丹尼爾•懷特在網際網路上的一個分形體論壇上引起了美國一位年輕計算機專家保羅•尼蘭德的注意。他接手懷特的研究，對演算法進行稍事改動，把反覆的平方操作換成更高次方（八次方），從而得到了一系列新的芒德球，指數越高，細節就越豐富。
這個芒德球引起了我的極大興趣，下決心要學學分形體，於是乎決定從最簡單的分形演算法學起，希望與各位共勉。
以下開始介紹幾例最簡單的分形演算法：
一、Cantor三分集的遞迴演算法
選取一個歐氏長度的直線段，將該線段三等分，去掉中間一段，剩下兩段。將剩下的兩段分別再三等分，各去掉中間一段，剩下四段。將這樣的操作繼續下去，直到無窮，則可得到一個離散的點集。點數趨於無窮多，而歐氏長度趨於零。經無限操作，達到極限時所得到的離散點集稱之為Cantor集。
1.給定初始直線兩個端點的座標（ax，ay）和（bx，by），按Cantor三分集的生成規則計算出個關鍵點的座標如下：
cx=ax+(bx-ax)/3
cy=ay-d
dx=bx-(bx-ax)/3
dy=by-d
ay=ay-d
by=by-d
2.利用遞迴演算法，將計算出來的新點分別對應於（ax，ay）和（bx，by），然後利用步驟1的計算關係計算出下一級新點（cx，cy）和（dx，dy），並壓入堆疊。
3.給定一個小量c，當（bx，by）<c時，被壓入堆疊中的值依次釋放完畢，同時繪製直線段（ax，ay）-（bx，by），然後程式結束。
下面給出matlab程式：
function f=cantor(ax,ay,bx,by)
c=0.005;d=0.005;
if (bx-ax)>c
x=[ax,bx];y=[ay,by];hold on;
plot(x,y,'LineWidth',2);hold off;
cx=ax+(bx-ax)/3;
cy=ay-d;
dx=bx-(bx-ax)/3;
dy=by-d;
ay=ay-d;
by=by-d;
cantor(ax,ay,cx,cy);
cantor(dx,dy,bx,by);
end
執行cantor（0,5,5,5），出現圖例如下：
 |
二、Koch曲線的遞迴演算法
在一單位長度的線段上對其三等分，將中間段直線換成一個去掉底邊的等邊三角形，再在每條直線上重複以上操作，如此進行下去直到無窮，就得到分形曲線Koch曲線。
1.給定初始直線（ax，ay）-（bx，by），按Koch曲線的構成原理計算出各關鍵點座標如下：
cx=ax+(bx-ax)/3
cy=ay+(by-ay)/3
ex=bx-(bx-ax)/3
ey=by-(by-ay)/3
l=sqrt((ex-cx)^2+(ey-cy)^2)
alpha=atan((ey-cy)/(ex-cx))
dy=cy+sin(alpha+pi/3)*l
dx=cx+cos(alpha+pi/3)*l
2.利用遞迴演算法，將計算出來的新點分別對應於（ax，ay）和（bx，by），然後利用步驟1中的計算公式計算出下一級新點（cx，cy），（dx，dy），（ex，ey），並壓入堆疊。
3.給定一個小量c，當l<c時，被壓入堆疊中的值依次釋放完畢，同時繪製直線段（ax，ay）-（bx，by），然後結束程式。
下面給出matlab程式：
function f=Koch(ax,ay,bx,by,c)
if (bx-ax)^2+(by-ay)^2<c
x=[ax,bx];y=[ay,by];
plot(x,y);hold on;
else
cx=ax+(bx-ax)/3; cy=ay+(by-ay)/3;
ex=bx-(bx-ax)/3; ey=by-(by-ay)/3;
l=sqrt((ex-cx)^2+(ey-cy)^2);
alpha=atan((ey-cy)/(ex-cx));
if (alpha>=0&(ex-cx)<0)|(alpha<=0&(ex-cx)<0)
alpha=alpha+pi;
end
dy=cy+sin(alpha+pi/3)*l;
dx=cx+cos(alpha+pi/3)*l;
Koch(ax,ay,cx,cy,c);
Koch(ex,ey,bx,by,c);
Koch(cx,cy,dx,dy,c);
Koch(dx,dy,ex,ey,c);
end
執行Koch(0,0,100,0,10)，出現圖例如下：
 |
三、生成填充Julia集
1.設定引數a，b以及一個最大的迭代步數N。
2.設定一個限界值R，即實數R≧max（2，sqrt（a^2+b^2）。
3.對於平面上以R為半徑的圓盤內的每一點進行迭代，如果對於所有的n≦N，都有|x^2+y^2|≦R，那麼，在螢幕上繪製出相應的起始點，否則不繪製。
下面給出matlab程式：
a=-0.11;b=0.65;r=2;
for x0=-1:0.01:1
for y0=-1:0.01:1
x=x0;y=y0;
if x0^2+y0^2<1
for n=1:80
x1=x*x-y*y+a;
y1=2*x*y+b;
x=x1;
y=y1;
end
if (x*x+y*y)<r
plot(x0,y0);
end
hold on;
end
end
end
 |
 |
 |
 |
四、牛頓迭代
牛頓迭代是在數值求解非線性方程（組）的時候經常使用的方法。有些牛頓迭代能夠繪製出漂亮的圖形來，所以現在也常用於設計圖形。
Matlab程式如下：
首先編寫newton函式：
function y=newton(z)
if (z==0)
y=0;
return;
end
for i=1:1:2000
y=z-(z^3-1)/(3*z^2);
if (abs(y-z)<1.0e-7)
break;
end
z=y;
end
接著進入主程式：
clear all;clc;
A=1;B=0;C=1;
for a=-1:0.005:1
for b=-1:0.005:1
x0=a+b*i;
y=newton(x0);
if abs(y-A)<1.0e-6
plot(a,b,'r');hold on;
elseif abs(y-B)<1.0e-6
plot(a,b,'g');hold on;
elseif abs(y-C)<1.0e-6
plot(a,b,'y');hold on;
end
end
end
 |
五、迭代函式系IFS
IFS是分形的重要分支。它是分形影象處理中最富生命力而且最具有廣闊應用前景的領域之一。這一工作最早可以追溯到Hutchinson於1981年對自相似集的研究。美國科學家M.F.Barnsley於1985年發展了這一分形構型系統，並命名為迭代函式系統（Iterated Function System，IFS），後來又由Stephen Demko等人將其公式化，並引入到影象合成領域中。IFS將待生成的影象看做是由許多與整體相似的（自相似）或經過一定變換與整體相似的（自仿射）小塊拼貼而成。
演算法：
1.設定一個起始點（x0，y0）及總的迭代步數。
2.以概率P選取仿射變換W，形式為
X1=a x0+b y0 +e
Y1=c x0+d y0+f
3.以W作用點（x0，y0），得到新座標（x1，y1）。
4.令x0=x1，y0=y1。
5.在螢幕上打出（x0，y0）。
6.重返第2步，進行下一次迭代，直到迭代次數大於總步數為止。
下面給出一些IFS植物形態的matlab程式：
a=[0.195 -0.488 0.344 0.433 0.4431 0.2452 0.25;
0.462 0.414 -0.252 0.361 0.2511 0.5692 0.25;
-0.058 -0.07 0.453 -0.111 0.5976 0.0969 0.25;
-0.035 0.07 -0.469 -0.022 0.4884 0.5069 0.2;
-0.637 0 0 0.501 0.8562 0.2513 0.05];
x0=1;y0=1;
for i=1:10000
r=rand;
if r<=0.25
x1=a(1,1)*x0+a(1,2)*y0+a(1,5);
y1=a(1,3)*x0+a(1,4)*y0+a(1,6);
end
if r>0.25 & r<=0.5
x1=a(2,1)*x0+a(2,2)*y0+a(2,5);
y1=a(2,3)*x0+a(2,4)*y0+a(2,6);
end
if r>0.5 & r<=0.75
x1=a(3,1)*x0+a(3,2)*y0+a(3,5);
y1=a(3,3)*x0+a(3,4)*y0+a(3,6);
end
if r>0.75 & r<=0.95
x1=a(4,1)*x0+a(4,2)*y0+a(4,5);
y1=a(4,3)*x0+a(4,4)*y0+a(4,6);
end
if r>0.95 & r<=1
x1=a(5,1)*x0+a(5,2)*y0+a(5,5);
y1=a(5,3)*x0+a(5,4)*y0+a(5,6);
end
x0=x1;y0=y1;
plot(x1,y1);hold on;
end
得到圖例如下：
 |
修改部分系數便可得到另一種形態：
 |
六、三角形分形
function triangles(n);
clc;close all;
if nargin==0;
n=4;
end
rand('state',2);
C=rand(n+4,3);
figure;
axis square equal;hold on;
a=-pi/6;
p=0;
r=1;
[p,r,n,a]=tritri(p,r,n,a,C);
function [p,r,n,a]=tritri(p,r,n,a,C);
% 畫一個三角形
% p 是三角形中心
% r是三角形半徑
% n是遞迴次數
% a是三角形角度
% C是顏色矩陣
z=p+r*exp(i*([0:3]*pi*2/3+a));
zr=p+r*exp(i*([0:3]*pi*2/3+a))/2;
pf=fill(real(z),imag(z),C(n+2,:));
set(pf,'EdgeColor',C(n+2,:));
if n>0;
[p,r,n,a]=tritri(p,r/2,n-1,a+pi/3,C);
n=n+1;r=r*2;a=a-pi/3;
[zr(1),r,n,a]=tritri(zr(1),r/4,n-1,a,C);
n=n+1;r=r*4;
[zr(2),r,n,a]=tritri(zr(2),r/4,n-1,a,C);
n=n+1;r=r*4;
[zr(3),r,n,a]=tritri(zr(3),r/4,n-1,a,C);
n=n+1;r=r*4;
end
 |
七、曼德布羅集合
Mandelbrot set是在複平面上組成分形的點的集合。Mandelbrot集合可以用復二次多項式f(z)=z^2+c來定義。其中c是一個復引數。對於每一個c，從z=0開始對f(z)進行迭代序列 (0, f(0), f(f(0)), f(f(f(0))), .......)的值或者延伸到無限大，或者只停留在有限半徑的圓盤內。曼德布羅集合就是使以上序列不延伸至無限大的所有c點的集合。從數學上來講，曼德布羅集合是一個複數的集合。一個給定的複數c或者屬於曼德布羅集合M，或者不是。
1.設定引數a，b，以及一個最大的迭代步數N。
2.設定一個限界值R，不妨設實數R=2。
3.對於引數平面上的每一點c（a，b），使用以R為半徑的圓盤內的每一點進行迭代，如果對於所有的n≤N，都有|x*x+y*y|≤R*R，那麼，在螢幕上繪製出相應的起始點c（a，b），否則不繪製。
下面給出matlab程式：
r=4;%限界值
for a=-2:0.002:1
for b=-2:0.002:1%引數a，b取到一個範圍
x=a;y=b;%初始的複數c
for n=1:20
x1=x*x-y*y+a;%複數平方加一個c的運算
y1=2*x*y+b;
x=x1;%迭代
y=y1;
end
if(x*x+y*y)<r%限界
plot(a,b);
end
hold on;
end
end
 |
八、腦分形
作為IFS的一種應用
a=[0.03 0 0 0.45 0 0 0.05;
-0.03 0 0 -0.45 0 0.4 0.15;
0.56 -0.56 0.56 0.56 0 0.4 0.4;
0.56 0.56 -0.56 0.56 0 0.4 0.4];
x0=1;y0=1;
for i=1:100000
r=rand;
if r<=0.05
x1=a(1,1)*x0+a(1,2)*y0+a(1,5);
y1=a(1,3)*x0+a(1,4)*y0+a(1,6);
end
if r>0.05 & r<=0.2
x1=a(2,1)*x0+a(2,2)*y0+a(2,5);
y1=a(2,3)*x0+a(2,4)*y0+a(2,6);
end
if r>0.2 & r<=0.6
x1=a(3,1)*x0+a(3,2)*y0+a(3,5);
y1=a(3,3)*x0+a(3,4)*y0+a(3,6);
end
if r>0.6 & r<=1
x1=a(4,1)*x0+a(4,2)*y0+a(4,5);
y1=a(4,3)*x0+a(4,4)*y0+a(4,6);
end
x0=x1;y0=y1;
plot(x1,y1);
hold on;
end