alice-sdf 5.0.0

//! Codec Bridge: ALICE-SDF × ALICE-Codec Integration
//!
//! Evaluate SDF fields on voxel grids and encode the resulting
//! 3D volume data using ALICE-Codec's wavelet compression.
//!
//! # Architecture
//!
//! ```text
//! SdfNode ──→ voxelize_sdf() ──→ SdfVolume ──→ encode_sdf_volume()
//!              (eval_grid)        (f32 data)    ├── f32 → i32 (fixed-point)
//!                                               ├── Wavelet3D forward
//!                                               ├── Quantization (RDO)
//!                                               └── rANS entropy coding
//! ```
//!
//! # Bitstream
//!
//! All integers little-endian:
//!
//! | bytes | field |
//! |---|---|
//! | 12 | width, height, depth (u32) |
//! | 12 | origin x, y, z (f32) |
//! | 4 + 4 | voxel size along x, fixed_point_scale (f32) |
//! | 1 | quality (u8) |
//! | 1 | flags: bit 0 CDF 5/3 wavelet, bit 1 rANS payload, bit 2 i32 coefficients, bit 3 per-axis voxel size; bits 4-7 reserved (must be 0, the decoder rejects them) |
//! | 4 | symbol count (u32, voxels) |
//! | 4 + 4 | voxel size along y, z (f32), only with bit 3 |
//! | 4 + 4 | quantizer step, dead zone (i32) |
//! | 1024 | byte histogram (256 x u32), only with bit 1 |
//! | 4 + n | payload length, payload |
//!
//! The payload (rANS-coded or raw) is the quantized coefficients as i16, or as
//! i32 when bit 2 is set. The encoder sets bit 2 only when a coefficient does
//! not fit in i16, so such a stream is otherwise identical to the i16 format.
//!
//! Bit 3 (5.0.0) carries a per-axis voxel size. The encoder sets it only when
//! the three steps differ (bit-wise), so a cubic-voxel stream is byte-identical
//! to the earlier format, and a stream without bit 3 (every stream written
//! before 5.0.0) decodes with the same step on all three axes. A decoder older
//! than 5.0.0 that checks the reserved bits rejects a bit-3 stream instead of
//! reading it with the wrong layout.
//!
//! # Example
//!
//! ```ignore
//! use alice_sdf::prelude::*;
//! use alice_sdf::codec_bridge::{voxelize_sdf, encode_sdf_volume, decode_sdf_volume, EncodeConfig};
//!
//! // Create SDF shape
//! let node = SdfNode::sphere(1.0)
//!     .smooth_union(SdfNode::box3d(0.5, 2.0, 0.5), 0.2);
//!
//! // Voxelize on 64^3 grid
//! let volume = voxelize_sdf(
//!     &node,
//!     Vec3::splat(-2.0),
//!     Vec3::splat(2.0),
//!     [64, 64, 64],
//! );
//!
//! // Encode with wavelet compression
//! let encoded = encode_sdf_volume(&volume, &EncodeConfig::default());
//!
//! // Decode back
//! let decoded = decode_sdf_volume(&encoded);
//! ```
//!
//! Author: Moroya Sakamoto

use alice_codec::quant::Quantizer;
use alice_codec::rans::{FrequencyTable, RansDecoder, RansEncoder};
use alice_codec::wavelet::Wavelet3D;

use crate::eval::eval_batch_parallel;
use crate::types::SdfNode;
use glam::Vec3;

// ============================================================================
// Volume Types
// ============================================================================

/// Voxelized SDF volume for codec encoding.
///
/// Distance values are stored in row-major order: `data[z * H * W + y * W + x]`
/// matching ALICE-Codec's 3D wavelet layout (x, y, z/depth).
#[derive(Clone, Debug)]
pub struct SdfVolume {
    /// Distance values in row-major order `[z][y][x]`.
    pub data: Vec<f32>,
    /// Volume width (X axis).
    pub width: usize,
    /// Volume height (Y axis).
    pub height: usize,
    /// Volume depth (Z axis).
    pub depth: usize,
    /// World-space minimum corner.
    pub origin: Vec3,
    /// Distance between adjacent voxels along each axis.
    ///
    /// Until 5.0.0 this was one `f32` (the smallest of the three steps), so
    /// on a grid whose steps differ [`Self::world_pos`] was wrong along the
    /// other axes. For a cubic voxel use `Vec3::splat(step)`.
    pub voxel_size: Vec3,
}

impl SdfVolume {
    /// Total number of voxels.
    #[inline]
    pub const fn len(&self) -> usize {
        self.width * self.height * self.depth
    }

    /// Check if volume is empty.
    #[inline]
    pub fn is_empty(&self) -> bool {
        self.data.is_empty()
    }

    /// Sample a voxel at (x, y, z) indices.
    #[inline]
    pub fn get(&self, x: usize, y: usize, z: usize) -> f32 {
        self.data[z * self.height * self.width + y * self.width + x]
    }

    /// World-space position for voxel indices:
    /// `origin + (x·voxel_size.x, y·voxel_size.y, z·voxel_size.z)`.
    #[inline]
    pub fn world_pos(&self, x: usize, y: usize, z: usize) -> Vec3 {
        self.origin
            + Vec3::new(
                x as f32 * self.voxel_size.x,
                y as f32 * self.voxel_size.y,
                z as f32 * self.voxel_size.z,
            )
    }
}

/// Statistics about a voxelized SDF volume.
#[derive(Clone, Debug)]
pub struct VolumeStats {
    /// Minimum distance value in the volume.
    pub min_distance: f32,
    /// Maximum distance value in the volume.
    pub max_distance: f32,
    /// Number of sign changes between adjacent voxels (surface crossings).
    pub zero_crossings: usize,
    /// Total number of voxels.
    pub total_voxels: usize,
    /// Ratio of surface voxels to total voxels.
    pub surface_ratio: f64,
    /// Root-mean-square of distance values.
    pub rms: f64,
}

