Back to Blog
MathPython

终端甜甜圈的数学

SP

Shinwoo PARK

4 min readTranslated
Info

本文是从我的 Velog 搬运过来的文章
原文可在这里查看。
这是修正了深色背景、LaTeX,以及内容错误等多个问题的版本
如果还有需要修改的部分,请联系 shinwoo.park@psw.kr

原文撰写日期:2022/12/4
更新日期:2025/7/3

这种东西,你至少见过一次吧?

donut.c-spinning-donut.gif

在黑色的屏幕上,用白色特殊字符做出一个旋转的 3D 甜甜圈。
看起来明明像是个简单的程序,却莫名很耐看。
这就是 Donut.c。
像这样用 ASCII 字符串绘图的方式称为 ASCII Art。

“哎呀,说到底不就是在控制台窗口里输出文字嘛。”

但如果这么轻视它可不行。
既然和图形有关,一旦深入挖掘,就会展开一个深奥的数学世界。

这篇文章想尝试理解其中的数学原理。

Donut.c

我对这个控制台甜甜圈的兴趣,始于一则 YouTube 视频。
那部视频是 "why you NEED math for programming",翻译过来就是“为什么编程需要数学”。

视频在结尾用简要的数学原理解释了这个甜甜圈是如何被渲染成 ASCII 码的。
听完之后我产生了兴趣,于是上网搜索了一下。

令人惊讶的是,关于这个的文章比想象中少得多。
不过只有一些简略的韩文介绍文章,以及讲解复杂英文数学原理的文章而已。

一开始因为英语和复杂数学混在一起,让人有点不想读,但下定决心读下去后,发现其实还是能读的。

那么现在就真正开始分析吧。
这里涉及的数学原理参考了这篇文章。

关于 ASCII Art

首先来看看,一般的 ASCII Art 是如何表达明暗的。

我们输出 ASCII Art 的地方是终端,而终端通常是深色背景、浅色文字。
因此,在一个字符所占据的空间里,如果像素密度高,就可以视为亮色;像素密度低,就可以视为暗色。
下面这些字符,用简单的 12 个字符表示了这种区别。
.,-~:;=!*#$@
左边在固定空间内密度最低,所以最暗;右边在固定空间内密度最高,所以最亮。

数学原理

先来想想,要如何实现这个东西。
需要两件事。

  1. 画出甜甜圈的形状
  2. 为甜甜圈的形状表现明暗

好,现在就从 3D 第一人称视角中使用的简单数学开始吧。

