こういうもの、一度くらいは見たことがあるのではないか？

黒い画面で、白い特殊文字だけを使って回転する3Dドーナツを作る。
たしかにシンプルなプログラムに見えるが、不思議と見入ってしまう魅力がある。
これがDonut.cだ。
このようにASCII文字列で絵を描くことをASCII Artという。

「いや、しょせんコンソール画面に文字を出力しているだけでしょ。」

しかし、簡単に見てはいけない。
グラフィックスに関わるだけあって、一度掘り下げてみると深遠な数学の世界が広がっているからだ。

この記事では、その数学を一度理解してみようと思う。

Donut.c

このコンソールドーナツへの私の関心は、あるYouTube動画から始まった。
その動画は**"why you NEED math for programming"**、訳すと「プログラミングに数学が必要な理由」という動画だ。

この記事の最後では、ドーナツがどのようにASCIIコードでレンダリングされるのかを、簡潔な数学的原理で説明してくれる。
それを聞いて興味が湧き、インターネットで検索してみることになった。

驚いたことに、思ったよりこれについて書かれた記事は多くなかった。
せいぜい簡単な韓国語の紹介記事と、複雑な英語の数学的原理を説明する記事があるだけだった。

最初は英語に複雑な数学まで混ざっていて読む気がしなかったが、覚悟を決めて読んでみると、意外と読めるものだった。

では、ここから本当に分析を始めてみよう。
ここで出てくる数学的原理はこの記事を参考にした。

ASCII Artについて

まず、一般的なASCII Artではどのように明るさを表現するのかを見てみよう。

私たちがASCII Artを出力する場所はターミナルだが、ターミナルは一般的に背景が暗い色で、文字が明るい色だ。
したがって、1文字が占める空間の中でピクセルの密度が高ければ明るい色、ピクセルの密度が低ければ暗い色と定めることができる。
これを簡単な12文字で表したのが以下の文字たちだ。
.,-~:;=!*#$@
左側は一定の空間内で密度が最も低いので最も暗い色、右側は一定の空間内で密度が最も高いので最も明るい色である。

数学的原理

まず、あれをどう実装すればよいのかを考えてみよう。
必要なのは2つある。

1. ドーナツの形状を描くこと
2. ドーナツの形状に明るさを表すこと

さあ、ここから3D一人称視点で使われる簡単な数学から始めてみよう。

