相対論的レイトレーシング

これは、以前ロシア語で公開した記事を少し洗練させたものです。
この記事では、一般相対性理論に完全に準拠して、物理的に正確なブラックホールをレンダリングする方法を示します。そのために、光線の運動方程式を導出し、GLSLでルンゲ＝クッタ積分器を実装し、最終的にカメラから過去への光線追跡を行うフラグメントシェーダーにこれらを統合します。

準備：計量、接続、測地線

光は直線を進みます。しかし、「直線」とは本当は何を意味するのでしょうか？線形方程式で記述されるもの？いいえ、そうではありません。極座標系を取ると、直線の方程式は線形とは程遠いものになります。そして、空間自体が平坦でない場合（このトーラスのように）、直線が線形になるような座標系の選択肢はありません。

［画像: thorus_geo.png|トーラス上の直線］

では、どのように曲線が「直線」であると考えるのでしょうか？ある点での曲線に沿ったベクトルを取り、それを曲線に沿って移動させます。それはまだ曲線に沿った方向を向いているでしょうか？もしあらゆる点でそうであれば、そのような曲線は測地線と呼ばれ、これは任意の多様体上で直線に最も近いものです。

4次元多様体上の測地線の方程式を導出するのは非常に骨の折れる作業です。幸いなことに、手作業で行う必要はありません。機械に骨の折れる作業をやらせましょう。Pythonパッケージのsympyを使用します。

from sympy import *
from sympy.diffgeom import *
import numpy as np

sympy.diffgeomモジュールはやや不便で扱いにくいですが、アトラスで記述された多様体上の方程式を導出する仕事は比較的うまくこなします。私たちの場合はアトラスは1つのパッチしか持ちませんが、いくつかの座標系を持ちます。1つは球面座標系（シュワルツシルト計量で作業する際に便利）、もう1つは将来のカメラ用のデカルト座標系です。座標シンボルを定義しましょう。

t_, x_, y_, z_ = symbols("t x y z")
rho_, theta_, phi_ = symbols("\\rho \\theta \\phi")

（私は終端シンボルを慣習として_サフィックスを使用しています）

球面座標からデカルト座標への変換方程式を構築しましょう。

x = rho_ * sin(theta_) * sin(phi_)
y = rho_ * sin(theta_) * cos(phi_)
z = rho_ * cos(theta_)
print(latex(x_) + " = " + latex(x))
print(latex(y_) + " = " + latex(y))
print(latex(z_) + " = " + latex(z))

$$ \begin{array}{l} x = \rho \sin{\left(\phi \right)} \sin{\left(\theta \right)} \ y = \rho \sin{\left(\theta \right)} \cos{\left(\phi \right)} \ z = \rho \cos{\left(\theta \right)} \ \end{array} $$

新しい4次元多様体spacetimeを作成し、そのアトラスにpatchを追加し、2つの座標系をそのパッチにアタッチします。

spacetime = Manifold("spacetime", 4)
patch = Patch("patch", spacetime)
relation_dict = {
    ("spherical", "cart"): [(x_, y_, z_, t_), (x, y, z, t_)]
}
coord_sh = CoordSystem("spherical", patch, (t_, rho_, theta_, phi_), relation_dict)
coord_cart = CoordSystem("cart", patch, (x_, y_, z_, t_), relation_dict)

sympy.diffgeomの不便な点の一つは、CoordSystem関数に渡された座標シンボルを直接使用できないことです。代わりに、そこから新しいシンボルを取得する必要があります。

c_t, c_rho, c_theta, c_phi = coord_sh.base_scalars()
print(",".join(map(latex, (c_t, c_rho, c_theta, c_phi))))

$$ \mathbf{t},\mathbf{\rho},\mathbf{\theta},\mathbf{\phi} $$

ご覧のとおり、シンボルは同じですが、sympy.diffgeomで使用できるのは新しい変数のセットだけです。

次に、シュワルツシルト計量を定義します。

Rs = symbols("r_s")
d_t, d_rho, d_theta, d_phi = coord_sh.base_oneforms()
TP = TensorProduct
metric = (
    + (1 - Rs/c_rho) * TP(d_t, d_t)
    - 1 / (1 - Rs / c_rho) * TP(d_rho, d_rho)
    - (c_rho**2) * TP(d_theta, d_theta)
    - (c_rho**2) * (sin(c_theta)**2) * TP(d_phi, d_phi)
)
print(latex(metric))

