Question

自然界的很多事物，例如树木、云彩、山脉、闪电、雪花以及海岸线等等都呈现出传统的几何学不能描述的形状。这些形状都有如下的特性：

• 有着十分精细的不规则的结构
• 整体与局部相似，例如一根树杈的形状和一棵树很像

分形几何学就是用来研究这样一类的几何形状的科学，借助计算机的高速计算和图像显示，使得我们可以更加深入地直观地观察分形几何。在本章中，让我们用 Python 绘制一些经典的分形图案。

Mandelbrot 集合

Mandelbrot(曼德布洛特) 集合是在复平面上组成分形的点的集合。

Mandelbrot 集合的定义（摘自维基百科)

Mandelbrot 集合可以用下面的复二次多项式定义：

其中 c 是一个复参数。对于每一个 c，从 z=0 开始对函数 进行迭代。

序列 的值或者延伸到无限大，或者只停留在有限半径的圆盘内。

Mandelbrot 集合就是使以上序列不发散的所有 c 点的集合。

从数学上来讲，Mandelbrot 集合是一个复数的集合。一个给定的复数 c 或者属于 Mandelbrot 集合，或者不是。

用程序绘制 Mandelbrot 集合时不能进行无限次迭代，最简单的方法是使用逃逸时间（迭代次数) 进行绘制，具体算法如下：

• 判断每次调用函数 得到的结果是否在半径 R 之内，即复数的模小于 R
• 记录下模大于 R 时的迭代次数
• 迭代最多进行 N 次
• 不同的迭代次数的点使用不同的颜色绘制

下面是完整的绘制 Mandelbrot 集合的程序：

# -*- coding: utf-8 -*-

import numpy as np
import pylab as pl
import time
from matplotlib import cm

def iter_point(c):
  z = c
  for i in xrange(1, 100): # 最多迭代 100 次
    if abs(z)>2: break # 半径大于 2 则认为逃逸
    z = z*z+c
  return i # 返回迭代次数

def draw_mandelbrot(cx, cy, d):
  """
 绘制点(cx, cy) 附近正负 d 的范围的 Mandelbrot
 """
  x0, x1, y0, y1 = cx-d, cx+d, cy-d, cy+d 
  y, x = np.ogrid[y0:y1:200j, x0:x1:200j]
  c = x + y*1j
  start = time.clock()
  mandelbrot = np.frompyfunc(iter_point,1,1)(c).astype(np.float)
  print "time=",time.clock() - start
  pl.imshow(mandelbrot, cmap=cm.Blues_r, extent=[x0,x1,y0,y1])
  pl.gca().set_axis_off()

x,y = 0.27322626, 0.595153338

pl.subplot(231)
draw_mandelbrot(-0.5,0,1.5)
for i in range(2,7):  
  pl.subplot(230+i)
  draw_mandelbrot(x, y, 0.2**(i-1))
pl.subplots_adjust(0.02, 0, 0.98, 1, 0.02, 0)
pl.show()

Mandelbrot 集合，以 5 倍的倍率放大点(0.273, 0.595) 附近

程序中的 iter_point 函数计算点 c 的逃逸时间，逃逸半径 R 为 2.0，最大迭代次数为 100。draw_mandelbrot 函数绘制以点(cx, cy) 为中心，边长为 2*d 的正方形区域内的 Mandelbrot 集合。

下面 3 行计算指定范围内的迭代公式的参数 c，c 是一个元素为复数的二维数组，大小为 200*200，注意 np.ogrid 不是函数：

x0, x1, y0, y1 = cx-d, cx+d, cy-d, cy+d
y, x = np.ogrid[y0:y1:200j, x0:x1:200j]
c = x + y*1j

下面一行程序通过调用 np.frompyfunc 将 iter_point 转换为 NumPy 的 ufunc 函数，这样它可以自动对 c 中的每个元素调用 iter_point 函数，由于结果的数组元素类型为 object，还需要调用 astype 方法将其元素类型转换为浮点类型：

mandelbrot = np.frompyfunc(iter_point,1,1)(c).astype(np.float)

