wyj/GR-raytracing
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

init

Browse Source
This commit is contained in:
wyj committed 2026-08-25 22:45:43 -04:00
commit d30f9ac56c
56 files changed
+40776
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

src/catalog.c
+107
View File
@@ -0,0 +1,107 @@
#include "catalog.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define PI 3.14159265358979323846
enum {
TEST_GRID_LINE_DEG = 10,
TEST_GRID_SAMPLE_DEG = 2,
TEST_GRID_RED_TEMPERATURE_K = 3000,
TEST_GRID_BLUE_TEMPERATURE_K = 12000,
};
/*
* With the renderer's 380--780 nm CIE integration and linear-sRGB luminance,
* B(3000 K) / B(12000 K) = 0.00141095580387. Red stars retain unit scale.
*/
#define TEST_GRID_BLUE_AMPLITUDE 0.00141095580387
static int test_grid_temperature_K(int longitude_deg, int latitude_deg)
{
/* Boundaries belong to the octant immediately east/north of them. */
const int longitude_sector = longitude_deg / 90;
const int hemisphere = latitude_deg < 0 ? 0 : 1;
/* Add the hemisphere bit to flip the color across the equator. */
const int octant = hemisphere + longitude_sector;
return octant % 2 == 0 ? TEST_GRID_RED_TEMPERATURE_K
: TEST_GRID_BLUE_TEMPERATURE_K;
}
static double test_grid_amplitude(int temperature_K)
{
return temperature_K == TEST_GRID_BLUE_TEMPERATURE_K
? TEST_GRID_BLUE_AMPLITUDE
: 1.0;
}
int catalog_write_octant_grid(const char *path)
{
FILE *file = fopen(path, "w");
if (file == NULL) return -1;
fputs("longitude_deg,latitude_deg,temperature_K,amplitude\n", file);
for (int latitude = -90; latitude <= 90; latitude += TEST_GRID_SAMPLE_DEG) {
for (int longitude = 0; longitude < 360;
longitude += TEST_GRID_SAMPLE_DEG) {
const int on_longitude_line = longitude % TEST_GRID_LINE_DEG == 0;
const int on_latitude_line = latitude % TEST_GRID_LINE_DEG == 0;
if ((!on_longitude_line && !on_latitude_line) ||
((latitude == -90 || latitude == 90) && longitude != 0))
continue;
const int temperature_K =
test_grid_temperature_K(longitude, latitude);
fprintf(file, "%d,%d,%d,%.12g\n", longitude, latitude,
temperature_K, test_grid_amplitude(temperature_K));
}
}
return fclose(file) == 0 ? 0 : -1;
}
int catalog_load_csv(StarCatalog *catalog, const char *path)
{
FILE *file = fopen(path, "r");
char line[256];
size_t capacity = 0;
catalog->stars = NULL;
catalog->count = 0;
if (file == NULL || fgets(line, sizeof line, file) == NULL) goto fail;
while (fgets(line, sizeof line, file) != NULL) {
double longitude, latitude, temperature, amplitude;
if (sscanf(line, "%lf,%lf,%lf,%lf", &longitude, &latitude,
&temperature, &amplitude) != 4)
goto fail;
if (catalog->count == capacity) {
size_t next = capacity == 0 ? 256 : capacity * 2;
Star *stars = realloc(catalog->stars, next * sizeof *stars);
if (stars == NULL) goto fail;
catalog->stars = stars;
capacity = next;
}
const double lon = longitude * PI / 180.0;
const double lat = latitude * PI / 180.0;
const double cos_lat = cos(lat);
Star *star = &catalog->stars[catalog->count++];
star->direction[0] = cos_lat * cos(lon);
star->direction[1] = sin(lat);
star->direction[2] = cos_lat * sin(lon);
star->temperature_K = temperature;
star->amplitude = amplitude;
}
fclose(file);
return 0;
fail:
if (file != NULL) fclose(file);
catalog_destroy(catalog);
return -1;
}
void catalog_destroy(StarCatalog *catalog)
{
free(catalog->stars);
catalog->stars = NULL;
catalog->count = 0;
}
src/catalog.h
+27
View File
@@ -0,0 +1,27 @@
#ifndef CATALOG_H
#define CATALOG_H
#include <stddef.h>
typedef struct {
double direction[3];
double temperature_K;
double amplitude;
} Star;
typedef struct {
Star *stars;
size_t count;
} StarCatalog;
/*
* Synthetic lensing fixture: stars lie on the union of 10-degree longitude
* and latitude lines, sampled every 2 degrees. The eight longitude/hemisphere
* octants alternate between red and blue blackbody temperatures. Blue stars'
* scale compensates for their greater visible blackbody radiance.
*/
int catalog_write_octant_grid(const char *path);
int catalog_load_csv(StarCatalog *catalog, const char *path);
void catalog_destroy(StarCatalog *catalog);
#endif
src/frame.c
+378
View File
@@ -0,0 +1,378 @@
#include "frame.h"
#include "optics.h"
#include <limits.h>
#include <math.h>
#include <omp.h>
#include <stdint.h>
#include <stdlib.h>
/* Private HDR buffers make splatting independent across workers. This is an
* allocation cap, not a rendering parameter: callers transparently fall back
* to the serial implementation when the frame is too large for two buffers. */
#define FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES ((size_t)512 * 1024 * 1024)
static const double pi = 3.14159265358979323846;
static double dot(const double a[3], const double b[3]) {
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
static void cross(const double a[3], const double b[3], double out[3]) {
out[0] = a[1] * b[2] - a[2] * b[1];
out[1] = a[2] * b[0] - a[0] * b[2];
out[2] = a[0] * b[1] - a[1] * b[0];
}
static double normalize(double vector[3]) {
const double length = sqrt(dot(vector, vector));
if (length > 0.0)
for (int i = 0; i < 3; ++i)
vector[i] /= length;
return length;
}
static size_t vertex_index(int column, int row, int columns) {
return (size_t)row * (columns + 1) + column;
}
int frame_lens_mesh_build_coarse(FrameLensMesh *mesh, int width, int height,
int cell_pixels, double horizontal_fov_deg) {
if (mesh == NULL || width <= 0 || height <= 0 || cell_pixels <= 0 ||
horizontal_fov_deg <= 0.0 || horizontal_fov_deg >= 179.0)
return -1;
const int columns = (width + cell_pixels - 1) / cell_pixels;
const int rows = (height + cell_pixels - 1) / cell_pixels;
const size_t vertex_count = (size_t)(columns + 1) * (rows + 1);
const size_t triangle_count = (size_t)columns * rows * 2;
LensVertex *vertices = calloc(vertex_count, sizeof *vertices);
LensTriangle *triangles = malloc(triangle_count * sizeof *triangles);
if (vertices == NULL || triangles == NULL) {
free(vertices);
free(triangles);
return -1;
}
const double tan_half_x = tan(horizontal_fov_deg * pi / 360.0);
const double tan_half_y = tan_half_x * (double)height / width;
for (int row = 0; row <= rows; ++row) {
const double image_y = (double)row * height / rows;
for (int column = 0; column <= columns; ++column) {
LensVertex *vertex = &vertices[vertex_index(column, row, columns)];
vertex->image_x = (double)column * width / columns;
vertex->image_y = image_y;
vertex->camera_direction[0] = 1.0;
vertex->camera_direction[1] = (0.5 - image_y / height) * 2.0 * tan_half_y;
vertex->camera_direction[2] =
(vertex->image_x / width - 0.5) * 2.0 * tan_half_x;
normalize(vertex->camera_direction);
}
}
size_t next_triangle = 0;
for (int row = 0; row < rows; ++row)
for (int column = 0; column < columns; ++column) {
const size_t top_left = vertex_index(column, row, columns);
const size_t top_right = vertex_index(column + 1, row, columns);
const size_t bottom_left = vertex_index(column, row + 1, columns);
const size_t bottom_right = vertex_index(column + 1, row + 1, columns);
triangles[next_triangle++] =
(LensTriangle){{top_left, bottom_left, bottom_right}};
triangles[next_triangle++] =
(LensTriangle){{top_left, bottom_right, top_right}};
}
*mesh = (FrameLensMesh){.vertices = vertices,
.triangles = triangles,
.vertex_count = vertex_count,
.triangle_count = triangle_count};
return 0;
}
int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime,
const ObserverState *observer,
const GeodesicTraceConfig *trace) {
if (mesh == NULL || spacetime == NULL || observer == NULL || trace == NULL)
return -1;
/* Each iteration exclusively owns one vertex. SpacetimeSource is shared
* read-only here; backends with mutable evaluation state must keep it
* thread-local. */
#pragma omp parallel for schedule(static)
for (size_t i = 0; i < mesh->vertex_count; ++i) {
LensVertex *vertex = &mesh->vertices[i];
RayEndpoint endpoint = geodesic_trace_past(spacetime, observer,
vertex->camera_direction, trace);
vertex->status = endpoint.status;
if (endpoint.status == RAY_ENDPOINT_ESCAPED) {
for (int axis = 0; axis < 3; ++axis)
vertex->n_infinity[axis] = endpoint.n_infinity[axis];
vertex->log_frequency_ratio = log(endpoint.frequency_ratio);
}
}
return 0;
}
static double spherical_area(const double a[3], const double b[3],
const double c[3]) {
double b_cross_c[3];
cross(b, c, b_cross_c);
return 2.0 * atan2(fabs(dot(a, b_cross_c)),
1.0 + dot(a, b) + dot(b, c) + dot(c, a));
}
static int spherical_barycentric(const double point[3], const double a[3],
const double b[3], const double c[3],
double weights[3]) {
const double area = spherical_area(a, b, c);
double edge_cross[3];
if (area < 1e-14)
return -1;
const double *corners[3] = {a, b, c};
for (int edge = 0; edge < 3; ++edge) {
const double *left = corners[edge];
const double *right = corners[(edge + 1) % 3];
const double *opposite = corners[(edge + 2) % 3];
cross(left, right, edge_cross);
/* This is a sign test, so its tolerance must scale with the source
* triangle. A fixed absolute threshold turns sufficiently fine triangles
* into near-all-sky queries. */
if (dot(edge_cross, point) * dot(edge_cross, opposite) <
-1e-14 * dot(edge_cross, edge_cross))
return -1;
}
weights[0] = spherical_area(point, b, c) / area;
weights[1] = spherical_area(point, c, a) / area;
weights[2] = spherical_area(point, a, b) / area;
return 0;
}
static int usable_triangle(const FrameLensMesh *mesh,
const LensTriangle *triangle,
const LensVertex *vertices[3]) {
for (int i = 0; i < 3; ++i) {
vertices[i] = &mesh->vertices[triangle->vertex[i]];
if (vertices[i]->status != RAY_ENDPOINT_ESCAPED)
return 0;
}
return spherical_area(vertices[0]->n_infinity, vertices[1]->n_infinity,
vertices[2]->n_infinity) >= 1e-14;
}
/* Give a source lying exactly on a shared source edge to one triangle only.
* Interior overlaps remain valid separate lens images. */
static int owns_source_boundary(const LensTriangle *triangle,
const double weights[3]) {
for (int opposite = 0; opposite < 3; ++opposite) {
if (weights[opposite] > 1e-11)
continue;
const size_t left = triangle->vertex[(opposite + 1) % 3];
const size_t right = triangle->vertex[(opposite + 2) % 3];
if (left > right)
return 0;
}
return 1;
}
static size_t splat_catalog_triangles(const FrameLensMesh *mesh,
const StarCatalog *catalog, double *hdr,
int width, int height, double exposure,
const PointSpreadFunction *psf,
size_t first_triangle,
size_t last_triangle) {
size_t images = 0;
for (size_t t = first_triangle; t < last_triangle; ++t) {
const LensVertex *vertex[3];
if (!usable_triangle(mesh, &mesh->triangles[t], vertex))
continue;
const double source_area = spherical_area(
vertex[0]->n_infinity, vertex[1]->n_infinity, vertex[2]->n_infinity);
const double image_area =
spherical_area(vertex[0]->camera_direction, vertex[1]->camera_direction,
vertex[2]->camera_direction);
const double magnification = image_area / source_area;
for (size_t s = 0; s < catalog->count; ++s) {
const Star *star = &catalog->stars[s];
double weights[3];
if (spherical_barycentric(star->direction, vertex[0]->n_infinity,
vertex[1]->n_infinity, vertex[2]->n_infinity,
weights))
continue;
if (!owns_source_boundary(&mesh->triangles[t], weights))
continue;
const double image_x = weights[0] * vertex[0]->image_x +
weights[1] * vertex[1]->image_x +
weights[2] * vertex[2]->image_x;
const double image_y = weights[0] * vertex[0]->image_y +
weights[1] * vertex[1]->image_y +
weights[2] * vertex[2]->image_y;
const double log_g = weights[0] * vertex[0]->log_frequency_ratio +
weights[1] * vertex[1]->log_frequency_ratio +
weights[2] * vertex[2]->log_frequency_ratio;
const LinearRgb color =
blackbody_to_linear_rgb(star->temperature_K * exp(log_g));
splat_moffat(hdr, width, height, image_x, image_y, color,
exposure * star->amplitude * magnification, psf);
++images;
}
}
return images;
}
size_t frame_splat_catalog(const FrameLensMesh *mesh,
const StarCatalog *catalog, double *hdr, int width,
int height, double exposure,
const PointSpreadFunction *psf) {
if (mesh == NULL || catalog == NULL || hdr == NULL || exposure <= 0.0 ||
psf == NULL || width <= 0 || height <= 0)
return 0;
const size_t pixel_count = (size_t)width * height * 3;
if (pixel_count > SIZE_MAX / sizeof(double) ||
pixel_count * sizeof(double) > FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES / 2)
return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure,
psf, 0, mesh->triangle_count);
const size_t buffer_bytes = pixel_count * sizeof(double);
size_t worker_count = FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES / buffer_bytes;
const int max_threads = omp_get_max_threads();
if (worker_count > (size_t)max_threads)
worker_count = (size_t)max_threads;
if (worker_count < 2 || worker_count > INT_MAX)
return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure,
psf, 0, mesh->triangle_count);
double **private_hdr = calloc(worker_count, sizeof *private_hdr);
if (private_hdr == NULL)
return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure,
psf, 0, mesh->triangle_count);
size_t allocated = 0;
for (; allocated < worker_count; ++allocated) {
private_hdr[allocated] = calloc(pixel_count, sizeof **private_hdr);
if (private_hdr[allocated] == NULL)
break;
}
if (allocated != worker_count) {
while (allocated > 0)
free(private_hdr[--allocated]);
free(private_hdr);
return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure,
psf, 0, mesh->triangle_count);
}
size_t images = 0;
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images)
{
const size_t worker = (size_t)omp_get_thread_num();
const size_t first = mesh->triangle_count * worker / worker_count;
const size_t last = mesh->triangle_count * (worker + 1) / worker_count;
images += splat_catalog_triangles(mesh, catalog, private_hdr[worker], width,
height, exposure, psf, first, last);
}
#pragma omp parallel for schedule(static)
for (size_t pixel = 0; pixel < pixel_count; ++pixel)
for (size_t worker = 0; worker < worker_count; ++worker)
hdr[pixel] += private_hdr[worker][pixel];
for (size_t worker = 0; worker < worker_count; ++worker)
free(private_hdr[worker]);
free(private_hdr);
return images;
}
static void blend_gray(double *hdr, int width, int height, int x, int y,
double gray, double alpha) {
if (x < 0 || x >= width || y < 0 || y >= height)
return;
double *pixel = &hdr[3 * (y * width + x)];
for (int channel = 0; channel < 3; ++channel)
pixel[channel] = (1.0 - alpha) * pixel[channel] + alpha * gray;
}
static double fractional_part(double value) { return value - floor(value); }
static void plot_aa(double *hdr, int width, int height, int steep, int x, int y,
double coverage, double gray, double opacity) {
if (coverage > 0.0)
blend_gray(hdr, width, height, steep ? y : x, steep ? x : y, gray,
coverage * opacity);
}
/* Xiaolin Wu line rasterization: a one-pixel line with coverage-based alpha. */
static void draw_line(double *hdr, int width, int height,
const LensVertex *from, const LensVertex *to, double gray,
double opacity) {
double x0 = from->image_x, y0 = from->image_y;
double x1 = to->image_x, y1 = to->image_y;
const int steep = fabs(y1 - y0) > fabs(x1 - x0);
if (steep) {
double swap = x0;
x0 = y0;
y0 = swap;
swap = x1;
x1 = y1;
y1 = swap;
}
if (x0 > x1) {
double swap = x0;
x0 = x1;
x1 = swap;
swap = y0;
y0 = y1;
y1 = swap;
}
const double dx = x1 - x0;
if (dx == 0.0) {
plot_aa(hdr, width, height, steep, (int)lround(x0), (int)floor(y0), 1.0,
gray, opacity);
return;
}
const double gradient = (y1 - y0) / dx;
double x_end = round(x0);
double y_end = y0 + gradient * (x_end - x0);
double x_gap = 1.0 - fractional_part(x0 + 0.5);
int x_pixel_start = (int)x_end;
int y_pixel = (int)floor(y_end);
plot_aa(hdr, width, height, steep, x_pixel_start, y_pixel,
(1.0 - fractional_part(y_end)) * x_gap, gray, opacity);
plot_aa(hdr, width, height, steep, x_pixel_start, y_pixel + 1,
fractional_part(y_end) * x_gap, gray, opacity);
double inter_y = y_end + gradient;
x_end = round(x1);
y_end = y1 + gradient * (x_end - x1);
x_gap = fractional_part(x1 + 0.5);
const int x_pixel_end = (int)x_end;
y_pixel = (int)floor(y_end);
plot_aa(hdr, width, height, steep, x_pixel_end, y_pixel,
(1.0 - fractional_part(y_end)) * x_gap, gray, opacity);
plot_aa(hdr, width, height, steep, x_pixel_end, y_pixel + 1,
fractional_part(y_end) * x_gap, gray, opacity);
for (int x = x_pixel_start + 1; x < x_pixel_end; ++x) {
y_pixel = (int)floor(inter_y);
plot_aa(hdr, width, height, steep, x, y_pixel,
1.0 - fractional_part(inter_y), gray, opacity);
plot_aa(hdr, width, height, steep, x, y_pixel + 1, fractional_part(inter_y),
gray, opacity);
inter_y += gradient;
}
}
void frame_draw_mesh(const FrameLensMesh *mesh, double *hdr, int width,
int height, double gray, double opacity) {
if (mesh == NULL || hdr == NULL || width <= 0 || height <= 0 || gray < 0.0 ||
opacity < 0.0 || opacity > 1.0)
return;
for (size_t i = 0; i < mesh->triangle_count; ++i) {
const LensTriangle *triangle = &mesh->triangles[i];
for (int edge = 0; edge < 3; ++edge) {
const size_t from_id = triangle->vertex[edge];
const size_t to_id = triangle->vertex[(edge + 1) % 3];
if (from_id < to_id)
draw_line(hdr, width, height, &mesh->vertices[from_id],
&mesh->vertices[to_id], gray, opacity);
}
}
}
void frame_lens_mesh_destroy(FrameLensMesh *mesh) {
if (mesh == NULL)
return;
free(mesh->vertices);
free(mesh->triangles);
*mesh = (FrameLensMesh){0};
}
src/frame.h
+46
View File
@@ -0,0 +1,46 @@
#ifndef FRAME_H
#define FRAME_H
#include "catalog.h"
#include "geodesic.h"
#include "optics.h"
#include "observer.h"
#include "spacetime.h"
#include <stddef.h>
typedef struct {
double image_x, image_y;
double camera_direction[3];
double n_infinity[3];
double log_frequency_ratio;
RayEndpointStatus status;
} LensVertex;
typedef struct {
size_t vertex[3];
} LensTriangle;
typedef struct {
LensVertex *vertices;
LensTriangle *triangles;
size_t vertex_count, triangle_count;
} FrameLensMesh;
int frame_lens_mesh_build_coarse(FrameLensMesh *mesh, int width, int height,
int cell_pixels, double horizontal_fov_deg);
int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime,
const ObserverState *observer,
const GeodesicTraceConfig *trace);
/* Each locally invertible escaped triangle contributes one image per contained
* star. */
size_t frame_splat_catalog(const FrameLensMesh *mesh,
const StarCatalog *catalog, double *hdr, int width,
int height, double exposure,
const PointSpreadFunction *psf);
void frame_draw_mesh(const FrameLensMesh *mesh, double *hdr, int width,
int height, double gray, double opacity);
void frame_lens_mesh_destroy(FrameLensMesh *mesh);
#endif
src/geodesic.c
+189
View File
@@ -0,0 +1,189 @@
#include "geodesic.h"
#include <math.h>
typedef struct {
double x[3], Pi[3], log_alpha_p0;
} State;
typedef struct {
double x[3], Pi[3], log_alpha_p0;
} Derivative;
static double dot(const double a[3], const double b[3]) {
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
static int invert(double a[3][3], double b[3][3]) {
double d = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
if (!isfinite(d) || fabs(d) < 1e-14)
return -1;
b[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / d;
b[0][1] = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) / d;
b[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / d;
b[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / d;
b[1][1] = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) / d;
b[1][2] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) / d;
b[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / d;
b[2][1] = (a[0][1] * a[2][0] - a[0][0] * a[2][1]) / d;
b[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) / d;
return 0;
}
/* Equation (4) and (5) of Bohn et al., arXiv:1410.7775. */
static int rhs(const SpacetimeSource *source, double t, const State *s,
Derivative *out) {
MetricData m;
double inv[3][3], up[3] = {0}, da_pi = 0, k_pi_pi = 0;
if (spacetime_eval(source, t, s->x, &m) || m.alpha <= 0 ||
invert(m.gamma, inv))
return -1;
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
up[i] += inv[i][j] * s->Pi[j];
for (int i = 0; i < 3; i++) {
da_pi += m.d_alpha[i] * up[i];
for (int j = 0; j < 3; j++)
k_pi_pi += m.K[i][j] * up[i] * up[j];
}
for (int i = 0; i < 3; i++) {
double db_pi = 0, dg_pi_pi = 0;
out->x[i] = m.alpha * up[i] - m.beta[i];
for (int j = 0; j < 3; j++) {
db_pi += m.d_beta[i][j] * s->Pi[j];
for (int k = 0; k < 3; k++) {
/* Equation (4) requires partial_i gamma^{jk}, while MetricData
* deliberately stores partial_i gamma_jk because that is what metric
* backends interpolate naturally. Differentiate gamma^{-1}:
* partial_i gamma^{jk} = -gamma^{ja}(partial_i gamma_ab)gamma^{bk}.
*/
double d_inverse_gamma = 0.0;
for (int a = 0; a < 3; ++a)
for (int b = 0; b < 3; ++b)
d_inverse_gamma -= inv[j][a] * m.d_gamma[i][a][b] * inv[b][k];
dg_pi_pi += d_inverse_gamma * s->Pi[j] * s->Pi[k];
}
}
out->Pi[i] = -m.d_alpha[i] + (da_pi - m.alpha * k_pi_pi) * s->Pi[i] +
db_pi - 0.5 * m.alpha * dg_pi_pi;
}
out->log_alpha_p0 = -da_pi + m.alpha * k_pi_pi;
return 0;
}
static State add(const State *s, const Derivative *d, double h) {
State r = *s;
for (int i = 0; i < 3; i++) {
r.x[i] += h * d->x[i];
r.Pi[i] += h * d->Pi[i];
}
r.log_alpha_p0 += h * d->log_alpha_p0;
return r;
}
static int rk4(const SpacetimeSource *source, double t, double h, State *s) {
Derivative a, b, c, d;
State q;
if (rhs(source, t, s, &a))
return -1;
q = add(s, &a, h / 2);
if (rhs(source, t + h / 2, &q, &b))
return -1;
q = add(s, &b, h / 2);
if (rhs(source, t + h / 2, &q, &c))
return -1;
q = add(s, &c, h);
if (rhs(source, t + h, &q, &d))
return -1;
for (int i = 0; i < 3; i++) {
s->x[i] += h * (a.x[i] + 2 * b.x[i] + 2 * c.x[i] + d.x[i]) / 6;
s->Pi[i] += h * (a.Pi[i] + 2 * b.Pi[i] + 2 * c.Pi[i] + d.Pi[i]) / 6;
}
s->log_alpha_p0 += h *
(a.log_alpha_p0 + 2 * b.log_alpha_p0 + 2 * c.log_alpha_p0 +
d.log_alpha_p0) /
6;
return 0;
}
static int initialize(const SpacetimeSource *source, const ObserverState *o,
const double n[3], State *s) {
MetricData m;
double k[4] = {o->tetrad[0][0], o->tetrad[0][1], o->tetrad[0][2],
o->tetrad[0][3]};
if (spacetime_eval(source, o->coordinate_time, o->coordinate_position, &m) ||
m.alpha <= 0)
return -1;
for (int a = 0; a < 3; a++)
for (int mu = 0; mu < 4; mu++)
k[mu] -= n[a] * o->tetrad[a + 1][mu];
if (k[0] <= 0)
return -1;
for (int i = 0; i < 3; i++) {
s->x[i] = o->coordinate_position[i];
s->Pi[i] = 0;
for (int j = 0; j < 3; j++)
s->Pi[i] += m.gamma[i][j] * (k[j + 1] + m.beta[j] * k[0]);
s->Pi[i] /= m.alpha * k[0];
}
s->log_alpha_p0 = log(m.alpha * k[0]);
return isfinite(s->log_alpha_p0) ? 0 : -1;
}
static int escaped_direction(const SpacetimeSource *source, double t,
const State *s, double n[3]) {
MetricData m;
double inv[3][3], norm = 0;
if (spacetime_eval(source, t, s->x, &m) || invert(m.gamma, inv))
return -1;
for (int i = 0; i < 3; i++) {
n[i] = 0;
for (int j = 0; j < 3; j++)
n[i] -= inv[i][j] * s->Pi[j];
norm += n[i] * n[i];
}
if (norm <= 0)
return -1;
for (int i = 0; i < 3; i++)
n[i] /= sqrt(norm);
return 0;
}
RayEndpoint geodesic_trace_past(const SpacetimeSource *source,
const ObserverState *observer,
const double n[3],
const GeodesicTraceConfig *config) {
RayEndpoint out = {.frequency_ratio = 0,
.magnification = 1,
.status = RAY_ENDPOINT_INTEGRATION_FAILURE};
State s;
double t;
if (!source || !observer || !config || config->coordinate_time_step <= 0 ||
!config->max_steps || fabs(dot(n, n) - 1) > 1e-10 ||
initialize(source, observer, n, &s))
return out;
t = observer->coordinate_time;
for (unsigned int i = 0; i < config->max_steps; i++) {
if (config->capture_log_alpha_p0 > 0.0 &&
s.log_alpha_p0 >= config->capture_log_alpha_p0) {
out.status = RAY_ENDPOINT_CAPTURED;
return out;
}
SpacetimeRayStatus status = spacetime_classify(source, t, s.x);
if (status != SPACETIME_RAY_ACTIVE) {
out.status = status == SPACETIME_RAY_ESCAPED ? RAY_ENDPOINT_ESCAPED
: RAY_ENDPOINT_CAPTURED;
if (out.status == RAY_ENDPOINT_ESCAPED &&
escaped_direction(source, t, &s, out.n_infinity) == 0)
out.frequency_ratio = exp(-s.log_alpha_p0);
else if (out.status == RAY_ENDPOINT_ESCAPED)
out.status = RAY_ENDPOINT_INTEGRATION_FAILURE;
return out;
}
if (rk4(source, t, -config->coordinate_time_step, &s))
return out;
t -= config->coordinate_time_step;
}
out.status = RAY_ENDPOINT_MAX_STEPS;
return out;
}
src/geodesic.h
+37
View File
@@ -0,0 +1,37 @@
#ifndef GEODESIC_H
#define GEODESIC_H
#include "observer.h"
#include "spacetime.h"
typedef enum {
RAY_ENDPOINT_ESCAPED,
RAY_ENDPOINT_CAPTURED,
RAY_ENDPOINT_MAX_STEPS,
RAY_ENDPOINT_INTEGRATION_FAILURE
} RayEndpointStatus;
typedef struct {
double n_infinity[3];
double frequency_ratio; /* E_camera / E_infinity */
double magnification; /* Filled by the future local inverse lens map. */
RayEndpointStatus status;
} RayEndpoint;
typedef struct {
double coordinate_time_step;
unsigned int max_steps;
/* A positive value terminates a backwards ray whose horizon redshift has
* made log(alpha p^0) reach this value. Zero disables this analytic/demo
* criterion; numerical moving-puncture backends use their AH-calibrated
* spatial cutoff instead. */
double capture_log_alpha_p0;
} GeodesicTraceConfig;
/* camera_direction is a unit vector in the observer's (forward, up, right)
* tetrad. */
RayEndpoint geodesic_trace_past(const SpacetimeSource *source,
const ObserverState *observer,
const double camera_direction[3],
const GeodesicTraceConfig *config);
#endif
src/main.c
+200
View File
@@ -0,0 +1,200 @@
#include "catalog.h"
#include "frame.h"
#include "optics.h"
#include "spacetime.h"
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
typedef struct {
int width, height;
int coarse_cell_pixels;
int draw_mesh;
double horizontal_fov_deg, look_ra_deg, look_dec_deg, exposure;
double observer_inward_speed;
PointSpreadFunction psf;
const char *catalog_path;
const char *output_path;
} Settings;
static int parse_int(const char *text, int *value) {
char *end;
errno = 0;
long parsed = strtol(text, &end, 10);
if (errno || *end || parsed <= 0 || parsed > 16384)
return -1;
*value = (int)parsed;
return 0;
}
static int parse_double(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value <= 0.0 || *value >= 179.0 ? -1 : 0;
}
static int parse_ra_deg(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value < 0.0 || *value >= 360.0 ? -1 : 0;
}
static int parse_dec_deg(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value < -90.0 || *value > 90.0 ? -1 : 0;
}
static int parse_positive(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value <= 0.0 ? -1 : 0;
}
static int parse_speed(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value < 0.0 || *value >= 1.0 ? -1 : 0;
}
static int parse_moffat_beta(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value <= 1.0 ? -1 : 0;
}
static int parse_args(int argc, char **argv, Settings *s,
const char **write_path) {
*s = (Settings){1280,
720,
16,
0,
30.0,
270.0,
0.0,
100.0,
0.0,
{2.7, 4.5},
"assets/sky_grid_5deg.csv",
"output/imgs/minkowski_sky.ppm"};
*write_path = NULL;
for (int i = 1; i < argc; ++i) {
if (!strcmp(argv[i], "--catalog") && i + 1 < argc)
s->catalog_path = argv[++i];
else if (!strcmp(argv[i], "--output") && i + 1 < argc)
s->output_path = argv[++i];
else if (!strcmp(argv[i], "--width") && i + 1 < argc &&
!parse_int(argv[++i], &s->width)) {
} else if (!strcmp(argv[i], "--height") && i + 1 < argc &&
!parse_int(argv[++i], &s->height)) {
} else if (!strcmp(argv[i], "--coarse-cell-pixels") && i + 1 < argc &&
!parse_int(argv[++i], &s->coarse_cell_pixels)) {
} else if (!strcmp(argv[i], "--draw-mesh")) {
s->draw_mesh = 1;
} else if (!strcmp(argv[i], "--fov-deg") && i + 1 < argc &&
!parse_double(argv[++i], &s->horizontal_fov_deg)) {
} else if (!strcmp(argv[i], "--look-ra-deg") && i + 1 < argc &&
!parse_ra_deg(argv[++i], &s->look_ra_deg)) {
} else if (!strcmp(argv[i], "--look-dec-deg") && i + 1 < argc &&
!parse_dec_deg(argv[++i], &s->look_dec_deg)) {
} else if (!strcmp(argv[i], "--exposure") && i + 1 < argc &&
!parse_positive(argv[++i], &s->exposure)) {
} else if (!strcmp(argv[i], "--observer-inward-speed") && i + 1 < argc &&
!parse_speed(argv[++i], &s->observer_inward_speed)) {
} else if (!strcmp(argv[i], "--psf-fwhm-pixels") && i + 1 < argc &&
!parse_positive(argv[++i], &s->psf.fwhm_pixels)) {
} else if (!strcmp(argv[i], "--psf-moffat-beta") && i + 1 < argc &&
!parse_moffat_beta(argv[++i], &s->psf.moffat_beta)) {
} else if (!strcmp(argv[i], "--write-catalog") && i + 1 < argc)
*write_path = argv[++i];
else
return -1;
}
return 0;
}
static int render_frame(const Settings *s, const StarCatalog *catalog,
const SpacetimeSource *spacetime) {
#ifdef SPACETIME_SCHWARZSCHILD
ObserverState observer;
if (observer_inward_schwarzschild_ks(1.0, 30.0,
s->observer_inward_speed, &observer))
return -1;
const GeodesicTraceConfig trace = {.coordinate_time_step = 0.1,
.max_steps = 4096,
.capture_log_alpha_p0 = 8.0};
#else
const ObserverState observer =
observer_fixed_at_origin_look_at(s->look_ra_deg, s->look_dec_deg);
const GeodesicTraceConfig trace = {.coordinate_time_step = 1.0,
.max_steps = 2048};
#endif
FrameLensMesh mesh = {0};
double *hdr = calloc((size_t)s->width * s->height * 3, sizeof *hdr);
if (hdr == NULL ||
frame_lens_mesh_build_coarse(&mesh, s->width, s->height,
s->coarse_cell_pixels,
s->horizontal_fov_deg) ||
frame_lens_mesh_trace(&mesh, spacetime, &observer, &trace)) {
frame_lens_mesh_destroy(&mesh);
free(hdr);
return -1;
}
size_t images = frame_splat_catalog(&mesh, catalog, hdr, s->width, s->height,
s->exposure, &s->psf);
if (s->draw_mesh)
frame_draw_mesh(&mesh, hdr, s->width, s->height, 0.5, 0.5);
int result = write_tonemapped_image(s->output_path, hdr, s->width, s->height);
fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s)\n",
images, catalog->count, s->output_path,
result == 0 ? "ok" : "write failed");
frame_lens_mesh_destroy(&mesh);
free(hdr);
return result;
}
int main(int argc, char **argv) {
Settings settings;
const char *write_path;
if (parse_args(argc, argv, &settings, &write_path)) {
fprintf(stderr,
"Usage: %s [--catalog PATH] [--output PATH] [--width N] [--height "
"N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] "
"[--exposure E] [--observer-inward-speed V] "
"[--psf-fwhm-pixels N] [--psf-moffat-beta N] "
"[--coarse-cell-pixels N] [--draw-mesh] "
"[--write-catalog PATH]\n",
argv[0]);
return 2;
}
if (write_path != NULL)
return catalog_write_octant_grid(write_path) == 0 ? 0
: (perror(write_path), 1);
StarCatalog catalog;
if (catalog_load_csv(&catalog, settings.catalog_path)) {
if (catalog_write_octant_grid(settings.catalog_path) ||
catalog_load_csv(&catalog, settings.catalog_path)) {
perror(settings.catalog_path);
return 1;
}
fprintf(stderr, "Created test catalog: %s\n", settings.catalog_path);
}
SpacetimeSource spacetime = {0};
if (spacetime_create_default(&spacetime)) {
fputs("Could not create spacetime source\n", stderr);
catalog_destroy(&catalog);
return 1;
}
int result = render_frame(&settings, &catalog, &spacetime);
spacetime_destroy(&spacetime);
catalog_destroy(&catalog);
return result == 0 ? 0 : 1;
}
src/observer.c
+67
View File
@@ -0,0 +1,67 @@
#include "observer.h"
#include <math.h>
#include <stddef.h>
static const double pi = 3.14159265358979323846;
ObserverState observer_fixed_at_origin(void) {
return (ObserverState){.coordinate_time = 0.0,
.coordinate_position = {0.0, 0.0, 0.0},
.tetrad = {{1.0, 0.0, 0.0, 0.0},
{0.0, 0.0, 0.0, -1.0},
{0.0, 0.0, 1.0, 0.0},
{0.0, 1.0, 0.0, 0.0}}};
}
ObserverState observer_fixed_at_origin_look_at(double ra_deg, double dec_deg) {
const double ra = ra_deg * pi / 180.0;
const double dec = dec_deg * pi / 180.0;
const double cos_ra = cos(ra), sin_ra = sin(ra);
const double cos_dec = cos(dec), sin_dec = sin(dec);
const double forward[3] = {cos_dec * cos_ra, sin_dec, cos_dec * sin_ra};
const double up[3] = {-sin_dec * cos_ra, cos_dec, -sin_dec * sin_ra};
const double right[3] = {-sin_ra, 0.0, cos_ra};
return (ObserverState){.coordinate_time = 0.0,
.coordinate_position = {0.0, 0.0, 0.0},
.tetrad = {{1.0, 0.0, 0.0, 0.0},
{0.0, forward[0], forward[1], forward[2]},
{0.0, up[0], up[1], up[2]},
{0.0, right[0], right[1], right[2]}}};
}
int observer_static_schwarzschild_ks(double mass, double radius,
ObserverState *out) {
if (out == NULL || mass <= 0.0 || radius <= 2.0 * mass)
return -1;
const double f = 2.0 * mass / radius;
const double normalization = sqrt(1.0 - f);
*out = (ObserverState){
.coordinate_time = 0.0,
.coordinate_position = {radius, 0.0, 0.0},
/* e_(0) is the static four-velocity. e_(1) points inward; its time
* component makes the tetrad orthonormal in the KS metric. */
.tetrad = {{1.0 / normalization, 0.0, 0.0, 0.0},
{-f / normalization, -normalization, 0.0, 0.0},
{0.0, 0.0, 1.0, 0.0},
{0.0, 0.0, 0.0, 1.0}}};
return 0;
}
int observer_inward_schwarzschild_ks(double mass, double radius,
double inward_speed, ObserverState *out) {
ObserverState static_observer;
if (inward_speed < 0.0 || inward_speed >= 1.0 ||
observer_static_schwarzschild_ks(mass, radius, &static_observer))
return -1;
const double gamma = 1.0 / sqrt(1.0 - inward_speed * inward_speed);
*out = static_observer;
for (int mu = 0; mu < 4; ++mu) {
const double e0 = static_observer.tetrad[0][mu];
const double forward = static_observer.tetrad[1][mu];
out->tetrad[0][mu] = gamma * (e0 + inward_speed * forward);
out->tetrad[1][mu] = gamma * (inward_speed * e0 + forward);
}
return 0;
}
src/observer.h
+24
View File
@@ -0,0 +1,24 @@
#ifndef OBSERVER_H
#define OBSERVER_H
typedef struct {
double coordinate_time;
double coordinate_position[3];
/* e_(0), then spatial (forward, up, right), in coordinate components. */
double tetrad[4][4];
} ObserverState;
ObserverState observer_fixed_at_origin(void);
/* Point the fixed inertial observer at an ICRS-style RA/Dec direction.
* The local spatial axes remain (forward, celestial north, increasing RA). */
ObserverState observer_fixed_at_origin_look_at(double ra_deg, double dec_deg);
/* Static camera at (radius, 0, 0) in Cartesian Kerr--Schild coordinates,
* directed toward the Schwarzschild black hole at the origin. */
int observer_static_schwarzschild_ks(double mass, double radius,
ObserverState *out);
/* Camera at (radius, 0, 0), moving inward at local speed inward_speed with
* respect to the static Schwarzschild observer. */
int observer_inward_schwarzschild_ks(double mass, double radius,
double inward_speed, ObserverState *out);
#endif
src/optics.c
+183
View File
@@ -0,0 +1,183 @@
#include "optics.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef ENABLE_PNG
#include <png.h>
#endif
static double clamp(double value, double low, double high)
{
return value < low ? low : value > high ? high : value;
}
/* Wyman, Sloan, and Shirley (2013), Eq. 4: analytic fits to the CIE 1931
* 2-degree color-matching functions. Wavelength is in nanometres. */
static void cie_1931_2deg(double wavelength_nm, double xyz[3])
{
const double x1 = (wavelength_nm - 442.0) *
(wavelength_nm < 442.0 ? 0.0624 : 0.0374);
const double x2 = (wavelength_nm - 599.8) *
(wavelength_nm < 599.8 ? 0.0264 : 0.0323);
const double x3 = (wavelength_nm - 501.1) *
(wavelength_nm < 501.1 ? 0.0490 : 0.0382);
const double y1 = (wavelength_nm - 568.8) *
(wavelength_nm < 568.8 ? 0.0213 : 0.0247);
const double y2 = (wavelength_nm - 530.9) *
(wavelength_nm < 530.9 ? 0.0613 : 0.0322);
const double z1 = (wavelength_nm - 437.0) *
(wavelength_nm < 437.0 ? 0.0845 : 0.0278);
const double z2 = (wavelength_nm - 459.0) *
(wavelength_nm < 459.0 ? 0.0385 : 0.0725);
xyz[0] = 0.362 * exp(-0.5 * x1 * x1) +
1.056 * exp(-0.5 * x2 * x2) - 0.065 * exp(-0.5 * x3 * x3);
xyz[1] = 0.821 * exp(-0.5 * y1 * y1) +
0.286 * exp(-0.5 * y2 * y2);
xyz[2] = 1.217 * exp(-0.5 * z1 * z1) +
0.681 * exp(-0.5 * z2 * z2);
}
static double planck_radiance_lambda(double wavelength_m, double temperature_K)
{
const double h = 6.62607015e-34;
const double c = 299792458.0;
const double k = 1.380649e-23;
const double exponent = h * c / (wavelength_m * k * temperature_K);
return 2.0 * h * c * c /
(pow(wavelength_m, 5.0) * expm1(exponent));
}
LinearRgb blackbody_to_linear_rgb(double temperature_K)
{
/* Integrate Planck spectral radiance from 380 to 780 nm into CIE XYZ,
* then transform XYZ to linear sRGB. Results are W m^-2 sr^-1 in each
* linear-primary channel, before the catalog amplitude and exposure. */
double xyz[3] = {0.0, 0.0, 0.0};
const double wavelength_step_m = 5e-9;
if (!isfinite(temperature_K) || temperature_K <= 0.0)
return (LinearRgb){0.0, 0.0, 0.0};
for (int wavelength_nm = 380; wavelength_nm <= 780; wavelength_nm += 5) {
double matching[3];
const double radiance =
planck_radiance_lambda(wavelength_nm * 1e-9, temperature_K);
cie_1931_2deg(wavelength_nm, matching);
for (int channel = 0; channel < 3; ++channel)
xyz[channel] += radiance * matching[channel] * wavelength_step_m;
}
return (LinearRgb){
clamp(3.24096994 * xyz[0] - 1.53738318 * xyz[1] - 0.49861076 * xyz[2],
0.0, INFINITY),
clamp(-0.96924364 * xyz[0] + 1.87596750 * xyz[1] + 0.04155506 * xyz[2],
0.0, INFINITY),
clamp(0.05563008 * xyz[0] - 0.20397696 * xyz[1] + 1.05697151 * xyz[2],
0.0, INFINITY)};
}
void splat_moffat(double *hdr, int width, int height, double x, double y,
LinearRgb color, double flux,
const PointSpreadFunction *psf)
{
/* The tail omitted outside this radius contains 1e-8 of the Moffat's
* total flux. Unlike the old fixed 3-sigma box, this is both circular and
* far below the displayed HDR precision for the chosen beta. */
const double tail_fraction = 1e-8;
if (hdr == NULL || psf == NULL || flux <= 0.0 ||
psf->fwhm_pixels <= 0.0 || psf->moffat_beta <= 1.0)
return;
const double beta = psf->moffat_beta;
const double alpha = psf->fwhm_pixels /
(2.0 * sqrt(pow(2.0, 1.0 / beta) - 1.0));
const double support_radius = alpha * sqrt(
pow(tail_fraction, 1.0 / (1.0 - beta)) - 1.0);
const double support_radius_squared = support_radius * support_radius;
const int min_x = (int)floor(x - support_radius);
const int max_x = (int)ceil(x + support_radius);
const int min_y = (int)floor(y - support_radius);
const int max_y = (int)ceil(y + support_radius);
const double normalization = flux * (beta - 1.0) /
(3.14159265358979323846 * alpha * alpha);
for (int py = min_y; py <= max_y; ++py) for (int px = min_x; px <= max_x; ++px) {
if (px < 0 || px >= width || py < 0 || py >= height) continue;
double dx = (px + 0.5) - x, dy = (py + 0.5) - y;
const double radius_squared = dx * dx + dy * dy;
if (radius_squared > support_radius_squared) continue;
const double w = normalization *
pow(1.0 + radius_squared / (alpha * alpha), -beta);
double *pixel = &hdr[3 * (py * width + px)];
pixel[0] += color.r * w; pixel[1] += color.g * w; pixel[2] += color.b * w;
}
}
static unsigned char tonemap_channel(double hdr_value)
{
/* Reinhard tone mapping followed by the sRGB display transfer curve. */
const double linear = hdr_value / (1.0 + hdr_value);
const double display = linear <= 0.0031308 ? 12.92 * linear
: 1.055 * pow(linear, 1.0 / 2.4) - 0.055;
return (unsigned char)lround(255.0 * clamp(display, 0.0, 1.0));
}
static int write_tonemapped_ppm(const char *path, const double *hdr, int width,
int height)
{
FILE *file = fopen(path, "wb");
if (file == NULL) return -1;
fprintf(file, "P6\n%d %d\n255\n", width, height);
for (int i = 0; i < width * height * 3; ++i) {
const unsigned char value = tonemap_channel(hdr[i]);
if (fwrite(&value, 1, 1, file) != 1) { fclose(file); return -1; }
}
return fclose(file) == 0 ? 0 : -1;
}
#ifdef ENABLE_PNG
static int write_tonemapped_png(const char *path, const double *hdr, int width,
int height)
{
FILE *file = fopen(path, "wb");
png_structp png = NULL;
png_infop info = NULL;
unsigned char *pixels = NULL;
int result = -1;
if (file == NULL) return -1;
png = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
if (png == NULL) goto done;
info = png_create_info_struct(png);
if (info == NULL || setjmp(png_jmpbuf(png))) goto done;
pixels = malloc((size_t)width * height * 3);
if (pixels == NULL) goto done;
for (int i = 0; i < width * height * 3; ++i)
pixels[i] = tonemap_channel(hdr[i]);
png_init_io(png, file);
png_set_IHDR(png, info, (png_uint_32)width, (png_uint_32)height, 8,
PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
png_write_info(png, info);
for (int row = 0; row < height; ++row)
png_write_row(png, &pixels[(size_t)row * width * 3]);
png_write_end(png, info);
result = 0;
done:
free(pixels);
png_destroy_write_struct(&png, &info);
if (fclose(file) != 0) result = -1;
return result;
}
#endif
int write_tonemapped_image(const char *path, const double *hdr, int width, int height)
{
const size_t path_length = strlen(path);
if (path_length >= 4 && strcmp(path + path_length - 4, ".png") == 0) {
#ifdef ENABLE_PNG
return write_tonemapped_png(path, hdr, width, height);
#else
fputs("PNG output is disabled; rebuild with make ENABLE_PNG=1.\n", stderr);
return -1;
#endif
}
return write_tonemapped_ppm(path, hdr, width, height);
}
src/optics.h
+20
View File
@@ -0,0 +1,20 @@
#ifndef OPTICS_H
#define OPTICS_H
typedef struct { double r, g, b; } LinearRgb;
typedef struct {
double fwhm_pixels;
double moffat_beta;
} PointSpreadFunction;
/* Integrate a Planck spectrum into absolute linear-sRGB spectral radiance
* (W m^-2 sr^-1), before catalog amplitude and display exposure. */
LinearRgb blackbody_to_linear_rgb(double temperature_K);
void splat_moffat(double *hdr, int width, int height, double x, double y,
LinearRgb color, double flux,
const PointSpreadFunction *psf);
/* Writes PPM by default; a .png path requires a build with ENABLE_PNG=1. */
int write_tonemapped_image(const char *path, const double *hdr, int width,
int height);
#endif
src/spacetime.h
+48
View File
@@ -0,0 +1,48 @@
#ifndef SPACETIME_H
#define SPACETIME_H
typedef struct {
double alpha;
double beta[3];
double gamma[3][3];
double K[3][3];
double d_alpha[3];
double d_beta[3][3]; /* d_beta[spatial derivative][component] */
double d_gamma[3][3][3]; /* d_gamma[spatial derivative][j][k] */
} MetricData;
typedef enum {
SPACETIME_RAY_ACTIVE,
SPACETIME_RAY_ESCAPED,
SPACETIME_RAY_CAPTURED
} SpacetimeRayStatus;
typedef struct SpacetimeSource SpacetimeSource;
typedef struct {
int (*eval)(const SpacetimeSource *source, double t, const double x[3],
MetricData *metric);
SpacetimeRayStatus (*classify)(const SpacetimeSource *source, double t,
const double x[3]);
void (*destroy)(SpacetimeSource *source);
} SpacetimeOps;
struct SpacetimeSource {
const SpacetimeOps *ops;
void *context;
};
/* The selected build provides spacetime_create_default(). Named constructors
* remain available to backend-specific tests and tools. */
int spacetime_create_default(SpacetimeSource *source);
int spacetime_create_minkowski(SpacetimeSource *source, double escape_radius);
int spacetime_create_schwarzschild_ks(SpacetimeSource *source, double mass,
double escape_radius,
double capture_radius);
void spacetime_destroy(SpacetimeSource *source);
int spacetime_eval(const SpacetimeSource *source, double t, const double x[3],
MetricData *metric);
SpacetimeRayStatus spacetime_classify(const SpacetimeSource *source, double t,
const double x[3]);
#endif
src/spacetime_common.c
+22
View File
@@ -0,0 +1,22 @@
#include "spacetime.h"
#include <stddef.h>
void spacetime_destroy(SpacetimeSource *source) {
if (source != NULL && source->ops != NULL)
source->ops->destroy(source);
}
int spacetime_eval(const SpacetimeSource *source, double t, const double x[3],
MetricData *metric) {
return source == NULL || source->ops == NULL
? -1
: source->ops->eval(source, t, x, metric);
}
SpacetimeRayStatus spacetime_classify(const SpacetimeSource *source, double t,
const double x[3]) {
return source == NULL || source->ops == NULL
? SPACETIME_RAY_CAPTURED
: source->ops->classify(source, t, x);
}
src/spacetime_minkowski.c
+56
View File
@@ -0,0 +1,56 @@
#include "spacetime.h"
#include <stdlib.h>
typedef struct {
double escape_radius;
} MinkowskiContext;
static int minkowski_eval(const SpacetimeSource *source, double t,
const double x[3], MetricData *metric) {
(void)source;
(void)t;
(void)x;
*metric = (MetricData){
.alpha = 1.0,
.gamma = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}};
return 0;
}
static SpacetimeRayStatus minkowski_classify(const SpacetimeSource *source,
double t, const double x[3]) {
const MinkowskiContext *context = source->context;
const double radius_squared = x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
(void)t;
return radius_squared >= context->escape_radius * context->escape_radius
? SPACETIME_RAY_ESCAPED
: SPACETIME_RAY_ACTIVE;
}
static void minkowski_destroy(SpacetimeSource *source) {
free(source->context);
source->context = NULL;
source->ops = NULL;
}
static const SpacetimeOps minkowski_ops = {
.eval = minkowski_eval,
.classify = minkowski_classify,
.destroy = minkowski_destroy,
};
int spacetime_create_minkowski(SpacetimeSource *source, double escape_radius) {
if (source == NULL || escape_radius <= 0.0)
return -1;
MinkowskiContext *context = malloc(sizeof *context);
if (context == NULL)
return -1;
context->escape_radius = escape_radius;
source->ops = &minkowski_ops;
source->context = context;
return 0;
}
int spacetime_create_default(SpacetimeSource *source) {
return spacetime_create_minkowski(source, 1024.0);
}
src/spacetime_schwarzschild.c
+131
View File
@@ -0,0 +1,131 @@
#include "spacetime.h"
#include <math.h>
#include <stdlib.h>
typedef struct {
double mass;
double escape_radius;
double capture_radius;
} SchwarzschildKsContext;
/* Schwarzschild in ingoing Cartesian Kerr--Schild coordinates:
* g_mu_nu = eta_mu_nu + (2 M / r) l_mu l_nu, l_mu = (1, x_i / r).
* These slices are regular at r = 2 M; only the physical r = 0 singularity
* is excluded by the conservative capture cutoff. */
static int schwarzschild_ks_eval(const SpacetimeSource *source, double t,
const double x[3], MetricData *metric) {
const SchwarzschildKsContext *context = source->context;
double r2 = 0.0;
(void)t;
for (int i = 0; i < 3; ++i)
r2 += x[i] * x[i];
if (!isfinite(r2) || r2 <= 0.0)
return -1;
const double r = sqrt(r2);
const double m = context->mass;
const double f = 2.0 * m / r;
const double alpha = 1.0 / sqrt(1.0 + f);
const double beta_scale = 2.0 * m / (r * (r + 2.0 * m));
const double dbeta_scale_dr =
-4.0 * m * (r + m) / (r2 * (r + 2.0 * m) * (r + 2.0 * m));
*metric = (MetricData){.alpha = alpha};
for (int i = 0; i < 3; ++i) {
metric->beta[i] = beta_scale * x[i];
metric->d_alpha[i] = m * alpha * alpha * alpha * x[i] / (r2 * r);
for (int j = 0; j < 3; ++j) {
metric->gamma[i][j] = (i == j ? 1.0 : 0.0) +
2.0 * m * x[i] * x[j] / (r2 * r);
metric->d_beta[i][j] = beta_scale * (i == j ? 1.0 : 0.0) +
dbeta_scale_dr * x[i] * x[j] / r;
for (int k = 0; k < 3; ++k)
metric->d_gamma[i][j][k] =
2.0 * m * (((i == j ? 1.0 : 0.0) * x[k] +
(i == k ? 1.0 : 0.0) * x[j]) /
(r2 * r) -
3.0 * x[i] * x[j] * x[k] / (r2 * r2 * r));
}
}
/* K_ij = (D_i beta_j + D_j beta_i)/(2 alpha) for these stationary slices.
* Build the connection from the analytic spatial-metric derivatives above;
* this avoids finite differences in the geodesic RHS. */
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j) {
double d_beta_cov_i_j =
2.0 * m * ((i == j ? 1.0 : 0.0) / r2 -
2.0 * x[i] * x[j] / (r2 * r2));
double d_beta_cov_j_i = d_beta_cov_i_j;
double connection_term_ij = 0.0;
double connection_term_ji = 0.0;
for (int ell = 0; ell < 3; ++ell) {
double gamma_inverse_ell_k[3];
for (int k = 0; k < 3; ++k)
gamma_inverse_ell_k[k] = (ell == k ? 1.0 : 0.0) -
f / (1.0 + f) * x[ell] * x[k] / r2;
double gamma_ell_ij = 0.0, gamma_ell_ji = 0.0;
for (int k = 0; k < 3; ++k) {
gamma_ell_ij += 0.5 * gamma_inverse_ell_k[k] *
(metric->d_gamma[i][k][j] +
metric->d_gamma[j][k][i] -
metric->d_gamma[k][i][j]);
gamma_ell_ji += 0.5 * gamma_inverse_ell_k[k] *
(metric->d_gamma[j][k][i] +
metric->d_gamma[i][k][j] -
metric->d_gamma[k][j][i]);
}
const double beta_cov_ell = 2.0 * m * x[ell] / r2;
connection_term_ij += gamma_ell_ij * beta_cov_ell;
connection_term_ji += gamma_ell_ji * beta_cov_ell;
}
metric->K[i][j] = (d_beta_cov_i_j - connection_term_ij +
d_beta_cov_j_i - connection_term_ji) /
(2.0 * alpha);
}
return 0;
}
static SpacetimeRayStatus schwarzschild_ks_classify(
const SpacetimeSource *source, double t, const double x[3]) {
const SchwarzschildKsContext *context = source->context;
double r2 = 0.0;
(void)t;
for (int i = 0; i < 3; ++i)
r2 += x[i] * x[i];
if (!isfinite(r2) || r2 <= context->capture_radius * context->capture_radius)
return SPACETIME_RAY_CAPTURED;
return r2 >= context->escape_radius * context->escape_radius
? SPACETIME_RAY_ESCAPED
: SPACETIME_RAY_ACTIVE;
}
static void schwarzschild_ks_destroy(SpacetimeSource *source) {
free(source->context);
source->context = NULL;
source->ops = NULL;
}
static const SpacetimeOps schwarzschild_ks_ops = {
.eval = schwarzschild_ks_eval,
.classify = schwarzschild_ks_classify,
.destroy = schwarzschild_ks_destroy,
};
int spacetime_create_schwarzschild_ks(SpacetimeSource *source, double mass,
double escape_radius,
double capture_radius) {
if (source == NULL || mass <= 0.0 || escape_radius <= 2.0 * mass ||
capture_radius <= 0.0 || capture_radius >= 2.0 * mass ||
capture_radius >= escape_radius)
return -1;
SchwarzschildKsContext *context = malloc(sizeof *context);
if (context == NULL)
return -1;
*context = (SchwarzschildKsContext){mass, escape_radius, capture_radius};
source->ops = &schwarzschild_ks_ops;
source->context = context;
return 0;
}
int spacetime_create_default(SpacetimeSource *source) {
return spacetime_create_schwarzschild_ks(source, 1.0, 256.0, 1.5);
}