grib-reader 0.7.0

//! Data Section (Section 7) decoding.

use crate::error::{Error, Result};
use grib_core::binary::decode_ieee_binary128;
use grib_core::bit::BitReader;
pub use grib_core::data::{
    CcsdsBlockSize, CcsdsFlags, CcsdsPackingParams, ComplexPackingParams, DataRepresentation,
    Grib1SecondOrderLayout, Grib1SecondOrderPackingParams, IeeePackingParams, IeeePrecision,
    Jpeg2000PackingParams, PngPackingParams, RunLengthPackingParams, ScaledPackingParams,
    SpatialDifferencingParams,
};
use grib_core::filled_vec;

#[cfg(feature = "ccsds")]
use grib_aec::{decode as aec_decode, AecFlags, AecParams};

/// Numeric target type for decoded field values.
pub trait DecodeSample: Copy + Sized {
    fn from_f64(value: f64) -> Self;
    fn nan() -> Self;
}

impl DecodeSample for f32 {
    fn from_f64(value: f64) -> Self {
        value as f32
    }

    fn nan() -> Self {
        f32::NAN
    }
}

impl DecodeSample for f64 {
    fn from_f64(value: f64) -> Self {
        value
    }

    fn nan() -> Self {
        f64::NAN
    }
}

/// Decode Section 7 payload into field values, applying Section 6 bitmap when present.
pub fn decode_field(
    data_section: &[u8],
    representation: &DataRepresentation,
    bitmap_section: Option<&[u8]>,
    num_grid_points: usize,
) -> Result<Vec<f64>> {
    if data_section.len() < 5 || data_section[4] != 7 {
        return Err(Error::InvalidSection {
            section: data_section.get(4).copied().unwrap_or(7),
            reason: "not a data section".into(),
        });
    }

    decode_payload(
        &data_section[5..],
        representation,
        bitmap_section,
        num_grid_points,
    )
}

/// Decode Section 7 payload into a caller-provided buffer, applying Section 6
/// bitmap when present.
pub fn decode_field_into<T: DecodeSample>(
    data_section: &[u8],
    representation: &DataRepresentation,
    bitmap_section: Option<&[u8]>,
    num_grid_points: usize,
    out: &mut [T],
) -> Result<()> {
    if data_section.len() < 5 || data_section[4] != 7 {
        return Err(Error::InvalidSection {
            section: data_section.get(4).copied().unwrap_or(7),
            reason: "not a data section".into(),
        });
    }

    decode_payload_into_with_group_lengths(
        &data_section[5..],
        representation,
        bitmap_section,
        num_grid_points,
        out,
        None,
    )
}

pub(crate) fn decode_field_into_with_group_lengths<T: DecodeSample>(
    data_section: &[u8],
    representation: &DataRepresentation,
    bitmap_section: Option<&[u8]>,
    num_grid_points: usize,
    out: &mut [T],
    row_group_lengths: Option<&[usize]>,
) -> Result<()> {
    if data_section.len() < 5 || data_section[4] != 7 {
        return Err(Error::InvalidSection {
            section: data_section.get(4).copied().unwrap_or(7),
            reason: "not a data section".into(),
        });
    }

    decode_payload_into_with_group_lengths(
        &data_section[5..],
        representation,
        bitmap_section,
        num_grid_points,
        out,
        row_group_lengths,
    )
}

pub(crate) fn decode_payload(
    payload: &[u8],
    representation: &DataRepresentation,
    bitmap_section: Option<&[u8]>,
    num_grid_points: usize,
) -> Result<Vec<f64>> {
    let mut values = filled_vec(num_grid_points, 0.0, "decoded field")?;
    decode_payload_into_with_group_lengths(
        payload,
        representation,
        bitmap_section,
        num_grid_points,
        &mut values,
        None,
    )?;
    Ok(values)
}

pub(crate) fn decode_payload_into_with_group_lengths<T: DecodeSample>(
    payload: &[u8],
    representation: &DataRepresentation,
    bitmap_section: Option<&[u8]>,
    num_grid_points: usize,
    out: &mut [T],
    row_group_lengths: Option<&[usize]>,
) -> Result<()> {
    if out.len() != num_grid_points {
        return Err(Error::DataLengthMismatch {
            expected: num_grid_points,
            actual: out.len(),
        });
    }

    let expected_values = match representation {
        DataRepresentation::SimplePacking(params) => params.encoded_values,
        DataRepresentation::Grib1SecondOrderPacking(params) => params.packing().encoded_values,
        DataRepresentation::ComplexPacking(params) => params.encoded_values,
        DataRepresentation::IeeePacking(params) => params.encoded_values,
        DataRepresentation::Jpeg2000Packing(params) => params.packing.encoded_values,
        DataRepresentation::PngPacking(params) => params.packing.encoded_values,
        DataRepresentation::CcsdsPacking(params) => params.packing().encoded_values,
        DataRepresentation::RunLengthPacking(params) => params.encoded_values,
        DataRepresentation::Unsupported(template) => {
            return Err(Error::UnsupportedDataTemplate(*template));
        }
        _ => {
            return Err(Error::InvalidMessage(
                "data representation variant is not supported by this reader".into(),
            ));
        }
    };
    match bitmap_section {
        Some(bitmap_payload) => {
            let present_values = count_bitmap_present_points(bitmap_payload, num_grid_points)?;
            if expected_values != present_values {
                return Err(Error::DataLengthMismatch {
                    expected: present_values,
                    actual: expected_values,
                });
            }
        }
        None if expected_values != num_grid_points => {
            return Err(Error::DataLengthMismatch {
                expected: num_grid_points,
                actual: expected_values,
            });
        }
        None => {}
    }

    let mut output = OutputCursor::new(out, bitmap_section);
    match representation {
        DataRepresentation::SimplePacking(params) => {
            unpack_simple_into(payload, params, expected_values, &mut output)?
        }
        DataRepresentation::Grib1SecondOrderPacking(params) => {
            unpack_grib1_second_order_into(payload, params, row_group_lengths, &mut output)?
        }
        DataRepresentation::ComplexPacking(params) => {
            unpack_complex_into(payload, params, row_group_lengths, &mut output)?
        }
        DataRepresentation::IeeePacking(params) => {
            unpack_ieee_into(payload, params, expected_values, &mut output)?
        }
        DataRepresentation::Jpeg2000Packing(params) => {
            unpack_jpeg2000_into(payload, params, expected_values, &mut output)?
        }
        DataRepresentation::PngPacking(params) => {
            unpack_png_into(payload, params, expected_values, &mut output)?
        }
        DataRepresentation::CcsdsPacking(params) => {
            unpack_ccsds_into(payload, params, expected_values, &mut output)?
        }
        DataRepresentation::RunLengthPacking(params) => {
            unpack_run_length_into(payload, params, expected_values, &mut output)?
        }
        DataRepresentation::Unsupported(_) => unreachable!(),
        _ => unreachable!("unsupported representation rejected before decoding"),
    }
    output.finish()
}

fn unpack_ieee_into<T: DecodeSample>(
    data_bytes: &[u8],
    params: &IeeePackingParams,
    num_values: usize,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    let expected_bytes = num_values
        .checked_mul(params.precision.bytes_per_value())
        .ok_or(Error::ArithmeticOverflow {
            operation: "computing IEEE-packed byte count",
        })?;
    if data_bytes.len() != expected_bytes {
        return Err(Error::DataLengthMismatch {
            expected: expected_bytes,
            actual: data_bytes.len(),
        });
    }

    macro_rules! decode_values {
        ($chunks:expr, $decode:expr) => {{
            if let Some(dense) = output.take_dense(num_values)? {
                for (target, bytes) in dense.iter_mut().zip($chunks) {
                    *target = T::from_f64($decode(bytes));
                }
            } else {
                for bytes in $chunks {
                    output.push_present(T::from_f64($decode(bytes)))?;
                }
            }
        }};
    }

    match params.precision {
        IeeePrecision::Single => decode_values!(data_bytes.chunks_exact(4), |bytes: &[u8]| {
            f64::from(f32::from_be_bytes(
                bytes.try_into().expect("four-byte IEEE chunk"),
            ))
        }),
        IeeePrecision::Double => decode_values!(data_bytes.chunks_exact(8), |bytes: &[u8]| {
            f64::from_be_bytes(bytes.try_into().expect("eight-byte IEEE chunk"))
        }),
        IeeePrecision::Quadruple => {
            decode_values!(data_bytes.chunks_exact(16), |bytes: &[u8]| {
                decode_ieee_binary128(bytes.try_into().expect("sixteen-byte IEEE chunk"))
            })
        }
        _ => {
            return Err(Error::UnsupportedPackingWidth(
                params.precision.bits_per_value(),
            ));
        }
    }

    Ok(())
}

fn unpack_run_length_into<T: DecodeSample>(
    data_bytes: &[u8],
    params: &RunLengthPackingParams,
    num_values: usize,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    if params.bits_per_value > 64 {
        return Err(Error::InvalidSection {
            section: 5,
            reason: format!(
                "template 5.200 bits per value must not exceed 64, got {}",
                params.bits_per_value
            ),
        });
    }
    if params.max_level_value == 0 {
        if !params.level_values.is_empty() {
            return Err(Error::InvalidSection {
                section: 5,
                reason: format!(
                    "template 5.200 zero maximum level requires an empty level table, got {} entries",
                    params.level_values.len()
                ),
            });
        }
        if !remaining_bits_are_zero(data_bytes, 0)? {
            return Err(invalid_run_length_data(
                "all-missing field contains nonzero packed data",
            ));
        }
        let missing = T::nan();
        if let Some(dense) = output.take_dense(num_values)? {
            dense.fill(missing);
        } else {
            for _ in 0..num_values {
                output.push_present(missing)?;
            }
        }
        return Ok(());
    }

    let bits = usize::from(params.bits_per_value);
    if !(1..=64).contains(&bits) {
        return Err(Error::InvalidSection {
            section: 5,
            reason: format!(
                "template 5.200 bits per value must be between 1 and 64, got {}",
                params.bits_per_value
            ),
        });
    }
    if usize::from(params.max_level_value) > params.level_values.len() {
        return Err(Error::InvalidSection {
            section: 5,
            reason: format!(
                "template 5.200 maximum level {} exceeds the {}-entry level table",
                params.max_level_value,
                params.level_values.len()
            ),
        });
    }
    let max_code = if bits == 64 {
        u64::MAX
    } else {
        (1_u64 << bits) - 1
    };
    let range = max_code
        .checked_sub(u64::from(params.max_level_value))
        .filter(|range| *range > 0)
        .ok_or_else(|| Error::InvalidSection {
            section: 5,
            reason: format!(
                "template 5.200 maximum level {} leaves no run-length codes at {bits} bits",
                params.max_level_value
            ),
        })?;
    let data_bits = data_bytes
        .len()
        .checked_mul(8)
        .ok_or(Error::BitOffsetOverflow)?;
    let decimal_factor = 10.0_f64.powi(-i32::from(params.decimal_scale));
    let mut reader = BitReader::new(data_bytes);
    let mut decoded_values = 0usize;

    while decoded_values < num_values {
        if reader
            .bit_offset()
            .checked_add(bits)
            .ok_or(Error::BitOffsetOverflow)?
            > data_bits
        {
            return Err(Error::DataLengthMismatch {
                expected: num_values,
                actual: decoded_values,
            });
        }

        let level = reader.read(bits)?;
        if level > u64::from(params.max_level_value) {
            return Err(invalid_run_length_data(format!(
                "run starts with code {level}, above maximum level {}",
                params.max_level_value
            )));
        }

        let remaining_values = num_values - decoded_values;
        let run_limit = remaining_values as u128;
        let mut run_length = 1_u128;
        let mut factor = 1_u128;
        loop {
            if reader
                .bit_offset()
                .checked_add(bits)
                .ok_or(Error::BitOffsetOverflow)?
                > data_bits
            {
                break;
            }

            let mut lookahead = reader;
            let code = lookahead.read(bits)?;
            if code <= u64::from(params.max_level_value) {
                break;
            }
            reader = lookahead;

            let digit = u128::from(code - u64::from(params.max_level_value) - 1);
            let addition = factor.checked_mul(digit).ok_or(Error::ArithmeticOverflow {
                operation: "decoding a template 5.200 run length",
            })?;
            run_length = run_length
                .checked_add(addition)
                .ok_or(Error::ArithmeticOverflow {
                    operation: "decoding a template 5.200 run length",
                })?;
            if run_length > run_limit {
                return Err(Error::DataLengthMismatch {
                    expected: num_values,
                    actual: num_values.saturating_add(1),
                });
            }
            factor = factor.saturating_mul(u128::from(range)).min(run_limit + 1);
        }

        let run_length = usize::try_from(run_length).map_err(|_| Error::ArithmeticOverflow {
            operation: "converting a template 5.200 run length",
        })?;
        let value = if level == 0 {
            T::nan()
        } else {
            let scaled = params.level_values[usize::try_from(level - 1).map_err(|_| {
                Error::ArithmeticOverflow {
                    operation: "indexing a template 5.200 level",
                }
            })?];
            T::from_f64(f64::from(scaled) * decimal_factor)
        };
        if let Some(dense) = output.take_dense(run_length)? {
            dense.fill(value);
        } else {
            for _ in 0..run_length {
                output.push_present(value)?;
            }
        }
        decoded_values += run_length;
    }

    if !remaining_bits_are_zero(data_bytes, reader.bit_offset())? {
        return Err(invalid_run_length_data(
            "nonzero bits remain after the declared number of values",
        ));
    }
    Ok(())
}

