C# 中的黑洞光线追踪

• 从 GitHub 下载源代码

注意

本文使用 MathJax 引擎渲染数学公式。

引言

那么，黑洞看起来是什么样的？

一个朴素的答案是：嗯，由于黑洞的极端引力无法逃脱任何光线，它必然是完全漆黑的。

这是真的，但关于光线在黑洞周围如何运动以及可以观察到什么效应的问题，比这有趣得多。这些效应可能非常壮观。

电影《星际穿越》极具启发性，因为它包含了对假设黑洞的非凡渲染，以及围绕它旋转的发光物质盘（称为“吸积盘”）。根据电影的一些“制作过程”材料，这可能是历史上第一次如此精确地生成了天体物理学家预测的黑洞渲染。

在普通 PC 上可以重现其中一些效果吗？在一定程度上可以，使用相当简单、标准的计算机图形技术、相当复杂的数学、鲁棒的数值分析方法和并行计算概念。

本文提供的代码是 [1] 中介绍的方法的扩展版本，增加了纹理映射、相机路径计算和基本的并行计算机制。

背景

为了理解本文，对以下概念有所了解至关重要。

光线追踪

光线追踪是一种标准计算机图形技术，用于模拟计算机内存中三维场景的图像。本质上，它涉及对从观察者射入场景的假想“光线”进行数学模拟。光线追踪算法模拟光线照射到场景中的物体时的行为，并考虑表面着色、反射和折射。它非常适合模拟各种物理现象，并且能够渲染具有精美真实感的图像。

黑洞理论

维基百科上有关于黑洞的详尽文章。在本文中，了解光线在黑洞附近传播的非直观行为——称为引力透镜效应——很有趣。引力会弯曲光线，因此当我们观察一个带有吸积盘的黑洞时，我们实际上同时看到了吸积盘的正面以及从上方和下方观察到的吸积盘的背面！

可视化技术需要计算图像每个像素的颜色，这些颜色来自被从观察者投射到被观察物体上的光线击中的物体的属性。在“普通”光线追踪中，这是（相对）简单的，因为光线是直线传播的。然而，在黑洞这种极端引力源的附近，光线不再沿直线传播。因此，有必要模拟到达观察者的光子的轨迹。

克尔度规

本文描述的代码基于爱因斯坦广义相对论场方程推导出的数学。该模拟基于称为克尔度规的解，该度规描述了围绕未带电、球对称且旋转的黑洞的时空几何。

光子的轨迹从以 Boyer-Lindquist 坐标表示的线元方程推导而来。

其中

\(r\),\(\theta\) 和 \(\phi\) 是物理学中使用的球坐标。

（来源：维基百科）

和

\(\Delta = r^2-2Mr+a^2\)
\(\Sigma = r^2 + a^2\cos^2 \theta\)

\(a = \) 角动量因子 - 非零，因为我们的黑洞正在旋转。

上述方程可用于构建一个常微分方程 (ODE) 系统，该系统描述了光子在黑洞周围场中传播的路径。

有了计算光子轨迹的能力，就可以为光线追踪算法计算光线。光线方程用球坐标 (\(r\), \(\theta\), \(\phi\)) 表示而不是常用的笛卡尔坐标 (\(x\), \(y\), \(z\)) 并不重要——因为需要的是光线击中点的颜色，而不是其坐标。

龙格-库塔法

为了模拟光子的轨迹，有必要采用数值分析。该代码使用 Cash-Karp 方法计算 ODE 系统的数值解。Cash-Karp 方法属于通用的龙格-库塔数值积分方法系列。

程序

实现克尔度规以确定光子轨迹

下面引用了用于计算光子轨迹的数学公式，仅供说明。理解它并非必要才能理解渲染引擎的原理。这些方程只是公式，当由通用的数值分析引擎处理时，它们能够实现光线相交检测，从而确定像素的颜色。

为简化计算，假设许多变量为“归一化”，即假设 M（黑洞质量）和 E（光子能量）为 1，光子质量假设为 0。L（轨道角动量）和 K（卡特常数）需要为每条光线计算。

以下八个微分方程系统描述了光子在克尔黑洞附近（如 [1] 中所示）的路径：

请注意，在代码中，所有系数都被取反。这是因为我们“反向”模拟光子的位置——从观察者回溯到光子的起源。

现在可以忽略 \(t'\)、\(p_\phi\) 和 \(p_t'\)，因为它们不影响光子轨迹的形状。最终，使用龙格-库塔算法对五个 ODE 进行积分。

实现这些方程的代码位于

