Science 版 (精华区)

发信人: qpcwth (独翅鸟), 信区: Science
标  题: 《分形艺术》41
发信站: 哈工大紫丁香 (2001年11月03日18:15:38 星期六), 站内信件

第五章 林氏系统与迭代函数系统
5.4 迭代函数系统方法
　
    魏尔(C.H.H.Weyl,1885-1955)在《对称》一书中就指出过，“Turritella duplica
ta ”的壳与“鹦鹉螺”的壳可用简单的伸缩变换刻划。不过，在80年代以前，只有少数
杰 出人物洞察到复杂的生物形态背后可能隐藏着简单的规则。分形几何学兴起后，通过
分形迭 代法模拟生物形态发生过程，已成为倍受科学界注目的一件事情。
    分形模拟自然形态所涉及的主要学科包括计算机图形学、仿射几何、形式语言等，
它们也都 因此而突然变得活跃起来。科学家早已可以轻而易举、维妙维肖地生成“人工
生物”了。只 要拥有一台286以上的微机，经过努力人人都可以尝试造一些生物，一方
面模拟已有的物种 ，另一方面可以任意创造新的物种。你的绘画技巧可能不甚高明，但
计算机可以助你一臂之 力，只需调系数、参数就可以生成各种不同的生物形态，其精细
程度不亚于素描大师的作品 。
    迭代函数系统(IFS，简称迭代函数系)方法是美国佐治亚理工学院的巴恩斯利教授首
创的。I FS方法的魅力在于它是分形迭代生成的“反问题”，根据拼接定理(collage t
heorem)，对 于一个给定的图形(比如一幅照片)，求得几个生成规则，就可以大幅度压
缩信息。
    这里采用确定性算法与随机性算法相结合的办法生成植物杆茎或叶片。“确定性”
指用以迭 代的规则是确定性的，它们由一组仿射变换(如R_1，R_2，R_3等等) 构成；“
随机性”指迭代过程是不确定的，每一次究竟迭代哪一个规则，即R_i中具体哪一 个，
不是预先定好的，而要靠掷骰子的办法来决定。设最终要生成的图形(植物形态图)为M 
，它要满足下述集合方程：
　
M=R_1∪R_2∪…∪R_N
　
上式的含义是，随机地从R_i(i=1,…，N)中挑选一个迭代规则迭代 一次，然后再随机地
在R_i(i=1,…，N)中选一个规则迭代一次，不 断重复此过程，最后生成的极限图形M就
是欲求的植物形态图。
    每个迭代规则R_i都是一个仿射变换。正交变换保持几何图形的度量性质(向量的夹
 角，点与点之间的距离，图形的面积等)不变；而仿射变换一般会改变图形中向量的夹
角、 点对间的距离、图形的面积等，但仿射变换不改变共线、平行、相交、共线点的顺
序、中心 对称、二次曲线的次数等。举例来说，原来平行的线段，经仿射变换后仍然是
平行的，但长 度可能发生了变化，之间的距离也可能改变了。一个图形经仿射变换后面
积变小了，则此变 换是收缩的，如果变大了则是扩张的，若保持不变则是恒等的。这里
只用到收缩性仿射变换 ，因为极限图形M应当是所有迭代R_i的吸引子，每个仿射变换是
收缩性的才 能保证迭代收敛到M上。
    仿射变换是一种线性变换，在二维平面上进行讨论，二维仿射变换的形式为：
R：
x′=ax+by+e
y′=cx+dy+f
上式中未知数有两个：x和y；系数有四个：a,b,c和d；常数有两个 ：e和f。
　
    设平面上一有面积区域为D，经仿射变换R变换后变成D′，则D与 D′的面积关系为
：
　
S(D′)=｜det(A)｜·S(D)
　
这里｜det(A)｜表示矩阵A的行列式的绝对值，它是小于1的正数，其中A 为系数矩阵。

对于不同的R_i(i=1，2，…，N)，有相应的A_i(由 a_i，b_i，c_i，d_i组成)，相应的
常数e_i和 f_i。通常N取2，3，4，有时高达16。
    剩下的一个问题是怎样实现掷骰子操作。在计算机上可以很容易地用伪随机数发生
器不断生 成随机数，用来代替人工抛掷有N个面的骰子。不失一般性，设N=4，每次用 
计算机生成一个随机数E∈(0，100)，设0＜β_1＜β_2＜β _3＜100，作如下规定：
若0＜E＜β_1，则选择规则R_1，
若β_1≤E＜β_2，则选择规则R_2，
若β_2≤E＜β_3，则选择规则R_3，
若β_3≤E＜100，则选择规则R_4.
指定β_i的过程，相当于为每一种迭代规则R_i指派一个概率p_i。 虽说具体每一次用哪
个规则迭代不确定，但从长期行为看，每种规则使用的频率是一定的。 当然这种概率可
以随意指定，但是为了使得迭代高效、迅速，概率的大小应与仿射变换的图 形压缩率成
正比。图形压缩率就是前文说的｜det(A)｜，对于N=4的情况， 它有4个。
    至此可以知道，每一种规则R_i都包括7个参数(或系数)：a_i，b_i， c_i，d_i，e