fn remaining_bits_are_zero(data: &[u8], start_bit: usize) -> Result<bool> {
    let data_bits = data.len().checked_mul(8).ok_or(Error::BitOffsetOverflow)?;
    if start_bit > data_bits {
        return Err(Error::BitOffsetOverflow);
    }
    let mut byte_index = start_bit / 8;
    let bit_index = start_bit % 8;
    if bit_index != 0 {
        let remaining_mask = u8::MAX >> bit_index;
        if data[byte_index] & remaining_mask != 0 {
            return Ok(false);
        }
        byte_index += 1;
    }
    Ok(data[byte_index..].iter().all(|byte| *byte == 0))
}

fn invalid_run_length_data(reason: impl Into<String>) -> Error {
    Error::InvalidSection {
        section: 7,
        reason: format!("invalid template 5.200 data: {}", reason.into()),
    }
}

/// Parse bitmap presence from Section 6.
pub fn bitmap_payload(section_bytes: &[u8]) -> Result<Option<&[u8]>> {
    if section_bytes.len() < 6 {
        return Err(Error::InvalidSection {
            section: 6,
            reason: format!("expected at least 6 bytes, got {}", section_bytes.len()),
        });
    }
    if section_bytes[4] != 6 {
        return Err(Error::InvalidSection {
            section: section_bytes[4],
            reason: "not a bitmap section".into(),
        });
    }

    match section_bytes[5] {
        255 => Ok(None),
        0 => Ok(Some(&section_bytes[6..])),
        indicator => Err(Error::UnsupportedBitmapIndicator(u16::from(indicator))),
    }
}

pub(crate) fn count_bitmap_present_points(
    bitmap_payload: &[u8],
    num_grid_points: usize,
) -> Result<usize> {
    let full_bytes = num_grid_points / 8;
    let remaining_bits = num_grid_points % 8;
    let required_bytes = full_bytes + usize::from(remaining_bits > 0);
    if bitmap_payload.len() < required_bytes {
        return Err(Error::DataLengthMismatch {
            expected: required_bytes,
            actual: bitmap_payload.len(),
        });
    }

    let mut present = bitmap_payload[..full_bytes]
        .iter()
        .map(|byte| byte.count_ones() as usize)
        .sum();
    if remaining_bits > 0 {
        let mask = u8::MAX << (8 - remaining_bits);
        present += (bitmap_payload[full_bytes] & mask).count_ones() as usize;
    }

    Ok(present)
}

/// Unpack simple-packed values.
pub fn unpack_simple(
    data_bytes: &[u8],
    params: &ScaledPackingParams,
    num_values: usize,
) -> Result<Vec<f64>> {
    let mut values = filled_vec(num_values, 0.0, "simple-packed field")?;
    let mut output = OutputCursor::new(&mut values, None);
    unpack_simple_into(data_bytes, params, num_values, &mut output)?;
    output.finish()?;
    Ok(values)
}

fn unpack_simple_into<T: DecodeSample>(
    data_bytes: &[u8],
    params: &ScaledPackingParams,
    num_values: usize,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    let bits = params.bits_per_value as usize;
    let binary_factor = 2.0_f64.powi(params.binary_scale as i32);
    let decimal_factor = 10.0_f64.powi(-(params.decimal_scale as i32));
    let reference = params.reference_value as f64;
    if bits == 0 {
        let constant = T::from_f64(scale_decoded_value(
            reference,
            0.0,
            binary_factor,
            decimal_factor,
        ));
        if let Some(dense) = output.take_dense(num_values)? {
            dense.fill(constant);
            return Ok(());
        }
        for _ in 0..num_values {
            output.push_present(constant)?;
        }
        return Ok(());
    }
    if bits > u64::BITS as usize {
        return Err(Error::UnsupportedPackingWidth(params.bits_per_value));
    }

    let required_bits = bits
        .checked_mul(num_values)
        .ok_or(Error::ArithmeticOverflow {
            operation: "computing packed bit count",
        })?;
    let required_bytes = required_bits.div_ceil(8);
    if data_bytes.len() < required_bytes {
        return Err(Error::Truncated {
            offset: data_bytes.len() as u64,
        });
    }

    if let Some(dense) = output.take_dense(num_values)? {
        return unpack_simple_dense(
            data_bytes,
            bits,
            reference,
            binary_factor,
            decimal_factor,
            dense,
        );
    }

    let mut reader = BitReader::new(data_bytes);
    for _ in 0..num_values {
        let packed = reader.read(bits)?;
        output.push_present(T::from_f64(scale_decoded_value(
            reference,
            packed as f64,
            binary_factor,
            decimal_factor,
        )))?;
    }

    Ok(())
}

fn unpack_simple_dense<T: DecodeSample>(
    data_bytes: &[u8],
    bits: usize,
    reference: f64,
    binary_factor: f64,
    decimal_factor: f64,
    output: &mut [T],
) -> Result<()> {
    let scale = |packed: u64| {
        T::from_f64(scale_decoded_value(
            reference,
            packed as f64,
            binary_factor,
            decimal_factor,
        ))
    };

    match bits {
        8 => {
            for (target, &packed) in output.iter_mut().zip(data_bytes) {
                *target = scale(u64::from(packed));
            }
        }
        16 => {
            for (target, bytes) in output.iter_mut().zip(data_bytes.chunks_exact(2)) {
                *target = scale(u64::from(u16::from_be_bytes(
                    bytes.try_into().expect("two-byte chunk"),
                )));
            }
        }
        32 => {
            for (target, bytes) in output.iter_mut().zip(data_bytes.chunks_exact(4)) {
                *target = scale(u64::from(u32::from_be_bytes(
                    bytes.try_into().expect("four-byte chunk"),
                )));
            }
        }
        64 => {
            for (target, bytes) in output.iter_mut().zip(data_bytes.chunks_exact(8)) {
                *target = scale(u64::from_be_bytes(
                    bytes.try_into().expect("eight-byte chunk"),
                ));
            }
        }
        _ => {
            let mut reader = BitReader::new(data_bytes);
            for target in output {
                *target = scale(reader.read(bits)?);
            }
        }
    }
    Ok(())
}

#[cfg(feature = "ccsds")]
fn unpack_ccsds_into<T: DecodeSample>(
    data_bytes: &[u8],
    params: &CcsdsPackingParams,
    num_values: usize,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    let packing = params.packing();
    let bits = packing.bits_per_value;

    if bits == 0 {
        if !data_bytes.is_empty() {
            return Err(Error::DataLengthMismatch {
                expected: 0,
                actual: data_bytes.len(),
            });
        }
        // ecCodes stores a constant CCSDS field directly in referenceValue and
        // ignores both scale factors when bitsPerValue is zero.
        let value = T::from_f64(f64::from(packing.reference_value));
        if let Some(dense) = output.take_dense(num_values)? {
            dense.fill(value);
        } else {
            for _ in 0..num_values {
                output.push_present(value)?;
            }
        }
        return Ok(());
    }

    if num_values == 0 {
        if data_bytes.is_empty() {
            return Ok(());
        }
        return Err(Error::DataLengthMismatch {
            expected: 0,
            actual: data_bytes.len(),
        });
    }

    let bytes_per_sample = usize::from(bits).div_ceil(8);
    let decoded_size =
        num_values
            .checked_mul(bytes_per_sample)
            .ok_or(Error::ArithmeticOverflow {
                operation: "computing CCSDS decoded byte count",
            })?;
    let codec_params = AecParams::new(
        bits,
        params.block_size().samples(),
        params.reference_sample_interval(),
        ccsds_decoder_flags(params.flags()),
    )
    .map_err(|error| Error::Codec {
        codec: "CCSDS/AEC",
        operation: "configure",
        reason: error.to_string(),
    })?;
    let decoded =
        aec_decode(data_bytes, num_values, codec_params).map_err(|error| Error::Codec {
            codec: "CCSDS/AEC",
            operation: "decode",
            reason: error.to_string(),
        })?;
    if decoded.len() != decoded_size {
        return Err(Error::DataLengthMismatch {
            expected: decoded_size,
            actual: decoded.len(),
        });
    }

    let value_mask = if bits == 32 {
        u32::MAX
    } else {
        (1u32 << bits) - 1
    };
    for bytes in decoded.chunks_exact(bytes_per_sample) {
        let container = bytes
            .iter()
            .fold(0u32, |value, byte| (value << 8) | u32::from(*byte));
        let raw = container & value_mask;
        let signed = params.flags().contains(CcsdsFlags::SIGNED);
        let high_mask = if bytes_per_sample == 4 {
            !value_mask
        } else {
            ((1u32 << (bytes_per_sample * 8)) - 1) & !value_mask
        };
        let expected_high = if signed && raw & (1u32 << (bits - 1)) != 0 {
            high_mask
        } else {
            0
        };
        if container & high_mask != expected_high {
            return Err(Error::ValueOutOfRange(format!(
                "CCSDS sample container 0x{container:08x} has invalid extension bits for its {bits}-bit width"
            )));
        }
        // SIGNED controls the AEC preprocessor. GRIB template 5.42 still
        // stores non-negative scaled differences, so reconstruction uses the
        // declared low bits as an unsigned value.
        output.push_present(scale_ccsds_value(packing, f64::from(raw)))?;
    }
    Ok(())
}

#[cfg(feature = "ccsds")]
fn ccsds_decoder_flags(flags: CcsdsFlags) -> AecFlags {
    // MSB and THREE_BYTE describe the caller's decoded sample containers, not
    // the CCSDS bitstream. Request one canonical, compact big-endian layout so
    // the reader never depends on host endianness or four-byte padding.
    let mut codec_flags = AecFlags::DATA_MSB | AecFlags::DATA_3BYTE;
    for (wire_flag, codec_flag) in [
        (CcsdsFlags::SIGNED, AecFlags::DATA_SIGNED),
        (CcsdsFlags::PREPROCESS, AecFlags::DATA_PREPROCESS),
        (CcsdsFlags::RESTRICTED, AecFlags::RESTRICTED),
        (CcsdsFlags::PAD_RSI, AecFlags::PAD_RSI),
        (CcsdsFlags::NOT_ENFORCE, AecFlags::NOT_ENFORCE),
    ] {
        if flags.contains(wire_flag) {
            codec_flags |= codec_flag;
        }
    }
    codec_flags
}

#[cfg(feature = "ccsds")]
fn scale_ccsds_value<T: DecodeSample>(params: &ScaledPackingParams, packed: f64) -> T {
    T::from_f64(scale_decoded_value(
        f64::from(params.reference_value),
        packed,
        2.0_f64.powi(i32::from(params.binary_scale)),
        10.0_f64.powi(-i32::from(params.decimal_scale)),
    ))
}

#[cfg(not(feature = "ccsds"))]
fn unpack_ccsds_into<T: DecodeSample>(
    _data_bytes: &[u8],
    _params: &CcsdsPackingParams,
    _num_values: usize,
    _output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    Err(Error::UnsupportedDataTemplate(42))
}