最后调用 matplotlib 的 imshow 函数将结果数组绘制成图，通过 cmap 关键字参数指定图的值和颜色的映射表：

pl.imshow(mandelbrot, cmap=cm.Blues_r, extent=[x0,x1,y0,y1])

使用 Python 绘制 Mandelbrot 集合最大的问题就是运算速度太慢，下面是上面每幅图的计算时间：

time= 0.88162629608
time= 1.53712748408
time= 1.71502160191
time= 1.8691174437
time= 3.03812691278

因为计算每个点的逃逸时间均不相同，因此每幅图的计算时间也不相同。

计算速度慢的最大的原因是因为 iter_point 函数的运算速度慢，如果将此函数用 C 语言重写的话将能显著地提高计算速度，下面使用 scipy.weave 库将 C++重写的 iter_point 函数转换为 Python 能调用的函数：

import scipy.weave as weave

def weave_iter_point(c):
  code = """
 std::complex<double> z;
 int i;
 z = c;
 for(i=1;i<100;i++)
 {
 if(std::abs(z) > 2) break;
 z = z*z+c;
 }
 return_val=i;
 """

  f = weave.inline(code, ["c"], compiler="gcc")
  return f

下面是使用 weave_iter_point 函数计算 Mandelbrot 集合的时间：

time= 0.285266982256
time= 0.271430028118
time= 0.293769180161
time= 0.308515188383
time= 0.411168179196

通过 NumPy 的数组运算也可以提高计算速度，前面的计算都是先对复数平面上的每个点进行循环，然后再循环迭代计算每个点的逃逸时间。如果要用 NumPy 的数组运算加速计算的话，可以将这两个循环的顺序颠倒过来，下面的程序演示这一算法：

# -*- coding: utf-8 -*-

import numpy as np
import pylab as pl
import time
from matplotlib import cm

def draw_mandelbrot(cx, cy, d, N=200):
  """
 绘制点(cx, cy) 附近正负 d 的范围的 Mandelbrot
 """
  global mandelbrot

  x0, x1, y0, y1 = cx-d, cx+d, cy-d, cy+d 
  y, x = np.ogrid[y0:y1:N*1j, x0:x1:N*1j]
  c = x + y*1j

  # 创建 X,Y 轴的坐标数组
  ix, iy = np.mgrid[0:N,0:N]

  # 创建保存 mandelbrot 图的二维数组，缺省值为最大迭代次数
  mandelbrot = np.ones(c.shape, dtype=np.int)*100

  # 将数组都变成一维的
  ix.shape = -1
  iy.shape = -1
  c.shape = -1
  z = c.copy() # 从 c 开始迭代，因此开始的迭代次数为 1

  start = time.clock()

  for i in xrange(1,100):
    # 进行一次迭代
    z *= z
    z += c
    # 找到所有结果逃逸了的点
    tmp = np.abs(z) > 2.0
    # 将这些逃逸点的迭代次数赋值给 mandelbrot 图
    mandelbrot[ix[tmp], iy[tmp]] = i

    # 找到所有没有逃逸的点
    np.logical_not(tmp, tmp)
    # 更新 ix, iy, c, z 只包含没有逃逸的点
    ix,iy,c,z = ix[tmp], iy[tmp], c[tmp],z[tmp]
    if len(z) == 0: break

  print "time=",time.clock() - start
  pl.imshow(mandelbrot, cmap=cm.Blues_r, extent=[x0,x1,y0,y1])
  pl.gca().set_axis_off()

x,y = 0.27322626, 0.595153338

pl.subplot(231)
draw_mandelbrot(-0.5,0,1.5)
for i in range(2,7):  
  pl.subplot(230+i)
  draw_mandelbrot(x, y, 0.2**(i-1))
pl.subplots_adjust(0.02, 0, 0.98, 1, 0.02, 0)
pl.show()

