Optics: parallelize PSF cache construction

psf_lookup_optimization_plan.md

@@ -116,7 +116,8 @@ case.  The direct evaluator remains a selectable regression reference.
 Implemented as one immutable `PsfKernelCache` per invocation, built after CLI
 parsing and passed read-only through the existing frame splat workers.  The
 cache uses 64 by 64 phase intervals (65 nodes on each axis), 4-point
-Gauss-Legendre pixel-area quadrature, and bilinear phase interpolation.  The
+Gauss-Legendre pixel-area quadrature, bilinear phase interpolation, and a
+static OpenMP phase-plane construction loop.  The
 `--psf-direct` regression mode uses an independent 8-point quadrature
 reference; events whose requested support exceeds the cache use that path
 automatically.

src/optics.c

@@ -2,11 +2,11 @@
 
 #include <limits.h>
 #include <math.h>
+#include <omp.h>
 #include <stdint.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
-#include <time.h>
 
 #ifdef ENABLE_PNG
 #include <png.h>
@@ -213,7 +213,7 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf)
     float *weights = malloc(count * sizeof *weights);
     if (weights == NULL)
         return -1;
-    const clock_t start = clock();
+    const double start = omp_get_wtime();
     PsfKernelCache building = {.fwhm_pixels = psf->fwhm_pixels,
                                .moffat_beta = psf->moffat_beta,
                                .alpha_pixels = alpha,
@@ -221,6 +221,11 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf)
                                .weights = weights,
                                .phase_resolution = PSF_PHASE_RESOLUTION,
                                .radius_pixels = radius};
+    int invalid_kernel = 0;
+    /* Each sub-pixel phase owns a disjoint kernel plane.  Parallelizing this
+     * regular construction needs no locks, and the finished kernel remains
+     * immutable to splat workers. */
+#pragma omp parallel for collapse(2) schedule(static)
     for (int phase_y = 0; phase_y <= PSF_PHASE_RESOLUTION; ++phase_y)
         for (int phase_x = 0; phase_x <= PSF_PHASE_RESOLUTION; ++phase_x) {
             const double star_x = (double)phase_x / PSF_PHASE_RESOLUTION;
@@ -234,10 +239,12 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf)
                         alpha, psf->moffat_beta, offset_x, offset_y, star_x, star_y);
                     weights[index] = (float)weight;
                     sum += weight;
-                }
+            }
             if (!(sum > 0.0) || !isfinite(sum)) {
-                free(weights);
-                return -1;
+                /* An error cannot return from an OpenMP structured block. */
+#pragma omp atomic write
+                invalid_kernel = 1;
+                continue;
             }
             for (int offset_y = -radius; offset_y <= radius; ++offset_y)
                 for (int offset_x = -radius; offset_x <= radius; ++offset_x) {
@@ -246,7 +253,11 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf)
                     weights[index] = (float)(weights[index] / sum);
                 }
         }
-    building.build_seconds = (double)(clock() - start) / CLOCKS_PER_SEC;
+    if (invalid_kernel) {
+        free(weights);
+        return -1;
+    }
+    building.build_seconds = omp_get_wtime() - start;
     building.ready = 1;
     *cache = building;
     return 0;