# 迭代函数系统的分形

```# -*- 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

# 通过函数概率，计算函数的选择序列
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.gcf().patch.set_facecolor("white")
pl.show()
```

## 迭代函数系统设计器

```# -*- coding: utf-8 -*-
from enthought.traits.ui.api import *
from enthought.traits.ui.menu import OKCancelButtons
from enthought.traits.api import *
from enthought.traits.ui.wx.editor import Editor

import matplotlib
# matplotlib采用WXAgg为后台，这样才能将绘图控件嵌入以wx为后台界面库的traitsUI窗口中
matplotlib.use("WXAgg")
from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas
from matplotlib.figure import Figure

import numpy as np
import time
import wx
import pickle

ITER_COUNT = 4000 # 一次ifs迭代的点数
ITER_TIMES = 10   # 总共调用ifs的次数

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

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

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

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

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

# 通过函数概率，计算函数的选择序列
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

class _MPLFigureEditor(Editor):
"""
使用matplotlib figure的traits编辑器
"""
scrollable = True

def init(self, parent):
self.control = self._create_canvas(parent)

def update_editor(self):
pass

def _create_canvas(self, parent):
panel = wx.Panel(parent, -1, style=wx.CLIP_CHILDREN)
sizer = wx.BoxSizer(wx.VERTICAL)
panel.SetSizer(sizer)
mpl_control = FigureCanvas(panel, -1, self.value)
sizer.Add(mpl_control, 1, wx.LEFT | wx.TOP | wx.GROW)
self.value.canvas.SetMinSize((10,10))
return panel

class MPLFigureEditor(BasicEditorFactory):
"""
相当于traits.ui中的EditorFactory，它返回真正创建控件的类
"""
klass = _MPLFigureEditor

class IFSTriangles(HasTraits):
"""
三角形编辑器
"""
version = Int(0) # 三角形更新标志

def __init__(self, ax):
super(IFSTriangles, self).__init__()
self.colors = ["r","g","b","c","m","y","k"]
self.points = np.array([(0,0),(2,0),(2,4),(0,1),(1,1),(1,3),(1,1),(2,1),(2,3)], dtype=np.float)
self.equations = self.get_eqs()
self.ax = ax
self.ax.set_ylim(-10,10)
self.ax.set_xlim(-10,10)
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)
self.canvas = canvas
self._ind = None
self.background = None
self.update_lines()

def refresh(self):
"""
重新绘制所有的三角形
"""
self.update_lines()
self.canvas.draw()
self.version += 1

def del_triangle(self):
"""
删除最后一个三角形
"""
self.points = self.points[:-3].copy()
self.refresh()

"""
添加一个三角形
"""
self.points = np.vstack((self.points, np.array([(0,0),(1,0),(0,1)],dtype=np.float)))
self.refresh()

def set_points(self, points):
"""
直接设置三角形定点
"""
self.points = points.copy()
self.refresh()

def get_eqs(self):
"""
计算所有的仿射方程
"""
eqs = []
for i in range(1,len(self.points)/3):
eqs.append( solve_eq( self.points[:3,:], self.points[i*3:i*3+3,:]) )
return eqs

def get_areas(self):
"""
通过三角形的面积计算仿射方程的迭代概率
"""
areas = []
for i in range(1, len(self.points)/3):
areas.append( triangle_area(self.points[i*3:i*3+3,:]) )
s = sum(areas)
return [x/s for x in areas]

def update_lines(self):
"""
重新绘制所有的三角形
"""
del self.ax.lines[:]
for i in xrange(0,len(self.points),3):
color = self.colors[i/3%len(self.colors)]
x0, x1, x2 = self.points[i:i+3, 0]
y0, y1, y2 = self.points[i:i+3, 1]
type = color+"%so"
if i==0:
linewidth = 3
else:
linewidth = 1
self.ax.plot([x0,x1],[y0,y1], type % "-", linewidth=linewidth)
self.ax.plot([x1,x2],[y1,y2], type % "--", linewidth=linewidth)
self.ax.plot([x0,x2],[y0,y2], type % ":", linewidth=linewidth)

self.ax.set_ylim(-10,10)
self.ax.set_xlim(-10,10)

def button_release_callback(self, event):
"""
鼠标按键松开事件
"""
self._ind = None

def button_press_callback(self, event):
"""
鼠标按键按下事件
"""
if event.inaxes!=self.ax: return
if event.button != 1: return
self._ind = self.get_ind_under_point(event.xdata, event.ydata)

def get_ind_under_point(self, mx, my):
"""
找到距离mx, my最近的顶点
"""
for i, p in enumerate(self.points):
if abs(mx-p[0]) < 0.5 and abs(my-p[1])< 0.5:
return i
return None

def motion_notify_callback(self, event):
"""
鼠标移动事件
"""
self.event = event
if self._ind is None: return
if event.inaxes != self.ax: return
if event.button != 1: return
x,y = event.xdata, event.ydata

#更新定点坐标
self.points[self._ind,:] = [x, y]

i = self._ind / 3 * 3
# 更新顶点对应的三角形线段
x0, x1, x2 = self.points[i:i+3, 0]
y0, y1, y2 = self.points[i:i+3, 1]
self.ax.lines[i].set_data([x0,x1],[y0,y1])
self.ax.lines[i+1].set_data([x1,x2],[y1,y2])
self.ax.lines[i+2].set_data([x0,x2],[y0,y2])

# 背景为空时，捕捉背景
if self.background == None:
self.ax.clear()
self.ax.set_axis_off()
self.canvas.draw()
self.background = self.canvas.copy_from_bbox(self.ax.bbox)
self.update_lines()

# 快速绘制所有三角形
self.canvas.restore_region(self.background) #恢复背景
# 绘制所有三角形
for line in self.ax.lines:
self.ax.draw_artist(line)
self.canvas.blit(self.ax.bbox)

self.version += 1

name = Str("")
view = View(
Item("name", label = u"名称"),
kind = "modal",
buttons = OKCancelButtons
)

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

class IFSDesigner(HasTraits):
figure = Instance(Figure) # 控制绘图控件的Figure对象
ifs_triangle = Instance(IFSTriangles)
del_button = Button(u"删除三角形")
save_button = Button(u"保存当前IFS")
unsave_button = Button(u"删除当前IFS")
clear = Bool(True)
exit = Bool(False)
ifs_names = List()
ifs_points = List()
current_name = Str

view = View(
VGroup(
HGroup(
Item("del_button"),
Item("current_name", editor = EnumEditor(name="object.ifs_names")),
Item("save_button"),
Item("unsave_button"),
show_labels = False
),
Item("figure", editor=MPLFigureEditor(), show_label=False, width=600),
),
resizable = True,
height = 350,
width = 600,
title = u"迭代函数系统设计器",
handler = IFSHandler()
)

def _current_name_changed(self):
self.ifs_triangle.set_points( self.ifs_points[ self.ifs_names.index(self.current_name) ] )

"""
添加三角形按钮事件处理
"""

def _del_button_fired(self):
self.ifs_triangle.del_triangle()

def _unsave_button_fired(self):
if self.current_name in self.ifs_names:
index = self.ifs_names.index(self.current_name)
del self.ifs_names[index]
del self.ifs_points[index]
self.save_data()

def _save_button_fired(self):
"""
保存按钮处理
"""
if ask.name not in self.ifs_names:
self.ifs_points.append( self.ifs_triangle.points.copy() )
else:
self.ifs_points[index] = self.ifs_triangle.points.copy()
self.save_data()

def save_data(self):
with file("IFS.data", "wb") as f:
pickle.dump(self.ifs_names[:], f) # ifs_names不是list，因此需要先转换为list
for data in self.ifs_points:
np.save(f, data) # 保存多个数组

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]
# 不绘制迭代的初始100个点
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

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

def _figure_default(self):
"""
figure属性的缺省值，直接创建一个Figure对象
"""
figure = Figure()
self.ax2.set_axis_off()
self.ax.set_axis_off()
figure.patch.set_facecolor("w")
return figure

def init_gui_component(self):
self.ifs_triangle = IFSTriangles(self.ax)
self.figure.canvas.draw()
try:
with file("ifs.data","rb") as f: