科学可视化：Python和Matplotlib（英文原名：Scientific Visualization: Python + Matplotlib）阅读记录【第3章】

本文是阅读《科学可视化：Python和Matplotlib》（英文原名：《Scientific Visualization: Python + Matplotlib》）的记录，书和代码均为开源。

第0章 引言（Introduction）

附有链接

链接

第1章 基础原理（Fundamentals）

附有链接

链接

第2章 图像设计（Figure design）

附有链接

链接

第3章 高级概念（Advanced concepts）

3.1 动画（Animation）

使用 matplotlib 的 animation framework 能很容易地创建动画。

从一个非常简单的动画开始。想制作一个在屏幕上逐步绘制正弦和余弦函数的动画。

为此，首先需要告诉 matplotlib 我们想制作一个动画。

然后，必须指定想在每一帧绘制什么。

一个常见的错误是在每一帧重新绘制所有内容，这使得整个过程非常缓慢。相反，只需更新必要的内容，因为此时很多东西在一帧到另一帧之间不会改变。

对于 line plot，使用 set_data 来更新绘图，matplotlib 将负责其余的工作。代码如下：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig = plt.figure(figsize=(7,2), dpi=100)
ax = plt.subplot()

X = np.linspace(-np.pi, np.pi, 256, endpoint=True)
C, S = np.cos(X), np.sin(X)
line1, = ax.plot(X, C, marker="o", markevery=[-1],
                       markeredgecolor="white")
line2, = ax.plot(X, S, marker="o", markevery=[-1],
                       markeredgecolor="white")

def update(frame):
    line1.set_data(X[:frame], C[:frame])
    line2.set_data(X[:frame], S[:frame])

ani = animation.FuncAnimation(fig, update, interval=10)
plt.show()

注意随着动画移动的终点标记。原因是在末尾指定了一个标记（markevery=[-1]），这样每次设置新数据时，标记都会自动更新并随动画移动。如下图。

若要保存这个动画，matplotlib 可以创建一个视频文件，但选项相当少。

更好的解决方案是使用大多数系统上都可用的 FFPeg 等外部库。安装后，就可以使用专用的 FFMpegWriter，如下代码：

writer = animation.FFMpegWriter(fps=30)
anim = animation.FuncAnimation(fig, update, 
                               interval=10, frames=len(X))
anim.save("sine-cosine.mp4", writer=writer, dpi=100)

当保存视频时，动画不会立即开始，因为存在与创建视频相对应的延迟。

对于正弦和余弦，延迟相当短，这实际上不是问题。

然而，对于长而复杂的动画，这种延迟可能会变得相当显著，因此有必要跟踪其进度。因此，使用 tqdm 库添加一些信息。

from tqdm.autonotebook import tqdm
bar = tqdm(total=len(X))
anim.save("../data/sine-cosine.mp4", writer=writer, dpi=300,
          progress_callback = lambda i, n: bar.update(1))
bar.close()

创建时间保持不变，但至少现在可以检查它有多慢或多快。以下是动画的一些输出：

3.1.1 像降雨一样的效果（Make it rain）

一个非常简单的降雨效果可以通过在一个图像上随机放置小的生长环来获得。当然，它们不会永远生长，因为涟漪应该会随着时间的推移而衰减。

为了模拟这种现象，可以使用一种随着环的生长而变得越来越透明的颜色，直到它不再可见。此时，删除环并创建一个新的环。

首先，创建一个空白图像

fig = plt.figure(figsize=(6,6), facecolor='white', dpi=50)

ax = fig.add_axes([0,0,1,1], frameon=False, aspect=1)
ax.set_xlim(0,1), ax.set_xticks([])
ax.set_ylim(0,1), ax.set_yticks([])

然后，创建一个空散点图，注意对任何新数据都设置线宽（linewidth）为 0.5 和填充色（facecolor）为空（“None”）

scatter = ax.scatter([], [], s=[], lw=0.5,
                     edgecolors=[], facecolors="None")

接下来，需要创建几个环。

为此，可以使用散点图对象，其通常用于可视化点云，但也可以通过指定其没有填充色来绘制环。

还需注意每个环的初始尺寸和颜色，以便有最小和最大尺寸之间的所有尺寸。

此外，需要确保最大的环几乎是透明的。

n = 50
R = np.zeros(n, dtype=[ ("position", float, (2,)),
                        ("size",     float, (1,)),
                        ("color",    float, (4,)) ])
R["position"] = np.random.uniform(0, 1, (n,2))
R["size"] = np.linspace(0, 1, n).reshape(n,1)
R["color"][:,3] = np.linspace(0, 1, n)

现在，需要为动画编写更新函数。

在每个时间步长，每个环都应该生长并变得更加透明，而最大的环应该完全透明并因此被移除。

当然，实际上不会删除最大的环，而是重新使用它，在一个新的随机位置设置一个新环，具有指定的尺寸和颜色。因此，保持环的数量不变。

