zenjpeg 0.8.2

//! Adaptive Quantization — C++ parity implementation.
//!
//! Full port of jpegli's adaptive quantization pipeline from
//! `lib/jpegli/adaptive_quantization.cc`. Produces per-block `aq_strength`
//! values (typically 0.0–0.2, mean ~0.08) used to modulate zero-bias
//! thresholds during coefficient quantization.
//!
//! ## Algorithm Steps
//!
//! 1. **ComputePreErosion()** — Initial quant field from pixel statistics.
//!    Uses `QuantMasking()` for spatial frequency masking and `MaskingSqrt()`
//!    to combine with butteraugli distance.
//!
//! 2. **FuzzyErosion()** — Spatial smoothing of the quant field.
//!    5×5 kernel with asymmetric weights, separate horizontal/vertical passes.
//!
//! 3. **PerBlockModulations()** — Final per-block modulations using spatial
//!    frequency analysis and AC energy distribution masking.
//!
//! 4. **Final Transform** — Converts quant_field to aq_strength:
//!    ```text
//!    aq_strength = max(0.0, (0.6 / quant_field) - 1.0)
//!    ```
//!
//! All four steps have SIMD implementations (AVX2/AVX-512) and a streaming
//! variant for strip-based encoding. See [`streaming::StreamingAQ`].
//!
//! See also:
//! - `docs/ADAPTIVE_QUANTIZATION.md` for detailed analysis
//! - `tests/aq_locked_tests.rs` for invariant tests

#![allow(dead_code)]

use std::f32::consts::PI;

use crate::foundation::aligned_alloc::{AlignedVec, AllocError, try_alloc_zeroed};

// SIMD implementations for hot paths
pub mod simd;

// Auto-vectorized implementations using archmage autoversion
pub mod autovec;

// Streaming AQ for strip-based encoding
pub mod streaming;

use simd::{compute_pre_erosion_simd, fuzzy_erosion_simd, per_block_modulations_simd};

// ============================================================================
// Constants ported from C++ adaptive_quantization.cc
// ============================================================================

/// Gamma-related constant from ComputePreErosion.
/// Note: C++ divides by kInputScaling (255.0), applied at usage.
const MATCH_GAMMA_OFFSET: f32 = 0.019;

/// Limit threshold for pre-erosion.
#[allow(dead_code)] // Kept for reference - used in C++ but handled differently in Rust
const LIMIT: f32 = 0.2;

/// AC quantization scaling constant.
#[allow(dead_code)] // Kept for reference - used in C++ but handled differently in Rust
const K_AC_QUANT: f32 = 0.841;

/// Input scaling factor (1/255 for 8-bit input).
const K_INPUT_SCALING: f32 = 1.0 / 255.0;

/// Gamma modulation bias (adjusted for scaling).
const K_GAMMA_MOD_BIAS: f32 = 0.16 * K_INPUT_SCALING;

/// Gamma modulation scale.
const K_GAMMA_MOD_SCALE: f32 = 1.0 / 64.0;

/// Inverse of ln(2) = ln(2).
const K_INV_LOG2E: f32 = 0.6931471805599453;

/// Gamma modulation gamma coefficient.
const K_GAMMA_MOD_GAMMA: f32 = -0.15526878023684174 * K_INV_LOG2E;

/// High-frequency modulation coefficient.
const K_HF_MOD_COEFF: f32 = -2.0052193233688884 / 112.0;

// Constants for ComputeMask (from C++)
// CRITICAL: kBase is NEGATIVE in C++ (-0.74174993)!
const K_MASK_BASE: f32 = -0.74174993;
#[allow(dead_code)] // Old constants kept for reference - replaced by K_MASK_* versions
const K_MUL4: f32 = 0.03879999369382858;
#[allow(dead_code)] // Old constants kept for reference - replaced by K_MASK_* versions
const K_MUL2: f32 = 0.17580001056194305;
const K_MASK_MUL4: f32 = 3.2353257320940401;
const K_MASK_MUL2: f32 = 12.906028311180409;
const K_MASK_OFFSET2: f32 = 305.04035728311436;
const K_MASK_MUL3: f32 = 5.0220313103171232;
#[allow(dead_code)] // Old constant kept for reference - replaced by K_MASK_* version
const K_MUL3: f32 = 0.30230000615119934;
const K_MASK_OFFSET3: f32 = 2.1925739705298404;
const K_MASK_OFFSET4: f32 = 0.25 * K_MASK_OFFSET3;
const K_MASK_MUL0: f32 = 0.74760422233706747;

