分形算法与程序设计

第1

章初识分形1.1Fractal的含义1.2分形的几何特征1.3分形的度量1.4分形维数1.5分形是一种方法论1.6分形与计算机图形学1参考书：《分形算法与程序设计》1.1Fractal的含义

英文单词Fractal，在大陆被译为“分形”，在台湾被译为“碎形”。它是由美籍法国数学家曼德勃罗（BenoitMandelbrot）创造出来的。其含义是不规则的、破碎的、分数的。曼德勃罗是想用此词来描述自然界中传统欧几里得几何学所不能描述的一大类复杂无规的几何对象。

2参考书：《分形算法与程序设计》1.2分形的几何特征自相似性

自相似，便是局部与整体的相似。自仿射性

自仿射性是自相似性的一种拓展。如果，将自相似性看成是局部到整体在各个方向上的等比例变换的结果的话，那么，自仿射性就是局部到整体在不同方向上的不等比例变换的结果。前者称为自相似变换，后者称为自仿射变换。

精细结构

任意小局部总是包含细致的结构。3参考书：《分形算法与程序设计》1.3分形的度量（1）长度的测量

Length(n=0)=1

Length(n=1)=4/3

Length(n=2)=16/9…………Length=lim(Length(n))=lim(4/3)=∞n→∞n→∞n4参考书：《分形算法与程序设计》1.3分形的度量（2）面积的测量

Area(n0)=(1╳√3/6)/2=√3/12

Area(n1)=√3/12╳(4/9)Area(n2)=√3/12╳(4/9)2…………

Area(n)=lim(√3/12╳(4/9)n)=0n→∞如上所述，koch曲线在一维欧氏空间中的度量为∞，在二维欧氏空间中的面积为0。如此看来，Koch曲线在传统欧氏空间中不可度量。5参考书：《分形算法与程序设计》1.4分形维数

分形维数是分形的很好的不变量，它一般是分数，用它可以把握住分形体的基本特征。

图a是边长为1的正方形，当边长变成原来的1∕2时，原正方形中包含4个小正方形，如图b，而4=22；图c是边长为1的正立方体，当边长变成原来的1∕2时，原正立方体中包含8个小正立方体，如图d，而8=23。则有N＝kD,D=log(N)/log(k)这样Koch曲线的分形维数D=log(4)log(3)=1.26186参考书：《分形算法与程序设计》1.4分形维数

对于实际的自然景物，我们可以用计盒维数的方法测量分维。7参考书：《分形算法与程序设计》1.5分形是一种方法论

沃尔夫奖（WolfPrize）在颁发给分形理论创始人曼德勃罗时的评语所说的，“分形几何改变了我们对世界的看法”。

分形理论至少会在三个方面改变我们对世界的认识。首先，自然界中许多不规则的形态其背后都有规则，都可以用分形的方法建立模型并在计算机上构造出以假乱真的景象来，显然利用这套方法我们可以把世界压缩到几个分形规则中，便于携带和传播。其次，许多以前被认为是随机的现象，从分形理论的角度看并不是随机的，比如布朗运动、股票价格的波动、传染病的流行传播等，这为我们控制这些貌似随机的现象奠定了理论基础。最后，分形理论中的分数维概念，为我们认识世界中的复杂形态提供了一个新的尺度。复杂性科学是现代科学的前沿，在这门科学的研究过程中，发现了许多符合分形规则的复杂形态，而分数维是测量这些形态复杂程度的一种度量。也就是说，我们找到了对复杂性做定量分析的工具。8参考书：《分形算法与程序设计》1.6分形与计算机图形学

分形理论的发展离不开计算机图形学的支持，如果一个分形构造的表达，不用计算机的帮助是很难让人理解的。不仅如此，分形算法与现有计算机图形学的其他算法相结合，还会产生出非常美丽的图形，而且可以构造出复杂纹理和复杂形状，从而产生非常逼真的物质形态和视觉效果。分形作为一种方法，在图形学领域主要是利用迭代、递归等技术来实现某一具体的分形构造。分形几何学与计算机图形学相结合，将会产生一门新的学科——分形图形学。它的主要任务是以分形几何学为数学基础，构造非规则的几何图素，从而实现分形体的可视化，以及对自然景物的逼真模拟。9参考书：《分形算法与程序设计》第2

章分形图的递归算法2.1Cantor三分集的递归算法2.2Koch曲线的递归算法

2.3Sierpinski垫片的递归算法2.4Hilbert-Peano曲线的算法2.5分支结构分形递归算法2.6分形树递归算法10参考书：《分形算法与程序设计》递归算法u

直接递归调用的例子如下：voidRecur(n){

……Recur(m);……}过程Recur的内部又调用了自身——Recur过程。