def rain_update(frame):
    global R, scatter
    
    # 每个环的透明度增加 (Transparency of each ring is increased)
    R["color"][:,3] = np.maximum(0, R["color"][:,3] - 1/len(R))
    
    # 每个环的尺寸增加 (Size of each rings is increased)
    R["size"] += 1/len(R)
    
    # 重置最后一个环 (Reset last ring)
    i = frame % len(R)
    R["position"][i] = np.random.uniform(0, 1, 2)
    R["size"][i] = 0
    R["color"][i,3] = 1
    
    # 相应地更新散点图对象 (Update scatter object accordingly)
    scatter.set_edgecolors(R["color"])
    scatter.set_sizes(1000*R["size"].ravel())
    scatter.set_offsets(R["position"])

最后，告诉 matplotlib 将此函数用作动画的更新函数并显示结果（或将其另存为视频）：

animation = animation.FuncAnimation(fig, rain_update,
                                    interval=10, frames=200)

3.1.2 可视化地球地震（Visualizing earthquakes on Earth）

现在使用降雨动画来可视化过去30天地球上的地震。

美国地质调查局（United States Geological Survey，USGS）地震灾害计划（Earthquake Hazards Program）是国家地震灾害减少计划（National Earthquake Hazards Reduction Program，NEHRP）的一部分，并在其网站上提供了一些数据。

这些数据根据地震震级进行排序，从仅重大地震到所有地震，无论是大地震还是小地震。

由于小地震数据太多，仅使用震级 >4.5 的地震。

首先，读取和转换数据。

使用 urllib 库，可以打开和读取远程数据。网站上的数据使用 CSV 格式，其内容由第一行给出：

time,latitude,longitude,depth,mag,magType,nst,gap,dmin,rms,...
2015-08-17T13:49:17.320Z,37.8365,-122.2321667,4.82,4.01,mw,...
2015-08-15T07:47:06.640Z,-10.9045,163.8766,6.35,6.6,mwp,...

只对纬度（latitude）、经度（longitude）和震级（magnitude）感兴趣，因此不会解析事件的时间。

import urllib
import numpy as np

# -> http://earthquake.usgs.gov/earthquakes/feed/v1.0/csv.php
feed = "http://earthquake.usgs.gov/" \
     + "earthquakes/feed/v1.0/summary/"

# 震级 > 4.5 (Magnitude > 4.5)
url = urllib.request.urlopen(feed + "4.5_month.csv")

# 震级 > 2.5 (Magnitude > 2.5)
# url = urllib.request.urlopen(feed + "2.5_month.csv")

# 震级 > 1.0 (Magnitude > 1.0)
# url = urllib.request.urlopen(feed + "1.0_month.csv")

# 读取和存储数据 (Reading and storage of data)
data = url.read().split(b'\n')[+1:-1]
E = np.zeros(len(data), dtype=[('position',  float, (2,)),
                               ('magnitude', float, (1,))])

for i in range(len(data)):
    row = data[i].split(b',')
    E['position'][i] = float(row[2]),float(row[1])
    E['magnitude'][i] = float(row[4])

需要绘制地球以精确显示地震中心的位置，并在 matplotlib 可以处理的一些坐标系中转换纬度/经度。幸运的是，有一个 cartopy 库很容易使用，但是其安装并不简单。

第一步是定义一个转换（projection），将地球绘制到屏幕上。存在许多不同的转换，但对非专业人士来说，使用等矩形转换（Equirectangular projection）是相当标准的。

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()

plt.show()

现在可以调整降雨动画来展示地震，只需在散点图中添加一个 transform，坐标就会自动变换（通过 cartopy 库）。

scatter = ax.scatter([], [], transform=ccrs.PlateCarree())

3.1.3 场景化动画（Scenarized animation）

我们已经了解了动画的基本原理。现在是时候为动画定义一个更详细的场景了。

为此，将用流体模拟（fluid simulation），因为它很有趣。

在本书代码 animation/fluid.py 中，可以看到 Gregory Johnson 基于 Joe Stam 的论文编写的稳定流体模拟的实现。

本书作者修改了原始脚本，并编写了一个 inflow 方法，在给定位置（角度）定义了一个源。在每一帧，希望定义活动源，以便整个动画显示一系列发射源。

在下面的场景中，随意地定义了一个旋转的源序列，以最大限度地提高中心的混合效果，但你也可以想象将这个动画与一些音乐同步。

import numpy as np
from fluid import Fluid, inflow
from scipy.special import erf
import matplotlib.pyplot as plt
import matplotlib.animation as animation

shape = 256, 256
duration = 500
fluid = Fluid(shape, 'dye')
inflows = [inflow(fluid, x)
           for x in np.linspace(-np.pi, np.pi, 8, endpoint=False)]

# 动画设置 (Animation setup)
fig = plt.figure(figsize=(5, 5), dpi=100)
ax = fig.add_axes([0, 0, 1, 1], frameon=False)
ax.set_xlim(0, 1); ax.set_xticks([]);
ax.set_ylim(0, 1); ax.set_yticks([]);
im = ax.imshow( np.zeros(shape), extent=[0, 1, 0, 1],
                vmin=0, vmax=1, origin="lower",
                interpolation='bicubic', cmap=plt.cm.RdYlBu)

# 动画场景 (Animation scenario)
scenario = []
for i in range(8):
    scenario.extend( [[i]]*20 )
