Frame: refine critical lens regions

src/frame.c

@@ -335,6 +335,37 @@ static double direction_angle(const double a[3], const double b[3]) {
   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.
+ * Captured and degenerate image triangles have no reliable value. */
+static int discrete_jacobian(const FrameLensMesh *mesh,
+                             const LensTriangle *triangle, double *value) {
+  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;
+  return 1;
+}
+
 static int probe_requires_split(const FrameLensMesh *mesh,
                                 const LensTriangle *triangle,
                                 unsigned int side,
@@ -441,6 +472,10 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh,
     if (allowed[i]) {
       const unsigned int side = longest_side(mesh, triangle);
       requested[3 * i + side] = probe_requires_split(mesh, triangle, side, config);
+      double jacobian;
+      if (discrete_jacobian(mesh, triangle, &jacobian) &&
+          fabs(jacobian) < config->jacobian_minimum)
+        requested[3 * i + side] = 1;
     }
   }
   qsort(edges, edge_count, sizeof *edges, compare_mesh_edge);

src/frame.h

@@ -28,11 +28,15 @@ typedef struct {
   unsigned int max_level;
   double angle_absolute_rad;
   double angle_relative;
+  double jacobian_minimum;
   double min_edge_pixels;
   double min_area_pixels2;
 } RefinementConfig;
 
-typedef enum { FRAME_SAMPLE_VERTEX, FRAME_SAMPLE_PROBE } FrameSampleKind;
+typedef enum {
+  FRAME_SAMPLE_VERTEX,
+  FRAME_SAMPLE_PROBE
+} FrameSampleKind;
 
 typedef struct {
   FrameSampleKind kind;
@@ -74,8 +78,8 @@ 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. */
+/* Builds one generation of vertex/probe ray requests.  Results are installed
+ * only after the caller has completed the generation's sweep. */
 int frame_lens_mesh_prepare_generation(FrameLensMesh *mesh,
                                        const RefinementConfig *config);
 const FrameSample *frame_lens_mesh_samples(const FrameLensMesh *mesh,

src/main.c

@@ -154,6 +154,7 @@ static int parse_args(int argc, char **argv, Settings *s,
                   .refinement = {.angle_absolute_rad =
                                      1e-3 * 3.14159265358979323846 / 180.0,
                                  .angle_relative = 0.1,
+                                 .jacobian_minimum = 1e-3,
                                  .min_edge_pixels = 0.5,
                                  .min_area_pixels2 = 0.25}};
   *write_path = NULL;
@@ -181,6 +182,8 @@ static int parse_args(int argc, char **argv, Settings *s,
       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-jacobian-min") && i + 1 < argc &&
+               !parse_positive(argv[++i], &s->refinement.jacobian_minimum)) {
     } 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 &&
@@ -528,6 +531,7 @@ int main(int argc, char **argv) {
 #endif
             "[--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] "
             "[--draw-mesh] [--write-catalog PATH] "
             "[--catalog-load-workers N] "