#[cfg(feature = "jpeg2000")]
fn unpack_jpeg2000_into<T: DecodeSample>(
    data_bytes: &[u8],
    params: &Jpeg2000PackingParams,
    num_values: usize,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    validate_jpeg2000_bits(params.packing.bits_per_value)?;

    // Parse the SIZ marker in safe Rust before OpenJPEG allocates or decodes
    // sample data. jpeg2k's header-only DumpImage path currently frees a null
    // OpenJPEG code-block pointer on some valid codestreams, so it cannot be
    // used as the validation boundary.
    let header = parse_jpeg2000_header(data_bytes)?;
    if header.components != 1 {
        return Err(invalid_data_section(format!(
            "JPEG 2000 GRIB packing requires one component, got {}",
            header.components
        )));
    }
    if header.is_signed {
        return Err(invalid_data_section(
            "JPEG 2000 GRIB packing requires unsigned component data",
        ));
    }
    if header.precision > params.packing.bits_per_value {
        return Err(invalid_data_section(format!(
            "JPEG 2000 component precision {} exceeds GRIB bits-per-value {}",
            header.precision, params.packing.bits_per_value
        )));
    }

    let sample_count = image_sample_count(header.width, header.height)?;
    if sample_count != num_values {
        return Err(Error::DataLengthMismatch {
            expected: num_values,
            actual: sample_count,
        });
    }

    let image = jpeg2k::Image::from_bytes(data_bytes).map_err(|err| Error::Codec {
        codec: "JPEG 2000",
        operation: "decode",
        reason: err.to_string(),
    })?;
    if image.num_components() != 1 {
        return Err(invalid_data_section(format!(
            "JPEG 2000 GRIB packing requires one component, got {}",
            image.num_components()
        )));
    }
    let component = image
        .components()
        .first()
        .ok_or_else(|| invalid_data_section("JPEG 2000 image has no components"))?;
    if component.is_signed() {
        return Err(invalid_data_section(
            "JPEG 2000 GRIB packing requires unsigned component data",
        ));
    }
    if component.precision() > u32::from(params.packing.bits_per_value) {
        return Err(invalid_data_section(format!(
            "JPEG 2000 component precision {} exceeds GRIB bits-per-value {}",
            component.precision(),
            params.packing.bits_per_value
        )));
    }
    let decoded_sample_count = image_sample_count(component.width(), component.height())?;
    if decoded_sample_count != num_values {
        return Err(Error::DataLengthMismatch {
            expected: num_values,
            actual: decoded_sample_count,
        });
    }

    for &sample in component.data() {
        let raw = jpeg2000_raw_value(sample, params.packing.bits_per_value)?;
        push_image_value(output, &params.packing, raw)?;
    }

    Ok(())
}

#[cfg(feature = "jpeg2000")]
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
struct Jpeg2000Header {
    width: u32,
    height: u32,
    components: u16,
    precision: u8,
    is_signed: bool,
}

#[cfg(feature = "jpeg2000")]
fn parse_jpeg2000_header(data: &[u8]) -> Result<Jpeg2000Header> {
    let codestream = jpeg2000_codestream(data)?;
    if codestream.get(..4) != Some(&[0xff, 0x4f, 0xff, 0x51]) {
        return Err(invalid_jpeg2000_header(
            "codestream does not begin with SOC followed by SIZ",
        ));
    }

    let segment_length = usize::from(u16::from_be_bytes(codestream[4..6].try_into().unwrap()));
    let segment_end = 4usize
        .checked_add(segment_length)
        .ok_or_else(|| invalid_jpeg2000_header("SIZ segment length overflow"))?;
    let segment = codestream
        .get(4..segment_end)
        .ok_or_else(|| invalid_jpeg2000_header("truncated SIZ segment"))?;
    if segment.len() < 38 {
        return Err(invalid_jpeg2000_header("SIZ segment is too short"));
    }

    let components = u16::from_be_bytes(segment[36..38].try_into().unwrap());
    let required_length = usize::from(components)
        .checked_mul(3)
        .and_then(|component_bytes| 38usize.checked_add(component_bytes))
        .ok_or_else(|| invalid_jpeg2000_header("SIZ component table length overflow"))?;
    if segment_length != required_length {
        return Err(invalid_jpeg2000_header(format!(
            "SIZ length {segment_length} does not match {components} components"
        )));
    }
    if components == 0 {
        return Err(invalid_jpeg2000_header("SIZ declares zero components"));
    }

    let xsiz = u32::from_be_bytes(segment[4..8].try_into().unwrap());
    let ysiz = u32::from_be_bytes(segment[8..12].try_into().unwrap());
    let xosiz = u32::from_be_bytes(segment[12..16].try_into().unwrap());
    let yosiz = u32::from_be_bytes(segment[16..20].try_into().unwrap());
    let xrsiz = u32::from(segment[39]);
    let yrsiz = u32::from(segment[40]);
    if xrsiz == 0 || yrsiz == 0 {
        return Err(invalid_jpeg2000_header(
            "SIZ component subsampling factors must be nonzero",
        ));
    }
    let width = ceil_div_u32(xsiz, xrsiz)
        .checked_sub(ceil_div_u32(xosiz, xrsiz))
        .filter(|width| *width > 0)
        .ok_or_else(|| invalid_jpeg2000_header("SIZ has an invalid horizontal extent"))?;
    let height = ceil_div_u32(ysiz, yrsiz)
        .checked_sub(ceil_div_u32(yosiz, yrsiz))
        .filter(|height| *height > 0)
        .ok_or_else(|| invalid_jpeg2000_header("SIZ has an invalid vertical extent"))?;
    let sample_format = segment[38];

    Ok(Jpeg2000Header {
        width,
        height,
        components,
        precision: (sample_format & 0x7f) + 1,
        is_signed: sample_format & 0x80 != 0,
    })
}

#[cfg(feature = "jpeg2000")]
fn jpeg2000_codestream(data: &[u8]) -> Result<&[u8]> {
    if data.starts_with(&[0xff, 0x4f]) {
        return Ok(data);
    }

    let mut offset = 0usize;
    while offset
        .checked_add(8)
        .is_some_and(|header_end| header_end <= data.len())
    {
        let short_length = u32::from_be_bytes(data[offset..offset + 4].try_into().unwrap());
        let box_type = &data[offset + 4..offset + 8];
        let (header_length, box_length) = match short_length {
            0 => (8usize, data.len() - offset),
            1 => {
                let extended_end = offset
                    .checked_add(16)
                    .ok_or_else(|| invalid_jpeg2000_header("JP2 box header overflow"))?;
                let extended = data
                    .get(offset + 8..extended_end)
                    .ok_or_else(|| invalid_jpeg2000_header("truncated extended JP2 box"))?;
                let length = usize::try_from(u64::from_be_bytes(extended.try_into().unwrap()))
                    .map_err(|_| invalid_jpeg2000_header("JP2 box length exceeds usize"))?;
                (16usize, length)
            }
            length => (8usize, length as usize),
        };
        if box_length < header_length {
            return Err(invalid_jpeg2000_header(
                "JP2 box is shorter than its header",
            ));
        }
        let box_end = offset
            .checked_add(box_length)
            .filter(|end| *end <= data.len())
            .ok_or_else(|| invalid_jpeg2000_header("JP2 box extends past the payload"))?;
        if box_type == b"jp2c" {
            return data
                .get(offset + header_length..box_end)
                .ok_or_else(|| invalid_jpeg2000_header("truncated JP2 codestream box"));
        }
        if short_length == 0 {
            break;
        }
        offset = box_end;
    }

    Err(invalid_jpeg2000_header(
        "payload is neither a raw codestream nor a JP2 file with a codestream box",
    ))
}

#[cfg(feature = "jpeg2000")]
const fn ceil_div_u32(value: u32, divisor: u32) -> u32 {
    value / divisor + (!value.is_multiple_of(divisor)) as u32
}

#[cfg(feature = "jpeg2000")]
fn invalid_jpeg2000_header(reason: impl Into<String>) -> Error {
    invalid_data_section(format!("invalid JPEG 2000 header: {}", reason.into()))
}

#[cfg(any(feature = "jpeg2000", feature = "png"))]
fn invalid_data_section(reason: impl Into<String>) -> Error {
    Error::InvalidSection {
        section: 7,
        reason: reason.into(),
    }
}

#[cfg(not(feature = "jpeg2000"))]
fn unpack_jpeg2000_into<T: DecodeSample>(
    _data_bytes: &[u8],
    _params: &Jpeg2000PackingParams,
    _num_values: usize,
    _output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    Err(Error::UnsupportedDataTemplate(40))
}

#[cfg(feature = "jpeg2000")]
fn validate_jpeg2000_bits(bits_per_value: u8) -> Result<()> {
    if !(1..=32).contains(&bits_per_value) {
        return Err(Error::UnsupportedPackingWidth(bits_per_value));
    }
    Ok(())
}

#[cfg(feature = "jpeg2000")]
fn jpeg2000_raw_value(sample: i32, bits_per_value: u8) -> Result<u64> {
    let raw = if bits_per_value == 32 {
        u64::from(u32::from_ne_bytes(sample.to_ne_bytes()))
    } else {
        u64::try_from(sample).map_err(|_| {
            Error::ValueOutOfRange("JPEG 2000 unsigned component yielded a negative sample".into())
        })?
    };
    validate_raw_value_fits(raw, bits_per_value)?;
    Ok(raw)
}

#[cfg(feature = "png")]
fn unpack_png_into<T: DecodeSample>(
    data_bytes: &[u8],
    params: &PngPackingParams,
    num_values: usize,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    validate_png_bits(params.packing.bits_per_value)?;

    let decoder = png::Decoder::new(std::io::Cursor::new(data_bytes));
    let mut reader = decoder.read_info().map_err(|err| Error::Codec {
        codec: "PNG",
        operation: "read header",
        reason: err.to_string(),
    })?;
    let header_sample_count = image_sample_count(reader.info().width, reader.info().height)?;
    if header_sample_count != num_values {
        return Err(Error::DataLengthMismatch {
            expected: num_values,
            actual: header_sample_count,
        });
    }
    let buffer_size = reader
        .output_buffer_size()
        .ok_or(Error::ArithmeticOverflow {
            operation: "computing PNG output buffer size",
        })?;
    let mut buffer = filled_vec(buffer_size, 0, "PNG output buffer")?;
    let info = reader.next_frame(&mut buffer).map_err(|err| Error::Codec {
        codec: "PNG",
        operation: "decode",
        reason: err.to_string(),
    })?;
    let data = &buffer[..info.buffer_size()];

    let sample_count = image_sample_count(info.width, info.height)?;
    if sample_count != num_values {
        return Err(Error::DataLengthMismatch {
            expected: num_values,
            actual: sample_count,
        });
    }

    match (
        info.color_type,
        info.bit_depth,
        params.packing.bits_per_value,
    ) {
        (png::ColorType::Grayscale, png::BitDepth::One, 1)
        | (png::ColorType::Grayscale, png::BitDepth::Two, 2)
        | (png::ColorType::Grayscale, png::BitDepth::Four, 4) => unpack_png_subbyte_grayscale(
            data,
            info.width,
            info.height,
            params.packing.bits_per_value,
            &params.packing,
            output,
        ),
        (png::ColorType::Grayscale, png::BitDepth::Eight, 8) => {
            unpack_png_bytes(data, 1, num_values, &params.packing, output)
        }
        (png::ColorType::Grayscale, png::BitDepth::Sixteen, 16) => {
            unpack_png_u16(data, num_values, &params.packing, output)
        }
        (png::ColorType::Rgb, png::BitDepth::Eight, 24) => {
            unpack_png_bytes(data, 3, num_values, &params.packing, output)
        }
        (png::ColorType::Rgba, png::BitDepth::Eight, 32) => {
            unpack_png_bytes(data, 4, num_values, &params.packing, output)
        }
        (color_type, bit_depth, bits_per_value) => Err(invalid_data_section(format!(
            "PNG image layout {color_type:?}/{bit_depth:?} is incompatible with GRIB bits-per-value {bits_per_value}"
        ))),
    }
}

#[cfg(not(feature = "png"))]
fn unpack_png_into<T: DecodeSample>(
    _data_bytes: &[u8],
    _params: &PngPackingParams,
    _num_values: usize,
    _output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    Err(Error::UnsupportedDataTemplate(41))
}

