Optimize Radial Blur sampling without changing output (#1499)

Radial Blur prepares its sample transforms once per filter (blur::radial::Plan) and shares them across output tiles, reading each bilinear tap once with channels kept together. Output is identical to the previous sampler at Draft, Good and Best, across colour models and depths (tests compare against the unchanged reference path); a failed table allocation falls back to that sampler. Measured on 24 MP RGBA8, Good, amount 1: Spin 66.1 -> 16.4 s, Zoom 19.2 -> 5.4 s (bench_radial example, docs/radial-blur-performance.md).

Co-Authored-By: Claude Opus 5.5 <[email protected]>

crates/algo/examples/bench_radial.rs

@@ -0,0 +1,116 @@
+//! Paired, end-to-end Radial Blur benchmark against the current upstream sampler.
+//! cargo run --release -j 12 -p photocraft-algo --example bench_radial
+//! Defaults: 6000x4000 RGBA8, Good, amount 100, 12 Rayon workers, three runs.
+//! Optional positional arguments: width height amount repeats.
+
+#[cfg(not(target_arch = "wasm32"))]
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    use photocraft_algo::{Ctx, FilterParams, Image, RadialMethod, RadialQuality, TILE, apply_in, kernel};
+    use photocraft_color::PixelFormat;
+    use photocraft_geom::Rect;
+    use photocraft_raster::Surface;
+    use rayon::prelude::*;
+    use std::time::Instant;
+
+    // Match apply_tiled_with's current upstream pipeline for this benchmark's
+    // full-canvas, unselected global filter. kernel retains the upstream sampler.
+    // Both paths read one normalized source, batch the same 256px tiles within
+    // the same 256MiB result budget, write the same surface depth, and prune.
+    fn upstream(source: &Surface, params: &FilterParams, rect: Rect) -> Surface {
+        let mut output = source.clone();
+        let image = Image::read(source, rect);
+        let format = source.format();
+        let ctx = Ctx { bounds: rect, mode: format.mode, alpha: format.alpha };
+        let tiles: Vec<Rect> = rect.tiles().map(|tile| tile.rect().intersect(&rect)).collect();
+        let per_tile = TILE as usize * TILE as usize * format.channels() * size_of::<f32>();
+        let group = ((256 << 20) / per_tile).max(rayon::current_num_threads());
+        for chunk in tiles.chunks(group) {
+            let results: Vec<_> = chunk.par_iter().map(|tile| (*tile, kernel(params, &image, *tile, &ctx))).collect();
+            for (tile, samples) in results {
+                output.write_region(tile, &samples);
+            }
+        }
+        output.prune();
+        output
+    }
+
+    let args: Vec<_> = std::env::args().skip(1).collect();
+    let parse = |i: usize, default: i32| -> Result<i32, Box<dyn std::error::Error>> {
+        Ok(match args.get(i) {
+            Some(value) => value.parse()?,
+            None => default,
+        })
+    };
+    let width = parse(0, 6000)?;
+    let height = parse(1, 4000)?;
+    let amount = parse(2, 100)?;
+    let repeats = parse(3, 3)?;
+    if !(1..=10000).contains(&width) || !(1..=10000).contains(&height) || !(0..=100).contains(&amount) || !(1..=7).contains(&repeats) {
+        return Err("expected dimensions 1..10000, amount 0..100, repeats 1..7".into());
+    }
+    let pool = rayon::ThreadPoolBuilder::new().num_threads(12).build()?;
+    let rect = Rect::new(0, 0, width, height);
+    let mut source = Surface::new(PixelFormat::RGBA8);
+    for tile in rect.tiles() {
+        let area = tile.rect().intersect(&rect);
+        let mut values = Vec::new();
+        values.try_reserve_exact(area.width() as usize * area.height() as usize * 4)?;
+        for y in area.y0..area.y1 {
+            for x in area.x0..area.x1 {
+                let seed = (x as u32).wrapping_mul(0x9e37_79b9) ^ (y as u32).wrapping_mul(0x85eb_ca6b);
+                for channel in 0..4 {
+                    let level = seed.rotate_left(channel * 7).wrapping_add(channel * 71) % 256;
+                    values.push(if channel == 3 && seed.is_multiple_of(5) { 0.0 } else { level as f32 / 255.0 });
+                }
+            }
+        }
+        source.write_region(area, &values);
+    }
+    pool.install(|| -> Result<(), Box<dyn std::error::Error + Send + Sync>> {
+        // Small untimed warmup initializes worker/code paths without doubling
+        // the cost of each 24MP measurement. Full-size runs include page faults.
+        let warm = Rect::new(0, 0, width.min(128), height.min(128));
+        for method in [RadialMethod::Spin, RadialMethod::Zoom] {
+            let params = FilterParams::RadialBlur { amount: amount as f32, method, quality: RadialQuality::Good, center_x: 0.5, center_y: 0.5 };
+            upstream(&source, &params, warm);
+            apply_in(&source, &params, warm, warm, None, warm);
+            let (mut before, mut after) = (Vec::new(), Vec::new());
+            for run in 0..repeats {
+                let execute = |original: bool| {
+                    let start = Instant::now();
+                    let output = if original { upstream(&source, &params, rect) } else { apply_in(&source, &params, rect, rect, None, rect) };
+                    (output, start.elapsed().as_secs_f64())
+                };
+                let ((original, old_time), (optimized, new_time)) = if run % 2 == 0 {
+                    let first = execute(true);
+                    (first, execute(false))
+                } else {
+                    let first = execute(false);
+                    (execute(true), first)
+                };
+                // Exact comparison after timing, including alpha and native U8
+                // writeback; benchmark errors never silently publish a speedup.
+                for tile in rect.tiles() {
+                    let area = tile.rect().intersect(&rect);
+                    if original.read_region(area) != optimized.read_region(area) {
+                        return Err(format!("output mismatch for {method:?} at {area:?}").into());
+                    }
+                }
+                before.push(old_time);
+                after.push(new_time);
+                println!("{method:?} {width}x{height} Good amount={amount} run={}: upstream={old_time:.6}s optimized={new_time:.6}s exact=true", run + 1);
+            }
+            before.sort_by(f64::total_cmp);
+            after.sort_by(f64::total_cmp);
+            let i = before.len() / 2;
+            let old = before.get(i).copied().ok_or("missing upstream timing")?;
+            let new = after.get(i).copied().ok_or("missing optimized timing")?;
+            println!("MEDIAN {method:?}: upstream={old:.6}s optimized={new:.6}s speedup={:.3}x exact=true workers=12 repeats={repeats}", old / new);
+        }
+        Ok(())
+    })
+    .map_err(|error| -> Box<dyn std::error::Error> { error })
+}
+
+#[cfg(target_arch = "wasm32")]
+fn main() {}

