python 几何库 计算曲面的交线

本章学习利用matplotlib画几何图形和分形（fractal）

使用MATPLOTLIB库PATCHES包画几何图形

回顾之前的画图：

plt.plot(...)
plt.show()

也可以拆解为多步骤：

>>> import matplotlib.pyplot as plt
>>> plt.figure()	# 建立一个空图表
<Figure size 640x480 with 0 Axes>
>>> plt.gcf()		# 返回当前的图表，会自动创建新图表如果之前没有创建过
<Figure size 640x480 with 0 Axes>
>>> plt.axes()		# 建立坐标轴
<matplotlib.axes._subplots.AxesSubplot object at 0x7ffbdad50310>
>>> plt.gca()		# 返回当前坐标轴，会自动创建带坐标轴的新图表如果之前没有创建过
<matplotlib.axes._subplots.AxesSubplot object at 0x7ffbdad50310>
>>> plt.show()

画圆

>>> circle = plt.Circle((0, 0), radius = 0.5)
>>> ax = plt.gca()
>>> ax.add_patch(circle)
<matplotlib.patches.Circle object at 0x7ffbc17f09d0>
>>> plt.axis('scaled')	# 设置自动调整坐标轴
(-0.55, 0.55, -0.55, 0.55)
>>> plt.show()	# 显示坐标轴和一个实心圆

创建动画

以下代码显示一个循环往复由小逐渐变大的圆。

import matplotlib.pyplot as plt
from matplotlib import animation

#这个函数书固定写法，第一个参数是表示帧序号，是自动传入的
def update_radius(i, circle):
       circle.radius = i*0.5	# 半径是帧序号 * 0.5, 因此半径会逐渐变大
       return circle
	   
fig = plt.gcf()
ax = plt.axes(xlim=(-10, 10), ylim=(-10, 10))		# 设定x, y坐标轴的范围从-10到10
ax.set_aspect('equal')		# 设定纵横比为1:1
circle = plt.Circle((0, 0), radius = 0.05)				# 初始圆，半径非常小，看上去是一小点
ax.add_patch(circle)
# update_radius是更新画面函数，其中帧序号自动传入，因此这里只传了circle
# frames是帧数量
# interval是两帧之间的间隔时间，单位为毫秒。将其改为1000（1秒）可看到慢动作
anim = animation.FuncAnimation( fig, update_radius, fargs = (circle,), frames=30, interval=50)
plt.title('Simple Circle Animation')
plt.show()

需要说明的是，fargs = (circle,)不能写成fargs = (circle)，因为这里需要传入一个tuple，而不是一个circle对象：

>>> type((1))
<class 'int'>
>>> type((1,))
<class 'tuple'>

另外，我们将animation.FuncAnimation赋给变量anim，但这个变量后续并没有引用。这是一个bug，目前必须这么做。
所有30帧播放完后，又会重新从第1帧开始播放，可以指定repeat=False不做循环。

制造弹道轨迹动画

这个轨迹是通过画点实现的，也就是在每一帧更新的是点的位置circle.center = x, y。
这里最难的其实还是计算x轴的宽度和y轴的高度。以及在某一时刻点的位置。
动画不太对，以后再更正。

画分形

分形是由简单公式生成的几何图形，生活中常见的分形有海岸线，树和雪花。

平面上点的变换

变换是指通过规则形成新的点的位置。规则可以是多种，每次应用的规则可以从中随机选择。

$P_{new}=Rule(P_{old})$

例如两个规则$Rule_1 = (x+1, y+1)$，$Rule_2 = (x+1, y-1)$，x总是向右移动，y则一个向上，一个向下。通过随机选择就可以形成锯齿(zigzag)形。

在作者代码中，有一点可以学习。就是将函数作为参数和返回值。

见以下代码片段，作者用类似的方法来选择规则：

import random

def f1():
	print('I am f1')

def f2():
	print('I am f2')

def call_f(f):
	f()

flist = [f1, f2]
for i in range(5):
	f = random.choice(flist)
	call_f(f)

输出如下：

$ p3 randomfunc.py
I am f1
I am f1
I am f2
I am f1
I am f2

这里函数没有参数，当然也可以加参数。

除list外, random.choice()也适用于tuple和字符串：

