Fix NVIDIA Dust & Scratches threshold rounding (#211)

Squashed from #211 without its before/after PNG (we don't commit binary images).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

crates/gpu/src/fx_noise.rs

@@ -156,14 +156,17 @@ pub(crate) fn fractal(e: &mut Enc, ctx: &EffectCtx, b: GBuf, turbulent: bool) ->
 /// noise::median_image (of the straight colour with alpha 1 when `straight`).
 pub(crate) fn median_image(e: &mut Enc, img: &GpuImage, r: usize, straight: bool) -> GpuImage {
     let q = run(e, "fxn_quant", &with_u0([straight as u32, 0, 0, 0]), img, None, None);
+    // Preserve the CPU's rounded level / 511 values before the strict Dust & Scratches
+    // threshold. GPU division (even residual-corrected) may differ by one ulp.
+    let levels = e.data(&(0..512).map(|level| level as f32 / 511.0).collect::<Vec<_>>());
     if r <= MEDIAN_BISECT_R {
-        return run(e, "fxn_median", &with_u0([r as u32, 0, 0, 0]), &q, None, None);
+        return run(e, "fxn_median", &with_u0([r as u32, 0, 0, 0]), &q, None, Some(&levels));
     }
     // Segments long enough that filling the first window (2r + 1)² costs about as much as
     // sliding along the rest.
     let seg = (4 * r as u32).clamp(32, 512).min(q.width);
     let out = e.scratch(q.width, q.height);
-    e.dispatch("fxn_median_huang", &with_u0([r as u32, seg, 0, 0]), &q, None, &out, None, (q.height.div_ceil(64), q.width.div_ceil(seg)));
+    e.dispatch("fxn_median_huang", &with_u0([r as u32, seg, 0, 0]), &q, None, &out, Some(&levels), (q.height.div_ceil(64), q.width.div_ceil(seg)));
     out
 }
 

crates/gpu/src/shaders/fx_noise.wgsl

@@ -14,9 +14,9 @@ fn fxn_pixel(gid: vec3<u32>) -> vec2<i32> {
     return p;
 }
 
-// a / b rounded like the CPU's IEEE division (GPU division may be approximate): one
-// residual correction with fma, followed by a second refinement for half-step cases.
-// Used where a strict threshold follows the quotient.
+// Refine approximate GPU division with two fma residual corrections before restoring
+// the median's straight colour. Driver float optimisations can still affect rounding;
+// the fixed median levels instead come from a CPU-rounded lookup table.
 fn fxn_div4(a: vec4<f32>, b: vec4<f32>) -> vec4<f32> {
     let q = a / b;
     let r = fma(-q, b, a);
@@ -53,6 +53,7 @@ fn fxn_quant(@builtin(global_invocation_id) gid: vec3<u32>) {
 
 // Per-channel median of the quantised levels over a (2r+1)² window (edges repeated): the
 // smallest level m with more than half the window at or below m. u[0].x = r.
+// data = the 512 CPU-rounded level / 511 values.
 @compute @workgroup_size(16, 16)
 fn fxn_median(@builtin(global_invocation_id) gid: vec3<u32>) {
     let p = fxn_pixel(gid);
@@ -77,7 +78,8 @@ fn fxn_median(@builtin(global_invocation_id) gid: vec3<u32>) {
         hi = select(hi, mid, above);
         lo = select(mid + 1.0, lo, above);
     }
-    textureStore(out, p, fxn_div4(lo, vec4<f32>(511.0)));
+    let level = vec4<u32>(lo);
+    textureStore(out, p, vec4<f32>(data[level.x], data[level.y], data[level.z], data[level.w]));
 }
 
 // The same median for large radii, as the CPU computes it: a sliding 512-bin histogram per
@@ -85,6 +87,7 @@ fn fxn_median(@builtin(global_invocation_id) gid: vec3<u32>) {
 // the segment's first window, then slides right one column at a time (remove the column
 // leaving, add the one entering) and moves each channel's median by the counts below it.
 // u[0] = (r, segment length); dispatched over (rows, segments) with 64 × 1 workgroups.
+// data = the 512 CPU-rounded level / 511 values.
 @compute @workgroup_size(64, 1)
 fn fxn_median_huang(@builtin(global_invocation_id) gid: vec3<u32>) {
     let dims = out_dims();
@@ -138,7 +141,7 @@ fn fxn_median_huang(@builtin(global_invocation_id) gid: vec3<u32>) {
                 med[c] += 1u;
             }
         }
-        textureStore(out, vec2<i32>(x, y), fxn_div4(vec4<f32>(med), vec4<f32>(511.0)));
+        textureStore(out, vec2<i32>(x, y), vec4<f32>(data[med.x], data[med.y], data[med.z], data[med.w]));
     }
 }
 

crates/gpu/src/tests_fx_noise.rs

@@ -203,6 +203,60 @@ fn echo_and_posterize_time() {
     }
 }
 
+#[test]
+fn median_levels_match_cpu_division() {
+    let Some(g) = crate::tests::gpu() else { return };
+    let mut image = effectcraft_raster::Image::new(512, 1);
+    for (i, pixel) in image.data.iter_mut().enumerate() {
+        *pixel = [i as f32 / 511.0; 4];
+    }
+    let input = g.ctx.upload_image(&image).unwrap();
+    for radius in [3, 5] {
+        // A monotone row's median is its centre, including repeated edges. Exercise both
+        // the bisection and sliding-histogram kernels at every quantised level.
+        let mut enc = crate::context::Enc::new(&g.ctx);
+        let median = crate::fx_noise::median_image(&mut enc, &input, radius, false);
+        let actual = enc.download(&median).unwrap();
+        for (level, (want, got)) in image.data.iter().zip(&actual.data).enumerate() {
+            assert_eq!(got, want, "radius {radius}, level {level}");
+        }
+    }
+}
+
+#[test]
+fn dust_strict_threshold_boundaries() {
+    let Some(g) = crate::tests::gpu() else { return };
+    let spec = effectcraft_effects::find("ec.noise.dustscratches").unwrap();
+    let size = [7.0, 7.0];
+    let mut params =
+        effectcraft_effects::Params { values: spec.params.iter().map(|p| (p.id.to_string(), effectcraft_effects::default_value(p, size))).collect() };
+    let exact = 63.75_f32;
+    let half_colour = [128.0 / 511.0, 128.0 / 511.0, 128.0 / 511.0, 256.0 / 511.0];
+    let cases = [
+        // The NVIDIA failure: rounding 154 / 511 down changes the keep/replace decision.
+        ([126.0 / 511.0, 60.0 / 511.0, 62.0 / 511.0, 154.0 / 511.0], [73.0 / 255.0, 32.0 / 255.0, 32.0 / 255.0, 77.0 / 255.0], 10.0, true),
+        // Candidate red = 0.25, original red = 0.5. Equality must keep the original;
+        // adjacent representable thresholds exercise both sides of the strict comparison.
+        (half_colour, [0.5, 0.25, 0.25, 0.5], exact.next_down(), true),
+        (half_colour, [0.5, 0.25, 0.25, 0.5], exact, false),
+        (half_colour, [0.5, 0.25, 0.25, 0.5], exact.next_up(), false),
+    ];
+    for radius in [3.0, 5.0] {
+        params.values.insert("radius".into(), n(radius));
+        for (median, original, threshold, replaced) in cases {
+            params.values.insert("threshold".into(), n(threshold as f64));
+            let ctx = || effectcraft_effects::EffectCtx { params: &params, time: 0.0, layer_size: size, seed: 0, adjustment: false, env: Default::default() };
+            let mut img = effectcraft_raster::Image::filled(7, 7, median);
+            img.set(3, 3, original);
+            let buf = effectcraft_effects::Buf { img, offset: [0.0; 2], scale: 1.0 };
+            let cpu = (spec.render)(&ctx(), buf.clone());
+            assert_eq!(cpu.img.data[24] != original, replaced, "CPU radius {radius}, threshold {threshold}");
+            let gpu = effectcraft_render::Accelerator::effects(g, &[effectcraft_render::FxStep { spec, ctx: ctx() }], &buf, None).unwrap();
+            assert_eq!(gpu.img.data[24], cpu.img.data[24], "radius {radius}, threshold {threshold}");
+        }
+    }
+}
+
 #[test]
 fn dust_threshold_on_quantized_input() {
     let Some(g) = crate::tests::gpu() else { return };