scenario.extend([[0,2,4,6]]*30)
scenario.extend([[1,3,5,7]]*30)

# 动画更新 (Animation update)
def update(frame):
    frame = frame % len(scenario)
    for i in scenario[frame]:
        inflow_velocity, inflow_dye = inflows[i]
        fluid.velocity += inflow_velocity
        fluid.dye += inflow_dye
    divergence, curl, pressure = fluid.step()
    Z = curl
    Z = (erf(Z * 2) + 1) / 4
    
    im.set_data(Z)
    im.set_clim(vmin=Z.min(), vmax=Z.max())

anim = animation.FuncAnimation(fig, update, interval=10, frames=duration)
plt.show()

注意，在 update 函数中，限制了颜色图的更新，这是必要的，因为显示的图像是动态的，最小值和最大值可能因帧而异。如果不这样做，可能会有一些闪烁。

也可以有更详细的场景，比如一个最初由黑马分析设计的动画的重制，见下图。

3.1.4 练习（Exercise）

创建一个动画，显示如何生成李萨如曲线（Lissajous curves），下图显示了动画中的静止图像。

3.2 走向 3D（Going 3D）

Matplotlib 有一个非常好的 3D 界面，有许多功能（和一些限制），在用户中非常受欢迎。

然而，对于某些用户（或大多数用户）来说，3D 仍然被认为是某种黑魔法。

但一旦理解了一些概念，3D 渲染真的很容易。为了证明这一点，将使用60行 Python 和一个 matplotlib 调用来渲染下图中的兔子，即不使用 3D axis。

3.2.1 加载兔子（Loading the bunny）

首先，需要加载模型。

使用斯坦福兔子的简化版本。该文件使用波阵面格式（wavefront format），这是最简单的格式之一。

于是，制作一个非常简单（但容易出错）的加载器，它只会为这个特定的模型完成工作。否则，可以使用 meshio 库

V, F = [], []
with open("bunny.obj") as f:
    for line in f.readlines():
        if line.startswith('#'): continue
        values = line.split()
        if not values: continue
        if values[0] == 'v':
            V.append([float(x) for x in values[1:4]])
        elif values[0] == 'f':
            F.append([int(x) for x in values[1:4]])
V, F = np.array(V), np.array(F)-1

V 现在是一组顶点（如果你愿意，可以是 3D 点），F 是一组面（即三角形）。每个三角形由相对于顶点数组的 3 个索引描述。

对顶点进行归一化，使整个兔子适合单位框：

V = (V-(V.max(0)+V.min(0))/2)/max(V.max(0)-V.min(0))

现在，可以通过只获取顶点的 x、y 坐标并去掉 z 坐标来初步查看模型。

为此，可以使用功能强大的 PolyCollection 对象，其允许高效渲染非规则多边形的集合。因为我们想渲染一堆三角形，所以这是一个完美的匹配。因此，首先提取三角形并去掉 z 坐标：

T = V[F][...,:2]

现在可以渲染：

fig = plt.figure(figsize=(6,6))
ax = fig.add_axes([0,0,1,1], xlim=[-1,+1], ylim=[-1,+1],
                  aspect=1, frameon=False)
collection = PolyCollection(T, closed=True, linewidth=0.1,
                            facecolor="None", edgecolor="black")
ax.add_collection(collection)
plt.show()

结果如下图

3.2.2 透视投影（Perspective Projection）

刚刚进行的渲染实际上是正交投影（orthographic projection），而之前的兔子使用透视投影（perspective projection）：

在这两种情况下，定义投影的正确方法都是首先定义一个观察体，即想要在屏幕上渲染的 3d 空间中的物体。

为此，需要考虑 6 个截面（左、右、上、下、远、近），它们相对于相机包围了观察体。如果定义一个相机位置和一个观察方向，每个平面都可以用一个标量来描述。一旦有了这个观察体，就可以使用正交投影或透视投影投影到屏幕上。

幸运的是，这些投影可以用 4×4 矩阵表示：

def frustum(left, right, bottom, top, znear, zfar):
    M = np.zeros((4, 4), dtype=np.float32)
    M[0, 0] = +2.0 * znear / (right - left)
    M[1, 1] = +2.0 * znear / (top - bottom)
    M[2, 2] = -(zfar + znear) / (zfar - znear)
    M[0, 2] = (right + left) / (right - left)
    M[2, 1] = (top + bottom) / (top - bottom)
    M[2, 3] = -2.0 * znear * zfar / (zfar - znear)
    M[3, 2] = -1.0
    return M

def perspective(fovy, aspect, znear, zfar):
    h = np.tan(0.5*radians(fovy)) * znear
    w = h * aspect
    return frustum(-w, w, -h, h, znear, zfar)

对于透视投影，还需要指定孔径角度（aperture angle），该角度（多或少）设置近平面相对于远平面的大小。因此，对于高孔径，会得到很多“变形”（deformation）。

但是，上面的两个函数都返回 4×4 矩阵，而坐标是 3d。那么，如何使用这些矩阵？答案是齐次坐标（homogeneous coordinates）。齐次坐标最适合处理 3D 中的变换和投影。