>>> import random
>>> random.choice('HelloWorld!')
'r'
>>> random.choice((2,3,4))
3

绘制Barnsley Fern分形

巴恩斯利蕨型(Barnsley Fern)是英国数学家Michael Barnsley发明的分形。

这里有4个转换规则：

规则1：

$X_{n+1}= 0.85X_n + 0.04Y_n$

$Y_{n+1}= -0.04X_n + 0.85Y_n+1.6$

规则2：

$X_{n+1}= 0.2X_n - 0.26Y_n$

$Y_{n+1}= 0.23X_n + 0.22Y_n+1.6$

规则3：

$X_{n+1}= -0.15X_n + 0.28Y_n$

$Y_{n+1}= 0.26X_n + 0.24Y_n+0.44$

规则4：

$X_{n+1}= 0$

$Y_{n+1}= -0.16Y_n$

同时，这4中转换规则是不均匀分布的，概率分别为0.85, 0.07, 0.07, 0.01

不均匀分布在上一章讲过了。

以下就是Barnsley Fern效果：

我的代码如下，这回是一次通过的：

import random
import matplotlib.pyplot as plt

def tr1(x, y):
	return (0.85*x + 0.04*y), (-0.04*x + 0.85*y + 1.6)


def tr2(x, y):
	return (0.2*x - 0.26*y), (0.23*x + 0.22*y + 1.6)


def tr3(x, y):
	return (-0.15*x + 0.28*y), (0.26*x + 0.24*y + 0.44)


def tr4(x, y):
	return (0), (0.16*y)


def calc_sum_props(props):
	sum_props=[]
	for i in range(len(props)):
		if i == 0:
			sum_props.append(props[0])
		else:
			sum_props.append(props[i] + sum_props[i-1])

	return sum_props


def next_point(x, y):
	r = random.random()
	for i in range(len(sum_props)):
		if r <= sum_props[i]:
			break;

	return rules[i](x, y)


rules = [tr1, tr2, tr3, tr4]
props = [.85, .07, .07, .01]
sum_props = calc_sum_props(props)
x = [0,]
y = [0,]

maxp = int(input('please input number of points:'))

for n in range(maxp):
	x_next, y_next = next_point(x[n], y[n])
	x.append(x_next)
	y.append(y_next)

plt.plot(x, y, 'go')
plt.title('Barnsley Fern Fractal')
plt.show()

编程挑战

在正方形中填充圆

学到两点。
第一，绘制多边形（此处是正方形）：

square = plt.Polygon([(1, 1), (5, 1), (5, 5), (1, 5)], closed=True)

第二，可多次调用add_patch然后一起显示：

ax = plt.gca()
ax.add_patch(square)
for i in range(len(circle)):
	ax.add_patch(circle[i])
plt.show()

绘制 Sierpiński 三角分形

3个转换规则，概率均匀分布：
规则1：
$X_{n+1}= 0.5X_n$
$Y_{n+1}= 0.5Y_n$
规则2：
$X_{n+1}= 0.5X_n + 0.5$
$Y_{n+1}= 0.5Y_n + 0.5$
规则3：
$X_{n+1}= 0.5X_n + 1$
$Y_{n+1}= 0.5Y_n$

这个不难，仿造Barnsley Fern分形即可。

输出如下：

代码如下：

import random
import matplotlib.pyplot as plt

def tr1(x, y):
	return (0.5*x), (0.5*y)


def tr2(x, y):
	return (0.5*x + 0.5), (0.5*y + 0.5)


def tr3(x, y):
	return (0.5*x + 1), (0.5*y)


def calc_sum_props(props):
	sum_props=[]
	for i in range(len(props)):
		if i == 0:
			sum_props.append(props[0])
		else:
			sum_props.append(props[i] + sum_props[i-1])

	return sum_props


def next_point(x, y):
	r = random.random()
	for i in range(len(sum_props)):
		if r <= sum_props[i]:
			break;

	return rules[i](x, y)


rules = [tr1, tr2, tr3]
props = [1/3, 1/3, 1/3]
sum_props = calc_sum_props(props)
x = [0,]
y = [0,]

maxp = int(input('please input number of points:'))

for n in range(maxp):
	x_next, y_next = next_point(x[n], y[n])
	x.append(x_next)
	y.append(y_next)