// Constants from RatioOfDerivatives
const K_EPSILON_RATIO: f32 = 1e-2;
const K_NUM_OFFSET_RATIO: f32 = K_EPSILON_RATIO / K_INPUT_SCALING / K_INPUT_SCALING;
const K_SG_MUL: f32 = 226.0480446705883;
const K_SG_MUL2: f32 = 1.0 / 73.377132366608819;
const K_SG_RET_MUL: f32 = K_SG_MUL2 * 18.6580932135 * K_INV_LOG2E;
const K_NUM_MUL_RATIO: f32 = K_SG_RET_MUL * 3.0 * K_SG_MUL;
const K_SG_VOFFSET: f32 = 7.14672470003;
const K_VOFFSET_RATIO: f32 = (K_SG_VOFFSET * K_INV_LOG2E + K_EPSILON_RATIO) / K_INPUT_SCALING;
const K_DEN_MUL_RATIO: f32 = K_INV_LOG2E * K_SG_MUL * K_INPUT_SCALING * K_INPUT_SCALING;

// ============================================================================
// Public API
// ============================================================================

/// Per-block adaptive quantization strength.
///
/// Values are in the range 0.0-0.2 (matching C++ output).
#[derive(Debug, Clone)]
pub struct AQStrengthMap {
    /// Width in 8x8 blocks
    pub width_blocks: usize,
    /// Height in 8x8 blocks
    pub height_blocks: usize,
    /// Per-block aq_strength values (0.0 to ~0.2)
    pub strengths: Vec<f32>,
}

impl AQStrengthMap {
    /// Creates a uniform AQ map with the given constant strength.
    ///
    /// The default C++ mean is ~0.08.
    #[must_use]
    pub fn uniform(width_blocks: usize, height_blocks: usize, strength: f32) -> Self {
        Self {
            width_blocks,
            height_blocks,
            strengths: vec![strength; width_blocks * height_blocks],
        }
    }

    /// Creates a uniform map with the C++ testdata mean (0.08).
    #[must_use]
    pub fn with_cpp_mean(width_blocks: usize, height_blocks: usize) -> Self {
        Self::uniform(width_blocks, height_blocks, 0.08)
    }

    /// Returns the aq_strength for a block.
    #[inline]
    #[must_use]
    pub fn get(&self, bx: usize, by: usize) -> f32 {
        let idx = by * self.width_blocks + bx;
        self.strengths.get(idx).copied().unwrap_or(0.08)
    }

    /// Returns the mean AQ strength across all blocks.
    ///
    /// This value can be used with [`crate::encode::trellis::should_use_hybrid`]
    /// to decide whether hybrid trellis will benefit this image.
    #[must_use]
    pub fn mean(&self) -> f32 {
        if self.strengths.is_empty() {
            return 0.0;
        }
        self.strengths.iter().sum::<f32>() / self.strengths.len() as f32
    }

    /// Returns the standard deviation of AQ strength.
    #[must_use]
    pub fn std(&self) -> f32 {
        if self.strengths.is_empty() {
            return 0.0;
        }
        let mean = self.mean();
        let variance = self
            .strengths
            .iter()
            .map(|&x| {
                let d = x - mean;
                d * d
            })
            .sum::<f32>()
            / self.strengths.len() as f32;
        variance.sqrt()
    }

    /// Returns (min, max, mean, std) statistics for the AQ map.
    ///
    /// Computes all statistics in a single pass over the data.
    #[must_use]
    pub fn stats(&self) -> (f32, f32, f32, f32) {
        if self.strengths.is_empty() {
            return (0.0, 0.0, 0.0, 0.0);
        }
        let mut min = f32::INFINITY;
        let mut max = f32::NEG_INFINITY;
        let mut sum = 0.0f32;
        for &x in &self.strengths {
            if x < min {
                min = x;
            }
            if x > max {
                max = x;
            }
            sum += x;
        }
        let n = self.strengths.len() as f32;
        let mean = sum / n;
        let mut var_sum = 0.0f32;
        for &x in &self.strengths {
            let d = x - mean;
            var_sum += d * d;
        }
        let std = (var_sum / n).sqrt();
        (min, max, mean, std)
    }

    /// Scale all AQ strengths by a factor.
    ///
    /// - `scale > 1.0` → more bits to complex regions → LARGER files
    /// - `scale < 1.0` → fewer bits to complex regions → SMALLER files
    ///
    /// This can be used to compensate for hybrid trellis size increase
    /// without changing the base quality level.
    pub fn scale(&mut self, factor: f32) {
        for s in &mut self.strengths {
            *s *= factor;
        }
    }

    /// Scale to achieve a target mean AQ strength.
    ///
    /// Useful for normalizing AQ maps to a consistent aggressiveness level.
    pub fn scale_to_mean(&mut self, target_mean: f32) {
        let current_mean = self.mean();
        if current_mean > 0.0 {
            self.scale(target_mean / current_mean);
        }
    }