在此例中，因为处理的是顶点（而不是向量），所以只需要将 1 作为所有顶点的第四个坐标（w）。然后，可以使用点积应用透视变换。

V = np.c_[V, np.ones(len(V))] @ perspective(25,1,1,100).T

最后，需要重新归一化齐次坐标。这意味着用最后一个分量（w）分割每个变换的顶点，这样每个顶点的 w=1。

V /= V[:,3].reshape(-1,1)

结果见下图，发现完全错误

其错误原因，是相机实际上在兔子体内。

为了正确渲染，需要将兔子从相机上移开，或者将相机从兔子上移开。以后再说。

相机当前位于 (0,0,0) 处，并在 z 方向上向上看。因此，在透视变换之前，需要在 z 负方向上稍微移动相机。

V = V - (0,0,3.5)
V = np.c_[V, np.ones(len(V))] @ perspective(25,1,1,100).T
V /= V[:,3].reshape(-1,1)

正确输出如下图

3.2.3 模型、视图和投影（Model, view, projection (MVP)）

不能明显看出，最后一次渲染实际上是一次透视转换。为了让其更明显，旋转兔子。

为此，需要一些旋转矩阵（4×4），同时也可以定义平移矩阵：

def translate(x, y, z):
    return np.array([[1, 0, 0, x],
                     [0, 1, 0, y],
                     [0, 0, 1, z],
                     [0, 0, 0, 1]], dtype=float)

def xrotate(theta):
    t = np.pi * theta / 180
    c, s = np.cos(t), np.sin(t)
    return np.array([[1, 0,  0, 0],
                     [0, c, -s, 0],
                     [0, s,  c, 0],
                     [0, 0,  0, 1]], dtype=float)

def yrotate(theta):
    t = np.pi * theta / 180
    c, s = np.cos(t), np.sin(t)
    return np.array([[ c, 0, s, 0],
                     [ 0, 1, 0, 0],
                     [-s, 0, c, 0],
                     [ 0, 0, 0, 1]], dtype=float)

现在，根据模型（局部变换）、视图（全局变换）和投影来分解要应用的变换，以便能计算一个全局 MVP 矩阵，该矩阵可以同时执行所有操作：

model = xrotate(20) @ yrotate(45)
view  = translate(0,0,-3.5)
proj  = perspective(25, 1, 1, 100)
MVP   = proj @ view @ model

现在可用：

V = np.c_[V, np.ones(len(V))] @ MVP.T
V /= V[:,3].reshape(-1,1)

可得如下结果

稍微调整一下孔径（aperture），可见下图的差异。还需要调整与相机的距离，以便兔子具有相同的表观尺寸。

3.2.4 深度排序（Depth sorting）

填充这些三角形，结果见下图，可以发现完全是错误的。

问题在于，PolyCollection 按照给定的顺序绘制三角形，而我们希望从后到前绘制它们。这意味着我们需要根据它们的深度对其进行排序。

幸运的是，当应用 MVP 转换时，已经计算出了这些信息，其存储在新的 z 坐标中。然而，这些 z 值是基于顶点的，而我们需要对三角形进行排序。

因此，取平均 z 值作为三角形深度的代表。如果三角形相对较小且不相交，则效果很好，所有内容都正确渲染，如下：

T =  V[:,:,:2]
Z = -V[:,:,2].mean(axis=1)
I = np.argsort(Z)
T = T[I,:]

使用存储的深度添加一些颜色。根据每个三角形的深度为其上色。PolyCollection 对象的美妙之处在于，可以使用 numpy 数组指定每个三角形的颜色：

zmin, zmax = Z.min(), Z.max()
Z = (Z-zmin)/(zmax-zmin)
C = plt.get_cmap("magma")(Z)
I = np.argsort(Z)
T, C = T[I,:], C[I,:]

结果见下图

3.2.5 练习（Exercise）

使用正交投影（orthographic projection）的定义，复现下图，其混合了透视（左上）和正交投影，每只兔子都被包含在一个 axes 中。

3.3 架构与优化（Architecture & optimization）

尽管使用 matplotlib 的架构不需要对其有深入的了解，但了解其架构以优化速度、内存甚至渲染仍然是有用的。

3.3.1 透明度级别（Transparency levels）

已经学过了如何在散点图中使用透明度来感知数据密度。

如果没有太多的数据，这表现很好，但究竟多少是太多？

这很难给出一个明确的限制，因为它取决于许多参数，如图像尺寸、标记的形状和大小以及 $\alpha$ 水平（即透明度）。

对于透明度，一种颜色的透明度实际上是有限的，正好是 0.002（=1/500）。

这意味着，如果用 0.002 的透明度绘制 500 个黑点，会得到一个准黑色的标记，如下图：

对于 500 个点来说，并不完全是黑色的，因为这还取决于 $\alpha$ 合成的内部计算方式，但这仍然提供了一个有用的近似值。

知道了这个有此限制，就解释了为什么当数据很多时，在密集区域会得到纯色，如下图。其中数据量分别为 10000、100000 和 1000000。对于 10000 和100000 个点，可以调整透明度级别，以显示密集区域的位置。此时是简单的正态分布，可以看到中心区域更暗。对于 1000000 个点，到达了透明度技巧的极限（$\alpha$=0.002），中心暗点隐藏了信息。

