Python数学编程 第六章 绘制几何图形与分形 第二节 绘制分形

6.2 绘制分形

分形势复杂的几何图案或形状，但却由简单的数学公式所生成。与圆形和矩形等几何图形相比，分形看上去并不规则，也没有明显的模式或描述。但如果仔细观察，你将发现模式的存在，即整个图形有无数个自身的副本构成。由于分形设计平面上的点的相同几何变换的重复引用，因此非常适合使用计算机程序来构建。这一章中，我们将学习如何绘制Barnsley蕨类植物、Sierpinski三角和Manderlbrot集合（后两个出现在编程挑战中），这些都是分形研究领域的常见例子。分形在自然界中比比皆是，常见的例子包括海岸线、树木和雪花。

6.2.1 平面上点的变化

) l6 c& f9 |& `! T) b7 |1 I

创建分形的基本思想是点的变换。给定x-y平面上的一个点P(x,y)，变幻的一个例子是P(x,y)->Q(x+1, y+1)，这意味着应用变换后，创建好了一个新点Q，它位于P向上和向右各一个单位处。如果以Q点作为起始点，将得到另一个点R，它位于点Q向上和向右各一个单位处。假设起始点P位于（1，1），下图展示了这些点的位置。

P、Q、R三个点的位置图

因此，变换是一种描述点在x-y平面上移动的规则，从起始位置开始，每一次迭代使得点移动到新的位置。我们可以把变换想象成点在平面上的轨迹。现在换个思路考虑两个变换规则，每次迭代将从这两个规则中随机选择一个。两个规则如下：

5 K7 p# a9 F4 Z

规则1：P1(x,y) -> P2(x+1, y-1)

x! S( l ^* P8 h

规则2：P1(x,y) -> P2(x+1, y+1)

6 y5 ]) f4 ?; G$ {' U, F

以P(1, 1)作为起始点，如果我们执行4次迭代，将得到以下一系列的点：

5 ^& t+ D$ H2 J: @6 N& c

P1(1, 1) -> P2(2, 0) (规则1)

P2(2, 0) -> P3(3, 1) (规则2)

$ ~9 c; U' B/ \' ?

P3(3, 1) -> P4(4, 2) (规则2)

P4(4, 2) -> P5(5, 1) (规则2)

...

( L( \3 H8 i) b4 W, `

两个变换规则是随机选择的，其中每个规则被选择的概率相同。无论选择哪一个规则，点都会右移，因为在两种规则中都会增加x坐标。当点向右移动时，他们会向上或向下移动，从而形成一条Z字路径。以下程序绘制了在指定的迭代次数内进行上述转换规则之一时点的路径：

8 V+ C% L1 c, b6 j% G Example of selecting a transformation from two equally probable transformations % n. T7 @2 z6 ]6 z @ import matplotlib.pyplot as plt import random 6 g" n0 p% }5 R 9 w3 ^# v5 U9 \0 i5 r" z( e def transformation_1(p): . Z# I. X' l9 ]/ |9 o x = p[0] ! O! M' ^, F8 u# \) o0 Z# Y4 M y = p[1] return x + 1, y - 1 def transformation_2(p): x = p[0] * d; Y8 a3 @# p# c y = p[1] + E& V4 V% b* d$ G k" o7 e return x + 1, y + 1 def transform(p): # List of transformation functions transformations = [transformation_1, transformation_2] # pick a random transformation function and call it t = random.choice(transformations) + T; | B3 A6 [ x, y = t(p) return x, y def build_trajectory(p, n): x = [p[0]] & Q. ~4 r2 h+ {1 J% l y = [p[1]] - @' q. }6 ]. ^2 s( v4 o: u* E for i in range(n): " S6 L9 L9 F4 e$ y T1 \ p = transform(p) * y0 o% _; O* Z) c4 |( v0 p x.append(p[0]) y.append(p[1]) return x, y . T9 m4 {6 p0 x& S: c+ _! i + b. v$ b/ L- m; S. c4 \& u' V if __name__ == __main__: # Initial point p = (1,1) ) e3 U; R, d, r3 [; b" q/ P' g& b6 Q9 ^ n = int(input(Enter the number of iterations: )) * j* ?* E; k4 Z x , y = build_trajectory(p, n) # PLot % C$ a& d( _! S0 y+ @. g9 X/ a; B Y plt.plot(x, y) & o- Q+ L$ ?% E# t0 o; P plt.xlabel(X) plt.ylabel(Y) plt.show()

我们定义了两个函数transformation_1()和transformation_2()，分别对应于之前的两种变换规则，在transform()函数中，我们创建了一个包含这两个函数名字的列表，然后使用random.choice()函数从该列表中选择去哦中一个进行变换。现在假设我们选择了其中一个变换，蒋点P作为输入参数调用该变换函数，然后分别使用x，y来存储变换后的坐标并且返回它们。

9 G' X7 q" V% K1 |. f& |

从列表中随机选择一个元素

我们在第一个分形程序中看到的random.choice()函数可以用来从列表中随机选择一个元素，每一个元素被选中的概率相同，下面是一个例子：

import random 5 G) j. ~# ?" \ >>> l = [1,2,3] >>> random.choice(l) 1 >>> random.choice(l) 3 >>> random.choice(l) 3 >>> random.choice(l) 1 >>> random.choice(l) 1 8 Q3 c' D( E$ I: F0 K/ C2 z

