广义相对论的可视化

作者：米凯洛·莫罗兹

当考虑曲速驱动器和黑洞等事物的渲染时，我们通常只期望看到一个简单的近似值或艺术家的再现，假设实现准确的数学结果需要至少拥有数学物理学博士学位的人。我不会说这是完全不真实的，但在这篇博文中，我将尝试解释一种在 100 行左右的代码内进行相当准确的可视化的方法，基本上适用于您可以实现的任何类型的时空将其指标写为代码。尽管如此，这种方法的详细数学推导可能有点数学繁重。

任何 GR 渲染的主要成分都是弄清楚光线如何移动。了解光线如何移动，我们可以追踪从相机进入场景的光线并查看光线来自哪里。因此，为了渲染没有对象的最简单场景，我们只需为每个像素跟踪一条光线，并将像素的颜色分配给光线最终指向的方向上的天空盒的颜色。

什么是测地线？
最短路径的数学描述
测地线的拉格朗日描述
测地线的哈密顿描述
用 GLSL 编写哈密顿测地线示踪器
结论
参考
什么是测地线？
我们到底如何在曲线空间中追踪光线？弯曲空间内的任何物体都遵循所谓的测地线。

从某种意义上说，测地线只是一个奇特的词，指空间内两点之间的最短长度的路径。

我应该注意到，可能有多个这样的路径，这些路径在局部是最小的，从某种意义上说，您无法微调路径以使其更短，而在全局范围内可能存在更短的路径。另外，在闵可夫斯基时空中，由于虚距离的原因（当 \( ds^2 < 0 \) 时），定义有点复杂。ds2< 0）。

在我们的例子中，我们感兴趣的不是两点之间的路径，而是在给定初始点和方向的情况下找出射线如何移动，但在推导描述测地线的方程时，上面的定义仍然有用，我们将在此处使用该方程。

最短路径的数学描述
我将尝试快速进行推导。如果您想跳过数学部分，请跳转到用 GLSL 编写哈密顿测地线示踪器。

从数学上来说，我们有一些坐标系、一条路径和一种计算两点之间距离的方法。

坐标系是一组标记空间中每个点的数字。路径是一个接受路径参数并输出坐标的函数。在广义相对论中，路径参数通常是原时间（就像随着物体移动的时钟，标记每个点），但它实际上可以是任何东西。计算距离的方法称为度量（它是这里可怕的数学的主要来源）。

在物理学或更普遍的微分几何中，度量被定义为称为度量张量的积分（“和”）。度量张量是双线性形式 \( g(a, b) \)，它本质上将向量对映射到实数，并且是弯曲空间点积的推广。因此，使用度量张量，我们可以找到空间中向量的长度，以及空间中无限接近的点之间的距离 \( ds \) 。克（一、二），它本质上是将向量对映射到实数，并且是弯曲空间的点积的推广。因此，使用度量张量，我们可以找到空间中向量的长度以及距离ds空间中无限接近的点之间。

ds2= g( d, d​x )
在我们的例子中，我们将向量描述为一组数字，度量只是某个矩阵 \( g_{\mu \nu} \) 乘以向量的矩阵乘积。对于我们的无穷小距离 \( ds \)，我们得到这个表达式：克μν​乘以向量。为了我们的无限小距离ds我们得到这个表达式：

ds2=Σμν​氮克μν​dxμdxν
通常，总和只是通过爱因斯坦符号[1] 隐式假设的。

ds2=克μν​dxμdxν
在这里我们实际上可以看到，对于 \( g_{\mu \nu} \) 的一些简单选择，我们可以通过毕达哥拉斯定理得到距离。特别针对度量张量矩阵是单位矩阵的情况。克μν​我们可以通过毕达哥拉斯定理得到距离。特别针对度量张量矩阵是单位矩阵的情况。

ds2= dx21+ dx22+ dx23
对于像空间这样的平坦时空，我们得到类似的结果，但时间坐标分量带有负号。