#[cfg(feature = "png")]
fn validate_png_bits(bits_per_value: u8) -> Result<()> {
    if !matches!(bits_per_value, 1 | 2 | 4 | 8 | 16 | 24 | 32) {
        return Err(Error::UnsupportedPackingWidth(bits_per_value));
    }
    Ok(())
}

#[cfg(any(feature = "jpeg2000", feature = "png"))]
fn image_sample_count(width: u32, height: u32) -> Result<usize> {
    let width = usize::try_from(width)
        .map_err(|_| Error::ValueOutOfRange("image width exceeds usize".into()))?;
    let height = usize::try_from(height)
        .map_err(|_| Error::ValueOutOfRange("image height exceeds usize".into()))?;
    width.checked_mul(height).ok_or(Error::ArithmeticOverflow {
        operation: "computing image sample count",
    })
}

#[cfg(feature = "png")]
fn unpack_png_subbyte_grayscale<T: DecodeSample>(
    data: &[u8],
    width: u32,
    height: u32,
    bits_per_sample: u8,
    params: &ScaledPackingParams,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    let width = usize::try_from(width)
        .map_err(|_| Error::ValueOutOfRange("PNG width exceeds usize".into()))?;
    let height = usize::try_from(height)
        .map_err(|_| Error::ValueOutOfRange("PNG height exceeds usize".into()))?;
    let bits = usize::from(bits_per_sample);
    let row_bytes = width
        .checked_mul(bits)
        .ok_or(Error::ArithmeticOverflow {
            operation: "computing PNG row width",
        })?
        .div_ceil(8);
    let expected_bytes = row_bytes
        .checked_mul(height)
        .ok_or(Error::ArithmeticOverflow {
            operation: "computing PNG data length",
        })?;
    if data.len() < expected_bytes {
        return Err(Error::DataLengthMismatch {
            expected: expected_bytes,
            actual: data.len(),
        });
    }

    let mask = (1u8 << bits) - 1;
    for row in data[..expected_bytes].chunks_exact(row_bytes) {
        for x in 0..width {
            let bit_offset = x.checked_mul(bits).ok_or(Error::BitOffsetOverflow)?;
            let shift = 8 - bits - (bit_offset % 8);
            let raw = u64::from((row[bit_offset / 8] >> shift) & mask);
            push_image_value(output, params, raw)?;
        }
    }

    Ok(())
}

#[cfg(feature = "png")]
fn unpack_png_bytes<T: DecodeSample>(
    data: &[u8],
    bytes_per_sample: usize,
    num_values: usize,
    params: &ScaledPackingParams,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    let expected_bytes =
        num_values
            .checked_mul(bytes_per_sample)
            .ok_or(Error::ArithmeticOverflow {
                operation: "computing PNG data length",
            })?;
    if data.len() < expected_bytes {
        return Err(Error::DataLengthMismatch {
            expected: expected_bytes,
            actual: data.len(),
        });
    }

    for sample in data[..expected_bytes].chunks_exact(bytes_per_sample) {
        let raw = sample
            .iter()
            .fold(0u64, |acc, byte| (acc << 8) | u64::from(*byte));
        push_image_value(output, params, raw)?;
    }

    Ok(())
}

#[cfg(feature = "png")]
fn unpack_png_u16<T: DecodeSample>(
    data: &[u8],
    num_values: usize,
    params: &ScaledPackingParams,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    let expected_bytes = num_values.checked_mul(2).ok_or(Error::ArithmeticOverflow {
        operation: "computing PNG data length",
    })?;
    if data.len() < expected_bytes {
        return Err(Error::DataLengthMismatch {
            expected: expected_bytes,
            actual: data.len(),
        });
    }

    for sample in data[..expected_bytes].chunks_exact(2) {
        let raw = u64::from(u16::from_be_bytes(sample.try_into().unwrap()));
        push_image_value(output, params, raw)?;
    }

    Ok(())
}

#[cfg(any(feature = "jpeg2000", feature = "png"))]
fn push_image_value<T: DecodeSample>(
    output: &mut OutputCursor<'_, T>,
    params: &ScaledPackingParams,
    raw: u64,
) -> Result<()> {
    validate_raw_value_fits(raw, params.bits_per_value)?;
    output.push_present(scale_image_value(params, raw))
}

#[cfg(any(feature = "jpeg2000", feature = "png"))]
fn validate_raw_value_fits(raw: u64, bits_per_value: u8) -> Result<()> {
    let max_value = if bits_per_value == u64::BITS as u8 {
        u64::MAX
    } else {
        (1u64 << bits_per_value) - 1
    };
    if raw > max_value {
        return Err(Error::ValueOutOfRange(format!(
            "decoded image sample {raw} exceeds GRIB bits-per-value {bits_per_value}"
        )));
    }
    Ok(())
}

#[cfg(any(feature = "jpeg2000", feature = "png"))]
fn scale_image_value<T: DecodeSample>(params: &ScaledPackingParams, raw: u64) -> T {
    let binary_factor = 2.0_f64.powi(params.binary_scale as i32);
    let decimal_factor = 10.0_f64.powi(-(params.decimal_scale as i32));
    T::from_f64(scale_decoded_value(
        params.reference_value as f64,
        raw as f64,
        binary_factor,
        decimal_factor,
    ))
}

#[cfg(test)]
fn unpack_complex(data_bytes: &[u8], params: &ComplexPackingParams) -> Result<Vec<f64>> {
    let mut values = filled_vec(params.encoded_values, 0.0, "complex-packed field")?;
    let mut output = OutputCursor::new(&mut values, None);
    unpack_complex_into(data_bytes, params, None, &mut output)?;
    output.finish()?;
    Ok(values)
}

fn unpack_grib1_second_order_into<T: DecodeSample>(
    data_bytes: &[u8],
    params: &Grib1SecondOrderPackingParams,
    row_group_lengths: Option<&[usize]>,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    let group_lengths = match params.layout() {
        Grib1SecondOrderLayout::RowByRow => {
            row_group_lengths.ok_or_else(|| Error::InvalidSection {
                section: 4,
                reason: "GRIB1 row-by-row second-order packing requires grid-derived group lengths"
                    .into(),
            })?
        }
        Grib1SecondOrderLayout::ConstantWidth | Grib1SecondOrderLayout::VariableWidth => {
            params.group_lengths()
        }
        _ => {
            return Err(Error::InvalidSection {
                section: 4,
                reason: "unsupported GRIB1 second-order group layout".into(),
            });
        }
    };
    validate_external_group_lengths(
        group_lengths,
        params.group_count(),
        params.packing().encoded_values,
        "GRIB1 second-order",
        4,
    )?;

    let first_order_offset = params
        .first_order_values_offset()
        .checked_mul(8)
        .ok_or(Error::BitOffsetOverflow)?;
    let second_order_offset = params
        .second_order_values_offset()
        .checked_mul(8)
        .ok_or(Error::BitOffsetOverflow)?;
    let mut first_order_reader = BitReader::with_offset(data_bytes, first_order_offset);
    let mut second_order_reader = BitReader::with_offset(data_bytes, second_order_offset);
    let packing = params.packing();
    let binary_factor = 2.0_f64.powi(i32::from(packing.binary_scale));
    let decimal_factor = 10.0_f64.powi(-i32::from(packing.decimal_scale));
    let reference = f64::from(packing.reference_value);

    for (group_index, group_length) in group_lengths.iter().copied().enumerate() {
        let first_order = first_order_reader.read(usize::from(packing.bits_per_value))?;
        let group_width = params
            .group_width(group_index)
            .ok_or_else(|| Error::InvalidSection {
                section: 4,
                reason: format!("missing width for GRIB1 second-order group {group_index}"),
            })?;
        for _ in 0..group_length {
            let second_order = second_order_reader.read(usize::from(group_width))?;
            let packed =
                first_order
                    .checked_add(second_order)
                    .ok_or(Error::ArithmeticOverflow {
                        operation: "reconstructing a GRIB1 second-order value",
                    })?;
            output.push_present(T::from_f64(scale_decoded_value(
                reference,
                packed as f64,
                binary_factor,
                decimal_factor,
            )))?;
        }
    }

    let payload_bits = data_bytes
        .len()
        .checked_mul(8)
        .ok_or(Error::BitOffsetOverflow)?;
    let usable_bits = payload_bits
        .checked_sub(usize::from(params.unused_bits()))
        .ok_or_else(|| Error::InvalidSection {
            section: 4,
            reason: "GRIB1 BDS unused-bit count exceeds its payload".into(),
        })?;
    if second_order_reader.bit_offset() != usable_bits {
        return Err(Error::InvalidSection {
            section: 4,
            reason: format!(
                "GRIB1 second-order values end at bit {}, expected {usable_bits}",
                second_order_reader.bit_offset()
            ),
        });
    }
    if output.values_written() != packing.encoded_values {
        return Err(Error::DataLengthMismatch {
            expected: packing.encoded_values,
            actual: output.values_written(),
        });
    }
    Ok(())
}

fn unpack_complex_into<T: DecodeSample>(
    data_bytes: &[u8],
    params: &ComplexPackingParams,
    row_group_lengths: Option<&[usize]>,
    output: &mut OutputCursor<'_, T>,
) -> Result<()> {
    if params.num_groups == 0 {
        return Err(Error::InvalidSection {
            section: 5,
            reason: "complex packing requires at least one group".into(),
        });
    }

    let mut reader = BitReader::new(data_bytes);
    let mut spatial = params
        .spatial_differencing
        .map(|spatial| read_spatial_descriptors(&mut reader, spatial))
        .transpose()?
        .map(SpatialRestoreState::new);

    let layout = GroupReaderLayout::new(reader.bit_offset(), params)?;
    let mut reference_reader = BitReader::with_offset(data_bytes, layout.reference_offset);
    let mut width_reader = BitReader::with_offset(data_bytes, layout.width_offset);
    let mut length_reader = BitReader::with_offset(data_bytes, layout.length_offset);
    let mut value_reader = BitReader::with_offset(data_bytes, layout.value_offset);
    let binary_factor = 2.0_f64.powi(params.binary_scale as i32);
    let decimal_factor = 10.0_f64.powi(-(params.decimal_scale as i32));
    let reference = params.reference_value as f64;
    let mut actual_total = 0usize;
    let row_group_lengths = match params.group_splitting_method {
        0 => {
            let lengths = row_group_lengths.ok_or_else(|| Error::InvalidSection {
                section: 5,
                reason:
                    "complex packing group splitting method 0 requires grid-derived group lengths"
                        .into(),
            })?;
            validate_external_group_lengths(
                lengths,
                params.num_groups,
                params.encoded_values,
                "complex-packing row",
                5,
            )?;
            Some(lengths)
        }
        1 => None,
        method => return Err(Error::UnsupportedGroupSplittingMethod(method)),
    };

    for group_index in 0..params.num_groups {
        let group_reference = reference_reader.read(params.group_reference_bits as usize)?;
        let width_delta = width_reader.read(params.group_width_bits as usize)?;
        let group_width = usize::from(params.group_width_reference)
            .checked_add(width_delta as usize)
            .ok_or(Error::ArithmeticOverflow {
                operation: "computing complex-packing group width",
            })?;
        let group_length = match row_group_lengths {
            Some(lengths) => lengths[group_index],
            None => read_group_length(&mut length_reader, params, group_index)?,
        };

        actual_total = actual_total
            .checked_add(group_length)
            .ok_or(Error::ArithmeticOverflow {
                operation: "summing complex-packing group lengths",
            })?;

        if group_width == 0 {
            let raw_value = decode_constant_group_value(
                group_reference,
                params.group_reference_bits as usize,
                params.missing_value_management,
            )?;
            for _ in 0..group_length {
                let value = spatial
                    .as_mut()
                    .map_or(Ok(raw_value), |state| state.restore(raw_value))?;
                output.push_present(scale_complex_value(
                    reference,
                    binary_factor,
                    decimal_factor,
                    value,
                ))?;
            }
            continue;
        }

        if group_width > u64::BITS as usize {
            return Err(Error::UnsupportedPackingWidth(group_width as u8));
        }

        let group_reference = i64::try_from(group_reference)
            .map_err(|_| Error::ValueOutOfRange("group reference exceeds i64".into()))?;
        for _ in 0..group_length {
            let packed = value_reader.read(group_width)?;
            let value = decode_group_value(
                group_reference,
                packed,
                group_width,
                params.missing_value_management,
            )?;
            let value = spatial
                .as_mut()
                .map_or(Ok(value), |state| state.restore(value))?;
            output.push_present(scale_complex_value(
                reference,
                binary_factor,
                decimal_factor,
                value,
            ))?;
        }
    }

    if actual_total != params.encoded_values {
        return Err(Error::DataLengthMismatch {
            expected: params.encoded_values,
            actual: actual_total,
        });
    }

    if output.values_written() != params.encoded_values {
        return Err(Error::DataLengthMismatch {
            expected: params.encoded_values,
            actual: output.values_written(),
        });
    }

    if let Some(spatial) = spatial {
        spatial.finish()?;
    }

    Ok(())
}

