..

甜甜圈

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

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

主要步骤如下2：

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

效果图：









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



sphere_like_degenerate_torus



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



torus_cross_section



在 xy 平面上，点 (x, y) 绕3圆心 \((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\)，则环面上每个点的坐标为4：


\[(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 坐标5，在绘制前判断当前位置是否绘制过。






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



light-equation



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



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



torus_cross_section_2.svg



单位圆上的起始点为 \((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}]\)7，取 L > 0），为了将 L 与表示亮度的 11 个 ASCII 字符一一映射，将 \(L \times 8\)，值域扩大为 [0, 11]。






齐活！开始写代码8：



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)
    }
}








  
    
    
  1. 
      
    原文是 Donut math: how donut.c works，Joma 用视频演绎了一遍。 ↩
    
    
    2. 
    
  2. 
      
    Donut-shaped C code that generates a 3D spinning donut ↩
    
    
    3. 
    
  3. 
      
    右手系 ↩
    
    
    4. 
    
  4. 
      
    旋转矩阵(rotation matrix)、绕任意轴的旋转矩阵 ↩
    
    
    5. 
    
  5. 
      
    实际存的是 \(\frac{1}{z}\)：\(\frac{1}{z}=0\) 可以表示无穷远，\(\frac{1}{z}\) 也可以复用。 ↩
    
    
    6. 
    
  6. 
      
    Introduction to Lighting ↩
    
    
    7. 
    
  7. 
      
    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]\) ↩
    
    
    8. 
    
  8. 
      
    代码从原文中的 C 重写成 Go，但 x’ 乘上了系数 0.5，不然甜甜圈不圆：二维数组的第一维表示列，即纵向，将纵向收缩。 ↩
    
    
    9.