wyj/GR-raytracing
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
d02ecb2181321c51b52a81b95b33fd02a6165b2a
GR-raytracing/src/frame.c
T
wyj d02ecb2181 Feat: Add dummy PSF chunk diagnostics 
Reuse production catalog mapping and per-worker chunk boundaries to report triangle and spatial distributions without initializing HIP or writing image outputs. Include a protected full-catalog runner and build documentation.
2026-09-11 23:27:42 -04:00

1810 lines
74 KiB
C
Raw Blame History

#include "frame.h"
#ifdef PSF_BACKEND_HIP
#include "hip_psf.h"
#endif
#ifdef PSF_BACKEND_DUMMY
#include "dummy_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
#ifdef PSF_BACKEND_DUMMY
DummyPsfSink *dummy;
DummyPsfChunk *dummy_chunk;
int borrowed_dummy;
#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) {
#ifdef PSF_BACKEND_DUMMY
if (dummy_psf_sink_create(&sink->dummy, width, height,
PSF_EVENT_SINK_CAPACITY)) {
fputs("Dummy PSF backend initialization failed\n", stderr);
dummy_psf_sink_destroy(sink->dummy);
sink->dummy = NULL;
sink->failed = 1;
return -1;
}
#else
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
#endif
}
return 0;
}
static int psf_event_sink_flush_unlocked(PsfEventSink *sink) {
if (sink->failed)
return -1;
#ifdef PSF_BACKEND_DUMMY
if (sink->dummy_chunk != NULL)
return dummy_psf_chunk_flush(sink->dummy_chunk);
#endif
#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;
}
#ifndef PSF_BACKEND_DUMMY
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;
}
#endif
static int psf_event_sink_destroy(PsfEventSink *sink) {
#ifdef PSF_BACKEND_DUMMY
int dummy_result = sink->dummy_chunk != NULL
? dummy_psf_chunk_flush(sink->dummy_chunk)
: 0;
dummy_psf_chunk_destroy(sink->dummy_chunk);
sink->dummy_chunk = NULL;
if (!sink->borrowed_dummy) {
dummy_psf_sink_report(sink->dummy);
dummy_psf_sink_destroy(sink->dummy);
sink->dummy = NULL;
}
return dummy_result;
#endif
#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) || defined(PSF_BACKEND_DUMMY)
static void psf_event_sink_emit(PsfEventSink *sink, const PsfCachedEvent *event,
size_t triangle, double triangle_center_x,
double triangle_center_y) {
if (sink->failed)
return;
#ifdef PSF_BACKEND_DUMMY
if (dummy_psf_chunk_emit(sink->dummy_chunk, event, triangle,
triangle_center_x, triangle_center_y))
sink->failed = 1;
return;
#else
(void)triangle;
(void)triangle_center_x;
(void)triangle_center_y;
#endif
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);
#ifdef FRAME_PSF_DIAGNOSTIC
/* Test-only consumer: normal query/mapping/colour/classification above.
* Never produces an HDR image; not compiled into renderer binaries. */
if (frame_psf_diagnostic_visit(&event, direct_fallback,
context->triangle_index, image_x, image_y,
color, flux))
return -1;
context->direct_fallbacks += direct_fallback == 1;
context->cached_wing_clipped += direct_fallback == 2;
context->discarded_below_min_y += direct_fallback == 3;
++context->images;
continue;
#endif
if (direct_fallback == 1) {
/* Exclude every other submit and fallback until HDR is back on device. */
#ifdef PSF_BACKEND_DUMMY
if (psf_event_sink_flush(context->event_sink)) return -1;
#else
#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;
#endif
} else if (direct_fallback != 3) {
#if FRAME_PSF_EVENT_SINK || defined(PSF_BACKEND_HIP) || defined(PSF_BACKEND_DUMMY)
const double triangle_center_x =
(context->vertex[0]->image_x + context->vertex[1]->image_x +
context->vertex[2]->image_x) / 3.0;
const double triangle_center_y =
(context->vertex[0]->image_y + context->vertex[1]->image_y +
context->vertex[2]->image_y) / 3.0;
psf_event_sink_emit(context->event_sink, &event, context->triangle_index,
triangle_center_x, triangle_center_y);
#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) {
#ifdef FRAME_PSF_DIAGNOSTIC
if (frame_psf_diagnostic_stopped()) break;
#endif
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;
#ifdef FRAME_PSF_DIAGNOSTIC
if (!frame_psf_diagnostic_select(raw_magnification)) continue;
#endif
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_DUMMY
PsfEventSink owner;
if (psf_event_sink_init(&owner, hdr, width, height, psf_cache)) return SIZE_MAX;
size_t dummy_images = 0, dummy_direct = 0, dummy_clipped = 0,
dummy_discarded = 0;
int dummy_failed = 0, dummy_workers = 0;
const double dummy_start = omp_get_wtime();
#ifdef GR_DEBUG
double dummy_max_magnification = 0.0;
size_t dummy_clamped = 0;
#pragma omp parallel reduction(+ : dummy_images, dummy_direct, dummy_clipped, dummy_discarded, dummy_clamped) reduction(max : dummy_failed, dummy_max_magnification)
#else
#pragma omp parallel reduction(+ : dummy_images, dummy_direct, dummy_clipped, dummy_discarded) reduction(max : dummy_failed)
#endif
{
#pragma omp single
dummy_workers = omp_get_num_threads();
const size_t worker = (size_t)omp_get_thread_num();
PsfEventSink local = {.hdr = hdr, .width = width, .height = height,
.cache = psf_cache, .dummy = owner.dummy, .borrowed_dummy = 1};
if (dummy_psf_chunk_create(&local.dummy_chunk, owner.dummy, worker)) {
fputs("Dummy PSF worker chunk initialization failed\n", stderr);
local.failed = 1;
}
size_t completed = 0, next_report = 8;
if (progress && progress->worker_callback)
progress->worker_callback(progress->context, worker,
(size_t)dummy_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);
dummy_images += stats.images;
dummy_direct += stats.direct_fallbacks;
dummy_clipped += stats.cached_wing_clipped;
dummy_discarded += stats.discarded_below_min_y;
dummy_failed |= stats.failed;
#ifdef GR_DEBUG
dummy_max_magnification = fmax(dummy_max_magnification,
stats.max_raw_magnification);
dummy_clamped += stats.magnification_clamped_triangles;
#endif
++completed;
if (progress && progress->worker_callback && completed == next_report) {
progress->worker_callback(progress->context, worker,
(size_t)dummy_workers, completed, 0);
if (next_report <= SIZE_MAX / 2) next_report *= 2;
}
}
if (psf_event_sink_destroy(&local)) dummy_failed = 1;
if (progress && progress->worker_callback)
progress->worker_callback(progress->context, worker,
(size_t)dummy_workers, completed, 1);
}
if (psf_event_sink_destroy(&owner)) dummy_failed = 1;
fprintf(stderr,
"Dummy PSF producers: %d workers; classification/chunk wall %.3f s\n",
dummy_workers, omp_get_wtime() - dummy_start);
copy_psf_splat_stats(psf_stats, (CatalogSplatStats){
.images = dummy_images, .direct_fallbacks = dummy_direct,
.cached_wing_clipped = dummy_clipped,
.discarded_below_min_y = dummy_discarded,
#ifdef GR_DEBUG
.max_raw_magnification = dummy_max_magnification,
.magnification_clamped_triangles = dummy_clamped,
#endif
});
if (dummy_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 dummy_images;
#endif
#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