Spektrafilm: More bugfixing

data/kernels/spektrafilm.cl

@@ -248,27 +248,64 @@ static float3 sf_pchip3d(__global const float *lut, __global const float *sx,
 
 /* ---- base-10 <-> base-2, matching the host's SF_LOG10F/SF_POW10F -------- */
 /* spektra_sim.c does NOT call log10f()/exp10f(). It defines
-     SF_LOG10F(x) = log2f(x) * 0.3010299956639812f
-     SF_POW10F(x) = exp2f(x * 3.321928094887362f)
-   so calling log10()/exp10() here was not a more-or-less accurate version of
-   the same computation, it was a different computation: the host's scaling
-   multiply carries its own rounding that exp10()/log10() never perform. Use
-   the identical formulation so only the library ULP gap remains.
-
-   That remaining gap is real but small: OpenCL specs log2/exp2 to <=3 ULP and
-   most GPUs implement both in hardware, while glibc's log2f/exp2f are
-   correctly rounded. If the integration test still shows a residue after
-   this, these two are the next candidates for the portable-polynomial
-   treatment sf_exp_neg already got -- log2 via exponent extraction plus a
-   mantissa polynomial, exp2 via sf_exp2i plus a fractional polynomial. */
+     SF_LOG10F(x) = sf_log2f(x) * 0.3010299956639812f
+     SF_POW10F(x) = sf_exp2f(x * 3.321928094887362f)
+   so calling log10()/exp10() here would not be a more-or-less accurate
+   version of the same computation, it would be a different one: the host's
+   scaling multiply carries its own rounding that exp10()/log10() never
+   perform.
+
+   log2()/exp2() are not usable either, for the reason that governs
+   sf_exp_neg: OpenCL specs both to <=3 ULP and most GPUs implement them in
+   hardware, while glibc rounds correctly. That gap lands in the film density
+   feeding sf_layer_particle, and sf_poisson's accept/reject loop below turns
+   it into whole-integer grain draws on isolated pixels. sf_exp2f/sf_log2f
+   below are the same portable polynomials as spektra_core.h's, built only
+   from correctly-rounded operations, so both sides agree bit-for-bit. Keep
+   them in lockstep with the host copies -- same constants, same order of
+   operations -- and do not substitute the library or native_ variants. */
+/* defined below, next to sf_exp_neg, which is the other user */
+static inline float sf_exp2i(int k);
+
+static inline float sf_exp2f(float x)
+{
+  const int k = (int)floor(x + 0.5f);
+  const float t = (x - (float)k) * 0.6931471824645996f; /* ln(2) */
+  float p = 0.000198412700f;                            /* 1/5040 */
+  p = p * t + 0.00138888892f;                           /* 1/720 */
+  p = p * t + 0.00833333377f;                           /* 1/120 */
+  p = p * t + 0.0416666679f;                            /* 1/24 */
+  p = p * t + 0.166666672f;                             /* 1/6 */
+  p = p * t + 0.5f;
+  p = p * t + 1.0f;
+  p = p * t + 1.0f;
+  return p * sf_exp2i(k);
+}
+
+static inline float sf_log2f(float x)
+{
+  const uint xu = as_uint(x);
+  int e = (int)((xu >> 23) & 0xffu) - 127;
+  float m = as_float((xu & 0x007fffffu) | 0x3f800000u);
+  if(m > 1.41421356f) { m *= 0.5f; e += 1; }
+  const float s = (m - 1.0f) / (m + 1.0f);
+  const float s2 = s * s;
+  float p = 0.222222224f;    /* 2/9 */
+  p = p * s2 + 0.285714298f; /* 2/7 */
+  p = p * s2 + 0.400000006f; /* 2/5 */
+  p = p * s2 + 0.666666687f; /* 2/3 */
+  p = p * s2 + 2.0f;
+  return (float)e + p * s * 1.4426950216293335f; /* log2(e) */
+}
+
 static inline float sf_log10f(float x)
 {
-  return log2(x) * 0.3010299956639812f;
+  return sf_log2f(x) * 0.3010299956639812f;
 }
 
 static inline float sf_pow10f(float x)
 {
-  return exp2(x * 3.321928094887362f);
+  return sf_exp2f(x * 3.321928094887362f);
 }
 
 /* ---- [gc] Reinhard knee + OkLCh output gamut compression ---------------- */

src/common/spektra_core.h

@@ -349,6 +349,68 @@ SPEKTRA_INLINE float sf_exp_neg(float lam)
   return p * sf_exp2i(k);
 }
 