crates/algo/src/blur.rs

@@ -18,6 +18,7 @@ mod motion_fft;
 mod motion_fft_tests;
 #[cfg(test)]
 mod motion_tests;
+pub(crate) mod radial;
 
 /// Normalized Gaussian kernel with standard deviation `sigma` (radius 3σ).
 pub(crate) fn gaussian_kernel(sigma: f32) -> Vec<f32> {

crates/algo/src/blur/radial.rs

@@ -0,0 +1,355 @@
+//! Direct Radial Blur, with immutable sampling tables shared by every output tile.
+//!
+//! Pixel positions and the quality-dependent interval ceilings match the upstream
+//! sampler. The hot loop keeps channels together, reads each bilinear tap once,
+//! and preserves the upstream alpha arithmetic for every sample.
+
+use photocraft_geom::Rect;
+use photocraft_raster::Interrupt;
+
+use crate::image::Image;
+use crate::{Ctx, RadialMethod, RadialQuality};
+
+/// Position-independent rotations (Spin) or ray scales (Zoom), prepared once.
+pub(crate) struct Plan {
+    center: (f32, f32),
+    extent: f32,
+    method: RadialMethod,
+    // Packed tables for interval counts 1..=max_intervals, each with n+1 samples.
+    // Best needs at most 8,394,752 pairs (about 64 MiB), reserved fallibly.
+    tables: Vec<(f32, f32)>,
+    max_intervals: usize,
+    valid: bool,
+}
+
+impl Plan {
+    #[cfg(test)]
+    pub(crate) fn new(bounds: Rect, amount: f32, method: RadialMethod, center: (f32, f32)) -> Self {
+        Self::with_quality(bounds, amount, method, RadialQuality::Draft, center, &Interrupt::NONE).unwrap_or(Self {
+            center: (0.0, 0.0),
+            extent: 0.0,
+            method,
+            tables: Vec::new(),
+            max_intervals: 64,
+            valid: false,
+        })
+    }
+
+    #[cfg(test)]
+    pub(crate) fn with_quality(
+        bounds: Rect,
+        amount: f32,
+        method: RadialMethod,
+        quality: RadialQuality,
+        center: (f32, f32),
+        ctl: &Interrupt<'_>,
+    ) -> Option<Self> {
+        Self::prepare(bounds, amount, method, quality, center, None, ctl)
+    }
+
+    pub(crate) fn for_area(
+        bounds: Rect,
+        amount: f32,
+        method: RadialMethod,
+        quality: RadialQuality,
+        center: (f32, f32),
+        area: Rect,
+        ctl: &Interrupt<'_>,
+    ) -> Option<Self> {
+        Self::prepare(bounds, amount, method, quality, center, Some(area), ctl)
+    }
+
+    #[allow(clippy::too_many_arguments)]
+    fn prepare(
+        bounds: Rect,
+        amount: f32,
+        method: RadialMethod,
+        quality: RadialQuality,
+        center: (f32, f32),
+        area: Option<Rect>,
+        ctl: &Interrupt<'_>,
+    ) -> Option<Self> {
+        if !amount.is_finite() || !center.0.is_finite() || !center.1.is_finite() || ctl.cancelled() {
+            return None;
+        }
+        let amount = amount.clamp(0.0, 100.0);
+        let center = (bounds.x0 as f32 + bounds.width() as f32 * center.0, bounds.y0 as f32 + bounds.height() as f32 * center.1);
+        if !center.0.is_finite() || !center.1.is_finite() {
+            return None;
+        }
+        let extent = match method {
+            RadialMethod::Spin => amount.to_radians(),
+            RadialMethod::Zoom => amount / 200.0,
+        };
+        // Only prepare interval counts reachable in the requested output. This
+        // avoids preparing all Best tables (64MiB) for a tiny image/selection.
+        // Distance from a fixed center is greatest at a rectangle's corners.
+        let path = area.map_or(f32::INFINITY, |area| {
+            let mut radius = 0.0f32;
+            for y in [area.y0, area.y1.saturating_sub(1)] {
+                for x in [area.x0, area.x1.saturating_sub(1)] {
+                    let dx = x as f32 + 0.5 - center.0;
+                    let dy = y as f32 + 0.5 - center.1;
+                    radius = radius.max((dx * dx + dy * dy).sqrt());
+                }
+            }
+            extent * radius
+        });
+        let max_intervals = super::radial_intervals(path, quality);
+        let count = max_intervals.checked_mul(max_intervals.checked_add(3)?)?.checked_div(2)?;
+        let mut tables = Vec::new();
+        tables.try_reserve_exact(count).ok()?;
+        for n in 1..=max_intervals {
+            if ctl.cancelled() {
+                return None;
+            }
+            for i in 0..=n {
+                tables.push(match method {
+                    RadialMethod::Spin => ((i as f32 / n as f32 - 0.5) * extent).sin_cos(),
+                    RadialMethod::Zoom => (0.0, 1.0 - extent * i as f32 / n as f32),
+                });
+            }
+        }
+        Some(Self { center, extent, method, tables, max_intervals, valid: true })
+    }
+
+    /// `None` on cancellation or invalid image storage or allocation failure. Sampling beyond
+    /// `src.rect` is transparent, including the individual taps at its edges.
+    pub(crate) fn filter(&self, src: &Image, out: Rect, ctx: &Ctx, ctl: &Interrupt<'_>) -> Option<Vec<f32>> {
+        if !self.valid || ctl.cancelled() || src.ch == 0 {
+            return None;
+        }
+        let source = Source::new(src)?;
+        let len = sample_len(out, src.ch)?;
+        let mut result = Vec::new();
+        result.try_reserve_exact(len).ok()?;
+        result.resize(len, 0.0);
+        if out.is_empty() {
+            return Some(result);
+        }
+        match src.ch {
+            1 => self.dispatch_fixed::<1>(&source, out, ctx.alpha, ctl, &mut result),
+            2 => self.dispatch_fixed::<2>(&source, out, ctx.alpha, ctl, &mut result),
+            3 => self.dispatch_fixed::<3>(&source, out, ctx.alpha, ctl, &mut result),
+            4 => self.dispatch_fixed::<4>(&source, out, ctx.alpha, ctl, &mut result),
+            5 => self.dispatch_fixed::<5>(&source, out, ctx.alpha, ctl, &mut result),
+            _ => self.dispatch_generic(&source, out, ctx.alpha, ctl, &mut result),
+        }?;
+        Some(result)
+    }
+
+    fn dispatch_fixed<const N: usize>(&self, src: &Source<'_>, out: Rect, alpha: bool, ctl: &Interrupt<'_>, dst: &mut [f32]) -> Option<()> {
+        match (alpha, self.method) {
+            (true, RadialMethod::Spin) => self.fixed::<N, true, true>(src, out, ctl, dst),
+            (false, RadialMethod::Spin) => self.fixed::<N, false, true>(src, out, ctl, dst),
+            (true, RadialMethod::Zoom) => self.fixed::<N, true, false>(src, out, ctl, dst),
+            (false, RadialMethod::Zoom) => self.fixed::<N, false, false>(src, out, ctl, dst),
+        }
+    }
+
+    fn fixed<const N: usize, const ALPHA: bool, const SPIN: bool>(&self, src: &Source<'_>, out: Rect, ctl: &Interrupt<'_>, dst: &mut [f32]) -> Option<()> {
+        let row_len = (out.width() as usize).checked_mul(N)?;
+        for (row_index, row) in dst.chunks_exact_mut(row_len).enumerate() {
+            if row_index.is_multiple_of(8) && ctl.cancelled() {
+                return None;
+            }
+            let dy = (i64::from(out.y0) + row_index as i64) as f32 + 0.5 - self.center.1;
+            for (column, pixel) in row.as_chunks_mut::<N>().0.iter_mut().enumerate() {
+                let dx = (i64::from(out.x0) + column as i64) as f32 + 0.5 - self.center.0;
+                let table = self.table(dx, dy)?;
+                let mut acc = [0.0; N];
+                for &tap in table {
+                    let (x, y) = self.position::<SPIN>(dx, dy, tap);
+                    sample_fixed::<N, ALPHA>(src, x, y, &mut acc);
+                }
+                finish(&mut acc, 1.0 / table.len() as f32, ALPHA);
+                pixel.copy_from_slice(&acc);
+            }
+        }
+        Some(())
+    }
+
+    fn dispatch_generic(&self, src: &Source<'_>, out: Rect, alpha: bool, ctl: &Interrupt<'_>, dst: &mut [f32]) -> Option<()> {
+        match (alpha, self.method) {
+            (true, RadialMethod::Spin) => self.generic::<true, true>(src, out, ctl, dst),
+            (false, RadialMethod::Spin) => self.generic::<false, true>(src, out, ctl, dst),
+            (true, RadialMethod::Zoom) => self.generic::<true, false>(src, out, ctl, dst),
+            (false, RadialMethod::Zoom) => self.generic::<false, false>(src, out, ctl, dst),
+        }
+    }
+
+    fn generic<const ALPHA: bool, const SPIN: bool>(&self, src: &Source<'_>, out: Rect, ctl: &Interrupt<'_>, dst: &mut [f32]) -> Option<()> {
+        let n = src.image.ch;
+        let row_len = (out.width() as usize).checked_mul(n)?;
+        let mut sample = Vec::new();
+        sample.try_reserve_exact(n).ok()?;
+        sample.resize(n, 0.0);
+        for (row_index, row) in dst.chunks_exact_mut(row_len).enumerate() {
+            if row_index.is_multiple_of(8) && ctl.cancelled() {
+                return None;
+            }
+            let dy = (i64::from(out.y0) + row_index as i64) as f32 + 0.5 - self.center.1;
+            for (column, pixel) in row.chunks_exact_mut(n).enumerate() {
+                let dx = (i64::from(out.x0) + column as i64) as f32 + 0.5 - self.center.0;
+                let table = self.table(dx, dy)?;
+                for &tap in table {
+                    let (x, y) = self.position::<SPIN>(dx, dy, tap);
+                    sample_generic::<ALPHA>(src, x, y, pixel, &mut sample);
+                }
+                finish(pixel, 1.0 / table.len() as f32, ALPHA);
+            }
+        }
+        Some(())
+    }
+
+    #[inline]
+    fn table(&self, dx: f32, dy: f32) -> Option<&[(f32, f32)]> {
+        let radius = (dx * dx + dy * dy).sqrt();
+        let n = ((self.extent * radius).ceil() as usize).clamp(1, self.max_intervals);
+        let start = (n - 1).checked_mul(n.checked_add(2)?)?.checked_div(2)?;
+        self.tables.get(start..start.checked_add(n.checked_add(1)?)?)
+    }
+
+    #[inline]
+    fn position<const SPIN: bool>(&self, dx: f32, dy: f32, (s, c): (f32, f32)) -> (f32, f32) {
+        if SPIN { (self.center.0 + dx * c - dy * s, self.center.1 + dx * s + dy * c) } else { (self.center.0 + dx * c, self.center.1 + dy * c) }
+    }
+}
+
+fn sample_len(rect: Rect, channels: usize) -> Option<usize> {
+    (rect.width() as usize).checked_mul(rect.height() as usize)?.checked_mul(channels)
+}
+
+/// Storage dimensions are checked once, before coordinate-derived offset arithmetic.
+struct Source<'a> {
+    image: &'a Image,
+    width: usize,
+}
+
+impl<'a> Source<'a> {
+    fn new(image: &'a Image) -> Option<Self> {
+        if image.ch == 0 || sample_len(image.rect, image.ch)? != image.data.len() {
+            return None;
+        }
+        Some(Self { image, width: image.rect.width() as usize })
+    }
+
+    #[inline]
+    fn pixel(&self, x: i64, y: i64) -> Option<&'a [f32]> {
+        let b = self.image.rect;
+        if x < i64::from(b.x0) || x >= i64::from(b.x1) || y < i64::from(b.y0) || y >= i64::from(b.y1) {
+            return None;
+        }
+        // The rectangular storage length was checked by Source::new, and the
+        // coordinates above are inside that rectangle, so these products fit.
+        let i = ((y - i64::from(b.y0)) as usize * self.width + (x - i64::from(b.x0)) as usize) * self.image.ch;
+        self.image.data.get(i..i + self.image.ch)
+    }
+}
+
+#[inline]
+fn taps(x: f32, y: f32) -> Option<[(i64, i64, f32); 4]> {
+    if !x.is_finite() || !y.is_finite() {
+        return None;
+    }
+    let (fx, fy) = (x - 0.5, y - 0.5);
+    let (x0, y0) = (fx.floor(), fy.floor());
+    let (ax, ay) = (fx - x0, fy - y0);
+    // f32 casts saturate, and i64 leaves room for the neighbouring tap even
+    // when an adversarial center places a sample outside the i32 image domain.
+    let (x0, y0) = (i64::from(x0 as i32), i64::from(y0 as i32));
+    Some([(x0, y0, (1.0 - ax) * (1.0 - ay)), (x0 + 1, y0, ax * (1.0 - ay)), (x0, y0 + 1, (1.0 - ax) * ay), (x0 + 1, y0 + 1, ax * ay)])
+}
+
+#[inline]
+fn sample_fixed<const N: usize, const ALPHA: bool>(src: &Source<'_>, x: f32, y: f32, acc: &mut [f32; N]) {
+    let Some(taps) = taps(x, y) else { return };
+    let mut sample = [0.0; N];
+    for (x, y, weight) in taps {
+        if weight <= 0.0 {
+            continue;
+        }
+        let Some(pixel) = src.pixel(x, y) else { continue };
+        let Ok(pixel) = <&[f32; N]>::try_from(pixel) else { continue };
+        if ALPHA {
+            let alpha = pixel.last().copied().unwrap_or(0.0);
+            for (sum, &value) in sample.iter_mut().zip(pixel).take(N.saturating_sub(1)) {
+                *sum += value * alpha * weight;
+            }
+            if let Some(sum) = sample.last_mut() {
+                *sum += alpha * weight;
+            }
+        } else {
+            for (sum, &value) in sample.iter_mut().zip(pixel) {
+                *sum += value * weight;
+            }
+        }
+    }
+    accumulate_sample(acc, &sample, ALPHA);
+}
+
+#[inline]
+fn sample_generic<const ALPHA: bool>(src: &Source<'_>, x: f32, y: f32, acc: &mut [f32], sample: &mut [f32]) {
+    let Some(taps) = taps(x, y) else { return };
+    sample.fill(0.0);
+    for (x, y, weight) in taps {
+        if weight <= 0.0 {
+            continue;
+        }
+        let Some(pixel) = src.pixel(x, y) else { continue };
+        if ALPHA {
+            let alpha = pixel.last().copied().unwrap_or(0.0);
+            let colours = sample.len().saturating_sub(1);
+            for (sum, &value) in sample.iter_mut().zip(pixel).take(colours) {
+                *sum += value * alpha * weight;
+            }
+            if let Some(sum) = sample.last_mut() {
+                *sum += alpha * weight;
+            }
+        } else {
+            for (sum, &value) in sample.iter_mut().zip(pixel) {
+                *sum += value * weight;
+            }
+        }
+    }
+    accumulate_sample(acc, sample, ALPHA);
+}
+
+#[inline]
+fn accumulate_sample(acc: &mut [f32], sample: &[f32], alpha: bool) {
+    // Retain Image::sample's unpremultiply followed by average_samples'
+    // premultiply, including its floating-point rounding and unusual alpha.
+    let a = sample.last().copied().unwrap_or(0.0);
+    let colours = sample.len().saturating_sub(1);
+    for (c, (sum, &value)) in acc.iter_mut().zip(sample).enumerate() {
+        *sum += if alpha && c < colours {
+            let straight = if a > 0.0 { value / a } else { value };
+            straight * a
+        } else {
+            value
+        };
+    }
+}
+
+#[inline]
+fn finish(pixel: &mut [f32], reciprocal: f32, alpha: bool) {
+    for value in pixel.iter_mut() {
+        *value *= reciprocal;
+    }
+    if alpha {
+        let a = pixel.last().copied().unwrap_or(0.0);
+        let colours = pixel.len().saturating_sub(1);
+        for value in pixel.iter_mut().take(colours) {
+            *value = if a > 1e-7 { *value / a } else { 0.0 };
+        }
+    }
+}
+
+/// Original implementation, retained only as a quality/performance oracle.
+#[cfg(test)]
+pub(crate) fn reference(src: &Image, out: Rect, ctx: &Ctx, amount: f32, method: RadialMethod, center: (f32, f32)) -> Vec<f32> {
+    super::radial(src, out, ctx, amount, method, RadialQuality::Draft, center)
+}
+
+#[cfg(test)]
+mod tests;

