Spacetime Raytracer

Overview

Real-time simulation of light-bending around black holes or wormholes, with relativity-accurate movements of the spaceship. Camera rotations are modelled using quaternions.

Technical highlights

• Runs at 60 fps at 800×800 resolution on the CPU.
• Templated math library (also using concepts) written from scratch.
• Adaptive Runge-Kutta-Fehlberg 45 solver written from scratch.
• Spacetimes templated on lambdas describing the Lorentzian metric.
• Parallelization via OpenMP.
• GUI written using the ImGui library.
• CMake build system.

Implementation

General

The project is built using the CMake build system. All of the third-party dependencies are fetched automatically at the CMake configuration step. High resolution textures that can be used in the simulation are also automatically downloaded provided that the FETCH_TEXTURES flag is set. For the release build, the -ffast-math flag is enabled (visual inspection with and without the flag showed no visible loss of accuracy; the same is true for using double-precision types rather than single-precision).

Structure

The project almost entirely consists of header-only modules, with a single main.cpp file that constructs the App and runs it. The modules are Image, Math, Physics, GUI, CLI, each wrapped in its own namespace.

• The Image namespace provides some basic implementations of contiguous pixel buffers, textures, and allows for loading/saving images from the filesystem.
• The Math namespace is a generic math library implementing tensors (and hence vectors/matrices), quaternions, as well as the adaptive integrators used in the simulation. It is built with performance in mind, so the “hot” functions and structures only use stack allocated variables. E.g. tensors are class templates with compile-time rank and dimensions, so index mismatches are caught at compile time, storage is fixed-size on the stack, and small loops can be unrolled by the compiler.
• The Physics namespace uses the Math and Image libraries, and contains the bulk of the computational details of the simulation.
• The GUI namespace sets up the interface using ImGui and controls the simulation generated by the render engine from the Physics library at runtime.
• The CLI namespace features a simple parser that allows the user to control the render engine directly from the terminal.

Raytracing optimization

Each frame of the simulation is generated pixel-by-pixel. For each pixel, the direction of the incoming light ray is determined based on the projection used for the camera. The light ray is then traced backwards to its “origin” by integrating (using an adaptive RKF45 scheme) the geodesic equation for the selected spacetime, which in particular describes how light bends due to the curvature of spacetime.

Note that the origin of the light ray is either:

• a point on the skybox, if the backward trace reaches some fixed large radius away from the center of the simulated body,
• a point on the body, if the backward trace reaches the physical radius of the simulated body.

The color of the pixel is then determined by mapping the origin of the light ray to the selected texture for the universe/body.

A naive per-pixel geodesic solve runs too slowly on the CPU, even when parallelized. Profiling with perf revealed that the CPU time is mainly dominated by floating-point operations, so the main goal is to reduce them. Two details are key for achieving this:

• Spherically symmetric reduction. All geodesics in spherically symmetric spacetimes are planar. The engine thus rotates the coordinate system in a way that forces the geodesic to stay on the equatorial plane, solves the equation constrained to this plane, and then rotates it back to get the correct solution. This is beneficial because the equatorial geodesic equation takes a highly reduced and simplified form (4-dimensional rather than 8-dimensional, with approximately 50 times fewer FP operations per Runge-Kutta step).

• One ODE solve per-angle instead of per-pixel. In that plane, a light ray’s trajectory depends only on one angle $\alpha \in [0, π]$ that parametrizes its initial direction. Each frame precomputes the end-state for $16 \max(W, H)$ angles into a lookup table, and the pixel loop does only a nearest lookup and a texture fetch. This turns $O(WH)$ integrations per frame into $O(\max(W, H))$: at 800×800, 12,800 solves instead of 640,000. Both the construction of the lookup table and the pixel loop are parallelized with dynamic scheduling. The lookup table depends on the position of the observer so it is recomputed at each frame. However, it only depends on the radius of the position, so one could also precompute an array of lookup tables by discretizing some reasonable radius interval, which would eliminate the need to run the integrator at runtime, at the cost of using a large amount of memory. The lookup table can be further optimized to reduce the number of necessary $\text{atan2}$ computations (work in progress).