//! Colour transforms between profiles through the PCS, with rendering intents, black point
//! compensation, soft-proofing chains, and fast lookup-table evaluation of pixel buffers.

use std::collections::HashMap;
use std::sync::{Arc, Mutex, OnceLock, PoisonError};

use crate::clut::{Clut, MAX_CHANNELS};
use crate::curve::Curve;
use crate::math;
use crate::pipeline::{Pipeline, Stage};
use crate::profile::{ColorSpace, Pcs, Profile, ProfileClass};
use crate::{CmsError, GpuMatrixShaper, Intent, check_len, invalid_arg};

/// Sample storage of a byte buffer for [`Transform::convert_bytes_many`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum SampleKind {
    U8,
    U16,
    F32,
}

/// Options for building a [`Transform`].
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub struct TransformOptions {
    pub intent: Intent,
    /// Adobe-style black point compensation (ignored for absolute colorimetric).
    pub bpc: bool,
    /// `true`: float buffers go through the exact pipeline (default). `false`: float buffers use
    /// the same lookup tables as integer buffers (faster, for previews).
    pub precise_float: bool,
}

impl Default for TransformOptions {
    fn default() -> Self {
        Self { intent: Intent::RelativeColorimetric, bpc: true, precise_float: true }
    }
}

/// Grid sizes of the device-link tables built for integer (and fast float) evaluation.
pub const GRID_3D: usize = 33;
pub const GRID_4D: usize = 17;
const TABLE_1D: usize = GpuMatrixShaper::INPUT_INTERVALS;
/// Resolution of output shaper tables (indexed by `x^(1/4)`).
const OUT_TABLE: usize = GpuMatrixShaper::OUTPUT_ENTRIES;
// The blocked matrix/TRC path wraps table indices with a mask instead of a bounds check.
const _: () = assert!(OUT_TABLE.is_power_of_two());

/// ICC v4 perceptual reference medium black (ICC.1:2010 §6.3.4.4), XYZ.
const PERCEPTUAL_BLACK: [f64; 3] = [0.00336, 0.0034731, 0.00287];

/// Accelerated representation of a transform.
#[derive(Debug)]
enum Core {
    /// Input curves (sampled per input channel over `[0, 1]`, 4096+1 entries; 256 exact
    /// entries for 8-bit input) followed by an optional matrix.
    Shaper { curves8: Vec<[f32; 256]>, curves: Vec<Vec<f32>>, matrix: Option<(usize, usize, Vec<f32>, Vec<f32>)> },
    /// Device link table over the input space.
    Grid(Clut),
    /// No acceleration possible: use the exact pipeline.
    Exact,
}

/// Output inverse-curve tables indexed by `x^(1/4)` (accurate near black even for pure gamma
/// curves, whose slope is infinite at zero).
#[derive(Debug)]
struct Post {
    tables: Vec<Vec<f32>>,
}

impl Post {
    #[inline]
    fn eval(&self, c: usize, x: f32) -> f32 {
        let t = &self.tables[c];
        let u = x.clamp(0.0, 1.0).sqrt().sqrt() * (OUT_TABLE - 1) as f32;
        let i = (u as usize).min(OUT_TABLE - 2);
        let f = u - i as f32;
        t[i] + (t[i + 1] - t[i]) * f
    }
}

/// A colour transform from one profile's device space to another's.
///
/// Build once (it precomputes lookup tables) and apply to any number of buffers from any
/// thread. Buffers are interleaved; colour channels come first in each pixel and any further
/// channels (alpha, spot) are copied or left alone.
#[derive(Debug)]
pub struct Transform {
    inputs: usize,
    outputs: usize,
    opts: TransformOptions,
    pipeline: Pipeline,
    core: Core,
    post: Option<Post>,
    /// Pre-located grid coordinates per 8-bit code: (offset contribution, fraction) per input.
    grid8: Vec<[(u32, f32); 256]>,
    /// For an exact core: the first stage's curves at every 8-bit code.
    exact8: Vec<[f32; 256]>,
    white: Option<White>,
}

/// Device white in and out, for transforms that map white to white exactly (every intent but
/// absolute colorimetric).
#[derive(Debug)]
struct White {
    src: Vec<f32>,
    dst: Vec<f32>,
    /// 8-bit codes of `src` when the 8-bit grid path must check for it (white is not a node of
    /// the grid, e.g. Lab's a* = b* = 0).
    src8: Option<Vec<u8>>,
}

impl Transform {
    /// Transform from `src` to `dst` with `intent`, optionally with black point compensation.
    pub fn new(src: &Profile, dst: &Profile, intent: Intent, bpc: bool) -> Result<Transform, CmsError> {
        Self::with_options(src, dst, TransformOptions { intent, bpc, ..Default::default() })
    }

    pub fn with_options(src: &Profile, dst: &Profile, opts: TransformOptions) -> Result<Transform, CmsError> {
        let stages = link_stages(src, dst, opts.intent, opts.bpc)?;
        let mut t = Self::from_pipeline(Pipeline::new(src.channels(), stages)?, opts)?;
        if opts.intent != Intent::AbsoluteColorimetric {
            t.keep_white(src.color_space, dst.color_space);
        }
        Ok(t)
    }

    /// Soft-proof transform `src → proof → dst` (Photoshop View › Proof Setup): `intent`/`bpc`
    /// apply to the document → proof step; the proof → display step is relative colorimetric
    /// with BPC, or absolute colorimetric without BPC when `simulate_paper` is set (paper
    /// colour and black ink simulated). Without paper simulation, and unless `intent` is
    /// absolute colorimetric, device white maps to device white.
    pub fn proof(src: &Profile, proof: &Profile, dst: &Profile, intent: Intent, bpc: bool, simulate_paper: bool) -> Result<Transform, CmsError> {
        let mut stages = link_stages(src, proof, intent, bpc)?;
        stages.push(Stage::Clamp01);
        let (i2, b2) = if simulate_paper { (Intent::AbsoluteColorimetric, false) } else { (Intent::RelativeColorimetric, true) };
        stages.extend(link_stages(proof, dst, i2, b2)?);
        let mut t = Self::from_pipeline(Pipeline::new(src.channels(), stages)?, TransformOptions { intent, bpc, precise_float: true })?;
        if !simulate_paper && intent != Intent::AbsoluteColorimetric {
            t.keep_white(src.color_space, dst.color_space);
        }
        Ok(t)
    }