fn validate_external_group_lengths(
    lengths: &[usize],
    expected_groups: usize,
    expected_values: usize,
    packing_name: &str,
    section: u8,
) -> Result<()> {
    if lengths.len() != expected_groups {
        return Err(Error::InvalidSection {
            section,
            reason: format!(
                "{packing_name} packing declares {expected_groups} groups, but the grid defines {} scan lines",
                lengths.len()
            ),
        });
    }
    let actual_values = lengths.iter().try_fold(0usize, |total, length| {
        total.checked_add(*length).ok_or(Error::ArithmeticOverflow {
            operation: "summing row-wise packing group lengths",
        })
    })?;
    if actual_values != expected_values {
        return Err(Error::DataLengthMismatch {
            expected: expected_values,
            actual: actual_values,
        });
    }
    Ok(())
}

fn read_group_length(
    reader: &mut BitReader<'_>,
    params: &ComplexPackingParams,
    group_index: usize,
) -> Result<usize> {
    if params.scaled_group_length_bits as usize > u64::BITS as usize {
        return Err(Error::UnsupportedPackingWidth(
            params.scaled_group_length_bits,
        ));
    }

    if group_index + 1 == params.num_groups {
        return Ok(params.true_length_last_group as usize);
    }

    let scaled = reader
        .read(params.scaled_group_length_bits as usize)?
        .checked_mul(u64::from(params.group_length_increment))
        .ok_or(Error::ArithmeticOverflow {
            operation: "scaling complex-packing group length",
        })?;
    let length = u64::from(params.group_length_reference)
        .checked_add(scaled)
        .ok_or(Error::ArithmeticOverflow {
            operation: "computing complex-packing group length",
        })?;
    usize::try_from(length)
        .map_err(|_| Error::ValueOutOfRange("complex-packing group length exceeds usize".into()))
}

fn read_spatial_descriptors(
    reader: &mut BitReader<'_>,
    params: SpatialDifferencingParams,
) -> Result<SpatialDescriptors> {
    if params.descriptor_octets == 0 {
        return Err(Error::InvalidSection {
            section: 5,
            reason: "spatial differencing requires at least one descriptor octet".into(),
        });
    }

    let bit_count = usize::from(params.descriptor_octets) * 8;
    if bit_count > u64::BITS as usize {
        return Err(Error::InvalidSection {
            section: 5,
            reason: "spatial differencing descriptors wider than 8 octets are not supported".into(),
        });
    }

    let first_value = reader.read_signed(bit_count)?;
    let second_value = if params.order == 2 {
        Some(reader.read_signed(bit_count)?)
    } else {
        None
    };
    for _ in usize::from(params.order.min(2))..params.order as usize {
        let _ = reader.read_signed(bit_count)?;
    }
    let overall_minimum = reader.read_signed(bit_count)?;

    Ok(SpatialDescriptors {
        order: params.order,
        first_value,
        second_value,
        overall_minimum,
    })
}

fn decode_constant_group_value(
    group_reference: u64,
    group_reference_bits: usize,
    missing_value_management: u8,
) -> Result<Option<i64>> {
    if is_missing_code(
        group_reference,
        group_reference_bits,
        missing_value_management,
        MissingKind::Primary,
    )? || is_missing_code(
        group_reference,
        group_reference_bits,
        missing_value_management,
        MissingKind::Secondary,
    )? {
        return Ok(None);
    }

    let value = i64::try_from(group_reference)
        .map_err(|_| Error::ValueOutOfRange("group reference exceeds i64".into()))?;
    Ok(Some(value))
}

fn decode_group_value(
    group_reference: i64,
    packed: u64,
    group_width: usize,
    missing_value_management: u8,
) -> Result<Option<i64>> {
    if is_missing_code(
        packed,
        group_width,
        missing_value_management,
        MissingKind::Primary,
    )? || is_missing_code(
        packed,
        group_width,
        missing_value_management,
        MissingKind::Secondary,
    )? {
        return Ok(None);
    }

    let packed = i64::try_from(packed)
        .map_err(|_| Error::ValueOutOfRange("packed value exceeds i64".into()))?;
    let value = group_reference
        .checked_add(packed)
        .ok_or(Error::ArithmeticOverflow {
            operation: "reconstructing complex-packing value",
        })?;
    Ok(Some(value))
}

fn is_missing_code(
    value: u64,
    bit_width: usize,
    missing_value_management: u8,
    kind: MissingKind,
) -> Result<bool> {
    let required_mode = match kind {
        MissingKind::Primary => 1,
        MissingKind::Secondary => 2,
    };
    if missing_value_management < required_mode {
        return Ok(false);
    }

    let Some(code) = missing_code(bit_width, kind)? else {
        return Ok(false);
    };
    Ok(value == code)
}

fn missing_code(bit_width: usize, kind: MissingKind) -> Result<Option<u64>> {
    if bit_width == 0 {
        return Ok(None);
    }
    if bit_width > u64::BITS as usize {
        return Err(Error::UnsupportedPackingWidth(bit_width as u8));
    }

    let max_value = if bit_width == u64::BITS as usize {
        u64::MAX
    } else {
        (1u64 << bit_width) - 1
    };

    let code = match kind {
        MissingKind::Primary => max_value,
        MissingKind::Secondary => max_value.saturating_sub(1),
    };
    Ok(Some(code))
}

fn scale_complex_value<T: DecodeSample>(
    reference: f64,
    binary_factor: f64,
    decimal_factor: f64,
    value: Option<i64>,
) -> T {
    match value {
        Some(value) => T::from_f64(scale_decoded_value(
            reference,
            value as f64,
            binary_factor,
            decimal_factor,
        )),
        None => T::nan(),
    }
}

fn scale_decoded_value(
    reference: f64,
    packed_delta: f64,
    binary_factor: f64,
    decimal_factor: f64,
) -> f64 {
    (reference + packed_delta * binary_factor) * decimal_factor
}

struct GroupReaderLayout {
    reference_offset: usize,
    width_offset: usize,
    length_offset: usize,
    value_offset: usize,
}

impl GroupReaderLayout {
    fn new(start_bit_offset: usize, params: &ComplexPackingParams) -> Result<Self> {
        if params.group_reference_bits as usize > u64::BITS as usize {
            return Err(Error::UnsupportedPackingWidth(params.group_reference_bits));
        }
        if params.group_width_bits as usize > u64::BITS as usize {
            return Err(Error::UnsupportedPackingWidth(params.group_width_bits));
        }
        if params.scaled_group_length_bits as usize > u64::BITS as usize {
            return Err(Error::UnsupportedPackingWidth(
                params.scaled_group_length_bits,
            ));
        }

        let reference_offset = start_bit_offset;
        let width_offset = align_bit_offset(add_group_bits(
            reference_offset,
            params.num_groups,
            params.group_reference_bits as usize,
        )?)?;
        let length_offset = align_bit_offset(add_group_bits(
            width_offset,
            params.num_groups,
            params.group_width_bits as usize,
        )?)?;
        let value_offset = align_bit_offset(add_group_bits(
            length_offset,
            params.num_groups,
            params.scaled_group_length_bits as usize,
        )?)?;

        Ok(Self {
            reference_offset,
            width_offset,
            length_offset,
            value_offset,
        })
    }
}

fn add_group_bits(start_bit_offset: usize, count: usize, bits_per_group: usize) -> Result<usize> {
    count
        .checked_mul(bits_per_group)
        .and_then(|total_bits| start_bit_offset.checked_add(total_bits))
        .ok_or(Error::BitOffsetOverflow)
}

fn align_bit_offset(bit_offset: usize) -> Result<usize> {
    let remainder = bit_offset % 8;
    if remainder == 0 {
        Ok(bit_offset)
    } else {
        bit_offset
            .checked_add(8 - remainder)
            .ok_or(Error::BitOffsetOverflow)
    }
}

struct OutputCursor<'a, T> {
    output: &'a mut [T],
    bitmap: Option<&'a [u8]>,
    next_index: usize,
    values_written: usize,
}

impl<'a, T: DecodeSample> OutputCursor<'a, T> {
    fn new(output: &'a mut [T], bitmap: Option<&'a [u8]>) -> Self {
        Self {
            output,
            bitmap,
            next_index: 0,
            values_written: 0,
        }
    }

    fn push_present(&mut self, value: T) -> Result<()> {
        match self.bitmap {
            Some(bitmap) => {
                if self.advance_to_present(bitmap)? {
                    self.output[self.next_index] = value;
                    self.next_index += 1;
                    self.values_written += 1;
                    return Ok(());
                }

                let expected = count_bitmap_present_points(bitmap, self.output.len())?;
                Err(Error::DataLengthMismatch {
                    expected,
                    actual: self.values_written + 1,
                })
            }
            None => {
                if self.next_index >= self.output.len() {
                    return Err(Error::DataLengthMismatch {
                        expected: self.output.len(),
                        actual: self.next_index + 1,
                    });
                }

                self.output[self.next_index] = value;
                self.next_index += 1;
                self.values_written += 1;
                Ok(())
            }
        }
    }

    fn take_dense(&mut self, count: usize) -> Result<Option<&mut [T]>> {
        if self.bitmap.is_some() {
            return Ok(None);
        }
        let start = self.next_index;
        let end = start.checked_add(count).ok_or(Error::ArithmeticOverflow {
            operation: "advancing dense decoded output",
        })?;
        if end > self.output.len() {
            return Err(Error::DataLengthMismatch {
                expected: self.output.len(),
                actual: end,
            });
        }
        self.next_index = end;
        self.values_written =
            self.values_written
                .checked_add(count)
                .ok_or(Error::ArithmeticOverflow {
                    operation: "counting dense decoded values",
                })?;
        Ok(Some(&mut self.output[start..end]))
    }

    fn advance_to_present(&mut self, bitmap: &[u8]) -> Result<bool> {
        while self.next_index < self.output.len() {
            let byte_index = self.next_index / 8;
            let bit_index = self.next_index % 8;
            let byte = *bitmap.get(byte_index).ok_or(Error::DataLengthMismatch {
                expected: byte_index + 1,
                actual: bitmap.len(),
            })?;
            let remaining = byte & (u8::MAX >> bit_index);
            if remaining != 0 {
                let present_index = byte_index
                    .checked_mul(8)
                    .and_then(|base| base.checked_add(remaining.leading_zeros() as usize))
                    .ok_or(Error::BitOffsetOverflow)?;
                if present_index < self.output.len() {
                    self.output[self.next_index..present_index].fill(T::nan());
                    self.next_index = present_index;
                    return Ok(true);
                }
            }

            let next_byte = byte_index
                .checked_add(1)
                .and_then(|index| index.checked_mul(8))
                .ok_or(Error::BitOffsetOverflow)?
                .min(self.output.len());
            self.output[self.next_index..next_byte].fill(T::nan());
            self.next_index = next_byte;
        }
        Ok(false)
    }

    fn finish(mut self) -> Result<()> {
        if let Some(bitmap) = self.bitmap {
            if self.advance_to_present(bitmap)? {
                return Err(Error::DataLengthMismatch {
                    expected: count_bitmap_present_points(bitmap, self.output.len())?,
                    actual: self.values_written,
                });
            }
        }

        if self.next_index != self.output.len() {
            return Err(Error::DataLengthMismatch {
                expected: self.output.len(),
                actual: self.next_index,
            });
        }

        Ok(())
    }
    fn values_written(&self) -> usize {
        self.values_written
    }
}

