#include "frame.h" #include "optics.h" #include #include #include #include #include #include /* 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; } typedef struct { size_t a, b, triangle; unsigned int side; size_t midpoint; } 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 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) return 0; /* 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 side = longest_side(mesh, triangle); 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 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 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 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; } 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); if (edges == NULL || requested == NULL || allowed == NULL || parity == NULL || jacobians == NULL) { free(edges); free(requested); free(allowed); free(parity); free(jacobians); return -1; } 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, SIZE_MAX}; } if (allowed[i]) { 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; } /* Two requested sides require red refinement. Add the third side, then * close the new shared edge requests before allocating any vertices. */ for (;;) { int changed = 0; for (size_t t = 0; t < mesh->triangle_count; ++t) { unsigned int count = 0; for (unsigned int side = 0; side < 3; ++side) count += requested[3 * t + side] != 0; if (count >= 2 && allowed[t]) for (unsigned int side = 0; side < 3; ++side) if (!requested[3 * t + side]) { requested[3 * t + side] = 1; changed = 1; } } 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) { any |= requested[3 * edges[i].triangle + edges[i].side] != 0; possible &= allowed[edges[i].triangle] != 0; } if (any && possible) for (size_t i = first; i < last; ++i) if (!requested[3 * edges[i].triangle + edges[i].side]) { requested[3 * edges[i].triangle + edges[i].side] = 1; changed = 1; } first = last; } if (!changed) break; } 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); 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); 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) edges[i].midpoint = 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); 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; size_t a = triangle->vertex[side], b = triangle->vertex[(side + 1) % 3]; if (a > b) { size_t swap = a; a = b; b = swap; } for (size_t e = 0; e < edge_count; ++e) if (edges[e].a == a && edges[e].b == b) { middle[side] = edges[e].midpoint; break; } } 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 { /* Two or three requested edges become conforming red refinement. */ size_t ab = middle[0], bc = middle[1], ca = middle[2]; if (ab == SIZE_MAX || bc == SIZE_MAX || ca == SIZE_MAX) { /* This cannot be made conforming from one-probe-per-triangle data. */ 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); return (int)midpoint_count; } int frame_lens_mesh_refine(FrameLensMesh *mesh, const SpacetimeSource *spacetime, const ObserverState *observer, const GeodesicTraceConfig *trace, const RefinementConfig *config) { if (mesh == NULL || spacetime == NULL || observer == NULL || trace == NULL || config == NULL) return -1; for (;;) { const int requested = frame_lens_mesh_prepare_generation(mesh, config); if (requested < 0) return -1; if (requested == 0) return 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); } if (frame_lens_mesh_finish_generation(mesh, config) < 0) return -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; 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; double *hdr; int width, height; double exposure, magnification; const PointSpreadFunction *psf; const PsfKernelCache *psf_cache; size_t images; size_t direct_fallbacks; } TriangleSplatContext; 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)) return -1; } 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)); context->direct_fallbacks += splat_moffat_cached( context->hdr, context->width, context->height, image_x, image_y, color, context->exposure * star->amplitude * context->magnification, context->psf, context->psf_cache); ++context->images; } return 0; } static size_t splat_catalog_triangles(const FrameLensMesh *mesh, StarCatalog *catalog, double *hdr, int width, int height, double exposure, const PointSpreadFunction *psf, const PsfKernelCache *psf_cache, size_t first_triangle, size_t last_triangle, size_t *direct_fallbacks) { size_t images = 0; for (size_t t = first_triangle; t < last_triangle; ++t) { const LensVertex *vertex[3]; if (!usable_triangle(mesh, &mesh->triangles[t], vertex)) continue; const double source_area = spherical_area( vertex[0]->n_infinity, vertex[1]->n_infinity, vertex[2]->n_infinity); const double image_area = spherical_area(vertex[0]->camera_direction, vertex[1]->camera_direction, vertex[2]->camera_direction); const double magnification = image_area / source_area; 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], .hdr = hdr, .width = width, .height = height, .exposure = exposure, .magnification = magnification, .psf = psf, .psf_cache = psf_cache}; if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile, &context) == 0) images += context.images; *direct_fallbacks += context.direct_fallbacks; } return images; } 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, 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) 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); const size_t pixel_count = (size_t)width * height * 3; if (pixel_count > SIZE_MAX / sizeof(double)) { size_t fallbacks = 0; const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure, psf, psf_cache, 0, mesh->triangle_count, &fallbacks); if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks}; return 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) { size_t fallbacks = 0; const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure, psf, psf_cache, 0, mesh->triangle_count, &fallbacks); if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks}; return images; } double **private_hdr = calloc(worker_count, sizeof *private_hdr); if (private_hdr == NULL) { size_t fallbacks = 0; const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure, psf, psf_cache, 0, mesh->triangle_count, &fallbacks); if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks}; return 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); size_t fallbacks = 0; const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure, psf, psf_cache, 0, mesh->triangle_count, &fallbacks); if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks}; return images; } size_t images = 0, direct_fallbacks = 0; #pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks) { const size_t worker = (size_t)omp_get_thread_num(); /* 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) images += splat_catalog_triangles(mesh, catalog, private_hdr[worker], width, height, exposure, psf, psf_cache, triangle, triangle + 1, &direct_fallbacks); } #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); if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - direct_fallbacks, direct_fallbacks}; 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}; }