为了减少计算次数，程序中每次迭代之后，都将已经逃逸的点剔除出去，这样就需要保存每个点的下标，程序中用 ix 和 iy 这两个数组来保存没有逃逸的点的下标，因为有额外的数组保存下标，因此数组 z 和 c 不需要是二维的。函数迭代部分的程序如下：

# 进行一次迭代
z *= z
z += c

使用 * =, += 这样的运算符能够让 NumPy 不分配额外的空间直接在数组 z 上进行运算。

下面的程序计算出逃逸点，tmp 是逃逸点在 z 中的下标，由于 z 和 ix 和 iy 等数组始终是同时更新的，因此 ix[tmp], iy[tmp]就是逃逸点在图像中的下标：

# 找到所有结果逃逸了的点
tmp = np.abs(z) > 2.0
# 将这些逃逸点的迭代次数赋值给 mandelbrot 图
mandelbrot[ix[tmp], iy[tmp]] = i

最后通过对 tmp 中的每个元素取逻辑反，更新所有没有逃逸的点的对应的 ix, iy, c, z：

# 找到所有没有逃逸的点
np.logical_not(tmp, tmp)
# 更新 ix, iy, c, z 只包含没有逃逸的点
ix,iy,c,z = ix[tmp], iy[tmp], c[tmp], z[tmp]

此程序的计算时间如下：

time= 0.186070576008
time= 0.327006365334
time= 0.372756034636
time= 0.410074464771
time= 0.681048289658
time= 0.878626752841

连续的逃逸时间

修改逃逸半径 R 和最大迭代次数 N，可以绘制出不同效果的 Mandelbrot 集合图案。但是前面所述的方法计算出的逃逸时间是大于逃逸半径时的迭代次数，因此所输出的图像最多只有 N 种不同的颜色值，有很强的梯度感。为了在不同的梯度之间进行渐变处理，使用下面的公式进行逃逸时间计算：

是迭代 n 次之后的结果，通过在逃逸时间的计算中引入迭代结果的模值，结果将不再是整数，而是平滑渐变的。

下面是计算此逃逸时间的程序：

def smooth_iter_point(c):
  z = c
  for i in xrange(1, iter_num):
    if abs(z)>escape_radius: break
    z = z*z+c
  absz = abs(z)
  if absz > 2.0:
    mu = i - log(log(abs(z),2),2)
  else:
    mu = i
  return mu # 返回正规化的迭代次数

如果你的逃逸半径设置得很小，例如 2.0，那么有可能结果不够平滑，这时可以在迭代循环之后添加几次迭代保证 z 能够足够逃逸，例如：

z = z*z+c
z = z*z+c
i += 2

下图是逃逸半径为 10，最大迭代次数为 20 时，绘制的结果：

逃逸半径=10，最大迭代次数=20 的平滑处理后的 Mandelbrot 集合

逃逸时间公式是如何得出的？

请参考： http://linas.org/art-gallery/escape/ray.html

完整的程序请参考绘制 Mandelbrot 集合

迭代函数系统(IFS)

迭代函数系统是一种用来创建分形图案的算法，它所创建的分形图永远是绝对自相似的。下面我们直接通过绘制一种蕨类植物的叶子来说明迭代函数系统的算法：

有下面 4 个线性函数将二维平面上的坐标进行线性映射变换：