    /// Estimate scale factor needed to reduce file size by a given percentage.
    ///
    /// Based on empirical model: higher AQ strength → more bits to complex
    /// regions → LARGER files. To reduce size, we need to DECREASE AQ strength.
    ///
    /// Model: ~20% size change per 0.1 AQ mean change (empirical from Kodak corpus).
    /// Note: There's also a baseline overhead (~3-10%) from trellis that can't be
    /// eliminated by AQ scaling alone.
    ///
    /// # Arguments
    /// * `size_reduction_pct` - Target size reduction (e.g., 15.0 for 15%)
    ///
    /// # Returns
    /// Scale factor to apply to the AQ map (will be < 1.0 for size reduction)
    #[must_use]
    pub fn scale_for_size_reduction(&self, size_reduction_pct: f32) -> f32 {
        // Empirical: ~20% size increase per 0.1 AQ mean increase
        // To REDUCE by X%, need to DECREASE mean by X/200
        // (i.e., to reduce by 10%, decrease mean by 0.05)
        let current_mean = self.mean();
        if current_mean <= 0.0 {
            return 1.0;
        }
        let mean_decrease_needed = size_reduction_pct / 200.0;
        let target_mean = (current_mean - mean_decrease_needed).max(0.01);
        target_mean / current_mean
    }
}

/// Computes per-block adaptive quantization strength.
///
/// # Status: RETURNS CONSTANT (C++ mean)
///
/// This function currently returns a uniform map with the C++ testdata mean.
/// The algorithm below is implemented but needs verification against C++ testdata.
///
/// # Arguments
///
/// * `y_plane` - Luminance plane (Y channel) as f32 values (0-255 range)
/// * `width` - Image width in pixels
/// * `height` - Image height in pixels
/// * `distance` - Quality distance parameter
///
/// # Returns
///
/// Per-block aq_strength map with values in 0.0-0.2 range.
///
/// # Arguments
/// * `y_plane` - Luminance plane (0-255 range)
/// * `width` - Image width in pixels
/// * `height` - Image height in pixels
/// * `y_quant_01` - Y quantization table value at position 1 (first AC coefficient)
///                  This is what C++ uses for the dampen calculation.
///
/// # Errors
///
/// Returns `AllocError` if internal buffer allocation fails.
pub fn compute_aq_strength_map(
    y_plane: &[f32],
    width: usize,
    height: usize,
    y_quant_01: u16,
) -> Result<AQStrengthMap, AllocError> {
    // Use per-block implementation
    // y_quant_01 is the actual quant table value, matching C++ behavior
    compute_aq_strength_map_impl(y_plane, width, height, y_quant_01 as f32)
}

/// Converts quant_field to aq_strength.
///
/// C++ formula: `aq_strength = max(0.0, (0.6 / quant_field) - 1.0)`
///
/// This is the CRITICAL missing piece from the failed port.
#[inline]
#[must_use]
pub fn quant_field_to_aq_strength(quant_field: f32) -> f32 {
    (0.6 / quant_field - 1.0).max(0.0)
}

/// SIMD-optimized quant_field to aq_strength conversion.
///
/// # Errors
///
/// Returns `AllocError` if result buffer allocation fails.
pub fn quant_field_to_aq_strength_simd(quant_field: &[f32]) -> Result<AlignedVec<f32>, AllocError> {
    let mut result = try_alloc_zeroed(quant_field.len())?;
    for i in 0..quant_field.len() {
        result[i] = quant_field_to_aq_strength(quant_field[i]);
    }
    Ok(result)
}

// ============================================================================
// Implementation (needs verification against C++ testdata)
// ============================================================================

