Frame: guard sliver-producing refinement splits

README.md

@@ -148,8 +148,8 @@ 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, minimum area `0.25` pixel-squared, and a provisional
-minimum discrete-Jacobian magnitude of `1e-3`.  Each value can be overridden
-independently:
+minimum discrete-Jacobian magnitude of `1e-3`, and minimum triangle quality
+`0.2`.  Each value can be overridden independently:
 
 ```text
 --refine-max-level N
@@ -158,6 +158,7 @@ independently:
 --refine-jacobian-min J
 --refine-min-edge-pixels P
 --refine-min-area-pixels2 A
+--refine-min-triangle-quality Q
 ```
 
 `N` caps the triangle refinement level.  Let `e` be the angle between the
@@ -170,7 +171,12 @@ 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.
+long-edge and area scales.  `Q` is the normalized image-plane triangle quality
+`4 sqrt(3) A / (l0^2 + l1^2 + l2^2)`: it is one for an equilateral triangle
+and approaches zero for a sliver.  Its default is `0.2`; `0` disables this
+quality promotion.  If a requested
+single-edge split would create a child below `Q`, that leaf is promoted to red
+refinement before the conforming shared-edge closure is applied.
 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.

src/frame.c

@@ -139,6 +139,20 @@ static double image_triangle_area(const LensVertex *a, const LensVertex *b,
                     (b->image_y - a->image_y) * (c->image_x - a->image_x));
 }
 