1.
  x(n+1）= 0
  y(n+1) = 0.16 * y(n)

2.
  x(n+1) = 0.2 * x(n) − 0.26 * y(n)
  y(n+1) = 0.23 * x(n) + 0.22 * y(n) + 1.6

3.
  x(n+1) = −0.15 * x(n) + 0.28 * y(n)
  y(n+1) = 0.26 * x(n) + 0.24 * y(n) + 0.44

4.
  x(n+1) = 0.85 * x(n) + 0.04 * y(n)
  y(n+1) = −0.04 * x(n) + 0.85 * y(n) + 1.6

所谓迭代函数是指将函数的输出再次当作输入进行迭代计算，因此上面公式都是通过坐标 x(n),y(n) 计算变换后的坐标 x(n+1),y(n+1)。现在的问题是有 4 个迭代函数，迭代时选择哪个函数进行计算呢？我们为每个函数指定一个概率值，它们依次为 1%, 7%, 7%和 85%。选择迭代函数时使用通过每个函数的概率随机选择一个函数进行迭代。上面的例子中，第四个函数被选择迭代的概率最高。

最后我们从坐标原点(0,0) 开始迭代，将每次迭代所得到的坐标绘制成图，就得到了叶子的分形图案。下面的程序演示这一计算过程：

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as pl
import time

# 蕨类植物叶子的迭代函数和其概率值
eq1 = np.array([[0,0,0],[0,0.16,0]])
p1 = 0.01

eq2 = np.array([[0.2,-0.26,0],[0.23,0.22,1.6]])
p2 = 0.07

eq3 = np.array([[-0.15, 0.28, 0],[0.26,0.24,0.44]])
p3 = 0.07

eq4 = np.array([[0.85, 0.04, 0],[-0.04, 0.85, 1.6]])
p4 = 0.85

def ifs(p, eq, init, n):
  """
 进行函数迭代
 p: 每个函数的选择概率列表
 eq: 迭代函数列表
 init: 迭代初始点
 n: 迭代次数

 返回值： 每次迭代所得的 X 坐标数组， Y 坐标数组， 计算所用的函数下标 
 """

  # 迭代向量的初始化
  pos = np.ones(3, dtype=np.float)
  pos[:2] = init

  # 通过函数概率，计算函数的选择序列
  p = np.add.accumulate(p)  
  rands = np.random.rand(n)
  select = np.ones(n, dtype=np.int)*(n-1)
  for i, x in enumerate(p[::-1]):
    select[rands<x] = len(p)-i-1

  # 结果的初始化
  result = np.zeros((n,2), dtype=np.float)
  c = np.zeros(n, dtype=np.float)

  for i in xrange(n):
    eqidx = select[i] # 所选的函数下标
    tmp = np.dot(eq[eqidx], pos) # 进行迭代
    pos[:2] = tmp # 更新迭代向量

    # 保存结果
    result[i] = tmp
    c[i] = eqidx

  return result[:,0], result[:, 1], c

start = time.clock()
x, y, c = ifs([p1,p2,p3,p4],[eq1,eq2,eq3,eq4], [0,0], 100000)
print time.clock() - start
pl.figure(figsize=(6,6))
pl.subplot(121)
pl.scatter(x, y, s=1, c="g", marker="s", linewidths=0)
pl.axis("equal")
pl.axis("off")
pl.subplot(122)
pl.scatter(x, y, s=1,c = c, marker="s", linewidths=0)
pl.axis("equal")
pl.axis("off")
pl.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)
pl.gcf().patch.set_facecolor("white")
pl.show()

程序中的 ifs 函数是进行函数迭代的主函数，我们希望通过矩阵乘法计算函数(numpy.dot) 的输出，因此需要将乘法向量扩充为三维的：

pos = np.ones(3, dtype=np.float)
pos[:2] = init

这样每次和迭代函数系数进行矩阵乘积运算的向量就变成了： x(n), y(n), 1.0 。

为了减少计算时间，我们不在迭代循环中计算随机数选择迭代方程，而是事先通过每个函数的概率，计算出函数选择数组 select，注意这里使用 accumulate 函数先将概率累加，然后产生一组 0 到 1 之间的随机数，通过判断随机数所在的概率区间选择不同的方程下标：

p = np.add.accumulate(p)
rands = np.random.rand(n)
select = np.ones(n, dtype=np.int)*(n-1)
for i, x in enumerate(p[::-1]):
  select[rands<x] = len(p)-i-1

最后我们通过调用 scatter 绘图函数将所得到的坐标进行散列图绘制：

pl.scatter(x, y, s=1, c="g", marker="s", linewidths=0)

其中每个关键字参数的含义如下：

• s : 散列点的大小，因为我们要绘制 10 万点，因此大小选择为 1
• c : 点的颜色，这里选择绿色
• marker : 点的形状，"s"表示正方形，方形的绘制是最快的
• linewidths : 点的边框宽度，0 表示没有边框