_i，f_i，p_i。要求p_1+ p_2+…+p_N=1。如果生成某一种植物形态需要4种规则，则共
有 4×7=28个参数。只要给定了这些参数，实际上极限图形M就完全确定了。开始迭代 
时可能还一时看不出M的面貌，但过了一会以后，M的轮廓就可以看清楚了， 最后M的精
细结构就展示出来了。浙江大学数学系陈刚(1963- )最近还研究了准仿射 变换的一些性
质。用准仿射变换也能生成各种分形图形。
    图5.6(上左)就是用IFS方法生成的图象，此图象的算法由美国佐治亚工学院的巴恩
斯利 首创，它像不像蕨类植物的叶子?生成此叶子用了4个仿射变换R_1，R_2， R_3和R
_4。7个参数的取值见表5.3。
表5.3 巴恩斯利蕨的参数表
R a b c d e f p
1 0 0 0 0.16 0 0 0.01
2 0.85 0.04 －0.04 0.85 0 1.6 0.85
3 0.2 －0.26 0.23 0.22 0 1.6 0.07
4 －0.15 0.28 0.26 0.24 0 0.44 0.07
　
表5.4 树叶B的参数表
R a b c d e f p
1 0 0 0 0.5 0 0 0.05
2 0.12 -0.82 0.42 0.42 0 0.2 0.4
3 0.12 0.82 -0.42 0.42 0 0.2 0.4
4 0.1 0 0 0.1 0 0.2 0.15
　
增减规则R_i，可以改变最终植物M的形态。即使不改变迭代规则，采用同样的程序 ，只
改变参数也可以生成完全不同的植物形态。为了使程序具有通用性，可以设计如下“过
 程”：
　
Procedure AFF(a,b,c,d,e,f,S,T:real);
Var lins:real;
Begin
  lins:=a*S+b*T+e;
  y:=c*S+d*y+f;
  x:=lins;
End;
　
    程序中不断调用此过程AFF，即可完成迭代。这样，每次修改参数，只需改变数组中
的数据 就可以了。由此得到一个启示：外表极不相同、极复杂的形态，其背后控制其生
长的规律可 能是相同(或相似)的、简单的。再推一步，分子水平上控制形态发生的DNA
遗传编码，也应 当是简单的。这一猜想有许多合理之处，直观上想象一下，DNA分子不
可能承载过多的形态 控制信息。
　
图5.6 用IFS方法生成的植物形态
    对于表5.3，如果把其中的第二行中表项(a,b,f)的值由(0.85，0.04,1.6)改 为(-0
.85, 0.09,2.6)，则会生成蒿草。对于表5.4，若把第二、第三行的表项c分别由0.42 改
为0.2(保留符号)，则会生成六角枫叶；若把(a,b)两列的0.12和0.82都改 成0.42(保留
符号)，则会生成树冠。读者也可以自己改变参数，试验着生成各种精美的分形 植物形
态。
　
｛IFS生成植物形态PASCAL程序｝
Program FractalPlant28;
Uses Graph,Crt;
var
x,y,E,u:real;
Gd,Gm:integer;
Begin
  Gd:=VGA;｛用VGA方式640×480｝
  Gm:=VGAHi;
  Initgraph(Gd,Gm,’D:\\PASCAL’);｛初始化图形｝
  if GraphResult <> grOK then Halt(1);｛若出错则停机｝
  x:=0;y:=0;｛赋初值｝
  Randomize;｛初始化随机数发生器｝
  Repeat
    E:=Random(100);｛产生一个介于0至100之间的随机数｝
    if E < 1 then
    Begin
      x:=0;y:=0.16 * y;｛用R_1迭代｝
    End;
  if (E >=1) and (E < 86) then
    Begin
      u:=0.85*x+0.04*y;y:=-0.04*x+0.85*y+1.6;｛用R_2迭代｝
      x:=u;
    End;
  if (E >= 86) and (E < 97) then
    Begin
      u:=0.2*x-0.26*y;y:=0.23*x+0.22*y+1.6;｛用R_3迭代｝
      x:=u;
    End;
  if E >=97 then
    Begin
      u:=-0.15*x+0.28*y;y:=0.26*x+0.24*y+0.44;｛用R_4迭代｝
      x:=u;
    End;
  PutPixel(Round(50*x)+300,500-Round(50*y),2);｛描点，绿色｝
  Until KeyPressed;｛按任意键退出｝
  CloseGraph;｛关闭图形方式｝
