Python 绘制分形图(曼德勃罗集、分形树叶、科赫曲线、分形龙、谢尔宾斯基三角等)附代码

lypr7463 8年前
   <p style="text-align:center"><img src="https://simg.open-open.com/show/f948af9904df92b647cf12887064b80f.jpg"></p>    <h3>1. 曼德勃罗集</h3>    <pre>  <code class="language-python">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.jet, 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()</code></pre>    <p style="text-align:center"><img src="https://simg.open-open.com/show/00d8f6099cea9215ef1da4e88a6eddf0.png"></p>    <h3>2. 分形树叶</h3>    <pre>  <code class="language-python">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 range(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)  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("#D3D3D3")  pl.show()</code></pre>    <p style="text-align:center"><img src="https://simg.open-open.com/show/69f23767aeb42c7016652bd1a1aa79c2.png"></p>    <h3>3. 其它分形图(科赫曲线、分形龙、谢尔宾斯基三角等)</h3>    <pre>  <code class="language-python">from math import sin, cos, pi  import matplotlib.pyplot as pl  from matplotlib import collections    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    rules = [  {  "F":"F+F--F+F", "S":"F",  "direct":180,  "angle":60,  "iter":5,  "title":"Koch"  },  {  "X":"X+YF+", "Y":"-FX-Y", "S":"FX",  "direct":0,  "angle":90,  "iter":13,  "title":"Dragon"  },  {  "f":"F-f-F", "F":"f+F+f", "S":"f",  "direct":0,  "angle":60,  "iter":7,  "title":"Triangle"  },  {  "X":"F-[[X]+X]+F[+FX]-X", "F":"FF", "S":"X",  "direct":-45,  "angle":25,  "iter":6,  "title":"Plant"  },  {  "S":"X", "X":"-YF+XFX+", "Y":"+XF-YFY-FX+",  "direct":0,  "angle":90,  "iter":6,  "title":"Hilbert"  },  {  "S":"L--F--L--F", "L":"+R-F-R+", "R":"-L+F+",  "direct":0,  "angle":45,  "iter":10,  "title":"Sierpinski"  },    ]    def draw(ax, rule, iter=None):  if iter!=None:  rule["iter"] = iter  lines = L_System(rule).get_lines()  linecollections = collections.LineCollection(lines)  ax.add_collection(linecollections, autolim=True)  ax.axis("equal")  ax.set_axis_off()  ax.set_xlim(ax.dataLim.xmin, ax.dataLim.xmax)  ax.invert_yaxis()    fig = pl.figure(figsize=(7,4.5))  fig.patch.set_facecolor("papayawhip")    for i in xrange(6):  ax = fig.add_subplot(231+i)  draw(ax, rules[i])    fig.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)  pl.show()</code></pre>    <p style="text-align:center"><img src="https://simg.open-open.com/show/713231ec1f45ae415511556aa469fc06.png"></p>    <p> </p>    <p> </p>    <p>来自:https://zhuanlan.zhihu.com/p/25792397</p>    <p> </p>