#[derive(Debug, Clone)]
struct SpatialDescriptors {
    order: u8,
    first_value: i64,
    second_value: Option<i64>,
    overall_minimum: i64,
}

#[derive(Debug, Clone)]
struct SpatialRestoreState {
    descriptors: SpatialDescriptors,
    previous: Option<i64>,
    previous_difference: Option<i64>,
    non_missing_seen: usize,
}

impl SpatialRestoreState {
    fn new(descriptors: SpatialDescriptors) -> Self {
        Self {
            descriptors,
            previous: None,
            previous_difference: None,
            non_missing_seen: 0,
        }
    }

    fn restore(&mut self, value: Option<i64>) -> Result<Option<i64>> {
        let Some(value) = value else {
            return Ok(None);
        };

        let restored = match self.descriptors.order {
            1 => self.restore_first_order(value)?,
            2 => self.restore_second_order(value)?,
            other => return Err(Error::UnsupportedSpatialDifferencingOrder(other)),
        };

        self.previous = Some(restored);
        self.non_missing_seen += 1;
        Ok(Some(restored))
    }

    fn finish(self) -> Result<()> {
        let expected = match self.descriptors.order {
            1 => 1,
            2 => 2,
            other => return Err(Error::UnsupportedSpatialDifferencingOrder(other)),
        };
        if self.non_missing_seen < expected {
            return Err(Error::DataLengthMismatch {
                expected,
                actual: self.non_missing_seen,
            });
        }
        Ok(())
    }

    fn restore_first_order(&mut self, value: i64) -> Result<i64> {
        if self.non_missing_seen == 0 {
            return Ok(self.descriptors.first_value);
        }

        let delta = value.checked_add(self.descriptors.overall_minimum).ok_or(
            Error::ArithmeticOverflow {
                operation: "restoring first-order spatial difference",
            },
        )?;
        self.previous
            .and_then(|previous| previous.checked_add(delta))
            .ok_or(Error::ArithmeticOverflow {
                operation: "restoring first-order spatial difference",
            })
    }

    fn restore_second_order(&mut self, value: i64) -> Result<i64> {
        match self.non_missing_seen {
            0 => Ok(self.descriptors.first_value),
            1 => {
                let second_value = self.descriptors.second_value.ok_or(Error::InvalidSection {
                    section: 5,
                    reason: "missing second-order spatial differencing descriptors".into(),
                })?;
                self.previous_difference = second_value.checked_sub(self.descriptors.first_value);
                self.previous_difference.ok_or(Error::ArithmeticOverflow {
                    operation: "initializing second-order spatial difference",
                })?;
                Ok(second_value)
            }
            _ => {
                let second_difference = value.checked_add(self.descriptors.overall_minimum).ok_or(
                    Error::ArithmeticOverflow {
                        operation: "restoring second-order spatial difference",
                    },
                )?;
                let difference = self
                    .previous_difference
                    .and_then(|previous_difference| {
                        previous_difference.checked_add(second_difference)
                    })
                    .ok_or(Error::ArithmeticOverflow {
                        operation: "restoring second-order spatial difference",
                    })?;
                let next = self
                    .previous
                    .and_then(|previous| previous.checked_add(difference))
                    .ok_or(Error::ArithmeticOverflow {
                        operation: "restoring second-order spatial difference",
                    })?;
                self.previous_difference = Some(difference);
                Ok(next)
            }
        }
    }
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum MissingKind {
    Primary,
    Secondary,
}

#[cfg(test)]
mod tests {
    use super::{
        bitmap_payload, count_bitmap_present_points, decode_field, decode_payload, unpack_complex,
        unpack_simple, CcsdsBlockSize, CcsdsFlags, CcsdsPackingParams, ComplexPackingParams,
        DataRepresentation, IeeePackingParams, IeeePrecision, PngPackingParams,
        RunLengthPackingParams, ScaledPackingParams, SpatialDifferencingParams,
    };
    use crate::error::Error;
    use grib_core::binary::encode_ieee_binary128;
    use grib_core::bit::BitWriter;

    #[test]
    fn unpack_simple_constant() {
        let params = ScaledPackingParams {
            encoded_values: 5,
            reference_value: 42.0,
            binary_scale: 0,
            decimal_scale: 0,
            bits_per_value: 0,
            original_field_type: 0,
        };
        let values = unpack_simple(&[], &params, 5).unwrap();
        assert_eq!(values, vec![42.0; 5]);
    }

    #[test]
    fn unpack_simple_basic() {
        let params = ScaledPackingParams {
            encoded_values: 5,
            reference_value: 0.0,
            binary_scale: 0,
            decimal_scale: 0,
            bits_per_value: 8,
            original_field_type: 0,
        };
        let values = unpack_simple(&[0, 1, 2, 3, 4], &params, 5).unwrap();
        assert_eq!(values, vec![0.0, 1.0, 2.0, 3.0, 4.0]);
    }

    #[test]
    fn unpack_simple_dense_matches_byte_aligned_and_bit_packed_values() {
        let cases: &[(usize, &[u64])] = &[
            (7, &[1, 2, 3, 127]),
            (8, &[1, 17, 255]),
            (16, &[1, 0x1234, u16::MAX as u64]),
            (32, &[1, 0x1234_5678, u32::MAX as u64]),
            (64, &[1, 1 << 40, u32::MAX as u64 + 1]),
        ];

        for &(bits, packed_values) in cases {
            let mut packed = BitWriter::new();
            for &value in packed_values {
                packed.write(value, bits).unwrap();
            }
            let params = ScaledPackingParams {
                encoded_values: packed_values.len(),
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: bits as u8,
                original_field_type: 0,
            };

            let decoded = unpack_simple(packed.as_bytes(), &params, packed_values.len()).unwrap();
            let expected: Vec<_> = packed_values.iter().map(|&value| value as f64).collect();
            assert_eq!(decoded, expected, "bits_per_value={bits}");
        }
    }

    #[test]
    fn unpack_simple_applies_decimal_scale_to_reference_and_values() {
        let params = ScaledPackingParams {
            encoded_values: 2,
            reference_value: 10.0,
            binary_scale: 0,
            decimal_scale: 1,
            bits_per_value: 8,
            original_field_type: 0,
        };

        let values = unpack_simple(&[0, 20], &params, 2).unwrap();
        assert!((values[0] - 1.0).abs() < 1e-12);
        assert!((values[1] - 3.0).abs() < 1e-12);
    }

    #[test]
    fn decodes_all_ieee_packing_precisions() {
        let values = [1.0, -2.5, 65_536.125];

        for precision in [
            IeeePrecision::Single,
            IeeePrecision::Double,
            IeeePrecision::Quadruple,
        ] {
            let payload = ieee_payload(&values, precision);
            let representation = DataRepresentation::IeeePacking(IeeePackingParams {
                encoded_values: values.len(),
                precision,
            });
            assert_eq!(
                decode_payload(&payload, &representation, None, values.len()).unwrap(),
                values
            );
        }
    }

    #[test]
    fn decodes_ieee_values_through_a_bitmap() {
        let values = [1.0, -2.5, 65_536.125];
        let payload = ieee_payload(&values, IeeePrecision::Double);
        let representation = DataRepresentation::IeeePacking(IeeePackingParams {
            encoded_values: values.len(),
            precision: IeeePrecision::Double,
        });

        let decoded = decode_payload(&payload, &representation, Some(&[0b1011_0000]), 4).unwrap();
        assert_eq!(decoded[0], 1.0);
        assert!(decoded[1].is_nan());
        assert_eq!(decoded[2..], values[1..]);
    }

    #[test]
    fn rejects_truncated_and_trailing_ieee_payload_bytes() {
        let representation = DataRepresentation::IeeePacking(IeeePackingParams {
            encoded_values: 1,
            precision: IeeePrecision::Double,
        });

        for payload in [&[0_u8; 7][..], &[0_u8; 9][..]] {
            assert!(matches!(
                decode_payload(payload, &representation, None, 1),
                Err(Error::DataLengthMismatch {
                    expected: 8,
                    actual
                }) if actual == payload.len()
            ));
        }
    }

    #[test]
    fn decodes_run_length_levels_and_multidigit_runs() {
        let params = RunLengthPackingParams {
            encoded_values: 9,
            bits_per_value: 3,
            max_level_value: 2,
            decimal_scale: 1,
            level_values: vec![10, 20],
        };
        // With range 5, digits 2 and 1 encode seven repetitions after the
        // first level value: 1 + 2 + 1*5 = 8 total values.
        let payload = packed_symbols(3, &[1, 5, 4, 2]);
        let representation = DataRepresentation::RunLengthPacking(params);

        assert_eq!(
            decode_payload(&payload, &representation, None, 9).unwrap(),
            [vec![1.0; 8], vec![2.0]].concat()
        );
    }

    #[test]
    fn decodes_run_length_level_zero_as_missing() {
        let params = RunLengthPackingParams {
            encoded_values: 4,
            bits_per_value: 3,
            max_level_value: 2,
            decimal_scale: 0,
            level_values: vec![10, 20],
        };
        let payload = packed_symbols(3, &[0, 5, 1]);
        let representation = DataRepresentation::RunLengthPacking(params);
        let decoded = decode_payload(&payload, &representation, None, 4).unwrap();

        assert!(decoded[..3].iter().all(|value| value.is_nan()));
        assert_eq!(decoded[3], 10.0);
    }

    #[test]
    fn decodes_all_missing_run_length_field_through_bitmap() {
        let params = RunLengthPackingParams {
            encoded_values: 2,
            bits_per_value: 0,
            max_level_value: 0,
            decimal_scale: 0,
            level_values: Vec::new(),
        };
        let representation = DataRepresentation::RunLengthPacking(params);
        let decoded = decode_payload(&[], &representation, Some(&[0b1010_0000]), 4).unwrap();

        assert!(decoded.iter().all(|value| value.is_nan()));
    }

    #[test]
    fn rejects_malformed_run_length_streams() {
        let representation = DataRepresentation::RunLengthPacking(RunLengthPackingParams {
            encoded_values: 2,
            bits_per_value: 8,
            max_level_value: 2,
            decimal_scale: 0,
            level_values: vec![10, 20],
        });

        let error = decode_payload(&[3], &representation, None, 2).unwrap_err();
        assert!(error.to_string().contains("run starts with code 3"));

        let error = decode_payload(&[1], &representation, None, 2).unwrap_err();
        assert!(matches!(
            error,
            Error::DataLengthMismatch {
                expected: 2,
                actual: 1
            }
        ));

        let error = decode_payload(&[1, 255], &representation, None, 2).unwrap_err();
        assert!(matches!(
            error,
            Error::DataLengthMismatch { expected: 2, .. }
        ));

        for (params, reason) in [
            (
                RunLengthPackingParams {
                    encoded_values: 1,
                    bits_per_value: 65,
                    max_level_value: 0,
                    decimal_scale: 0,
                    level_values: Vec::new(),
                },
                "bits per value must not exceed 64",
            ),
            (
                RunLengthPackingParams {
                    encoded_values: 1,
                    bits_per_value: 0,
                    max_level_value: 0,
                    decimal_scale: 0,
                    level_values: vec![10],
                },
                "zero maximum level requires an empty level table",
            ),
        ] {
            let error = decode_payload(&[], &DataRepresentation::RunLengthPacking(params), None, 1)
                .unwrap_err();
            assert!(error.to_string().contains(reason), "{error}");
        }

        let one_value = DataRepresentation::RunLengthPacking(RunLengthPackingParams {
            encoded_values: 1,
            bits_per_value: 3,
            max_level_value: 2,
            decimal_scale: 0,
            level_values: vec![10, 20],
        });
        let error = decode_payload(&[0b0010_0001], &one_value, None, 1).unwrap_err();
        assert!(error.to_string().contains("nonzero bits remain"));
    }

    #[test]
    fn decodes_bitmap_masked_field() {
        let data_section = [0, 0, 0, 8, 7, 10, 20, 30];
        let bitmap_section = [0, 0, 0, 7, 6, 0, 0b1011_0000];
        let representation = DataRepresentation::SimplePacking(ScaledPackingParams {
            encoded_values: 3,
            reference_value: 0.0,
            binary_scale: 0,
            decimal_scale: 0,
            bits_per_value: 8,
            original_field_type: 0,
        });

        let bitmap = bitmap_payload(&bitmap_section).unwrap();
        let decoded = decode_field(&data_section, &representation, bitmap, 4).unwrap();
        assert_eq!(decoded[0], 10.0);
        assert!(decoded[1].is_nan());
        assert_eq!(decoded[2], 20.0);
        assert_eq!(decoded[3], 30.0);
    }