此外，关键字参数 c 还可以传入一个数组，作为每个点的颜色值，我们将计算坐标的函数下标传入，这样可以直观地看出哪个点是哪个函数迭代产生的：

pl.scatter(x, y, s=1,c = c, marker="s", linewidths=0)

下图是程序的输出：

函数迭代系统所绘制的蕨类植物的叶子

观察右图的 4 种颜色的部分可以发现概率为 1%的函数 1 所计算的是叶杆部分（深蓝色)，概率为 7%的两个函数计算的是左右两片子叶，而概率为 85%的函数计算的是整个叶子的迭代：即最下面的三种颜色的点通过此函数的迭代产生上面的所有的深红色的点。

我们可以看出整个叶子呈现出完美的自相似特性，任意取其中的一个子叶，将其旋转放大之后都和整个叶子相同。

2D 仿射变换

上面所介绍的 4 个变换方程的一般形式如下：

x(n+1) = A * x(n) + B * y(n) + C
y(n+1) = D * x(n) + E * y(n) + F

这种变换被称为 2D 仿射变换，它是从 2D 坐标到其他 2D 坐标的线性映射，保留直线性和平行性。即原来是直线上的坐标，变换之后仍然成一条直线，原来是平行的直线，变换之后仍然是平行的。这种变换我们可以看作是一系列平移、缩放、翻转和旋转变换构成的。

为了直观地显示仿射变换，我们可以使用平面上的两个三角形来表示。因为仿射变换公式中有 6 个未知数：A, B, C, D, E, F，而每两个点之间的变换决定两个方程，因此一共需要 3 组点来决定六个变换方程，正好是两个三角形，如下图所示：

两个三角形决定一个 2D 仿射变换的六个参数

从红色三角形的每个顶点变换到绿色三角形的对应顶点，正好能够决定仿射变换中的六个参数。这样我们可是使用 N+1 个三角形，决定 N 个仿射变换，其中的每一个变换的参数都是由第 0 个三角形和其它的三角形决定的。这第 0 个三角形我们称之为基础三角形，其余的三角形称之为变换三角形。

为了绘制迭代函数系统的图像，我们还需要给每个仿射变换方程指定一个迭代概率的参数。此参数也可以使用三角形直观地表达出来：迭代概率和变换三角形的面积成正比。即迭代概率为变换三角形的面积除以所有变换三角形的面积之和。

如下图所示，前面介绍的蕨类植物的分形图案的迭代方程可以由 5 个三角形决定，可以很直观地看出紫色的小三角形决定了叶子的茎；而两个蓝色的三角形决定了左右两片子叶；绿色的三角形将茎和两片子叶往上复制，形成整片叶子。

5 个三角形的仿射方程绘制蕨类植物的叶子

迭代函数系统设计器

按照上节所介绍的三角形法，我们可以设计一个迭代函数系统的设计工具，如下图所示：

具体的程序请参照迭代函数系统的分形 ，这里简单地介绍一下程序的几个核心组成部分：

首先通过两个三角形求解仿射方程的系数相当于求六元线性方程组的解，solve_eq 函数完成这一工作，它先计算出线性方程组的矩阵 a 和 b, 然后调用 NumPy 的 linalg.solve 对线性方程组 a*X = b 求解：

def solve_eq(triangle1, triangle2):
  """
 解方程，从 triangle1 变换到 triangle2 的变换系数
 triangle1,2 是二维数组：
 x0,y0
 x1,y1
 x2,y2
 """
  x0,y0 = triangle1[0]
  x1,y1 = triangle1[1]
  x2,y2 = triangle1[2]

  a = np.zeros((6,6), dtype=np.float)
  b = triangle2.reshape(-1)
  a[0, 0:3] = x0,y0,1
  a[1, 3:6] = x0,y0,1
  a[2, 0:3] = x1,y1,1
  a[3, 3:6] = x1,y1,1
  a[4, 0:3] = x2,y2,1
  a[5, 3:6] = x2,y2,1

  c = np.linalg.solve(a, b)
  c.shape = (2,3)
  return c