/// Configuration for SDF volume encoding.
#[derive(Clone, Debug)]
pub struct EncodeConfig {
    /// Quality factor (0-100). Higher = better quality, larger output.
    pub quality: u8,
    /// Fixed-point scale factor: `i32_value = (f32_value * scale).round()`.
    /// Higher scale preserves more precision but increases coefficient magnitude.
    pub fixed_point_scale: f32,
    /// Use CDF 5/3 wavelet (lossless integers) instead of CDF 9/7.
    pub lossless_wavelet: bool,
}

impl Default for EncodeConfig {
    fn default() -> Self {
        Self {
            quality: 75,
            fixed_point_scale: 1024.0,
            lossless_wavelet: false,
        }
    }
}

impl EncodeConfig {
    /// High quality preset (quality 95).
    pub const fn high_quality() -> Self {
        Self {
            quality: 95,
            fixed_point_scale: 4096.0,
            lossless_wavelet: true,
        }
    }

    /// Fast / small preset (quality 50).
    pub const fn fast() -> Self {
        Self {
            quality: 50,
            fixed_point_scale: 512.0,
            lossless_wavelet: false,
        }
    }

    /// Lossless wavelet preset (CDF 5/3).
    pub const fn lossless() -> Self {
        Self {
            quality: 100,
            fixed_point_scale: 4096.0,
            lossless_wavelet: true,
        }
    }
}

/// Header stored at the beginning of an encoded SDF volume bitstream.
#[derive(Clone, Debug)]
struct EncodedHeader {
    width: u32,
    height: u32,
    depth: u32,
    origin: [f32; 3],
    /// Voxel size along x (along every axis without `FLAG_ANISOTROPIC_VOXEL`).
    voxel_size: f32,
    fixed_point_scale: f32,
    quality: u8,
    /// Flags byte: bit 0 = lossless_wavelet, bit 1 = rans_compressed,
    /// bit 2 = wide (i32) coefficients. Bits 3-7 are reserved and must be 0.
    flags: u8,
    /// Number of quantized coefficients (voxel count).
    symbol_count: u32,
}

impl EncodedHeader {
    const SIZE: usize = 4 * 3 + 4 * 3 + 4 + 4 + 1 + 1 + 4; // 34 bytes

    const FLAG_LOSSLESS_WAVELET: u8 = 0x01;
    const FLAG_RANS_COMPRESSED: u8 = 0x02;
    /// Coefficients are stored as i32 LE (4 bytes) instead of i16 LE (2 bytes).
    /// Set only when a quantized coefficient does not fit in i16, so a stream
    /// whose coefficients all fit is byte-identical to the original format.
    const FLAG_WIDE_COEFFS: u8 = 0x04;
    /// The voxel sizes along y and z follow the header (2 × f32 LE). Set only
    /// when the three steps differ, so a cubic-voxel stream is byte-identical
    /// to the original format.
    const FLAG_ANISOTROPIC_VOXEL: u8 = 0x08;
    const KNOWN_FLAGS: u8 = Self::FLAG_LOSSLESS_WAVELET
        | Self::FLAG_RANS_COMPRESSED
        | Self::FLAG_WIDE_COEFFS
        | Self::FLAG_ANISOTROPIC_VOXEL;

    fn to_bytes(&self) -> Vec<u8> {
        let mut buf = Vec::with_capacity(Self::SIZE);
        buf.extend_from_slice(&self.width.to_le_bytes());
        buf.extend_from_slice(&self.height.to_le_bytes());
        buf.extend_from_slice(&self.depth.to_le_bytes());
        buf.extend_from_slice(&self.origin[0].to_le_bytes());
        buf.extend_from_slice(&self.origin[1].to_le_bytes());
        buf.extend_from_slice(&self.origin[2].to_le_bytes());
        buf.extend_from_slice(&self.voxel_size.to_le_bytes());
        buf.extend_from_slice(&self.fixed_point_scale.to_le_bytes());
        buf.push(self.quality);
        buf.push(self.flags);
        buf.extend_from_slice(&self.symbol_count.to_le_bytes());
        buf
    }

    fn from_bytes(data: &[u8]) -> Self {
        let read_u32 = |offset: usize| -> u32 {
            u32::from_le_bytes([
                data[offset],
                data[offset + 1],
                data[offset + 2],
                data[offset + 3],
            ])
        };
        let read_f32 = |offset: usize| -> f32 {
            f32::from_le_bytes([
                data[offset],
                data[offset + 1],
                data[offset + 2],
                data[offset + 3],
            ])
        };

        Self {
            width: read_u32(0),
            height: read_u32(4),
            depth: read_u32(8),
            origin: [read_f32(12), read_f32(16), read_f32(20)],
            voxel_size: read_f32(24),
            fixed_point_scale: read_f32(28),
            quality: data[32],
            flags: data[33],
            symbol_count: read_u32(34),
        }
    }

    const fn lossless_wavelet(&self) -> bool {
        self.flags & Self::FLAG_LOSSLESS_WAVELET != 0
    }

    const fn rans_compressed(&self) -> bool {
        self.flags & Self::FLAG_RANS_COMPRESSED != 0
    }

    const fn wide_coeffs(&self) -> bool {
        self.flags & Self::FLAG_WIDE_COEFFS != 0
    }

    const fn anisotropic_voxel(&self) -> bool {
        self.flags & Self::FLAG_ANISOTROPIC_VOXEL != 0
    }
}

// ============================================================================
// Voxelization
// ============================================================================