/// Full implementation of AQ strength computation.
///
/// # Status: UNVERIFIED - DO NOT USE IN PRODUCTION
///
/// This implementation needs to be verified against C++ testdata before use.
/// Exposed as pub for testing against C++ testdata.
///
/// # Errors
///
/// Returns `AllocError` if any internal buffer allocation fails.
pub fn compute_aq_strength_map_impl(
    y_plane: &[f32],
    width: usize,
    height: usize,
    y_quant_01: f32,
) -> Result<AQStrengthMap, AllocError> {
    let width_blocks = (width + 7) / 8;
    let height_blocks = (height + 7) / 8;
    let num_blocks = width_blocks * height_blocks;

    if width == 0 || height == 0 {
        return Ok(AQStrengthMap::uniform(0, 0, 0.08));
    }

    // 1. ComputePreErosion (downsamples 4x) - SIMD accelerated
    // NOTE: y_plane is in 0-255 range - ratio_of_derivatives expects this range
    // (its constants have kInputScaling = 1/255 baked in)
    let pre_erosion = compute_pre_erosion_simd(y_plane, width, height)?;

    // 2. FuzzyErosion - SIMD accelerated
    let pre_erosion_w = (width + 3) / 4;
    let pre_erosion_h = (height + 3) / 4;
    let mut quant_field = try_alloc_zeroed(num_blocks)?;
    fuzzy_erosion_simd(
        &pre_erosion,
        pre_erosion_w,
        pre_erosion_h,
        width_blocks,
        height_blocks,
        &mut quant_field,
    )?;

    // 3. PerBlockModulations - SIMD accelerated
    // C++ uses y_quant_01 = quant_table[1] (first AC coefficient)
    // This is passed in directly from the quant table, NOT computed from distance
    per_block_modulations_simd(
        y_quant_01,
        y_plane, // Use original unscaled input, not input_scaled
        width,
        height,
        width_blocks,
        height_blocks,
        &mut quant_field,
    );

    // 4. Final transform: quant_field -> aq_strength (SIMD optimized)
    let strengths = quant_field_to_aq_strength_simd(&quant_field)?;

    Ok(AQStrengthMap {
        width_blocks,
        height_blocks,
        strengths: strengths.to_vec(),
    })
}

// ============================================================================
// Helper functions (ported from previous attempt)
// ============================================================================

/// Generates a 1D Gaussian kernel.
#[allow(dead_code)]
fn gaussian_kernel(sigma: f32, radius: usize) -> Vec<f32> {
    let mut kernel = vec![0.0; 2 * radius + 1];
    let sigma_sq = sigma * sigma;
    let norm_factor = 1.0 / (2.0 * PI * sigma_sq).sqrt();
    let mut sum = 0.0;

    for i in 0..=radius {
        let dist_sq = (i * i) as f32;
        let val = norm_factor * (-dist_sq / (2.0 * sigma_sq)).exp();
        kernel[radius + i] = val;
        kernel[radius - i] = val;
        sum += if i == 0 { val } else { 2.0 * val };
    }

    // Normalize
    if sum > 1e-6 {
        for val in &mut kernel {
            *val /= sum;
        }
    }

    kernel
}

/// Calculates the ratio of derivatives for psychovisual modulation.
/// Ported from `RatioOfDerivativesOfCubicRootToSimpleGamma`.
#[allow(dead_code)]
fn ratio_of_derivatives(val: f32, invert: bool) -> f32 {
    let v = val.max(0.0);
    let v2 = v * v;

    let num = v2.mul_add(K_NUM_MUL_RATIO, K_NUM_OFFSET_RATIO);
    let den = (v * K_DEN_MUL_RATIO).mul_add(v2, K_VOFFSET_RATIO);

    let safe_den = if den == 0.0 { 1e-9 } else { den };

    if invert {
        num / safe_den
    } else {
        safe_den / num
    }
}

/// Updates the 4 minimum values in a sorted array.
#[inline]
#[allow(dead_code)]
fn update_min4(val: f32, mins: &mut [f32; 4]) {
    if val < mins[3] {
        if val < mins[2] {
            mins[3] = mins[2];
            if val < mins[1] {
                mins[2] = mins[1];
                if val < mins[0] {
                    mins[1] = mins[0];
                    mins[0] = val;
                } else {
                    mins[1] = val;
                }
            } else {
                mins[2] = val;
            }
        } else {
            mins[3] = val;
        }
    }
}

/// MaskingSqrt from C++.
/// `return 0.25 * sqrt(v * sqrt(211.50759899638012e8) + 28)`
#[inline]
fn masking_sqrt(v: f32) -> f32 {
    const K_LOG_OFFSET: f32 = 28.0;
    const K_MUL: f32 = 211.50759899638012;
    0.25 * v.mul_add((K_MUL * 1e8_f32).sqrt(), K_LOG_OFFSET).sqrt()
}

/// FastPow2f from C++ jpegli - bit manipulation approximation of 2^x.
/// Based on Highway FastPow2f implementation.
/// Note: This doesn't match C++ exactly due to different polynomial coefficients.
/// Kept for potential future use/optimization.
#[allow(dead_code)]
#[inline]
fn fast_pow2f(x: f32) -> f32 {
    const LN2: f32 = 0.6931471805599453;

    // Clamp to avoid overflow/underflow
    let x = x.clamp(-126.0, 126.0);

    // Split into integer and fractional parts
    let xi = x.floor();
    let xf = x - xi;

    // Polynomial approximation for 2^xf (minimax polynomial for [0, 1])
    let p = 1.0
        + xf * (LN2
            + xf * (0.240226507 + xf * (0.0558325133 + xf * (0.00898934009 + xf * 0.00187757667))));

    // Multiply by 2^xi using bit manipulation
    let bits = p.to_bits();
    let exp_offset = ((xi as i32) << 23) as u32;
    f32::from_bits(bits.wrapping_add(exp_offset))
}