+/* sf_exp2f / sf_log2f: 2^x and log2(x) built from +, -, *, / and the exact
+   floor/exponent manipulation above, for the reason sf_exp_neg is -- and
+   this pair is what actually reaches the grain sampler. SF_POW10F and
+   SF_LOG10F are defined in terms of these rather than the platform
+   exp2f/log2f, which OpenCL specifies only to <=3 ULP while glibc rounds
+   correctly: that slack lands in the film density arriving at
+   sf_layer_particle, and sf_poisson's accept/reject loop turns a one-ULP
+   density difference into a whole-integer grain count difference wherever a
+   partial product happens to sit near limit. The result is isolated pixels,
+   scattered evenly and uncorrelated with image structure, each off by a full
+   grain quantum rather than by a rounding.
+
+   Every operation here is correctly rounded on both sides by IEEE-754 and by
+   the OpenCL spec, so the same x yields the same bits. The same caveat as
+   sf_exp_neg applies: none of these multiply-adds may be contracted into an
+   FMA on one side and not the other, which is what -ffp-contract=off on the
+   host and #pragma OPENCL FP_CONTRACT OFF in the kernel are for.
+
+   Both are ~1 ULP against glibc over the domains this module uses, and are
+   not general-purpose replacements outside them: sf_exp2f assumes the result
+   stays normal, sf_log2f assumes x is positive and normal. SF_LOG10F floors
+   its argument at SF_LOG_EPS, which keeps it there. */
+SPEKTRA_INLINE float sf_exp2f(float x)
+{
+  /* x = k + r, k integer and |r| <= 0.5, so 2^x = 2^k * e^(r ln2) with the
+     exponential taken over |t| <= 0.347 by a degree-7 Taylor polynomial in
+     Horner form and 2^k injected exactly. */
+  const int k = (int)floorf(x + 0.5f);
+  const float t = (x - (float)k) * 0.6931471824645996f; /* ln(2) */
+  float p = 0.000198412700f;                            /* 1/5040 */
+  p = p * t + 0.00138888892f;                           /* 1/720 */
+  p = p * t + 0.00833333377f;                           /* 1/120 */
+  p = p * t + 0.0416666679f;                            /* 1/24 */
+  p = p * t + 0.166666672f;                             /* 1/6 */
+  p = p * t + 0.5f;
+  p = p * t + 1.0f;
+  p = p * t + 1.0f;
+  return p * sf_exp2i(k);
+}
+
+SPEKTRA_INLINE float sf_log2f(float x)
+{
+  /* Take the binary exponent off by hand, then fold the mantissa into
+     [1/sqrt2, sqrt2] so that s = (m-1)/(m+1) stays inside +-0.1716, where
+     log(m) = 2(s + s^3/3 + s^5/5 + s^7/7 + s^9/9) is good to well under an
+     ULP. Both steps of the fold are exact. */
+  union { float f; uint32_t u; } v;
+  v.f = x;
+  int e = (int)((v.u >> 23) & 0xffu) - 127;
+  v.u = (v.u & 0x007fffffu) | 0x3f800000u;
+  float m = v.f;
+  if(m > 1.41421356f) { m *= 0.5f; e += 1; }
+  const float s = (m - 1.0f) / (m + 1.0f);
+  const float s2 = s * s;
+  float p = 0.222222224f;    /* 2/9 */
+  p = p * s2 + 0.285714298f; /* 2/7 */
+  p = p * s2 + 0.400000006f; /* 2/5 */
+  p = p * s2 + 0.666666687f; /* 2/3 */
+  p = p * s2 + 2.0f;
+  return (float)e + p * s * 1.4426950216293335f; /* log2(e) */
+}
+
 SPEKTRA_INLINE float sf_poisson(float lam,
                                 uint32_t seed)
 {

src/common/spektra_fetch.c

@@ -98,6 +98,7 @@ static struct
   double progress;
   char message[256];
   guint generation; /* bumped whenever an install changes what is on disk */
+  guint finished_idle; /* pending _finished_idle source, 0 when none */
   gboolean inited;
 } _sf = { .state = SF_FETCH_IDLE };
 