/// Voxelize an SDF node onto a regular 3D grid.
///
/// Evaluates the SDF at every grid point using parallel batch evaluation.
///
/// # Arguments
/// * `sdf` - The SDF tree to evaluate
/// * `origin` - World-space minimum corner of the volume
/// * `extent` - World-space maximum corner of the volume
/// * `dims` - Grid dimensions `[width, height, depth]`
///
/// # Returns
/// An `SdfVolume` with distance values at each grid point.
pub fn voxelize_sdf(sdf: &SdfNode, origin: Vec3, extent: Vec3, dims: [usize; 3]) -> SdfVolume {
    let [w, h, d] = dims;
    let size = extent - origin;
    let voxel_size_vec = Vec3::new(
        if w > 1 {
            size.x / (w - 1) as f32
        } else {
            size.x
        },
        if h > 1 {
            size.y / (h - 1) as f32
        } else {
            size.y
        },
        if d > 1 {
            size.z / (d - 1) as f32
        } else {
            size.z
        },
    );
    let total = w * h * d;
    let mut points = Vec::with_capacity(total);

    // Build point list in z-major, y-minor, x-innermost order
    // matching Wavelet3D layout: volume[z * H * W + y * W + x]
    for z in 0..d {
        for y in 0..h {
            for x in 0..w {
                points.push(Vec3::new(
                    (x as f32).mul_add(voxel_size_vec.x, origin.x),
                    (y as f32).mul_add(voxel_size_vec.y, origin.y),
                    (z as f32).mul_add(voxel_size_vec.z, origin.z),
                ));
            }
        }
    }

    let data = eval_batch_parallel(sdf, &points);

    SdfVolume {
        data,
        width: w,
        height: h,
        depth: d,
        origin,
        voxel_size: voxel_size_vec,
    }
}

/// Voxelize an SDF using a uniform cubic grid at the given resolution.
///
/// Convenience wrapper that computes dims from a single resolution value.
pub fn voxelize_sdf_uniform(
    sdf: &SdfNode,
    origin: Vec3,
    extent: Vec3,
    resolution: usize,
) -> SdfVolume {
    voxelize_sdf(sdf, origin, extent, [resolution, resolution, resolution])
}

// ============================================================================
// Encoding / Decoding
// ============================================================================

/// Encode a voxelized SDF volume using 3D wavelet transform + rANS.
///
/// Pipeline: f32 -> fixed-point i32 -> Wavelet3D forward -> Quantize -> rANS
///
/// # Arguments
/// * `volume` - The SDF volume to encode
/// * `config` - Encoding parameters (quality, scale, wavelet type)
///
/// # Returns
/// Compressed byte stream including header for decoding.
pub fn encode_sdf_volume(volume: &SdfVolume, config: &EncodeConfig) -> Vec<u8> {
    let total = volume.len();
    let scale = config.fixed_point_scale;

    // 1. Convert f32 distances to fixed-point i32
    let mut coeffs: Vec<i32> = volume
        .data
        .iter()
        .map(|&d| (d * scale).round() as i32)
        .collect();

    // 2. Apply 3D wavelet transform
    let wavelet = if config.lossless_wavelet {
        Wavelet3D::cdf53()
    } else {
        Wavelet3D::cdf97()
    };
    wavelet.forward(&mut coeffs, volume.width, volume.height, volume.depth);

    // 3. Quantize
    // Compute step size from quality: quality 100 -> step 1, dead_zone 0
    // (near-lossless). quality 0 -> aggressive quantization.
    // SDF-specific: wavelet detail coefficients are small but crucial for
    // sign preservation, so dead_zone must be kept very small at high quality.
    let quantizer = {
        if config.quality >= 95 {
            // Near-lossless: no quantization loss, only wavelet rounding
            Quantizer::with_dead_zone(1, 0)
        } else {
            let max_abs = coeffs.iter().map(|c| c.abs()).max().unwrap_or(1).max(1) as f32;
            // Cubic curve: high quality stays near step=1
            let inv_quality = (100 - config.quality.min(100)) as f32 / 100.0;
            let quality_factor = inv_quality * inv_quality * inv_quality;
            let step = (1.0 + quality_factor * (max_abs / 4.0)).round() as i32;
            let step = step.max(1);
            // Dead zone proportional to step but capped to avoid zeroing detail
            let dead_zone = step / 2;
            Quantizer::with_dead_zone(step, dead_zone)
        }
    };

    let mut quantized = vec![0i32; total];
    // `InvalidBufferSize` is impossible: both buffers are `total` long.
    quantizer
        .quantize_buffer(&coeffs, &mut quantized)
        .expect("quantized buffer is sized to the coefficient buffer");

    // 4. Serialize quantized coefficients as little-endian bytes, then
    //    compress those bytes with rANS. This avoids the u8-symbol clamping
    //    issue of alice_codec::quant::to_symbols which is designed for video
    //    pixel data (0-255 range), not for wide-range SDF wavelet coefficients.
    //    i16 when every coefficient fits (the original format), otherwise i32
    //    with FLAG_WIDE_COEFFS: narrowing would silently change the values.
    let wide = quantized.iter().any(|&v| i16::try_from(v).is_err());
    let raw_bytes: Vec<u8> = if wide {
        quantized.iter().flat_map(|&v| v.to_le_bytes()).collect()
    } else {
        quantized
            .iter()
            .flat_map(|&v| (v as i16).to_le_bytes())
            .collect()
    };

    // 5. Build byte-level histogram and compress with rANS.
    //    If the histogram is too skewed (any symbol dominates > 75% of data),
    //    store raw bytes to avoid u32 overflow in the rANS encoder.
    let mut histogram = [0u32; 256];
    for &b in &raw_bytes {
        histogram[b as usize] += 1;
    }

    let max_freq = *histogram.iter().max().unwrap_or(&0);
    // rANS encoder overflows when x_max = ((RANS32_L >> PROB_BITS) << 8) * freq
    // where freq is the *normalized* frequency (max = PROB_SCALE = 4096).
    // The product limit is u32::MAX. With RANS32_L=1<<23, PROB_BITS=12:
    //   ((1<<23 >> 12) << 8) * 4096 = 2^31 which barely overflows.
    // Safe limit: if max_freq dominates > 50% of data, normalized freq approaches
    // PROB_SCALE and can overflow. Skip rANS for small or highly-skewed data.
    let use_rans = raw_bytes.len() >= 512 && max_freq < (raw_bytes.len() as u32 / 2);

    let payload: Vec<u8> = if use_rans {
        let freq_table = FrequencyTable::from_histogram(&histogram);
        let mut encoder = RansEncoder::with_capacity(raw_bytes.len());
        encoder.encode_symbols(&raw_bytes, &freq_table);
        encoder.finish()
    } else {
        // Store raw (rANS can overflow with extreme histograms on small data)
        raw_bytes
    };

    // 6. Assemble output: header + quantizer params + histogram + payload
    let mut flags = 0u8;
    if config.lossless_wavelet {
        flags |= EncodedHeader::FLAG_LOSSLESS_WAVELET;
    }
    if use_rans {
        flags |= EncodedHeader::FLAG_RANS_COMPRESSED;
    }
    if wide {
        flags |= EncodedHeader::FLAG_WIDE_COEFFS;
    }
    let vs = volume.voxel_size;
    let anisotropic = vs.y.to_bits() != vs.x.to_bits() || vs.z.to_bits() != vs.x.to_bits();
    if anisotropic {
        flags |= EncodedHeader::FLAG_ANISOTROPIC_VOXEL;
    }

    let header = EncodedHeader {
        width: volume.width as u32,
        height: volume.height as u32,
        depth: volume.depth as u32,
        origin: [volume.origin.x, volume.origin.y, volume.origin.z],
        voxel_size: vs.x,
        fixed_point_scale: scale,
        quality: config.quality,
        flags,
        symbol_count: total as u32,
    };

    let mut output = header.to_bytes();
    if anisotropic {
        output.extend_from_slice(&vs.y.to_le_bytes());
        output.extend_from_slice(&vs.z.to_le_bytes());
    }

    // Write quantizer step + dead_zone (8 bytes)
    output.extend_from_slice(&quantizer.step.to_le_bytes());
    output.extend_from_slice(&quantizer.dead_zone.to_le_bytes());

    if use_rans {
        // Write histogram (256 * 4 = 1024 bytes) - only needed for rANS decoding
        for &count in &histogram {
            output.extend_from_slice(&count.to_le_bytes());
        }
    }

    // Write payload length + data
    let payload_len = payload.len() as u32;
    output.extend_from_slice(&payload_len.to_le_bytes());
    output.extend(payload);

    output
}