/// Ported from ComputePreErosion (scalar version).
///
/// C++ algorithm:
/// 1. For each pixel: base = 0.25 * (left + right + top + bottom)
/// 2. diff = ratio_of_derivatives(pixel + gamma_offset) * (pixel - base)
/// 3. diff_squared = diff * diff, clamped to 0.2
/// 4. masked = MaskingSqrt(diff_squared)
/// 5. Sum over 4x4 blocks to produce 4x downsampled output
///
/// NOTE: input must be in 0-255 range (not 0-1 scaled) because
/// ratio_of_derivatives has kInputScaling=1/255 baked into its constants.
#[allow(dead_code)]
fn compute_pre_erosion_scalar(input: &[f32], width: usize, height: usize) -> Vec<f32> {
    let pre_erosion_w = (width + 3) / 4;
    let pre_erosion_h = (height + 3) / 4;
    let mut pre_erosion = vec![0.0f32; pre_erosion_w * pre_erosion_h];

    const LIMIT: f32 = 0.2;
    // gamma_offset is in 0-255 range (0.019 * 255 = 4.85)
    let gamma_offset = MATCH_GAMMA_OFFSET / K_INPUT_SCALING;

    // Helper to get pixel with clamped bounds (returns 0-255 range)
    let get = |x: isize, y: isize| -> f32 {
        let x = x.clamp(0, width as isize - 1) as usize;
        let y = y.clamp(0, height as isize - 1) as usize;
        input[y * width + x]
    };

    // Temporary buffer for diff values at full resolution
    let mut diff_buffer = vec![0.0f32; width];

    for y_block in 0..pre_erosion_h {
        // Accumulate over 4 rows
        diff_buffer.fill(0.0);

        for iy in 0..4 {
            let y = y_block * 4 + iy;
            if y >= height {
                continue;
            }

            for x in 0..width {
                let ix = x as isize;
                let iy_s = y as isize;

                let pixel = get(ix, iy_s);

                // Base is average of 4 neighbors
                let base = 0.25
                    * (get(ix - 1, iy_s)
                        + get(ix + 1, iy_s)
                        + get(ix, iy_s - 1)
                        + get(ix, iy_s + 1));

                // Gamma-corrected ratio as weight
                let gammacv = ratio_of_derivatives(pixel + gamma_offset, false);

                // Weighted difference
                let diff = gammacv * (pixel - base);

                // Square and clamp
                let diff_sq = (diff * diff).min(LIMIT);

                // Apply MaskingSqrt
                let masked = masking_sqrt(diff_sq);

                // Accumulate
                diff_buffer[x] += masked;
            }
        }

        // Downsample 4x in x direction
        // C++ multiplies by 0.25 after summing 4 values
        for x_block in 0..pre_erosion_w {
            let x_start = x_block * 4;
            let mut sum = 0.0f32;
            for ix in 0..4 {
                let x = x_start + ix;
                if x < width {
                    sum += diff_buffer[x];
                }
            }
            pre_erosion[y_block * pre_erosion_w + x_block] = sum * 0.25;
        }
    }

    pre_erosion
}