plt.plot(x, y, 'go')
plt.title('Sierpiński Trangle Fractal')
plt.show()

探索Hénon’s Function函数

扩展阅读：Fractals: A Very Short Introduction by Kenneth Falconer (Oxford University Press, 2013)

转换规则：
$x_{n+1}=y+1-1.4x^2$
$y=0.3x$

我的代码：

import matplotlib.pyplot as plt

x = [1,]
y = [1,]

maxp = int(input('please input number of points:'))

for i in range(maxp):
	x_next = y[i] + 1 -  1.4 * x[i]**2
	y_next = 0.3 * x[i]
	x.append(x_next)
	y.append(y_next)

plt.title('Henon’s Function Fractal')
plt.plot(x, y, 'go')
plt.show()

挑战升级：制作轨迹动画。作者做了个视频。

绘制Mandelbrot Set

首先介绍了matplotlib’s imshow() 函数。参考这里。

及这里。

用此函数画出了灰度色块：

灰度是由随机数产生的，取值为0到20。

代码如下：

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import random

x_p = 20
y_p = 20

image = []
for i in range(y_p):
	image.append([0] * x_p)

for i in range(y_p):
	for j in range(x_p):
		image[i][j] = random.randint(0, 10)

# origin='lower'表示image[0][0]对应坐标(0,0)的颜色
# extent=(0, 8, 0, 8)表示区域为0,0到8,8
# interpolation='nearest'表示按最近（之前）的颜色插补，因为在image中只指定了
# 64个点的灰度，其它的点并无对应的值，需要按指定规则插补，这也是为何会形成方块的原因
plt.imshow(image, origin='lower', extent=(0, 8, 0, 8), cmap=cm.Greys_r, interpolation='nearest')
plt.colorbar()	# 右侧显示色带
plt.show()

在初始化image时，不能写成[[0] * x_p] *y_p，看下面的示例：

>>> a = [0] * 4
>>> b = [a] * 4
>>> b
[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
>>> b[0][0] = 5
>>> b
[[5, 0, 0, 0], [5, 0, 0, 0], [5, 0, 0, 0], [5, 0, 0, 0]]

很奇怪，我赋一个值，结果4个地方全改了，估计这4个都是指针。

利用imshow，就可以实现Mandelbrot Set的可视化了。应用灰度，数值越大，颜色越深。

效果如图，运行时间在几十秒左右：

类似的算法称为逃逸时间算法（Escape Time Algorithm)。图中白色的部分属于Mandelbrot集，黑色的部分就是escape的部分。

我的代码如下，参照了作者的代码。很多可以学习的地方。

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import random

x_p = 400
y_p = 400
x0, x1 = -2.5, 1
y0, y1 = -1, 1
max_iteration = 1000


def init_image():
	image = []

	for i in range(y_p):
		image.append([0] * x_p)
	return image

def calc(x, y):
	iteration = 0
	z1 = complex(0 ,0) 
	c = complex(x, y)

	while iteration < max_iteration:
		z1 = z1**2 + c
		iteration += 1
		if (abs(z1) < 2):
			continue
		else:
			break
	
	return iteration


image = init_image()
dx = (x1-x0)/(x_p-1)
dy = (y1-y0)/(y_p-1)
x_coords = [x0 + i*dx for i in range(x_p)]
y_coords = [y0 + i*dy for i in range(y_p)]

for i,x in enumerate(x_coords):
	for k,y in enumerate(y_coords):
		image[k][i] = calc(x, y)

plt.imshow(image, origin='lower', extent=(x0, x1, y0, y1), cmap=cm.Greys_r, interpolation='nearest')
plt.colorbar()
plt.show()

值得学习的几点：

1. 常数都用变量表示，且放在文件最前。
2. 简洁的赋值，例如x0, y0 = 1, 1用于坐标的赋值是非常好的；其实也可以(x,y) = (1,1)
3. 简洁的序列赋值，例如x = [ 1 + i for i in range(5)]，x为[1, 2, 3, 4, 5]
4. 以上是小节，最重要的是理解了坐标系。绘图和计算灰度用的绝对坐标系，而image是相对坐标系，坐标都是整数

<未完，还有许多扩展阅读>