crates/algo/src/blur/radial/tests.rs

@@ -0,0 +1,260 @@
+use super::*;
+use photocraft_color::{ColorMode, PixelFormat, SampleType};
+use photocraft_raster::Surface;
+
+fn fixture(format: PixelFormat, rect: Rect, hdr: bool) -> Surface {
+    let mut surface = Surface::new(format);
+    for tile in rect.tiles() {
+        let area = tile.rect().intersect(&rect);
+        let mut values = Vec::with_capacity(area.width() as usize * area.height() as usize * format.channels());
+        for y in area.y0..area.y1 {
+            for x in area.x0..area.x1 {
+                let seed = (x as u32).wrapping_mul(0x9e37_79b9) ^ (y as u32).wrapping_mul(0x85eb_ca6b);
+                for channel in 0..format.channels() {
+                    let level = seed.rotate_left(channel as u32 * 7).wrapping_add(channel as u32 * 71) % 256;
+                    let alpha = format.alpha && channel + 1 == format.channels();
+                    let value = if alpha && seed.is_multiple_of(5) { 0.0 } else { level as f32 / 255.0 };
+                    values.push(if hdr && !alpha { value * 4.0 - 0.75 } else { value });
+                }
+            }
+        }
+        surface.write_region(area, &values);
+    }
+    surface
+}
+
+fn premultiplied_max_error(mut a: Vec<f32>, mut b: Vec<f32>, channels: usize, alpha: bool) -> f32 {
+    crate::image::premultiply(&mut a, channels, alpha);
+    crate::image::premultiply(&mut b, channels, alpha);
+    a.iter().zip(b).map(|(a, b)| (a - b).abs()).fold(0.0, f32::max)
+}
+
+#[test]
+fn direct_matches_legacy_depth_models_edges_centers_and_hdr() {
+    let source_rect = Rect::new(-17, -13, 45, 39);
+    let bounds = Rect::new(-41, -32, 51, 48);
+    let out = Rect::new(-22, -19, 47, 43);
+    let mut worst = 0.0f32;
+    for sample in [SampleType::U8, SampleType::U16, SampleType::F32] {
+        for (mode, alpha) in [
+            (ColorMode::Grayscale, false),
+            (ColorMode::Grayscale, true),
+            (ColorMode::Rgb, false),
+            (ColorMode::Rgb, true),
+            (ColorMode::Cmyk, false),
+            (ColorMode::Cmyk, true),
+            (ColorMode::Lab, true),
+        ] {
+            let format = PixelFormat::new(mode, sample, alpha);
+            let src = Image::read(&fixture(format, source_rect, sample == SampleType::F32), source_rect);
+            let ctx = Ctx { bounds, mode, alpha };
+            for method in [RadialMethod::Spin, RadialMethod::Zoom] {
+                for (amount, center) in [(0.0, (0.5, 0.5)), (1.0, (0.5, 0.5)), (20.0, (0.53, -0.17)), (100.0, (0.0, 1.0))] {
+                    let direct = Plan::new(bounds, amount, method, center).filter(&src, out, &ctx, &Interrupt::NONE).expect("valid input");
+                    let legacy = reference(&src, out, &ctx, amount, method, center);
+                    assert_eq!(direct, legacy, "preserve finite upstream samples exactly");
+                    let error = premultiplied_max_error(direct, legacy, src.ch, alpha);
+                    worst = worst.max(error);
+                    assert!(error <= 2e-5, "{format:?}, {method:?}, amount {amount}, center {center:?}: {error}");
+                }
+            }
+        }
+    }
+    eprintln!("direct vs legacy: maximum premultiplied float error {worst}");
+}
+
+#[test]
+fn direct_generic_channels_and_unusual_float_alpha_match_legacy() {
+    let rect = Rect::new(-4, -3, 21, 16);
+    for channels in [1, 2, 4, 5, 7, 12] {
+        let mut src = Image::new(rect, channels);
+        for (index, pixel) in src.data.chunks_exact_mut(channels).enumerate() {
+            for (channel, value) in pixel.iter_mut().enumerate() {
+                *value = ((index * 13 + channel * 3) % 37) as f32 / 9.0 - 0.75;
+            }
+            *pixel.last_mut().expect("nonempty pixel") = [-0.7, 0.0, 1e-9, 0.4, 1.7, f32::NAN, f32::INFINITY, f32::NEG_INFINITY][index % 8];
+        }
+        for alpha in [false, true] {
+            let ctx = Ctx { bounds: rect, mode: ColorMode::Rgb, alpha };
+            for method in [RadialMethod::Spin, RadialMethod::Zoom] {
+                let fast = Plan::new(rect, 53.0, method, (0.51, 0.47)).filter(&src, rect, &ctx, &Interrupt::NONE).expect("valid input");
+                let slow = reference(&src, rect, &ctx, 53.0, method, (0.51, 0.47));
+                for (a, b) in fast.iter().zip(&slow) {
+                    if b.is_nan() {
+                        assert!(a.is_nan(), "match legacy NaN classification");
+                    } else if b.is_infinite() {
+                        assert_eq!(a, b, "match legacy infinity sign");
+                    } else {
+                        assert_eq!(a, b, "finite upstream output is preserved exactly");
+                    }
+                }
+                let error = fast.iter().zip(&slow).map(|(a, b)| (a - b).abs() / b.abs().max(1.0)).fold(0.0f32, f32::max);
+                assert!(error < 2e-5, "channels {channels}, alpha {alpha}, {method:?}: relative error {error}");
+            }
+        }
+    }
+}
+
+#[test]
+fn direct_invalid_storage_extreme_coordinates_and_cancellation_are_bounded() {
+    let rect = Rect::new(0, 0, 23, 19);
+    let ctx = Ctx { bounds: rect, mode: ColorMode::Rgb, alpha: true };
+    let plan = Plan::new(rect, 20.0, RadialMethod::Spin, (0.5, 0.5));
+    let mut src = Image::new(rect, 4);
+    let cancel = || true;
+    assert!(plan.filter(&src, rect, &ctx, &Interrupt::cancel_only(&cancel)).is_none());
+    src.data.pop();
+    assert!(plan.filter(&src, rect, &ctx, &Interrupt::NONE).is_none());
+    src.ch = 0;
+    assert!(plan.filter(&src, rect, &ctx, &Interrupt::NONE).is_none());
+    let src = Image::new(Rect::EMPTY, 4);
+    assert!(plan.filter(&src, Rect::new(i32::MIN, i32::MIN, i32::MAX, i32::MAX), &ctx, &Interrupt::NONE).is_none());
+    assert!(Plan::new(rect, f32::NAN, RadialMethod::Spin, (0.5, 0.5)).filter(&src, rect, &ctx, &Interrupt::NONE).is_none());
+    assert!(Plan::new(rect, 20.0, RadialMethod::Spin, (f32::INFINITY, 0.5)).filter(&src, rect, &ctx, &Interrupt::NONE).is_none());
+    assert_eq!(plan.filter(&src, Rect::EMPTY, &ctx, &Interrupt::NONE), Some(Vec::new()));
+    let edge = Rect::new(i32::MAX - 2, i32::MAX - 2, i32::MAX, i32::MAX);
+    let src = Image::new(edge, 4);
+    assert!(Plan::new(edge, 100.0, RadialMethod::Spin, (-100.0, 100.0)).filter(&src, edge, &ctx, &Interrupt::NONE).is_some());
+}
+
+#[test]
+fn direct_observes_cancellation_inside_an_output_tile() {
+    use std::sync::atomic::{AtomicUsize, Ordering};
+    let rect = Rect::new(0, 0, 16, 24);
+    let src = Image::new(rect, 4);
+    let ctx = Ctx { bounds: rect, mode: ColorMode::Rgb, alpha: true };
+    let calls = AtomicUsize::new(0);
+    let cancel = || calls.fetch_add(1, Ordering::Relaxed) >= 2;
+    assert!(Plan::new(rect, 100.0, RadialMethod::Spin, (0.5, 0.5)).filter(&src, rect, &ctx, &Interrupt::cancel_only(&cancel)).is_none());
+    assert_eq!(calls.load(Ordering::Relaxed), 3, "entry, first row, eighth row");
+}
+
+#[test]
+fn direct_apply_preserves_selection_and_canvas_edges_at_all_depths() {
+    let rect = Rect::new(0, 0, 37, 29);
+    let bounds = Rect::new(3, 5, 31, 23);
+    let selection = Surface::with_default(PixelFormat::new(ColorMode::Grayscale, SampleType::F32, false), &[0.37]);
+    for sample in [SampleType::U8, SampleType::U16, SampleType::F32] {
+        let format = PixelFormat::new(ColorMode::Rgb, sample, true);
+        let src = fixture(format, rect, sample == SampleType::F32);
+        for method in [RadialMethod::Spin, RadialMethod::Zoom] {
+            let params = crate::FilterParams::RadialBlur { quality: crate::RadialQuality::Draft, amount: 100.0, method, center_x: 0.2, center_y: 0.8 };
+            let fast = crate::apply_in_with(&src, &params, rect, bounds, Some(&selection), rect, &Interrupt::NONE).expect("direct input");
+            let legacy = crate::apply_radial_reference(&src, &params, rect, bounds, Some(&selection), rect, &Interrupt::NONE).expect("valid input");
+            assert_eq!(fast.read_region(rect), legacy.read_region(rect));
+            let error = premultiplied_max_error(fast.read_region(rect), legacy.read_region(rect), 4, true);
+            let tolerance = match sample {
+                SampleType::U8 => 1.0 / 255.0,
+                SampleType::U16 => 1.0 / 65535.0,
+                _ => 0.0,
+            } + 2e-5;
+            assert!(error <= tolerance, "{sample:?}, {method:?}: {error}");
+        }
+    }
+}
+
+#[test]
+fn direct_preserves_upstream_quality_ceilings_and_discrete_samples() {
+    let bounds = Rect::new(0, 0, 1000, 1000);
+    let out = Rect::new(999, 500, 1000, 501);
+    let mut source = Image::new(bounds, 1);
+    for (i, v) in source.data.iter_mut().enumerate() {
+        *v = ((i * 37 + i / 1000 * 13) % 101) as f32 / 100.0;
+    }
+    let ctx = Ctx { bounds, mode: ColorMode::Grayscale, alpha: false };
+    for (quality, limit) in [(crate::RadialQuality::Draft, 64), (crate::RadialQuality::Good, 256), (crate::RadialQuality::Best, 4096)] {
+        for method in [RadialMethod::Spin, RadialMethod::Zoom] {
+            let plan = Plan::with_quality(bounds, 100.0, method, quality, (0.5, 0.5), &Interrupt::NONE).expect("quality plan");
+            assert_eq!(plan.max_intervals, limit);
+            assert_eq!(plan.table(10000.0, 0.0).expect("capped table").len(), limit + 1);
+            let actual = plan.filter(&source, out, &ctx, &Interrupt::NONE).expect("direct output");
+            let expected = super::super::average_samples(&source, out, &ctx, |x, y, samples| {
+                let (dx, dy) = (x - 500.0, y - 500.0);
+                let extent = match method {
+                    RadialMethod::Spin => 100.0f32.to_radians(),
+                    RadialMethod::Zoom => 0.5,
+                };
+                let intervals = (extent * (dx * dx + dy * dy).sqrt()).ceil().clamp(1.0, limit as f32) as usize;
+                for i in 0..=intervals {
+                    samples.push(match method {
+                        RadialMethod::Spin => {
+                            let (sin, cos) = ((i as f32 / intervals as f32 - 0.5) * extent).sin_cos();
+                            (500.0 + dx * cos - dy * sin, 500.0 + dx * sin + dy * cos)
+                        }
+                        RadialMethod::Zoom => {
+                            let k = 1.0 - extent * i as f32 / intervals as f32;
+                            (500.0 + dx * k, 500.0 + dy * k)
+                        }
+                    });
+                }
+            });
+            assert_eq!(actual, expected, "{quality:?}, {method:?}: preserve discrete positions");
+        }
+    }
+}
+
+#[test]
+fn quality_plan_cancellation_stops_table_preparation() {
+    use std::sync::atomic::{AtomicUsize, Ordering};
+    let calls = AtomicUsize::new(0);
+    let cancel = || calls.fetch_add(1, Ordering::Relaxed) >= 2;
+    assert!(
+        Plan::with_quality(Rect::new(0, 0, 100, 100), 100.0, RadialMethod::Spin, crate::RadialQuality::Best, (0.5, 0.5), &Interrupt::cancel_only(&cancel))
+            .is_none()
+    );
+    assert_eq!(calls.load(Ordering::Relaxed), 3);
+}
+
+#[test]
+fn all_qualities_preserve_samples_at_large_radii_across_models_and_depths() {
+    let bounds = Rect::new(0, 0, 6000, 4000);
+    let source_rect = Rect::new(5970, 1980, 6000, 2020);
+    let out = Rect::new(5990, 2000, 5992, 2002);
+    for quality in [crate::RadialQuality::Draft, crate::RadialQuality::Good, crate::RadialQuality::Best] {
+        for method in [RadialMethod::Spin, RadialMethod::Zoom] {
+            let plan = Plan::with_quality(bounds, 100.0, method, quality, (0.5, 0.5), &Interrupt::NONE).unwrap();
+            for sample in [SampleType::U8, SampleType::U16, SampleType::F32] {
+                for mode in [ColorMode::Grayscale, ColorMode::Rgb, ColorMode::Cmyk, ColorMode::Lab] {
+                    for alpha in [false, true] {
+                        let format = PixelFormat::new(mode, sample, alpha);
+                        let src = Image::read(&fixture(format, source_rect, sample == SampleType::F32), source_rect);
+                        let ctx = Ctx { bounds, mode, alpha };
+                        let actual = plan.filter(&src, out, &ctx, &Interrupt::NONE).unwrap();
+                        let expected = super::super::radial(&src, out, &ctx, 100.0, method, quality, (0.5, 0.5));
+                        assert_eq!(actual, expected, "{quality:?} {method:?} {format:?}");
+                    }
+                }
+            }
+        }
+    }
+}
+
+#[test]
+fn unavailable_optional_plan_uses_upstream_instead_of_cancellation() {
+    let rect = Rect::new(0, 0, 7, 5);
+    let source = fixture(PixelFormat::RGBA8, rect, false);
+    // Nonfinite centers decline preparation. The unchanged upstream sampler is
+    // still used; a declined optimization must not look like engine cancellation.
+    let params =
+        crate::FilterParams::RadialBlur { amount: 20.0, method: RadialMethod::Spin, quality: crate::RadialQuality::Good, center_x: f32::NAN, center_y: 0.5 };
+    let actual = crate::apply_in_with(&source, &params, rect, rect, None, rect, &Interrupt::NONE).unwrap();
+    let expected = crate::apply_radial_reference(&source, &params, rect, rect, None, rect, &Interrupt::NONE).unwrap();
+    assert_eq!(actual.read_region(rect), expected.read_region(rect));
+}
+
+#[test]
+fn small_best_jobs_prepare_only_reachable_tables() {
+    let rect = Rect::new(-30, -20, 30, 20);
+    let format = PixelFormat::RGBA32F;
+    let source = fixture(format, rect, true);
+    let ctx = Ctx { bounds: rect, mode: format.mode, alpha: true };
+    let image = Image::read(&source, rect);
+    for method in [RadialMethod::Spin, RadialMethod::Zoom] {
+        let plan = Plan::for_area(rect, 100.0, method, crate::RadialQuality::Best, (0.5, 0.5), rect, &Interrupt::NONE).unwrap();
+        assert!(plan.max_intervals < 64, "small Best image does not need 4096 intervals");
+        assert!(plan.tables.len() < 2200);
+        let actual = plan.filter(&image, rect, &ctx, &Interrupt::NONE).unwrap();
+        let expected = super::super::radial(&image, rect, &ctx, 100.0, method, crate::RadialQuality::Best, (0.5, 0.5));
+        assert_eq!(actual, expected);
+    }
+}

