Frame: localize shadow boundary refinement

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 unsigned int first_side = longest_side(mesh, triangle);
-    const size_t a = triangle->vertex[side];
+    const unsigned int side_count = terminal_mismatch(
-    const size_t b = triangle->vertex[(side + 1) % 3];
+                                        &mesh->vertices[triangle->vertex[0]],
-    if (find_probe(mesh, a, b) != SIZE_MAX)
+                                        &mesh->vertices[triangle->vertex[1]],
-      continue;
+                                        &mesh->vertices[triangle->vertex[2]])
-    FrameSample probe = {.kind = FRAME_SAMPLE_PROBE, .edge_vertex = {a, b}};
+                                    ? 3
-    const LensVertex *left = &mesh->vertices[a];
+                                    : 1;
-    const LensVertex *right = &mesh->vertices[b];
+    for (unsigned int offset = 0; offset < side_count; ++offset) {
-    probe.vertex.image_x = 0.5 * (left->image_x + right->image_x);
+      const unsigned int side = (first_side + offset) % 3;
-    probe.vertex.image_y = 0.5 * (left->image_y + right->image_y);
+      const size_t a = triangle->vertex[side];
-    for (int axis = 0; axis < 3; ++axis)
+      const size_t b = triangle->vertex[(side + 1) % 3];
-      probe.vertex.camera_direction[axis] =
+      if (find_probe(mesh, a, b) != SIZE_MAX)
-          left->camera_direction[axis] + right->camera_direction[axis];
+        continue;
-    if (normalize(probe.vertex.camera_direction) == 0.0)
+      FrameSample probe = {.kind = FRAME_SAMPLE_PROBE, .edge_vertex = {a, b}};
-      return -1;
+      const LensVertex *left = &mesh->vertices[a];
-    if (add_sample(mesh, &probe))
+      const LensVertex *right = &mesh->vertices[b];
-      return -1;
+      probe.vertex.image_x = 0.5 * (left->image_x + right->image_x);
-    index_probe(mesh, mesh->sample_count - 1);
+      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);
+      if (terminal_mismatch(&mesh->vertices[triangle->vertex[0]],
-      requested[3 * i + side] = probe_requires_split(mesh, triangle, side, config);
+                            &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,
-        append_triangle(children, &child_count, old_count * 4, a, b, c, triangle->level, 1);
+                        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);

tests/test_frame.c

@@ -7,6 +7,56 @@
 #include <stdlib.h>
 #include <string.h>
 
+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};