/// Why an encoded SDF volume could not be decoded.
///
/// `#[non_exhaustive]` since 5.0.0: later versions may add reasons.
#[derive(Clone, Debug, PartialEq, Eq)]
#[non_exhaustive]
pub enum DecodeError {
    /// The stream ends before the field that starts at `offset` (`needed`
    /// bytes in total, `len` available).
    Truncated {
        /// Bytes the stream must have to hold the next field.
        needed: usize,
        /// Bytes the stream has.
        len: usize,
    },
    /// The flags byte has bits this decoder does not know (`flags & !0x0f`).
    /// Such a stream was written by a newer encoder; reading it with this
    /// layout would produce wrong values.
    UnknownFlags(u8),
    /// `width · height · depth` overflows or differs from the symbol count in
    /// the header (the encoder writes the voxel count there).
    InvalidDimensions {
        /// Width in the header.
        width: u32,
        /// Height in the header.
        height: u32,
        /// Depth in the header.
        depth: u32,
        /// Symbol count in the header.
        symbol_count: u32,
    },
    /// The quantizer step is below 1 or the dead zone is negative (the encoder
    /// writes `step ≥ 1`, `0 ≤ dead_zone`).
    InvalidQuantizer {
        /// Quantizer step in the stream.
        step: i32,
        /// Dead zone in the stream.
        dead_zone: i32,
    },
    /// The rANS histogram cannot be the one the encoder wrote for this many
    /// bytes: its counts do not sum to the coefficient byte count, a byte
    /// value holds half of the counts or more, there are fewer than 512
    /// bytes (the encoder stores such data raw), or the payload is too short
    /// to hold that many symbols (`bytes > 16 · payload length`).
    InvalidHistogram,
    /// A coefficient leaves the `i32` range while it is dequantized or while
    /// the inverse wavelet transform runs (only a corrupted stream does this;
    /// the reference decoder would overflow there).
    CoefficientOverflow,
}

impl std::fmt::Display for DecodeError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::Truncated { needed, len } => {
                write!(
                    f,
                    "encoded SDF volume truncated: need {needed} bytes, have {len}"
                )
            }
            Self::UnknownFlags(flags) => {
                write!(f, "encoded SDF volume has unknown flag bits 0x{flags:02x}")
            }
            Self::InvalidDimensions {
                width,
                height,
                depth,
                symbol_count,
            } => write!(
                f,
                "encoded SDF volume has dimensions {width}x{height}x{depth} but symbol count {symbol_count}"
            ),
            Self::InvalidQuantizer { step, dead_zone } => write!(
                f,
                "encoded SDF volume has quantizer step {step}, dead zone {dead_zone} (need step >= 1, dead zone >= 0)"
            ),
            Self::InvalidHistogram => {
                write!(f, "encoded SDF volume has an rANS histogram that does not match its payload")
            }
            Self::CoefficientOverflow => {
                write!(f, "encoded SDF volume has a coefficient outside the i32 range")
            }
        }
    }
}

impl std::error::Error for DecodeError {}

/// Decode a compressed SDF volume back to an `SdfVolume`.
///
/// Pipeline: rANS decode -> un-quantize -> Wavelet3D inverse -> i32 -> f32
///
/// # Arguments
/// * `data` - Compressed byte stream produced by `encode_sdf_volume`.
///
/// # Returns
/// Reconstructed `SdfVolume` (with lossy reconstruction error from quantization).
///
/// # Panics
///
/// Panics when the stream cannot be decoded (see [`try_decode_sdf_volume`]
/// for the non-panicking form and the reasons).
pub fn decode_sdf_volume(data: &[u8]) -> SdfVolume {
    try_decode_sdf_volume(data).unwrap_or_else(|e| panic!("decode_sdf_volume: {e}"))
}

