甜甜圈
之前偶尔刷到 Joma 的视频 why you NEED math for programming：在终端渲染了一个甜甜圈1，当时只是不明觉厉。最近回忆起来，有点好奇就去研究了下。
他所做的事是：以固定视角，在甜甜圈表面绘制大量对应每个点亮度的像素点，让其看起来是立体的。这些像素点是 ASCII 字符：.,-~:;=!*#$@ 依次代表最暗到最亮。
主要步骤如下2：
- 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)
效果图：
首先我们要有一个甜甜圈（也就是环面），其本质上是一个旋转体。我们以 3D 空间中的点为圆心画一个 2D 的圆，然后让圆绕着环面的中心轴旋转就得到了一个旋转体。
这是穿过环面中心的横切面：
在 xy 平面上，点 (x, y) 绕3圆心 $(R_2, 0, 0)$ 旋转一圈，得到一个半径为 $R_1$ 的圆，假设旋转角度为 $\theta$，则圆上每个点的坐标为：
再让圆绕 y 轴旋转另一个角度 $\phi$，则环面上每个点的坐标为4：
如果我们想要动画中整个甜甜圈绕着两个轴转动，假设绕 x 轴转动的角度是 A，绕 y 轴转动的角度是 B，则环面上每个点的坐标为：
把上面的矩阵乘开，我们就获得了，以原点为中心，绕着 x、y 轴旋转的甜甜圈上所有的点：
获得甜甜圈后，我们要把它显示在屏幕上。这是侧视图，一个人在 2D 屏幕前观看后面的 3D 物体：
假设人到屏幕的距离是 $K_1$，人到物体的距离是 $K_2$，为了在 2D 屏幕上渲染 3D 物体，我们将 3D 空间中的每个点 (x, y, z)，映射到屏幕 (x’, y’, z’)上。人眼、屏幕和人眼、Y 轴形成了两个三角形，它们的关系是：
推导出：
同理可得 x’（想象俯视图），所以屏幕上每个点的坐标为：
当我们绘制大量点时，可能在相同位置 (x’, y’) 上绘制甜甜圈上不同位置的点。它们离屏幕的距离（深度）不同，所以我们要维护一个深度缓冲，来存储我们绘制的点的 Z 坐标5，在绘制前判断当前位置是否绘制过。
知道像素点的绘制位置后，我们还需要根据每个点的亮度绘制不同的像素点。假设有一个光源6：
N 是需要计算亮度的点的法线向量，L 是光到达表面方向的反方向向量，它们都是长度为 1 的单位向量，所以 N 与 L 的点积就是这两个向量夹角的余弦值：如果点积大于 0，表面朝向光源；如果点积小于 0，表面背向光源。点积的值越大，表面越亮。
甜甜圈表面法线方向的推导和我们获取甜甜圈的推导类似，当 $\theta、\phi、A、B$ 相同时，圆上法线方向与以原点为中心的单位圆（半径为 1）相同。
单位圆上的起始点为 $(cos\theta, sin\theta, 0)$，进行相同的转动，表面法线 $(N_x, N_y, N_z)$ 为：
这也是甜甜圈的表面法线，再将光源放置在人的后上方 (0, 1, -1)，就可以计算出亮度：
我们只考虑表面朝向光源的情况（$L \in [-\sqrt{2}, \sqrt{2}]$7，取 L > 0），为了将 L 与表示亮度的 11 个 ASCII 字符一一映射，将 $L \times 8$，值域扩大为 [0, 11]。
齐活！开始写代码8：
1package main
2
3import (
4 "fmt"
5 "math"
6 "time"
7)
8
9const (
10 screenWidth = 80
11 screenHeight = 80
12
13 thetaSpacing = 0.07
14 phiSpacing = 0.02
15
16 R1 = 1
17 R2 = 2
18 // Calculate K1 based on screen size: the maximum x-distance
19 // occurs roughly at the edge of the torus, which is at x=R1+R2, z=0. we
20 // want that to be displaced 3/8ths of the width of the screen, which
21 // is 3/4th of the way from the center to the side of the screen.
22 // screen_width*3/8 = K1*(R1+R2)/(K2+0)
23 // screen_width*K2*3/(8*(R1+R2)) = K1
24 K1 = screenWidth * K2 * 3 / (8 * (R1 + R2))
25 K2 = 5
26)
27
28func RenderFrame(A, B float64) {
29 // precompute sines and cosines of A and B
30 sinA, cosA := math.Sin(A), math.Cos(A)
31 sinB, cosB := math.Sin(B), math.Cos(B)
32
33 output := [screenWidth][screenHeight]byte{}
34 zBuffer := [screenWidth][screenHeight]float64{}
35 for i := 0; i < screenWidth; i++ {
36 for j := 0; j < screenWidth; j++ {
37 output[i][j] = ' '
38 zBuffer[i][j] = 0
39 }
40 }
41
42 // theta goes around the cross-sectional circle of a torus
43 for theta := 0.0; theta < 2*math.Pi; theta += thetaSpacing {
44 // precompute sines and cosines of theta
45 sinTheta, cosTheta := math.Sin(theta), math.Cos(theta)
46
47 // phi goes around the center of revolution of a torus
48 for phi := 0.0; phi < 2*math.Pi; phi += phiSpacing {
49 // precompute sines and cosines of phi
50 sinPhi, cosPhi := math.Sin(phi), math.Cos(phi)
51
52 // the x,y coordinate of the circle, before revolving
53 // (factored out of the above equations)
54 circleX := R2 + R1*cosTheta
55 circleY := R1 * sinTheta
56
57 // final 3D (x,y,z) coordinate after rotations, directly from
58 // our math above
59 x := (circleX*(cosB*cosPhi+sinA*sinB*sinPhi) - circleY*cosA*sinB)
60 y := (circleX*(cosPhi*sinB-cosB*sinA*sinPhi) + circleY*cosA*cosB)
61 z := K2 + cosA*circleX*sinPhi + circleY*sinA // 实际是 z + K2
62 ooz := 1 / z // "one over z"
63
64 // 在 x，y 坐标系统，x‘，y' 可能为负值，向右上方平移，将坐标全部转化为
65 // 正值，方便存储在二维数组中
66 xp := int(screenWidth/2 + 0.5*K1*ooz*x) // 这里乘上系数 0.5
67 yp := int(screenWidth/2 + K1*ooz*y)
68
69 // calculate luminance. ugly, but correct.
70 L := cosPhi*cosTheta*sinB - cosA*cosTheta*sinPhi - sinA*sinTheta + cosB*(cosA*sinTheta-cosTheta*sinA*sinPhi)
71 if L > 0 {
72 // test against the z-buffer. larger 1/z means the pixel is
73 // closer to the viewer than what's already plotted.
74 if ooz > zBuffer[xp][yp] {
75 zBuffer[xp][yp] = ooz
76 luminanceIndex := int(L * 8)
77 // luminance_index is now in the range 0..11
78 // now we lookup the character corresponding to the
79 // luminance and plot it in our output:
80 output[xp][yp] = ".,-~:;=!*#$@"[luminanceIndex]
81 }
82 }
83 }
84 }
85
86 // now, dump output[] to the screen.
87 // bring cursor to "home" location, in just about any currently-used
88 // terminal emulation mode
89 fmt.Printf("\x1b[H") // 清空屏幕
90 for i := 0; i < screenWidth; i++ {
91 for j := 0; j < screenHeight; j++ {
92 fmt.Printf("%c", output[i][j])
93 }
94 fmt.Println()
95 }
96}
97
98func main() {
99 var A, B float64 = 1.0, 1.0
100 for {
101 A += 0.07
102 B += 0.03
103 RenderFrame(A, B)
104 time.Sleep(100 * time.Millisecond)
105 }
106}
原文是 Donut math: how donut.c works，Joma 用视频演绎了一遍。 ↩︎
右手系 ↩︎
实际存的是 $\frac{1}{z}$：$\frac{1}{z}=0$ 可以表示无穷远，$\frac{1}{z}$ 也可以复用。 ↩︎
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，不然甜甜圈不圆：二维数组的第一维表示列，即纵向，将纵向收缩。 ↩︎