    /// Builds a transform from an arbitrary pipeline over normalised values.
    pub fn from_pipeline(mut pipeline: Pipeline, opts: TransformOptions) -> Result<Transform, CmsError> {
        pipeline.optimize();
        let (inputs, outputs) = (pipeline.inputs(), pipeline.outputs());
        // Split trailing inverse curves into accurate output tables.
        let mut head = pipeline.stages().to_vec();
        let mut post = None;
        if let Some(Stage::InvCurves(cs)) = head.last()
            && cs.len() == outputs
        {
            let tables = cs
                .iter()
                .map(|c| {
                    let inverse = c.inverse64_many();
                    (0..OUT_TABLE)
                        .map(|i| {
                            let u = i as f64 / (OUT_TABLE - 1) as f64;
                            inverse(u * u * u * u).clamp(0.0, 1.0) as f32
                        })
                        .collect()
                })
                .collect();
            post = Some(Post { tables });
            head.pop();
            if matches!(head.last(), Some(Stage::Clamp01)) {
                head.pop();
            }
        }
        let head_pipe = Pipeline::new(inputs, head.clone())?;
        let mid = head_pipe.outputs();
        let sampled = |grid: Vec<usize>| match Clut::sample(grid, mid, |i, o| head_pipe.run(i, o)) {
            Ok(c) => Core::Grid(c),
            // More outputs than a table holds: evaluate exactly.
            Err(_) => Core::Exact,
        };
        let core = match head.as_slice() {
            [] => Core::Shaper { curves8: Vec::new(), curves: Vec::new(), matrix: None },
            [Stage::Curves(cs)] if cs.len() == inputs => shaper(cs, None),
            [Stage::Matrix { rows, cols, m, offset }] => shaper(&[], Some((*rows, *cols, m, offset))),
            [Stage::Curves(cs), Stage::Matrix { rows, cols, m, offset }] if cs.len() == inputs => shaper(cs, Some((*rows, *cols, m, offset))),
            // Purely analytic heads (Lab/XYZ maths, matrices, curves) are evaluated exactly:
            // a grid would amplify interpolation error through steep output curves.
            _ if !head.iter().any(|s| matches!(s, Stage::Clut(_))) => Core::Exact,
            _ => match inputs {
                1 => sampled(vec![TABLE_1D]),
                3 => sampled(vec![GRID_3D; 3]),
                4 => sampled(vec![GRID_4D; 4]),
                2 => sampled(vec![GRID_3D; 2]),
                _ => Core::Exact,
            },
        };
        let grid8 = match &core {
            Core::Grid(c) if c.inputs >= 2 => (0..c.inputs)
                .map(|d| {
                    let n = c.grid[d];
                    let stride = c.data.len() / c.outputs / c.grid[..=d].iter().product::<usize>() * c.outputs;
                    std::array::from_fn(|v| {
                        let p = v as f32 / 255.0 * (n - 1) as f32;
                        let i = (p as usize).min(n - 2);
                        ((i * stride) as u32, p - i as f32)
                    })
                })
                .collect(),
            _ => Vec::new(),
        };
        let exact8 = match (&core, pipeline.stages().first()) {
            (Core::Exact, Some(Stage::Curves(cs))) if cs.len() == inputs && inputs <= 4 => {
                cs.iter().map(|c| std::array::from_fn(|i| c.eval(i as f32 / 255.0))).collect()
            }
            _ => Vec::new(),
        };
        Ok(Transform { inputs, outputs, opts, pipeline, core, post, grid8, exact8, white: None })
    }

    /// White preservation: the source's device white maps to exactly the destination's device
    /// white. LUT profiles rarely put white on a grid node (the legacy Lab encoding puts L* 100
    /// at 0.996 and a* = 0 at 0.50001), so interpolation would give a trace of ink or an
    /// additive white slightly off 1. Every other colour is left as the profiles give it.
    ///
    /// The device-link node at white is set to the destination white (`fix_white`),
    /// and the per-pixel paths that do not read that node exactly (the exact pipeline, float
    /// interpolation, 8-bit lookups when white is no node) compare each pixel with white.
    fn keep_white(&mut self, src: ColorSpace, dst: ColorSpace) {
        let (Some(sw), Some(dw)) = (device_white(src), device_white(dst)) else { return };
        if sw.len() != self.inputs || dw.len() != self.outputs {
            return;
        }
        let on_node = self.fix_white(&sw, &dw);
        let src8 = (!on_node && matches!(self.core, Core::Grid(_))).then(|| sw.iter().map(|v| (v * 255.0).round() as u8).collect());
        self.white = Some(White { src: sw, dst: dw, src8 });
    }

    /// Sets the device-link node at `sw` to `dw`; `false` when `sw` is no grid node.
    fn fix_white(&mut self, sw: &[f32], dw: &[f32]) -> bool {
        if self.post.is_some() {
            // Output curves map 1 → 1 for additive spaces; the grid holds pre-curve values.
            if dw.iter().any(|v| *v != 1.0) {
                return false;
            }
        }
        let Core::Grid(c) = &mut self.core else { return false };
        if c.inputs != sw.len() || c.outputs != dw.len() {
            return false;
        }
        let mut off = 0usize;
        let mut stride = c.outputs;
        for d in (0..c.inputs).rev() {
            let n = c.grid[d];
            let pos = sw[d] * (n - 1) as f32;
            if (pos - pos.round()).abs() > 1e-6 {
                return false;
            }
            off += pos.round() as usize * stride;
            stride *= n;
        }
        c.data[off..off + c.outputs].copy_from_slice(dw);
        true
    }