11参考书：《分形算法与程序设计》递归算法u

间接递归调用的例子如下：voidRecur_A(n){……Recur_B(m);……}voidRecur_B(n){……Recur_A(m);……}

12参考书：《分形算法与程序设计》2.1Cantor三分集的递归算法算法：Cantor(ax,ay,bx,by)标题：Cantor三分集的递归算法参数：c（终止递归的小量）

d（不同层次线之间的距离）变量：ax,ay（曲线端点坐标）

bx,by（曲线端点坐标）

cx,xy

（曲线端点坐标）

dx,dy（曲线端点坐标）函数：plot(x1,y1)–(x2,y2)

（画直线函数）BEGINIF((bx-ax)<c)THENBEGIN

plot(ax,ay)-(bx,by)ENDELSEBEGIN

plot(ax,ay)-(bx,by)

cx=ax+(bx-ax)/3cy=ay+d

dx=bx-(bx-ax)/3

dy=by+day=ay+dby=by+d

cantor(ax,ay,cx,cy)

cantor(dx,dy,bx,by)ENDEND13参考书：《分形算法与程序设计》2.2Koch曲线的递归算法算法：Koch(ax,ay,bx,by,c)标题：Koch曲线的递归算法参数：c（终止递归的小量）

PI（π值）变量：ax,ay（线段端点坐标）

bx,by（线段端点坐标）

cx,xy

（线段端点坐标）

dx,dy（线段端点坐标）

ex,ey

（线段端点坐标）

L

（线段长度）

alpha（基线与水平线正方向夹角）函数：plot(x1,y1)–(x2,y2)

（画直线函数）

sin()

（正弦函数）

cos()

（余弦函数）

ArcTan()

（反正切函数）

sqrt()

（开平方函数）14参考书：《分形算法与程序设计》2.2Koch曲线的递归算法BEGIN

IF((bx-ax)*(bx-ax)+(by-ay)*(by-ay))<cTHEN

BEGIN

plot(ax,ay)-(bx,by)

END

ELSEBEGIN

cx=ax+(bx-ax)/3cy=ay+(by-ay)/3ex=bx-(bx-ax)/3

ey=by-(by-ay)/3L=sqrt((ex-cx)*(ex-cx)+(ey-cy)*(ey-cy))alpha=ArcTan((ey-cy)/(ex-cx))

IF((ex-cx)<0)THEN15参考书：《分形算法与程序设计》2.2Koch曲线的递归算法

BEGINalpha=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)ENDEND16参考书：《分形算法与程序设计》2.3Sierpinski垫片的递归算法算法：Sierpinski(x,y,L,n)标题：Sierpinski垫片递归算法参数：n（递归深度）变量：x,y

（三角形中心点坐标）

x1,y1

（三角形顶点坐标）

x2,y2（三角形顶点坐标）

x3,y3

（三角形顶点坐标）

x01,y01（小三角形中心点坐标）

x02,y02（小三角形中心点坐标）

x03,y03（小三角形中心点坐标）

L

（三角形的边长）函数：plot(x1,y1)–(x2,y2)

（画直线函数）

sin()

（正弦函数）

cos()

（余弦函数）

sqrt()

（开平方函数）17参考书：《分形算法与程序设计》2.3Sierpinski垫片的递归算法BEGINIF(n=1)THENBEGIN

x1=x-L/2y1=y+L*(sin(PI/6)/cos(PI/6))/2x2=x+L/2y2=y+L*(sin(PI/6)/cos(PI/6))/2x3=xy3=y-L*(sin(PI/6)/cos(PI/6))

plot(x1,y1)-(x1,y1)plot(x2,y2)-(x2,y2)plot(x3,y3)-(x3,y3)

ENDELSEBEGIN

x01=x-L/4y01=y+L*(sin(PI/6)/cos(PI/6))/4x02=x-L/4y02=y+L*(sin(PI/6)/cos(PI/6))/4x03=xy03=y-L*(sin(PI/6)/cos(PI/6))/2Sierpinski(x01,y01,L/2,n-1)Sierpinski(x02,y02,L/2,n-1)Sierpinski(x03,y03,L/2,n-1)ENDEND18参考书：《分形算法与程序设计》2.4Hilbert-Peano曲线的算法算法：Peano(n)标题：Hilbert-Peano曲线递归算法变量：n（递归深度）

len（线的长度）

x,y（端点坐标）函数：lineto(x,y)（画直线函数）过程：a(n)（基本元素构型）

b(n)（基本元素构型）

c(n)（基本元素构型）

d(n)（基本元素构型）

19参考书：《分形算法与程序设计》2.4Hilbert-Peano曲线的算法BEGIN