/// Ported from FuzzyErosion (scalar version).
///
/// The C++ algorithm:
/// 1. For each pixel in pre_erosion, compute weighted sum of 4 smallest values in 3x3 window
/// 2. Weights: 0.125 * min0 + 0.075 * min1 + 0.06 * min2 + 0.05 * min3
/// 3. Sum 2x2 blocks to get final quant_field values
#[allow(dead_code)]
fn fuzzy_erosion_scalar(
    pre_erosion: &[f32],
    pre_erosion_w: usize,
    pre_erosion_h: usize,
    block_w: usize,
    block_h: usize,
    aq_map: &mut [f32],
) {
    assert_eq!(aq_map.len(), block_w * block_h);

    // Weights from C++ (sum = 0.31)
    const MUL0: f32 = 0.125;
    const MUL1: f32 = 0.075;
    const MUL2: f32 = 0.06;
    const MUL3: f32 = 0.05;

    // Temporary buffer for weighted min values
    let mut tmp = vec![0.0f32; pre_erosion_w * pre_erosion_h];

    // Helper to get value with clamped bounds
    let get = |x: isize, y: isize| -> f32 {
        let x = x.clamp(0, pre_erosion_w as isize - 1) as usize;
        let y = y.clamp(0, pre_erosion_h as isize - 1) as usize;
        pre_erosion[y * pre_erosion_w + x]
    };

    // Process each pixel - find 4 smallest in 3x3 window, compute weighted sum
    for y in 0..pre_erosion_h {
        for x in 0..pre_erosion_w {
            let ix = x as isize;
            let iy = y as isize;

            // Collect all 9 values in 3x3 window
            let mut vals = [
                get(ix - 1, iy - 1),
                get(ix, iy - 1),
                get(ix + 1, iy - 1),
                get(ix - 1, iy),
                get(ix, iy),
                get(ix + 1, iy),
                get(ix - 1, iy + 1),
                get(ix, iy + 1),
                get(ix + 1, iy + 1),
            ];

            // Partial sort to get 4 smallest values
            // (We only need the 4 smallest, not full sort)
            for i in 0..4 {
                for j in (i + 1)..9 {
                    if vals[j] < vals[i] {
                        vals.swap(i, j);
                    }
                }
            }

            // Weighted sum of 4 smallest (using FMA for consistency with SIMD path)
            let weighted = MUL0.mul_add(
                vals[0],
                MUL1.mul_add(vals[1], MUL2.mul_add(vals[2], MUL3 * vals[3])),
            );
            tmp[y * pre_erosion_w + x] = weighted;
        }
    }

    // Sum 2x2 blocks from tmp to get final aq_map values
    // C++ code (lines 476-481):
    //   aq_out[bx] = (row_out[x] + row_out[x + 1] + row_out0[x] + row_out0[x + 1]);
    // where x = bx * 2 and row_out0 = tmp.Row(y-1), row_out = tmp.Row(y)
    //
    // pre_erosion is at 4x4 pixel resolution (width/4, height/4)
    // blocks are at 8x8 pixel resolution (width/8, height/8)
    // So pre_erosion_w ≈ 2 * block_w, and we sum 2x2 pre_erosion cells per block

    // Helper to get tmp value with clamped bounds
    let get_tmp = |x: usize, y: usize| -> f32 {
        let x = x.min(pre_erosion_w.saturating_sub(1));
        let y = y.min(pre_erosion_h.saturating_sub(1));
        tmp[y * pre_erosion_w + x]
    };

    for by in 0..block_h {
        for bx in 0..block_w {
            // Each block sums 2x2 from tmp at (bx*2, by*2)
            let px = bx * 2;
            let py = by * 2;

            // Sum 2x2 values (NO multiplication factor - just sum)
            let sum = get_tmp(px, py)
                + get_tmp(px + 1, py)
                + get_tmp(px, py + 1)
                + get_tmp(px + 1, py + 1);

            aq_map[by * block_w + bx] = sum;
        }
    }
}

/// Ported from ComputeMask (scalar version).
#[allow(dead_code)]
fn compute_mask_scalar(out_val: f32) -> f32 {
    let v1 = (out_val * K_MASK_MUL0).max(1e-3);
    let v2 = 1.0 / (v1 + K_MASK_OFFSET2);
    let v3 = 1.0 / (v1 * v1 + K_MASK_OFFSET3);
    let v4 = 1.0 / (v1 * v1 + K_MASK_OFFSET4);
    // Use K_MASK_MUL* constants (3.24, 12.9, 5.02), NOT K_MUL* (0.04, 0.18, 0.30)
    // FMA chain for consistency with SIMD path
    v4.mul_add(
        K_MASK_MUL4,
        v2.mul_add(K_MASK_MUL2, v3.mul_add(K_MASK_MUL3, K_MASK_BASE)),
    )
}

/// Ported from HFModulation (scalar version).
#[allow(dead_code)]
fn hf_modulation_scalar(
    x: usize,
    y: usize,
    input_scaled: &[f32],
    width: usize,
    height: usize,
    current_val: f32,
) -> f32 {
    let center_idx = y * width + x;
    let center_val = input_scaled[center_idx];

    let left_idx = y * width + x.saturating_sub(1);
    let right_idx = y * width + (x + 1).min(width - 1);
    let top_idx = y.saturating_sub(1) * width + x;
    let bottom_idx = (y + 1).min(height - 1) * width + x;

    let diff_h =
        (input_scaled[left_idx] - center_val).abs() + (input_scaled[right_idx] - center_val).abs();
    let diff_v =
        (input_scaled[top_idx] - center_val).abs() + (input_scaled[bottom_idx] - center_val).abs();

    let diff_sum = diff_h + diff_v;
    diff_sum.mul_add(K_HF_MOD_COEFF, current_val)
}

/// Ported from GammaModulation (scalar version).
#[allow(dead_code)]
fn gamma_modulation_scalar(
    x: usize,
    y: usize,
    input_scaled: &[f32],
    width: usize,
    _height: usize,
    current_val: f32,
) -> f32 {
    let val = input_scaled[y * width + x];
    let log_arg = val.mul_add(K_GAMMA_MOD_SCALE, K_GAMMA_MOD_BIAS).max(1e-9);
    let modulation = K_GAMMA_MOD_GAMMA * super::ln_poly(log_arg);
    current_val + modulation
}