+/* Four times sqrt(3) times the area divided by the squared-edge sum.  This
+ * is scale-independent, equals one for an equilateral triangle, and tends to
+ * zero for a sliver. */
+static double image_triangle_quality(const LensVertex *a, const LensVertex *b,
+                                     const LensVertex *c) {
+  const double ab = image_edge_length(a, b);
+  const double bc = image_edge_length(b, c);
+  const double ca = image_edge_length(c, a);
+  const double squared_edges = ab * ab + bc * bc + ca * ca;
+  if (squared_edges == 0.0)
+    return 0.0;
+  return 4.0 * sqrt(3.0) * image_triangle_area(a, b, c) / squared_edges;
+}
+
 static unsigned int longest_side(const FrameLensMesh *mesh,
                                  const LensTriangle *triangle) {
   unsigned int best = 0;
@@ -430,6 +444,24 @@ static LensVertex midpoint_vertex(const LensVertex *a, const LensVertex *b) {
   return result;
 }
 
+/* A conformity request can require splitting a neighbor along a non-longest
+ * edge.  Unlike longest-edge bisection, that can create a sliver.  Promote
+ * that leaf to red refinement before allocating vertices when either of the
+ * two prospective children violates the requested image-plane quality. */
+static int single_edge_split_requires_red(const FrameLensMesh *mesh,
+                                          const LensTriangle *triangle,
+                                          unsigned int side,
+                                          const RefinementConfig *config) {
+  if (config->min_triangle_quality <= 0.0)
+    return 0;
+  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]];
+  const LensVertex midpoint = midpoint_vertex(a, b);
+  return image_triangle_quality(a, &midpoint, c) < config->min_triangle_quality ||
+         image_triangle_quality(&midpoint, b, c) < config->min_triangle_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) {
@@ -526,20 +558,34 @@ 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. */
+  /* Two requested sides require red refinement.  A single requested edge
+   * whose bisection would form a sliver is also promoted to red refinement.
+   * Close every new shared-edge request before allocating any vertices. */
   for (;;) {
     int changed = 0;
     for (size_t t = 0; t < mesh->triangle_count; ++t) {
       unsigned int count = 0;
+      unsigned int requested_side = 0;
       for (unsigned int side = 0; side < 3; ++side)
-        count += requested[3 * t + side] != 0;
-      if (count >= 2 && allowed[t])
+        if (requested[3 * t + side]) {
+          ++count;
+          requested_side = side;
+        }
+      if (count == 1 && allowed[t] &&
+          single_edge_split_requires_red(mesh, &mesh->triangles[t],
+                                         requested_side, config)) {
         for (unsigned int side = 0; side < 3; ++side)
           if (!requested[3 * t + side]) {
             requested[3 * t + side] = 1;
             changed = 1;
           }
+      } else 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;

src/frame.h

@@ -31,6 +31,9 @@ typedef struct {
   double jacobian_minimum;
   double min_edge_pixels;
   double min_area_pixels2;
+  /* Normalized image-plane quality, 1 for equilateral and 0 for degenerate.
+   * A non-positive value disables quality-driven red refinement. */
+  double min_triangle_quality;
 } RefinementConfig;
 
 typedef enum {

src/main.c

@@ -91,6 +91,13 @@ static int parse_positive(const char *text, double *value) {
   return errno || *end || *value <= 0.0 ? -1 : 0;
 }
 
+static int parse_unit_interval(const char *text, double *value) {
+  char *end;
+  errno = 0;
+  *value = strtod(text, &end);
+  return errno || *end || *value < 0.0 || *value > 1.0 ? -1 : 0;
+}
+
 static int parse_nonnegative(const char *text, double *value) {
   char *end;
   errno = 0;
@@ -156,7 +163,8 @@ static int parse_args(int argc, char **argv, Settings *s,
                                  .angle_relative = 0.1,
                                  .jacobian_minimum = 1e-3,
                                  .min_edge_pixels = 0.5,
-                                 .min_area_pixels2 = 0.25}};
+                                 .min_area_pixels2 = 0.25,
+                                 .min_triangle_quality = 0.2}};
   *write_path = NULL;
   for (int i = 1; i < argc; ++i) {
     if (!strcmp(argv[i], "--catalog") && i + 1 < argc)
@@ -188,6 +196,8 @@ static int parse_args(int argc, char **argv, Settings *s,
                !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], "--refine-min-triangle-quality") && i + 1 < argc &&
+               !parse_unit_interval(argv[++i], &s->refinement.min_triangle_quality)) {
     } else if (!strcmp(argv[i], "--draw-mesh")) {
       s->draw_mesh = 1;
     } else if (!strcmp(argv[i], "--fov-deg") && i + 1 < argc &&
@@ -532,7 +542,8 @@ int main(int argc, char **argv) {
             "[--coarse-cell-pixels N] [--refine-max-level N "
             "--refine-angle-abs-deg D --refine-angle-rel R "
             "--refine-jacobian-min J "
-            "--refine-min-edge-pixels P --refine-min-area-pixels2 A] "
+            "--refine-min-edge-pixels P --refine-min-area-pixels2 A "
+            "--refine-min-triangle-quality Q] "
             "[--draw-mesh] [--write-catalog PATH] "
             "[--catalog-load-workers N] "
             "[--observer-track PATH --frames-dir DIR --frames-prefix NAME "

tests/test_frame.c

@@ -193,6 +193,92 @@ int main(void) {
     goto done;
   }
   frame_lens_mesh_destroy(&adaptive_mesh);
+  /* A shared fold edge can be shorter than either triangle's longest edge.
+   * Splitting it alone would produce low-quality children, so the quality
+   * guard must promote both leaves to conforming red refinement. */
+  LensVertex *quality_vertices = calloc(4, sizeof *quality_vertices);
+  LensTriangle *quality_triangles = calloc(2, sizeof *quality_triangles);
+  if (quality_vertices == NULL || quality_triangles == NULL) {
+    free(quality_vertices);
+    free(quality_triangles);
+    goto done;
+  }
+  quality_vertices[0].image_x = 0.0;
+  quality_vertices[0].image_y = 0.0;
+  quality_vertices[1].image_x = 1.0;
+  quality_vertices[1].image_y = 0.0;
+  quality_vertices[2].image_x = 0.5;
+  quality_vertices[2].image_y = 1.0;
+  quality_vertices[3].image_x = 0.5;
+  quality_vertices[3].image_y = -1.0;
+  const double quality_source_directions[4][3] = {
+      {1.0, 0.0, 0.0}, {1.0, 0.1, 0.0},
+      {1.0, 0.0, 0.1}, {1.0, 0.0, 0.1}};
+  for (size_t i = 0; i < 4; ++i) {
+    quality_vertices[i].traced = 1;
+    quality_vertices[i].status = RAY_ENDPOINT_ESCAPED;
+    quality_vertices[i].camera_direction[0] = 1.0;
+    quality_vertices[i].camera_direction[1] = 0.01 * quality_vertices[i].image_x;
+    quality_vertices[i].camera_direction[2] = 0.01 * quality_vertices[i].image_y;
+    double length = sqrt(quality_vertices[i].camera_direction[0] *
+                             quality_vertices[i].camera_direction[0] +
+                         quality_vertices[i].camera_direction[1] *
+                             quality_vertices[i].camera_direction[1] +
+                         quality_vertices[i].camera_direction[2] *
+                             quality_vertices[i].camera_direction[2]);
+    for (size_t axis = 0; axis < 3; ++axis)
+      quality_vertices[i].camera_direction[axis] /= length;
+    memcpy(quality_vertices[i].n_infinity, quality_source_directions[i],
+           sizeof quality_source_directions[i]);
+    length = sqrt(quality_vertices[i].n_infinity[0] *
+                      quality_vertices[i].n_infinity[0] +
+                  quality_vertices[i].n_infinity[1] *
+                      quality_vertices[i].n_infinity[1] +
+                  quality_vertices[i].n_infinity[2] *
+                      quality_vertices[i].n_infinity[2]);
+    for (size_t axis = 0; axis < 3; ++axis)
+      quality_vertices[i].n_infinity[axis] /= length;
+  }
+  quality_triangles[0] = (LensTriangle){{0, 1, 2}, 0, 0};
+  quality_triangles[1] = (LensTriangle){{1, 0, 3}, 0, 0};
+  adaptive_mesh = (FrameLensMesh){.vertices = quality_vertices,
+                                  .triangles = quality_triangles,
+                                  .vertex_count = 4,
+                                  .vertex_capacity = 4,
+                                  .triangle_count = 2,
+                                  .triangle_capacity = 2};
+  refine.angle_absolute_rad = 3.14159265358979323846;
+  refine.angle_relative = 1e6;
+  refine.jacobian_minimum = 1000.0;
+  refine.min_edge_pixels = 0.1;
+  refine.min_area_pixels2 = 0.01;
+  refine.min_triangle_quality = 0.8;
+  const RayEndpoint quality_probe = {.n_infinity = {1.0, 0.0, 0.0},
+                                     .frequency_ratio = 1.0,
+                                     .status = RAY_ENDPOINT_ESCAPED};
+  const int quality_prepared =
+      frame_lens_mesh_prepare_generation(&adaptive_mesh, &refine);
+  if (quality_prepared != 2) {
+    fprintf(stderr, "adaptive quality probe setup regression failed (%d)\n",
+            quality_prepared);
+    goto done;
+  }
+  for (size_t i = 0; i < adaptive_mesh.sample_count; ++i)
+    if (frame_lens_mesh_install_sample(&adaptive_mesh, i, &quality_probe))
+      goto done;
+  const int quality_midpoints =
+      frame_lens_mesh_finish_generation(&adaptive_mesh, &refine);
+  if (quality_midpoints != 5 ||
+      adaptive_mesh.vertex_count != 9 || adaptive_mesh.triangle_count != 8) {
+    fprintf(stderr,
+            "adaptive quality red-refinement regression failed (%d midpoints, "
+            "%zu vertices, %zu triangles)\n",
+            quality_midpoints, adaptive_mesh.vertex_count,
+            adaptive_mesh.triangle_count);
+    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};