这个函数也适用于元组和字符串。在后一种情形中，函数从字符串中随即返回一个字符。

当运行程序，他首先会询问迭代次数n，即应用变化的次数。然后调用build_trajectory()函数，该函数的输入参数为迭代次数n和起始点坐标p，此处设置初始点坐标为（1，1）。build_trajectory()函数重复调用n次transform()函数，并使用两个列表x和y分别存储变换后的x和y坐标。最后使用这两个列表来绘制图形。

1 [) `/ R$ ~; w6 |8 @ F

下面两张图分别显示了100次和10000次迭代的点的轨迹图，这两幅图中的Z字形运动都很明显，这种Z字形路径通常被称为一条线上的随机游走。

迭代100次的Z字路径

) D1 ]! C/ t4 ^, q' e' d

- Q3 I3 [4 u2 I5 `* w* |

$ y3 L) a1 K2 D0 {

迭代10000次的Z字路径

) ?* }0 V4 _/ T7 ]

, H8 V" ~' z8 x; i+ y% t

这个例子展示了创建分形的基本思想，即从一个点开始并对其重复应用某种变换。接下来，我们将看到一个用相同想法绘制Barnsley蕨类植物的例子。

3 ^7 }% }$ N4 o

6.2.2 绘制Barnsley蕨类植物

英国数学家Michael Barnsley描述了如何对一个点进行重复的简单变换，从而创建蕨类植物的结构（如下图）。

# m# m' p% H3 V

蕨类植物

' N% _+ V7 U- u0 t {$ T" }

# A2 [1 E6 W" R( I/ e" c

他提出了以下步骤来创建类似于蕨类植物的结构，以（0，0）为初始点，按事先分配的概率随机选择下述某种变换。

& o6 h3 \& u1 M$ C- h# x# j

（1）变换1（概率为0.85）：

