甜甜圈 | Kohath's Reflections

之前偶尔刷到 Joma 的视频 why you NEED math for programming：在终端渲染了一个甜甜圈¹，当时只是不明觉厉。最近回忆起来，有点好奇就去研究了下。

他所做的事是：以固定视角，在甜甜圈表面绘制大量对应每个点亮度的像素点，让其看起来是立体的。这些像素点是 ASCII 字符：.,-~:;=!*#$@ 依次代表最暗到最亮。

主要步骤如下²：

Create a circle of radius R1 centered at R2
Create a donut by rotating about the Y axes
Spin the donut around the X and Z axes
Project donut onto 2D screen
Determine illumination by calculating surface normal (given a light source)

效果图：

donut

首先我们要有一个甜甜圈（也就是环面），其本质上是一个旋转体。我们以 3D 空间中的点为圆心画一个 2D 的圆，然后让圆绕着环面的中心轴旋转就得到了一个旋转体。

sphere_like_degenerate_torus

这是穿过环面中心的横切面：

torus_cross_section

在 xy 平面上，点 (x, y) 绕³圆心 $(R_2, 0, 0)$ 旋转一圈，得到一个半径为 $R_1$ 的圆，假设旋转角度为 $\theta$，则圆上每个点的坐标为：

$$ (x, y, z) = (R_2, 0, 0) + (R_1 cos\theta, R_1 sin\theta, 0) $$

再让圆绕 y 轴旋转另一个角度 $\phi$，则环面上每个点的坐标为⁴：

$$ (R_2 + R_1 \cos\theta, R_1 \sin\theta, 0) \cdot \begin{pmatrix} \cos\phi & 0 & \sin\phi \\ 0 & 1 & 0 \\ -\sin\phi & 0 & \cos\phi \\ \end{pmatrix} $$

如果我们想要动画中整个甜甜圈绕着两个轴转动，假设绕 x 轴转动的角度是 A，绕 y 轴转动的角度是 B，则环面上每个点的坐标为：

$$ (R_2 + R_1 cos\theta, R_1 sin\theta, 0) \cdot \begin{pmatrix} cos\phi & 0 & sin\phi \\ 0 & 1 & 0 \\ -sin\phi & 0 & cos\phi \\ \end{pmatrix} \cdot \begin{pmatrix} 1 & 0 & 0 \\ 0 & cosA & sinA \\ 0 & -sinA & cosA \\ \end{pmatrix} \cdot \begin{pmatrix} cosB & sinB & 0 \\ -sinB & cosB & 0 \\ 0 & 0 & 1\\ \end{pmatrix} $$

把上面的矩阵乘开，我们就获得了，以原点为中心，绕着 x、y 轴旋转的甜甜圈上所有的点：

$$ \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) $$

获得甜甜圈后，我们要把它显示在屏幕上。这是侧视图，一个人在 2D 屏幕前观看后面的 3D 物体：

perspective

假设人到屏幕的距离是 $K_1$，人到物体的距离是 $K_2$，为了在 2D 屏幕上渲染 3D 物体，我们将 3D 空间中的每个点 (x, y, z)，映射到屏幕 (x’, y’, z’)上。人眼、屏幕和人眼、Y 轴形成了两个三角形，它们的关系是：

$$ \frac {y’}{y} = \frac {K_1}{K_2 + z} $$

推导出：

$$ y’ = \frac {K_1 y}{K_2 + z} $$

同理可得 x’（想象俯视图），所以屏幕上每个点的坐标为：

$$ (x’, y’) = (\frac {K_1 x}{K_2 + z}, \frac {K_1 y}{K_2 + z}) $$

当我们绘制大量点时，可能在相同位置 (x’, y’) 上绘制甜甜圈上不同位置的点。它们离屏幕的距离（深度）不同，所以我们要维护一个深度缓冲，来存储我们绘制的点的 Z 坐标⁵，在绘制前判断当前位置是否绘制过。

知道像素点的绘制位置后，我们还需要根据每个点的亮度绘制不同的像素点。假设有一个光源⁶：

light-equation

N 是需要计算亮度的点的法线向量，L 是光到达表面方向的反方向向量，它们都是长度为 1 的单位向量，所以 N 与 L 的点积就是这两个向量夹角的余弦值：如果点积大于 0，表面朝向光源；如果点积小于 0，表面背向光源。点积的值越大，表面越亮。

甜甜圈表面法线方向的推导和我们获取甜甜圈的推导类似，当 $\theta、\phi、A、B$ 相同时，圆上法线方向与以原点为中心的单位圆（半径为 1）相同。

单位圆上的起始点为 $(cos\theta, sin\theta, 0)$，进行相同的转动，表面法线 $(N_x, N_y, N_z)$ 为：

$$ (N_x, N_y, N_z) = (cos\theta, sin\theta, 0) \cdot \begin{pmatrix} cos\phi & 0 & sin\phi \\ 0 & 1 & 0 \\ -sin\phi & 0 & cos\phi \\ \end{pmatrix} \cdot \begin{pmatrix} 1 & 0 & 0 \\ 0 & cosA & sinA \\ 0 & -sinA & cosA \\ \end{pmatrix} \cdot \begin{pmatrix} cosB & sinB & 0 \\ -sinB & cosB & 0 \\ 0 & 0 & 1\\ \end{pmatrix} $$

这也是甜甜圈的表面法线，再将光源放置在人的后上方 (0, 1, -1)，就可以计算出亮度：

$$ L = (N_x, N_y, N_z) \cdot (0, 1, -1) $$