crates/algo/src/lib.rs

@@ -983,6 +983,24 @@ pub fn apply_tiled_with(
     extent: Option<Rect>,
     ctl: &photocraft_raster::Interrupt,
 ) -> Option<Surface> {
+    apply_tiled_with_impl(surface, params, area, bounds, selection, tile, extent, ctl, true)
+}
+
+#[allow(clippy::too_many_arguments)]
+fn apply_tiled_with_impl(
+    surface: &Surface,
+    params: &FilterParams,
+    area: Rect,
+    bounds: Rect,
+    selection: Option<&Surface>,
+    tile: i32,
+    extent: Option<Rect>,
+    ctl: &photocraft_raster::Interrupt,
+    prepare_radial: bool,
+) -> Option<Surface> {
+    if ctl.cancelled() {
+        return None;
+    }
     let mut out = surface.clone();
     if area.is_empty() {
         return Some(out);
@@ -998,6 +1016,18 @@ pub fn apply_tiled_with(
     let tile = if matches!(params, FilterParams::ShapeBlur { .. }) { tile.clamp(1, 2048) } else { tile };
     let fmt = surface.format();
     let ctx = Ctx { bounds, mode: fmt.mode, alpha: fmt.alpha };
+    // Prepare the immutable sample transforms once for all output tiles. If the
+    // optional table allocation fails, use the existing sampler; only a real
+    // cancellation returns None. Source reads keep their upstream semantics.
+    let radial = match params {
+        FilterParams::RadialBlur { amount, method, quality, center_x, center_y } if prepare_radial => {
+            blur::radial::Plan::for_area(bounds, *amount, *method, *quality, (*center_x, *center_y), area, ctl)
+        }
+        _ => None,
+    };
+    if ctl.cancelled() {
+        return None;
+    }
     let halo = params.halo_for(bounds);
     let shared = (halo == Halo::Bounds).then(|| Image::read(surface, bounds.union(&area)));
     // Shape rasterization depends only on params, not on a tile's pixels.
@@ -1049,9 +1079,14 @@ pub fn apply_tiled_with(
                 &owned
             }
         };
-        let Some(mut data) = (match &shape {
-            Some(shape) => shape.checked_run(src, *t, &ctx),
-            None => Some(kernel(params, src, *t, &ctx)),
+        let Some(mut data) = (match (&shape, &radial) {
+            (Some(shape), _) => shape.checked_run(src, *t, &ctx),
+            (None, Some(plan)) => match plan.filter(src, *t, &ctx, ctl) {
+                Some(data) => Some(data),
+                None if ctl.cancelled() => None,
+                None => Some(kernel(params, src, *t, &ctx)),
+            },
+            (None, None) => Some(kernel(params, src, *t, &ctx)),
         }) else {
             return (*t, Vec::new());
         };
@@ -1189,3 +1224,18 @@ mod tests_ext;
 
 #[cfg(test)]
 mod tests_mosaic;
+
+/// The current upstream sampler through the identical source/selection/tile pipeline.
+#[cfg(test)]
+#[allow(clippy::too_many_arguments)]
+fn apply_radial_reference(
+    surface: &Surface,
+    params: &FilterParams,
+    area: Rect,
+    bounds: Rect,
+    selection: Option<&Surface>,
+    extent: Rect,
+    ctl: &photocraft_raster::Interrupt,
+) -> Option<Surface> {
+    apply_tiled_with_impl(surface, params, area.intersect(&extent), bounds, selection, auto_tile(params, bounds), Some(extent), ctl, false)
+}

docs/radial-blur-performance.md

@@ -0,0 +1,56 @@
+# Radial Blur sampler optimization
+
+The tiled filter prepares Spin rotations or Zoom scales once per operation and
+reuses them across output tiles. Bilinear taps load all channels together, with
+specialized loops for one through five channels and a generic fallback. Rayon
+continues to parallelize output tiles on native builds.
+
+Draft, Good and Best keep upstream's 64, 256 and 4096 interval ceilings, sample
+positions, transparent edges and alpha arithmetic. The positive-alpha
+unpremultiply/premultiply sequence is retained to preserve floating-point rounding.
+Tests compare against the unchanged upstream sampler at large radii across
+Gray, RGB, CMYK and Lab, with and without alpha, at U8, U16 and F32 depths.
+Selection blending and unusual float alpha are covered as well.
+
+Only interval counts reachable within the output area are prepared; small Best
+images avoid allocating/preparing the full roughly 64 MiB transform table.
+Preparation and output allocations use checked sizes and fallible reservation.
+If the optional optimization cannot prepare its tables, the existing sampler
+runs instead. Cancellation is checked during preparation and within output tiles.
+
+This revision adds no source-size rejection: it retains upstream's shared
+normalized source read and therefore still requires memory proportional to source
+size. It does not claim bounded source memory or unlimited physical capacity.
+The polar approximation and background preview scheduler are separate work.
+
+## Local comparison (2026-10-09)
+
+| Method | Upstream | Optimized | Speedup | Output |
+|---|---:|---:|---:|---|
+| Spin | 66.10 s | 16.43 s | 4.02× | Exact match |
+| Zoom | 19.19 s | 5.35 s | 3.59× | Exact match |
+
+
+Baseline: current upstream sampler at main `578bc905`, through `kernel`, using
+the same shared source read, 256px output tiles, 256 MiB result batches, surface
+writeback and pruning as the tiled filter. Optimized: `apply_in`.
+The benchmark's baseline wrapper covers full-canvas, unselected filtering only;
+selection behavior is checked separately in tests.
+
+Machine: Intel Core i7-9750H (6 cores / 12 logical workers), Linux, Rust 1.98.1,
+release profile with opt-level 3, thin LTO and one codegen unit. Source: deterministic
+textured RGBA8 with transparent pixels, 6000×4000 (24 MP). Good quality, amount 1,
+center (0.5, 0.5), 12 Rayon workers. One paired run per method, upstream first,
+after a small untimed warmup. Entire filtering calls are timed; fixture generation
+and exact comparison of every output tile are outside the timer.
+
+These are single local measurements, not a cross-machine baseline or an
+amount-100 performance claim. The lower amount keeps this review benchmark
+bounded on the local laptop. Quality ceilings are verified separately at large
+radii. The example defaults to three alternating-order runs and amount 100 for
+longer measurements. Historical handoff polar/default timings do not describe
+this output-preserving implementation.
+
+```sh
+cargo run --release -j 12 -p photocraft-algo --example bench_radial -- 6000 4000 1 1
+```

docs/roadmap.md

@@ -58,6 +58,11 @@ almost entirely missing or partial.
 of display and canvas zoom (#532). The current preference audit drops from 59 to 58 unread
 settings out of 135; see the regenerated scorecard.
 
+2026-10-09: Output-preserving Radial Blur direct sampling measured on a local
+i7-9750H, 12 workers, 24 MP RGBA8, Good quality, amount 1: Spin 66.10 → 16.43 s
+(4.02×), Zoom 19.19 → 5.35 s (3.59×), one paired run each, exact output equality.
+These are local measurements at amount 1; see [method and limits](radial-blur-performance.md).
+
 **Bottom line.** Two days after 0.2.0 we had merged ~96 PRs and closed ~48 issues, but **real
 Photoshop parity is still well below 50%**. The biggest gaps are AI, missing tools, professional
 workflow depth and the plug-in ecosystem. Most fixes since 0.2.0 have passed our tests but have