Matrix 首页推荐 

Matrix 是少数派的写作社区，我们主张分享真实的产品体验，有实用价值的经验与思考。我们会不定期挑选 Matrix 最优质的文章，展示来自用户的最真实的体验和观点。 

文章代表作者个人观点，少数派仅对标题和排版略作修改。

最近尝试用 Python 和 Matplotlib 从零手写复刻了一下 Pluribus 的片头。先看看效果：

1. 前言

最近看了 Apple TV 的一部剧叫 Pluribus。我很喜欢这部剧，原因有二：

• 它核心概念里的 "Joining" 和 《EVA》里的 「人类补完计划」 非常像，很对我的胃口；
• 剧情探讨了人类和 AI 的关系，也是我最近一直在深度思考的问题 <(")

除去剧情，我特别喜欢它的片头。极简但非常抓眼球，完全就是我的菜。Apple TV 的片头通常都很复杂且暗示剧情走向（比如《人生切割术》或者《羊毛战记》），但这一个很特别。这也是我第一次觉得「哎，这个我好像能用代码写出来」的片头 :>

2. 粒子系统 (Particle System)

因为我从来没碰过粒子系统，对计算机视觉也知之甚少，所以上手第一步就是先读几篇文章。下面这两个资源对拆解概念非常有帮助：

• [文章] Particle Systems From the Ground Up
• [视频] How Particle Life emerges from simplicity

简单来说，我只需要一堆点，然后追踪它们的物理状态：位置、速度和加速度。

class Particle:
    def __init__(self, pos: (int, int), 
                 velocities: (int, int), 
                 accelerations: (int, int)):
        self.pos = pos
        self.vel = velocities
        self.acc = accelerations

套用高中物理学过的标准公式：

写个函数来更新这些值：

def pos_update(dot, dt):
    dot.pos = (
        dot.pos[0] + dot.vel[0] * dt,
        dot.pos[1] + dot.vel[1] * dt
        )
    dot.vel = (dot.vel[0] + dot.acc[0] * dt,
                dot.vel[1] + dot.acc[1] * dt)

对每个点跑这个循环，就能得到一个基础的粒子系统（渲染代码略过不表，不过这里有个很好的 matplotlib 动画教程）。

最后，给每个点加点随机力。假设质量(m)为 1，根据 F=ma，我们可以直接把随机值加到加速度上：

def force_apply(p: Particle):
    p.acc = (
        p.acc[0] + random.randint(-2, 2), 
        p.acc[1] + random.randint(-2, 2)
        )

def dots_update(dots, dt):
    for dot in dots:
        pos_update(dot, dt)
        force_apply(dot)
    return

初始化网格里的点之后，大概长这样：

3. 背景点 (Background-dots)

把片头看了五遍以后，我发现里面的点可以分为三类，各个击破：

• 背景点 (Background-dots)
• 圆圈点 (Circle-dots)
• 文字点 (Text-dots)

对于背景点，简单的随机运动看着不自然。如果你仔细看（现在是第六遍了 :D），会发现它们之间是有交互的。基本上就是太近了会推开，太远了会拉近。我发现 Lennard-Jones 势能完美描述了这个行为：

简单说就是距离太近会排斥，距离远了（但在范围内）会吸引。就像下图这个曲线。（我是从这个博客学来的）。

实现起来也很简单，遍历每一对点应用这个力就行，复杂度是 O(n^2)。

def lj_force(p1, p2):
    dx = p1.pos[0] - p2.pos[0]
    dy = p1.pos[1] - p2.pos[1]
    dis = (dx**2 + dy**2) ** 0.5

    dx_dir = dx / dis
    dy_dir = dy / dis

    u = min(10, 4 * EPI * ((SIGMA/dis)**12 - (SIGMA/dis)**6))

    dx_acc = u * dx_dir / 1
    dy_acc = u * dy_dir / 1

    p1.acc = (p1.acc[0]+dx_acc, p1.acc[1]+dy_acc)
    p2.acc = (p2.acc[0]-dx_acc, p2.acc[1]-dy_acc)

加上 LJ 势能后的效果如下。能明显看到点之间相互作用产生的复杂运动。

4. 圆圈点 (Circle-dots)

加圆圈点之前，先快速复习一下如何在粒子系统中定义方向和距离。（记得的同学可以跳过 :O）

基本上给定一个角度 θ∈[0,2π) 我们可以得到方向的单位向量 ​dir_x=cos(θ) dir_y=sin(θ)​。给定两个点，我们可以得到从 p1​ 到 p2​ 的方向：

要得到方向（单位向量），我们用差值除以距离：

加圆圈点很容易。给个初始速度，按 2π（360度）均匀分布方向就行。

def add_wave(dots):
    for i in range(WAVE_DOTS_NUM):
        angle = 2 * math.pi * i / WAVE_DOTS_NUM
        
        pos = (WAVE_ORIGIN[0] + math.cos(angle)*5, 
            WAVE_ORIGIN[1] + math.sin(angle)*5)
        
        vx = WAVE_SPEED * math.cos(angle)
        vy = WAVE_SPEED * math.sin(angle)
            
        dots.append(Particle(pos, velocities=(vx, vy)))

• 碰撞问题： 但这里有个坑。因为我们加了 LJ 力，背景点会和圆圈点互怼。圆圈扩大的时候，撞到背景点会被推歪，形状就散了。
• 解决方案： 我的解法简单粗暴：给 Particle 类加个 mass（质量）属性。让圆圈点比背景点重得多，它们惯性就大，不容易被推跑。

更新物理计算遵循牛顿第二定律 (a=F/m)。基本就是更新速度的时候，把累计的力（加速度）除以质量：

def pos_update(dot, dt):
    dot.pos = (
        dot.pos[0] + dot.vel[0] * dt,
        dot.pos[1] + dot.vel[1] * dt
        )
    dot.vel = (
        dot.vel[0] + dot.acc[0] * dt / dot.mass,
        dot.vel[1] + dot.acc[1] * dt / dot.mass
        )

对比一下（左：无质量，右：有质量）。

加了质量以后看着舒服多了吧？能明显看到圆圈点把背景点推开，自己还能保持队形。

5. 文字点 (Text-dots)

用点渲染文字不难。找个字体（我用的 Arial）画出来，然后提取像素位置就行。

def get_text_draw(text = TEXT, font_path = FONT_PATH):
    mask_img = Image.new("L", (WIDTH, LENGTH), 0)
    draw = ImageDraw.Draw(mask_img)
    font = ImageFont.truetype(font_path, 35)

    bbox = draw.textbbox((0, 0), text, font=font)
    text_w, text_h = bbox[2] - bbox[0], bbox[3] - bbox[1]
    draw.text(((WIDTH - text_w) // 2, (LENGTH - text_h) // 2 - 5), text, fill=255, font=font)
    y_coords, x_coords = np.where(np.array(mask_img)[::-1] > 128)
    return x_coords, y_coords

难点在于做那个「指纹」图案。仔细看原片，它像个波浪，稍微有点不规则。为了简单，我用 sine wave 模拟：

基本上就是根据距离中心的远近推拉这些点。调整频率能搞出不同的环形图案。下图是 freq={1,4,7} 的效果。

def set_fingerprint(x, y, freq = RADIAL_FREQ, strengh = RADIAL_STRENGTH):
    dx = x_coords - RADIAL_ORIGIN[0]
    dy = y_coords - RADIAL_ORIGIN[1]

    dist = np.sqrt(dx**2 + dy**2)
    angle = np.arctan2(dy, dx)

    push = np.sin(dist * freq) * strengh

    x_new = x_coords + (np.cos(angle) * push)
    y_new = y_coords + (np.sin(angle) * push)
    return x_new, y_new

如下是从点 P(25,42) 发起正弦波应用到文字的效果。

其实调这个波的参数花了我好久。试了各种组合，最后选了个看着最舒服的。^_^

把所有东西合在一起，就有了第一版片头！8)

6. 性能优化 (Performance Optimization)

先停一下。目前渲染60帧要跑6分钟。感觉我在浪费生命等它跑完 :( 是时候做点优化了。

6.1 空间哈希 (Spatial Hashing)

前面说了，瓶颈在物理交互计算，复杂度 O(n2)。加上文字点和不断生成的圆圈点，数量轻松上千，意味着每帧要做 10^6 次距离检测。

我的解法是用空间哈希（分桶），把空间划成网格，只计算相邻网格里粒子的 LJ 力。灵感来自第 3 节的公式：距离 ≥3σ 时势能几乎归零。

我用哈希表记录每个点属于哪个格子：

def _bin_coords(pos):
    return int(pos[0]) // BIN_SIZE, int(pos[1]) // BIN_SIZE

def _build_bins(dots):
    bins = {}
    for idx, p in enumerate(dots):
        bx, by = _bin_coords(p.pos)
        if 0 <= bx < BIN_XNUM and 0 <= by < BIN_YNUM:
            bins.setdefault((bx, by), []).append(idx)
    return bins

这一改，速度提升了 5 倍，渲染时间从 6 分 10 秒降到了 1 分 06 秒。

（虽然我知道用树结构——类似二叉索引树——动态维护位置能把复杂度降到 O(nlogn)，毕竟最近在刷 LeetCode。但网格法目前够用了。）

6.2 生命周期管理

另一个优化是控制点的生命周期。圆圈点飞出屏幕（「越界」）后就不用算了。我加了个定期清理。这对减少内存占用很有效，之前内存都飙到 10GB 了。

def prune_dots(dots, circles, margin=50):
    alive_dots = []
    alive_circles = []

    for dot, circle in zip(dots, circles):
        x, y = dot.pos
        if -margin < x < WIDTH + margin and -margin < y < LENGTH + margin:
            alive_dots.append(dot)
            alive_circles.append(circle)
        else:
            circle.remove()

    dots[:] = alive_dots
    circles[:] = alive_circles

我很确定用内存池（链表+哈希表）能做到 O(1) 的插入删除，但对于这个项目有点杀鸡用牛刀了 :/

7. 视觉打磨 (Visual Optimization)

接下来打磨一下视觉效果。

7.1 文字形状

第一个问题是文字时间长了会「糊」掉或者散架。因为点挤得太紧，LJ 势能把它们推开了，导致我们（搞了半天的）指纹纹理丢了。

解决办法很简单：加个锚点力 (Anchor Force)。就像个弹簧，点飘太远了就把它拽回原位。我还加了点阻尼（摩擦力）防止它震荡个没完。

def anchor_force(p):
    dx = p.anchor[0] - p.pos[0]
    dy = p.anchor[1] - p.pos[1]
    dis = (dx**2 + dy**2) ** 0.5
    dx_dir = dx / dis
    dy_dir = dy / dis

    f = dis * ANCHOR_STRENGH

    damping_fx = -p.vel[0] * DAMPING
    damping_fy = -p.vel[1] * DAMPING

    p.acc = (
        p.acc[0] + (f * dx_dir + damping_fx) * random.randrange(5, 10) / 10, 
        p.acc[1] + (f * dy_dir + damping_fy) * random.randrange(5, 10) / 10
        )

7.2 呼吸与循环

另一个改进是给背景点加个「呼吸」效果，大小有节奏地缩放。给每个粒子加个相位属性，用正弦波更新就行。

最后，为了防止背景点飞出屏幕，我做了个屏幕循环 (Screen wrapping)。点从右边出去，就从左边回来。

def pos_update(dot, dt):
    dot.pos = (
        dot.pos[0] + dot.vel[0] * dt,
        dot.pos[1] + dot.vel[1] * dt
    )
    dot.vel = (
        dot.vel[0] + dot.acc[0] * dt / dot.mass,
        dot.vel[1] + dot.acc[1] * dt / dot.mass
    )
    dot.acc = (0, 0)

    dot.phase = (dot.phase + PHASE_INCREMENT) % (2 * math.pi)
    sine_wave = (math.sin(dot.phase) + 1) / 2

    if dot.type == 0:
        ## Keep background dots
        dot.vel = (dot.vel[0] * DECAY_RATIO, dot.vel[1] * DECAY_RATIO)
        dot.pos = (dot.pos[0] % WIDTH, dot.pos[1] % LENGTH)
        ## Change their size periodically
        dot.radius = 0.5 * (0.4 + 0.6 * sine_wave)
    else:
        dot.radius = 0.5 * (0.9 + 0.1 * sine_wave)

效果图解：

当然你也可以把文字换成任何你想要的：

8. 总结

这其实是我第一次尝试写粒子系统。本来计划在剧终（圣诞节）前搞定，但我高估了自己旅行时的精力和专注度。说实话，理解原理并实现它确实花了我不少时间。

相比之下，我看很多人用 Gemini 生成那种酷炫的 web 端 3D 粒子系统。跟那些比，我这个可能显得简陋甚至有点「丑」。但对我来说，从零构建的这个过程要更 enjoyable，虽然这肯定不是最高效的方法。最后，我觉得这种感觉大概也就是《Pluribus》想表达的东西吧。 :V