triangle_area 函数计算三角形的面积，它使用 NumPy 的 cross 函数计算三角形的两个边的向量的叉积：

def triangle_area(triangle):
  """
 计算三角形的面积
 """
  A = triangle[0]
  B = triangle[1]
  C = triangle[2]
  AB = A-B
  AC = A-C
  return np.abs(np.cross(AB,AC))/2.0

整个程序的界面使用 TraitsUI 库生成，将 matplotlib 的 Figure 控件通过 MPLFigureEditor 和 MPLFigureEditor 类嵌入到 TraitsUI 生成的界面中，请参考： [设计自己的 Trait 编辑器_](traitsui_manual_custom_editor.html)

IFSDesigner._figure_default 创建 Figure 对象，并且添加两个并排的子图 ax 和 ax2，ax 用于三角形编辑，而 ax2 用于分形图案显示。

def _figure_default(self):
  """
 figure 属性的缺省值，直接创建一个 Figure 对象
 """
  figure = Figure()
  self.ax = figure.add_subplot(121)
  self.ax2 = figure.add_subplot(122)
  self.ax2.set_axis_off()
  self.ax.set_axis_off()
  figure.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)
  figure.patch.set_facecolor("w")
  return figure

IFSTriangles 类完成三角形的编辑工作，其中通过如下的语句绑定 Figure 控件的 canvas 的鼠标事件

canvas = ax.figure.canvas
# 绑定 canvas 的鼠标事件
canvas.mpl_connect('button_press_event', self.button_press_callback)
canvas.mpl_connect('button_release_event', self.button_release_callback)
canvas.mpl_connect('motion_notify_event', self.motion_notify_callback)

由于 canvas 只有在真正显示 Figure 时才会创建，因此不能在创建 Figure 控件时创建 IFSTriangles 对象，而需要在界面生成之后，显示之前创建它。这里我们通过给 IFSDesigner 类的 view 属性指定其 handler 为 IFSHandler 对象，重载 Handler 的 init 方法，此方法在界面生成之后，显示之前被调用：

class IFSHandler(Handler):
  """
 在界面显示之前需要初始化的内容
 """
  def init(self, info):
    info.object.init_gui_component()
    return True

然后 IFSDesigner 类的 init_gui_component 方法完成实际和 canvas 相关的初始工作：

def init_gui_component(self):
  self.ifs_triangle = IFSTriangles(self.ax)
  self.figure.canvas.draw()
  thread.start_new_thread( self.ifs_calculate, ())
  ...

由于通过函数迭代计算分形图案比较费时，因此在另外一个线程中执行 ifs_calculate 方法进行运算，每计算 ITER_COUNT 个点，就调用 ax2.scatter 将产生的点添加进 ax2 中，由于随着 ax2 中的点数增加，界面重绘将越来越慢，因此在 draw_points 函数中限制最多只调用 ITER_TIMES 次 scatter 函数。因为在别的线程中不能更新界面，因此通过调用 wx.CallAfter 在管理 GUI 的线程中调用 draw_points 进行界面刷新。:

def ifs_calculate(self):
  """
 在别的线程中计算
 """
  def draw_points(x, y, c):
    if len(self.ax2.collections) < ITER_TIMES:
      try:
        self.ax2.scatter(x, y, s=1, c=c, marker="s", linewidths=0)
        self.ax2.set_axis_off()
        self.ax2.axis("equal")
        self.figure.canvas.draw()
      except:
        pass

  def clear_points():
    self.ax2.clear()

  while 1:
    try:
      if self.exit == True:
        break
      if self.clear == True:
        self.clear = False
        self.initpos = [0, 0]
        x, y, c = ifs( self.ifs_triangle.get_areas(),
          self.ifs_triangle.get_eqs(), self.initpos, 100)
        self.initpos = [x[-1], y[-1]]
        self.ax2.clear()

      x, y, c = ifs( self.ifs_triangle.get_areas(),
        self.ifs_triangle.get_eqs(), self.initpos, ITER_COUNT)
      if np.max(np.abs(x)) < 1000000 and np.max(np.abs(y)) < 1000000:
        self.initpos = [x[-1], y[-1]]
        wx.CallAfter( draw_points, x, y, c )
      time.sleep(0.05)
    except:
      pass