    fn ieee_payload(values: &[f64], precision: IeeePrecision) -> Vec<u8> {
        values
            .iter()
            .flat_map(|value| match precision {
                IeeePrecision::Single => Vec::from((*value as f32).to_be_bytes()),
                IeeePrecision::Double => Vec::from(value.to_be_bytes()),
                IeeePrecision::Quadruple => Vec::from(encode_ieee_binary128(*value)),
                _ => unreachable!("test covers all supported IEEE precisions"),
            })
            .collect()
    }

    fn packed_symbols(bits: usize, symbols: &[u64]) -> Vec<u8> {
        let mut writer = BitWriter::new();
        for symbol in symbols {
            writer.write(*symbol, bits).unwrap();
        }
        writer.align_to_byte().unwrap();
        writer.into_bytes()
    }

    #[test]
    fn decodes_bitmap_across_empty_bytes_and_a_partial_final_byte() {
        let data_section = [0, 0, 0, 7, 7, 11, 22];
        let bitmap_section = [0, 0, 0, 9, 6, 0, 0b0000_0000, 0b0010_0000, 0b1000_0000];
        let representation = DataRepresentation::SimplePacking(ScaledPackingParams {
            encoded_values: 2,
            reference_value: 0.0,
            binary_scale: 0,
            decimal_scale: 0,
            bits_per_value: 8,
            original_field_type: 0,
        });

        let bitmap = bitmap_payload(&bitmap_section).unwrap();
        let decoded = decode_field(&data_section, &representation, bitmap, 19).unwrap();
        assert!(decoded[..10].iter().all(|value| value.is_nan()));
        assert_eq!(decoded[10], 11.0);
        assert!(decoded[11..16].iter().all(|value| value.is_nan()));
        assert_eq!(decoded[16], 22.0);
        assert!(decoded[17..].iter().all(|value| value.is_nan()));
    }

    #[cfg(feature = "ccsds")]
    #[test]
    fn decodes_ccsds_unsigned_samples_across_container_widths() {
        for (bits, values) in [
            (8, vec![0, 1, 17, 255]),
            (16, vec![0, 1, 0x1234, 0xffff]),
            (20, vec![0, 1, 0x54321, 0xfffff]),
            (24, vec![0, 1, 0x654321, 0xffffff]),
            (32, vec![0, 1, 0x76543210, u32::MAX as i64]),
        ] {
            let flags = CcsdsFlags::PREPROCESS;
            let (payload, representation) = encode_ccsds_test_values(&values, bits, flags, 16);
            let decoded = decode_payload(&payload, &representation, None, values.len()).unwrap();
            let expected = values.iter().map(|value| *value as f64).collect::<Vec<_>>();
            assert_eq!(decoded, expected, "bits_per_value={bits}");
        }
    }

    #[cfg(feature = "ccsds")]
    #[test]
    fn decodes_ccsds_signed_preprocessing_as_unsigned_grib_differences() {
        let values = [0, 1, 2047, 2048, 4095];
        let flags = CcsdsFlags::SIGNED | CcsdsFlags::PREPROCESS;
        let (payload, mut representation) = encode_ccsds_test_values(&values, 12, flags, 16);
        let DataRepresentation::CcsdsPacking(params) = &mut representation else {
            unreachable!();
        };
        *params = CcsdsPackingParams::new(
            ScaledPackingParams {
                encoded_values: values.len(),
                reference_value: 10.0,
                binary_scale: 1,
                decimal_scale: 1,
                bits_per_value: 12,
                original_field_type: 0,
            },
            flags,
            CcsdsBlockSize::SIXTEEN,
            128,
        )
        .unwrap();

        let decoded = decode_payload(&payload, &representation, None, values.len()).unwrap();
        let expected = values
            .iter()
            .map(|value| (10.0 + *value as f64 * 2.0) / 10.0)
            .collect::<Vec<_>>();
        assert_float_values(&decoded, &expected);
    }

    #[cfg(feature = "ccsds")]
    #[test]
    fn decodes_ccsds_restricted_stream_with_wire_to_codec_flag_mapping() {
        let values = [0, 1, 2, 3, 15, 7, 4, 0];
        let flags = CcsdsFlags::RESTRICTED | CcsdsFlags::PREPROCESS;
        let (payload, representation) = encode_ccsds_test_values(&values, 4, flags, 16);

        let decoded = decode_payload(&payload, &representation, None, values.len()).unwrap();
        assert_eq!(
            decoded,
            values.iter().map(|value| *value as f64).collect::<Vec<_>>()
        );
        assert!(super::ccsds_decoder_flags(flags).contains(grib_aec::AecFlags::RESTRICTED));
    }

    #[cfg(feature = "ccsds")]
    #[test]
    fn decodes_ccsds_nonstandard_even_block_when_not_enforced() {
        let values = (0..24).collect::<Vec<_>>();
        let flags = CcsdsFlags::PREPROCESS | CcsdsFlags::NOT_ENFORCE;
        let (payload, representation) = encode_ccsds_test_values(&values, 8, flags, 12);

        let decoded = decode_payload(&payload, &representation, None, values.len()).unwrap();
        assert_eq!(
            decoded,
            values.iter().map(|value| *value as f64).collect::<Vec<_>>()
        );
    }

    #[cfg(feature = "ccsds")]
    #[test]
    fn decodes_constant_ccsds_field_with_bitmap() {
        let bitmap_section = [0, 0, 0, 7, 6, 0, 0b1011_0000];
        let representation = DataRepresentation::CcsdsPacking(
            CcsdsPackingParams::new(
                ScaledPackingParams {
                    encoded_values: 3,
                    reference_value: 42.0,
                    binary_scale: 0,
                    decimal_scale: 1,
                    bits_per_value: 0,
                    original_field_type: 0,
                },
                CcsdsFlags::DEFAULT,
                CcsdsBlockSize::THIRTY_TWO,
                128,
            )
            .unwrap(),
        );

        let decoded = decode_payload(
            &[],
            &representation,
            bitmap_payload(&bitmap_section).unwrap(),
            4,
        )
        .unwrap();
        assert_eq!(decoded[0], 42.0);
        assert!(decoded[1].is_nan());
        assert_eq!(decoded[2], 42.0);
        assert_eq!(decoded[3], 42.0);
    }

    #[cfg(not(feature = "ccsds"))]
    #[test]
    fn ccsds_packing_requires_ccsds_feature() {
        let representation = DataRepresentation::CcsdsPacking(
            CcsdsPackingParams::new(
                ScaledPackingParams {
                    encoded_values: 1,
                    reference_value: 0.0,
                    binary_scale: 0,
                    decimal_scale: 0,
                    bits_per_value: 8,
                    original_field_type: 0,
                },
                CcsdsFlags::DEFAULT,
                CcsdsBlockSize::THIRTY_TWO,
                128,
            )
            .unwrap(),
        );

        let error = decode_payload(&[], &representation, None, 1).unwrap_err();
        assert!(matches!(error, Error::UnsupportedDataTemplate(42)));
    }

    #[cfg(feature = "png")]
    #[test]
    fn decodes_png_packing_grayscale8() {
        let payload = encode_png(
            2,
            2,
            png::ColorType::Grayscale,
            png::BitDepth::Eight,
            &[1, 2, 3, 4],
        );
        let representation = DataRepresentation::PngPacking(PngPackingParams {
            packing: ScaledPackingParams {
                encoded_values: 4,
                reference_value: 10.0,
                binary_scale: 1,
                decimal_scale: 1,
                bits_per_value: 8,
                original_field_type: 0,
            },
        });

        let decoded = decode_payload(&payload, &representation, None, 4).unwrap();
        assert_float_values(&decoded, &[1.2, 1.4, 1.6, 1.8]);
    }

    #[cfg(feature = "png")]
    #[test]
    fn decodes_png_packing_grayscale16() {
        let payload = encode_png(
            2,
            1,
            png::ColorType::Grayscale,
            png::BitDepth::Sixteen,
            &[0x01, 0x00, 0x02, 0x00],
        );
        let representation = DataRepresentation::PngPacking(PngPackingParams {
            packing: ScaledPackingParams {
                encoded_values: 2,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 16,
                original_field_type: 0,
            },
        });

        let decoded = decode_payload(&payload, &representation, None, 2).unwrap();
        assert_eq!(decoded, vec![256.0, 512.0]);
    }

    #[cfg(feature = "png")]
    #[test]
    fn decodes_png_packing_grayscale_subbyte() {
        let payload = encode_png(
            4,
            1,
            png::ColorType::Grayscale,
            png::BitDepth::Four,
            &[0x12, 0x34],
        );
        let representation = DataRepresentation::PngPacking(PngPackingParams {
            packing: ScaledPackingParams {
                encoded_values: 4,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 4,
                original_field_type: 0,
            },
        });

        let decoded = decode_payload(&payload, &representation, None, 4).unwrap();
        assert_eq!(decoded, vec![1.0, 2.0, 3.0, 4.0]);
    }

    #[cfg(feature = "png")]
    #[test]
    fn decodes_png_packing_rgb8() {
        let payload = encode_png(
            2,
            1,
            png::ColorType::Rgb,
            png::BitDepth::Eight,
            &[0x01, 0x02, 0x03, 0x04, 0x05, 0x06],
        );
        let representation = DataRepresentation::PngPacking(PngPackingParams {
            packing: ScaledPackingParams {
                encoded_values: 2,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 24,
                original_field_type: 0,
            },
        });

        let decoded = decode_payload(&payload, &representation, None, 2).unwrap();
        assert_eq!(decoded, vec![0x01_02_03 as f64, 0x04_05_06 as f64]);
    }

    #[cfg(feature = "png")]
    #[test]
    fn decodes_png_packing_rgba8() {
        let payload = encode_png(
            1,
            1,
            png::ColorType::Rgba,
            png::BitDepth::Eight,
            &[0x01, 0x02, 0x03, 0x04],
        );
        let representation = DataRepresentation::PngPacking(PngPackingParams {
            packing: ScaledPackingParams {
                encoded_values: 1,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 32,
                original_field_type: 0,
            },
        });

        let decoded = decode_payload(&payload, &representation, None, 1).unwrap();
        assert_eq!(decoded, vec![0x01_02_03_04 as f64]);
    }

    #[cfg(feature = "png")]
    #[test]
    fn decodes_png_packing_with_bitmap() {
        let payload = encode_png(
            3,
            1,
            png::ColorType::Grayscale,
            png::BitDepth::Eight,
            &[10, 20, 30],
        );
        let bitmap_section = [0, 0, 0, 7, 6, 0, 0b1011_0000];
        let representation = DataRepresentation::PngPacking(PngPackingParams {
            packing: ScaledPackingParams {
                encoded_values: 3,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 8,
                original_field_type: 0,
            },
        });

        let decoded = decode_payload(
            &payload,
            &representation,
            bitmap_payload(&bitmap_section).unwrap(),
            4,
        )
        .unwrap();
        assert_eq!(decoded[0], 10.0);
        assert!(decoded[1].is_nan());
        assert_eq!(decoded[2], 20.0);
        assert_eq!(decoded[3], 30.0);
    }

    #[cfg(feature = "png")]
    #[test]
    fn rejects_png_header_dimension_mismatch_before_allocating_frame_buffer() {
        let mut payload = encode_png(1, 1, png::ColorType::Grayscale, png::BitDepth::Eight, &[1]);
        payload[16..20].copy_from_slice(&100_000u32.to_be_bytes());
        payload[20..24].copy_from_slice(&100_000u32.to_be_bytes());
        let crc = test_crc32(&payload[12..29]);
        payload[29..33].copy_from_slice(&crc.to_be_bytes());

        let representation = DataRepresentation::PngPacking(PngPackingParams {
            packing: ScaledPackingParams {
                encoded_values: 1,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 8,
                original_field_type: 0,
            },
        });

        let err = decode_payload(&payload, &representation, None, 1).unwrap_err();
        assert!(matches!(
            err,
            Error::DataLengthMismatch {
                expected: 1,
                actual: 10_000_000_000,
            }
        ));
    }