! b0 I3 O$ U. O" {" {0 R1 ~" i6 S/ ~

（2）变换2（概率为0.07）：

% |' }" z( T4 A( \3 b

（3）变换3（概率为0.07）：

% u8 b# J1 [" G# @/ ]2 v

（4）变换4（概率为0.01）：

上述每一个变换对应于蕨类植物的一部分。第一个变换被选中的概率最大，因此被执行的次数最多，从而产生了蕨类植物的茎和底部的叶子。第二个和第三个变换分别对应于左边和右边底部的叶子，第四个变换绘制了蕨类植物的茎。

这是一个非均匀概率选择的例子，我们在第五章中已经学习过，以下程序为指定的点数绘制Barnsley蕨类植物。

Draw a Barnsley Fern import random import matplotlib.pyplot as plt 1 b& X% z6 g$ { def transformation_1(p): x = p[0] y = p[1] x1 = 0.85*x + 0.04*y b! T9 I8 U. g- C# z! K, L y1 = -0.04*y + 0.85*y + 1.6 9 ]: i6 b9 a5 a' i: i4 A4 P" \ return x1, y1 def transformation_2(p): x = p[0] y = p[1] 1 Z1 O7 z% ]3 s' Y1 k* e, ? x1 = 0.2*x - 0.26*y 0 B! w+ q! n( _ y1 = 0.23*y + 0.22*y + 1.6 6 P! h; p( }6 I2 n; Y x( Z9 @/ @ return x1, y1 # \0 X( ?% {/ D6 T& X, T4 Z def transformation_3(p): x = p[0] ! {3 H: q" D+ J$ c1 J y = p[1] 7 ]; h$ N. `) V; J# s w x1 = -0.15*x + 0.28*y " y% f: `# H- K7 E y1 = 0.26*y + 0.24*y + 0.44 9 z" m6 Y4 z+ \. Q( z return x1, y1 U% T. n+ n& L$ U) h def transformation_4(p): x = p[0] y = p[1] x1 = 0 & e* l! C0 ^4 `- G9 W- ]1 n y1 = 0.16*y return x1, y1 def get_index(probability): r = random.random() ( X0 l$ N% y5 ~- i& N" X7 u c_probability = 0 sum_probability = [] : d3 E. d7 y4 C6 i2 {' B2 i7 n+ s for p in probability: & z9 L6 W8 q+ a& {# x* @, n c_probability += p sum_probability.append(c_probability) " _! x2 R% [/ f# v for item, sp in enumerate(sum_probability): # l7 \0 g) s& a if r < sp: return item return len(probability) - 1 + U6 M# Z: d2 C$ l, s def transform(p): 3 k' F# d7 \3 U+ k: M! j # List of transformation functions 4 H$ u6 N0 z7 }0 E" C# w6 H1 E transformation = [transformation_1, transformation_2, transformation_3, transformation_4] # C2 b# m+ u3 l, X2 |" i8 j* z2 _: h probability = [0.85, 0.07, 0.07, 0.01] # Pick a random transformation function and call it - n1 o, w0 l, ~ tindex = get_index(probability) 8 _: S7 \3 t6 r6 P/ p4 R$ L4 k( v t = transformation[tindex] 5 t' r* N0 b/ x' q; w1 p" q* q x, y = t(p) return x, y ! I6 W6 J3 V3 @3 v3 w; | def draw_fern(n): + p( [( p# v" ~8 D8 l* E3 w1 W- d # We start with (0,0) x = [0] 1 l1 s& S/ H x y = [0] 8 G! W0 {- G. Y, U; s8 r1 O. _7 d x1, y1 = 0, 0 for i in range(n): 2 b& \2 _0 o6 e/ K. L7 l x1, y1 = transform((x1, y1)) x.append(x1) y.append(y1) 3 K$ [/ E# _- B& `/ y6 B' p return x, y - Y- i }$ R5 W0 M" x2 V* U if __name__ == __main__: ! r% p7 Q. ]2 t5 I n = int(input(Enter the number of points in the Ferns: )) x, y = draw_fern(n) H5 d! f$ b7 C7 o$ T plt.plot(x, y, o) & ^, z8 A( n8 n/ u& e( @ plt.title(Fern with {0} points.format(n)) plt.show() * \( P! g8 N ^# u* I

运行此程序，他首先询问蕨类植物图中指定的点数，然后开始绘制该图。下面两张图分别展示了具有1000个点和10000个点的蕨类植物图。

& H4 V# i; X* C9 y4 ]4 l

制定了1000个点的蕨类植物叶片

f' [# ?1 M9 q5 Q1 S- _8 u6 n

; w9 @6 k7 j6 g4 Y: f$ ~+ L

. }! }4 |& q' ~: _5 |. r3 I; u9 R

指定了10000个点的蕨类植物叶片

6 ?* K6 @+ i8 y& L

这4个变换规则分别在transformation_1()，transformation_2()，transformation_3()，transformation_4()函数中定义。get_index()函数返回按照给定概率下的一种随机变化的索引值，这个概率我们之前已经讲过了，不清楚的可以返回去在第五章非均匀概率那里再次学习。返回索引之后即可选择某一个具体的变换来应用。

& y5 b5 {: s3 b9 o5 w0 _2 G

初始点（0，0）变换的次数即为程序输入中指定的蕨类植物点的个数。

) r$ }6 E U1 U( V. Y% A; D: o* F ( e1 v5 k: U [; c+ s 5 _$ k) ~6 Q6 v. n6 b+ q( A; ` + P7 b: J! Q Z+ l( F8 Z& n