用户修改三角形之后，需要重新迭代，并绘制分形图案，三角形的改变通过 IFSTriangles.version 属性通知给 IFSDesigner，在 IFSTriangles 中，三角形改变之后，将运行：

self.version += 1

在 IFSDesigner 中监听 version 属性的变化：

@on_trait_change("ifs_triangle.version")
def on_ifs_version_changed(self):
  """
 当三角形更新时，重新绘制所有的迭代点
 """
  self.clear = True

当 IFSDesigner.clear 为 True 时，真正进行迭代运算的 ifs_calculate 方法就知道需要重新计算了。

L-System 分形

前面所绘制的分形图案都是都是使用数学函数的迭代产生，而 L-System 分形则是采用符号的递归迭代产生。首先如下定义几个有含义的符号：

• F : 向前走固定单位
• + : 正方向旋转固定单位
• - : 负方向旋转固定单位

使用这三个符号我们很容易描述下图中由 4 条线段构成的图案：

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

如此替换迭代下去，并根据字串进行绘图（符号+和-分别正负旋转 60 度)，可得到如下的分形图案：

使用 F+F--F+F 迭代的分形图案

除了 F, +, - 之外我们再定义如下几个符号：

• f : 向前走固定单位，为了定义不同的迭代公式
• [ : 将当前的位置入堆栈
• ] : 从堆栈中读取坐标，修改当前位置
• S : 初始迭代符号

所有的符号（包括上面未定义的) 都可以用来定义迭代，通过引入两个方括号符号，使得我们能够描述分岔的图案。

例如下面的符号迭代能够绘制出一棵植物：

S -> X
X -> F-[[X]+X]+F[+FX]-X
F -> FF

我们用一个字典定义所有的迭代公式和其它的一些绘图信息：

{
  "X":"F-[[X]+X]+F[+FX]-X", "F":"FF", "S":"X",
  "direct":-45,
  "angle":25,
  "iter":6,
  "title":"Plant"
}

其中：

• direct : 是绘图的初始角度，通过指定不同的值可以旋转整个图案
• angle : 定义符号+,-旋转时的角度，不同的值能产生完全不同的图案
• iter : 迭代次数

下面的程序将上述字典转换为需要绘制的线段坐标：

class L_System(object):
  def __init__(self, rule):
    info = rule['S']
    for i in range(rule['iter']):
      ninfo = []
      for c in info:
        if c in rule:
          ninfo.append(rule[c])
        else:
          ninfo.append(c)
      info = "".join(ninfo)
    self.rule = rule
    self.info = info

  def get_lines(self):
    d = self.rule['direct']
    a = self.rule['angle']
    p = (0.0, 0.0)
    l = 1.0
    lines = []
    stack = []
    for c in self.info:
      if c in "Ff":
        r = d * pi / 180
        t = p[0] + l*cos(r), p[1] + l*sin(r)
        lines.append(((p[0], p[1]), (t[0], t[1])))
        p = t
      elif c == "+":
        d += a
      elif c == "-":
        d -= a
      elif c == "[":
        stack.append((p,d))
      elif c == "]":
        p, d = stack[-1]
        del stack[-1]
    return lines

我们使用 matplotlib 的 LineCollection 绘制所有的直线：

import matplotlib.pyplot as pl
from matplotlib import collections

# rule = {...}  此处省略 rule 的定义

lines = L_System(rule).get_lines()

fig = pl.figure()
ax = fig.add_subplot(111)
linecollections = collections.LineCollection(lines)
ax.add_collection(linecollections, autolim=True)
pl.show()

下面是几种 L-System 的分形图案，绘制此图的完整程序请参照绘制 L-System 的分形图 。

几种 L-System 的迭代图案