上の図は、人がスクリーンの前に座って、スクリーンの向こう側にある物体を見ている図である。 3次元の物体を2次元に描くために、3次元の各点 $(x, y, z)$(x,y,z) を、観察者から $z'$z′ だけ離れた平面に投影する。
すると、その点は $(x', y')$(x′,y′) になる。
3次元物体の各点の座標が分かっている状況で、2次元のどこに投影すべきかを知る必要があるので、$x$x, $y$y を使って $x'$x′, $y'$y′ を求めなければならない。 ($z'$z′ は定数なので $K_1$K1​ としよう。)

数式で表すと $\dfrac{y'}{K_1} = \dfrac{y}{z}$K1​y′​=zy​
$z'$z′ を移項すると $y' = \dfrac{yK_1}{z}$y′=zyK1​​
この式は $x$x にも同じように適用できる: $x' = \dfrac{xK_1}{z}$x′=zxK1​​
したがって投影方程式は $(x', y')=(\dfrac{K_1x}{z}, \dfrac{K_1y}{z})$(x′,y′)=(zK1​x​,zK1​y​) となる。
$K_1$K1​ は2Dウィンドウに表示したい視野角に応じて任意に設定できる。たとえば100x100サイズのウィンドウがあるなら、視野は(50, 50)に合わせられているはずだ。このとき3D空間で視野を基準に5だけ離れていて幅が10の物体を見たいなら、点 $x=10$x=10, $z=5$z=5 の投影が画面内に現れるようにしなければならないので:

• $x'<50$x′<50 でなければならない
• つまり $\dfrac{10K_1}{z}<50$z10K1​​<50 でなければならない
• つまり $K_1 < 25$K1​<25 でなければならない

また、多くの点を投影していると、$x$x と $y$y は同じだが $z$z が異なる点を投影しなければならないこともあるので、私たちが描くすべてのものについての $z$z 座標を保存する z-buffer を作る。z-bufferがあれば、ある位置に点を表示するとき、その位置にすでにあったものより前にあるかどうかを確認できる。
さらにz-bufferは $z^{-1}=\dfrac{1}{z}$z−1=z1​ を計算して深度をバッファリングするのにも役立つが、その理由は次の通りだ。

1. $z^{-1} = 0$z−1=0 なら $z$z は無限大に近いので、z-bufferの初期値を0にして背景を無限の深さに設定できる。
2. $x'$x′ と $y'$y′ を求めるときに再利用できる。 ($z^{-1}=\dfrac{1}{z}$z−1=z1​、あらかじめ $z$z で割っておくことで、わざわざ $z$z を2回割らずに少し効率よく計算できる。)

ドーナツ (solid of revolution)

では次に、ドーナツの形状をどう作るかを考えなければならない。
幸いなことに、ドーナツは回転体(solid of revolution) なので、次のような原理で各点を求めることができる。
見ての通り、$(R_2, 0, 0)$(R2​,0,0) を中心として、半径 $R_1$R1​ の円がある。
$R_1$R1​ を半径とする円は0〜360度回転させることで描ける。
0〜360度まで変化する変数を $\theta$θ としよう。

これを式で表すと次のようになる。
$(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)
この式で1つの円を求めることができる。

次にこの円を2y軸に沿って回転させてみよう。すぐにドーナツ形が現れることに気づいただろう。
この円をy軸に沿って回転させるには、回転変換行列を用いれば求められる。

y軸に沿って回転させるので、3次元のy軸回転変換行列を持ってきて掛けると、次の式になる。

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

しかし、私たちはこのドーナツが 2つ以上の軸で回転するように 実装したいので、さらに2つの回転変換行列を持ってきて掛ける。
また、持ってきた2つの回転変換行列で使われる角度はAとBと呼ぶ。
式は次の通りだ。

$$(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) を基準として、角度 $A$A と $B$B によって一定に回転するドーナツの座標を手に入れた。

次にこれを実際に画面へ投影するには、観察者の座標が $(0, 0, 0)$(0,0,0) だと仮定して、ドーナツを観察者の前方へ移動させなければならない。
ドーナツを観察者の前方へ移すために、$z$z に 観察者とドーナツの間の距離 を足す必要がある。
この 観察者とドーナツの間の距離 を $K_2$K2​ としよう。

さて、実際の $(x', y')$(x′,y′) を求める式は次のようになる。

$(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) を求める式は次の通りだ。

$$\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θ​​

できるだけ最適化するために、プログラミング段階では行列積を使わず、上の計算式を使うことにする。また、$(R_2 + R_1 * \cos\theta)$(R2​+R1​∗cosθ) のような値をあらかじめ計算して、式の中で再利用するようにする。

明るさ

これでどこに点を置くべきかは分かったが、それぞれの点の明るさはまだ分からない。
輝度を計算するのに必要な数学的概念は 法線surface normal だ。

この法線を求めると、法線と光の方向の内積を求められるようになり、この内積には法線と光の方向のなす角のコサイン値が含まれている。

ベクトル $a$a とベクトル $b$b の内積は $a \cdot b = |a||b|\cos\theta$a⋅b=∣a∣∣b∣cosθ なので
もしこの内積が0より大きければ $\cos\theta > 0$cosθ>0 であり、その位置の点が光を向いているので明るく、0より小さければその位置の点は光を向いていないので暗い。
つまり、内積が大きいほど、その点にはより多くの光が当たる。

ところで、考えてみるとこの法線ベクトルは非常に簡単に求められる。
なぜなら、$(0, 0, 0)$(0,0,0) から最終的な $(x, y, z)$(x,y,z) までのベクトルが、$(x, y, z)$(x,y,z) 点の法線ベクトルになるからだ。 したがって、上の回転変換行列と同じように式を書き直すことができる。

$$(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\theta, \sin\theta, 0)$(cosθ,sinθ,0) だということだ。
これは正確な $x$x, $y$y, $z$z の値ではなく、各 $\theta$θ に対するコサインとサインによって方向だけを表したものだ。 (大きさは正確ではない。)

さて、そうなると今度は光の方向を決めなければならない。
観察者の後ろ側の上方から光を当てるのはどうだろうか？
そうすると光のベクトルは $(0, 1, -1)$(0,1,−1) になる。
では先ほどの法線ベクトルとこの光ベクトルの内積を求めてみよう。

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

さて、最後にやるべきことは定数を決めることだ。
まだ決めていない定数には、ドーナツ全体の大きさとリングの太さを決める $R_1$R1​, $R_2$R2​ があり、さらに大きさを表す $K_1$K1​、観察者とドーナツの距離を表す $K_2$K2​ がある。

まず $R_1$R1​ と $R_2$R2​ はそれぞれ1と2に決め、$K_2$K2​ は5に決めた。
$K_1$K1​ は大きさを表すが、画面の大きさ(1画面に表示できるASCII文字数)によって変わるので、次のように式を立てることができる。

まずドーナツが取りうる最大の横幅(最大x座標)は $R_1+R_2$R1​+R2​ であり、$z$z は0である。
その点が画面の横方向の $\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の数学的な実装原理を説明します。3次元座標の2次元への投影、回転変換行列によるドーナツ形状の生成、そして光との内積を用いた明るさ表現の方法を扱います。