public class KerrBlackHoleEquation : IODESystem
{
	...
	/// <summary>
	/// Function that returns the equations that are included in the 
        /// coupled differential equation set.
	/// Coupled differential equations describing motion of photon.
	/// </summary>
	/// <param name="y">Vector describing current state of the ODE system</param>
	/// <param name="dydx">Coefficient vector of the differential equations</param>
	public unsafe void Function(double* y, double* dydx) {...}
	...
}

初始条件

为了完成数值积分，每条光线都需要一组初始条件，即描述计算初始状态的 ODE 系统系数的值。这对应于观察点处的模拟状态，因此初始条件需要考虑光线的“目标”。换句话说，有必要能够将光线“导向”图像像素空间的选定坐标，这必须反映在以球坐标表示的 ODE 系统的初始状态值中。因此，光线计算的输入将是（\(r_0\), \(\theta_0\), \(\phi_0\)）——观察者位置和（\(x_0\), \(y_0\)）——像素空间中的坐标。像素坐标会转换为 \(p_r'\) 和 \(p_\theta'\) 的初始值。

其中 \(\alpha\) 是孔径系数。

此外，还需要为光线计算（如 [1] 所示）计算两个常数（L 和 K）。

以上已在此实现。

public class KerrBlackHoleEquation : IODESystem
{
	...
        /// <summary>
        /// Set initial conditions for a starting point of the ray that is being simulated.
        /// </summary>
        /// <param name="y0">Vector of coefficients describing state of the system</param>
        /// <param name="ydot0">Vector of ODE coefficients that is initialized by this function call</param>
        /// <param name="x">x-coordinate of the ray</param>
        /// <param name="y">y-coordinate of the ray</param>
        public unsafe void SetInitialConditions(double* y0, double* ydot0, double x, double y) {...}
	...
}

此外，还将计算黑洞事件视界的半径：

吸积盘

黑洞周围的吸积盘从任意选定的半径延伸到最内层稳定轨道，该轨道也根据黑洞参数计算得出（如 [4] 所示）。

其中

光线路径示例

下图说明了程序计算的光线路径示例。这些路径是通过绘制几条穿过由黑洞中心（奇点）和观察者位置确定的垂直平面上的光线生成的。球坐标已重新计算为笛卡尔坐标，并沿着 XZ 轴（x - 水平，z - 垂直）绘制了图形。

这里光线击中了吸积盘，因此它们结束在坐标 z = 0 上。

请注意光线是如何弯曲的——这意味着黑洞的观察者实际上将看到来自场景周围宇宙的光线。这对外星宇宙的纹理提出了特殊的要求——纹理必须是“全方位的”！

光线相交检测

随着光线路径的计算，光线追踪器需要确定光线是否未击中物体，从而结束计算。可能出现三种情况：

1. 光线击中黑洞——当 \(r < r_{hor}\) 时发生。
2. 光线逃逸到无穷远——实际上，当 \(r > r_0\) 时假设如此。
3. 光线击中吸积盘——当 \(\theta = \pm {\pi \over 2}\) 且 \(r_s < r < r_{disc}\) 时发生。

为了检测光线是否击中水平面，龙格-库塔步长会在计算确定光线“穿过”水平面 \((\theta = \pm{\pi \over 2})\) 时，通过使用收缩的正负步长值来交替进行。这个过程类似于“二分查找”交点。

驱动光线追踪算法的逻辑实现在 RayTracer 类中。

/// <summary>
/// Ray tracer class that is responsible for all image calculations.
/// </summary>
public class RayTracer
{
	...
}

图像质量

光线追踪器类能够以可变质量渲染图像。有以下几种质量级别可用：

• Low - 计算目标图像的每 4x4 像素块一条光线。（适用于测试相机路径等。）
• Medium - 对图像相邻的 4x4 像素块执行自适应抗锯齿——如果两个相邻 4x4 像素块之间的颜色差异大于阈值，则追踪每个像素一条光线。
• High - 对相邻像素执行自适应抗锯齿——如果两个相邻像素之间的颜色差异大于阈值，则对像素进行超采样。
• UltraHigh - 在任何条件下都进行超采样——每个像素 16 条光线。

纹理映射

场景中的纹理有两种需求——吸积盘和场景周围的宇宙。

宇宙通过将高分辨率星空纹理映射到一个围绕场景的球体来建模，以黑洞的奇点为中心，以观察者位置的 R0 为半径。因此，算法后续步骤中积分的任何具有 r 坐标大于 R0 的光线，都被假定为“逃逸到无穷远”，其当前坐标根据以下公式映射到分辨率为 \((x_{res}, y_{res})\) 的矩形位图中：

请注意，球坐标如何使映射变得简单，并且 \(r\) 对于此映射无关紧要！
高质量、高分辨率的星空纹理来自此处。

吸积盘的纹理也需要高分辨率。当光线追踪器检测到光线与盘面相交时，会计算映射。映射到矩形位图如下：

