diff --git a/src/frame.c b/src/frame.c index 7c8e368..0034a1f 100644 --- a/src/frame.c +++ b/src/frame.c @@ -190,6 +190,17 @@ static int all_vertices_traced(const FrameLensMesh *mesh) { 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; @@ -271,24 +282,33 @@ int frame_lens_mesh_prepare_generation(FrameLensMesh *mesh, 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); + 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; @@ -319,17 +339,6 @@ int frame_lens_mesh_install_sample(FrameLensMesh *mesh, size_t sample_id, 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); @@ -430,6 +439,30 @@ static LensVertex midpoint_vertex(const LensVertex *a, const LensVertex *b) { 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) { @@ -439,6 +472,31 @@ static int append_triangle(LensTriangle *triangles, size_t *count, 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) @@ -483,8 +541,19 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, edges[3 * i + side] = (MeshEdge){a, b, i, side}; } if (allowed[i]) { - const unsigned int side = longest_side(mesh, triangle); - requested[3 * i + side] = probe_requires_split(mesh, triangle, side, config); + 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]); } } @@ -533,41 +602,6 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, 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; @@ -646,11 +680,39 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, 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. */ + } 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) { - /* 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); + 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); diff --git a/tests/test_frame.c b/tests/test_frame.c index a108914..51e95df 100644 --- a/tests/test_frame.c +++ b/tests/test_frame.c @@ -7,6 +7,56 @@ #include #include +static int mesh_has_hanging_vertex(const FrameLensMesh *mesh) { + for (size_t triangle = 0; triangle < mesh->triangle_count; ++triangle) + for (size_t side = 0; side < 3; ++side) { + const LensVertex *a = &mesh->vertices[mesh->triangles[triangle].vertex[side]]; + const LensVertex *b = + &mesh->vertices[mesh->triangles[triangle].vertex[(side + 1) % 3]]; + const double dx = b->image_x - a->image_x; + const double dy = b->image_y - a->image_y; + const double length_squared = dx * dx + dy * dy; + for (size_t vertex = 0; vertex < mesh->vertex_count; ++vertex) { + if (vertex == mesh->triangles[triangle].vertex[side] || + vertex == mesh->triangles[triangle].vertex[(side + 1) % 3]) + continue; + const LensVertex *p = &mesh->vertices[vertex]; + const double px = p->image_x - a->image_x; + const double py = p->image_y - a->image_y; + const double cross = px * dy - py * dx; + const double position = (px * dx + py * dy) / length_squared; + if (fabs(cross) <= 1e-12 * length_squared && position > 1e-12 && + position < 1.0 - 1e-12) + return 1; + } + } + return 0; +} + +static int mesh_has_same_winding_shared_edge(const FrameLensMesh *mesh) { + for (size_t left_triangle = 0; left_triangle < mesh->triangle_count; + ++left_triangle) + for (size_t left_side = 0; left_side < 3; ++left_side) { + const size_t from = mesh->triangles[left_triangle].vertex[left_side]; + const size_t to = + mesh->triangles[left_triangle].vertex[(left_side + 1) % 3]; + for (size_t right_triangle = left_triangle + 1; + right_triangle < mesh->triangle_count; ++right_triangle) + for (size_t right_side = 0; right_side < 3; ++right_side) { + const size_t other_from = + mesh->triangles[right_triangle].vertex[right_side]; + const size_t other_to = + mesh->triangles[right_triangle].vertex[(right_side + 1) % 3]; + if ((from == other_from && to == other_to) || + (from == other_to && to == other_from)) { + if (from == other_from && to == other_to) + return 1; + } + } + } + return 0; +} + int main(void) { const int width = 100, height = 100; const double test_exposure = 1e-3; @@ -213,6 +263,72 @@ int main(void) { goto done; } frame_lens_mesh_destroy(&adaptive_mesh); + /* A capture/escape discontinuity is a shadow boundary, not a smooth map + * error: request all three midpoint rays and red-refine in one generation. */ + if (frame_lens_mesh_build_coarse(&adaptive_mesh, width, height, 100, 30.0)) + goto done; + adaptive_mesh.triangle_count = 1; + 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; + } + adaptive_mesh.vertices[0].status = RAY_ENDPOINT_CAPTURED; + refine.max_level = 1; + refine.angle_absolute_rad = 3.14159265358979323846; + refine.angle_relative = 1e6; + refine.jacobian_minimum = 1e-12; + if (frame_lens_mesh_prepare_generation(&adaptive_mesh, &refine) != 3) { + fputs("shadow-boundary red-probe setup regression failed\n", stderr); + frame_lens_mesh_destroy(&adaptive_mesh); + goto done; + } + for (size_t i = 0; i < adaptive_mesh.sample_count; ++i) + if (frame_lens_mesh_install_sample(&adaptive_mesh, i, &bent_probe)) { + fputs("shadow-boundary red-probe installation regression failed\n", stderr); + frame_lens_mesh_destroy(&adaptive_mesh); + goto done; + } + if (frame_lens_mesh_finish_generation(&adaptive_mesh, &refine) != 3 || + adaptive_mesh.vertex_count != 7 || adaptive_mesh.triangle_count != 4) { + fputs("shadow-boundary red-refinement regression failed\n", stderr); + frame_lens_mesh_destroy(&adaptive_mesh); + goto done; + } + frame_lens_mesh_destroy(&adaptive_mesh); + /* Two shadow leaves can force two edges of an escaped neighbour. That + * neighbour must use a local three-child blue split, not create a third + * requested edge that spreads red refinement farther outward. */ + if (frame_lens_mesh_build_coarse(&adaptive_mesh, 200, 100, 100, 30.0)) + goto done; + adaptive_mesh.triangle_count = 3; + 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; + } + adaptive_mesh.vertices[3].status = RAY_ENDPOINT_CAPTURED; + adaptive_mesh.vertices[5].status = RAY_ENDPOINT_CAPTURED; + if (frame_lens_mesh_prepare_generation(&adaptive_mesh, &refine) != 6) { + fputs("shadow-boundary blue-neighbour probe setup regression failed\n", stderr); + frame_lens_mesh_destroy(&adaptive_mesh); + goto done; + } + for (size_t i = 0; i < adaptive_mesh.sample_count; ++i) + if (frame_lens_mesh_install_sample(&adaptive_mesh, i, &bent_probe)) { + fputs("shadow-boundary blue-neighbour probe installation regression failed\n", stderr); + frame_lens_mesh_destroy(&adaptive_mesh); + goto done; + } + if (frame_lens_mesh_finish_generation(&adaptive_mesh, &refine) != 6 || + adaptive_mesh.vertex_count != 12 || adaptive_mesh.triangle_count != 11 || + mesh_has_hanging_vertex(&adaptive_mesh) || + mesh_has_same_winding_shared_edge(&adaptive_mesh)) { + fputs("shadow-boundary blue-neighbour refinement regression failed\n", stderr); + frame_lens_mesh_destroy(&adaptive_mesh); + goto done; + } + frame_lens_mesh_destroy(&adaptive_mesh); const RayEndpoint flat_probe = {.n_infinity = {1.0, 0.0, 0.0}, .frequency_ratio = 1.0, .status = RAY_ENDPOINT_ESCAPED};