$$ =cos\phi cos\theta sinB - cosA cos\theta sin\phi - sinA sin\theta + cosB (cosA sin\theta - cos\theta sinA sin\phi) $$

我们只考虑表面朝向光源的情况（$L \in [-\sqrt{2}, \sqrt{2}]$⁷，取 L > 0），为了将 L 与表示亮度的 11 个 ASCII 字符一一映射，将 $L \times 8$，值域扩大为 [0, 11]。

齐活！开始写代码⁸：

package main

import (
    "fmt"
    "math"
    "time"
)

const (
    screenWidth  = 80
    screenHeight = 80

    thetaSpacing = 0.07
    phiSpacing   = 0.02

    R1 = 1
    R2 = 2
    // Calculate K1 based on screen size: the maximum x-distance
    // occurs roughly at the edge of the torus, which is at x=R1+R2, z=0.  we
    // want that to be displaced 3/8ths of the width of the screen, which
    // is 3/4th of the way from the center to the side of the screen.
    // screen_width*3/8 = K1*(R1+R2)/(K2+0)
    // screen_width*K2*3/(8*(R1+R2)) = K1
    K1 = screenWidth * K2 * 3 / (8 * (R1 + R2))
    K2 = 5
)

func RenderFrame(A, B float64) {
    // precompute sines and cosines of A and B
    sinA, cosA := math.Sin(A), math.Cos(A)
    sinB, cosB := math.Sin(B), math.Cos(B)

    output := [screenWidth][screenHeight]byte{}
    zBuffer := [screenWidth][screenHeight]float64{}
    for i := 0; i < screenWidth; i++ {
    for j := 0; j < screenWidth; j++ {
        output[i][j] = ' '
        zBuffer[i][j] = 0
        }
    }

    // theta goes around the cross-sectional circle of a torus
    for theta := 0.0; theta < 2*math.Pi; theta += thetaSpacing {
    // precompute sines and cosines of theta
        sinTheta, cosTheta := math.Sin(theta), math.Cos(theta)

        // phi goes around the center of revolution of a torus
        for phi := 0.0; phi < 2*math.Pi; phi += phiSpacing {
            // precompute sines and cosines of phi
            sinPhi, cosPhi := math.Sin(phi), math.Cos(phi)

            // the x,y coordinate of the circle, before revolving 
            // (factored out of the above equations)
            circleX := R2 + R1*cosTheta
            circleY := R1 * sinTheta

            // final 3D (x,y,z) coordinate after rotations, directly from
            // our math above
            x := (circleX*(cosB*cosPhi+sinA*sinB*sinPhi) - circleY*cosA*sinB)
            y := (circleX*(cosPhi*sinB-cosB*sinA*sinPhi) + circleY*cosA*cosB)
            z := K2 + cosA*circleX*sinPhi + circleY*sinA // 实际是 z + K2
            ooz := 1 / z                                 // "one over z"

            // 在 x，y 坐标系统，x‘，y' 可能为负值，向右上方平移，将坐标全部转化为
            // 正值，方便存储在二维数组中
            xp := int(screenWidth/2 + 0.5*K1*ooz*x) // 这里乘上系数 0.5
            yp := int(screenWidth/2 + K1*ooz*y)

            // calculate luminance.  ugly, but correct.
            L := cosPhi*cosTheta*sinB - cosA*cosTheta*sinPhi - sinA*sinTheta + cosB*(cosA*sinTheta-cosTheta*sinA*sinPhi)
            if L > 0 {
            // test against the z-buffer.  larger 1/z means the pixel is
            // closer to the viewer than what's already plotted.
            if ooz > zBuffer[xp][yp] {
                zBuffer[xp][yp] = ooz
                luminanceIndex := int(L * 8)
                // luminance_index is now in the range 0..11
                // now we lookup the character corresponding to the
                // luminance and plot it in our output:
                output[xp][yp] = ".,-~:;=!*#$@"[luminanceIndex]
                }
            }
        }
    }

    // now, dump output[] to the screen.
    // bring cursor to "home" location, in just about any currently-used
    // terminal emulation mode
    fmt.Printf("\x1b[H") // 清空屏幕
    for i := 0; i < screenWidth; i++ {
        for j := 0; j < screenHeight; j++ {
            fmt.Printf("%c", output[i][j])
        }
        fmt.Println()
    }
}

func main() {
    var A, B float64 = 1.0, 1.0
    for {
        A += 0.07
        B += 0.03
        RenderFrame(A, B)
        time.Sleep(100 * time.Millisecond)
    }
}

原文是 Donut math: how donut.c works，Joma 用视频演绎了一遍。 ↩︎
Donut-shaped C code that generates a 3D spinning donut ↩︎
右手系 ↩︎
旋转矩阵(rotation matrix)、绕任意轴的旋转矩阵 ↩︎
实际存的是 $\frac{1}{z}$：$\frac{1}{z}=0$ 可以表示无穷远，$\frac{1}{z}$ 也可以复用。 ↩︎
Introduction to Lighting ↩︎
L $=\vec{N}\cdot\vec{L}=\lvert\vec{N}\rvert\lvert\vec{L}\rvert cos\alpha，其中\lvert\vec{N}\rvert = 1，\lvert\vec{L}\rvert = \sqrt{2}，cos\alpha \in [-1, 1]$ ↩︎
代码从原文中的 C 重写成 Go，但 x’ 乘上了系数 0.5，不然甜甜圈不圆：二维数组的第一维表示列，即纵向，将纵向收缩。 ↩︎