ds2= − dx20+ dx21+ dx22+ dx23
这里我使用了 \( (- + + +) \) 签名，但实际上可以在不改变测地线的情况下翻转符号，并且在某些情况下，例如对于粒子物理学，使用相反的 \( (+ - - -） \） 签名。( − + + + )签名，但实际上可以在不改变测地线的情况下翻转符号，并且在某些情况下，例如粒子物理学，使用相反的符号更有意义( + --- )​​签名。

回到计算距离的主要问题，要计算 A 和 B 两个点之间沿某条路径 \( x^i(t) \) 的长度，我们只需使用积分将无穷小距离加在一起：x我( t )我们只需要使用积分将无穷小的距离加在一起：

l =∫乙一个克μν​dxμdxν--------√=∫乙一个克μν​dxμdxν--------√dtdt=∫乙一个克μν​dxμdtdxνdt----------√dt
其中 \( \frac{dx^i}{dt} \) 只是坐标 x 相对于路径参数（“时钟”）变化的速度，在某种意义上可以解释为速度。dx我dt只是坐标 x 相对于路径参数（“时钟”）变化的速度，在某种意义上可以解释为速度。

现在我们的主要问题是如何最小化路径长度？在这里，我们引入一种称为变分法的东西，粗略地说，这是一种通过其输入函数（路径）的无限小的变化来找出函数（距离）如何变化的方法。此类导数具有与普通函数导数类似的性质。事实上，与微积分类似，要找到函数的极值（最小值、最大值或驻点），我们只需将变化等于 0。

测地线的拉格朗日描述
物理学中有一个完整的分支与变分原理相关，它指出任何类型的物理系统都有一个它喜欢最小化的值（或者更一般地在路径的微小变化下保持不变，即“静止”）。该值称为作用，积分下的函数称为系统的拉格朗日函数。研究系统拉格朗日函数的物理学分支称为拉格朗日力学。

在我们的例子中，拉格朗日量可以写成这样：

L =克μν​dxμdtdxνdt----------√
事实证明，我们不需要函数最小值的平方根即可成为测地线，我们可以简单地将其从测地线拉格朗日中删除：

L =克μν​dxμdtdxνdt
您可以在此处找到证明[2]。这种简化造成的唯一区别是路径的参数化可能不同，但路径本身是相同的。另外值得注意的是，这种特定情况的方程是相同的，这意味着参数化也是相同的。

此外，我们希望使用 1/2 因子来简化方程。

L =12克μν​dxμdtdxνdt
所以，我们现在的目标是最小化这个函数：

S=∫乙一个12克μν​dxμdtdxνdtdt
一般来说，这样的函数的最小值可以通过应用欧拉-拉格朗日方程[3] 找到：

∂L∂x我-ddt∂L∂dx我dt= 0
欧拉-拉格朗日方程推导（点击展开）
您可以在这里找到更详细的推导[4]

让我们导出测地拉格朗日的欧拉-拉格朗日方程（请记住，每个坐标都有一个方程 \( x^i \) ）：x我）：

∂L∂dx我dt=12∂∂dx我dt克μν​dxμdtdxνdt=12克ν​dxνdt+12克μi​dxμdt=克ν​dxνdt
然后我们对路径参数求导：

ddt（克ν​dxνdt） =d克ν​dtdxνdt+克ν​d2xνdt2=d克ν​dxμdxμdtdxνdt+克ν​d2xνdt2
最后：

∂L∂x我=12d克μν​dx我dxμdtdxνdt
将 \eqref{el1} 和 \eqref{el2} 代入欧拉-拉格朗日方程 \eqref{el} 得出测地线方程：（？？？）和（？？？）代入欧拉-拉格朗日方程（？？？）引出了测地线方程：

12d克μν​dx我dxμdtdxνdt-d克ν​dxμdxμdtdxνdt-克ν​d2xνdt2= 0
乘以度量张量逆 \( - g^{i \nu} \) 我们得到：-克ν​我们得到：

d2x我dt2+克ν​（d克ν​dxμ-12d克μν​dx我）dxμdtdxνdt= 0
这就是我们最终的测地线方程组，我们当然也可以在这里替换 Christoffel 符号，但对于我们的应用来说没有区别。当然，我们可以只使用这些方程来追踪测地射线，然后就到此为止了，但不幸的是，这需要计算大量的导数（在 4d 时空中，具体来说有 64 个导数），无论是手动还是通过使用数值微分。值得庆幸的是，有一种方法可以避免这种情况，并且实际上简化了整个算法！ （以轻微的性能成本为代价）

测地线的哈密顿描述
这场秀的明星来了——哈密顿力学。哈密​​顿运动方程有一个非常好的形式，可以轻松编写一个计算机程序，通过使用欧拉积分对它们进行积分。

dp我dt= -∂H∂x我

\begin{方程}
dx我dt=∂H∂p我
其中 \( p \) 是所谓的广义动量，它是拉格朗日关于坐标路径参数导数的导数：p是所谓的广义动量，它是拉格朗日关于坐标路径参数导数的导数：

p我=∂L∂dx我dt
要获得哈密顿量本身，您需要对拉格朗日应用勒让德变换[6]：

H=Σ我氮p我dx我dt- L
哈密​​顿运动方程推导（点击展开）
汉密尔顿运动方程的不同推导可以在此处找到[5]。

对于我们的情况，动量如下，我们在写下欧拉-拉格朗日方程 \eqref{el0} 时已经计算出，并给出了广义动量 \eqref{momentumdef} 的定义：（？？？），并给出广义动量的定义（？？？）:

p我=克我jdxjdt
要获得“时间”导数，您只需将两边乘以度量张量逆：

dx我dt=克我jpj
以及哈密顿量本身：

H=Σ我氮dx我dtpi − L = g​我jdx我dtdxjdt-12克我jdx我dtdxjdt=12克我jdx我dtdxjdt= L
事实证明，对于测地拉格朗日的这个简单选择，哈密顿量等于拉格朗日！

另外，我们想通过将 \eqref{dxdt} 代入哈密顿方程来了解哈密顿量作为广义动量的函数：（？？？）代入哈密顿方程：

H=12克我jdx我dtdxjdt=12克我jp我pj
虽然运动方程简单地是：

dp我dt= -∂H∂x我
dx我dt=克我jpj
这就是我们编写数值测地积分器所需的全部内容！

用 GLSL 编写哈密顿测地线示踪器
你可能注意到了，在最后的汉密尔顿运动方程中我没有写出 \(\frac{\partial H}{\partial x^i} \)，这其实很重要！我们希望保持哈密顿量的导数不变，因为这样我们就不需要计算度量张量的 64 个导数，而只需要 4 个即可找到哈密顿量梯度。这是测地线追踪算法的主要简化。∂H∂x我，这个其实很重要！我们希望保持哈密顿量的导数不变，因为这样我们就不需要计算度量张量的 64 个导数，而只需要 4 个即可找到哈密顿量梯度。这是测地线追踪算法的主要简化。

这里我们将使用 GLSL 着色语言，因为它的变量和函数可以很好地映射到我们将在这里执行的数学运算。最重要的是，我们可以轻松制作实时 GR 可视化着色器。

首先，我们需要一个函数来计算空间和时间中 4d 点的度量张量。让我们使用Alcubierre 曲速驱动[7] 度量作为示例，因为它非常简单。

mat4 Metric ( vec4 x )
{
  //Alcubierre 度量
  const  float R =  1.0 ;
  常量 浮点西格玛=  35.0 ；
  常量 浮点v =  1.1 ；

  浮动x = v * x 。x ;
  float r =  sqrt ( sqr ( x . y  - x )  + x . z * x . z  + x . w * x . w ) ;
  浮点数f =  0.5 * ( tanh ( sigma * ( r + R ) )  - tanh( sigma * ( r - R ) ) ) / tanh ( sigma * R ) ;
  浮动gtt = v * v * f * f -  1.0 ;
  浮动gxt =  - v * f ;

  返回 mat4 ( gtt , gxt ,   0 ,   0 ,
              gxt ,    1 ,   0 ,   0 ,
                0 ,    0 ,   1 ,   0 ,
                0 ,    0 ,   0 ,   1 ) ;
}
在我们的例子中，x 是表示位置的 4D 向量。第一个组成部分x.x是x[0]时间。作为输出，我们得到一个 4 x 4 矩阵，用mat4GLSL 表示。

然后我们需要写下哈密顿量 \eqref{hamiltonian}。哈密​​顿量是一个函数，它接受两个因素：时空位置和 4d 动量向量，并输出一个标量。（？？？）。哈密​​顿量是一个函数，它接受两个因素：时空位置和 4d 动量向量，并输出一个标量。

浮点哈密顿量( vec4 x ,  vec4 p )
{
  mat4 g_inv =逆( Metric ( x ) ) ;
  返回 0.5 *点( g_inv * p , p ) ;
}
作为奖励，这里是拉格朗日

浮点拉格朗日( vec4 x ,  vec4 dxdt )
{
  返回 0.5 *点(公制( x ) * dxdt , dxdt ) ;
}
令人惊讶的是，GLSL 已经有了一个矩阵逆函数inverse()，在其之上，哈密顿量只是 和 的点积（在 GLSL 意义上）g_inv*p，p它们分别是逆变动量向量和协变动量向量。逆变动量实际上就是坐标的时间导数dxdt，即dot(dxdt,p)。

之后我们需要计算哈密顿量的 4D 梯度。我们可以通过在所有 4 个空间方向上使用前向数值差（使用一些小值）来做到这一点eps：

vec4 HamiltonianGradient ( vec4 x ,  vec4 p )
{
  const  float eps =  0.001 ;
  返回 ( vec4 (哈密顿量( x +  vec4 ( eps , 0 , 0 , 0 ) , p ) ,
               哈密顿量( x +  vec4 ( 0 , eps , 0 , 0 ) , p) ,
               哈密顿量( x +  vec4 ( 0 , 0 , eps , 0 ) , p ) ,
               哈密顿量( x +  vec4 ( 0 , 0 , 0 , eps ) , p ) )  -哈密顿量( x , p ) ) / eps ;
}
现在我们有了哈密顿梯度，我们终于可以写出运动方程 \eqref{eqmotion1} \eqref{eqmotion2} 积分代码（？？？） （？？？）集成代码

vec4 IntegrationStep ( inout  vec4 x ,  inout  vec4 p )
{
  const  float TimeStep =  0.1 ;
  p = p -时间步长*哈密顿梯度( x , p ) ;
  x = x +时间步长* inverse ( Metric ( x ) )  * p ;
}
您可能会问，“等等，就这些了？”，事实上这就是积分测地线所需的全部。当然，它是相当慢的，因为我们进行了多达 6 个矩阵逆评估，通过用没有逆矩阵的拉格朗日函数替换大多数哈密顿量，可以将其优化为 1，因为它们是相等的。更好的是已经通过分析计算了度量逆，但是不可能对每个度量，特别是对于隐式定义的度量。

当然还有最后一个问题，虽然初始化时空位置很容易，但是p在开始追踪时我们如何初始化动量向量的值呢？

在追踪测地线之前，您可以使用方程 \eqref{momentum}（？？？）
p =公制( x )  * dxdt ；
但 dxdt 是什么？它只不过是光线在时空中移动的 4D 方向。方向可分为 3 类：

类时间，当 \( A < 0 \)A < 0
空，当 \( A = 0 \)A = 0
类空间，当 \( A > 0 \)A > 0
其中 \(A\) 是 \begin{方程}一个是
A =克μν​dxμdtdxνdt
或者在 GLSL 中，A = dot(Metric(x) * dxdt, dxdt)

在模拟光的传播方式时，我们只需要零方向，这会导致零测地线解。另一方面，如果您想对运动速度慢于光速的物体进行建模，则需要类似时间的测地线。还有快子物质的类空间测地线，这在现实生活中不会发生，所以我们忽略它。

因此，假设平坦空间度量 \eqref{flat} 和空间中的某个 3D 方向，我们的光线 p 将为（？？？）空间中的某个 3D 方向我们的光线 p 是

vec4 GetNullMomentum ( vec3 dir )
{
  return Metric ( x )  *  vec4 ( 1.0 , 标准化( dir ) ) ;
}
以及这个操作的逆操作

vec3 GetDirection ( vec4 p )
{
  vec4 dxdt = inverse ( Metric ( x ) )  * p ;
  返回 标准化( dxdt . yzw ) ;
}
所以，最终的简单算法将如下所示：

void TraceGeodesic ( inout  vec3 pos ,  inout  vec3 dir ,  inout  float time )
{
  vec4 x =  vec4 (时间,位置) ;
  vec4 p = GetNullMomentum ( dir ) ;

  常量 int步骤=  256 ；
  for ( int i =  0 ; i <步骤; i ++ )
  {
    IntegrationStep ( x , p ) ;
    ​​ //例如，当 x 低于事件视界时，您可以在此处添加停止条件
  }

  位置= x 。yzw ;
  时间= x 。x ;
  dir = GetDirection ( p ) ;
}
本质上，这只是一个 4D 射线行进算法，其中射线的方向每一步都会改变。在这种特定情况下，步长的大小也会发生变化，这可以通过标准化动量来避免p = normalize(p)。这只会改变步长，不会改变测地线路径，即，它的工作原理就像路径的动态重新参数化一样。积分的时间步长也可以根据所使用的度量而变化。例如，在黑洞的情况下，我根据到事件视界的距离成比例地更改时间步长，以便测地线的精度大致与空间曲率成正比。这是获得准确结果的重要优化，同时保持计算成本相对较小。

示例着色器代码：

mat4 diag ( vec4 a )
{
    return  mat4 ( a . x , 0 , 0 , 0 ,
                0 , a . y , 0 , 0 ,
                0 , 0 , a . z , 0 ,
                0 , 0 , 0 , a . w ) ;
}

mat4 Metric ( vec4 x )
{
    //Kerr-Schild 坐标中的 Kerr-Newman 度量
    const  float a =  0.8 ;
    常量 浮点数m =  1.0 ；
    常量 浮点Q =  0.0 ；
    vec3 p = x 。yzw ;
    浮点rho = 点( p , p )  - a * a ;
    浮动r2 =  0.5 *( rho +  sqrt ( rho * rho +  4.0 * a * a * p .z * p .z ) ) ;​​浮点数r = sqrt ( r2 ) ; vec4 k = vec4 ( 1 , ( r * p . x + a * p . y ) / ( r2 + a
     
       * a ) ,  ( r * p . y  - a * p . x ) / ( r2 + a * a ) , p . z / r ) ;
    浮点数f = r2 * ( 2.0 * m * r - Q * Q ) / ( r2 * r2 + a * a *p 。z * p 。）；​
    返回f * mat4 ( k . x * k , k . y * k , k . z * k , k . w * k ) + diag ( vec4 ( - 1 , 1 , 1 , 1 ) ) ;
}

浮点哈密顿量( vec4 x ,  vec4 p )
{
    mat4 g_inv =逆( Metric ( x ) ) ;
    返回 0.5 *点( g_inv * p , p ) ;
}

/*
浮点拉格朗日(vec4 x, vec4 dxdt)
{
    返回 0.5*dot(公制(x)*dxdt,dxdt);
}
*/

vec4 HamiltonianGradient ( vec4 x ,  vec4 p )
{
    const  float eps =  0.001 ;
    返回 ( vec4 (哈密顿量( x +  vec4 ( eps , 0 , 0 , 0 ) , p ) ,
                 哈密顿量( x +  vec4 ( 0 , eps , 0 , 0 ) ,p ) ,
                 哈密顿量( x +  vec4 ( 0 , 0 , eps , 0 ) , p ) ,
                 哈密顿量( x +  vec4 ( 0 , 0 , 0 , eps ) , p ) )  -哈密顿量( x , p ) ) / eps ;
}

void IntegrationStep ( inout  vec4 x ,  inout  vec4 p )
{
    const  float TimeStep =  0.15 ;
    p = p -时间步长*哈密顿梯度( x , p ) ;
    x = x +时间步长* inverse ( Metric ( x ) )  * p ;
}

vec4 GetNullMomentum ( vec4 x ,  vec3 dir )
{
    return Metric ( x )  *  vec4 ( 1.0 , 标准化( dir ) ) ;
}

vec3 GetDirection ( vec4 x ,  vec4 p )
{
    vec4 dxdt = inverse ( Metric ( x ) )  * p ;
    返回 标准化( dxdt . yzw ) ;
}

void TraceGeodesic ( inout  vec3 pos ,  inout  vec3 dir ,  inout  float time )
{
    vec4 x =  vec4 (时间,位置) ;
    vec4 p = GetNullMomentum ( x , dir ) ;

    常量 int步骤=  256 ；
    for ( int i =  0 ; i <步骤; i ++ )
    {
        IntegrationStep ( x , p ) ;
        ​​ //例如，当 x 低于事件视界时，您可以在此处添加停止条件
    }

    位置= x 。yzw ;
    时间= x 。x ;
    dir = GetDirection ( x , p ) ;
}

void mainImage ( out  vec4 fragColor ,  in  vec2 fragCoord )  {
    vec2 uv =  2.0  *  ( fragCoord -  0.5  * iResolution . xy )  /  max ( iResolution . x , iResolution . y ) ;

    vec3 RayPos =  vec3 (  0.0,0.0,32.0 )   ;​​  ​​ vec3 RayDir = vec3 ( uv.x , uv.y , -1.0 ) ;​​​​​浮动  时间= 0.0 ;
      
     

    RayDir = 归一化( RayDir ) ;

    TraceGeodesic （ RayPos 、 RayDir 、时间）；

    fragColor =  vec4 (纹理( iChannel0 , RayDir ) .rgb , 1.0 ) ;​ }
您可以查看此 Shadertoy 实现来查看一些优化，例如可变时间步长、用拉格朗日量替换哈密顿量、使用对称矩阵求逆函数（更快一点）、重用一些计算值（如果 Shadertoy 为黑色则重新启动）：


上面的着色器在 Kerr-Schild 坐标[8]（本质上是笛卡尔坐标）中实现了 Alcubierre 度量和Kerr-Newman 度量。

mat4 diag ( vec4 a )
{
    return  mat4 ( a . x , 0 , 0 , 0 ,
                0 , a . y , 0 , 0 ,
                0 , 0 , a . z , 0 ,
                0 , 0 , 0 , a . w ) ;
}

mat4 Metric ( vec4 x )
{
  //Kerr-Schild 坐标中的 Kerr-Newman 度量
  const  float a =  0.8 ;
  常量 浮点数m =  1.0 ；
  常量 浮点Q =  0.0 ；
  vec3 p = x 。yzw ;
  浮点rho = 点( p , p )  - a * a ;
  浮动r2 =  0.5 *( rho +  sqrt ( rho * rho +  4.0 * a * a * p .z * p .z ) ) ;​​浮点数r = sqrt ( r2 ) ; vec4 k = vec4 ( 1 , ( r * p . x + a * p . y ) / ( r2 + a
   
     * a ) ,  ( r * p . y  - a * p . x ) / ( r2 + a * a ) , p . z / r ) ;
  浮点数f = r2 * ( 2.0 * m * r - Q * Q ) / ( r2 * r2 + a * a *p 。z * p 。）；​
  返回f * mat4 ( k . x * k , k . y * k , k . z * k , k . w * k ) + diag ( vec4 ( - 1 , 1 , 1 , 1 ) ) ;     
}
当然，这不是优化的限制。其他主要优化是计算度量的解析逆。对于一大类指标，您可以使用谢尔曼-莫里森公式[9]

一些需要记住的有用的事情
由于我们使用数值有限差分，因此结果实际上在很大程度上取决于浮点数的相对值。例如，远离坐标系中心，数值导数的精度会低很多，因此您需要改变 的大小eps以避免过多的数值噪声。此外，指标通常具有应该避免的数值奇点，除非您想获得 NaN 结果。

我通常避免使用球坐标中的度量，因为它们的极轴奇异性，它具有强烈的视觉效果，即使时间步长变化很小，也很难避免，尽管此类度量通常在数学上更简单，并且允许更大的时间步长而不破坏外观黑洞。对于球对称度量，例如非旋转黑洞和虫洞，有一个技巧可以完全避免极奇点！关于球对称的事情是，测地线总是在 2d 平面内移动，可以映射到坐标系的赤道平面，基本上将 3d + 时间问题简化为 2d + 时间。 （斯科特·曼利有一段视频解释了他如何渲染虫洞，在那里他使用了这个技巧+预先计算一个查找表来大大简化计算）

我还在我的 Shaderoy 虫洞中使用了降维技巧：


还有一种不同的方法可用于以数值方式计算导数，并且更加准确。本质上它是基于对偶数的前向自动微分，有一些Shadertoy示例使用了这种方法。

最后，您始终可以通过分析方式推导方程，虽然这是最烦人的方法，但它通常是性能最快的方法。一个折衷的解决方案是自动导出方程，这种方法被James Berrow制作的geodesic_raytracing使用（你应该在 Twitter 上关注他，他在这个主题上有很多很酷的东西）。

确定光线是否落入事件视界内实际上并不简单，并且没有通用的方法，虽然如果光线位于事件视界表面下方，您可以将颜色设置为 0，但从以下方向查看事物时，这是不正确的黑洞内部。追踪光线也应该及时向后完成，因为我们追踪来自相机的光线，而不是追踪到相机的光线，这对最终的渲染有明显的影响，如果不这样做，也会导致黑洞内部完全黑暗的渲染，即使光确实存在于视界之下，并且可以从外部到达。

广义相对论中的红移可以简单地根据物体发射/吸收速度的点积乘以光子动量的比率来计算，其中光子动量平行于光子运动方向。动量需要沿着测地线从发射到吸收平行传输，好消息是我们可以使用测地线 4 速度来代替，因为它在技术上也是沿着测地线平行传输并指向运动方向。 （以防万一，这也意味着我们在积分方程时需要避免重新规范化广义动量）

另外值得注意的是，广义相对论中的点积是简单地根据度量张量 \( g(u,v) \) 定义的克（紫外线）​​
u ⋅ v =克μν​你μvν
最后，如果您只想计算史瓦西黑洞的测地线，您可以简单地使用某个经典力场中质量为 1 的粒子的方程，详细信息请参见
这篇博文。