End.
　
　
IFS的规则实际上可以自由设计，下面再给出一个稍复杂的例子。令Γ=(x, y)^T=R^2，
取S_i，i=1,2,3，为下列函数：
　
　
| 0   0 |           |   1/2   |
S_1(Γ) = |          | Γ  +  |           |
                   | 0   c |           |      0   |
                   |rcosφ  -rsinφ|          |1/2-r/2cosφ|
S_2(Γ) = |                           | Γ  + |                    |
                   |rsinφ  rcosφ |          |c-r/2sinφ    |
                   |qcosψ  -rsinψ |         |1/2-q/2cosψ |
S_3(Γ)=  |                             | Γ+  |                      |
                   |qsinψ   rcosψ |          |3/5c-q/2sinψ|
其中c=0.255, r=0.75, q=0.625,φ=-π/8，ψ= π/5, p_1=p_2=p_3=1/3.
用PASCAL实现这个IFS，则有如下程序：
　
Program FractalTree9611;
uses Graph,Crt;
var
   x,y,ran,u,aa,bb,mm,nn,c,r,q,phi,psi:real;
Gd,Gm:integer;
begin
  Gd:=VGA;
  Gm:=VGAHi;
  InitGraph(Gd,Gm,'D:\PASCAL');
  if GraphResult <> grOK then Halt(1);
  x:=0;
  y:=0;
  c:=0.255;
  r:=0.75;
  q:=0.625;
  phi:=-pi/8;
  psi:=pi/5;
  aa:=sin(phi);
  bb:=cos(phi);
  mm:=sin(psi);
  nn:=cos(psi);
  Randomize;
  Repeat
    ran:=Random(100);
    if ran < 33 then
       begin{ if p1}
         x:=c*x+2/3;{x=ax+by+e}
         y:=c*y;{y=cx+dy+f}
       end;
    if (ran >33) and (ran < 66) then
       begin
         u:=r*bb*x-r*aa*y+1/2-r/2*bb;
         y:=r*aa*x+r*bb*y+c-r/2*aa;
         x:=u;
       end;
    if ran > 66 then
       begin
         u:=q*nn*x-r*mm*y+1/2-q/2*nn;
         y:=q*mm*x+r*nn*y+3/5*c-q/2*mm;
         x:=u;
       end;
    PutPixel(round(445*x+50),400-round(445*y),2);
  until KeyPressed;
  CloseGraph;
end.
　
    可以修改函数S_i(Γ)，得到许多意想不到的分形树。在修改迭代函数的过 程中，
PutPixel语句中的放大倍数也应作适当调整。
通过IFS方法还可以绘制朱丽亚集J。设z→z^2+c，求出两个逆变换 ：
　
w_1(z)=SQRT(z-c)，w_2(z) =-SQRT(z-c)
, 取概率(p_1,p_2)=(1/2,1/2)，迭代后生成的吸引子实际上就是各种朱丽亚 集!您觉得
有些不可思议?是的，细想一下，IFS方法、L系统、混沌动力学以及有关 M集和J集的迭
代，事实上都是相互关联着的，这种关联体现了非线性科学的 内在逻辑。
    最后我们给出一个用IFS方法生成山脉地形的小程序Mount.PAS：
　
{Mount.PAS}
Program IFSMountain;
uses
Graph,Crt;
var
Gd,Gm,ErrorCode : integer;
x,y: real;
k,MaxY:integer;
i:longint;
d: array［1..4,1..6］ of real;
begin
Gd:= Detect;
InitGraph(Gd,Gm,'d:\pascal');
ErrorCode := GraphResult;
if ErrorCode <> grOk then exit;
MaxY := GetMaxY;
d［1,1］:=0.5; d［1,2］:=0; d［1,3］:=0;
d［1,4］:=0.5; d［1,5］:=0;d［1,6］:=0;
d［2,1］:=0.5;  d［2,2］:=0;  d［2,3］:=0;
d［2,4］:=0.5;  d［2,5］:=2; d［2,6］:=0;
d［3,1］:=-0.4;  d［3,2］:=0;  d［3,3］:=1;
d［3,4］:=0.4;  d［3,5］:=0; d［3,6］:=1;
d［4,1］:=-0.5;  d［4,2］:=0;  d［4,3］:=0;
d［4,4］:=0.5;  d［4,5］:=2; d［4,6］:=1;
randomize;
x := 0; y := 0;
repeat
i:=i+1;
k := random(4) + 1;
x := d［k,1］*x + d［k,2］*y + d［k,5］;
y := d［k,3］*x + d［k,4］*y + d［k,6］;
if i > 10 then
putpixel(round(MaxY*x/2),MaxY-round(MaxY*y/2),15)
until keypressed;
closegraph;
end.

--
心事浩茫连广宇,于无声处听惊雷

※ 来源:·哈工大紫丁香 bbs.hit.edu.cn·[FROM: 202.118.229.154]