相机路径

该程序使用光线追踪类来生成动画的连续帧。为了使动画有趣，每一帧都需要不同，即场景设置需要遵循预定义的路径。场景设置通过 SceneDescription 类反映。

/// <summary>
/// Scene description class
/// </summary>
public class SceneDescription
{
	/// <summary>
	/// Camera position - Distance from black hole
	/// </summary>
	public double ViewDistance { get; set; }

	/// <summary>
	/// Camera position - Inclination (vertical angle) in degrees
	/// </summary>
	public double ViewInclination { get; set; }

	/// <summary>
	/// Camera position - Angle (horizontal) in degrees
	/// </summary>
	public double ViewAngle { get; set; }

	/// <summary>
	/// Camera tilt - in degrees
	/// </summary>
	public double CameraTilt { get; set; }

	/// <summary>
	/// Camera aperture - need to manipulate the camera angle.
	/// </summary>
	public double CameraAperture { get; set; }

	/// <summary>
	/// Camera yaw - if we want to look sideways.
	/// Note: this is expressed in % of image width.
	/// </summary>
	public double CameraYaw { get; set; }
}

需要有一个参数化函数，该函数为动画的给定时间点计算场景参数。这是 DefaultPathSceneGenerator 类的职责。

public class DefaultPathSceneGenerator : ISceneGenerator
{
	public SceneDescription GetScene(int frame, double fps)
	{
		SceneDescription result = new SceneDescription();
		result.CameraAperture = 2.0;

		double t = frame / fps;

		double r = 600 - t * 2.53;

		// factor of attenuation of camera's sinusoidal motions (the closer to black hole - the calmer the flight is)
		double calmFactor = Math.Pow((600 - r) / 575, 20); 
		
		double phi = t*3;
		double theta = 84
			+ 8 * Math.Sin(phi * Math.PI / 180) * (1 - calmFactor) // precession
			+ 3 * calmFactor;

		result.ViewAngle = phi;
		result.ViewDistance = r;
		result.ViewInclination = theta;
		result.CameraAperture = 24.00/500.0*r + 3.2;
		result.CameraTilt = 8.0 * Math.Cos(phi * Math.PI / 180) * (1 - calmFactor);
		result.CameraYaw =  calmFactor * 1.0; // we will be 'landing' on the accretion disc...

		return result;
	}
}

上面编码的相机路径描述了一个螺旋，始于（r = 600）单位，并围绕黑洞缓慢旋转，在旅程结束时“坠落”速度迅速加快。在路径结束时，相机略微转动，似乎“降落”在吸积盘上，以便更好地观察。相机路径还会定期倾斜相机，以模拟轨道略微偏离吸积盘平面的情况。整个旅程大约需要 4 分 30 秒，然后相机撞上事件视界。

渲染控制文件

为了定义长批次的计算（例如，用于多帧动画），程序使用一个控制文件。控制文件的后续行具有以下含义：

• 指定一个类作为每一帧的场景定义源。
• 输出路径，帧图像将存储在该路径下。
• 帧分辨率。
• FPS - 每秒帧数。
• 质量
• 帧行：格式为 {帧号} {帧状态：1 - 待处理，2 - 正在处理，3 - 已完成}
• ...

控制文件的示例如下：

GraviRayTraceSharp.Scene.DefaultPathSceneGenerator, GraviRayTraceSharp
.
640 480
25
High
1 1
2 1
3 1
...

分布式计算

程序的 Main() 方法接受一个参数：要使用的控制文件的名称。该文件以独占模式打开，如果失败，则在特定时间内重试。这意味着，控制文件成为了在多台机器上分布式计算的一个简单同步机制，例如在局域网中（如果控制文件和输出路径通过 UNC 路径发布）。

使用此机制，本文的代码已被用于计算一系列数千帧 1600x1200 的图像（约 4 分钟的 25 FPS 动画）。使用了多达 20 台 Intel I7 级 PC 的伪集群。总计算时间约为 200 小时。

未实现

本文提出的渲染引擎使用一个非常简单的吸积盘模型——假定吸积盘的厚度为 0，其外观通过纹理模拟。没有模拟红移效应。

关注点

该项目为探索 IT 和物理的多学科方面提供了有趣的机会。

• 计算机图形学，光线追踪。
• 数值分析——ODE 计算的龙格-库塔法。
• 分布式、并行计算。

……最终还附带了相当有趣的图像 :-)。

参考文献

• [1] 曲时空光线追踪
• [2] 用 Python 编写的黑洞光线追踪器
• [3] 克尔度规中的轨道（Sean Marshall，康奈尔大学）
• [4] 黑洞吸积盘理论基础（Marek Abramowicz，P. Chris Fragile）
• [5] 伪复数广义相对论中的光线追踪