/// Decode a compressed SDF volume, reporting malformed input as an error.
///
/// # Errors
///
/// [`DecodeError::Truncated`] when the stream is shorter than its header and
/// payload length say, [`DecodeError::UnknownFlags`] when the flags byte has a
/// reserved bit set, [`DecodeError::InvalidDimensions`] /
/// [`DecodeError::InvalidQuantizer`] / [`DecodeError::InvalidHistogram`] when
/// a header field is outside what the encoder writes, and
/// [`DecodeError::CoefficientOverflow`] when a coefficient would leave the
/// `i32` range. It does not panic on any input, and it allocates at most a
/// constant multiple of `data.len()` (raw payload: the payload holds every
/// coefficient; rANS payload: the histogram bounds the symbols per byte).
pub fn try_decode_sdf_volume(data: &[u8]) -> Result<SdfVolume, DecodeError> {
    let need = |needed: usize| {
        if data.len() < needed {
            Err(DecodeError::Truncated {
                needed,
                len: data.len(),
            })
        } else {
            Ok(())
        }
    };
    let read_u32 = |offset: usize| {
        u32::from_le_bytes([
            data[offset],
            data[offset + 1],
            data[offset + 2],
            data[offset + 3],
        ])
    };

    // 1. Read header + quantizer params (8 bytes)
    need(EncodedHeader::SIZE + 8)?;
    let header = EncodedHeader::from_bytes(data);
    let unknown = header.flags & !EncodedHeader::KNOWN_FLAGS;
    if unknown != 0 {
        return Err(DecodeError::UnknownFlags(unknown));
    }
    let (w, h, d) = (
        header.width as usize,
        header.height as usize,
        header.depth as usize,
    );
    let total = w
        .checked_mul(h)
        .and_then(|wh| wh.checked_mul(d))
        .filter(|&t| t as u64 == u64::from(header.symbol_count))
        .ok_or(DecodeError::InvalidDimensions {
            width: header.width,
            height: header.height,
            depth: header.depth,
            symbol_count: header.symbol_count,
        })?;
    let scale = header.fixed_point_scale;

    let mut offset = EncodedHeader::SIZE;

    // Voxel size: per-axis with FLAG_ANISOTROPIC_VOXEL, otherwise the one
    // step on every axis (every stream written before 5.0.0).
    let voxel_size = if header.anisotropic_voxel() {
        need(offset + 8 + 8)?;
        let read_f32 = |o: usize| f32::from_bits(read_u32(o));
        let v = Vec3::new(header.voxel_size, read_f32(offset), read_f32(offset + 4));
        offset += 8;
        v
    } else {
        Vec3::splat(header.voxel_size)
    };

    // 2. Read quantizer params
    let step = read_u32(offset) as i32;
    offset += 4;
    let dead_zone = read_u32(offset) as i32;
    offset += 4;
    if step < 1 || dead_zone < 0 {
        return Err(DecodeError::InvalidQuantizer { step, dead_zone });
    }
    let quantizer = Quantizer::with_dead_zone(step, dead_zone);

    // 3. Read histogram if rANS compressed
    let use_rans = header.rans_compressed();

    let mut histogram = [0u32; 256];
    if use_rans {
        need(offset + 256 * 4)?;
        for h in &mut histogram {
            *h = read_u32(offset);
            offset += 4;
        }
    }

    // 4. Read payload
    need(offset + 4)?;
    let payload_len = read_u32(offset) as usize;
    offset += 4;
    need(offset.saturating_add(payload_len))?;
    let payload_bytes = &data[offset..offset + payload_len];

    // 5. Decode raw bytes (i16 LE, or i32 LE with FLAG_WIDE_COEFFS)
    let coeff_bytes = if header.wide_coeffs() { 4 } else { 2 };
    let raw_byte_count = total
        .checked_mul(coeff_bytes)
        .ok_or(DecodeError::InvalidDimensions {
            width: header.width,
            height: header.height,
            depth: header.depth,
            symbol_count: header.symbol_count,
        })?;
    if use_rans {
        check_rans_histogram(&histogram, raw_byte_count, payload_len)?;
    }
    let raw_bytes: Vec<u8> = if use_rans {
        let freq_table = FrequencyTable::from_histogram(&histogram);
        let mut decoder = RansDecoder::new(payload_bytes);
        decoder.decode_n(raw_byte_count, &freq_table)
    } else {
        if payload_len < raw_byte_count {
            return Err(DecodeError::Truncated {
                needed: offset + raw_byte_count,
                len: data.len(),
            });
        }
        payload_bytes[..raw_byte_count].to_vec()
    };

    // 6. Convert the LE bytes back to quantized i32 coefficients
    let quantized: Vec<i32> = if header.wide_coeffs() {
        raw_bytes
            .chunks_exact(4)
            .map(|c| i32::from_le_bytes([c[0], c[1], c[2], c[3]]))
            .collect()
    } else {
        raw_bytes
            .chunks_exact(2)
            .map(|chunk| i32::from(i16::from_le_bytes([chunk[0], chunk[1]])))
            .collect()
    };

    // 7. Dequantize. The magnitude `dz + (|q| - 1)·step + step/2` grows with
    // |q| (step >= 1), so checking the largest |q| keeps every value in i32.
    let max_q = quantized
        .iter()
        .map(|q| q.unsigned_abs())
        .max()
        .unwrap_or(0);
    if max_q > 0 {
        let dz = i64::from(dead_zone.max(1));
        let mag = dz + (i64::from(max_q) - 1) * i64::from(step) + i64::from(step / 2);
        if mag > i64::from(i32::MAX) {
            return Err(DecodeError::CoefficientOverflow);
        }
    }
    let mut coeffs = vec![0i32; total];
    // `InvalidBufferSize` is impossible: both buffers are `total` long.
    quantizer
        .dequantize_buffer(&quantized, &mut coeffs)
        .expect("coefficient buffer is sized to the quantized buffer");

    // 8. Inverse wavelet transform, with every step checked for i32 overflow
    let steps = if header.lossless_wavelet() {
        CDF53_STEPS
    } else {
        CDF97_STEPS
    };
    checked_inverse_3d(&mut coeffs, w, h, d, steps)?;

    // 9. Convert fixed-point i32 back to f32 distances
    let inv_scale = 1.0 / scale;
    let float_data: Vec<f32> = coeffs.iter().map(|&c| c as f32 * inv_scale).collect();

    Ok(SdfVolume {
        data: float_data,
        width: w,
        height: h,
        depth: d,
        origin: Vec3::new(header.origin[0], header.origin[1], header.origin[2]),
        voxel_size,
    })
}

