wyj/GR-raytracing
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
ab7d879f6cdb8a0f5d2d8ec4a4c2b40c0e45e070
GR-raytracing/src/frame.c
T
wyj e1ec480669 HIP: accelerate PSF accumulation and restore parallel producers 
Cooperate across 32 lanes per PSF and use two completion-protected staging slots with complete batch timing. Restore coarse OpenMP event production while serializing shared GPU submissions and direct fallback boundaries.

Add bounded benchmarks, streaming and renderer regressions, and preserve validation evidence and ownership documentation.
2026-09-06 21:38:53 -04:00

1653 lines
68 KiB
C
Raw Blame History

#include "frame.h"
#ifdef PSF_BACKEND_HIP
#include "hip_psf.h"
#endif
#include "optics.h"
#include <limits.h>
#include <math.h>
#include <omp.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifndef FRAME_PSF_EVENT_SINK
#define FRAME_PSF_EVENT_SINK 1
#endif
/* Numerical metric backends may reserve substantial memory for slabs and
* thread-local evaluators, so they retain this private-HDR allocation budget.
* Analytic backends deliberately use all OpenMP render workers instead. */
#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}, 0, 0};
triangles[next_triangle++] =
(LensTriangle){{top_left, bottom_right, top_right}, 0, 0};
}
*mesh = (FrameLensMesh){.vertices = vertices,
.triangles = triangles,
.vertex_count = vertex_count,
.vertex_capacity = vertex_count,
.triangle_count = triangle_count,
.triangle_capacity = 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;
vertex->traced = 1;
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;
}
enum { PSF_EVENT_SINK_CAPACITY = 16384 };
typedef struct {
PsfCachedEvent *events;
size_t count;
double *hdr;
int width, height;
const PsfKernelCache *cache;
#ifdef PSF_BACKEND_HIP
HipPsfSink *hip;
omp_lock_t *hip_lock; /* Shared sink lock; acquired only per chunk/fallback. */
int borrowed_hip;
double submit_seconds;
char hip_message[256];
HipPsfTiming hip_timing;
#endif
int failed;
} PsfEventSink;
#ifdef PSF_BACKEND_HIP
static void psf_event_sink_mark_failed(PsfEventSink *sink, const char *stage) {
sink->failed = 1;
fprintf(stderr, "HIP PSF backend failed during %s: %s\n", stage,
sink->hip_message[0] == '\0' ? "unknown HIP error" : sink->hip_message);
}
#endif
static int psf_event_sink_init(PsfEventSink *sink, double *hdr, int width,
int height, const PsfKernelCache *cache) {
*sink = (PsfEventSink){.hdr = hdr, .width = width, .height = height,
.cache = cache};
if (cache != NULL) {
sink->events = malloc(PSF_EVENT_SINK_CAPACITY * sizeof *sink->events);
if (sink->events == NULL) {
#ifdef PSF_BACKEND_HIP
snprintf(sink->hip_message, sizeof sink->hip_message, "host event allocation failed");
psf_event_sink_mark_failed(sink, "initialization");
return -1;
#else
return 0;
#endif
}
#ifdef PSF_BACKEND_HIP
if (hip_psf_sink_create(&sink->hip, width, height, cache,
PSF_EVENT_SINK_CAPACITY, sink->hip_message,
sizeof sink->hip_message)) {
psf_event_sink_mark_failed(sink, "initialization");
free(sink->events);
sink->events = NULL;
return -1;
}
#endif
}
return 0;
}
static int psf_event_sink_flush_unlocked(PsfEventSink *sink) {
if (sink->failed)
return -1;
#ifdef PSF_BACKEND_HIP
if (sink->hip != NULL) {
if (hip_psf_sink_submit(sink->hip, sink->events, sink->count,
sink->hip_message, sizeof sink->hip_message)) {
psf_event_sink_mark_failed(sink, "event submission");
return -1;
}
sink->count = 0;
return 0;
}
#endif
for (size_t i = 0; i < sink->count; ++i)
splat_prepared_cached_event(sink->hdr, sink->width, sink->height,
&sink->events[i], sink->cache);
sink->count = 0;
return 0;
}
static int psf_event_sink_flush(PsfEventSink *sink) {
#ifdef PSF_BACKEND_HIP
const double start = omp_get_wtime();
if (sink->hip_lock) omp_set_lock(sink->hip_lock);
#endif
const int result = psf_event_sink_flush_unlocked(sink);
#ifdef PSF_BACKEND_HIP
if (sink->hip_lock) omp_unset_lock(sink->hip_lock);
sink->submit_seconds += omp_get_wtime() - start;
#endif
return result;
}
/* A borrowed sink calls this only while holding hip_lock across the entire
* download -> CPU fallback -> upload transaction. */
static int psf_event_sink_finish_for_cpu(PsfEventSink *sink) {
#ifdef PSF_BACKEND_HIP
if (psf_event_sink_flush_unlocked(sink))
#else
if (psf_event_sink_flush(sink))
#endif
return -1;
#ifdef PSF_BACKEND_HIP
if (sink->hip != NULL && hip_psf_sink_finish(sink->hip, sink->hdr,
sink->hip_message,
sizeof sink->hip_message)) {
psf_event_sink_mark_failed(sink, "HDR download");
return -1;
}
#endif
return 0;
}
static int psf_event_sink_resume_gpu(PsfEventSink *sink) {
#ifdef PSF_BACKEND_HIP
if (sink->hip != NULL && hip_psf_sink_load_hdr(sink->hip, sink->hdr,
sink->hip_message,
sizeof sink->hip_message)) {
psf_event_sink_mark_failed(sink, "HDR upload");
return -1;
}
#else
(void)sink;
#endif
return 0;
}
static int psf_event_sink_destroy(PsfEventSink *sink) {
#ifdef PSF_BACKEND_HIP
if (sink->borrowed_hip) {
const int result = psf_event_sink_flush(sink);
free(sink->events);
return result;
}
#endif
int result = psf_event_sink_finish_for_cpu(sink);
#ifdef PSF_BACKEND_HIP
if (sink->hip != NULL && hip_psf_sink_get_timing(sink->hip, &sink->hip_timing)) {
psf_event_sink_mark_failed(sink, "timing collection");
result = -1;
}
hip_psf_sink_destroy(sink->hip);
#endif
free(sink->events);
return result;
}
#if FRAME_PSF_EVENT_SINK || defined(PSF_BACKEND_HIP)
static void psf_event_sink_emit(PsfEventSink *sink, const PsfCachedEvent *event) {
if (sink->failed)
return;
if (sink->events == NULL) {
splat_prepared_cached_event(sink->hdr, sink->width, sink->height, event, sink->cache);
return;
}
if (sink->count == PSF_EVENT_SINK_CAPACITY)
(void)psf_event_sink_flush(sink);
if (sink->failed)
return;
sink->events[sink->count++] = *event;
}
#endif
typedef struct {
size_t a, b, triangle;
unsigned int side;
} MeshEdge;
static int compare_mesh_edge(const void *left, const void *right) {
const MeshEdge *a = left, *b = right;
if (a->a != b->a)
return a->a < b->a ? -1 : 1;
if (a->b != b->b)
return a->b < b->b ? -1 : 1;
return 0;
}
static double image_edge_length(const LensVertex *a, const LensVertex *b) {
return hypot(a->image_x - b->image_x, a->image_y - b->image_y);
}
static double image_triangle_area(const LensVertex *a, const LensVertex *b,
const LensVertex *c) {
return 0.5 * fabs((b->image_x - a->image_x) * (c->image_y - a->image_y) -
(b->image_y - a->image_y) * (c->image_x - a->image_x));
}
static unsigned int longest_side(const FrameLensMesh *mesh,
const LensTriangle *triangle) {
unsigned int best = 0;
double best_length = -1.0;
for (unsigned int side = 0; side < 3; ++side) {
const double length = image_edge_length(
&mesh->vertices[triangle->vertex[side]],
&mesh->vertices[triangle->vertex[(side + 1) % 3]]);
if (length > best_length) {
best_length = length;
best = side;
}
}
return best;
}
static int ensure_samples(FrameLensMesh *mesh, size_t count) {
if (count <= mesh->sample_capacity)
return 0;
size_t capacity = mesh->sample_capacity ? mesh->sample_capacity : 16;
while (capacity < count)
capacity *= 2;
FrameSample *samples = realloc(mesh->samples, capacity * sizeof *samples);
if (samples == NULL)
return -1;
mesh->samples = samples;
mesh->sample_capacity = capacity;
return 0;
}
static int ensure_vertices(FrameLensMesh *mesh, size_t count) {
if (count <= mesh->vertex_capacity)
return 0;
size_t capacity = mesh->vertex_capacity ? mesh->vertex_capacity : 16;
while (capacity < count)
capacity *= 2;
LensVertex *vertices = realloc(mesh->vertices, capacity * sizeof *vertices);
if (vertices == NULL)
return -1;
mesh->vertices = vertices;
mesh->vertex_capacity = capacity;
return 0;
}
static int all_vertices_traced(const FrameLensMesh *mesh) {
for (size_t i = 0; i < mesh->vertex_count; ++i)
if (!mesh->vertices[i].traced)
return 0;
return 1;
}
static int terminal_mismatch(const LensVertex *a, const LensVertex *b,
const LensVertex *c) {
int escaped = 0, captured = 0;
const LensVertex *vertices[] = {a, b, c};
for (size_t i = 0; i < 3; ++i) {
escaped |= vertices[i]->status == RAY_ENDPOINT_ESCAPED;
captured |= vertices[i]->status == RAY_ENDPOINT_CAPTURED;
}
return escaped && captured;
}
static int add_sample(FrameLensMesh *mesh, const FrameSample *sample) {
if (ensure_samples(mesh, mesh->sample_count + 1))
return -1;
mesh->samples[mesh->sample_count++] = *sample;
return 0;
}
static size_t probe_hash(size_t a, size_t b) {
uint64_t value = (uint64_t)a * UINT64_C(0x9e3779b185ebca87) ^
(uint64_t)b * UINT64_C(0xc2b2ae3d27d4eb4f);
value ^= value >> 33;
return (size_t)value;
}
static int prepare_probe_index(FrameLensMesh *mesh) {
size_t capacity = 16;
while (capacity < mesh->triangle_count * 2)
capacity *= 2;
if (capacity > mesh->probe_slot_capacity) {
size_t *slots = realloc(mesh->probe_slots, capacity * sizeof *slots);
if (slots == NULL)
return -1;
mesh->probe_slots = slots;
mesh->probe_slot_capacity = capacity;
}
memset(mesh->probe_slots, 0, mesh->probe_slot_capacity * sizeof *mesh->probe_slots);
return 0;
}
static size_t find_probe(const FrameLensMesh *mesh, size_t a, size_t b) {
if (a > b) { size_t swap = a; a = b; b = swap; }
if (mesh->probe_slot_capacity == 0)
return SIZE_MAX;
size_t slot = probe_hash(a, b) & (mesh->probe_slot_capacity - 1);
while (mesh->probe_slots[slot] != 0) {
const size_t sample_id = mesh->probe_slots[slot] - 1;
const FrameSample *sample = &mesh->samples[sample_id];
size_t x = sample->edge_vertex[0], y = sample->edge_vertex[1];
if (x > y) { size_t swap = x; x = y; y = swap; }
if (x == a && y == b)
return sample_id;
slot = (slot + 1) & (mesh->probe_slot_capacity - 1);
}
return SIZE_MAX;
}
static void index_probe(FrameLensMesh *mesh, size_t sample_id) {
FrameSample *sample = &mesh->samples[sample_id];
size_t a = sample->edge_vertex[0], b = sample->edge_vertex[1];
if (a > b) { size_t swap = a; a = b; b = swap; }
size_t slot = probe_hash(a, b) & (mesh->probe_slot_capacity - 1);
while (mesh->probe_slots[slot] != 0)
slot = (slot + 1) & (mesh->probe_slot_capacity - 1);
mesh->probe_slots[slot] = sample_id + 1;
}
int frame_lens_mesh_prepare_generation(FrameLensMesh *mesh,
const RefinementConfig *config) {
if (mesh == NULL || config == NULL || mesh->sample_count != 0)
return -1;
mesh->samples_include_probes = 0;
if (config->max_level > 0 && prepare_probe_index(mesh))
return -1;
const int has_untraced_vertices = !all_vertices_traced(mesh);
if (has_untraced_vertices) {
for (size_t i = 0; i < mesh->vertex_count; ++i)
if (!mesh->vertices[i].traced &&
add_sample(mesh, &(FrameSample){.kind = FRAME_SAMPLE_VERTEX,
.vertex_id = i,
.vertex = mesh->vertices[i]}))
return -1;
}
if (config->max_level == 0)
/* Coarse vertices still need tracing when refinement is disabled. */
return (int)mesh->sample_count;
/* Every generation may batch newly inserted vertices with probes for its
* new leaves: probe positions depend only on image-plane geometry. Their
* endpoints are considered only after this complete generation finishes. */
for (size_t i = 0; i < mesh->triangle_count; ++i) {
const LensTriangle *triangle = &mesh->triangles[i];
if (triangle->level >= config->max_level || triangle->evaluated)
continue;
const unsigned int first_side = longest_side(mesh, triangle);
const unsigned int side_count = terminal_mismatch(
&mesh->vertices[triangle->vertex[0]],
&mesh->vertices[triangle->vertex[1]],
&mesh->vertices[triangle->vertex[2]])
? 3
: 1;
for (unsigned int offset = 0; offset < side_count; ++offset) {
const unsigned int side = (first_side + offset) % 3;
const size_t a = triangle->vertex[side];
const size_t b = triangle->vertex[(side + 1) % 3];
if (find_probe(mesh, a, b) != SIZE_MAX)
continue;
FrameSample probe = {.kind = FRAME_SAMPLE_PROBE, .edge_vertex = {a, b}};
const LensVertex *left = &mesh->vertices[a];
const LensVertex *right = &mesh->vertices[b];
probe.vertex.image_x = 0.5 * (left->image_x + right->image_x);
probe.vertex.image_y = 0.5 * (left->image_y + right->image_y);
for (int axis = 0; axis < 3; ++axis)
probe.vertex.camera_direction[axis] =
left->camera_direction[axis] + right->camera_direction[axis];
if (normalize(probe.vertex.camera_direction) == 0.0)
return -1;
if (add_sample(mesh, &probe))
return -1;
index_probe(mesh, mesh->sample_count - 1);
}
}
mesh->samples_include_probes = mesh->sample_count != 0;
return (int)mesh->sample_count;
}
const FrameSample *frame_lens_mesh_samples(const FrameLensMesh *mesh,
size_t *count) {
if (count != NULL)
*count = mesh == NULL ? 0 : mesh->sample_count;
return mesh == NULL ? NULL : mesh->samples;
}
int frame_lens_mesh_install_sample(FrameLensMesh *mesh, size_t sample_id,
const RayEndpoint *endpoint) {
if (mesh == NULL || endpoint == NULL || sample_id >= mesh->sample_count)
return -1;
FrameSample *sample = &mesh->samples[sample_id];
LensVertex *vertex = sample->kind == FRAME_SAMPLE_VERTEX
? &mesh->vertices[sample->vertex_id]
: &sample->vertex;
vertex->status = endpoint->status;
vertex->traced = 1;
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 direction_angle(const double a[3], const double b[3]) {
const double product = fmax(-1.0, fmin(1.0, dot(a, b)));
return acos(product);
}
static double spherical_signed_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(dot(a, b_cross_c),
1.0 + dot(a, b) + dot(b, c) + dot(c, a));
}
/* Returns whether the discrete source/image solid-angle ratio is available.
* A zero output parity is a valid, critical (zero-Jacobian) result; captured
* and degenerate image triangles have no reliable parity. */
static int discrete_jacobian(const FrameLensMesh *mesh,
const LensTriangle *triangle, double *value,
signed char *parity) {
const LensVertex *a = &mesh->vertices[triangle->vertex[0]];
const LensVertex *b = &mesh->vertices[triangle->vertex[1]];
const LensVertex *c = &mesh->vertices[triangle->vertex[2]];
if (a->status != RAY_ENDPOINT_ESCAPED || b->status != RAY_ENDPOINT_ESCAPED ||
c->status != RAY_ENDPOINT_ESCAPED)
return 0;
const double image_area = spherical_signed_area(
a->camera_direction, b->camera_direction, c->camera_direction);
if (!isfinite(image_area) || fabs(image_area) <= 1e-15)
return 0;
const double source_area = spherical_signed_area(
a->n_infinity, b->n_infinity, c->n_infinity);
const double jacobian = source_area / image_area;
if (!isfinite(jacobian))
return 0;
*value = jacobian;
*parity = jacobian > 0.0 ? 1 : jacobian < 0.0 ? -1 : 0;
return 1;
}
static int probe_requires_split(const FrameLensMesh *mesh,
const LensTriangle *triangle,
unsigned int side,
const RefinementConfig *config) {
const LensVertex *a = &mesh->vertices[triangle->vertex[side]];
const LensVertex *b = &mesh->vertices[triangle->vertex[(side + 1) % 3]];
const LensVertex *c = &mesh->vertices[triangle->vertex[(side + 2) % 3]];
if (terminal_mismatch(a, b, c))
return 1;
const size_t probe_id = find_probe(mesh, triangle->vertex[side],
triangle->vertex[(side + 1) % 3]);
if (probe_id == SIZE_MAX)
return 0;
const LensVertex *probe = &mesh->samples[probe_id].vertex;
if ((probe->status == RAY_ENDPOINT_ESCAPED) !=
(a->status == RAY_ENDPOINT_ESCAPED) ||
(probe->status == RAY_ENDPOINT_ESCAPED) !=
(b->status == RAY_ENDPOINT_ESCAPED))
return (probe->status == RAY_ENDPOINT_ESCAPED ||
a->status == RAY_ENDPOINT_ESCAPED || b->status == RAY_ENDPOINT_ESCAPED) &&
(probe->status == RAY_ENDPOINT_CAPTURED ||
a->status == RAY_ENDPOINT_CAPTURED || b->status == RAY_ENDPOINT_CAPTURED);
if (a->status != RAY_ENDPOINT_ESCAPED || b->status != RAY_ENDPOINT_ESCAPED ||
probe->status != RAY_ENDPOINT_ESCAPED)
return 0;
double predicted[3] = {a->n_infinity[0] + b->n_infinity[0],
a->n_infinity[1] + b->n_infinity[1],
a->n_infinity[2] + b->n_infinity[2]};
if (normalize(predicted) == 0.0)
return 1;
const double error = direction_angle(predicted, probe->n_infinity);
/* The relative error is normalized by the image triangle's own angular
* scale (the camera-direction span of its longest edge), not by the
* source-side lens mapping. */
const double scale = fmax(direction_angle(a->camera_direction,
b->camera_direction), 1e-15);
return error > config->angle_absolute_rad &&
error / scale > config->angle_relative;
}
static int triangle_allows_children(const FrameLensMesh *mesh,
const LensTriangle *triangle,
const RefinementConfig *config) {
const LensVertex *a = &mesh->vertices[triangle->vertex[0]];
const LensVertex *b = &mesh->vertices[triangle->vertex[1]];
const LensVertex *c = &mesh->vertices[triangle->vertex[2]];
const double edge = fmax(image_edge_length(a, b),
fmax(image_edge_length(b, c), image_edge_length(c, a)));
const double area = image_triangle_area(a, b, c);
return edge > config->min_edge_pixels && area > config->min_area_pixels2;
}
static LensVertex midpoint_vertex(const LensVertex *a, const LensVertex *b) {
LensVertex result = {.image_x = 0.5 * (a->image_x + b->image_x),
.image_y = 0.5 * (a->image_y + b->image_y)};
for (int axis = 0; axis < 3; ++axis)
result.camera_direction[axis] = a->camera_direction[axis] + b->camera_direction[axis];
normalize(result.camera_direction);
return result;
}
static double image_triangle_quality(const FrameLensMesh *mesh, size_t a,
size_t b, size_t c) {
const LensVertex *va = &mesh->vertices[a];
const LensVertex *vb = &mesh->vertices[b];
const LensVertex *vc = &mesh->vertices[c];
const double ab = image_edge_length(va, vb);
const double bc = image_edge_length(vb, vc);
const double ca = image_edge_length(vc, va);
const double denominator = ab * ab + bc * bc + ca * ca;
return denominator > 0.0
? 4.0 * sqrt(3.0) * image_triangle_area(va, vb, vc) / denominator
: 0.0;
}
static double minimum_child_quality(const FrameLensMesh *mesh,
size_t children[3][3]) {
double quality = INFINITY;
for (size_t i = 0; i < 3; ++i)
quality = fmin(quality, image_triangle_quality(mesh, children[i][0],
children[i][1],
children[i][2]));
return quality;
}
static int append_triangle(LensTriangle *triangles, size_t *count,
size_t capacity, size_t a, size_t b, size_t c,
unsigned int level, int evaluated) {
if (*count >= capacity)
return -1;
triangles[(*count)++] = (LensTriangle){{a, b, c}, level, evaluated};
return 0;
}
static int append_triangle_with_parent_winding(
LensTriangle *triangles, size_t *count, size_t capacity,
const FrameLensMesh *mesh, const LensTriangle *parent, size_t a, size_t b,
size_t c, unsigned int level, int evaluated) {
const LensVertex *p0 = &mesh->vertices[parent->vertex[0]];
const LensVertex *p1 = &mesh->vertices[parent->vertex[1]];
const LensVertex *p2 = &mesh->vertices[parent->vertex[2]];
const LensVertex *v0 = &mesh->vertices[a];
const LensVertex *v1 = &mesh->vertices[b];
const LensVertex *v2 = &mesh->vertices[c];
const double parent_winding =
(p1->image_x - p0->image_x) * (p2->image_y - p0->image_y) -
(p1->image_y - p0->image_y) * (p2->image_x - p0->image_x);
const double child_winding =
(v1->image_x - v0->image_x) * (v2->image_y - v0->image_y) -
(v1->image_y - v0->image_y) * (v2->image_x - v0->image_x);
if (parent_winding * child_winding < 0.0) {
const size_t swap = b;
b = c;
c = swap;
}
return append_triangle(triangles, count, capacity, a, b, c, level,
evaluated);
}
int frame_lens_mesh_finish_generation(FrameLensMesh *mesh,
const RefinementConfig *config) {
if (mesh == NULL || config == NULL || mesh->sample_count == 0)
return -1;
for (size_t i = 0; i < mesh->sample_count; ++i)
if (!mesh->samples[i].vertex.traced &&
mesh->samples[i].kind == FRAME_SAMPLE_PROBE)
return -1;
if (!mesh->samples_include_probes) {
mesh->sample_count = 0;
mesh->samples_include_probes = 0;
return 0;
}
for (size_t i = 0; i < mesh->triangle_count; ++i)
if (mesh->triangles[i].level < config->max_level)
mesh->triangles[i].evaluated = 1;
const size_t edge_count = mesh->triangle_count * 3;
MeshEdge *edges = calloc(edge_count, sizeof *edges);
unsigned char *requested = calloc(edge_count, sizeof *requested);
unsigned char *allowed = calloc(mesh->triangle_count, sizeof *allowed);
signed char *parity = calloc(mesh->triangle_count, sizeof *parity);
double *jacobians = calloc(mesh->triangle_count, sizeof *jacobians);
/* Sorted edge groups assign one midpoint to every incident triangle side.
* Keep that relation in the original triangle-side order so child emission
* stays O(T), rather than scanning every sorted edge for every child side. */
size_t *side_midpoints = malloc(edge_count * sizeof *side_midpoints);
if (edges == NULL || requested == NULL || allowed == NULL || parity == NULL ||
jacobians == NULL || side_midpoints == NULL) {
free(edges); free(requested); free(allowed); free(parity); free(jacobians);
free(side_midpoints);
return -1;
}
for (size_t i = 0; i < edge_count; ++i)
side_midpoints[i] = SIZE_MAX;
for (size_t i = 0; i < mesh->triangle_count; ++i) {
const LensTriangle *triangle = &mesh->triangles[i];
allowed[i] = triangle->level < config->max_level &&
triangle_allows_children(mesh, triangle, config);
for (unsigned int side = 0; side < 3; ++side) {
size_t a = triangle->vertex[side], b = triangle->vertex[(side + 1) % 3];
if (a > b) { size_t swap = a; a = b; b = swap; }
edges[3 * i + side] = (MeshEdge){a, b, i, side};
}
if (allowed[i]) {
if (terminal_mismatch(&mesh->vertices[triangle->vertex[0]],
&mesh->vertices[triangle->vertex[1]],
&mesh->vertices[triangle->vertex[2]])) {
/* Capture is discontinuous across the shadow boundary. Red-refine
* directly so its image-plane scale halves every generation while
* preserving the parent triangle's shape. */
for (unsigned int side = 0; side < 3; ++side)
requested[3 * i + side] = 1;
} else {
const unsigned int side = longest_side(mesh, triangle);
requested[3 * i + side] =
probe_requires_split(mesh, triangle, side, config);
}
(void)discrete_jacobian(mesh, triangle, &jacobians[i], &parity[i]);
}
}
qsort(edges, edge_count, sizeof *edges, compare_mesh_edge);
/* A fold is selected only when its adjacent discrete parities disagree and
* at least one of those leaves is close enough to the critical curve. */
for (size_t first = 0; first < edge_count;) {
size_t last = first + 1;
while (last < edge_count && edges[last].a == edges[first].a &&
edges[last].b == edges[first].b)
++last;
if (last - first == 2) {
const size_t left = edges[first].triangle;
const size_t right = edges[first + 1].triangle;
if (parity[left] != 0 && parity[right] != 0 &&
parity[left] != parity[right] &&
fmin(fabs(jacobians[left]), fabs(jacobians[right])) <
config->jacobian_minimum) {
if (allowed[left] && allowed[right]) {
requested[3 * left + edges[first].side] = 1;
requested[3 * right + edges[first + 1].side] = 1;
}
}
}
first = last;
}
/* A requested interior edge is split by both incident leaves, preserving a
* conforming mesh. If either side has reached its geometric limit, reject
* the whole edge instead of introducing a T-junction. */
for (size_t first = 0; first < edge_count;) {
size_t last = first + 1;
while (last < edge_count && edges[last].a == edges[first].a &&
edges[last].b == edges[first].b)
++last;
int any = 0, possible = 1;
for (size_t i = first; i < last; ++i) {
const MeshEdge *edge = &edges[i];
any |= requested[3 * edge->triangle + edge->side] != 0;
possible &= allowed[edge->triangle] != 0;
}
if (any && possible)
for (size_t i = first; i < last; ++i)
requested[3 * edges[i].triangle + edges[i].side] = 1;
else if (any)
for (size_t i = first; i < last; ++i)
requested[3 * edges[i].triangle + edges[i].side] = 0;
first = last;
}
size_t split_edges = 0;
for (size_t i = 0; i < edge_count; ++i)
split_edges += requested[3 * edges[i].triangle + edges[i].side] != 0;
if (split_edges == 0) {
mesh->sample_count = 0;
mesh->samples_include_probes = 0;
free(edges); free(requested); free(allowed); free(parity); free(jacobians);
free(side_midpoints);
return 0;
}
/* Allocate a single stable midpoint vertex for each requested edge group. */
size_t midpoint_count = 0;
for (size_t first = 0; first < edge_count;) {
size_t last = first + 1;
while (last < edge_count && edges[last].a == edges[first].a &&
edges[last].b == edges[first].b)
++last;
int any = 0;
for (size_t i = first; i < last; ++i)
any |= requested[3 * edges[i].triangle + edges[i].side] != 0;
if (any) ++midpoint_count;
first = last;
}
if (ensure_vertices(mesh, mesh->vertex_count + midpoint_count)) {
free(edges); free(requested); free(allowed); free(parity); free(jacobians);
free(side_midpoints);
return -1;
}
size_t next_vertex = mesh->vertex_count;
for (size_t first = 0; first < edge_count;) {
size_t last = first + 1;
while (last < edge_count && edges[last].a == edges[first].a &&
edges[last].b == edges[first].b)
++last;
int any = 0;
for (size_t i = first; i < last; ++i)
any |= requested[3 * edges[i].triangle + edges[i].side] != 0;
if (any) {
const size_t probe_id = find_probe(mesh, edges[first].a, edges[first].b);
LensVertex midpoint = midpoint_vertex(&mesh->vertices[edges[first].a],
&mesh->vertices[edges[first].b]);
if (probe_id != SIZE_MAX)
midpoint = mesh->samples[probe_id].vertex;
mesh->vertices[next_vertex++] = midpoint;
for (size_t i = first; i < last; ++i)
side_midpoints[3 * edges[i].triangle + edges[i].side] = next_vertex - 1;
}
first = last;
}
const size_t old_count = mesh->triangle_count;
LensTriangle *children = calloc(old_count * 4, sizeof *children);
if (children == NULL) {
free(edges); free(requested); free(allowed); free(parity); free(jacobians);
free(side_midpoints);
return -1;
}
size_t child_count = 0;
for (size_t t = 0; t < old_count; ++t) {
const LensTriangle *triangle = &mesh->triangles[t];
size_t middle[3] = {SIZE_MAX, SIZE_MAX, SIZE_MAX};
unsigned int count = 0;
for (unsigned int side = 0; side < 3; ++side) {
if (!requested[3 * t + side]) continue;
++count;
middle[side] = side_midpoints[3 * t + side];
}
const size_t a = triangle->vertex[0], b = triangle->vertex[1], c = triangle->vertex[2];
const unsigned int level = triangle->level + 1;
if (count == 0)
append_triangle(children, &child_count, old_count * 4, a, b, c, triangle->level,
triangle->evaluated);
else if (count == 1) {
unsigned int side = middle[0] != SIZE_MAX ? 0 : middle[1] != SIZE_MAX ? 1 : 2;
const size_t v0 = triangle->vertex[side];
const size_t v1 = triangle->vertex[(side + 1) % 3];
const size_t other = triangle->vertex[(side + 2) % 3];
append_triangle(children, &child_count, old_count * 4, v0, middle[side], other, level, 0);
append_triangle(children, &child_count, old_count * 4, middle[side], v1, other, level, 0);
} else if (count == 2) {
size_t ab = middle[0], bc = middle[1], ca = middle[2];
const unsigned int missing = ab == SIZE_MAX ? 0 : bc == SIZE_MAX ? 1 : 2;
size_t first[3][3], second[3][3];
if (missing == 0) {
memcpy(first, (size_t[3][3]){{c, bc, ca}, {a, b, ca}, {b, bc, ca}},
sizeof first);
memcpy(second, (size_t[3][3]){{c, bc, ca}, {a, b, bc}, {a, bc, ca}},
sizeof second);
} else if (missing == 1) {
memcpy(first, (size_t[3][3]){{a, ab, ca}, {b, c, ab}, {c, ca, ab}},
sizeof first);
memcpy(second, (size_t[3][3]){{a, ab, ca}, {b, c, ca}, {b, ca, ab}},
sizeof second);
} else {
memcpy(first, (size_t[3][3]){{b, ab, bc}, {a, ab, c}, {ab, bc, c}},
sizeof first);
memcpy(second, (size_t[3][3]){{b, ab, bc}, {a, ab, bc}, {a, bc, c}},
sizeof second);
}
size_t (*chosen)[3] =
minimum_child_quality(mesh, first) >= minimum_child_quality(mesh, second)
? first
: second;
for (size_t child = 0; child < 3; ++child)
append_triangle_with_parent_winding(
children, &child_count, old_count * 4, mesh, triangle,
chosen[child][0], chosen[child][1], chosen[child][2], level, 0);
} else { /* Three requested edges: red refinement. */
const size_t ab = middle[0], bc = middle[1], ca = middle[2];
if (ab == SIZE_MAX || bc == SIZE_MAX || ca == SIZE_MAX) {
append_triangle(children, &child_count, old_count * 4, a, b, c,
triangle->level, 1);
} else {
append_triangle(children, &child_count, old_count * 4, a, ab, ca, level, 0);
append_triangle(children, &child_count, old_count * 4, ab, b, bc, level, 0);
append_triangle(children, &child_count, old_count * 4, ca, bc, c, level, 0);
append_triangle(children, &child_count, old_count * 4, ab, bc, ca, level, 0);
}
}
}
free(mesh->triangles);
mesh->triangles = children;
mesh->triangle_count = child_count;
mesh->triangle_capacity = old_count * 4;
mesh->vertex_count = next_vertex;
mesh->sample_count = 0;
mesh->samples_include_probes = 0;
free(edges); free(requested); free(allowed); free(parity); free(jacobians);
free(side_midpoints);
return (int)midpoint_count;
}
int frame_lens_mesh_refine(FrameLensMesh *mesh,
const SpacetimeSource *spacetime,
const ObserverState *observer,
const GeodesicTraceConfig *trace,
const RefinementConfig *config) {
return frame_lens_mesh_refine_with_progress(mesh, spacetime, observer, trace,
config, NULL, NULL);
}
int frame_lens_mesh_refine_with_progress(
FrameLensMesh *mesh, const SpacetimeSource *spacetime,
const ObserverState *observer, const GeodesicTraceConfig *trace,
const RefinementConfig *config, FrameRefinementProgressCallback callback,
void *context) {
if (mesh == NULL || spacetime == NULL || observer == NULL || trace == NULL ||
config == NULL)
return -1;
for (size_t generation = 0;; ++generation) {
const int requested = frame_lens_mesh_prepare_generation(mesh, config);
if (requested < 0) return -1;
if (requested == 0) return 0;
if (callback != NULL)
callback(context, generation, mesh->sample_count, mesh->vertex_count,
mesh->triangle_count, 0, 0);
#pragma omp parallel for schedule(static)
for (size_t i = 0; i < mesh->sample_count; ++i) {
const FrameSample *sample = &mesh->samples[i];
const RayEndpoint endpoint = geodesic_trace_past(
spacetime, observer, sample->vertex.camera_direction, trace);
/* Each request has a distinct destination vertex or probe slot. */
(void)frame_lens_mesh_install_sample(mesh, i, &endpoint);
}
const int added = frame_lens_mesh_finish_generation(mesh, config);
if (added < 0) return -1;
if (callback != NULL)
callback(context, generation, 0, mesh->vertex_count, mesh->triangle_count,
added, 1);
}
}
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_weights(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);
if (area < 1e-14)
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;
/* A thin source triangle can admit an exterior point through the tolerant
* edge test. Unsigned subareas then do not partition the total area, and
* interpolating with their ratios can move an image outside its triangle.
* Reject that case before normalizing roundoff in valid convex weights.
* The 1e-8 bound is dimensionless; see the inverse-map design notes. */
const double weight_sum = weights[0] + weights[1] + weights[2];
if (!isfinite(weight_sum) || weight_sum <= 0.0 ||
fabs(weight_sum - 1.0) > 1e-8)
return -1;
for (int i = 0; i < 3; ++i)
weights[i] /= weight_sum;
return 0;
}
static int spherical_barycentric(const double point[3], const double a[3],
const double b[3], const double c[3],
double weights[3]) {
double edge_cross[3];
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;
}
return spherical_barycentric_weights(point, a, b, c, weights);
}
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;
}
typedef struct {
const LensVertex *vertex[3];
const LensTriangle *triangle;
size_t triangle_index;
double *hdr;
int width, height;
double exposure, magnification;
const PointSpreadFunction *psf;
const PsfKernelCache *psf_cache;
double max_cache_psf_flux;
double psf_relative_tail;
double psf_min_y;
PsfEventSink *event_sink;
size_t images;
size_t direct_fallbacks;
size_t cached_wing_clipped;
size_t discarded_below_min_y;
} TriangleSplatContext;
typedef struct {
size_t images;
size_t direct_fallbacks;
size_t cached_wing_clipped;
size_t discarded_below_min_y;
size_t gpu_event_count, gpu_batch_count, gpu_timed_batch_count;
double gpu_upload_seconds, gpu_kernel_seconds, gpu_download_seconds;
int failed;
#ifdef GR_DEBUG
double max_raw_magnification;
size_t magnification_clamped_triangles;
#endif
} CatalogSplatStats;
static void copy_psf_splat_stats(PsfSplatStats *destination,
CatalogSplatStats source)
{
if (destination == NULL)
return;
*destination = (PsfSplatStats){.cached_splats =
source.images - source.direct_fallbacks -
source.discarded_below_min_y,
.cached_wing_clipped = source.cached_wing_clipped,
.direct_fallbacks = source.direct_fallbacks,
.discarded_below_min_y = source.discarded_below_min_y,
.gpu_event_count = source.gpu_event_count,
.gpu_batch_count = source.gpu_batch_count,
.gpu_timed_batch_count = source.gpu_timed_batch_count,
.gpu_upload_seconds = source.gpu_upload_seconds,
.gpu_kernel_seconds = source.gpu_kernel_seconds,
.gpu_download_seconds = source.gpu_download_seconds,
#ifdef GR_DEBUG
.max_raw_magnification = source.max_raw_magnification,
.magnification_clamped_triangles =
source.magnification_clamped_triangles,
#endif
};
}
static int splat_catalog_tile(const Star *stars, size_t count,
int fully_contained, void *opaque) {
TriangleSplatContext *context = opaque;
for (size_t s = 0; s < count; ++s) {
const Star *star = &stars[s];
double weights[3];
if (!fully_contained &&
spherical_barycentric(star->direction, context->vertex[0]->n_infinity,
context->vertex[1]->n_infinity,
context->vertex[2]->n_infinity, weights))
continue;
if (fully_contained) {
/* Only inverse-map weights remain: no per-star containment test. */
if (spherical_barycentric_weights(star->direction,
context->vertex[0]->n_infinity,
context->vertex[1]->n_infinity,
context->vertex[2]->n_infinity, weights))
continue;
}
if (!owns_source_boundary(context->triangle, weights))
continue;
const double image_x = weights[0] * context->vertex[0]->image_x +
weights[1] * context->vertex[1]->image_x +
weights[2] * context->vertex[2]->image_x;
const double image_y = weights[0] * context->vertex[0]->image_y +
weights[1] * context->vertex[1]->image_y +
weights[2] * context->vertex[2]->image_y;
const double log_g = weights[0] * context->vertex[0]->log_frequency_ratio +
weights[1] * context->vertex[1]->log_frequency_ratio +
weights[2] * context->vertex[2]->log_frequency_ratio;
const LinearRgb color = blackbody_to_linear_rgb(star->temperature_K * exp(log_g));
const double flux = context->exposure * star->amplitude * context->magnification;
PsfCachedEvent event;
const int direct_fallback = psf_prepare_cached_event(
&event, image_x, image_y, color, flux, context->psf, context->psf_cache,
context->max_cache_psf_flux, context->psf_relative_tail, context->psf_min_y);
if (direct_fallback == 1) {
/* Exclude every other submit and fallback until HDR is back on device. */
#ifdef PSF_BACKEND_HIP
const double submit_start = omp_get_wtime();
if (context->event_sink->hip_lock)
omp_set_lock(context->event_sink->hip_lock);
#endif
int failed = psf_event_sink_finish_for_cpu(context->event_sink);
if (!failed) {
splat_moffat_direct(context->hdr, context->width, context->height, image_x,
image_y, color, flux, context->psf,
context->psf_relative_tail, context->psf_min_y);
failed = psf_event_sink_resume_gpu(context->event_sink);
}
#ifdef PSF_BACKEND_HIP
if (context->event_sink->hip_lock)
omp_unset_lock(context->event_sink->hip_lock);
context->event_sink->submit_seconds += omp_get_wtime() - submit_start;
#endif
if (failed) return -1;
} else if (direct_fallback != 3) {
#if FRAME_PSF_EVENT_SINK || defined(PSF_BACKEND_HIP)
psf_event_sink_emit(context->event_sink, &event);
#else
splat_prepared_cached_event(context->hdr, context->width, context->height,
&event, context->psf_cache);
#endif
}
if (context->event_sink->failed)
return -1;
context->direct_fallbacks += direct_fallback == 1;
context->cached_wing_clipped += direct_fallback == 2;
context->discarded_below_min_y += direct_fallback == 3;
#ifdef GR_DEBUG
if (direct_fallback == 1)
/* This is deliberately emitted by the active splat worker: a direct
* fallback can be the long-running work a Debug render is waiting on.
* Do not add a critical section here; interleaved Debug lines are more
* useful than stalling the other workers. */
fprintf(stderr,
"Debug: %.3f s splat worker %d triangle %zu star %zu uses direct "
"PSF fallback (image %.3f, %.3f; flux %.6g).\n",
omp_get_wtime(), omp_get_thread_num(), context->triangle_index,
s, image_x, image_y, flux);
#endif
++context->images;
}
return 0;
}
static CatalogSplatStats splat_catalog_triangles(
const FrameLensMesh *mesh, StarCatalog *catalog, double *hdr, int width,
int height, double exposure, const PointSpreadFunction *psf,
const PsfKernelCache *psf_cache, double max_magnification,
double max_cache_psf_flux, double psf_relative_tail,
double psf_min_y,
size_t first_triangle, size_t last_triangle, PsfEventSink *event_sink) {
CatalogSplatStats stats = {0};
PsfEventSink owned_sink;
const int owns_sink = event_sink == NULL;
if (owns_sink) {
if (psf_event_sink_init(&owned_sink, hdr, width, height, psf_cache))
return (CatalogSplatStats){.failed = 1};
event_sink = &owned_sink;
}
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 raw_magnification = image_area / source_area;
double magnification = raw_magnification;
#ifdef GR_DEBUG
stats.max_raw_magnification = fmax(stats.max_raw_magnification,
raw_magnification);
#endif
if (raw_magnification > max_magnification) {
magnification = max_magnification;
#ifdef GR_DEBUG
++stats.magnification_clamped_triangles;
#endif
}
const double direction[3][3] = {
{vertex[0]->n_infinity[0], vertex[0]->n_infinity[1], vertex[0]->n_infinity[2]},
{vertex[1]->n_infinity[0], vertex[1]->n_infinity[1], vertex[1]->n_infinity[2]},
{vertex[2]->n_infinity[0], vertex[2]->n_infinity[1], vertex[2]->n_infinity[2]}};
TriangleSplatContext context = {.vertex = {vertex[0], vertex[1], vertex[2]},
.triangle = &mesh->triangles[t],
.triangle_index = t, .hdr = hdr,
.width = width, .height = height,
.exposure = exposure, .magnification = magnification,
.psf = psf, .psf_cache = psf_cache,
.max_cache_psf_flux = max_cache_psf_flux,
.psf_relative_tail = psf_relative_tail,
.psf_min_y = psf_min_y,
.event_sink = event_sink};
const int visit_result = catalog_visit_source_triangle(
catalog, direction, 0, splat_catalog_tile, &context);
if (visit_result == 0)
stats.images += context.images;
else {
stats.failed = event_sink->failed;
if (stats.failed)
break;
}
stats.direct_fallbacks += context.direct_fallbacks;
stats.cached_wing_clipped += context.cached_wing_clipped;
stats.discarded_below_min_y += context.discarded_below_min_y;
}
if (owns_sink) {
if (psf_event_sink_destroy(&owned_sink))
stats.failed = 1;
#ifdef PSF_BACKEND_HIP
stats.gpu_event_count = owned_sink.hip_timing.event_count;
stats.gpu_batch_count = owned_sink.hip_timing.batch_count;
stats.gpu_timed_batch_count = owned_sink.hip_timing.timed_batch_count;
stats.gpu_upload_seconds = owned_sink.hip_timing.upload_seconds;
stats.gpu_kernel_seconds = owned_sink.hip_timing.kernel_seconds;
stats.gpu_download_seconds = owned_sink.hip_timing.download_seconds;
#endif
}
return stats;
}
static void prefetch_catalog_for_mesh(const FrameLensMesh *mesh,
StarCatalog *catalog,
int worker_count,
CatalogPrefetchStats *stats) {
if (catalog->kind != STAR_CATALOG_ALL_SKY)
return;
unsigned char requested[CATALOG_ALL_SKY_TILE_COUNT] = {0};
for (size_t t = 0; t < mesh->triangle_count; ++t) {
const LensVertex *vertex[3];
if (!usable_triangle(mesh, &mesh->triangles[t], vertex))
continue;
const double direction[3][3] = {
{vertex[0]->n_infinity[0], vertex[0]->n_infinity[1], vertex[0]->n_infinity[2]},
{vertex[1]->n_infinity[0], vertex[1]->n_infinity[1], vertex[1]->n_infinity[2]},
{vertex[2]->n_infinity[0], vertex[2]->n_infinity[1], vertex[2]->n_infinity[2]}};
(void)catalog_mark_source_triangle_tiles(direction, requested);
}
(void)catalog_prefetch_marked_tiles(catalog, requested, worker_count, stats);
}
size_t frame_splat_catalog(const FrameLensMesh *mesh,
StarCatalog *catalog, double *hdr, int width,
int height, double exposure,
const PointSpreadFunction *psf,
const PsfKernelCache *psf_cache,
double max_magnification,
double max_cache_psf_flux,
double psf_relative_tail,
double psf_min_y,
int limit_workers_by_memory,
int catalog_load_workers,
CatalogPrefetchStats *prefetch_stats,
PsfSplatStats *psf_stats,
const FrameSplatProgress *progress) {
if (mesh == NULL || catalog == NULL || hdr == NULL || exposure <= 0.0 ||
psf == NULL || width <= 0 || height <= 0 || catalog_load_workers <= 0 ||
isnan(max_magnification) || max_magnification <= 0.0 ||
!isfinite(max_cache_psf_flux) || max_cache_psf_flux < 1.0 ||
!isfinite(psf_relative_tail) || psf_relative_tail <= 0.0 ||
psf_relative_tail >= 1.0 || !isfinite(psf_min_y) || psf_min_y < 0.0)
return 0;
/* A bounded parallel read phase completes before splatting. Its serial cache
* commit leaves immutable tile data for the OpenMP splat workers. */
if (progress != NULL && progress->callback != NULL)
progress->callback(progress->context, FRAME_SPLAT_PROGRESS_PREFETCH_BEGIN,
0, mesh->triangle_count);
prefetch_catalog_for_mesh(mesh, catalog, catalog_load_workers, prefetch_stats);
if (progress != NULL && progress->callback != NULL)
progress->callback(progress->context, FRAME_SPLAT_PROGRESS_PREFETCH_END,
prefetch_stats == NULL ? 0 : prefetch_stats->requested_tiles,
prefetch_stats == NULL ? 0 : prefetch_stats->requested_tiles);
if (psf_stats != NULL)
*psf_stats = (PsfSplatStats){0};
if (progress != NULL && progress->callback != NULL)
progress->callback(progress->context, FRAME_SPLAT_PROGRESS_BEGIN, 0,
mesh->triangle_count);
#ifdef PSF_BACKEND_HIP
/* CPU workers own only bounded event chunks. A single device cache/HDR is
* borrowed under a coarse submission lock, never duplicated per worker. */
PsfEventSink owner;
if (psf_event_sink_init(&owner, hdr, width, height, psf_cache)) return SIZE_MAX;
omp_lock_t submit_lock;
omp_init_lock(&submit_lock);
size_t hip_images = 0, hip_direct = 0, hip_clipped = 0, hip_discarded = 0;
int hip_failed = 0, hip_workers = 0;
double hip_generate_seconds = 0.0, hip_submit_seconds = 0.0;
const double hip_start = omp_get_wtime();
#ifdef GR_DEBUG
double hip_max_magnification = 0.0;
size_t hip_clamped = 0;
#pragma omp parallel reduction(+ : hip_images, hip_direct, hip_clipped, hip_discarded, hip_clamped, hip_generate_seconds, hip_submit_seconds) reduction(max : hip_failed, hip_max_magnification)
#else
#pragma omp parallel reduction(+ : hip_images, hip_direct, hip_clipped, hip_discarded, hip_generate_seconds, hip_submit_seconds) reduction(max : hip_failed)
#endif
{
const double worker_start = omp_get_wtime();
#pragma omp single
hip_workers = omp_get_num_threads();
PsfEventSink local = {.hdr = hdr, .width = width, .height = height,
.cache = psf_cache, .hip = owner.hip, .hip_lock = &submit_lock,
.borrowed_hip = 1};
local.events = malloc(PSF_EVENT_SINK_CAPACITY * sizeof *local.events);
if (!local.events) {
snprintf(local.hip_message, sizeof local.hip_message, "worker event allocation failed");
psf_event_sink_mark_failed(&local, "producer initialization");
}
const size_t worker = (size_t)omp_get_thread_num();
const size_t workers = (size_t)omp_get_num_threads();
size_t completed = 0, next_report = 8;
if (progress && progress->worker_callback)
progress->worker_callback(progress->context, worker, workers, 0, 0);
#pragma omp for schedule(dynamic, 1) nowait
for (size_t t = 0; t < mesh->triangle_count; ++t) {
if (local.failed) continue;
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y,
t, t + 1, &local);
hip_images += stats.images;
hip_direct += stats.direct_fallbacks;
hip_clipped += stats.cached_wing_clipped;
hip_discarded += stats.discarded_below_min_y;
hip_failed |= stats.failed;
#ifdef GR_DEBUG
hip_max_magnification = fmax(hip_max_magnification, stats.max_raw_magnification);
hip_clamped += stats.magnification_clamped_triangles;
#endif
++completed;
if (progress && progress->worker_callback && completed == next_report) {
progress->worker_callback(progress->context, worker, workers, completed, 0);
if (next_report <= SIZE_MAX / 2) next_report *= 2;
}
}
if (psf_event_sink_destroy(&local)) hip_failed = 1;
hip_submit_seconds += local.submit_seconds;
hip_generate_seconds += omp_get_wtime() - worker_start - local.submit_seconds;
if (progress && progress->worker_callback)
progress->worker_callback(progress->context, worker, workers, completed, 1);
}
omp_destroy_lock(&submit_lock);
if (psf_event_sink_destroy(&owner)) hip_failed = 1;
fprintf(stderr, "HIP producers: %d workers; summed generation %.3f s, submission/fallback %.3f s; wall %.3f s\n",
hip_workers, hip_generate_seconds, hip_submit_seconds, omp_get_wtime() - hip_start);
copy_psf_splat_stats(psf_stats, (CatalogSplatStats){
.images = hip_images, .direct_fallbacks = hip_direct,
.cached_wing_clipped = hip_clipped, .discarded_below_min_y = hip_discarded,
.gpu_event_count = owner.hip_timing.event_count,
.gpu_batch_count = owner.hip_timing.batch_count,
.gpu_timed_batch_count = owner.hip_timing.timed_batch_count,
.gpu_upload_seconds = owner.hip_timing.upload_seconds,
.gpu_kernel_seconds = owner.hip_timing.kernel_seconds,
.gpu_download_seconds = owner.hip_timing.download_seconds,
#ifdef GR_DEBUG
.max_raw_magnification = hip_max_magnification,
.magnification_clamped_triangles = hip_clamped,
#endif
});
if (hip_failed) return SIZE_MAX;
if (progress != NULL && progress->callback != NULL)
progress->callback(progress->context, FRAME_SPLAT_PROGRESS_END,
mesh->triangle_count, mesh->triangle_count);
return hip_images;
#endif
const size_t pixel_count = (size_t)width * height * 3;
if (pixel_count > SIZE_MAX / sizeof(double))
{
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y,
0, mesh->triangle_count, NULL);
copy_psf_splat_stats(psf_stats, stats);
return stats.images;
}
const size_t buffer_bytes = pixel_count * sizeof(double);
const int max_threads = omp_get_max_threads();
size_t worker_count = max_threads > 0 ? (size_t)max_threads : 0;
if (limit_workers_by_memory) {
worker_count = FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES / buffer_bytes;
if (worker_count > (size_t)max_threads)
worker_count = (size_t)max_threads;
}
if (worker_count < 2 || worker_count > INT_MAX)
{
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y,
0, mesh->triangle_count, NULL);
copy_psf_splat_stats(psf_stats, stats);
return stats.images;
}
double **private_hdr = calloc(worker_count, sizeof *private_hdr);
if (private_hdr == NULL)
{
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y,
0, mesh->triangle_count, NULL);
copy_psf_splat_stats(psf_stats, stats);
return stats.images;
}
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);
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y,
0, mesh->triangle_count, NULL);
copy_psf_splat_stats(psf_stats, stats);
return stats.images;
}
size_t images = 0, direct_fallbacks = 0, cached_wing_clipped = 0,
discarded_below_min_y = 0;
#ifdef GR_DEBUG
double max_raw_magnification = 0.0;
size_t magnification_clamped_triangles = 0;
#endif
/* Keep the ordinary render loop byte-for-byte free of progress checks. */
if (progress != NULL && progress->worker_callback != NULL) {
#ifdef GR_DEBUG
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y, magnification_clamped_triangles) reduction(max : max_raw_magnification)
#else
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y)
#endif
{
const size_t worker = (size_t)omp_get_thread_num();
PsfEventSink event_sink;
psf_event_sink_init(&event_sink, private_hdr[worker], width, height, psf_cache);
size_t local_triangles = 0, next_report = 8;
progress->worker_callback(progress->context, worker, worker_count, 0, 0);
/* Source density and lens magnification can vary by orders of magnitude
* between neighboring image triangles. Dynamic single-triangle chunks
* prevent a small sky region from leaving the other private HDR workers
* idle. Each worker still owns its HDR buffer exclusively. */
#pragma omp for schedule(dynamic, 1)
for (size_t triangle = 0; triangle < mesh->triangle_count; ++triangle) {
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, private_hdr[worker], width, height, exposure, psf,
psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail,
psf_min_y,
triangle, triangle + 1, &event_sink);
images += stats.images;
direct_fallbacks += stats.direct_fallbacks;
cached_wing_clipped += stats.cached_wing_clipped;
discarded_below_min_y += stats.discarded_below_min_y;
#ifdef GR_DEBUG
max_raw_magnification = fmax(max_raw_magnification, stats.max_raw_magnification);
magnification_clamped_triangles += stats.magnification_clamped_triangles;
#endif
++local_triangles;
if (local_triangles == next_report) {
progress->worker_callback(progress->context, worker, worker_count,
local_triangles, 0);
if (next_report <= SIZE_MAX / 2)
next_report *= 2;
}
}
psf_event_sink_destroy(&event_sink);
progress->worker_callback(progress->context, worker, worker_count,
local_triangles, 1);
}
} else {
#ifdef GR_DEBUG
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y, magnification_clamped_triangles) reduction(max : max_raw_magnification)
#else
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y)
#endif
{
const size_t worker = (size_t)omp_get_thread_num();
PsfEventSink event_sink;
psf_event_sink_init(&event_sink, private_hdr[worker], width, height, psf_cache);
#pragma omp for schedule(dynamic, 1)
for (size_t triangle = 0; triangle < mesh->triangle_count; ++triangle) {
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, private_hdr[worker], width, height, exposure, psf,
psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail,
psf_min_y,
triangle, triangle + 1, &event_sink);
images += stats.images;
direct_fallbacks += stats.direct_fallbacks;
cached_wing_clipped += stats.cached_wing_clipped;
discarded_below_min_y += stats.discarded_below_min_y;
#ifdef GR_DEBUG
max_raw_magnification = fmax(max_raw_magnification, stats.max_raw_magnification);
magnification_clamped_triangles += stats.magnification_clamped_triangles;
#endif
}
psf_event_sink_destroy(&event_sink);
}
}
#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);
copy_psf_splat_stats(psf_stats, (CatalogSplatStats){
.images = images,
.direct_fallbacks = direct_fallbacks,
.cached_wing_clipped = cached_wing_clipped,
.discarded_below_min_y = discarded_below_min_y,
#ifdef GR_DEBUG
.max_raw_magnification = max_raw_magnification,
.magnification_clamped_triangles =
magnification_clamped_triangles,
#endif
});
if (progress != NULL && progress->callback != NULL)
progress->callback(progress->context, FRAME_SPLAT_PROGRESS_END,
mesh->triangle_count, mesh->triangle_count);
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);
free(mesh->samples);
free(mesh->probe_slots);
*mesh = (FrameLensMesh){0};
}
Reference in new issue View Git Blame Copy Permalink