$$ \left(- \frac{r
_s}{\mathbf{\rho}} + 1\right) \operatorname{d}t \otimes \operatorname{d}t - \sin^{2}{\left(\mathbf{\theta} \right)} \mathbf{\rho}^{2} \operatorname{d}\phi \otimes \operatorname{d}\phi - \mathbf{\rho}^{2} \operatorname{d}\theta \otimes \operatorname{d}\theta - \frac{\operatorname{d}\rho \otimes \operatorname{d}\rho}{- \frac{r
_s}{\mathbf{\rho}} + 1} $$

ここで1形式のテンソル積を使用しなければならないのは、もう一つの不便な点です。g_ij成分の2次元配列があればより簡単だと想像するかもしれませんが、残念ながらそうはいきません。結局、微分形式を使います。いずれにせよ、計量をこの形式で入手できれば、それを使ってレヴィ＝チヴィタ接続（別名、2種クリストッフェル記号）の成分を計算できます。多様体上の接続は、ある点の接空間のベクトルを別の点の接空間に連続的に移動させたときにどのように変換されるかを定義します。そして、レヴィ＝チヴィタ接続は、そのような移動の下で計量テンソルを保存する特別な種類の接続です。つまり、計量テンソルの共変微分はゼロですが、その成分は変化する可能性があります。

sympy.diffgeomは、計量テンソル保存条件に基づいて、クリストッフェル記号を導出する関数metric_to_Christoffel_2ndを提供しています。

