From 48dcf4e7077faad9318a3a6609ddc235d76823d1 Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Fri, 28 Aug 2026 23:06:53 -0400 Subject: [PATCH] Frame: add adaptive mesh refinement --- README.md | 52 ++++ src/frame.c | 525 ++++++++++++++++++++++++++++++++++++- src/frame.h | 45 +++- src/main.c | 162 +++++++----- tests/test_frame.c | 32 +++ tests/test_schwarzschild.c | 20 ++ video_rendering_plan.md | 4 +- 7 files changed, 768 insertions(+), 72 deletions(-) diff --git a/README.md b/README.md index e67cf57..461be75 100644 --- a/README.md +++ b/README.md @@ -142,6 +142,58 @@ metric or its infinity criterion; frame, observer, and integrator use only as `--coarse-cell-pixels`; it is a Phase-0 sampling knob, not a settled production refinement threshold. +### Adaptive image mesh refinement + +Adaptive refinement is disabled by default (`--refine-max-level 0`), so the +existing coarse-mesh renders remain unchanged. When enabled, its defaults are +an absolute direction error of `1e-3` degrees, relative error `0.1`, minimum +long edge `0.5` pixels, and minimum area `0.25` pixel-squared. Each value can +be overridden independently: + +```text +--refine-max-level N +--refine-angle-abs-deg D +--refine-angle-rel R +--refine-min-edge-pixels P +--refine-min-area-pixels2 A +``` + +`N` caps the triangle refinement level. Let `e` be the angle between the +traced longest-edge midpoint direction and the normalized endpoint +interpolation, and let `s` be the angle between those two endpoint **camera +directions**. Both are evaluated internally in radians; the absolute CLI +threshold `D` is specified in degrees and converted before comparison. `s` +is the angular geometric size of the image triangle's test edge, not a +source-sky/lens-map length. A locally escaped triangle is split only when +**both** `e > D_rad` (the converted `--refine-angle-abs-deg D`) and +`e / max(s, 1e-15) > --refine-angle-rel`. `P` and `A` +prevent selecting a leaf already at or below the requested image-plane +long-edge and area scales. +Triangles whose three vertices disagree between capture and escape are split +independently of the direction-error thresholds, allowing the mesh to follow a +shadow boundary. This first implementation deliberately does not evaluate +orientation or Jacobian criteria. + +For a short Schwarzschild diagnostic that permits at most one actual split +generation, for example: + +```sh +./build/schwarzschild_sky --catalog assets/sky_grid_5deg.csv \ + --width 48 --height 48 --coarse-cell-pixels 24 --fov-deg 40 \ + --refine-max-level 1 --refine-angle-abs-deg 0.001 \ + --refine-angle-rel 0.001 --refine-min-edge-pixels 1 \ + --refine-min-area-pixels2 1 --draw-mesh \ + --output output/imgs/schwarzschild_refinement.png +``` + +For movies, each refinement generation completes the full newest-to-oldest +time-slab sweep before any probe becomes a mesh vertex. Newly added vertices +are therefore traced only by the next generation; the renderer never returns +to a slab that has already been released. At the start of every generation, +newly inserted vertices and geometry-only longest-edge probes for its new +leaves are collected together, so both ray sets use the same parallel +`RayPool` pass. + `--look-ra-deg` and `--look-dec-deg` rotate that fixed tetrad so its forward axis is the corresponding catalog direction; their defaults reproduce the original `-Z` view. `--exposure` converts a catalog's physical flux diff --git a/src/frame.c b/src/frame.c index df4c515..26ed786 100644 --- a/src/frame.c +++ b/src/frame.c @@ -7,6 +7,7 @@ #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. @@ -76,14 +77,16 @@ int frame_lens_mesh_build_coarse(FrameLensMesh *mesh, int width, int height, 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}}; + (LensTriangle){{top_left, bottom_left, bottom_right}, 0, 0}; triangles[next_triangle++] = - (LensTriangle){{top_left, bottom_right, top_right}}; + (LensTriangle){{top_left, bottom_right, top_right}, 0, 0}; } *mesh = (FrameLensMesh){.vertices = vertices, .triangles = triangles, .vertex_count = vertex_count, - .triangle_count = triangle_count}; + .vertex_capacity = vertex_count, + .triangle_count = triangle_count, + .triangle_capacity = triangle_count}; return 0; } @@ -101,6 +104,7 @@ int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime, 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]; @@ -110,6 +114,519 @@ int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime, 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 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); + if (edges == NULL || requested == NULL || allowed == NULL) { + free(edges); free(requested); free(allowed); + 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); + } + } + qsort(edges, edge_count, sizeof *edges, compare_mesh_edge); + /* 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); + 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); + 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); + 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); + 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]; @@ -500,5 +1017,7 @@ void frame_lens_mesh_destroy(FrameLensMesh *mesh) { return; free(mesh->vertices); free(mesh->triangles); + free(mesh->samples); + free(mesh->probe_slots); *mesh = (FrameLensMesh){0}; } diff --git a/src/frame.h b/src/frame.h index d9ce27e..6c3d499 100644 --- a/src/frame.h +++ b/src/frame.h @@ -15,16 +15,42 @@ typedef struct { double n_infinity[3]; double log_frequency_ratio; RayEndpointStatus status; + int traced; } LensVertex; typedef struct { size_t vertex[3]; + unsigned int level; + int evaluated; } LensTriangle; +typedef struct { + unsigned int max_level; + double angle_absolute_rad; + double angle_relative; + double min_edge_pixels; + double min_area_pixels2; +} RefinementConfig; + +typedef enum { FRAME_SAMPLE_VERTEX, FRAME_SAMPLE_PROBE } FrameSampleKind; + +typedef struct { + FrameSampleKind kind; + size_t vertex_id; + size_t edge_vertex[2]; + LensVertex vertex; +} FrameSample; + typedef struct { LensVertex *vertices; LensTriangle *triangles; - size_t vertex_count, triangle_count; + size_t vertex_count, vertex_capacity; + size_t triangle_count, triangle_capacity; + FrameSample *samples; + size_t sample_count, sample_capacity; + int samples_include_probes; + size_t *probe_slots; + size_t probe_slot_capacity; } FrameLensMesh; typedef enum { @@ -48,6 +74,23 @@ int frame_lens_mesh_build_coarse(FrameLensMesh *mesh, int width, int height, int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime, const ObserverState *observer, const GeodesicTraceConfig *trace); +/* Builds exactly one generation of requests. Probe results must be installed + * only after the caller has completed the generation's tracing sweep. */ +int frame_lens_mesh_prepare_generation(FrameLensMesh *mesh, + const RefinementConfig *config); +const FrameSample *frame_lens_mesh_samples(const FrameLensMesh *mesh, + size_t *count); +int frame_lens_mesh_install_sample(FrameLensMesh *mesh, size_t sample_id, + const RayEndpoint *endpoint); +/* Applies the completed generation. Returns the number of newly added + * vertices, zero when no topology change was made, or -1 on failure. */ +int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, + const RefinementConfig *config); +int frame_lens_mesh_refine(FrameLensMesh *mesh, + const SpacetimeSource *spacetime, + const ObserverState *observer, + const GeodesicTraceConfig *trace, + const RefinementConfig *config); /* Each locally invertible escaped triangle contributes one image per contained * star. */ diff --git a/src/main.c b/src/main.c index 926f30e..a42d348 100644 --- a/src/main.c +++ b/src/main.c @@ -40,6 +40,7 @@ typedef struct { double slab_duration; double minkowski_proper_acceleration; int catalog_load_workers; + RefinementConfig refinement; } Settings; static int parse_int(const char *text, int *value) { @@ -52,6 +53,16 @@ static int parse_int(const char *text, int *value) { return 0; } +static int parse_nonnegative_int(const char *text, unsigned int *value) { + char *end; + errno = 0; + unsigned long parsed = strtoul(text, &end, 10); + if (errno || *end || parsed > 1024) + return -1; + *value = (unsigned int)parsed; + return 0; +} + static int parse_double(const char *text, double *value) { char *end; errno = 0; @@ -139,7 +150,12 @@ static int parse_args(int argc, char **argv, Settings *s, .movie_fps = 30.0, .slab_duration = 64.0, .minkowski_proper_acceleration = 1.52, - .catalog_load_workers = 4}; + .catalog_load_workers = 4, + .refinement = {.angle_absolute_rad = + 1e-3 * 3.14159265358979323846 / 180.0, + .angle_relative = 0.1, + .min_edge_pixels = 0.5, + .min_area_pixels2 = 0.25}}; *write_path = NULL; for (int i = 1; i < argc; ++i) { if (!strcmp(argv[i], "--catalog") && i + 1 < argc) @@ -158,6 +174,17 @@ static int parse_args(int argc, char **argv, Settings *s, !parse_int(argv[++i], &s->height)) { } else if (!strcmp(argv[i], "--coarse-cell-pixels") && i + 1 < argc && !parse_int(argv[++i], &s->coarse_cell_pixels)) { + } else if (!strcmp(argv[i], "--refine-max-level") && i + 1 < argc && + !parse_nonnegative_int(argv[++i], &s->refinement.max_level)) { + } else if (!strcmp(argv[i], "--refine-angle-abs-deg") && i + 1 < argc && + !parse_positive(argv[++i], &s->refinement.angle_absolute_rad)) { + s->refinement.angle_absolute_rad *= 3.14159265358979323846 / 180.0; + } else if (!strcmp(argv[i], "--refine-angle-rel") && i + 1 < argc && + !parse_positive(argv[++i], &s->refinement.angle_relative)) { + } else if (!strcmp(argv[i], "--refine-min-edge-pixels") && i + 1 < argc && + !parse_positive(argv[++i], &s->refinement.min_edge_pixels)) { + } else if (!strcmp(argv[i], "--refine-min-area-pixels2") && i + 1 < argc && + !parse_positive(argv[++i], &s->refinement.min_area_pixels2)) { } else if (!strcmp(argv[i], "--draw-mesh")) { s->draw_mesh = 1; } else if (!strcmp(argv[i], "--fov-deg") && i + 1 < argc && @@ -257,19 +284,6 @@ static void report_splat_progress(void *context, FrameSplatProgressStage stage, } } -static void ray_pool_status_counts(const RayPool *rays, size_t *pending, - size_t *active, size_t *terminated, - size_t *failed) { - *pending = *active = *terminated = *failed = 0; - for (size_t i = 0; i < rays->count; ++i) - switch (rays->status[i]) { - case RAY_POOL_PENDING: ++*pending; break; - case RAY_POOL_ACTIVE: ++*active; break; - case RAY_POOL_TERMINATED: ++*terminated; break; - case RAY_POOL_FAILED: ++*failed; break; - } -} - static GeodesicTraceConfig trace_config(void) { #ifdef SPACETIME_SCHWARZSCHILD return (GeodesicTraceConfig){.coordinate_time_step = 0.1, @@ -303,7 +317,9 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog, frame_lens_mesh_build_coarse(&mesh, s->width, s->height, s->coarse_cell_pixels, s->horizontal_fov_deg) || - frame_lens_mesh_trace(&mesh, spacetime, observer, &trace)) { + frame_lens_mesh_trace(&mesh, spacetime, observer, &trace) || + frame_lens_mesh_refine(&mesh, spacetime, observer, &trace, + &s->refinement)) { frame_lens_mesh_destroy(&mesh); free(hdr); return -1; @@ -370,11 +386,65 @@ static int frame_output_path(char path[PATH_MAX], const Settings *s, return written < 0 || written >= PATH_MAX ? -1 : 0; } +static int trace_movie_generation(Movie *movie, const Settings *s, + const SpacetimeSource *spacetime, + const GeodesicTraceConfig *trace, + size_t generation) { + RayPool rays = {0}; + size_t ray_count = 0; + for (size_t f = 0; f < movie->frame_count; ++f) { + const int prepared = frame_lens_mesh_prepare_generation( + &movie->frames[f].mesh, &s->refinement); + if (prepared < 0) return -1; + ray_count += (size_t)prepared; + } + if (ray_count == 0) return 0; + if (ray_pool_init(&rays, ray_count)) return -1; + for (size_t f = 0; f < movie->frame_count; ++f) { + size_t count = 0; + const FrameSample *samples = frame_lens_mesh_samples(&movie->frames[f].mesh, &count); + for (size_t sample = 0; sample < count; ++sample) + if (ray_pool_append(&rays, &movie->frames[f].observer, + samples[sample].vertex.camera_direction, f, sample)) { + ray_pool_destroy(&rays); + return -1; + } + } + double slab_hi = movie->frames[movie->frame_count - 1].coordinate_time; + size_t slab_id = 0; + while (ray_pool_has_live(&rays)) { + const double slab_lo = slab_hi - s->slab_duration; + MetricSlab *slab = NULL; + if (spacetime_load_slab(spacetime, slab_hi, slab_lo, &slab)) { + ray_pool_destroy(&rays); + return -1; + } + ray_pool_activate_in_time_range(&rays, slab); + ray_pool_advance_active(&rays, slab, trace); + spacetime_free_slab(slab); + if (s->verbose) + fprintf(stderr, "Refinement generation %zu, slab %zu completed.\n", + generation, ++slab_id); + slab_hi = slab_lo; + } + for (size_t i = 0; i < rays.count; ++i) + if (frame_lens_mesh_install_sample(&movie->frames[rays.frame_id[i]].mesh, + rays.vertex_id[i], &rays.endpoint[i])) { + ray_pool_destroy(&rays); + return -1; + } + ray_pool_destroy(&rays); + for (size_t f = 0; f < movie->frame_count; ++f) + if (frame_lens_mesh_finish_generation(&movie->frames[f].mesh, + &s->refinement) < 0) + return -1; + return 1; +} + static int render_movie(const Settings *s, StarCatalog *catalog, const SpacetimeSource *spacetime) { ObserverTrack track = {0}; Movie movie = {0}; - RayPool rays = {0}; const GeodesicTraceConfig trace = trace_config(); int result = -1; if (s->observer_track_path == NULL || @@ -384,55 +454,13 @@ static int render_movie(const Settings *s, StarCatalog *catalog, movie_build_coarse_meshes(&movie, s->width, s->height, s->coarse_cell_pixels, s->horizontal_fov_deg)) goto done; - size_t ray_count = 0; - for (size_t i = 0; i < movie.frame_count; ++i) - ray_count += movie.frames[i].mesh.vertex_count; - if (ray_count == 0 || ray_pool_init(&rays, ray_count)) - goto done; - for (size_t f = 0; f < movie.frame_count; ++f) - for (size_t v = 0; v < movie.frames[f].mesh.vertex_count; ++v) - if (ray_pool_append(&rays, &movie.frames[f].observer, - movie.frames[f].mesh.vertices[v].camera_direction, - f, v)) - goto done; - double slab_hi = movie.frames[movie.frame_count - 1].coordinate_time; - size_t slab_id = 0; - while (ray_pool_has_live(&rays)) { - const double slab_lo = slab_hi - s->slab_duration; - MetricSlab *slab = NULL; - size_t pending_before, active_before, terminated_before, failed_before; - ray_pool_status_counts(&rays, &pending_before, &active_before, - &terminated_before, &failed_before); - if (s->verbose) - fprintf(stderr, "Time slab %zu: loading [%.6g, %.6g] with %zu pending and %zu active rays.\n", - slab_id + 1, slab_hi, slab_lo, pending_before, active_before); - if (spacetime_load_slab(spacetime, slab_hi, slab_lo, &slab)) + for (size_t generation = 0;; ++generation) { + const int traced = trace_movie_generation(&movie, s, spacetime, &trace, + generation); + if (traced < 0) goto done; - ray_pool_activate_in_time_range(&rays, slab); - size_t pending_active, active_active, terminated_active, failed_active; - ray_pool_status_counts(&rays, &pending_active, &active_active, - &terminated_active, &failed_active); - ray_pool_advance_active(&rays, slab, &trace); - spacetime_free_slab(slab); - size_t pending_after, active_after, terminated_after, failed_after; - ray_pool_status_counts(&rays, &pending_after, &active_after, - &terminated_after, &failed_after); - fprintf(stderr, - "Time slab %zu [%.6g, %.6g]: activated %zu; live %zu -> %zu, " - "terminated %zu, failed %zu\n", - ++slab_id, slab_hi, slab_lo, active_active - active_before, - pending_before + active_before, pending_after + active_after, - terminated_after, failed_after); - slab_hi = slab_lo; - } - for (size_t i = 0; i < rays.count; ++i) { - LensVertex *vertex = &movie.frames[rays.frame_id[i]].mesh.vertices[rays.vertex_id[i]]; - vertex->status = rays.endpoint[i].status; - if (vertex->status == RAY_ENDPOINT_ESCAPED) { - for (int axis = 0; axis < 3; ++axis) - vertex->n_infinity[axis] = rays.endpoint[i].n_infinity[axis]; - vertex->log_frequency_ratio = log(rays.endpoint[i].frequency_ratio); - } + if (traced == 0) + break; } for (size_t i = 0; i < movie.frame_count; ++i) { char output_path[PATH_MAX]; @@ -467,7 +495,6 @@ static int render_movie(const Settings *s, StarCatalog *catalog, } result = 0; done: - ray_pool_destroy(&rays); movie_destroy(&movie); observer_track_destroy(&track); return result; @@ -499,7 +526,10 @@ int main(int argc, char **argv) { #ifdef ENABLE_HDR_DEBUG "[--hdr-output PATH] " #endif - "[--coarse-cell-pixels N] [--draw-mesh] [--write-catalog PATH] " + "[--coarse-cell-pixels N] [--refine-max-level N " + "--refine-angle-abs-deg D --refine-angle-rel R " + "--refine-min-edge-pixels P --refine-min-area-pixels2 A] " + "[--draw-mesh] [--write-catalog PATH] " "[--catalog-load-workers N] " "[--observer-track PATH --frames-dir DIR --frames-prefix NAME " "--start-time T --duration T --fps N] " diff --git a/tests/test_frame.c b/tests/test_frame.c index e63c4e9..be8d520 100644 --- a/tests/test_frame.c +++ b/tests/test_frame.c @@ -160,6 +160,38 @@ int main(void) { goto done; } frame_lens_mesh_destroy(&fine_mesh); + /* Refinement probes are temporary until their generation is complete. A + * shared diagonal probe must produce one stable midpoint and conforming + * children only after its endpoint has been installed. */ + FrameLensMesh adaptive_mesh = {0}; + const RefinementConfig refine = {.max_level = 1, + .angle_absolute_rad = 1e-4, + .angle_relative = 1e-4, + .min_edge_pixels = 1.0, + .min_area_pixels2 = 1.0}; + if (frame_lens_mesh_build_coarse(&adaptive_mesh, width, height, 100, 30.0)) + goto done; + for (size_t i = 0; i < adaptive_mesh.vertex_count; ++i) { + adaptive_mesh.vertices[i].traced = 1; + adaptive_mesh.vertices[i].status = RAY_ENDPOINT_ESCAPED; + adaptive_mesh.vertices[i].n_infinity[0] = 1.0; + } + if (frame_lens_mesh_prepare_generation(&adaptive_mesh, &refine) != 1) { + fputs("adaptive shared-edge probe setup regression failed\n", stderr); + frame_lens_mesh_destroy(&adaptive_mesh); + goto done; + } + const RayEndpoint bent_probe = {.n_infinity = {0.0, 1.0, 0.0}, + .frequency_ratio = 1.0, + .status = RAY_ENDPOINT_ESCAPED}; + if (frame_lens_mesh_install_sample(&adaptive_mesh, 0, &bent_probe) || + frame_lens_mesh_finish_generation(&adaptive_mesh, &refine) != 1 || + adaptive_mesh.vertex_count != 5 || adaptive_mesh.triangle_count != 4) { + fputs("adaptive shared-edge split regression failed\n", stderr); + frame_lens_mesh_destroy(&adaptive_mesh); + goto done; + } + frame_lens_mesh_destroy(&adaptive_mesh); result = 0; done: frame_lens_mesh_destroy(&mesh); diff --git a/tests/test_schwarzschild.c b/tests/test_schwarzschild.c index 7afcef3..784db15 100644 --- a/tests/test_schwarzschild.c +++ b/tests/test_schwarzschild.c @@ -1,4 +1,5 @@ #include "geodesic.h" +#include "frame.h" #include #include @@ -56,6 +57,25 @@ int main(void) { central.status, inside_shadow.status, outside_shadow.status); goto done; } + /* A coarse field covering the shadow must genuinely refine: its initial + * capture/escape-discontinuous triangles are a separate trigger from the + * smooth direction-error criterion. */ + FrameLensMesh mesh = {0}; + const RefinementConfig refinement = {.max_level = 1, + .angle_absolute_rad = 1e-5, + .angle_relative = 1e-5, + .min_edge_pixels = 1.0, + .min_area_pixels2 = 1.0}; + if (frame_lens_mesh_build_coarse(&mesh, 48, 48, 24, 40.0) || + frame_lens_mesh_trace(&mesh, &spacetime, &observer, &trace) || + frame_lens_mesh_refine(&mesh, &spacetime, &observer, &trace, + &refinement) || + mesh.vertex_count <= 9 || mesh.triangle_count <= 8) { + fputs("Schwarzschild adaptive-refinement regression failed\n", stderr); + frame_lens_mesh_destroy(&mesh); + goto done; + } + frame_lens_mesh_destroy(&mesh); result = 0; done: spacetime_destroy(&spacetime); diff --git a/video_rendering_plan.md b/video_rendering_plan.md index 534fe01..8bb17ce 100644 --- a/video_rendering_plan.md +++ b/video_rendering_plan.md @@ -45,7 +45,7 @@ SampleRequest -> RayPool (SoA, inactive / active / terminated) | A | 基本完成 | canonical 21 列 observer-track CSV、插值、movie frame schedule、编号 PNG/PPM 序列和相关 CLI 已实现。 | 输出仍由 `optics` 写入，尚未拆出计划中的 `output.h/.c`；CSV 非有限数、非单调时间和 PNG 序号冲突的独立回归仍需补齐。 | | B | 部分完成 | Minkowski 恒 proper-acceleration 轨迹生成器、`0..2`/30 fps/61 frame 基准和前向解析 Doppler regression 已实现。 | 尚无逐 frame 诊断 CSV；a=0 全序列、后向红移、tetrad 正交性/有限差分速度和图像连续性回归仍需补齐。 | | C | 核心路径完成，验收未完成 | `MetricSlab` 接口、SoA `RayPool`、slab 到达时的 lazy ray initialization、倒序 slab sweep、OpenMP 批量推进和 `(frame_id, vertex_id)` endpoint 回填已实现。 | 仍未 compact terminated rays 或支持 RayPool 增长；未记录/测试每 generation 每 slab 一次 load/free；scheduler 与 wrapper 的完整 endpoint/image 对比尚未成为当前测试。 | -| D | 未开始 | 无。 | multi-pass adaptive refinement、请求队列、nmesh slab I/O/ghost-slice/memory contract。 | +| D | 部分完成 | 多代 request/probe mesh、最长边二分、共享边 conforming 闭包、相对+绝对天空方向误差、capture/escape 边界细分，以及单帧/movie 共用 generation 调度已实现。 | nmesh slab I/O/ghost-slice/memory contract、slab instrumentation 和强透镜收敛基准仍未完成。 | | E | 未开始 | 仅有静态 Schwarzschild 单帧和局部 inward boost。 | `30M -> 20M -> 30M` observer track、movie regression 和连续性检查。 | | F | 未开始 | PNG sequence 是默认 movie 输出。 | FFmpeg 编码命令/目标、解码核验和可选 direct-video backend。 | @@ -88,7 +88,7 @@ SampleRequest -> RayPool (SoA, inactive / active / terminated) ## Phase D：补全多 pass adaptive mesh 与 out-of-core 准备（未开始） -1. 给 `FrameLensMesh` 增加稳定 vertex ID、请求队列、triangle level/flags 和“缺失 endpoint”状态。每 pass 只为需要的 edge midpoint/center 产生 `SampleRequest`；安装结果后再判定 mapping interpolation error、orientation consistency 与局部 Jacobian 奇异性。 +1. `FrameLensMesh` 采用 append-only stable vertex ID、generation-local request/probe 队列、triangle level 与“缺失 endpoint”状态。每 pass 的 probe 只放在最长 image-plane 边中点；完整 sweep 后才安装 endpoint 并改变拓扑。细分要求同时超过绝对与相对 `n_infinity` 方向误差阈值，并受最大层数、最小长边和最小面积约束；capture/escape 不一致强制细分。首版不计算 orientation 或 Jacobian。 2. 以迭代 queue（可按 frame/root tile 并行、线程本地 request buffer 后 sort/deduplicate）替代递归 task。每一个 pass 完整执行 Phase C sweep；只有所有 frame 都无新请求才进行 catalog splat。 3. 在 movie 生命周期中及时释放已完成的 RayPool、临时 request 和单帧 HDR；保留最终 mesh/endpoints，或在渲染 PNG 后按明确策略释放，避免视频时无界增长。 4. 为 nmesh 预留并实现 source-side slab overlap / temporal ghost-slice 契约、可配置 memory budget、slab coverage 日志和线程本地 `MetricWorkspace` ownership。此 phase 不重采样为 Cartesian grid，也不假定相邻时间 slice 的 AMR tree 相同。