    /// The destination white when `input` is exactly the source's device white.
    #[inline]
    fn white_for(&self, input: &[f32]) -> Option<&[f32]> {
        self.white.as_ref().filter(|w| input[..self.inputs] == w.src[..]).map(|w| &w.dst[..])
    }

    /// The exact pipeline, with device white kept exact.
    #[inline]
    fn exact(&self, input: &[f32], out: &mut [f32]) {
        match self.white_for(input) {
            Some(w) => out[..self.outputs].copy_from_slice(w),
            None => self.pipeline.run(input, out),
        }
    }

    pub fn inputs(&self) -> usize {
        self.inputs
    }
    pub fn outputs(&self) -> usize {
        self.outputs
    }
    pub fn options(&self) -> TransformOptions {
        self.opts
    }
    /// The exact evaluation pipeline.
    pub fn pipeline(&self) -> &Pipeline {
        &self.pipeline
    }

    /// A separable RGB transform for accurate GPU display conversion, when available.
    ///
    /// Reuses the input/matrix/output acceleration of the integer conversion path instead
    /// of interpolating final encoded RGB over a coarse 3D grid. Returns `None` for other
    /// channel counts, CLUTs, proofing stages, or invalid/unavailable table data, so callers
    /// can keep their existing display-LUT fallback. This is an SDR representation; see
    /// [`GpuMatrixShaper`] for the precise sampling and clamping contract.
    pub fn gpu_matrix_shaper(&self) -> Option<GpuMatrixShaper> {
        if self.inputs != 3 || self.outputs != 3 {
            return None;
        }
        let Core::Shaper { curves, matrix, .. } = &self.core else { return None };
        let mut rows = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]];
        if let Some((nr, nc, coefficients, offsets)) = matrix {
            if *nr != 3 || *nc != 3 || coefficients.len() != 9 || offsets.len() != 3 {
                return None;
            }
            for (row, (coefficients, offset)) in rows.iter_mut().zip(coefficients.as_chunks::<3>().0.iter().zip(offsets)) {
                let [r, g, b] = coefficients;
                *row = [*r, *g, *b, *offset];
            }
        }
        GpuMatrixShaper::from_parts(rows, curves, self.post.as_ref().map(|p| p.tables.as_slice()))
    }

    /// Exact evaluation of one colour: `input` holds at least `inputs` normalised values, `out`
    /// receives `outputs` values.
    pub fn eval(&self, input: &[f32], out: &mut [f32]) -> Result<(), CmsError> {
        self.check_io(input, out)?;
        self.exact(input, out);
        Ok(())
    }

    /// Lookup-table evaluation of one colour (what integer buffers use); slices as in
    /// [`Transform::eval`].
    pub fn eval_fast(&self, input: &[f32], out: &mut [f32]) -> Result<(), CmsError> {
        self.check_io(input, out)?;
        self.fast(input, out);
        Ok(())
    }

    fn check_io(&self, input: &[f32], out: &[f32]) -> Result<(), CmsError> {
        check_len("transform input", input.len(), self.inputs)?;
        check_len("transform output", out.len(), self.outputs)
    }

    #[inline]
    fn fast(&self, input: &[f32], out: &mut [f32]) {
        let mut mid = [0.0f32; 16];
        match &self.core {
            Core::Exact => return self.exact(input, out),
            Core::Grid(c) => {
                if let Some(w) = self.white_for(input) {
                    return out[..self.outputs].copy_from_slice(w);
                }
                let mut inp = [0.0f32; 16];
                for (d, s) in inp.iter_mut().zip(&input[..self.inputs]) {
                    *d = s.clamp(0.0, 1.0);
                }
                c.interp(&inp[..self.inputs], &mut mid[..c.outputs]);
            }
            Core::Shaper { curves, matrix, .. } => {
                let mut v = [0.0f32; 16];
                for k in 0..self.inputs {
                    v[k] = if curves.is_empty() { input[k] } else { table(&curves[k], input[k]) };
                }
                apply_matrix(matrix, &v, &mut mid, self.inputs);
            }
        }
        self.finish(&mid, out);
    }

    #[inline]
    fn finish(&self, mid: &[f32], out: &mut [f32]) {
        match &self.post {
            Some(p) => {
                for k in 0..self.outputs {
                    out[k] = p.eval(k, mid[k]);
                }
            }
            None => out[..self.outputs].copy_from_slice(&mid[..self.outputs]),
        }
    }

    /// Converts float pixels. `src`/`dst` hold whole pixels of `src_stride`/`dst_stride`
    /// values (strides at least `inputs`/`outputs`); colour channels come first. `dst` holds at
    /// least as many pixels as `src`; pixels past that are left alone. Trailing channels are
    /// copied when `copy_extra`. Buffers that do not fit are refused with
    /// [`CmsError::InvalidArgument`] before any pixel is written.
    pub fn convert_f32(&self, src: &[f32], src_stride: usize, dst: &mut [f32], dst_stride: usize, copy_extra: bool) -> Result<(), CmsError> {
        self.check_buffers(src.len(), src_stride, dst.len(), dst_stride)?;
        let extra = extra_channels(self, src_stride, dst_stride, copy_extra);
        par_rows(src, src_stride, dst, dst_stride, |s, d| self.run_f32(s, src_stride, d, dst_stride, extra));
        Ok(())
    }

    fn run_f32(&self, s: &[f32], ss: usize, d: &mut [f32], ds: usize, extra: usize) {
        let precise = self.opts.precise_float || matches!(self.core, Core::Exact);
        for (sp, dp) in s.chunks_exact(ss).zip(d.chunks_exact_mut(ds)) {
            if precise {
                self.exact(sp, dp);
            } else {
                self.fast(sp, dp);
            }
            for e in 0..extra {
                dp[self.outputs + e] = sp[self.inputs + e];
            }
        }
    }

    /// In-place float conversion of whole pixels of `stride ≥ max(inputs, outputs)` values.
    pub fn apply(&self, buf: &mut [f32], stride: usize) -> Result<(), CmsError> {
        let need = self.inputs.max(self.outputs);
        if stride < need {
            return Err(invalid_arg(format!("pixel stride {stride} is below the channel count {need}")));
        }
        if !buf.len().is_multiple_of(stride) {
            return Err(invalid_arg(format!("buffer of {} values is not whole pixels of {stride}", buf.len())));
        }
        let precise = self.opts.precise_float || matches!(self.core, Core::Exact);
        let body = |chunk: &mut [f32]| {
            let mut out = [0.0f32; 16];
            for p in chunk.chunks_exact_mut(stride) {
                if precise {
                    self.exact(p, &mut out);
                } else {
                    self.fast(p, &mut out);
                }
                p[..self.outputs].copy_from_slice(&out[..self.outputs]);
            }
        };
        par_inplace(buf, stride, body);
        Ok(())
    }

    /// Converts 8-bit pixels (see [`Transform::convert_f32`] for the layout).
    pub fn convert_u8(&self, src: &[u8], src_stride: usize, dst: &mut [u8], dst_stride: usize, copy_extra: bool) -> Result<(), CmsError> {
        self.check_buffers(src.len(), src_stride, dst.len(), dst_stride)?;
        let extra = extra_channels(self, src_stride, dst_stride, copy_extra);
        par_rows(src, src_stride, dst, dst_stride, |s, d| self.run_u8(s, src_stride, d, dst_stride, extra));
        Ok(())
    }

    fn run_u8(&self, s: &[u8], ss: usize, d: &mut [u8], ds: usize, extra: usize) {
        let q = |v: f32| (v.clamp(0.0, 1.0) * 255.0 + 0.5) as u8;
        let mut mid = [0.0f32; 16];
        let mut out = [0.0f32; 16];
        match &self.core {
            Core::Grid(c) if c.inputs == 3 => {
                let (g0, g1, g2) = (&self.grid8[0], &self.grid8[1], &self.grid8[2]);
                let s1 = c.outputs * c.grid[2];
                let s0 = s1 * c.grid[1];
                let white8 = self.white.as_ref().and_then(|w| Some((w.src8.as_deref()?, &w.dst[..])));
                for (sp, dp) in s.chunks_exact(ss).zip(d.chunks_exact_mut(ds)) {
                    match white8 {
                        Some((w8, w)) if sp[..3] == *w8 => out[..self.outputs].copy_from_slice(w),
                        _ => {
                            let (a, b, cc) = (g0[sp[0] as usize], g1[sp[1] as usize], g2[sp[2] as usize]);
                            c.tetra_pub((a.0 + b.0 + cc.0) as usize, [s0, s1, c.outputs], [a.1, b.1, cc.1], &mut mid);
                            self.finish(&mid, &mut out);
                        }
                    }
                    for k in 0..self.outputs {
                        dp[k] = q(out[k]);
                    }
                    for e in 0..extra {
                        dp[self.outputs + e] = sp[self.inputs + e];
                    }
                }
            }
            Core::Grid(c) if c.inputs == 4 => {
                let g = &self.grid8;
                let s3 = c.outputs;
                let s2 = s3 * c.grid[3];
                let s1 = s2 * c.grid[2];
                let s0 = s1 * c.grid[1];
                let mut hi = [0.0f32; 16];
                for (sp, dp) in s.chunks_exact(ss).zip(d.chunks_exact_mut(ds)) {
                    let (a, b, cc, k4) = (g[0][sp[0] as usize], g[1][sp[1] as usize], g[2][sp[2] as usize], g[3][sp[3] as usize]);
                    let base = (a.0 + b.0 + cc.0 + k4.0) as usize;
                    c.tetra_pub(base, [s1, s2, s3], [b.1, cc.1, k4.1], &mut mid);
                    if a.1 > 0.0 {
                        c.tetra_pub(base + s0, [s1, s2, s3], [b.1, cc.1, k4.1], &mut hi);
                        for k in 0..c.outputs {
                            mid[k] += (hi[k] - mid[k]) * a.1;
                        }
                    }
                    self.finish(&mid, &mut out);
                    for k in 0..self.outputs {
                        dp[k] = q(out[k]);
                    }
                    for e in 0..extra {
                        dp[self.outputs + e] = sp[self.inputs + e];
                    }
                }
            }
            Core::Shaper { curves8, matrix: Some((3, 3, m, off)), .. } if self.shaper3_ok(curves8.len(), m.len(), off.len()) => {
                let load = |sp: &[u8]| [curves8[0][sp[0] as usize], curves8[1][sp[1] as usize], curves8[2][sp[2] as usize]];
                self.shaper3_blocks(s, ss, d, ds, extra, m, off, load, q);
            }
            Core::Shaper { curves8, matrix, .. } => {
                let mut v = [0.0f32; 16];
                for (sp, dp) in s.chunks_exact(ss).zip(d.chunks_exact_mut(ds)) {
                    for k in 0..self.inputs {
                        v[k] = if curves8.is_empty() { sp[k] as f32 / 255.0 } else { curves8[k][sp[k] as usize] };
                    }
                    apply_matrix(matrix, &v, &mut mid, self.inputs);
                    self.finish(&mid, &mut out);
                    for k in 0..self.outputs {
                        dp[k] = q(out[k]);
                    }
                    for e in 0..extra {
                        dp[self.outputs + e] = sp[self.inputs + e];
                    }
                }
            }
            Core::Exact if !self.exact8.is_empty() => {
                let mut inp = [0.0f32; 16];
                for (sp, dp) in s.chunks_exact(ss).zip(d.chunks_exact_mut(ds)) {
                    let mut v = [0.0f32; 16];
                    for k in 0..self.inputs {
                        inp[k] = sp[k] as f32 / 255.0;
                        v[k] = self.exact8[k][sp[k] as usize];
                    }
                    // Device white stays exact, as on the per-pixel exact path.
                    match self.white_for(&inp) {
                        Some(w) => out[..self.outputs].copy_from_slice(w),
                        None => self.pipeline.eval_after_first(&mut v, &mut out),
                    }
                    for k in 0..self.outputs {
                        dp[k] = q(out[k]);
                    }
                    for e in 0..extra {
                        dp[self.outputs + e] = sp[self.inputs + e];
                    }
                }
            }
            _ => {
                let mut inp = [0.0f32; 16];
                for (sp, dp) in s.chunks_exact(ss).zip(d.chunks_exact_mut(ds)) {
                    for k in 0..self.inputs {
                        inp[k] = sp[k] as f32 / 255.0;
                    }
                    self.fast(&inp, &mut out);
                    for k in 0..self.outputs {
                        dp[k] = q(out[k]);
                    }
                    for e in 0..extra {
                        dp[self.outputs + e] = sp[self.inputs + e];
                    }
                }
            }
        }
    }

    fn shaper3_ok(&self, curves: usize, m: usize, off: usize) -> bool {
        self.inputs == 3
            && self.outputs == 3
            && curves == 3
            && m == 9
            && off == 3
            && self.post.as_ref().is_some_and(|p| p.tables.len() == 3 && p.tables.iter().all(|t| t.len() == OUT_TABLE))
    }

    /// Matrix/TRC transform of whole pixels in blocks: the matrix and the output-table index
    /// maths run over a block of pixels at a time (same operations, same order, per pixel).
    #[allow(clippy::too_many_arguments)]
    fn shaper3_blocks<S: Copy>(
        &self,
        s: &[S],
        ss: usize,
        d: &mut [S],
        ds: usize,
        extra: usize,
        m: &[f32],
        off: &[f32],
        load: impl Fn(&[S]) -> [f32; 3],
        store: impl Fn(f32) -> S,
    ) {
        const B: usize = 64;
        let Some(post) = self.post.as_ref() else { return };
        let tabs = [&post.tables[0][..OUT_TABLE], &post.tables[1][..OUT_TABLE], &post.tables[2][..OUT_TABLE]];
        let m: [f32; 9] = std::array::from_fn(|i| m[i]);
        let off: [f32; 3] = std::array::from_fn(|i| off[i]);
        let n_px = (s.len() / ss).min(d.len() / ds);
        let mut x = [[0.0f32; B]; 3];
        let mut start = 0;
        while start < n_px {
            let n = B.min(n_px - start);
            let sb = &s[start * ss..(start + n) * ss];
            let db = &mut d[start * ds..(start + n) * ds];
            for (j, sp) in sb.chunks_exact(ss).enumerate() {
                let [a, b, c] = load(sp);
                for r in 0..3 {
                    let mut acc = off[r];
                    acc += m[r * 3] * a;
                    acc += m[r * 3 + 1] * b;
                    acc += m[r * 3 + 2] * c;
                    x[r][j] = acc;
                }
            }
            for (ch, t) in tabs.iter().enumerate() {
                let mut u = [0.0f32; B];
                let mut ix = [0i32; B];
                for j in 0..B {
                    let v = x[ch][j].clamp(0.0, 1.0).sqrt().sqrt() * (OUT_TABLE - 1) as f32;
                    let i = (v as i32).min(OUT_TABLE as i32 - 2);
                    ix[j] = i;
                    u[j] = v - i as f32;
                }
                for (j, dp) in db.chunks_exact_mut(ds).enumerate() {
                    let i = ix[j] as usize & (OUT_TABLE - 1);
                    let (v0, v1) = (t[i], t[(i + 1) & (OUT_TABLE - 1)]);
                    dp[ch] = store(v0 + (v1 - v0) * u[j]);
                }
            }
            for (sp, dp) in sb.chunks_exact(ss).zip(db.chunks_exact_mut(ds)) {
                dp[3..3 + extra].copy_from_slice(&sp[3..3 + extra]);
            }
            start += n;
        }
    }

    /// Converts 16-bit pixels (see [`Transform::convert_f32`] for the layout).
    pub fn convert_u16(&self, src: &[u16], src_stride: usize, dst: &mut [u16], dst_stride: usize, copy_extra: bool) -> Result<(), CmsError> {
        self.check_buffers(src.len(), src_stride, dst.len(), dst_stride)?;
        let extra = extra_channels(self, src_stride, dst_stride, copy_extra);
        par_rows(src, src_stride, dst, dst_stride, |s, d| self.run_u16(s, src_stride, d, dst_stride, extra));
        Ok(())
    }

    fn run_u16(&self, s: &[u16], ss: usize, d: &mut [u16], ds: usize, extra: usize) {
        if let Core::Shaper { curves, matrix: Some((3, 3, m, off)), .. } = &self.core
            && self.shaper3_ok(curves.len(), m.len(), off.len())
        {
            let load = |sp: &[u16]| {
                [table(&curves[0], sp[0] as f32 / 65535.0), table(&curves[1], sp[1] as f32 / 65535.0), table(&curves[2], sp[2] as f32 / 65535.0)]
            };
            return self.shaper3_blocks(s, ss, d, ds, extra, m, off, load, |v: f32| (v.clamp(0.0, 1.0) * 65535.0 + 0.5) as u16);
        }
        let mut inp = [0.0f32; 16];
        let mut out = [0.0f32; 16];
        for (sp, dp) in s.chunks_exact(ss).zip(d.chunks_exact_mut(ds)) {
            for k in 0..self.inputs {
                inp[k] = sp[k] as f32 / 65535.0;
            }
            self.fast(&inp, &mut out);
            for k in 0..self.outputs {
                dp[k] = (out[k].clamp(0.0, 1.0) * 65535.0 + 0.5) as u16;
            }
            for e in 0..extra {
                dp[self.outputs + e] = sp[self.inputs + e];
            }
        }
    }

    /// Converts many independent buffers of native-endian samples in one parallel pass (e.g.
    /// all tiles of a surface): far cheaper than one parallel call per small buffer. Each job
    /// is `(src bytes, dst bytes)` holding whole pixels of `src_stride`/`dst_stride` samples,
    /// as in [`Transform::convert_f32`]. Every job is checked before any is converted.
    pub fn convert_bytes_many(
        &self,
        sample: SampleKind,
        jobs: Vec<(&[u8], &mut [u8])>,
        src_stride: usize,
        dst_stride: usize,
        copy_extra: bool,
    ) -> Result<(), CmsError> {
        let size = match sample {
            SampleKind::U8 => 1,
            SampleKind::U16 => 2,
            SampleKind::F32 => 4,
        };
        for (s, d) in &jobs {
            if !s.len().is_multiple_of(size) || !d.len().is_multiple_of(size) {
                return Err(invalid_arg(format!("byte buffers of {}/{} bytes are not whole {size}-byte samples", s.len(), d.len())));
            }
            self.check_buffers(s.len() / size, src_stride, d.len() / size, dst_stride)?;
        }
        let extra = extra_channels(self, src_stride, dst_stride, copy_extra);
        let run = |(s, d): (&[u8], &mut [u8])| match sample {
            SampleKind::U8 => self.run_u8(s, src_stride, d, dst_stride, extra),
            SampleKind::U16 => {
                let a: Vec<u16> = s.as_chunks::<2>().0.iter().map(|b| u16::from_ne_bytes([b[0], b[1]])).collect();
                let mut o: Vec<u16> = d.as_chunks::<2>().0.iter().map(|b| u16::from_ne_bytes([b[0], b[1]])).collect();
                self.run_u16(&a, src_stride, &mut o, dst_stride, extra);
                for (b, v) in d.as_chunks_mut::<2>().0.iter_mut().zip(o) {
                    b.copy_from_slice(&v.to_ne_bytes());
                }
            }
            SampleKind::F32 => {
                let a: Vec<f32> = s.as_chunks::<4>().0.iter().map(|b| f32::from_ne_bytes([b[0], b[1], b[2], b[3]])).collect();
                let mut o: Vec<f32> = d.as_chunks::<4>().0.iter().map(|b| f32::from_ne_bytes([b[0], b[1], b[2], b[3]])).collect();
                self.run_f32(&a, src_stride, &mut o, dst_stride, extra);
                for (b, v) in d.as_chunks_mut::<4>().0.iter_mut().zip(o) {
                    b.copy_from_slice(&v.to_ne_bytes());
                }
            }
        };
        #[cfg(not(target_arch = "wasm32"))]
        {
            use rayon::prelude::*;
            jobs.into_par_iter().for_each(run);
        }
        #[cfg(target_arch = "wasm32")]
        jobs.into_iter().for_each(run);
        Ok(())
    }

    /// Checks a conversion's buffers once per call, so the per-pixel loops need no checks.
    fn check_buffers(&self, src_len: usize, ss: usize, dst_len: usize, ds: usize) -> Result<(), CmsError> {
        if ss < self.inputs || ds < self.outputs {
            return Err(invalid_arg(format!("pixel strides {ss}/{ds} are below the channel counts {}/{}", self.inputs, self.outputs)));
        }
        if !src_len.is_multiple_of(ss) || !dst_len.is_multiple_of(ds) {
            return Err(invalid_arg(format!("buffers of {src_len}/{dst_len} samples are not whole pixels of {ss}/{ds} samples")));
        }
        if dst_len / ds < src_len / ss {
            return Err(invalid_arg(format!("destination holds {} pixels, source {}", dst_len / ds, src_len / ss)));
        }
        Ok(())
    }
}