/// The rANS histogram as the encoder writes it: one count per byte value of
/// the `bytes` coefficient bytes, used only when `bytes >= 512` and no value
/// holds half of them (otherwise the payload is stored raw)
///
/// With the largest count below half, the normalized frequency of every
/// symbol is at most `2048 + 256` of 4096, so decoding a symbol shrinks the
/// rANS state (at least 2^23 before every step) by at least
/// `log2(4096 / 2304) > 0.8` bit. The state starts with 32 bits, ends at
/// 2^23, and gains 8 bits per payload byte read, so a stream of `bytes`
/// symbols needs `bytes · 0.8 <= 8 · payload_len + 9` bits: `bytes <= 16 ·
/// payload_len` holds for every stream the encoder writes and bounds the
/// allocation by the input size.
fn check_rans_histogram(
    histogram: &[u32; 256],
    bytes: usize,
    payload_len: usize,
) -> Result<(), DecodeError> {
    let sum: u64 = histogram.iter().map(|&c| u64::from(c)).sum();
    let max = histogram.iter().copied().max().map_or(0, u64::from);
    let ok = sum == bytes as u64
        && bytes >= 512
        && 2 * max < sum
        && (bytes as u64) <= 16 * payload_len as u64;
    if ok {
        Ok(())
    } else {
        Err(DecodeError::InvalidHistogram)
    }
}

/// Lifting steps `(coeff, predict)` of the CDF 5/3 wavelet of alice-codec
/// (`Wavelet1D::cdf53`; coefficients in units of 1/4096)
const CDF53_STEPS: &[(i32, bool)] = &[(-2048, true), (1024, false)];
/// Lifting steps of the CDF 9/7 wavelet of alice-codec (`Wavelet1D::cdf97`)
const CDF97_STEPS: &[(i32, bool)] = &[(-6497, true), (-217, false), (3616, true), (1817, false)];

/// `Wavelet3D::inverse` of alice-codec with every addition checked
///
/// Same order (the third axis, then per slice the columns and the rows), the
/// same interleave and the same lifting arithmetic (`odd -= round(coeff ·
/// (even_l + even_r) / 4096)`, mirror boundary, round half up), so for every
/// input on which the reference does not overflow the result is bit-identical
/// (pinned by a unit test against the reference). Where the reference would
/// overflow `i32` (a panic with overflow checks, a wrapped value without
/// them) this returns [`DecodeError::CoefficientOverflow`].
fn checked_inverse_3d(
    volume: &mut [i32],
    width: usize,
    height: usize,
    depth: usize,
    steps: &[(i32, bool)],
) -> Result<(), DecodeError> {
    let frame = width * height;
    let mut line = Vec::new();
    let mut run = |idx: &mut dyn Iterator<Item = usize>, volume: &mut [i32]| {
        let idx: Vec<usize> = idx.collect();
        line.clear();
        line.extend(idx.iter().map(|&i| volume[i]));
        checked_inverse_1d(&mut line, steps)?;
        for (&i, &v) in idx.iter().zip(&line) {
            volume[i] = v;
        }
        Ok::<(), DecodeError>(())
    };
    for y in 0..height {
        for x in 0..width {
            run(&mut (0..depth).map(|t| t * frame + y * width + x), volume)?;
        }
    }
    for t in 0..depth {
        let base = t * frame;
        for x in 0..width {
            run(&mut (0..height).map(|y| base + y * width + x), volume)?;
        }
        for y in 0..height {
            run(&mut (0..width).map(|x| base + y * width + x), volume)?;
        }
    }
    Ok(())
}

/// One line of [`checked_inverse_3d`] (`Wavelet1D::inverse` of alice-codec)
fn checked_inverse_1d(signal: &mut [i32], steps: &[(i32, bool)]) -> Result<(), DecodeError> {
    let n = signal.len();
    if n < 2 {
        return Ok(());
    }
    let half = n / 2;
    let low = n - half;
    // interleave: [e0, e1, .., o0, o1, ..] -> [e0, o0, e1, o1, ..]
    let src = signal.to_vec();
    for i in 0..low {
        signal[i * 2] = src[i];
    }
    for i in 0..half {
        signal[i * 2 + 1] = src[low + i];
    }
    let fit = |v: i64| i32::try_from(v).map_err(|_| DecodeError::CoefficientOverflow);
    for &(coeff, predict) in steps.iter().rev() {
        for i in 0..half {
            // predict: odd[i] -= delta(even_l + even_r); update: even[i] -=
            // delta(odd_l + odd_r), the neighbour mirrored at the boundary
            let (target, a, b) = if predict {
                let l = i * 2;
                let r = if l + 2 < n { l + 2 } else { l };
                (i * 2 + 1, l, r)
            } else {
                let r = i * 2 + 1;
                let l = if i > 0 { i * 2 - 1 } else { r };
                (i * 2, l, r)
            };
            let sum = fit(i64::from(signal[a]) + i64::from(signal[b]))?;
            let delta = fit((i64::from(sum) * i64::from(coeff) + 2048) >> 12)?;
            signal[target] = fit(i64::from(signal[target]) - i64::from(delta))?;
        }
    }
    Ok(())
}

// ============================================================================
// Volume Statistics
// ============================================================================