上图表示一个人坐在屏幕前,观看屏幕后方物体的情景。 donut.c-eye-screen-object.png 为了把三维物体画到二维平面上,需要将三维中的每个点 (x,y,z)(x, y, z)(x,y,z) 投影到距离观察者 z′z'z′ 的平面上。
那么这个点就会变成 (x′,y′)(x', y')(x′,y′)。
因为我们已知三维物体每个点的坐标,所以需要知道它应当投影到二维的哪个位置,也就是要通过 xxx、yyy 求出 x′x'x′、y′y'y′。 (由于 z′z'z′ 是常数,我们记作 K1K_1K1​。)

用数学式表示就是 y′K1=yz\dfrac{y'}{K_1} = \dfrac{y}{z}K1​y′​=zy​
移项得到 y′=yK1zy' = \dfrac{yK_1}{z}y′=zyK1​​
这个式子同样适用于 xxx:x′=xK1zx' = \dfrac{xK_1}{z}x′=zxK1​​
因此投影方程就是 (x′,y′)=(K1xz,K1yz)(x', y')=(\dfrac{K_1x}{z}, \dfrac{K_1y}{z})(x′,y′)=(zK1​x​,zK1​y​)。
K1K_1K1​ 可以根据想要显示在 2D 窗口中的视角任意设定。例如,如果有一个 100x100 大小的窗口,那么视野会对准 (50, 50)。此时,如果你想看到一个在 3D 空间中相对于视野距离为 5、宽度为 10 的物体,那么必须让点 x=10x=10x=10、z=5z=5z=5 的投影出现在屏幕上,因此:

  • x′<50x'<50x′<50
  • 即 10K1z<50\dfrac{10K_1}{z}<50z10K1​​<50
  • 即 K1<25K_1 < 25K1​<25

另外,在投影很多点时,也可能会遇到 xxx 和 yyy 相同但 zzz 不同的点,因此我们要为绘制的所有内容保存一个 z-buffer。有了 z-buffer,就可以在某个位置显示一个点时,判断它是否比该位置原本已有的内容更靠前。
此外,z-buffer 还能帮助我们计算并缓冲 z−1=1zz^{-1}=\dfrac{1}{z}z−1=z1​,原因如下。

  1. 当 z−1=0z^{-1} = 0z−1=0 时,说明 zzz 接近无穷大,因此可以把 z-buffer 的默认值设为 0,并把背景设置为无限深度。
  2. 在求 x′x'x′ 和 y′y'y′ 时还可以再次利用。 (因为 z−1=1zz^{-1}=\dfrac{1}{z}z−1=z1​,提前把 zzz 除掉,就不必再对 zzz 做两次除法,运算会更高效一些。)

甜甜圈(solid of revolution)

好,那么现在要考虑如何构造甜甜圈的形状。
幸运的是,甜甜圈是一个 旋转体(solid of revolution),因此可以通过下面这样的原理求出它的每一个点。
donut.c-eye-screen-object.png 如图所示,以点 (R2,0,0)(R_2, 0, 0)(R2​,0,0) 为圆心,有一个半径为 R1R_1R1​ 的圆。
围绕 R1R_1R1​ 的这个圆,通过从 0360 度旋转就可以画出来。
把从 0
360 度变化的变量记作 θ\thetaθ。

把它写成式子如下。
(x,y,z)=(R2,0,0)+(R1cos⁡θ,R1sin⁡θ,0)(x, y, z) = (R_2, 0, 0) + (R_1\cos\theta, R_1\sin\theta, 0)(x,y,z)=(R2​,0,0)+(R1​cosθ,R1​sinθ,0)
这个式子可以求出一个圆。

现在再让这个圆沿着 2y 轴旋转。你应该已经注意到,这样就会得到甜甜圈的形状。
要让这个圆绕 y 轴旋转,可以使用 旋转变换矩阵 来求得。

y 轴旋转的三维旋转变换矩阵如下,乘起来就会得到下面的式子。

(R2+R1cos⁡θ,R1sin⁡θ,0)⋅(cos⁡ϕ0sin⁡ϕ010−sin⁡ϕ0cos⁡ϕ)=((R2+R1cos⁡θ)cos⁡ϕ,R1sin⁡θ,−(R2+R1cos⁡θ)sin⁡ϕ)\begin{align} & (R_2+R_1\cos\theta,R_1\sin\theta,0)\cdot\left(\begin{matrix}\cos\phi & 0 & \sin\phi \\ 0 & 1 & 0 \\ -\sin\phi & 0 & \cos\phi\end{matrix}\right)\\= & ((R_2+R_1\cos\theta)\cos\phi, R_1\sin\theta, -(R_2+R_1\cos\theta)\sin\phi)\end{align}=​(R2​+R1​cosθ,R1​sinθ,0)⋅​cosϕ0−sinϕ​010​sinϕ0cosϕ​​((R2​+R1​cosθ)cosϕ,R1​sinθ,−(R2​+R1​cosθ)sinϕ)​​

不过,我们希望实现这个甜甜圈 绕两个以上的轴旋转,因此再取两个旋转变换矩阵相乘。
另外,这两个新增旋转变换矩阵中使用的角度分别记作 A 和 B。
式子如下。

(R2+R1cos⁡θ,R1sin⁡θ,0)⋅(cos⁡ϕ0sin⁡ϕ010−sin⁡ϕ0cos⁡ϕ)⋅(1000cos⁡Asin⁡A0−sin⁡Acos⁡A)⋅(cos⁡Bsin⁡B0−sin⁡Bcos⁡B0001)(R_2+R_1\cos\theta,R_1\sin\theta,0)\cdot\left(\begin{matrix}\cos\phi & 0 & \sin\phi \\ 0 & 1 & 0 \\ -\sin\phi & 0 & \cos\phi\end{matrix}\right)\cdot\left(\begin{matrix}1&0&0\\0&\cos A&\sin A\\0&-\sin A&\cos A\end{matrix}\right)\cdot\left(\begin{matrix}\cos B&\sin B&0\\-\sin B&\cos B & 0\\0 & 0 & 1\end{matrix}\right)(R2​+R1​cosθ,R1​sinθ,0)⋅​cosϕ0−sinϕ​010​sinϕ0cosϕ​​⋅​100​0cosA−sinA​0sinAcosA​​⋅​cosB−sinB0​sinBcosB0​001​​

现在,我们已经得到了一个以原点 (0,0,0)(0, 0, 0)(0,0,0) 为基准、由角度 AAA 和 BBB 持续旋转的甜甜圈坐标。

接下来,要把它真正投影到屏幕上,需要假设观察者坐标为 (0,0,0)(0, 0, 0)(0,0,0),并把甜甜圈移动到观察者前方。
为了把甜甜圈移到观察者前方,需要给 zzz 加上 观察者与甜甜圈之间的距离。
把这个 观察者与甜甜圈之间的距离 记作 K2K_2K2​。

这样一来,用于求实际 (x′,y′)(x', y')(x′,y′) 的式子如下。

(x′,y′)=(K1xK2+z,K1yK2+z)(x', y') = \left(\dfrac{K_1x}{K2 + z}, \dfrac{K_1y}{K_2+z}\right)(x′,y′)=(K2+zK1​x​,K2​+zK1​y​)

而用于求 (x,y,z)(x, y, z)(x,y,z) 的式子如下。

(xyz)=((R2+R1cos⁡θ)(cos⁡Bcos⁡ϕ+sin⁡Asin⁡Bsin⁡ϕ)−R1cos⁡Asin⁡Bsin⁡θ(R2+R1cos⁡θ)(cos⁡ϕsin⁡B−cos⁡Bsin⁡Asin⁡ϕ)+R1cos⁡Acos⁡Bsin⁡θcos⁡A(R2+R1cos⁡θ)sin⁡ϕ+R1sin⁡Asin⁡θ)\left(\begin{matrix}x\\y\\z\end{matrix}\right)=\left(\begin{matrix}(R_2+R_1\cos\theta)(\cos B\cos\phi+\sin A\sin B\sin\phi)-R_1\cos A\sin B\sin\theta\\(R_2+R_1\cos\theta)(\cos\phi\sin B-\cos B\sin A\sin\phi)+R_1\cos A\cos B\sin\theta\\\cos A(R_2+R_1\cos\theta)\sin\phi+R_1\sin A\sin\theta\end{matrix}\right)​xyz​​=​(R2​+R1​cosθ)(cosBcosϕ+sinAsinBsinϕ)−R1​cosAsinBsinθ(R2​+R1​cosθ)(cosϕsinB−cosBsinAsinϕ)+R1​cosAcosBsinθcosA(R2​+R1​cosθ)sinϕ+R1​sinAsinθ​​

为了尽可能优化,在编程阶段不会使用矩阵乘法,而是直接使用上面的计算式;另外,像 (R2+R1∗cos⁡θ)(R_2 + R_1 * \cos\theta)(R2​+R1​∗cosθ) 这样的值也会预先计算,并在表达式中重复利用。

明暗

现在我们已经知道该把点放在哪里了,但还不知道每个点的明暗。
计算亮度所需要的数学概念是 法线surface normal。

求出这个法线之后,就可以计算它与光照方向之间的点积,而这个点积中包含了法线与光照方向夹角的余弦值。

因为向量 aaa 与向量 bbb 的点积为 a⋅b=∣a∣∣b∣cos⁡θa \cdot b = |a||b|\cos\thetaa⋅b=∣a∣∣b∣cosθ,所以
如果这个点积大于 0,那么 cos⁡θ>0\cos\theta > 0cosθ>0,说明该位置的点正朝向光源,因此会更亮;如果小于 0,则说明该位置的点没有朝向光源,因此会更暗。
也就是说,点积越大,落在该点上的光就越多。

不过仔细想想,这个法线向量其实非常容易求。
因为从 (0,0,0)(0, 0, 0)(0,0,0) 到最终 (x,y,z)(x, y, z)(x,y,z) 的向量,就会成为点 (x,y,z)(x, y, z)(x,y,z) 的法线向量。 因此,可以像上面的旋转变换矩阵一样重新写出式子。

(Nx,Ny,Nz)=(cos⁡θ,sin⁡θ,0)⋅(cos⁡ϕ0sin⁡ϕ010−sin⁡ϕ0cos⁡ϕ)⋅(1000cos⁡Asin⁡A0−sin⁡Acos⁡A)⋅(cos⁡Bsin⁡B0−sin⁡Bcos⁡B0001)(N_x,N_y,N_z)=(\cos\theta,\sin\theta,0)\cdot\left(\begin{matrix}\cos\phi&0&\sin\phi\\0&1&0\\-\sin\phi&0&\cos\phi\end{matrix}\right)\cdot\left(\begin{matrix}1&0&0\\0&\cos A&\sin A\\0 & -\sin A&\cos A\end{matrix}\right)\cdot\left(\begin{matrix}\cos B&\sin B&0\\-\sin B&\cos B&0 \\0&0&1\end{matrix}\right)(Nx​,Ny​,Nz​)=(cosθ,sinθ,0)⋅​cosϕ0−sinϕ​010​sinϕ0cosϕ​​⋅​100​0cosA−sinA​0sinAcosA​​⋅​cosB−sinB0​sinBcosB0​001​​

这个式子里不同的一点在于,起点是 (cos⁡θ,sin⁡θ,0)(\cos\theta, \sin\theta, 0)(cosθ,sinθ,0)。
这表示的不是精确的 xxx、yyy、zzz 值,而只是各个 θ\thetaθ 对应的余弦和正弦所表示的方向。 (大小并不精确。)

好,那么接下来就该确定光照方向了。
让光从观察者后上方照过来怎么样?
那样的话,光向量就是 (0,1,−1)(0, 1, -1)(0,1,−1)。
现在来求刚才那个法线向量与这个光向量的点积。

L=(Nx,Ny,Nz)⋅(0,1,−1)=cos⁡ϕcos⁡θsin⁡B−cos⁡Acos⁡θsin⁡ϕ−sin⁡Asin⁡θ+cos⁡B(cos⁡Asin⁡θ−cos⁡θsin⁡Asin⁡ϕ)\begin{align}L&=(N_x,N_y,N_z)\cdot(0,1,-1)\\&=\cos\phi\cos\theta\sin B-\cos A\cos\theta\sin\phi-\sin A\sin\theta+\cos B(\cos A\sin\theta-\cos\theta\sin A\sin\phi)\end{align}L​=(Nx​,Ny​,Nz​)⋅(0,1,−1)=cosϕcosθsinB−cosAcosθsinϕ−sinAsinθ+cosB(cosAsinθ−cosθsinAsinϕ)​​

好,最后要做的就是确定常数。
尚未确定的常数有决定甜甜圈整体大小与圆环厚度的 R1R_1R1​、R2R_2R2​,以及表示尺寸的 K1K_1K1​、表示观察者与甜甜圈之间距离的 K2K_2K2​。

首先,将 R1R_1R1​ 和 R2R_2R2​ 分别设为 1 和 2,K2K_2K2​ 设为 5。
至于表示大小的 K1K_1K1​,它会随着屏幕大小(即一屏能显示的 ASCII 字符数)而变化,因此可以写出如下公式。

先考虑甜甜圈可能达到的最大横向宽度(最大 x 坐标)是 R1+R2R_1+R_2R1​+R2​,而此时 zzz 为 0。
如果想让这个点出现在屏幕横向大约 38\dfrac{3}{8}83​ 的位置,那么:

screen_width * 3/8 = K1 * (R1+R2) / (K2 + 0)
screen_width * K2 * 3 / (8 * (R1+R2)) = K1

可以写成上面的形式。

最终代码

因此,用 Python 和 numpy 表示的最终代码如下。

import numpy as np
screen_size = 40
theta_spacing = 0.07
phi_spacing = 0.02
illumination = np.fromiter(".,-~:;=!*#$@", dtype="<U1") 
A = 1
B = 1
R1 = 1
R2 = 2
K2 = 5
K1 = screen_size * K2 * 3 / (8 * (R1 + R2)) 
def render_frame(A: float, B: float) -> np.ndarray: 
    cos_A = np.cos(A) 
    sin_A = np.sin(A) 
    cos_B = np.cos(B) 
    sin_B = np.sin(B) 
    
    output = np.full((screen_size, screen_size), " ")
    zbuffer = np.zeros((screen_size, screen_size))
    
    cos_phi = np.cos(phi := np.arange(0, 2 * np.pi, phi_spacing)) # (315,) 
    sin_phi = np.sin(phi)
    cos_theta = np.cos(theta := np.arange(0, 2 * np.pi, theta_spacing)) # (90,)
    sin_theta = np.sin(theta)
    
    circle_x = R2 + R1 * cos_theta
    circle_y = R1 * sin_theta
    
    x = (np.outer(cos_B * cos_phi + sin_A * sin_B * sin_phi, circle_x) - circle_y * cos_A * sin_B).T
    y = (np.outer(sin_B * cos_phi - sin_A * cos_B * sin_phi, circle_x) + circle_y * cos_A * cos_B).T
    z = ((K2 + cos_A * np.outer(sin_phi, circle_x)) + circle_y * sin_A).T
    ooz = np.reciprocal(z) 
    xp = (screen_size / 2 + K1 * ooz * x).astype(int)
    yp = (screen_size / 2 - K1 * ooz * y).astype(int)
    L1 = (((np.outer(cos_phi, cos_theta) * sin_B) - cos_A * np.outer(sin_phi, cos_theta)) - sin_A * sin_theta)
    L2 = cos_B * (cos_A * sin_theta - np.outer(sin_phi, cos_theta * sin_A))
    L = np.around(((L1 + L2) * 8)).astype(int).T
    mask_L = L >= 0 
    chars = illumination[L] 
    for i in range(90): 
        mask = mask_L[i] & (ooz[i] > zbuffer[xp[i], yp[i]]) 
        zbuffer[xp[i], yp[i]] = np.where(mask, ooz[i], zbuffer[xp[i], yp[i]])     
        output[xp[i], yp[i]] = np.where(mask, chars[i], output[xp[i], yp[i]]) 
        return output 
def pprint(array: np.ndarray) -> None:
    print(*[" ".join(row) for row in array], sep="\n") 
    
if __name__ == "__main__": 
    for _ in range(screen_size * screen_size): 
        A += theta_spacing 
        B += phi_spacing 
        print("\x1b[H") 
        pprint(render_frame(A, B))

Summary: 这篇文章解释了在终端中渲染 3D ASCII 甜甜圈的 donut.c 的数学实现原理。内容涵盖了三维坐标到二维的投影、通过旋转变换矩阵生成甜甜圈形状,以及利用光照点积来表现明暗的方法。