/// Normalised device values of white: full additive signal, no ink, or Lab L* 100,
/// a* = b* = 0 in the v4 encoding Lab device values use.
fn device_white(cs: ColorSpace) -> Option<Vec<f32>> {
    Some(match cs {
        ColorSpace::Rgb => vec![1.0; 3],
        ColorSpace::Gray => vec![1.0],
        ColorSpace::Cmyk => vec![0.0; 4],
        ColorSpace::Cmy => vec![0.0; 3],
        ColorSpace::Lab => vec![1.0, 128.0 / 255.0, 128.0 / 255.0],
        _ => return None,
    })
}

fn extra_channels(t: &Transform, ss: usize, ds: usize, copy: bool) -> usize {
    if copy { (ss - t.inputs).min(ds - t.outputs) } else { 0 }
}

fn shaper(cs: &[Curve], matrix: Option<(usize, usize, &Vec<f64>, &Vec<f64>)>) -> Core {
    let curves8 = cs.iter().map(|c| std::array::from_fn(|i| c.eval64(i as f64 / 255.0) as f32)).collect();
    let curves = cs.iter().map(|c| (0..=TABLE_1D).map(|i| c.eval64(i as f64 / TABLE_1D as f64) as f32).collect()).collect();
    let matrix = matrix.map(|(r, c, m, o)| (r, c, m.iter().map(|v| *v as f32).collect(), o.iter().map(|v| *v as f32).collect()));
    Core::Shaper { curves8, curves, matrix }
}