christoffel = simplify(metric_to_Christoffel_2nd(metric))
print("\\begin{cases}") # v v v intentional reordering here, to sort by lower indices first
for (j,k,i) in np.ndindex(*np.shape(christoffel)):
    if not christoffel[i,j,k].is_zero:
        print(f"\\Gamma^{{{i}}}
_{{{j}{k}}} =" + latex(christoffel[i,j,k]) + " \\\\ ")
print("\\end{cases}")

$$ \begin{cases} \Gamma^{1}
_{00} =\frac{r
_s \left(- r
_s + \mathbf{\rho}\right)}{2 \mathbf{\rho}^{3}} \ \Gamma^{0}
_{01} =\frac{r
_s}{2 \left(- r
_s + \mathbf{\rho}\right) \mathbf{\rho}} \ \Gamma^{0}
_{10} =\frac{r
_s}{2 \left(- r
_s + \mathbf{\rho}\right) \mathbf{\rho}} \ \Gamma^{1}
_{11} =\frac{r
_s}{2 \left(r
_s - \mathbf{\rho}\right) \mathbf{\rho}} \ \Gamma^{2}
_{12} =\frac{1}{\mathbf{\rho}} \ \Gamma^{3}
_{13} =\frac{1}{\mathbf{\rho}} \ \Gamma^{2}
_{21} =\frac{1}{\mathbf{\rho}} \ \Gamma^{1}
_{22} =r
_s - \mathbf{\rho} \ \Gamma^{3}
_{23} =\frac{1}{\tan{\left(\mathbf{\theta} \right)}} \ \Gamma^{3}
_{31} =\frac{1}{\mathbf{\rho}} \ \Gamma^{3}
_{32} =\frac{1}{\tan{\left(\mathbf{\theta} \right)}} \ \Gamma^{1}
_{33} =\left(r
_s - \mathbf{\rho}\right) \sin^{2}{\left(\mathbf{\theta} \right)} \ \Gamma^{2}
_{33} =- \frac{\sin{\left(2 \mathbf{\theta} \right)}}{2} \ \end{cases} $$

接続成分が手に入ったので、測地線の方程式を導出できます。これは手動で行う必要があります。幸いなことに、それほど難しくはありません。ある曲線がパラメータ方程式 $x^\mu(t)$ で定義されているとします。その接線ベクトルは $\dot{x}^\mu(t) = \frac{d}{dt}x^\mu(t)$ の形をしています。測地線の定義によれば、このベクトルを自身に沿って共変微分したものはゼロになるはずです：

$ \nabla
_{\dot{x}} \dot{x} = 0 $

クリストッフェル記号を使って展開すると：

$ \nabla
_{\dot{x}} \dot{x} = \left( \frac{\partial}{\partial{x^\mu}} \dot{x}^\nu + \Gamma^{\nu}
_{\kappa\mu}\dot{x}^\kappa \right) \dot{x}^\mu = \left( \frac{\partial}{\partial{x^\mu}} \frac{dx^\nu}{dt} \right) \frac{dx^\mu}{dt} + \Gamma^{\nu}
_{\kappa\mu} \frac{dx^\kappa}{dt} \frac{dx^\mu}{dt} $

最初の項で、因子の順序を入れ替え、連鎖律を使用します。

$ \left( \frac{\partial}{\partial{x^\mu}} \frac{dx^\nu}{dt} \right) \frac{dx^\mu}{dt} = \frac{dx^\mu}{dt}\frac{\partial}{\partial{x^\mu}} \left( \frac{dx^\nu}{dt} \right) = \frac{d}{dt}\left( \frac{dx^\nu}{dt} \right) = \frac{d^2 x^\nu}{dt^2} $

すると、方程式は次のようになります。

$ \frac{d^2 x^\nu}{dt^2} + \Gamma^{\nu}
_{\kappa\mu} \frac{dx^\kappa}{dt} \frac{dx^\mu}{dt} = 0 $

または、単純に

$ \boxed{ \ddot{x}^\nu = -\Gamma^{\nu}
_{\kappa\mu} \dot{x}^\kappa \dot{x}^\mu } $

Pythonでは：

def geo_second_derivatives(christoffel, derivatives):
    ''' Parameters:
        christoffel: tensor of rank 3 - analytical expression of christoffel symbol in given coordinates
        derivatives: iterable consisting of sympy symbols corresponding to each first derivative of a coordinate variable
    '''
    d2u = np.zeros(len(derivatives), dtype='object')
    for n,k,m in np.ndindex(*np.shape(christoffel)):
        d2u[n] -= christoffel[n,k,m] * derivatives[k] * derivatives[m]
    return d2u

d2_t, d2_rho, d2_theta, d2_phi = geo_second_derivatives(christoffel, symbols("t' \\rho' \\theta' \\phi'"))
print("\\begin{array}{l}")
print(f"{{{latex(c_t)}}}'' = {latex(d2_t)} \\\\ ")
print(f"{{{latex(c_rho)}}}'' = {latex(d2_rho)} \\\\ ")
print(f"{{{latex(c_theta)}}}'' = {latex(d2_theta)} \\\\ ")
print(f"{{{latex(c_phi)}}}'' = {latex(d2_phi)} \\\\ ")
print("\\end{array}")

$$ \begin{array}{l}
{\mathbf{t}}'' = - \frac{\rho' r
_s t'}{\left(- r
_s + \mathbf{\rho}\right) \mathbf{\rho}} \
{\mathbf{\rho}}'' = - \phi'^{2} \left(r
_s - \mathbf{\rho}\right) \sin^{2}{\left(\mathbf{\theta} \right)} - \frac{\rho'^{2} r
_s}{2 \left(r
_s - \mathbf{\rho}\right) \mathbf{\rho}} - \theta'^{2} \left(r
_s - \mathbf{\rho}\right) - \frac{r
_s t'^{2} \left(- r
_s + \mathbf{\rho}\right)}{2 \mathbf{\rho}^{3}} \
{\mathbf{\theta}}'' = \frac{\phi'^{2} \sin{\left(2 \mathbf{\theta} \right)}}{2} - \frac{2 \rho' \theta'}{\mathbf{\rho}} \
{\mathbf{\phi}}'' = - \frac{2 \phi' \theta'}{\tan{\left(\mathbf{\theta} \right)}} - \frac{2 \rho' \phi'}{\mathbf{\rho}} \
\end{array} $$

このようにして、測地線の方程式の微分方程式系に到達しました。

レイトレーシング

［画像: sphere_pole.png|球の極の重複］

光線（または他のあらゆる物体）は、2階微分方程式で定義される「直線」に沿って進みます。Pythonのscipy.integrateを使って、これらの経路を数値的に簡単に計算できます。しかし、ブラックホールを見たいと思いませんか？画面の各ピクセルに対してレイを送信し、それらをブラックホールの周りを回らせて周囲に到達させることができます。そして、レイを過去に送り返せば、過去に放出された光を受信したのと同じになります！しかし、百万個ものレイを送信する必要があるでしょう。Pythonスクリプトにとっては容易なタスクではありません。並列で実行できればどれほど良いか！GPUはまさにそのためにあります。各ピクセルで同じ計算を実行するために、GPUプログラムの一種であるフラグメントシェーダーがあります。独自のシェーダーを作成するために、OpenGLとそのシェーディング言語GLSLを使用しましょう。計画は非常にシンプルです：ブラックホールの周りにカラフルな球殻を配置し、その近くにカメラを置き、カメラから光線を送信し、それらが球殻に当たって対応するピクセルに色を付けるか、事象の地平面に失われるようにします。

GLSLでの測地線方程式…

Sympyから得られた方程式を別の言語に打ち直したい人はいません。幸いなことに、ccode関数を使用してコードを生成し、その後はわずかなクリーンアップだけで済みます。コードを整形すると、測地線方程式 $p'' = f(p, p')$ の次の関数が得られます。

vec4 d2p(vec4 p, vec4 dp) {
    float t = p[0];
    float rho = p[1];
    float theta = p[2];
    float phi = p[3];

    float D_t = dp[0];
    float D_rho = dp[1];
    float D_theta = dp[2];
    float D_phi = dp[3];

    // Geodesics equation:
    float D2_t = - (D_rho * rs * D_t) / (rho * (rho - rs));
    float D2_rho = (rho - rs) * (
        + pow(D_phi * sin(theta), 2)
        + pow(D_theta, 2)
        - (rs * pow(D_t, 2)) / (2 * pow(rho, 3))
    ) + (rs * pow(D_rho, 2)) / (2 * rho * (rho - rs));
    float D2_theta = pow(D_phi, 2) * sin(theta) * cos(theta) - (2 * D_rho * D_theta) / rho;
    float D2_phi = -(2 * D_phi * D_theta) / tan(theta) - (2 * D_rho * D_phi) / rho;

    vec4 d2p = {D2_t, D2_rho, D2_theta, D2_phi};
    return d2p;
}

この関数は、光線の上の点の座標と、その点での4元速度を受け取ります。それは、曲線に沿って移動する際の座標の2階微分を返します。ここで、すぐに問題点が見つかります。tan(theta)による除算を見てください。この角度がゼロになると、2階微分の値は未定義になります。$ heta$が0に近い場合でも、精度が失われます。幸いなことに、これは大きな問題ではありません。適応的積分ステップを使用することで、精度問題に対処します。あるいは、非回転ブラックホールの場合、光線はその平面から離れることがないため、ゼロの問題を回避できることに気づくかもしれませんが、これは読者の演習として残しておきます。

これで、任意の測地線が点pを4元速度dpで与えられた方向に通過する場合の2階微分を計算できます。この2階微分は、4元速度がどのように変化するかを決定し、4元速度は次の点が何になるかを定義します。古典的な数値積分問題です。

...そしてその数値解

精度とパフォーマンスのバランスを見つける探求において、いくつかの異なる積分方法を試しましたが、最終的に3次と4次の適応型ルンゲ＝クッタ法に落ち着きました。適応型ルンゲ＝クッタ法のアイデアは、異なる次数（この場合は3次と4次）の2つの方法が並列で実行されることです。それらの結果の差は誤差の推定値となり、それを使用して積分ステップを調整します。誤差が大きすぎる場合はステップが小さくなり、十分に小さい場合はステップが大きくなります。二重に便利なのは、一方の次数法の中間結果を、もう一方の次数法内で異なる係数で再利用してパフォーマンスを向上させることができることです。これらの係数を選択することはそれ自体が科学であり、このトピックに関する何百もの研究論文が書かれています。当然ながら、私はこれらの何百もの論文を読みませんでした。私は古典的なRK4を主要な方法として採用し、できるだけ多くの計算済み部分を再利用するために3次法をハッキングしました。

[!info]

• RK4法について
  RK4は、初期値問題 $y'(t) = f(t, y(t))$ を満たし、点 $(t_0, y_0)$ を通る曲線 $y(t)$ を見つけるための方法です。それを数値的に解く最も簡単な方法はオイラー法です。
  $y_{i+1} = y_i + f(t_i, y_i) \Delta t$
  残念ながら、2階微分の値が大きいほど、この方法は悪くなります。特に曲線が同じ方向に曲がり続ける場合は非常に悪く、まさに私たちのケースです。
  しかし、それにいくつかの修正を加えることができます。導関数を $(t_i, y_i)$ の点だけでなく、その周りでも計算します。これを考えてみてください。

1. $a = f(t_i,y_i)$ - オイラー法の係数
2. $b=f\left( t_i + \frac{\Delta t}{2}, y_i + a \frac{\Delta t}{2} \right)$ - オイラー法の結果点まで半分進み、そこで導関数値を計算します。
3. この$b$値を使用して、初期点から半ステップ移動します：$y_b = y_i + b \frac{\Delta t}{2}$
4. そこで導関数値を計算します：
   $c=f\left( t_i + \frac{\Delta t}{2}, y_i + b \frac{\Delta t}{2} \right)$
5. 最後に、$c$を使用して初期点からフルステップを取り、そこで別の導関数を計算します：$d = f(t_i + \Delta t, y_i + c \Delta t)$
6. これらの4つの値を使用して、$a,b,c,d$の加重平均を $(t_i, y_i)$ における実際の導関数値とします：$y_i'(t) = \frac{1}{6}(a + 2b + 2c + d)$

そして、この値のみを使用して、実際のステップが取られます。
$y_{i+1} = y_i + \frac{1}{6}(a + 2b + 2c + d) \Delta t$
これが古典的なRK4です。当然ながら、異なる重みとステップを選択することで、他のルンゲ＝クッタ法の巨大なファミリーを作成できます。

メインのRK4と並行して実行される3次法については、単純に次を選択しました。
$y_{i+1} = y_i + \frac{1}{6}(a + 2b + 3c) \Delta t$

しかし、RK4は1階微分方程式を解くため、2階微分方程式をどのように解くのでしょうか？問題ありません！$p'' = f(\tau, p, p')$ の形の任意の方程式は、$p'$を独立変数と考えることで、追加の次元を追加して1階微分方程式に還元できます。

$ y = \begin{bmatrix} p \ p' \end{bmatrix} $

すると、$ y' = \begin{bmatrix} p \ p' \end{bmatrix}' = \begin{bmatrix} p' \ p'' \end{bmatrix} = \begin{bmatrix} p' \ f(\tau,p,p') \end{bmatrix} =g(\tau,y) $

この方程式は現在1階なので、解くことができます。この方法で使用するためにd2p関数のラッパーを作成しましょう。

void rk4_diff(vec4 y[2], out vec4 dy[2]) {
    dy[0] = y[1];
    dy[1] = d2p(y[0], y[1]);
}

このラッパーは、点の位置と速度を含む8次元ベクトルを受け取り、その導関数を返します。位置の導関数としての速度、およびd2pで計算された「加速度」です。結果のRK34積分器は最も洗練されていませんが、ここにあります。

[!example]

• rk34_step

/**
 * 適応型ルンゲ＝クッタ法（3次および4次）、移動する4次元点の2階微分を計算するために調整されています。
 * [in] epsilon - 誤差推定値の目標値。ステップはこの値に基づいて更新されます。
 * [in] max_step - ステップが取りうる最大のハードリミット。
 * [in] dtau - 積分ステップ。
 * [inout] p0 - 点の位置。
 * [inout] p1 - 点の速度。
 * [out] error_estimate - 推定誤差値（RK3とRK4の結果の差）。
 *
 * 戻り値：
 * true - 誤差は大きすぎず、処理を続行できます。
 * false - 誤差が大きすぎるため、更新されたdtauで再実行する必要があります。
 */

bool rk34_step(
    float epsilon, float max_step, inout float dtau, inout vec4 p0, inout vec4 p1, out float error_estimate)
{
    vec4 a[2], b[2], c[2], d[2];

    vec4 y0[2] = {p0, p1};
    rk4_diff(y0, a);

    vec4 y1[2] = {p0 + a[0] * dtau/2, p1 + a[1] * dtau/2};
    rk4_diff(y1, b);

    vec4 y2[2] = {p0 + b[0] * dtau/2, p1 + b[1] * dtau/2};
    rk4_diff(y2, c);

    vec4 y3[2] = {p0 + c[0] * dtau, p1 + c[1] * dtau};
    rk4_diff(y3, d);

    vec4 rk3_0 = (a[0] + 2 * b[0] + 3 * c[0]) * dtau / 6;
    vec4 rk3_1 = (a[1] + 2 * b[1] + 3 * c[1]) * dtau / 6;

    vec4 rk4_0 = (a[0] + 2 * b[0] + 2 * c[0] + d[0]) * dtau / 6;
    vec4 rk4_1 = (a[1] + 2 * b[1] + 2 * c[1] + d[1]) * dtau / 6;

    vec4 diff0 = rk3_0 - rk4_0;
    vec4 diff1 = rk3_1 - rk4_1;
    error_estimate = sqrt(abs(dot(diff0, diff0) + dot(diff1, diff1)));

    float new_dtau = dtau*sqrt(epsilon/error_estimate);
    if (new_dtau < dtau / 2) {
        dtau = new_dtau;
        return false;
    }
    dtau = new_dtau;
    if (dtau > max_step) {
        dtau = max_step;
    }
    p0 += rk4_0;
    p1 += rk4_1;

    return true;
}

座標

［画像: sphere_with_back.png|事象の地平面の近くにある球殻］

これで、点の位置と速度を入力としてrk34_step積分器に渡す方法がわかりました。これは測地線に沿って小さなステップを取り、誤差を範囲内に保つために次のステップのサイズを調整します。これで、光線をどこにでも送信するために必要なすべてのツールが揃いました。しかし、どこに実際に送信すべきかを知るためのツールはありません。そしてそのためには、以前と同様に、2つの座標系が必要です。計算用の球面座標系と、カメラを扱うためのデカルト座標系です。そして、それらを切り替える方法が必要です。

点座標の変換

球面座標系からデカルト座標系への変換方法はすでにわかっています。これをGLSLコードに変換しましょう。

vec4 sh_to_cart(vec4 sh) {
    float t = sh[0];
    float r = sh[1];
    float th = sh[2];
    float phi = sh[3];

    float x = r*sin(th)*sin(phi);
    float y = r*sin(th)*cos(phi);
    float z = r*cos(th);

    vec4 ret = {x, y, z, t};
    return ret;
}

デカルト座標では、時間は最後のものとして選択されるため、通常のGLSL 4元ベクトル成分のように$x, y, z$にアクセスできます。point.wを介して$t$にアクセスできます。光線をどこに送信するかを決定するために、逆変換も必要になります。デカルト座標（画面上のピクセル）から球面座標への変換です。これはSympyが苦労する部分です。しかし、人間にとっては pretty easy geometry problem です。

vec4 cart_to_sh(vec4 pos) {
    float r = sqrt(pos.x*pos.x + pos.y*pos.y + pos.z*pos.z);
    float th = PI / 2. - atan(sqrt(pos.x*pos.x + pos.y*pos.y), pos.z);
    float phi = atan(pos.y, pos.x);
    float t = pos.w;

    vec4 ret = {t, r, th, phi};
    return ret;
}

atan関数に2つの引数を渡していることに注意してください。GLSLでは、これは三角形の辺が与えられたときの角度を見つけるためのショートカットです（多くの他の言語のatan2関数と同じです）。Sympyはこれを実際には使用できず、通常のタンジェントを使用すると符号に関する情報が失われるため、機械は逆変換を理解できません。

ベクトル座標の変換

繰り返しますが、測地線（つまり光線）は、自身に沿って輸送されるベクトルの経路です。しかし、曲線座標では、同じベクトルが異なる点から始まると、座標線自体が変化するため、異なる成分を持つことになります。したがって、各点には、その接空間内の異なる座標系間を移動する一意の方法があります。この変換はヤコビ行列によって捉えられます。そして再び、Sympyが私たちを助けてくれます。

J = coord_sh.jacobian(coord_cart)
IJ = simplify(J.inv())
print(f"J_{{spherical \\rightarrow cart}} = {latex(J)}")
print(f"J_{{cart \\rightarrow spherical}} = {latex(IJ)}")

$$ J
_{spherical \rightarrow cart} = \left\lbrack\begin{matrix}0 & \sin{\left(\phi \right)} \sin{\left(\theta \right)} & \rho \sin{\left(\phi \right)} \cos{\left(\theta \right)} & \rho \sin{\left(\theta \right)} \cos{\left(\phi \right)}\0 & \sin{\left(\theta \right)} \cos{\left(\phi \right)} & \rho \cos{\left(\phi \right)} \cos{\left(\theta \right)} & - \rho \sin{\left(\phi \right)} \sin{\left(\theta \right)}\0 & \cos{\left(\theta \right)} & - \rho \sin{\left(\theta \right)} & 0\1 & 0 & 0 & 0\end{matrix}\right\rbrack $$

$$ J
_{cart \rightarrow spherical} = \left\lbrack\begin{matrix}0 & 0 & 0 & 1\\sin{\left(\phi \right)} \sin{\left(\theta \right)} & \sin{\left(\theta \right)} \cos{\left(\phi \right)} & \cos{\left(\theta \right)} & 0\\frac{\sin{\left(\phi \right)} \cos{\left(\theta \right)}}{\rho} & \frac{\cos{\left(\phi \right)} \cos{\left(\theta \right)}}{\rho} & - \frac{\sin{\left(\theta \right)}}{\rho} & 0\\frac{\cos{\left(\phi \right)}}{\rho \sin{\left(\theta \right)}} & - \frac{\sin{\left(\phi \right)}}{\rho \sin{\left(\theta \right)}} & 0 & 0\end{matrix}\right\rbrack $$

これにより、球面座標で定義された任意の点に対して、ベクトル座標がどのように変換されるかを決定する行列が存在します。ただし、$ heta = 0$の点ではゼロ除算が発生するため、例外です。しかし、詳細については心配しないでおきましょう。GLSLで使用するために、これらの行列は転置する必要があります。なぜなら、それらは列優先形式で格納されるからです。

// sh_pos (点座標はシュワルツシルト座標で与えられる) におけるヤコビ行列の値
// cartesian_vector = J * sh_vector
mat4 jacobian_inner_to_cart_at(vec4 sh_pos) {
    float rho = sh_pos[1];
    float theta = sh_pos[2];
    float phi = sh_pos[3];

    // OpenGLは行列を列優先形式で格納します。したがって、これは転置されたヤコビ行列です。
    mat4 ret = {
        {0, 0, 0, 1},
        {sin(phi)*sin(theta), sin(theta)*cos(phi), cos(theta), 0},
        {rho*sin(phi)*cos(theta), rho*cos(phi)*cos(theta), -rho*sin(theta), 0},
        {rho*sin(theta)*cos(phi), -rho*sin(phi)*sin(theta), 0, 0},
    };
    return ret;
}

// sh_pos (点座標はシュワルツシルト座標で与えられる) における逆ヤコビ行列の値
// sh_vector = IJ * cartesian_vector
mat4 jacobian_cart_to_inner_at(vec4 sh_pos) {
    float rho = sh_pos[1];
    float theta = sh_pos[2];
    float phi = sh_pos[3];

    // OpenGLは行列を列優先形式で格納します。したがって、これは転置されたヤコビ行列です。
    mat4 ret = {
        {0, sin(theta)*sin(phi), sin(phi)*cos(theta)/rho, cos(phi)/(rho*sin(theta))},
        {0, sin(theta)*cos(phi), cos(phi)*cos(theta)/rho, -sin(phi)/(rho*sin(theta))},
        {0, cos(theta), -sin(theta)/rho, 0},
        {1, 0, 0, 0},
    };
    return ret;
}

フクロウの残りの部分を描画する

［画像: outer_shell.png|事象の地平面のすぐ上にある球殻。片側から見ると、その表面全体が見えます。］

フラグメントシェーダーは、フレームバッファの各ピクセルに対して実行され、そのピクセルの色を出力します。

layout(location = 0) out vec4 diffuseColor;

#define PI 3.1415926535897932384626433832795

uniform vec2 window_size;
uniform vec4 camera_pos; // in Schwarzschild coordinates
uniform mat3 camera_rotation; // in Cartesian coordinates
uniform float viewport_dist = 1.0;
uniform float rs = 0.25;
uniform float shell_radius = 2.0;

入力は、次のuniform変数として提供されます。

• フレームバッファのサイズ
• カメラパラメータ：球面座標で定義された位置、デカルト座標でのカメラ回転行列、およびビューポート距離
• シュワルツシルト半径
• 球殻の半径

ピクセル自体の座標は、gl_FragCoord変数で利用できます。残りは pretty straightforward です。

• 各ピクセルを通して、過去へレイを送信する
• ループで、rk34_step関数を使用して測地線に沿って移動させる
• 球殻に当たったら、ピクセルをその色で塗りつぶす。レイが球殻から離れる場合は、終了する。

さあ、始めましょう！

[!example]
void main()

void main()
{
    // (-1, 1) の範囲での現在のピクセルの位置
    vec2 window_center = {window_size.x / 2, window_size.y / 2};
    vec2 pixel_pos = (vec2(gl_FragCoord) - window_center) / max(window_center.x, window_center.y);

    // 光線。
    // ビューポートからスクリーンピクセルへの方向：
    vec3 direction3 = {viewport_dist, pixel_pos.x, pixel_pos.y};
    direction3 = normalize(direction3);
    // カメラの向きに合わせた調整：
    direction3 = camera_rotation * direction3;
    vec4 direction4 = {direction3.x, direction3.y, direction3.z, -length(direction3)};

    // 現在の位置と4元速度（シュワルツシルト座標）：
    vec4 pos_s = camera_pos;
    vec4 dp = jacobian_cart_to_inner_at(pos_s) * direction4;

    // 初期ステップ
    float dtau = 0.01;

    vec3 color = {0,0,0};
    float prev_pho = pos_s[1];
    float max_step = 0.25;
    int i;
    for (i = 0; i < 2000; i++) {
        float r = pos_s[1];

        if (abs(r - shell_radius) < 1) {
            // 球殻の近くで精度を上げる
            max_step = 0.01 + abs(r - shell_radius)/20.;
        } else {
            max_step = 0.5;
        }

        float error_estimate;
        if (rk34_step(0.01, max_step, dtau, pos_s, dp, error_estimate) == false) {
            // このイテレーションをステップサイズの調整に費やす。
            continue;
        }
        if ((r > shell_radius * 2) && (r - prev_pho > 0)) {
            // 外側から見ている場合の最適化。
            // このレイは球殻を完全に外し、離れていくだけです。
            break;
        }

        // 最適化。これほど近い場合、もう離れることはありません。
        if (r < rs*1.01) {
            break;
        }

        bool entering_shell = (prev_pho >= shell_radius) && (r < shell_radius);
        bool exiting_shell = (prev_pho < shell_radius) && (r >= shell_radius);
        if ( entering_shell || exiting_shell) {
            float lat = pos_s[2]/PI*180; // 度単位
            float lon = pos_s[3]/PI*180; // 度単位
            float line_freq = 15; // 格子線間の度数

            // グリッド線
            float shine = 1.0;
            float thicc = 0.5;

            float prev_lat_line = floor(lat/line_freq)*line_freq;
            float next_lat_line = prev_lat_line + line_freq;
            color.b += shine * exp(-pow((lat - prev_lat_line)/thicc, 2));
            color.b += shine * exp(-pow((lat - next_lat_line)/thicc, 2));

            float prev_lon_line = floor(lon/line_freq)*line_freq;
            float next_lon_line = prev_lon_line + line_freq;
            color.b += shine * exp(-pow((lon - prev_lon_line)/thicc, 2));
            color.b += shine * exp(-pow((lon - next_lon_line)/thicc, 2));

            color.g = color.b/2;

            if (exiting_shell)
                break;
        }
        prev_pho = r;
    }

    diffuseColor = vec4(color, 1.0);
}

結果：

［画像: sphere_1.png|］

このシェーダーは実際には最適ではありません。ゼロに近い値での除算は、精度を維持するために積分器に非常に小さなステップを選択させるため、極の近くでパフォーマンスを低下させます。この邪悪な赤い輝きを見てください！これは、各ピクセルのレンダリングに費やされたイテレーション数に比例します。

［画像: sphere_with_iteration_count.png|］

そして、これが積分器からの誤差推定値です。

［画像: sphere_with_max_error.png|］

この記事の他の画像（トーラスを除く）は、球殻のサイズと色を調整するだけで作成されました。これを数倍高速化する余地があると思います。しかし、それは読者の演習として残しておきます。

<div class=custom-page-footer><hr/> <div class=spaced-lr> #SoME4 June 2025 </div> </div>