因此需要一种新的策略来显示数据。幸运的是，matplotlib 提供了 hist2d 和 hexbin，都会将点聚合到箱体（bin）（正方形或六边形）中，最终根据 bin 中的数量着色。这允许可视化任何数量的数据点的密度，并且不需要操纵标记的大小和/或透明度。

现在，准备好复现 Todd W.Schneider 对纽约出租车旅行的惊人可视化（分析 11 亿次纽约出租车和 Uber 旅行）。

$\alpha$ 合成会导致其他类型的折线图问题，特别是当折线图是自我覆盖的时候，如下图所示的一个高频信号。

该信号是两个不同频率的正弦波和读数的乘积：

xmin, xmax = 0*np.pi, 5*np.pi
ymin, ymax = -1.1, +1.1
def f(x): return np.sin(np.power(x,3)) * np.sin(x)
X = np.linspace(xmin,xmax,10000)
Y = f(X)

当绘制这个信号时，可以看到线条的密度从左到右变得越来越高。在图的右侧附近，频率最高，实际上高于屏幕分辨率，使得连续波之间没有空白。

然而，当使用具有一定透明度的常规绘图（上图第一行）时，看不到颜色的变化（虽然可以预料到绘图会自我覆盖）。原因是 matplotlib 渲染引擎会注意避免过度绘制属于同一绘图的区域，如下图所示：

这解释了为什么上上图中没有颜色变化。

为了抵消这种影响，可以使用由单个线段组成的线集合来渲染相同的图。

此时，每个部分都是单独考虑的，并将影响其他部分。这对应于上上图中的第二行，可以观察到颜色的变化，右侧颜色较深，表明频率较高。

还可以对信号进行多重采样（multisampling），这是信号处理中的一种标准技术。不绘制信号，而是创建了一个具有足够分辨率的空图像，对于这张图像的每个点（像素），随机考虑8个样本点，但要求这些样本点分布紧密，以确定其值。与常规绘图相比，这当然较慢，但渲染更忠实于信号，如上上图第三行所示。

3.3.2 加快渲染速度（Speed-up rendering）

渲染给定图像的总体速度取决于许多需要了解的 matplotlib 内部因素。

尽管在大多数情况下渲染速度相当不错，但当有大量对象时，情况会非常明显地恶化，而在之前的散点图示例中已经经历过这种缓慢。

注意到，渲染散点图有两种方法，可以使用仅带有标记的 plot 命令，也可以使用专用的 scatter 命令。

这两种方法既相似又不同。如果一个散点图的标记的大小、形状和颜色相同，那么可以使用更快的 plot 命令（大约是原来的两倍）。对于其他情况，应使用 scatter 命令。

可用以下代码来测量绘制一百万点的散点图的时间：

import numpy as np
import matplotlib.pyplot as plt
from timeit import default_timer as timer

n_points = 100_000
np.random.seed(1)
X = np.random.normal(0, 1, n_points)
Y = np.random.normal(0, 1, n_points)

fig = plt.figure(figsize=(9, 3.5))

# -----------------------------
ax = fig.add_subplot(
    1, 3, 2, aspect=1, xlim=[-5, 5], xticks=[], ylim=[-5, 5], yticks=[]
)
start = timer()
ax.scatter(X, Y, s=4, color="black", alpha=0.008, linewidth=0)
end = timer()
ax.set_title("Scatter: %.4fs" % (end - start))

# -----------------------------
ax = fig.add_subplot(
    1, 3, 3, aspect=1, xlim=[-5, 5], xticks=[], ylim=[-5, 5], yticks=[]
)
start = timer()
ax.plot(X, Y, "o", markersize=2, color="black", alpha=0.008, markeredgewidth=0)
end = timer()
ax.set_title("Plot: %.4fs" % (end - start))

# -----------------------------
ax = fig.add_subplot(
    1, 3, 1, aspect=1, xlim=[-5, 5], xticks=[], ylim=[-5, 5], yticks=[]
)
start = timer()
ax.hist2d(X, Y, 128, cmap="gray_r")
end = timer()
ax.set_title("Hist2d: %.4fs" % (end - start))

plt.show()

顺便一提，注意到 plot（markersize=2）和 sactter（s=2**2）之间的大小差异。原因是 plot 中标记的尺寸以点为单位进行测量，而 scatter 中标记的尺寸以平方点为单位。

对于折线图，如下图，不同方法之间的渲染速度差异可能很大。下图分别进行 1000 次调用 plot 方法（左）、一次单独调用 plot 方法（中）（其中单个线段坐标由 None 分隔）和一个线集合（右）绘制了 1000 条线段。

import numpy as np
import matplotlib.pyplot as plt
from timeit import default_timer as timer

n_lines, n_points = 1_000, 2
X = [np.random.uniform(0, 1, n_points) for i in range(n_lines)]
Y = [np.random.uniform(0, 1, n_points) for i in range(n_lines)]

fig = plt.figure(figsize=(9, 3.5))