#[inline]
fn table(t: &[f32], x: f32) -> f32 {
    // Sampled over [0, 1] with TABLE_1D segments; values outside extrapolate linearly from
    // the end segments so HDR floats stay usable in the fast path.
    let n = t.len() - 1;
    let p = x * n as f32;
    let i = (p.floor().max(0.0) as usize).min(n - 1);
    let f = p - i as f32;
    t[i] + (t[i + 1] - t[i]) * f
}

#[inline]
fn apply_matrix(m: &Option<(usize, usize, Vec<f32>, Vec<f32>)>, v: &[f32; 16], out: &mut [f32; 16], n: usize) {
    match m {
        None => out[..n].copy_from_slice(&v[..n]),
        Some((rows, cols, m, off)) => {
            for r in 0..*rows {
                let mut s = off[r];
                for c in 0..*cols {
                    s += m[r * cols + c] * v[c];
                }
                out[r] = s;
            }
        }
    }
}

// ------------------------------------------------------------ parallel helpers

#[cfg(not(target_arch = "wasm32"))]
const PAR_MIN_PIXELS: usize = 16 * 1024;
#[cfg(not(target_arch = "wasm32"))]
const PAR_CHUNK_PIXELS: usize = 8 * 1024;

#[cfg(not(target_arch = "wasm32"))]
fn par_rows<S: Sync, D: Send>(src: &[S], ss: usize, dst: &mut [D], ds: usize, f: impl Fn(&[S], &mut [D]) + Sync) {
    use rayon::prelude::*;
    let n = (src.len() / ss).min(dst.len() / ds);
    let (src, dst) = (&src[..n * ss], &mut dst[..n * ds]);
    if n < PAR_MIN_PIXELS {
        return f(src, dst);
    }
    src.par_chunks(PAR_CHUNK_PIXELS * ss).zip(dst.par_chunks_mut(PAR_CHUNK_PIXELS * ds)).for_each(|(s, d)| f(s, d));
}