/// Compute statistics about an SDF volume.
///
/// Useful for estimating compression characteristics before encoding.
pub fn volume_stats(volume: &SdfVolume) -> VolumeStats {
    let mut min = f32::MAX;
    let mut max = f32::MIN;
    let mut sum_sq: f64 = 0.0;

    for &v in &volume.data {
        if v < min {
            min = v;
        }
        if v > max {
            max = v;
        }
        sum_sq = (v as f64).mul_add(v as f64, sum_sq);
    }

    // Count zero-crossings along X axis (sign changes = surface voxels)
    let mut zero_crossing_count = 0usize;
    for z in 0..volume.depth {
        for y in 0..volume.height {
            for x in 1..volume.width {
                let idx = z * volume.height * volume.width + y * volume.width + x;
                if volume.data[idx] * volume.data[idx - 1] < 0.0 {
                    zero_crossing_count += 1;
                }
            }
        }
    }

    let n = volume.data.len();
    let rms = if n > 0 {
        (sum_sq / n as f64).sqrt()
    } else {
        0.0
    };

    VolumeStats {
        min_distance: min,
        max_distance: max,
        zero_crossings: zero_crossing_count,
        total_voxels: n,
        surface_ratio: if n > 0 {
            zero_crossing_count as f64 / n as f64
        } else {
            0.0
        },
        rms,
    }
}

/// Compute the compression ratio for an encoded SDF volume.
///
/// Returns `uncompressed_size / compressed_size`.
pub fn compression_ratio(volume: &SdfVolume, encoded: &[u8]) -> f64 {
    let uncompressed = volume.len() * std::mem::size_of::<f32>();
    uncompressed as f64 / encoded.len() as f64
}

// ============================================================================
// One-call Pipeline API
// ============================================================================

/// SDF 圧縮パイプライン結果。
#[derive(Clone, Debug)]
pub struct CompressResult {
    /// 圧縮されたバイトストリーム。
    pub encoded: Vec<u8>,
    /// ボリューム統計情報。
    pub stats: VolumeStats,
    /// 圧縮率（非圧縮サイズ / 圧縮サイズ）。
    pub ratio: f64,
    /// ボクセル化に使用した解像度。
    pub resolution: [usize; 3],
}

/// ワンコールで SDF ノードを圧縮する。
///
/// `SdfNode` → ボクセル化 → ウェーブレット変換 → 量子化 → rANS 符号化
/// の全パイプラインを 1 回の呼び出しで実行する。
///
/// # Arguments
/// * `sdf` - 圧縮対象の SDF ノード
/// * `origin` - ワールド空間の最小コーナー
/// * `extent` - ワールド空間の最大コーナー
/// * `resolution` - グリッド解像度（各軸）
/// * `config` - エンコード設定
pub fn compress_sdf(
    sdf: &SdfNode,
    origin: Vec3,
    extent: Vec3,
    resolution: usize,
    config: &EncodeConfig,
) -> CompressResult {
    let volume = voxelize_sdf_uniform(sdf, origin, extent, resolution);
    let stats = volume_stats(&volume);
    let encoded = encode_sdf_volume(&volume, config);
    let ratio = compression_ratio(&volume, &encoded);

    CompressResult {
        encoded,
        stats,
        ratio,
        resolution: [resolution, resolution, resolution],
    }
}

/// 圧縮された SDF を復元する。
///
/// `compress_sdf` で生成されたバイトストリームから `SdfVolume` を復元する。
pub fn decompress_sdf(data: &[u8]) -> SdfVolume {
    decode_sdf_volume(data)
}

// ============================================================================
// Tests
// ============================================================================

#[cfg(test)]
mod tests {
    use super::*;
    use crate::types::SdfNode;

    /// The checked inverse is bit-identical to alice-codec's
    /// `Wavelet3D::inverse` on inputs where the reference does not overflow
    /// (both wavelets, odd / even / unit lengths on every axis)
    #[test]
    fn checked_inverse_matches_the_reference_wavelet() {
        let mut state = 0x9e37_79b9_7f4a_7c15u64;
        let mut compared = 0;
        for (steps, reference) in [
            (CDF53_STEPS, Wavelet3D::cdf53()),
            (CDF97_STEPS, Wavelet3D::cdf97()),
        ] {
            for dims in [
                [1, 1, 1],
                [2, 1, 1],
                [1, 3, 1],
                [5, 4, 3],
                [8, 8, 8],
                [7, 2, 9],
            ] {
                let [w, h, d] = dims;
                let n = w * h * d;
                let data: Vec<i32> = (0..n)
                    .map(|_| {
                        state = state
                            .wrapping_mul(6_364_136_223_846_793_005)
                            .wrapping_add(1_442_695_040_888_963_407);
                        // +-2^20: far from overflow through 3 axes of lifting
                        ((state >> 43) as i32) - (1 << 20)
                    })
                    .collect();
                let mut want = data.clone();
                reference.inverse(&mut want, w, h, d);
                let mut got = data;
                checked_inverse_3d(&mut got, w, h, d, steps).unwrap();
                assert_eq!(got, want, "{dims:?}");
                compared += n;
            }
        }
        assert!(compared > 1000);
    }

    #[test]
    fn checked_inverse_reports_overflow() {
        let mut v = vec![i32::MAX, i32::MAX];
        assert_eq!(
            checked_inverse_3d(&mut v, 2, 1, 1, CDF53_STEPS),
            Err(DecodeError::CoefficientOverflow)
        );
    }

    #[test]
    fn test_voxelize_sphere() {
        let sphere = SdfNode::sphere(1.0);
        let volume = voxelize_sdf(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), [8, 8, 8]);

        assert_eq!(volume.len(), 512);
        assert_eq!(volume.width, 8);
        assert_eq!(volume.height, 8);
        assert_eq!(volume.depth, 8);

        // Center voxel (3,3,3) should be inside the sphere (negative distance)
        let center_d = volume.get(3, 3, 3);
        assert!(
            center_d < 0.0,
            "Center should be inside sphere, got {}",
            center_d
        );

