用 python matplotlib 画分形图
sicp 中 picture language 的 matplotlib 实现
[2023-09-25 Mon] 原题 "sicp 中 picture language 的 matplotlib 实现" 改为 subtile
picture language 来自 sicp 的 2.2.4 节, 个人发现网络上能找到的 python 实现大多都是基于 python turtle 包, 这不适合在 jupyter notebook (有 jupyter turle 包，但我没有尝试过)或者 emacs org babel 里运行，但我更喜欢在这类文学编程的环境下试验代码， 因此尝试用 matplotlib 作为绘图的底层实现。
注：本文并不是对 picture-langage 的介绍，更多是备注性质，其中很多用词不规范（自造的），理解原理还请阅读原文。
import numpy as np
import matplotlib.pyplot as plt
plt.plot 绘图的分解动作
python - How to draw a line with matplotlib? - Stack Overflow
plt.plot 基本用法就是 plot(x, y), 分别对应所有点的 x 坐标和 y 坐标 例如，下图在点(1,3) 和 (2,4) 之间连线
plt.rcParams["figure.figsize"] = (3,3)
plt.plot([1,2], [3,4]);
默认情况下，相邻的点会按传入顺序用线段自动连接（一笔划），因此，如果要 plot 出一个形状，需要考虑这个图是否能一笔划出，还要设计好点的传入顺序，以下是一个传入顺序混乱的例子
构造一个小人，其坐标是（用 lisp 的写法）：
(let ((head (make-vect 0.5 1))
(larm (make-vect 0 0.7)) ;left arm
(rarm (make-vect 1 0.3)) ;right arm
(body (make-vect 0.5 0.5))
(ltoe (make-vect 0.3 0)) ;left toe
(rtoe (make-vect 0.7 0))) ;right toe
把 x 和 y 坐标分开传入 plt.plot 的结果是：
x = [0.5, 0.0, 1.0, 0.5, 0.3, 0.7]
y = [1.0, 0.7, 0.3, 0.5, 0.0, 0.0]
plt.plot(x, y);
这些点的坐标是正确的，然而连接的顺序有问题，我们希望点之间是按以下代码里声明的关系相互连接的，但上图则是 按照 head-larm,larm-rarm,rarm-body,body-ltoe,ltoe-rtoe 的连接方式
(list (make-segment head body)
(make-segment larm body)
(make-segment rarm body)
(make-segment ltoe body)
(make-segment rtoe body))
plot 可以传入多组 x 和 y 坐标列表，这种语法能够比较方便地处理不同的点之间的特定连接（不是一笔划，从而避免传入过多重复的点）：
plt.plot([1, 2], [3, 4], [0, 1, 1.4], [1, 2, 3]);
点 (1,3) 与 (2,4) 之间画上了线段（右上），而点 (0, 1), (1, 2), (1.4, 3) 依次构成另外一组连接（左下）
这等价于：
plt.plot([1,2],[3,4])
plt.plot([0, 1, 1.4], [1, 2, 3])
plt.show()
以下提供对点（vect），线的构造函数，选择函数以及需要的基本运算
于是，小人的组成如下：
head = (0.5, 1)
larm = (0, 0.7)
rarm = (1, 0.3)
body = (0.5, 0.5)
ltoe = (0.3, 0)
rtoe = (0.7, 0)
stick_man = [(head, body), (larm, body), (rarm, body), (ltoe, body),(rtoe, body)]
为了把 stick_man 形式转成 plot 能够识别的方式，以下函数将其中每一个线段转为向量形式，转置后拼接成一个 2nx2 的 array ，然后解包成 2n 个参数输入到 plot 中。以下 arrow 参数会突出线段的第二个点，用于模拟向量的展示。
def segments2plot(segments, c=None, arrow=False):
pairs = np.concatenate(list(map(lambda x: np.array(x).T, segments)))
plt.plot(*pairs, c=c)
if arrow:
ends = np.array([s[1] for s in segments]).T
plt.plot(*ends, "bo")
注意由于之后要经常调用这个函数，所有的点都是要画在同一个画布上，因此 fig 和 ax 是全局变量
segments2plot(stick_man)
给连线加上一个端点：
segments2plot(stick_man, None, True)
来自 SICP 2.2.4 Example: A Picture LanguageをPythonでやってみたらデキタ━━━━(ﾟ∀ﾟ)━━━━ｯ!! - 牌語備忘録 -pygo 的例子（后文很多函数也来自此文）
meo = [
[[0.0, 0.0], [0.2, 1.0]],
[[0.2, 1.0], [0.3, 0.8]],
[[0.3, 0.8], [0.7, 0.8]],
[[0.7, 0.8], [0.8, 1.0]],
[[0.8, 1.0], [1.0, 0.0]],
[[0.4, 0.6], [0.3, 0.6]],
[[0.3, 0.6], [0.3, 0.5]],
[[0.3, 0.5], [0.4, 0.5]],
[[0.8, 0.6], [0.7, 0.6]],
[[0.7, 0.6], [0.7, 0.5]],
[[0.7, 0.5], [0.8, 0.5]],
[[0.7, 0.4], [0.4, 0.4]],
[[0.4, 0.4], [0.4, 0.3]]
]
segments2plot(meo, c='r')
vect(point), segment,draw-line 封装
有了以上的尝试，本节开始按照从低层到高层的方向构造 picture language（SICP 更多是从上往下，但也不是绝对），这样可以测试每一次抽象层的构建，所有的接口的命名和书本保持一致
以下是对点和线段的构造、选择和绘图函数
def make_vect(x, y):
return (x, y)
def xcor_vect(v):
return v[0]
def ycor_vect(v):
return v[1]
def add_vect(v1, v2):
return make_vect(xcor_vect(v1) + xcor_vect(v2),
ycor_vect(v1) + ycor_vect(v2))
def sub_vect(v1, v2):
return make_vect(xcor_vect(v1) - xcor_vect(v2),
ycor_vect(v1) - ycor_vect(v2))
def scale_vect(a, v):
return make_vect(a * xcor_vect(v), a * ycor_vect(v))
def draw_line(p1, p2):
segments2plot([[p1,p2]])
draw_line(make_vect(0,0), make_vect(1,1))
segment 的封装，
def make_segment(p1, p2):
return (p1, p2)
def start_segment(s):
return s[0]
def end_segment(s):
return s[1]
Frame 的构造
frame 的核心是构造坐标系，以及按照该坐标系对点进行转换，它接收的是三个向量：原点以及两个主轴。
frame_coord_map 函数返回一个闭包（坐标变换函数），该闭包接收一个向量，对这个向量进行坐标变换后返回新的向量。
def make_frame(origin, edge1, edge2):
return (origin, edge1, edge2)
def origin_frame(frame):
return frame[0]
def edge1_frame(frame):
return frame[1]
def edge2_frame(frame):
return frame[2]
def frame_coord_map(frame):
return lambda v:add_vect(
origin_frame(frame),
add_vect(scale_vect(xcor_vect(v),
edge1_frame(frame)),
scale_vect(ycor_vect(v),
edge2_frame(frame))))
可视化 frame, 底层还是调用 segements2plot 来展现
def plot_frame(frame):
o = origin_frame(frame)
v1 = add_vect(o, edge1_frame(frame))
v2 = add_vect(o, edge2_frame(frame))
segments2plot((((0, 0), o), (o, v1), (o, v2)), None, True)
standard_frame = make_frame(make_vect(0, 0),
make_vect(1, 0),
make_vect(0, 1))
print_frame = standard_frame
dia_frame = make_frame(make_vect(0.5,0.5), make_vect(1,1), make_vect(0,1))
plot_frame(dia_frame)
plot_frame(standard_frame)
传入 frame_coord_map 的向量是在源坐标系里的一个点，例如以下把标准坐标里的 (1,1) 转到了 dia_frame 坐标系里的 (1.5, 2.5)
frame_coord_map(dia_frame)(make_vect(1,1))
| 1.5 | 2.5 |
Painter 构造
painter 是一个闭包，它包含这一系列原始坐标（标准坐标系的坐标），当要绘制 painter 时，接收一个新的坐标系，将原始坐标转换到新的坐标系中得到最终坐标，然后将其 plot.
由于只从 segment 进行构造，因此其构造函数名称就称为 segment_list：
def null_list(lst):
return None
def segments2painter(segment_list):
return lambda frame: null_list([draw_line(
frame_coord_map(frame)(start_segment(segment)),
frame_coord_map(frame)(end_segment(segment)))
for segment in segment_list])
meo_painter = segments2painter(meo)
meo_painter(standard_frame)
meo_painter(dia_frame)
transform painter 的各类基础操作抽象层
transform_painter 函数接收一个 painter 以及新的坐标系的三个点的绝对坐标，不同于 frame 组成中的向量，用绝对坐标可以更直观地理解坐标变换，但由于是绝对坐标，需要转为相对于原点的向量之后后输入到 make_frame 中
def transform_painter(painter, origin, corner1, corner2):
return lambda frame:painter(make_frame(
frame_coord_map(frame)(origin),
sub_vect(frame_coord_map(frame)(corner1),
frame_coord_map(frame)(origin)),
sub_vect(frame_coord_map(frame)(corner2),
frame_coord_map(frame)(origin))))
可以看到，该函数返回一个 painter 闭包，闭包的输入参数 frame 直接对 origin, corner1, cornern2 进行坐标变化，将变换的结果用于构建一个新的 frame ，并把该 frame 喂入 painter 里。
一般来说， frame 始终是标准坐标系，即我们不用 painter 的参数 frame 来转换坐标，而是使用 transform_painter 函数，因此 lambda 内的 sub_vect 的结果就是新的 edge 基向量。如果对同一个 painter 加两层 transform_painter, 转换的 frame 分别是 f1 和 f2, 那么最终相当于得到 lambda f: painter(f2*f1), 这里 f2*f1 理解函数的复合，即先对 painter 里所有点做 f1 变换然后做 f2 变换，这个过程纯粹是函数的组合，只有在最后执行 painter(frame) 的时候才会真正对点进行变换并且打印出来，类似惰性求值的感觉。
def flip_vert(painter):
return transform_painter(painter,
make_vect(0.0, 1.0),
make_vect(1.0, 1.0),
make_vect(0.0, 0.0))
def shrink_to_upper_right(painter):
return transform_painter(painter,
make_vect(0.5, 0.5),
make_vect(1.0, 0.5),
make_vect(0.5, 1.0))
def rotate90(painter):
return transform_painter(painter,
make_vect(1.0, 0.0),
make_vect(1.0, 1.0),
make_vect(0.0, 0.0))
def squash_inwards(painter):
return transform_painter(painter,
make_vect(0.0, 0.0),
make_vect(0.65, 0.35),
make_vect(0.35, 0.65))
def beside(painter1, painter2):
split_point = make_vect(0.5, 0.0)
paint_left = transform_painter(painter1,
make_vect(0.0, 0.0),
split_point,
make_vect(0.0, 1.0))
paint_right = transform_painter(painter2,
split_point,
make_vect(1.0, 0.0),
make_vect(0.5, 1.0))
return lambda frame:null_list([paint_left(frame), paint_right(frame)])
def flip_horiz(painter):
return transform_painter(painter,
make_vect(1.0, 0.0),
make_vect(0.0, 0.0),
make_vect(1.0, 1.0))
def rotate180(painter):
return rotate90(rotate90(painter))
def rotate270(painter):
return transform_painter(painter,
make_vect(0.0, 1.0),
make_vect(0.0, 0.0),
make_vect(1.0, 1.0))
def below(painter1, painter2):
return rotate90(beside(rotate270(painter1), rotate270(painter2)))
below(meo_painter, meo_painter)(print_frame)
rotate180(meo_painter)(print_frame)
更高层的操作
def up_split(painter, n):
if n == 0:
return painter
smaller = up_split(painter, n - 1)
return below(painter, beside(smaller, smaller))
def right_split(painter, n):
if n == 0:
return painter
smaller = right_split(painter, n - 1)
return beside(painter, below(smaller, smaller))
def corner_split(painter, n):
if n == 0:
return painter
up = up_split(painter, n - 1)
right = right_split(painter, n -1)
top_left = beside(up, up)
bottom_right = below(right, right)
corner = corner_split(painter, n - 1)
return beside(below(painter, top_left),
below(bottom_right, corner))
def square_limit(painter, n):
quarter = corner_split(painter, n)
half = beside(flip_horiz(quarter), quarter)
return below(flip_vert(half), half)
plt.rcParams["figure.figsize"] = (5,5)
square_limit(meo_painter, 3)(print_frame)
stick_man_painter = segments2painter(stick_man)
square_limit(stick_man_painter, 4)(print_frame)
还可以用这个语言来绘制其他的分形图像，比如：
def triple(painter):
up0 = make_vect(0.25, 0.5)
up1 = make_vect(0.75, 0.5)
up2 = make_vect(0.25, 1)
left0 = make_vect(0, 0)
left1 = make_vect(0.5, 0)
left2 = make_vect(0, 0.5)
right0 = left1
right1 = make_vect(1, 0)
right2 = make_vect(0.5, 0.5)
paint_up = transform_painter(painter, up0, up1, up2)
paint_left = transform_painter(painter, left0, left1, left2)
paint_right = transform_painter(painter, right0, right1, right2)
return lambda frame:null_list([paint_up(frame), paint_left(frame), paint_right(frame)])
def triple_split(painter, n):
if n == 0:
return painter
new = triple_split(painter, n - 1)
return triple(new)
segments2plot([[[0,0], [1,0]], [[1,0],[0.5,np.sqrt(3)/2]], [[0,0],[0.5,np.sqrt(3)/2]]])
tri_painter = segments2painter([[[0,0], [1,0]], [[1,0],[0.5,1]], [[0,0],[0.5,1]]])
triple_split(tri_painter, 6)(standard_frame)