a(n)

BEGIN

IFn>0THEN

BEGINd(n-1)x=x-len

lineto(x,y)a(n-1)y=y+len

lineto(x,y)a(n-1)x=x+len

lineto(x,y)b(n-1)

END

END

b(n)

BEGIN

IFn>0THEN

BEGINc(n-1)y=y-len

lineto(x,y)b(n-1)x=x+len

lineto(x,y)b(n-1)y=y+len

lineto(x,y)a(n-1)

END

END

c(n)

BEGIN

IFn>0THEN

BEGINb(n-1)x=x-len

lineto(x,y)c(n-1)y=y+len

lineto(x,y)c(n-1)x=x+len

lineto(x,y)d(n-1)

END

END

d(n)

BEGIN

IFn>0THEN

BEGINa(n-1)y=y-len

lineto(x,y)d(n-1)x=x+len

lineto(x,y)d(n-1)y=y+len

lineto(x,y)c(n-1)

END

ENDEND20参考书：《分形算法与程序设计》2.5分支结构分形递归算法算法：Ramus

(x,y,alpha,L,n)标题：分支结构递归算法参数：PI（π值）变量：n（递归深度）

L（线段长度）

x,y

（线段起点坐标）

x1,y1（线段终点坐标）

alpha

（主干生成角度）

alpha_L（左支干生成角度）

alpha_R（右支干生成角度）函数：plot(x1,y1)-(x2,y2)

（画直线函数）

sin()（正弦函数）

cos()（余弦函数）21参考书：《分形算法与程序设计》2.5分支结构分形递归算法BEGIN

IF(n>0)THEN

BEGINx1=x+cos(alpha)*Ly1=y-sin(alpha)*Lplot(x,y)-(x1,y1)

alpha_L=alpha+PI/8

alpha_R=alpha-PI/8L=2*L/3Ramus(x2,y2,alpha_L,L,n-1)Ramus(x2,y2,alpha_R,L,n-1)

ENDEND22参考书：《分形算法与程序设计》2.6分形树递归算法算法：tree

(x,y,L,alpha)标题：分形树递归算法参数：PI（π值）

A（主干生长方向）

B（侧干与主干的夹角）

C（主干偏转角度）

s1（长度小量，控制递归深度）

s2（主干与侧干之比）

s3（上一级主干与下一级主干之比）变量：L

（主干的长度）

x,y

（主干起点坐标）

x1,y1（中间分叉点坐标）

x2,y2（主干终点坐标）

x1L,y1L（第一左分支终点坐标）

x1R,y1R（第一右分支终点坐标）

x2L,y2L（第二左分支终点坐标）

x2R,y2R（第二右分支终点坐标）函数：plot(x1,y1)-(x2,y2)

（画直线函数）

sin()（正弦函数）

cos()（余弦函数）23参考书：《分形算法与程序设计》2.6分形树递归算法BEGIN

IF(L>s1)THEN

BEGINx2=x+L*cos(A*PI/180)y2=y-L*sin(A*PI/180)x2L=x2+L/s2*cos((A+B)*PI/180)y2L=y2-L/s2*sin((A+B)*PI/180)x2R=x2+L/s2*cos((A-B)*PI/180)y2R=y2-L/s2*sin((A-B)*PI/180)x1=x+L/s2*cos(A*PI/180)y1=y-L/s2*sin(A*PI/180)x1L=x1+L/s2*cos((A+B)*PI/180)y1L=y1-L/s2*sin((A+B)*PI/180)x1R=x1+L/s2*cos((A-B)*PI/180)y1R=y1+L/s2*sin((A-B)*PI/180)plot(x,y)-(x2,y2)plot(x2,y2)-(x2L,y2L)plot(x2,y2)-(x2R,y2R)plot(x1,y1)-(x1R,y1R)plot(x1,y1)-(x1L,y1L)tree(x2,y2,L/s3,A-C)tree(x2L,y2L,L/s2,A+B)tree(x2R,y2R,L/s2,A-B)tree(x1R,y1R,L/s2,A-B)tree(x1L,y1L,L/s2,A-B)

ENDEND24参考书：《分形算法与程序设计》第3

章文法构图算法3.1LS文法3.2单一规则的LS文法生成

3.3多规则的LS文法生成3.4随机LS文法25参考书：《分形算法与程序设计》文法构图算法字母表：L，R生成规则：L→R，R→LR初始字母：R则有：R→LR→RLR→LRRLR→RLRLRRLR→LRRLRRLRLRRLR→……26参考书：《分形算法与程序设计》3.1LS文法二维LS是字母表的绘图规则如下：F：以当前方向前进一步，并画线；f：以当前方向前进一步，不画线；＋：逆时针旋转角；－：顺时针旋转角；[：将海龟当前信息压栈；]：将“[”时刻的海龟信息出栈。27参考书：《分形算法与程序设计》3.2单一规则的LS文法生成Koch曲线的LS文法如下：

w：Fa：60oP：F→F＋F－－F＋F

步骤0：F

步骤1：F＋

F－－F＋F

步骤2：F＋

F－－F＋

F＋

F＋

F－－F＋

F－－F＋

F－－F＋

F＋

F＋

F－－F＋

F

步骤3：F＋

F－－F＋

F＋

F＋

F－－F＋

F－－F＋

F－－F＋

F＋

F＋

F－－F＋

F＋

F＋

F－－F＋

F＋

F＋

F－－F＋

F－－F＋

F－－F＋

F＋

F＋

F－－F＋

F－－F＋F－－F＋F＋

F＋

F－－F＋

F－－F＋

F－－F＋

F＋

F＋

F－－F＋

F＋

F＋

F－－F＋

F＋

F＋

F－－F＋

F－－F＋

F－－F＋F＋

F＋

F－－F＋

F

步骤4：……

28参考书：《分形算法与程序设计》3.3单一规则的LS算法实现算法：LS标题：L_System生成算法参数：Axiom（公理）

n（字符串替换次数）变量：a[]（被替换字符）

x[]，y[]（替换字符串）

len（线段长度）

Aipha（基线与水平线夹角）

Dalta（生成角度）

ox,oy（起始坐标点）

n（递归深度）

bef（前一节点）

aft（后一节点）函数：lineto(x,y)（画直线函数）

GetLength(x)（字符串长度）BEGIN//初始化

x[0]=“F”a[1]=“F”x[1]=“F+F--F+F”

len=10

Aipha=0Delta=60ox=100

oy=400n=429参考书：《分形算法与程序设计》3.3单一规则的LS算法实现

y[]=x[0]

FORi=1TOn

yL=GetLength(y)j=0

WHILE(j<yL)

IF

y[j]==a[1]THENtree[]=tree[]+x[1]j=j+1

ELSEtree[]=tree[]+y[j]j=j+1

ENDIFENDWHILE

y[]=tree[]

tree.Empty()//清空

ENDFOR

tree[]=y[]

//按字符绘图

WHILE(i<xL)

SWITCH(tree[i])

CASE“F”

aft.x=bef.x+len*cos(Aipha*PI/180)

aft.y=bef.y-len*sin(Aipha*PI/180)

aft.d=bef.d//方向

LineTO(aft.x,aft.y)

bef=aft

BREAK

CASE“[”

stack[k]=befk=k+1

BREAK

CASE“]”

bef=stack[k-1]30参考书：《分形算法与程序设计》3.3单一规则的LS算法实现

k=k-1

LineTO(bef.x,bef.y)

BREAK

CASE“+”

bef=bef+DeltaBREAK

CASE“-”

bef=bef-Delta

BREAK

ENDSWITCHi=i+1

ENDWHILEEND31参考书：《分形算法与程序设计》3.3多规则的LS文法生成为了生成更复杂的图形，可将字母表增加字母元素Sierpinski垫片的LS文法如下：

w：La：600P1：L→＋R－L－R＋

P2：R→－L＋

R＋

L－相应的改造1：x[0]=“L”a[1]=“L”x[1]=“+R-L-R+”a[2]=“R”x[2]=“-L+R+L-”Delta=6032参考书：《分形算法与程序设计》3.3多规则的LS文法生成相应的改造2：

WHILE(j<xL)

IF

y[j]==a[1]THENtree[]=tree[]+x[1]j=j+1ELSEIFy[j]==a[2]THEN

tree[]=tree[]+x[2]j=j+1

ELSEtree[]=tree[]+y[i]j=j+1

ENDIFy[]=tree[]

ENDWHILE

33参考书：《分形算法与程序设计》3.3多规则的LS文法生成相应的改造3：

CASE“L”

aft.x=bef.x+len*cos(Aipha*PI/180)

aft.y=bef.y-len*sin(Aipha*PI/180)

aft.d=bef.d//方向

LineTO(aft.x,aft.y)

bef=aft

BREAK

CASE“R”

aft.x=bef.x+len*cos(Aipha*PI/180)

aft.y=bef.y-len*sin(Aipha*PI/180)

aft.d=bef.d//方向

LineTO(aft.x,aft.y)

bef=aft

BREAK34参考书：《分形算法与程序设计》3.4随机LS文法w:Fa:25P1：F————→F[＋F]F[－F]FP2：F————→F[＋F]F[－F[＋F]]P3：F————→FF[－F＋F＋F]＋[＋F－F－F]为了更好的模拟自然景物，需要随机使用多个变换规则。p1p2p3相应的改造1：x[0]=“F[+F]F[-F]F”x[1]=“F[+F]F[-F[+F]]”x[2]=“FF-[-F+F+F]+[+F-F-F]”相应的改造2：

WHILE(j<xL)

IF

y[j]==“F”THEN

r=Rand()%3tree[]=tree[]+x[r]j=j+1

ELSEtree[]=tree[]+y[r]j=j+1

ENDIFENDWHILE35参考书：《分形算法与程序设计》3.4随机LS文法36参考书：《分形算法与程序设计》第4

章迭代函数系统算法4.1混沌游戏4.2迭代函数系统4.3相似变换与仿射变换4.4IFS码4.5Sierpinski垫片的IFS生成4.6拼贴与IFS码的确定

4.7IFS植物形态实例4.8复平面上的IFS算法37参考书：《分形算法与程序设计》混沌游戏4.1给定平面上三点A,B,C。再任意给定初始点Z0,做下列迭代。ïîïíì+++=+,2/)(,2/)(,2/)(1CZBZAZZnnnn当掷出的硬币呈正面当掷出的硬币呈反面当掷出的硬币呈侧面38参考书：《分形算法与程序设计》迭代函数系统4.2迭代函数系统（IteratedFunctionSystem，IFS）是分形理论的重要分支。它将待生成的图像看成是由许多与整体相似的（自相似）或经过一定变换与整体相似的（自仿射）小块拼贴而成。39参考书：《分形算法与程序设计》相似变换与仿射变换直观上看：相似变换是指在各个方向上变换的比率必须相同的一种比例变换，仿射变换是指在不同的方向上变化的比率可以不同的一种比例变换。

4.3相似变换：如果对于任意两点A、B，以及对应点A’、B’，总有A’B’=k·AB（k为正实数），那么，这个变换叫做相似变换，实数k叫做相似比。仿射变换：x’=ax+by+ey’=cx+dy+f其中a,b,c,d,e,f为仿射变换系数。40参考书：《分形算法与程序设计》4.4IFS码用多个仿射变换式表达一个图象w1,w2,w3,……，使用每一个仿射变换式的概率p可以不同，一般面积越大，p值越大。于是，只要获得a,b,c,d,e,f,p（IFS码）的值便可以得到要表达的图形。x’=a1·x+b1·y+e1y’=c1·x+d1·y+f1w1x’=a2·x+b2·y+e2y’=c2·x+d2·y+f2w2x’=a3·x+b3·y+e3y’=c3·x+d3·y+f3w3p1p2p3p1+p2+p3=141参考书：《分形算法与程序设计》4.5Sierpinski垫片的IFS生成

由于生成的三个小三角形的面积相等，所以我们可以让w1、w2、w3出现的概率相同或相近。

x’=0.5xy’=0.5yx’=0.5x+0.5y’=0.5yx’=0.5x+0.25y’=0.5y-0.5x’=0.5·x+0·y+0y’=0·x+0.5·y+0x’=0.5·x+0·y+0.5y’=0·x+0.5·y+0x’=0.5·x+0·y+0.25y’=0·x+0.5·y-0.5w1w2w342参考书：《分形算法与程序设计》4.5Sierpinski垫片的IFS生成

算法：IFS标题：IFS生成算法参数：

n（迭代次数）

m[]（存储IFS码值）变量：a,b,c,d,e,f,p（IFS码变量）函数：Pset(x,y)（画点函数）

Rand（随机函数）BEGIN//IFS码赋值

m[0,0]=0.5;m[0,1]=0;m[0,2]=0m[0,3]=0;m[0,4]=0;m[0,5]=0m[0,6]=0.333m[1,0]=0.5;m[1,1]=0;m[1,2]=0m[1,3]=0.5;m[1,4]=0.5;m[1,5]=0m[1,6]=0.333m[2,0]=0.5;m[2,1]=0;m[2,2]=0m[2,3]=0.5;m[2,4]=0.25;m[2,5]=0.5m[2,6]=0.334n=6000

WHILEn>0p=Rand

SWITCH(P)

CASEIS<=m[0,6]a=m[0,0];b=m[0,1];c=m[0,2]d=m[0,3];e=m[0,4];f=m[0,5]

BREAK

CASEIS<=(m[0,6]+m[1,6])a=m[1,0];b=m[1,1];c=m[1,2]d=m[1,3];e=m[1,4];f=m[1,5]

BREAK43参考书：《分形算法与程序设计》4.5Sierpinski垫片的IFS生成

CASEIS<=(m[0,6]+m[1,6]+m[2,6])a=m[2,0];b=m[2,1];c=m[2,2]d=m[2,3];e=m[2,4];f=m[2,5]BREAKCASEIS<=(m[0,6]+m[1,6]+m[2,6]+m[3,6])a=m[3,0];b=m[3,1];c=m[3,2]d=m[3,3];e=m[3,4];f=m[3,5]

ENDSWITCH

x=a*x+b*y+ex=c*x+d*y+fPset(100+600*x,500-600*y)n=n-1

ENDWHILEEND（源代码：书中程序4.1）44参考书：《分形算法与程序设计》4.6拼贴与IFS码的确定

此时四个子图分别是目标图的1/4、1/5、1/4、1/2大小的复制品。然后按顺序交互式地在屏幕上调节每一个子图的仿射变换参数ai、bi、ci、di、ei，使得平移、旋转后基本覆盖住目标图。

45参考书：《分形算法与程序设计》4.7IFS植物形态实例

IFS码在书中表4.18IFS码在书中表4.20IFS码在书中表4.1946参考书：《分形算法与程序设计》4.8复平面上的IFS算法

47参考书：《分形算法与程序设计》4.8复平面上的IFS算法

算法：Julia_IFS标题：Julia集的IFS算法参数：z（迭代次数）

PI（π值）

RAND_MAX（随机最大值）变量：k（概率变量）

x,y（z的实部和虚部）

cx,cy（c的实部和虚部）

r（极距）

a,b,e,f,m,n,wx,wy,theta

函数：Pset(x,y)（画点函数）

Rand（随机函数）

sin（正弦函数）

cos（余弦函数）

atan（反正切函数）BEGINFORi=1TOzm=a*x+en=b*y+fIFi>10THEN

Pset(m,n)ENDIF

wx=x-cx

wy=y-cyIFwx>0THENtheta=atan(wy/wx)ENDIFIFwx<0THENtheta=PI+atan(wy/wx)ENDIFIFwx=0THENtheta=PI/2ENDIF48参考书：《分形算法与程序设计》4.8复平面上的IFS算法

theta=theta/2r=sqit(wx*wx+wy*wy)k=rand

rnd=k/RAND_MAXIFrnd<0.5THENr=sqrt(r)ELSEr=-sqrt(r)ENDIFx=r*cos(theta)y=r*sin(theta)ENDFOREND49参考书：《分形算法与程序设计》第5

章逃逸时间算法5.1基本思想

5.2Julia集的逃逸时间算法

5.3Mandelbrot集的逃逸时间算法5.4基于牛顿迭代的Julia集的逃逸时间算法50参考书：《分形算法与程序设计》逃逸时间算法的基本思想F(z)=z2+c当c=0时，由于z是复数，即z=x+yi，则有z2=z×z=(x+yi)×(x+yi)=x2+y2i2+2xyi=(x2-y2)+(2xy)i设复数z=x+yi的绝对值，即

|z|=SQR(x2+y2)|F(z0)|=|x02-y02+2x0y0i|=SQR((x02-y02)2+(2x0y0)2)=SQR(x04+y04-2x02y02+4x02y02)=SQR((x02+y02)2)=|z0|2若0<|z0|<1,|F(z0)|<|z0|，对于每一次迭代z趋向0，即z向0收敛。若|z0|>1,经过迭代z会趋向无穷，z向无穷逃逸。若|z0|>1,

z是平面上的单位圆。5.151参考书：《分形算法与程序设计》逃逸时间算法的基本思想当c≠0时，其吸引子不再是0，而是一个区域，称混沌区。如图，假设有一个充分大的整数N，当未逃逸区域M中的初始点a经过小于N次迭代就达到未逃逸区域M的边界，甚至超出了边界，我们就认为点a逃逸出去了；而如果经过N次迭代后a的轨迹仍未到达M的边界，我们就认为，a是A上的点。用这样的方法，描绘出A的边界图形，这便是逃逸时间算法的基本思想。

5.152参考书：《分形算法与程序设计》5.2Julia集的逃逸时间算法

53参考书：《分形算法与程序设计》5.2Julia集的逃逸时间算法

算法：Julia标题：Julia集的逃逸时间算法参数：K（逃逸时间）

m（逃逸半径）

Mx,My（绘图范围）

xs,xl,ys,yl（窗口范围）

p,q（复平面上C的坐标）变量：x0,y0（坐标变量）

r（膜变量）函数：SetP(x,y,color)（画点函数）

BEGIN//初始化

K=100m=500

Mx=800My=600

xs=-1.5xl=1.5

ys=-1.5

yl=1.5p=0.23q=0.043

xb=(xl-xs)/Mx

yb=(yl-ys)/My

54参考书：《分形算法与程序设计》5.2Julia集的逃逸时间算法

FOR

nx=0TO

Mx

FOR

ny=0TOMyx0=xs+nx*xby0=ys+ny*ybk=0Loop1:

xk=x0*x0-y0*y0+p

yk=2*x0*y0+qk=k+1r=xk*xk+yk*ykx0=xky0=yk

IFr>mTHENH=k;

gotoLoop2

ENDIF

IFk==KTHENH=r

GOTOLoop2

ENDIF

IFr<=m&&k<KTHEN

GOTOLoop1

ENDIF

Loop2:

SetP(nx,ny,H)

ENDFOR

ENDFOREND55参考书：《分形算法与程序设计》56参考书：《分形算法与程序设计》5.3Mandelbrot集的逃逸时间算法

算法：Mandelbrot标题：Mandelbrot集的逃逸时间算法参数：K（逃逸时间）

m（逃逸半径）

Mx,My（绘图范围）

xs,xl,ys,yl（窗口范围）

p,q（复平面上C的坐标）变量：x0,y0（坐标变量）

r（膜变量）函数：SetP(x,y,color)（画点函数）

BEGINpl=0.9

ps=2.3

ql=1.2

qs=-1.2K=100m=500

Mx=800My=600p=(pl-ps)/Mxq=(ql-qs)/My

FOR

np=0TO

Mx

FOR

nq=0TOMyp0=ps+np*pq0=qs+nq*qk=057参考书：《分形算法与程序设计》5.3Mandelbrot集的逃逸时间算法

x0=0y0=0Loop1:

xk=x0*x0-y0*y0+p0

yk=2*x0*y0+q0k=k+1r=xk*xk+yk*ykx0=xky0=yk

IFr>mTHENH=kGOTOLoop2

ENDIF

IFr<=mANDk<KTHEN

GOTOLoop1

ENDIFLoop2:

SetP(np,nq,H)

ENDFOR

ENDFOREND58参考书：《分形算法与程序设计》59参考书：《分形算法与程序设计》5.4基于牛顿迭代的Julia集的逃逸时间算法

牛顿迭代法求根公式：zn+1=zn-f(zn)/f’(zn)

其中，f’(zn)是f(zn)的导数。考虑f(z)=z3-1=0的情况,那么相应的牛顿变换是f(z)=(2z3+1)/3z2则z的三个根分别是w1=1,w2=ei2π/3,w3=ei4π/3，三个根的吸引域A(w1)，A(w2)，A(w3)的交界便是牛顿函数的Julia集。经过迭代，在A(wi)上的点都会被吸引到点wi上。设一个较大的迭代次数N，以及一个距离小量r，当迭代次数达到N，其与根点的距离小于r的被认为是收敛到某根上了，否则被认为是逃逸了。60参考书：《分形算法与程序设计》5.4基于牛顿迭代的Julia集的逃逸时间算法

算法：Julia_Newton标题：Julia集的牛顿迭代逃逸时间算法参数：N（迭代次数）

r（距离小量）

Mx,My（绘图范围）

xs,xl,ys,yl（窗口范围）

p,q（复平面上C的坐标）变量：x0,y0（坐标变量）

r（膜变量）函数：SetP(x,y,color)（画点函数）

BEGIN

fnm(x,y)=x*x+y*y

FORi=-150TO150

FORj=-150TO150x=i/75y=i/75

FORk=1TOn

IFi=0ANDj=0THEN

GOTOLoop

ENDIF

IFfnm(x-1,y)<rTHENSetP(i,j)

GOTOLoop

ENDIFfm=3*fnm(x,y)*fnm(x,y)61参考书：《分形算法与程序设计》5.4基于牛顿迭代的Julia集的逃逸时间算法

nx=2*x/3+(x*x-y*y)/fm

ny=2*y/3-2*x*y/fmx=nxy=ny

ENDFORLoop

ENDFOR

ENDFOREND62参考书：《分形算法与程序设计》第6

章分形显微镜6.1逃逸时间算法的放缩原理6.2Mandelbrot集的局部放大6.3Julia集的局部放大

6.4牛顿迭代法的局部放大6.5作为Julia集字典的Mandelbrot集63参考书：《分形算法与程序设计》逃逸时间算法的放缩原理6.164参考书：《分形算法与程序设计》6.2Mandelbrot集的局部放大

BEGIN

dp=(pmax-pmin)/(w_right-w_left)

dq=(qmax-qmin)/(w_bottom-w_top)FORi=w_leftTOw_rightp=pmin+dp*i;FORj=w_topTOw_bottomx0=0y0=0q=qmin+dq*jFORk=0TOm_times x=x0*x0-y0*y0+py=2*x0*y0+qr=x*x+y*yx0=xy0=y算法：Mandelbrot_zoom标题：Mandelbrot集的局部放大参数：m_times（逃逸时间）

m_deep（逃逸半径）变量：pmin，qmin

（参数窗口左上角坐标）

pmax，qmax

（参数窗口右下角坐标）

w_left，w_top

（绘图窗口左上角坐标）

w_right，w_bottom

（绘图窗口右下角坐标）

p,q（参数坐标）

x,y（绘图坐标）

x0,y0（绘图坐标）

r（极距）函数：SPost(x,y,color)

（画点）65参考书：《分形算法与程序设计》6.2Mandelbrot集的局部放大

IF(r>m_deep)THENbreakIF(k>=m_times)THENk=r

SPost(i,j,k);ENDFORENDFORENDFOREND 66参考书：《分形算法与程序设计》6.3Julia集的局部放大BEGIN

dp=(pmax-pmin)/(w_right-w_left)

dq=(qmax-qmin)/(w_bottom-w_top)FORi=w_leftTOw_rightp=pmin+dp*i;FORj=w_topTOw_bottomx0=pmin+dp*i;y0=qmin+dq*j;q=qmin+dq*jFORk=0TOm_times x=x0*x0-y0*y0+py=2*x0*y0+qr=x*x+y*yx0=xy0=y算法：Julia_zoom标题：Julia集的局部放大参数：m_times（逃逸时间）

m_deep（逃逸半径）变量：pmin，qmin

（参数窗口左上角坐标）

pmax，qmax

（参数窗口右下角坐标）

w_left，w_top

（绘图窗口左上角坐标）

w_right，w_bottom

（绘图窗口右下角坐标）

p,q（参数坐标）

x,y（绘图坐标）

x0,y0（绘图坐标）

r（极距）函数：SPost(x,y,color)

（画点）67参考书：《分形算法与程序设计》6.3Julia集的局部放大

IF(r>m_deep)THENbreakIF(k>=m_times)THENk=r

SPost(i,j,k);ENDFORENDFORENDFOREND 68参考书：《分形算法与程序设计》6.4牛顿迭代法的局部放大BEGIN

dp=(pmax-pmin)/(w_right-w_left)

dq=(qmax-qmin)/(w_bottom-w_top)FORi=w_leftTOw_rightFORj=w_topTOw_bottomx=pmin+dp*i;y=qmin+dq*j;count=0;WHILE(count<N)xx=x*x

yy=y*y d=3*((xx-yy)*(xx-yy)+4*xx*yy)

tmp=x; x=(2/3)*x+(xx-yy)/dy=(2/3)*y-2*tmp*y/d count=count+1; ENDWHILE算法：Newton_zoom标题：Newton法的局部放大参数：m_times（逃逸时间）

m_deep（逃逸半径）变量：pmin，qmin

（参数窗口左上角坐标）

pmax，qmax

（参数窗口右下角坐标）

w_left，w_top

（绘图窗口左上角坐标）

w_right，w_bottom

（绘图窗口右下角坐标）

p,q（参数坐标）

x,y（绘图坐标）

xx,yy（绘图坐标）

d（极距）函数：SPost(x,y,color)

（画点）69参考书：《分形算法与程序设计》6.4牛顿迭代法的局部放大IF(x>0)THENcolor=countELSEIF(x<-0.3)AND(y>0.0)THENcolor=count+100ELSEcolor=count+200ENDIFENDFORENDFOREND 70参考书：《分形算法与程序设计》6.5作为Julia集字典的Mandelbrot集

Julia集是一个固定的值的图形展现，而Mandelbrot集却要走遍所有的值，显然，每一个值都对应一个Julia集，所以我们说Mandelbrot集是Julia集微缩字典。我们可以在Mandelbrot集的绘图空间中任意取一点，并将其还原成相应的值，再将此值作为Julia集程序中选定的值，来绘制Julia集的图形。dp=(pl-ps)/Mxdq=(ql-qs)/MyDrawJulia(ps+dp*point.x,qs+dq*point.y)71参考书：《分形算法与程序设计》第7

章分形演化算法7.1从逻辑运算谈起7.2一维元胞自动机7.3二维元胞自动机7.4分形演化的DLA模型7.5用DLA模型模拟植物的生长7.6不同初始条件的DLA形态72参考书：《分形算法与程序设计》73参考书：《分形算法与程序设计》从逻辑运算谈起7.1逻辑异或

本行：00011011下一行：0110

74参考书：《分形算法与程序设计》7.2一维元胞自动机元胞按等间隔方式分布在一条向两侧无限延伸的直线中，称为一维元胞自动机。

本行：001