# -----------------------------
ax = fig.add_subplot(1, 3, 1, aspect=1, xlim=[0, 1], xticks=[], ylim=[0, 1], yticks=[])
start = timer()
for x, y in zip(X, Y):
    ax.plot(x, y, color="blue", alpha=0.1, linewidth=0.5)
end = timer()
ax.set_title("Individual plots: %.4fs" % (end - start))


# -----------------------------
ax = fig.add_subplot(1, 3, 2, aspect=1, xlim=[0, 1], xticks=[], ylim=[0, 1], yticks=[])
start = timer()
X_, Y_ = [], []
for x, y in zip(X, Y):
    X_.extend(x), X_.extend([None])
    Y_.extend(y), Y_.extend([None])
ax.plot(X_, Y_, color="blue", alpha=0.1, linewidth=0.5)
end = timer()
ax.set_title("Unified plot: %.4fs" % (end - start))

# -----------------------------
from matplotlib.collections import LineCollection

ax = fig.add_subplot(1, 3, 3, aspect=1, xlim=[0, 1], xticks=[], ylim=[0, 1], yticks=[])
start = timer()
V = [np.stack([x, y], axis=1) for x, y in zip(X, Y)]
lines = LineCollection(V, color="blue", alpha=0.1, linewidth=0.5)
ax.add_collection(lines)
end = timer()
ax.set_title("Line collection: %.4fs" % (end - start))

plt.show()

此时，不同的渲染方法会有很大不同。对于大量的线条，渲染可能需要几秒钟或几分钟。注意，由于没有如前所述的自我覆盖，最快的渲染（上图中间的图）并不完全等同于其他渲染。

3.3.3 文件大小（File sizes）

根据保存图形的格式，生成的文件可能相对较小，但也可能很大，高达几兆字节（MB），

这与脚本的复杂性无关，而是与细节的数量或元素的数量有关。考虑以下脚本：