#[cfg(target_arch = "wasm32")]
fn par_rows<S: Sync, D: Send>(src: &[S], ss: usize, dst: &mut [D], ds: usize, f: impl Fn(&[S], &mut [D]) + Sync) {
    let n = (src.len() / ss).min(dst.len() / ds);
    f(&src[..n * ss], &mut dst[..n * ds]);
}

#[cfg(not(target_arch = "wasm32"))]
fn par_inplace(buf: &mut [f32], stride: usize, f: impl Fn(&mut [f32]) + Sync) {
    use rayon::prelude::*;
    if buf.len() / stride < PAR_MIN_PIXELS {
        return f(buf);
    }
    buf.par_chunks_mut(PAR_CHUNK_PIXELS * stride).for_each(&f);
}

#[cfg(target_arch = "wasm32")]
fn par_inplace(buf: &mut [f32], _stride: usize, f: impl Fn(&mut [f32]) + Sync) {
    f(buf)
}

// ------------------------------------------------------------ linking

/// Stages from `src` device values to `dst` device values.
pub(crate) fn link_stages(src: &Profile, dst: &Profile, intent: Intent, bpc: bool) -> Result<Vec<Stage>, CmsError> {
    for p in [src, dst] {
        // A Lab DeviceLink stands in for an abstract profile only when it also ends in Lab.
        let lab_to_lab =
            p.color_space == ColorSpace::Lab && p.a2b.iter().flatten().all(|l| l.outputs == 3) && p.b2a.iter().flatten().all(|l| l.inputs == 3);
        if matches!(p.class, ProfileClass::DeviceLink | ProfileClass::NamedColor | ProfileClass::Abstract) && !lab_to_lab {
            return Err(CmsError::Unsupported(format!("{:?} profiles cannot be used as a source or destination", p.class)));
        }
    }
    let (mut stages, s_pcs) = src.device_to_pcs(intent)?;
    let (d_stages, d_pcs) = dst.pcs_to_device(intent)?;
    let adjust = connection(src, dst, intent, bpc)?;
    match adjust {
        Some((scale, offset)) => {
            if s_pcs == Pcs::Lab {
                stages.push(Stage::LabToXyz);
            }
            stages.push(Stage::Matrix {
                rows: 3,
                cols: 3,
                m: vec![scale[0], 0.0, 0.0, 0.0, scale[1], 0.0, 0.0, 0.0, scale[2]],
                offset: offset.to_vec(),
            });
            if d_pcs == Pcs::Lab {
                stages.push(Stage::XyzToLab);
            }
        }
        None => match (s_pcs, d_pcs) {
            (Pcs::Lab, Pcs::Xyz) => stages.push(Stage::LabToXyz),
            (Pcs::Xyz, Pcs::Lab) => stages.push(Stage::XyzToLab),
            _ => {}
        },
    }
    stages.extend(d_stages);
    Ok(stages)
}