        // Corner voxel (0,0,0) at (-2,-2,-2) should be outside
        let corner_d = volume.get(0, 0, 0);
        assert!(
            corner_d > 0.0,
            "Corner should be outside sphere, got {}",
            corner_d
        );
    }

    #[test]
    fn test_voxelize_uniform() {
        let box_sdf = SdfNode::box3d(1.0, 1.0, 1.0);
        let volume = voxelize_sdf_uniform(&box_sdf, Vec3::splat(-2.0), Vec3::splat(2.0), 4);
        assert_eq!(volume.len(), 64);
    }

    #[test]
    fn test_volume_stats() {
        let sphere = SdfNode::sphere(1.0);
        let volume = voxelize_sdf_uniform(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), 8);

        let stats = volume_stats(&volume);

        assert!(
            stats.min_distance < 0.0,
            "Sphere should have negative interior"
        );
        assert!(
            stats.max_distance > 0.0,
            "Sphere should have positive exterior"
        );
        assert!(
            stats.zero_crossings > 0,
            "Sphere should have surface crossings"
        );
        assert!(stats.surface_ratio > 0.0 && stats.surface_ratio < 1.0);
        assert!(stats.rms > 0.0);
        assert_eq!(stats.total_voxels, 512);
    }

    #[test]
    fn test_encode_decode_roundtrip() {
        // Use a small 8^3 sphere volume for fast testing
        let sphere = SdfNode::sphere(1.0);
        let volume = voxelize_sdf_uniform(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), 8);

        let config = EncodeConfig {
            quality: 99,
            fixed_point_scale: 4096.0,
            lossless_wavelet: true, // CDF 5/3 for better roundtrip
        };

        let encoded = encode_sdf_volume(&volume, &config);
        assert!(!encoded.is_empty(), "Encoded data should not be empty");

        let decoded = decode_sdf_volume(&encoded);

        // Verify dimensions match
        assert_eq!(decoded.width, volume.width);
        assert_eq!(decoded.height, volume.height);
        assert_eq!(decoded.depth, volume.depth);
        assert_eq!(decoded.data.len(), volume.data.len());

        // Verify origin preserved
        assert!((decoded.origin - volume.origin).length() < 1e-6);

        // Verify distances are approximately correct
        // (Quantization introduces some error)
        let mut max_error: f32 = 0.0;
        for (orig, recon) in volume.data.iter().zip(decoded.data.iter()) {
            let err = (orig - recon).abs();
            if err > max_error {
                max_error = err;
            }
        }

        // With quality 99 and scale 4096, max error should be very small
        assert!(
            max_error < 0.5,
            "Max reconstruction error {} is too large",
            max_error
        );
    }

    #[test]
    fn test_encode_preserves_sign() {
        // The sign of the SDF (inside/outside) must be preserved after encoding
        let sphere = SdfNode::sphere(1.0);
        let volume = voxelize_sdf_uniform(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), 8);

        let config = EncodeConfig::high_quality();
        let encoded = encode_sdf_volume(&volume, &config);
        let decoded = decode_sdf_volume(&encoded);

        // Count sign mismatches (excluding near-zero values where sign is ambiguous)
        let threshold = 0.1;
        let mut sign_mismatches = 0;
        let mut checked = 0;
        for (&orig, &recon) in volume.data.iter().zip(decoded.data.iter()) {
            if orig.abs() > threshold {
                checked += 1;
                if orig.signum() != recon.signum() {
                    sign_mismatches += 1;
                }
            }
        }

        let mismatch_rate = if checked > 0 {
            sign_mismatches as f64 / checked as f64
        } else {
            0.0
        };

        assert!(
            mismatch_rate < 0.05,
            "Sign mismatch rate {} is too high ({}/{} voxels)",
            mismatch_rate,
            sign_mismatches,
            checked
        );
    }

    #[test]
    fn test_compression_ratio_positive() {
        let sphere = SdfNode::sphere(1.0);
        let volume = voxelize_sdf_uniform(&sphere, Vec3::splat(-2.0), Vec3::splat(2.0), 8);

        let encoded = encode_sdf_volume(&volume, &EncodeConfig::default());
        let ratio = compression_ratio(&volume, &encoded);

        // Ratio should be positive (always is by definition)
        assert!(ratio > 0.0, "Compression ratio should be positive");
    }

    #[test]
    fn test_world_pos() {
        let volume = SdfVolume {
            data: vec![0.0; 8],
            width: 2,
            height: 2,
            depth: 2,
            origin: Vec3::new(1.0, 2.0, 3.0),
            voxel_size: Vec3::splat(0.5),
        };

        let p = volume.world_pos(1, 1, 1);
        assert!((p.x - 1.5).abs() < 1e-6);
        assert!((p.y - 2.5).abs() < 1e-6);
        assert!((p.z - 3.5).abs() < 1e-6);
    }

    #[test]
    fn test_compress_sdf_pipeline() {
        let sphere = SdfNode::sphere(1.0);
        let result = compress_sdf(
            &sphere,
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            8,
            &EncodeConfig::default(),
        );

        assert!(!result.encoded.is_empty());
        assert!(result.ratio > 0.0);
        assert_eq!(result.resolution, [8, 8, 8]);
        assert_eq!(result.stats.total_voxels, 512);
        assert!(result.stats.zero_crossings > 0);
    }

    #[test]
    fn test_compress_decompress_roundtrip() {
        let sphere = SdfNode::sphere(1.0);
        let result = compress_sdf(
            &sphere,
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            8,
            &EncodeConfig::high_quality(),
        );

        let volume = decompress_sdf(&result.encoded);
        assert_eq!(volume.width, 8);
        assert_eq!(volume.height, 8);
        assert_eq!(volume.depth, 8);
        assert_eq!(volume.data.len(), 512);
    }

    #[test]
    fn test_compress_result_fields() {
        let sphere = SdfNode::sphere(1.0);
        let result = compress_sdf(
            &sphere,
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            4,
            &EncodeConfig::fast(),
        );

        assert_eq!(result.resolution, [4, 4, 4]);
        assert_eq!(result.stats.total_voxels, 64);
    }

    #[test]
    fn test_encode_config_presets() {
        let fast = EncodeConfig::fast();
        assert_eq!(fast.quality, 50);

        let hq = EncodeConfig::high_quality();
        assert_eq!(hq.quality, 95);

        let lossless = EncodeConfig::lossless();
        assert!(lossless.lossless_wavelet);
        assert_eq!(lossless.quality, 100);
    }
}