    #[cfg(feature = "jpeg2000")]
    #[test]
    fn rejects_jpeg2000_header_dimension_mismatch_before_decoding_codestream() {
        let mut payload = vec![0u8; 45];
        payload[..4].copy_from_slice(&[0xff, 0x4f, 0xff, 0x51]);
        payload[4..6].copy_from_slice(&41u16.to_be_bytes());
        payload[8..12].copy_from_slice(&100_000u32.to_be_bytes());
        payload[12..16].copy_from_slice(&100_000u32.to_be_bytes());
        payload[40..42].copy_from_slice(&1u16.to_be_bytes());
        payload[42] = 7;
        payload[43] = 1;
        payload[44] = 1;

        let representation = DataRepresentation::Jpeg2000Packing(super::Jpeg2000PackingParams {
            packing: ScaledPackingParams {
                encoded_values: 1,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 8,
                original_field_type: 0,
            },
            compression_type: 0,
            target_compression_ratio: 0,
        });

        let err = decode_payload(&payload, &representation, None, 1).unwrap_err();
        assert!(matches!(
            err,
            Error::DataLengthMismatch {
                expected: 1,
                actual: 10_000_000_000,
            }
        ));
    }

    #[cfg(feature = "jpeg2000")]
    #[test]
    fn reports_negative_unsigned_jpeg2000_samples_as_out_of_range() {
        let error = super::jpeg2000_raw_value(-1, 8).unwrap_err();
        assert!(matches!(error, Error::ValueOutOfRange(_)));
    }

    #[cfg(feature = "png")]
    #[test]
    fn reports_incompatible_png_layout_as_invalid_section_seven() {
        let payload = encode_png(1, 1, png::ColorType::Rgb, png::BitDepth::Eight, &[1, 2, 3]);
        let representation = DataRepresentation::PngPacking(PngPackingParams {
            packing: ScaledPackingParams {
                encoded_values: 1,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 8,
                original_field_type: 0,
            },
        });

        let error = decode_payload(&payload, &representation, None, 1).unwrap_err();
        assert!(matches!(error, Error::InvalidSection { section: 7, .. }));
    }

    #[cfg(not(feature = "png"))]
    #[test]
    fn png_packing_requires_png_feature() {
        let representation = DataRepresentation::PngPacking(PngPackingParams {
            packing: ScaledPackingParams {
                encoded_values: 1,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 8,
                original_field_type: 0,
            },
        });

        let err = decode_payload(&[], &representation, None, 1).unwrap_err();
        assert!(matches!(err, Error::UnsupportedDataTemplate(41)));
    }

    #[cfg(not(feature = "jpeg2000"))]
    #[test]
    fn jpeg2000_packing_requires_jpeg2000_feature() {
        let representation = DataRepresentation::Jpeg2000Packing(super::Jpeg2000PackingParams {
            packing: ScaledPackingParams {
                encoded_values: 1,
                reference_value: 0.0,
                binary_scale: 0,
                decimal_scale: 0,
                bits_per_value: 8,
                original_field_type: 0,
            },
            compression_type: 1,
            target_compression_ratio: 255,
        });

        let err = decode_payload(&[], &representation, None, 1).unwrap_err();
        assert!(matches!(err, Error::UnsupportedDataTemplate(40)));
    }

    #[test]
    fn rejects_simple_packing_wider_than_u64() {
        let params = ScaledPackingParams {
            encoded_values: 1,
            reference_value: 0.0,
            binary_scale: 0,
            decimal_scale: 0,
            bits_per_value: 65,
            original_field_type: 0,
        };
        let err = unpack_simple(&[0; 9], &params, 1).unwrap_err();
        assert!(matches!(err, Error::UnsupportedPackingWidth(65)));
    }

    #[test]
    fn rejects_encoded_value_count_mismatch_without_bitmap() {
        let data_section = [0, 0, 0, 8, 7, 10, 20, 30];
        let representation = DataRepresentation::SimplePacking(ScaledPackingParams {
            encoded_values: 3,
            reference_value: 0.0,
            binary_scale: 0,
            decimal_scale: 0,
            bits_per_value: 8,
            original_field_type: 0,
        });

        let err = decode_field(&data_section, &representation, None, 4).unwrap_err();
        assert!(matches!(
            err,
            Error::DataLengthMismatch {
                expected: 4,
                actual: 3,
            }
        ));
    }

    #[test]
    fn rejects_bitmap_present_count_mismatch() {
        let data_section = [0, 0, 0, 7, 7, 10, 20];
        let bitmap_section = [0, 0, 0, 7, 6, 0, 0b1011_0000];
        let representation = DataRepresentation::SimplePacking(ScaledPackingParams {
            encoded_values: 2,
            reference_value: 0.0,
            binary_scale: 0,
            decimal_scale: 0,
            bits_per_value: 8,
            original_field_type: 0,
        });

        let bitmap = bitmap_payload(&bitmap_section).unwrap();
        let err = decode_field(&data_section, &representation, bitmap, 4).unwrap_err();
        assert!(matches!(
            err,
            Error::DataLengthMismatch {
                expected: 3,
                actual: 2,
            }
        ));
    }

    #[test]
    fn counts_bitmap_present_points_with_partial_bytes() {
        let present = count_bitmap_present_points(&[0b1011_1111], 3).unwrap();
        assert_eq!(present, 2);
    }

    #[cfg(feature = "ccsds")]
    fn encode_ccsds_test_values(
        values: &[i64],
        bits_per_value: u8,
        flags: CcsdsFlags,
        block_size: u8,
    ) -> (Vec<u8>, DataRepresentation) {
        let bytes_per_sample = usize::from(bits_per_value).div_ceil(8);
        let mask = if bits_per_value == 32 {
            u32::MAX
        } else {
            (1u32 << bits_per_value) - 1
        };
        let mut unpacked = Vec::with_capacity(values.len() * bytes_per_sample);
        for &value in values {
            let raw = (value as u32) & mask;
            let bytes = raw.to_be_bytes();
            unpacked.extend_from_slice(&bytes[4 - bytes_per_sample..]);
        }

        let codec_params = grib_aec::AecParams::new(
            bits_per_value,
            block_size,
            128,
            super::ccsds_decoder_flags(flags),
        )
        .unwrap();
        let payload = grib_aec::encode(&unpacked, codec_params).unwrap();
        let representation = DataRepresentation::CcsdsPacking(
            CcsdsPackingParams::new(
                ScaledPackingParams {
                    encoded_values: values.len(),
                    reference_value: 0.0,
                    binary_scale: 0,
                    decimal_scale: 0,
                    bits_per_value,
                    original_field_type: 0,
                },
                flags,
                CcsdsBlockSize::new(block_size).unwrap(),
                128,
            )
            .unwrap(),
        );
        (payload, representation)
    }

    #[test]
    fn unpacks_complex_packing_with_explicit_missing() {
        let params = ComplexPackingParams {
            encoded_values: 4,
            reference_value: 0.0,
            binary_scale: 0,
            decimal_scale: 0,
            group_reference_bits: 4,
            original_field_type: 0,
            group_splitting_method: 1,
            missing_value_management: 1,
            primary_missing_substitute: 0x6655_4433,
            secondary_missing_substitute: 0x2211_00ff,
            num_groups: 2,
            group_width_reference: 0,
            group_width_bits: 2,
            group_length_reference: 2,
            group_length_increment: 1,
            true_length_last_group: 2,
            scaled_group_length_bits: 0,
            spatial_differencing: None,
        };

        assert_eq!(params.primary_missing_substitute(), Some(0x6655_4433));
        assert_eq!(params.secondary_missing_substitute(), None);

        let values = unpack_complex(&[0x79, 0x90, 0x34], &params).unwrap();
        assert_eq!(values[0], 7.0);
        assert!(values[1].is_nan());
        assert_eq!(values[2], 9.0);
        assert!(values[3].is_nan());
    }

    #[test]
    fn unpacks_complex_packing_applies_decimal_scale_to_reference_and_values() {
        let params = ComplexPackingParams {
            encoded_values: 2,
            reference_value: 10.0,
            binary_scale: 0,
            decimal_scale: 1,
            group_reference_bits: 4,
            original_field_type: 0,
            group_splitting_method: 1,
            missing_value_management: 0,
            primary_missing_substitute: u32::MAX,
            secondary_missing_substitute: u32::MAX,
            num_groups: 1,
            group_width_reference: 4,
            group_width_bits: 0,
            group_length_reference: 2,
            group_length_increment: 1,
            true_length_last_group: 2,
            scaled_group_length_bits: 0,
            spatial_differencing: None,
        };

        let values = unpack_complex(&[0x10, 0x02], &params).unwrap();
        assert!((values[0] - 1.1).abs() < 1e-12);
        assert!((values[1] - 1.3).abs() < 1e-12);
    }

    #[test]
    fn unpacks_complex_packing_with_second_order_spatial_differencing() {
        let params = ComplexPackingParams {
            encoded_values: 4,
            reference_value: 0.0,
            binary_scale: 0,
            decimal_scale: 0,
            group_reference_bits: 1,
            original_field_type: 0,
            group_splitting_method: 1,
            missing_value_management: 0,
            primary_missing_substitute: u32::MAX,
            secondary_missing_substitute: u32::MAX,
            num_groups: 1,
            group_width_reference: 1,
            group_width_bits: 0,
            group_length_reference: 4,
            group_length_increment: 1,
            true_length_last_group: 4,
            scaled_group_length_bits: 0,
            spatial_differencing: Some(SpatialDifferencingParams {
                order: 2,
                descriptor_octets: 2,
            }),
        };

        let values =
            unpack_complex(&[0x00, 0x0A, 0x00, 0x0D, 0x00, 0x03, 0x00, 0x10], &params).unwrap();
        assert_eq!(values, vec![10.0, 13.0, 19.0, 29.0]);
    }

    #[test]
    fn unpacks_complex_packing_with_spatial_differencing_and_missing_values() {
        let params = ComplexPackingParams {
            encoded_values: 4,
            reference_value: 0.0,
            binary_scale: 0,
            decimal_scale: 0,
            group_reference_bits: 1,
            original_field_type: 0,
            group_splitting_method: 1,
            missing_value_management: 1,
            primary_missing_substitute: u32::MAX,
            secondary_missing_substitute: u32::MAX,
            num_groups: 1,
            group_width_reference: 2,
            group_width_bits: 0,
            group_length_reference: 4,
            group_length_increment: 1,
            true_length_last_group: 4,
            scaled_group_length_bits: 0,
            spatial_differencing: Some(SpatialDifferencingParams {
                order: 1,
                descriptor_octets: 2,
            }),
        };

        let values = unpack_complex(&[0x00, 0x0A, 0x00, 0x03, 0x00, 0x32], &params).unwrap();
        assert_eq!(values[0], 10.0);
        assert!(values[1].is_nan());
        assert_eq!(values[2], 13.0);
        assert_eq!(values[3], 18.0);
    }

    #[cfg(feature = "png")]
    fn encode_png(
        width: u32,
        height: u32,
        color_type: png::ColorType,
        bit_depth: png::BitDepth,
        data: &[u8],
    ) -> Vec<u8> {
        let mut payload = Vec::new();
        {
            let mut encoder = png::Encoder::new(&mut payload, width, height);
            encoder.set_color(color_type);
            encoder.set_depth(bit_depth);
            let mut writer = encoder.write_header().unwrap();
            writer.write_image_data(data).unwrap();
        }
        payload
    }

    #[cfg(feature = "png")]
    fn test_crc32(bytes: &[u8]) -> u32 {
        let mut crc = u32::MAX;
        for &byte in bytes {
            crc ^= u32::from(byte);
            for _ in 0..8 {
                let mask = 0u32.wrapping_sub(crc & 1);
                crc = (crc >> 1) ^ (0xedb8_8320 & mask);
            }
        }
        !crc
    }

    #[cfg(feature = "png")]
    fn assert_float_values(actual: &[f64], expected: &[f64]) {
        assert_eq!(actual.len(), expected.len());
        for (actual, expected) in actual.iter().zip(expected) {
            assert!(
                (actual - expected).abs() < 1e-12,
                "expected {expected}, got {actual}"
            );
        }
    }
}