/// Per-component XYZ scale and offset applied in the PCS: absolute colorimetric media-white
/// scaling or black point compensation. `None` = identity.
/// Per-component XYZ `(scale, offset)`.
type Affine3 = ([f64; 3], [f64; 3]);

fn connection(src: &Profile, dst: &Profile, intent: Intent, bpc: bool) -> Result<Option<Affine3>, CmsError> {
    if intent == Intent::AbsoluteColorimetric {
        let (ws, wd) = (src.media_white(), dst.media_white());
        // `wtpt` comes straight from the file. A zero, negative or non-finite component would
        // make the scale NaN or infinite and poison every colour the transform produces, so a
        // profile with such a media white is refused for this intent instead.
        let usable = |w: &[f64; 3]| w.iter().all(|v| v.is_finite() && *v > 0.0);
        if !usable(&ws) || !usable(&wd) {
            return Err(CmsError::Invalid(format!("absolute colorimetric needs a positive media white point (source {ws:?}, destination {wd:?})")));
        }
        let s = [ws[0] / wd[0], ws[1] / wd[1], ws[2] / wd[2]];
        if s.iter().all(|v| (v - 1.0).abs() < 1e-6) {
            return Ok(None);
        }
        return Ok(Some((s, [0.0; 3])));
    }
    if !bpc {
        return Ok(None);
    }
    let bs = black_point(src, intent, true)?;
    let bd = black_point(dst, intent, false)?;
    if (bs - bd).abs() < 1e-6 {
        return Ok(None);
    }
    // Neutral black points (luminance only, as in Adobe's BPC), mapped linearly in XYZ so that
    // the source black lands on the destination black and the white stays fixed.
    let w = math::D50;
    let mut scale = [1.0; 3];
    let mut offset = [0.0; 3];
    for i in 0..3 {
        let (s, d) = (bs * w[i], bd * w[i]);
        let a = (w[i] - d) / (w[i] - s);
        scale[i] = a;
        offset[i] = w[i] * (1.0 - a);
    }
    Ok(Some((scale, offset)))
}