plt.scatter(np.random.rand(n=int(1e6), np.random.rand(n=int(1e6))
plt.savefig("vector.pdf")

生成的文件大小约为 15MB，原因是 pdf 格式是矢量格式，意味着需要对每个点的坐标进行编码。

上述代码有一百万个点，每个点有两个浮点坐标。一个浮点数由 4 个字节表示，那就要 8000000 个字节来存储坐标。若添加单个颜色（4字节，RGBA）和大小（1个浮点数，4字节），可以很容易地达到 16MB。

对上述代码稍微修改一下：

plt.scatter(np.random.rand(n=int(1e6), 
            np.random.rand(n=int(1e6),
            rasterized=True)
plt.savefig("vector.pdf", dpi=600)

新文件的大小约为 50KB，即使它不再是纯矢量格式，质量也大致相当。

事实上，rasterized 意味着 maplotlib 将创建散点图的光栅化（rasterized）（即位图（bitmap））表示，可以节省大量硬盘空间，还会使图像的渲染速度更快，因为 pdf 查看器不需要渲染单个元素。

然而，矢量格式与光栅化元素的组合并不总是最好的选择。

例如，如果要制作一个具有很高清晰度的巨大图像（例如海报），要是没有太多元素，纯矢量格式可能是最好的格式。

如何进行选择并没有明确的诀窍，主要取决于经验。

3.3.4 多线程渲染（Multithread rendering）

matplotlib 本身不支持多线程渲染，但无论如何都可以实现。

最明显的情况发生在需要渲染几个不同的图时，此时只需同时启动几个线程即可。

更有趣的是使用多线程渲染生成单个图像。为此，需要将图像拆分为不同且不重叠的部分，以便每个部分都可以独立渲染。

例如，考虑一个完整范围为 xlim=[0,9] 和 ylim=[0.9] 的图像。此时，很容易定义 9 个不重叠的部分：

X = np.random.normal(4.5, 2, 5_000_000)
Y = np.random.normal(4.5, 2, 5_000_000)

extents = [[x,x+3,y,y+3] for x in range(0,9,3)
                         for y in range(0,9,3)]

对于每个部分，可以用 Figure Canvas 绘制一个离线图像（offline figure），并将结果保存为图片：

def plot(extent):
    
    xmin, xmax, ymin, ymax = extent
    fig = Figure(figsize=(2,2))
    canvas = FigureCanvas(fig)
    ax = fig.add_axes([0,0,1,1], frameon=False,
                      xlim = [xmin,xmax], xticks=[],
                      ylim = [ymin,ymax], yticks=[])
    epsilon = 0.1
    I = np.argwhere((X >= (xmin-epsilon)) &
                    (X <= (xmax+epsilon)) &
                    (Y >= (ymin-epsilon)) &
                    (Y <= (ymax+epsilon)))
    ax.scatter(X[I], Y[I], 3, clip_on=False,
        color="black", edgecolor="None", alpha=.0025)
    canvas.draw()
    return np.array(canvas.renderer.buffer_rgba())

注意，要认真选择所提供范围内的 X 和 Y（模（modulo）epsilon）。这非常重要，因为不想在每个子部分中绘制所有数据，否则会减慢速度。

现在使用几个 imshow 将每个部分放在一起：

from multiprocessing import Pool

extents = [[x,x+3,y,y+3] for x in range(0,9,3)
                         for y in range(0,9,3)]
pool = Pool()
images = pool.map(plot, extents)
pool.close()

fig = plt.figure(figsize=(6,6))
ax = plt.subplot(xlim=[0,9], ylim=[0,9])
for img, extent in zip(images, extents):
    ax.imshow(img, extent=extent, interpolation="None")
plt.show()

结果如下

你可以看到不同片段的完美蒙太奇。如果将 epsilon 值设置为零，则会观察到不同部分之间出现空白。原因是，如果强制执行非常严格的剪裁，则中心在范围外的标记将不会被绘制，而由于其大小，它可能会重叠。

3.4 图像库（Graphic library）

除了用于科学可视化之外，matplotlib 还是一个综合库，像其官网上写的那样，用 Python 创建静态、动画和交互式可视化。

换句话说，matplotlib 是一个图形库，几乎可以用于任何目的，即使性能可能因渲染的复杂性和不同应用程序而异。

这种多功能性可以用许多低级对象的存在来解释，这些对象允许产生几乎任何渲染，并由下图所示的许多标准操作来支持：

上图是展示 matplotlib 在多边形（polygon）和剪裁（clipping）方面的能力的第一个示例。可以看到，剪裁允许渲染两个多边形的任意组合。

这种剪裁也可以在常规图形中使用，以产生一些有趣的效果，如下图：

3.4.1 Matplotlib地牢（Matplotlib dungeon）

如果你玩过角色扮演游戏（尤其是《龙与地下城》），可能会遇到如下图所示的“戴森阴影线（Dyson hatching）”（看墙的外边界）。

这种阴影线（hatching）样式非常独特，可以立即将该计划识别为某种地牢。这个阴影线最初是由 Dyson Logos 设计的，他很友好地解释了他是如何（手工）绘制 的。那么问题是，如何使用 matplotlib 复现它？

这实际上并不太难，但也不完全简单，因为必须处理好很多细节才能得到一个好的结果。

起点是一个随机的二维分布，其中点之间不需要太靠近。

为了实现这样的结果，可以使用 Bridson 算法，这是一种非常流行的方法来产生这种蓝噪声采样点分布，保证没有两个点比给定的距离更近。

如下图，可以看到，与纯均匀分布或具有一些正常抖动的规则网格相比，该算法确实有所不同。

从这个蓝噪声分布中，可以在每个位置以随机方向插入阴影图案。填充图案是一组 n 条平行线，其中包含一些噪声：

def hatch(n=4, theta=None):
    theta = theta or np.random.uniform(0,np.pi)
    P = np.zeros((n,2,2))
    X = np.linspace(-0.5,+0.5,n,endpoint=True)
    P[:,0,1] = -0.5 + np.random.normal(0, 0.05, n)
    P[:,1,1] = +0.5 + np.random.normal(0, 0.05, n)
    P[:,1,0] = X + np.random.normal(0, 0.025, n)
    P[:,0,0] = X + np.random.normal(0, 0.025, n)
    c, s = np.cos(theta), np.sin(theta)
    Z = np.array([[ c, s],[-s, c]])
    return P @ Z.T

结果为下图中间的子图，乍一看像戴森阴影线，但还不令人满意，因为阴影线相互覆盖。

为了避免这种情况，需要使用相应的 voronoi 单元格来裁剪阴影线。

最简单的方法是使用 shapey 库，其提供了计算与通用多边形（generic polygon）相交的方法。

结果如下图的右图，看着好多了。

下一部分是生成一个地牢。如果在互联网上搜索地牢生成器，会发现许多，从最基本的到非常复杂的。此例只用 inkscape 设计地牢，并从 svg 文件中提取了墙壁的坐标：

Walls = np.array([
    [1,1],[5,1],[5,3],[8,3],[8,2],[11,2],[11,5],[10,5],
    [10,6],[12,6],[12,8],[13,8],[13,10],[11,10],[11,12],
    [2,12],[2,10],[1,10],[1,7],[4,7],[4,10],[3,10],
    [3,11],[10,11],[10,10],[9,10],[9,8],[11,8],[11,7],
    [9,7],[9,5],[8,5],[8,4],[5,4],[5,6],[1,6], [1,1]])
walls = Polygon(Walls, closed=True, zorder=10,
                facecolor="white", edgecolor="None",
                lw=3, joinstyle="round")
ax.add_patch(walls)

下一步是将阴影线限制在墙壁附近。由于阴影线对应于初始点分布，因此只需过滤中心足够靠近任何墙壁的阴影线即可。

因此，只需计算点到线段的距离。此时，不在乎阴影线是在地牢内部还是外部，因为内部阴影线被地牢内部隐藏了。

接着在走廊里添加虚线方块，使用一组垂直和水平线以及一些随机的“岩石”，这些岩石实际上是小椭圆的集合。

最后，用一种看起来很旧的字体添加一个漂亮的标题，用的是Dieter Steffmann 的 Morris Roman 字体。

结果见下图（在本小节开头已经展示过）

结果看着不错，但还能进一步改进。例如，可以在墙上引入一些噪音来模仿手绘，可以通过添加噪音来改进岩石，等等。Matplotlib 提供了所需的一切，唯一的限制是你的想象力。

3.4.2 微型机器人模拟器（Tiny bot simulator）

使用相同的方法，可以设计一个微型机器人模拟器，如下图（模拟的截屏）。

为了设计这个模拟器，首先使用 gridspec 分割图形，如下代码：

fig = plt.figure(figsize=(10,5), frameon=False)
G = GridSpec(8, 2, width_ratios=(1,2))
ax = plt.subplot( G[:,0], aspect=1, frameon=False)
...

for i in range(8): # 8 个传感器 (8 sensors)
    sax = plt.subplot( G[i,1])
    ...

ax 是左侧显示迷宫和机器人的 axes，而 sax 是右侧显示传感器值的 axes。

迷宫墙壁使用线集合渲染，而机器人使用圆（用于身体）、线（用于“头部”，即指示方向的线）和传感器的线集合渲染。

整个模拟是一个 matplotlib 动画，其中更新函数负责更新机器人位置和传感器值。

除了传感器和墙壁交叉点的计算外，没有真正的困难，可以使用 Numpy 进行矢量化，使其加速：

def line_intersect(p1, p2, P3, P4):

    p1 = np.atleast_2d(p1)
    p2 = np.atleast_2d(p2)
    P3 = np.atleast_2d(P3)
    P4 = np.atleast_2d(P4)
    
    x1, y1 = p1[:,0], p1[:,1]
    x2, y2 = p2[:,0], p2[:,1]
    X3, Y3 = P3[:,0], P3[:,1]
    X4, Y4 = P4[:,0], P4[:,1]
    
    D = (Y4-Y3)*(x2-x1) - (X4-X3)*(y2-y1)
    
    # 共线性检验 (Colinearity test)
    C = (D != 0)
    UA = ((X4-X3)*(y1-Y3) - (Y4-Y3)*(x1-X3))
    UA = np.divide(UA, D, where=C)
    UB = ((x2-x1)*(y1-Y3) - (y2-y1)*(x1-X3))
    UB = np.divide(UB, D, where=C)
    
    # 测试交叉口是否在每个路段内
    # (Test if intersections are inside each segment)
    C = C * (UA > 0) * (UA < 1) * (UB > 0) * (UB < 1)
    X = np.where(C, x1 + UA*(x2-x1), np.inf)
    Y = np.where(C, y1 + UA*(y2-y1), np.inf)
    return np.stack([X,Y],axis=1)

这个模拟器可以很容易地一个相机进行扩展，此相机使用在第 3.2 节中介绍的渲染器，并以 3D 显示环境。

3.4.3 实例（Real example）

当这些图形元件（Graphic primitives）放在一起时，可以绘制出如下图所示的非常精细的图形。该图来自论文《一种用于放置和连接生物细胞的图形化、可扩展和直观的方法》（A graphical, scalable and intuitive method for the placement and the connection of biological cells），该论文介绍了一种源自计算机图形领域的图形方法，用于在二维流形上任意放置细胞。

上图表示基底神经节（basal ganglia）（纹状体（striatum）和苍白球（globus pallidus））的示意性切片，该切片分为四个不同的子图：

• 子图 A：由显示任意密度神经元的位图图像组成。使用位图图像，是因为还不可能使用 matplotlib 渲染如此任意的渐变。
• 子图 B：表示用于定位任意数量的神经元以强制颜色梯度表示的密度的实际方法。使用了一个简单的散点图，并根据它们的输入/输出状态给一些神经元着色。
• 子图 C：表示神经元活动的插值，并使用 2D 直方图绘制。为此，构建了一个代表整个图像的大数组，并使用一个半径恒定的圆盘来设置神经元周围的活动。这只是将神经元的二维坐标转换为图像阵列内的二维索引的问题。然后，使用 imshow 显示结果，并在每个结构的边界上绘制。这种渲染有助于查看结构内部的整体活动。
• 子图 D：可能是最复杂的，因为它涉及 Voronoi 单元的计算及其与结构边界的交点。再一次，shapely 库对于实现这样的结果非常有用。一旦计算出细胞，只需根据它们的活动为它们绘制一个颜色图即可。为了提高效率，这是使用 poly collection 制成的。

这实际上是一个相当复杂的例子，但一旦编写了代码，它就可以适应任何输入（在这种情况下是 SVG 文件），这样最终结果就完全自动化了。

当然，这代表的工作量应该与实际需求相平衡。如果只需要一次这个图像，要是能手动完成，可能就不值得付出这么多努力。

3.4.4 练习（Exercise）

1. 邮票般的效果：
   Fancy boxes 提供了几种风格，可用于实现不同的效果，如下图，复现它。

2. 放射迷宫（Radial Maze）：
   复现下图，其中显示了一个放射迷宫（在神经科学中经常用于研究小鼠或大鼠行为）和一个代表大鼠探索迷宫的模拟路径（这是通过记录（计算机）鼠标运动生成的）。每个块的颜色表示占用率，即块内记录点的数量。

第4章 案例展示（Showcase）

附有链接

链接