@@ -114,6 +115,7 @@ static void _set_status(const sf_fetch_state_t state,
                         const double progress,
                         const char *msg)
 {
+  if(!_sf.inited) return;
   g_mutex_lock(&_sf.lock);
   _sf.state = state;
   if(progress >= 0.0) _sf.progress = CLAMP(progress, 0.0, 1.0);
@@ -123,6 +125,11 @@ static void _set_status(const sf_fetch_state_t state,
 
 static gboolean _cancelled(void)
 {
+  /* Not started, or already torn down: nothing is in flight to cancel, and the
+     mutex below may have been cleared. Every entry point that takes _sf.lock
+     tests this first -- the lock only exists between sf_fetch_init() and
+     sf_fetch_cleanup(), and callers on the GUI side outlive neither reliably. */
+  if(!_sf.inited) return FALSE;
   g_mutex_lock(&_sf.lock);
   const gboolean c = _sf.cancel;
   g_mutex_unlock(&_sf.lock);
@@ -150,6 +157,17 @@ void sf_fetch_cleanup(void)
   _sf.thread = NULL;
   g_mutex_unlock(&_sf.lock);
   if(t) g_thread_join(t);
+
+  /* The worker's last act is to post _finished_idle, so joining guarantees it
+     has been posted, not that it has run. It calls sf_fetch_status(), so a
+     dispatch after the clear below would take a cleared mutex; drop it while
+     the lock is still valid. Zero means it already ran. */
+  g_mutex_lock(&_sf.lock);
+  const guint idle = _sf.finished_idle;
+  _sf.finished_idle = 0;
+  g_mutex_unlock(&_sf.lock);
+  if(idle) g_source_remove(idle);
+
   g_mutex_clear(&_sf.lock);
   _sf.inited = FALSE;
 }
@@ -158,6 +176,12 @@ sf_fetch_state_t sf_fetch_status(char *msg,
                                  size_t msgsz,
                                  double *progress)
 {
+  if(!_sf.inited)
+  {
+    if(msg && msgsz) msg[0] = '\0';
+    if(progress) *progress = 0.0;
+    return SF_FETCH_IDLE;
+  }
   g_mutex_lock(&_sf.lock);
   const sf_fetch_state_t s = _sf.state;
   if(msg && msgsz) g_strlcpy(msg, _sf.message, msgsz);
@@ -168,6 +192,7 @@ sf_fetch_state_t sf_fetch_status(char *msg,
 
 void sf_fetch_cancel(void)
 {
+  if(!_sf.inited) return;
   g_mutex_lock(&_sf.lock);
   _sf.cancel = TRUE;
   g_mutex_unlock(&_sf.lock);
@@ -878,6 +903,15 @@ typedef struct sf_worker_args_t
 static gboolean _finished_idle(gpointer user_data)
 {
   const gboolean ok = GPOINTER_TO_INT(user_data);
+
+  /* Taken before anything else, and posted under the same lock: if the main
+     loop reaches this before the worker has recorded the source id, that store
+     is still in progress and this blocks until it lands, so the id can never be
+     left behind pointing at a source that has already run. */
+  g_mutex_lock(&_sf.lock);
+  _sf.finished_idle = 0;
+  g_mutex_unlock(&_sf.lock);
+
   char msg[256] = { 0 };
   sf_fetch_status(msg, sizeof(msg), NULL);
 
@@ -901,6 +935,7 @@ static gpointer _fetch_worker(gpointer data)
   gboolean success = FALSE;
   char *repo = NULL, *ref = NULL, *manifest_url = NULL, *manifest = NULL;
   char *base = NULL, *tmpdir = NULL, *destdir = NULL, *profdir = NULL;
+  char *olddir = NULL;
   GPtrArray *files = NULL;
   CURL *curl = NULL;
 
@@ -1034,11 +1069,45 @@ static gpointer _fetch_worker(gpointer data)
     goto out;
   }
 
-  _rmdir_recursive(destdir); /* replacing an older copy of the same hash */
+  /* Replacing an older copy of the same hash: move it aside rather than delete
+     it, so there is never a moment with no pack at this path. Deleting first
+     and then failing the rename -- out of space, a permission change, a handle
+     held open on the file the user is mid-render against -- leaves the user
+     with nothing, because the out: block below then removes the download too.
+     The name is dotted like .incoming, and lookups address packs by their exact
+     %08x name rather than scanning, so neither is ever mistaken for a pack. */
+  if(g_file_test(destdir, G_FILE_TEST_IS_DIR))
+  {
+    olddir = g_strdup_printf("%s%s.replaced-%08x", packs, G_DIR_SEPARATOR_S, got_hash);
+    _rmdir_recursive(olddir);
+    if(g_rename(destdir, olddir) != 0)
+    {
+      /* Cannot move it aside, so it could not be restored either. Stop here
+         and keep it: a pack that already renders is worth more than the one
+         being installed, and deleting it to make room would risk ending up
+         with neither. */
+      dt_print(DT_DEBUG_ALWAYS,
+               "[spektrafilm] cannot move the installed pack aside at %s: %s",
+               destdir, strerror(errno));
+      g_free(olddir);
+      olddir = NULL;
+      _set_status(SF_FETCH_FAILED, -1.0,
+                  _("could not replace the installed pack -- the existing one "
+                    "has been kept"));
+      goto out;
+    }
+  }
+
   if(g_rename(tmpdir, destdir) != 0)
   {
     dt_print(DT_DEBUG_ALWAYS, "[spektrafilm] cannot install pack into %s: %s",
              destdir, strerror(errno));
+    /* Put the working pack back before reporting the failure. */
+    if(olddir && g_rename(olddir, destdir) == 0)
+    {
+      g_free(olddir);
+      olddir = NULL;
+    }
     _set_status(SF_FETCH_FAILED, -1.0, _("could not install the downloaded pack"));
     goto out;
   }
@@ -1055,6 +1124,9 @@ static gpointer _fetch_worker(gpointer data)
 
 out:
   if(!success && tmpdir) _rmdir_recursive(tmpdir);
+  /* Reached with olddir set only once the new pack is in place, or after a
+     restore that itself failed; either way the copy it names is superseded. */
+  if(olddir) _rmdir_recursive(olddir);
   if(curl) curl_easy_cleanup(curl);
   _files_free(files);
   g_free(manifest);
@@ -1065,13 +1137,13 @@ static gpointer _fetch_worker(gpointer data)
   g_free(tmpdir);
   g_free(destdir);
   g_free(profdir);
+  g_free(olddir);
 
   g_mutex_lock(&_sf.lock);
   _sf.cancel = FALSE;
-  g_mutex_unlock(&_sf.lock);
-
   if(darktable.gui)
-    g_idle_add(_finished_idle, GINT_TO_POINTER(success ? 1 : 0));
+    _sf.finished_idle = g_idle_add(_finished_idle, GINT_TO_POINTER(success ? 1 : 0));
+  g_mutex_unlock(&_sf.lock);
 
   return NULL;
 }

src/common/spektra_sim.c

@@ -58,10 +58,6 @@
 #include "spektra_core.h"
 
 #include <math.h>
-/* ensure C99 math functions for SF_POW10F/SF_LOG10F (exp2f, log2f) */
-#if !defined(exp2f) && !defined(_GNU_SOURCE)
-/* exp2f and log2f are C99; every compiler since GCC 4.x / Clang 3.x has them */
-#endif
 #include <stdint.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -93,11 +89,15 @@ static inline void neon_mat3_mulv_batch(const float m[9],
 }
 #endif /* __ARM_NEON */
 
-/* Fast pow10 / log10 via exp2f/log2f. Using compiler builtins gives the
-   optimizer a better chance to inline/reduce them vs libm powf(10, x)
-   which internally computes exp2(x * log2(10)) with extra overhead. */
-#define SF_POW10F(x) __builtin_exp2f((x) * 3.321928094887362f)  /* x * log2f(10) */
-#define SF_LOG10F(x) (__builtin_log2f(x) * 0.3010299956639812f)  /* log2f(x) * log10(2) */
+/* pow10 / log10 through spektra_core.h's own exp2/log2, NOT the platform
+   exp2f/log2f. The kernel has to compute the same thing, and OpenCL specifies
+   exp2/log2 only to <=3 ULP where glibc rounds correctly, so a library call
+   here is a library call the GPU cannot match. sf_exp2f/sf_log2f are built
+   from correctly-rounded operations alone and agree bit-for-bit on both
+   sides; see their comment for why a ULP here is not a ULP by the time it
+   reaches the grain sampler. */
+#define SF_POW10F(x) sf_exp2f((x) * 3.321928094887362f)  /* x * log2f(10) */
+#define SF_LOG10F(x) (sf_log2f(x) * 0.3010299956639812f)  /* log2f(x) * log10(2) */
 #define SF_TC_KNEE_T 0.0 /* [gc] InputGamutCompressSpec.knee */
 #define SF_TC_KNEE_L 1.0
 #define SF_TC_KNEE_P 6.0
@@ -1320,6 +1320,7 @@ void sf_sim_params_defaults(sf_sim_params_t *p)
   p->grain_rms_scale = -1.0;
   p->grain_uniformity_scale = -1.0;
   p->grain_particle_scale = -1.0;
+  p->grain_density_min_scale = -1.0;
   p->coupler_diffusion_um = -1.0;
   p->coupler_tail_um = -1.0;
   p->coupler_tail_weight = -1.0;
@@ -3447,8 +3448,26 @@ sf_sim_t *sf_sim_build(const sf_pack_t *pack,
     if(p->grain_uniformity_scale >= 0.0)
       for(int c = 0; c < 3; c++)
         s->grain_uniformity[c] = fmin(s->grain_uniformity[c] * p->grain_uniformity_scale, 0.999);
+    /* A scale and not an absolute value, for the same reason rms and uniformity
+       are: the pack's floors are per channel -- kodak_vision3_500t is
+       0.12/0.10/0.35 -- and one number replacing all three would flatten a
+       shape that carries real colour information. Scaling keeps the stock's own
+       ratios and still reaches any overall floor. Applied after the pack read,
+       so it lands on the film's own value; sf_pack_film_grain() leaves
+       p->grain_density_min alone for a stock it does not characterise, and the
+       scale then multiplies the caller's fallback instead. */
+    if(p->grain_density_min_scale >= 0.0)
+      for(int c = 0; c < 3; c++) p->grain_density_min[c] *= p->grain_density_min_scale;
+    /* Sub-layer 0 is the coarsest and stays the reference at 1.0: this control
+       moves the FINER sub-layers relative to it, so it starts at i == 1.
+       Scaling the whole array instead is an exact no-op -- _sf_build_grain_layers
+       derives a_coarsest as sig^2*A48/peak[c], peak[c] is linear in
+       particle_scale[], and particle_area = a_coarsest * particle_scale[l], so a
+       common factor k cancels and every layer_npart comes out unchanged. The
+       only value that did anything was exactly 0, where the 1e-9 floors on peak
+       and particle_area take over and npart explodes, i.e. grain disappears. */
     if(p->grain_particle_scale >= 0.0)
-      for(int i = 0; i < n_scale; i++) particle_scale[i] *= p->grain_particle_scale;
+      for(int i = 1; i < n_scale; i++) particle_scale[i] *= p->grain_particle_scale;
     _sf_build_grain_layers(s, film, p->grain_density_min, s->grain_uniformity,
                            s->grain_rms, particle_scale, n_scale);
   }