/// Luminance (Y, D50-relative) of the profile's black point for `intent`.
///
/// Follows the structure of Adobe's published BPC algorithm: matrix/TRC profiles use device
/// zero; v4 LUT profiles in perceptual/saturation use the perceptual reference black; other LUT
/// profiles use the darkest colour the device reaches (Lab 0 round-tripped through BToA/AToB,
/// or device black for additive spaces). The destination-side curve fitting of the paper is
/// not implemented; the round-trip black is used directly.
pub fn black_point(p: &Profile, intent: Intent, _input: bool) -> Result<f64, CmsError> {
    if matches!(p.color_space, ColorSpace::Lab | ColorSpace::Xyz) {
        return Ok(0.0);
    }
    let lut_based = p.a2b.iter().any(Option::is_some);
    if lut_based && p.version.0 >= 4 && matches!(intent, Intent::Perceptual | Intent::Saturation) {
        return Ok(PERCEPTUAL_BLACK[1]);
    }
    let pcs_y = |st: Vec<Stage>, pcs: Pcs, dev: &[f32]| -> Result<f64, CmsError> {
        let pipe = Pipeline::new(p.channels(), st)?;
        let mut o = [0.0f32; 16];
        pipe.eval(dev, &mut o)?;
        let y = match pcs {
            Pcs::Xyz => o[1] as f64,
            Pcs::Lab => math::lab_to_xyz([o[0] as f64, 0.0, 0.0], math::D50)[1],
        };
        Ok(y.clamp(0.0, 0.5))
    };
    let (to_pcs, pcs) = p.device_to_pcs(intent)?;
    let additive = matches!(p.color_space, ColorSpace::Rgb | ColorSpace::Gray);
    if additive || !lut_based {
        let zero = vec![0.0f32; p.channels()];
        return pcs_y(to_pcs, pcs, &zero);
    }
    // Subtractive LUT profile: darkest reachable colour.
    let colorimetric = if intent == Intent::Saturation { Intent::Saturation } else { Intent::RelativeColorimetric };
    match p.pcs_to_device(colorimetric) {
        Ok((from_pcs, fpcs)) => {
            let mut st = Vec::new();
            if fpcs == Pcs::Xyz {
                st.push(Stage::LabToXyz);
            }
            st.extend(from_pcs);
            let pipe = Pipeline::new(3, st)?;
            let mut dev = [0.0f32; 16];
            pipe.eval(&[0.0, 0.0, 0.0], &mut dev)?;
            pcs_y(to_pcs, pcs, dev.get(..p.channels()).unwrap_or(&dev))
        }
        Err(_) => {
            let ones = vec![1.0f32; p.channels()];
            pcs_y(to_pcs, pcs, &ones)
        }
    }
}

// ------------------------------------------------------------ cache

type Key = (u64, u64, TransformOptions);

/// A cached transform with the ICC bytes of the profiles it was built from.
struct Entry {
    src: Arc<Vec<u8>>,
    dst: Arc<Vec<u8>>,
    transform: Arc<Transform>,
}

/// Process-wide cache of transforms keyed by profile content and options. The key is the
/// profiles' 64-bit hashes, which a crafted profile can share with another, so a hit counts only
/// when the bytes match too; a colliding pair gets its own, uncached transform.
pub fn cached(src: &Profile, dst: &Profile, opts: TransformOptions) -> Result<Arc<Transform>, CmsError> {
    static CACHE: OnceLock<Mutex<HashMap<Key, Entry>>> = OnceLock::new();
    let key = (src.content_hash(), dst.content_hash(), opts);
    let (src_bytes, dst_bytes) = (src.to_bytes(), dst.to_bytes());
    let cache = CACHE.get_or_init(Default::default);
    let hit =
        cache.lock().unwrap_or_else(|e| e.into_inner()).get(&key).map(|e| (e.src == src_bytes && e.dst == dst_bytes).then(|| e.transform.clone()));
    match hit {
        Some(Some(t)) => return Ok(t),
        Some(None) => return Ok(Arc::new(Transform::with_options(src, dst, opts)?)),
        None => {}
    }
    let t = Arc::new(Transform::with_options(src, dst, opts)?);
    let mut c = cache.lock().unwrap_or_else(PoisonError::into_inner);
    if c.len() > 64 {
        c.clear();
    }
    c.insert(key, Entry { src: src_bytes, dst: dst_bytes, transform: t.clone() });
    Ok(t)
}

impl Clut {
    /// Tetrahedral interpolation in one cube (crate-internal fast paths).
    #[inline]
    pub(crate) fn tetra_pub(&self, base: usize, s: [usize; 3], r: [f32; 3], out: &mut [f32; 16]) {
        let d = &self.data;
        let [rx, ry, rz] = r;
        let [sx, sy, sz] = s;
        let (p1, p2, p3, w1, w2, w3) = if rx >= ry {
            if ry >= rz {
                (sx, sx + sy, sx + sy + sz, rx, ry, rz)
            } else if rx >= rz {
                (sx, sx + sz, sx + sy + sz, rx, rz, ry)
            } else {
                (sz, sx + sz, sx + sy + sz, rz, rx, ry)
            }
        } else if rx >= rz {
            (sy, sx + sy, sx + sy + sz, ry, rx, rz)
        } else if ry >= rz {
            (sy, sy + sz, sx + sy + sz, ry, rz, rx)
        } else {
            (sz, sy + sz, sx + sy + sz, rz, ry, rx)
        };
        let o = self.outputs.min(MAX_CHANNELS);
        for k in 0..o {
            let v0 = d[base + k];
            let v1 = d[base + p1 + k];
            let v2 = d[base + p2 + k];
            let v3 = d[base + p3 + k];
            out[k] = v0 + (v1 - v0) * w1 + (v2 - v1) * w2 + (v3 - v2) * w3;
        }
    }
}