/// Ported from PerBlockModulations (scalar version).
///
/// C++ algorithm:
/// 1. ComputeMask on fuzzy erosion value
/// 2. HfModulation: sum of abs diffs for 8x8 block * kSumCoeff + out_val
/// 3. GammaModulation: kGamma * log2(sum_ratio * scale) + out_val
/// 4. Final: 2^(out_val * 1.442695) * mul + add
#[allow(dead_code)]
fn per_block_modulations_scalar(
    y_quant_01: f32,
    input_scaled: &[f32],
    width: usize,
    height: usize,
    block_w: usize,
    block_h: usize,
    aq_map: &mut [f32],
) {
    const K_AC_QUANT: f32 = 0.841;
    const K_DAMPEN_RAMP_START: f32 = 9.0;
    const K_DAMPEN_RAMP_END: f32 = 65.0;

    let base_level = 0.48 * K_AC_QUANT;

    // Compute dampen based on y_quant_01
    let dampen = if y_quant_01 >= K_DAMPEN_RAMP_START {
        let d =
            1.0 - (y_quant_01 - K_DAMPEN_RAMP_START) / (K_DAMPEN_RAMP_END - K_DAMPEN_RAMP_START);
        d.max(0.0)
    } else {
        1.0
    };

    let mul = K_AC_QUANT * dampen;
    let add = (1.0 - dampen) * base_level;

    // Helper to get pixel with clamped bounds
    let get = |x: usize, y: usize| -> f32 {
        let x = x.min(width.saturating_sub(1));
        let y = y.min(height.saturating_sub(1));
        input_scaled[y * width + x]
    };

    for by in 0..block_h {
        let y_start = by * 8;
        for bx in 0..block_w {
            let x_start = bx * 8;
            let block_idx = by * block_w + bx;

            // Get value from fuzzy erosion
            let fuzzy_val = aq_map[block_idx];

            // 1. ComputeMask
            let mut out_val = compute_mask_scalar(fuzzy_val);

            // 2. HfModulation: sum of absolute differences with right and below neighbors
            // kSumCoeff = -2.0052193233688884 * kInputScaling / 112.0
            const K_SUM_COEFF: f32 = -2.0052193233688884 * K_INPUT_SCALING / 112.0;
            let mut sum = 0.0f32;
            for dy in 0..8 {
                let y = y_start + dy;
                for dx in 0..8 {
                    let x = x_start + dx;
                    let p = get(x, y);

                    // Right neighbor (skip for rightmost column in block)
                    if dx < 7 {
                        let pr = get(x + 1, y);
                        sum += (p - pr).abs();
                    }

                    // Below neighbor (skip for bottom row)
                    if dy < 7 {
                        let pd = get(x, y + 1);
                        sum += (p - pd).abs();
                    }
                }
            }
            out_val += sum * K_SUM_COEFF;

            // 3. GammaModulation: sum ratio_of_derivatives (inverted) for 8x8 block
            // kGamma = -0.15526878023684174 * kInvLog2e
            const K_GAMMA: f32 = -0.15526878023684174 * K_INV_LOG2E;
            const K_BIAS: f32 = 0.16 / K_INPUT_SCALING;
            const K_SCALE: f32 = K_INPUT_SCALING / 64.0;

            let mut overall_ratio = 0.0f32;
            for dy in 0..8 {
                let y = y_start + dy;
                for dx in 0..8 {
                    let x = x_start + dx;
                    let iny = get(x, y) + K_BIAS;
                    let ratio_g = ratio_of_derivatives(iny, true);
                    overall_ratio += ratio_g;
                }
            }
            overall_ratio *= K_SCALE;

            // C++ uses FastLog2f which has ~3e-7 relative error
            // For parity, use std log2 (difference is negligible at ~0.0006% of AQ threshold)
            let log_ratio = if overall_ratio > 0.0 {
                overall_ratio.log2()
            } else {
                0.0
            };
            out_val += K_GAMMA * log_ratio;

            // 4. Final transform: 2^(out_val * 1.442695041) * mul + add
            // C++ uses FastPow2f, but standard exp2 produces equivalent results
            // 1.442695041 = 1/ln(2) = log2(e)
            const LOG2_E: f32 = 1.442695041;
            let quant_field = (out_val * LOG2_E).exp2().mul_add(mul, add);

            aq_map[block_idx] = quant_field;
        }
    }
}

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

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

    #[test]
    fn test_aq_strength_map_uniform() {
        let map = AQStrengthMap::uniform(10, 10, 0.08);
        assert_eq!(map.strengths.len(), 100);
        assert!((map.get(5, 5) - 0.08).abs() < 1e-6);
    }

    #[test]
    fn test_aq_strength_map_cpp_mean() {
        let map = AQStrengthMap::with_cpp_mean(10, 10);
        assert!((map.get(0, 0) - 0.08).abs() < 1e-6);
    }

    #[test]
    fn test_quant_field_to_aq_strength() {
        // quant_field = 0.6 → aq_strength = 0.0
        assert!((quant_field_to_aq_strength(0.6) - 0.0).abs() < 1e-6);

        // quant_field = 0.3 → aq_strength = 1.0
        assert!((quant_field_to_aq_strength(0.3) - 1.0).abs() < 1e-6);

        // quant_field = 6.0 → aq_strength = 0.0 (clamped)
        assert!(quant_field_to_aq_strength(6.0) >= 0.0);
    }

    #[test]
    fn test_compute_returns_uniform() {
        let plane = vec![128.0f32; 64 * 64];
        let map = compute_aq_strength_map(&plane, 64, 64, 1).unwrap();
        assert_eq!(map.width_blocks, 8);
        assert_eq!(map.height_blocks, 8);

        // For uniform input, all output values should be identical
        let first = map.strengths[0];
        for &s in &map.strengths {
            assert!(
                (s - first).abs() < 1e-6,
                "Uniform input should produce uniform output, got {} vs {}",
                s,
                first
            );
        }

        // The value should be positive and reasonable (0 < aq_strength < 1)
        assert!(first > 0.0, "aq_strength should be positive, got {}", first);
        assert!(first < 1.0, "aq_strength should be < 1, got {}", first);

        // NOTE: Current implementation produces ~0.047 for uniform 128 input.
        // C++ testdata shows mean ~0.08. This discrepancy needs investigation.
        // See CLAUDE.md section on adaptive quantization debugging.
    }

    #[test]
    fn test_ratio_of_derivatives() {
        // Basic sanity checks
        let r1 = ratio_of_derivatives(0.5, false);
        let r2 = ratio_of_derivatives(0.5, true);
        assert!(r1.is_finite());
        assert!(r2.is_finite());
        assert!(r1 > 0.0);
        assert!(r2 > 0.0);
    }

    #[test]
    fn test_update_min4() {
        let mut mins = [f32::INFINITY; 4];
        update_min4(5.0, &mut mins);
        assert!((mins[0] - 5.0).abs() < 1e-6);

        update_min4(3.0, &mut mins);
        assert!((mins[0] - 3.0).abs() < 1e-6);
        assert!((mins[1] - 5.0).abs() < 1e-6);

        update_min4(7.0, &mut mins);
        assert!((mins[0] - 3.0).abs() < 1e-6);
        assert!((mins[1] - 5.0).abs() < 1e-6);
        assert!((mins[2] - 7.0).abs() < 1e-6);
    }

    #[test]
    fn test_compute_mask_scalar() {
        let mask = compute_mask_scalar(1.0);
        assert!(mask.is_finite());
        assert!(mask > 0.0);
    }

    #[test]
    fn test_impl_runs_without_panic() {
        // Test that the full implementation runs without panicking
        let plane = vec![128.0f32; 64 * 64];
        let map = compute_aq_strength_map_impl(&plane, 64, 64, 1.0).unwrap();
        assert_eq!(map.width_blocks, 8);
        assert_eq!(map.height_blocks, 8);
        assert_eq!(map.strengths.len(), 64);

        // Values should be finite
        for &s in &map.strengths {
            assert!(s.is_finite(), "aq_strength should be finite");
        }
    }

    #[test]
    fn test_impl_output_range() {
        // Test with a smooth gradient image (more realistic than high-contrast pattern)
        // Each row has the same value, creating horizontal gradient
        let plane: Vec<f32> = (0..64 * 64)
            .map(|i| {
                let row = i / 64;
                (row * 4) as f32 // 0, 4, 8, ... 252 - smooth gradient
            })
            .collect();

        let distance = 1.0 / 3.0; // y_quant_01 = 3.0
        let map = compute_aq_strength_map_impl(&plane, 64, 64, distance).unwrap();

        for &s in &map.strengths {
            assert!(s >= 0.0, "aq_strength {} should be non-negative", s);
            assert!(s.is_finite(), "aq_strength {} should be finite", s);
        }

        // For smooth gradients, aq_strength should be relatively low
        // C++ produces 0.0-0.2 for typical images, but synthetic patterns can vary
        let mean: f32 = map.strengths.iter().sum::<f32>() / map.strengths.len() as f32;
        println!("aq_strength mean for gradient: {}", mean);

        // Just ensure it's bounded - the exact range depends on image content
        assert!(
            mean <= 2.0,
            "aq_strength mean {} is unexpectedly high even for synthetic image",
            mean
        );
    }
}