fits-well 0.1.4

use super::*;
use crate::reader::{FitsReader, StreamReader};
use std::fs::File;

fn open(name: &str) -> StreamReader<File> {
    FitsReader::open(File::open(format!("tests/data/fits/{name}")).unwrap()).unwrap()
}

/// The fixtures encode value(x, y) = x*7 − y*5 over a 24×16 i16 image.
fn expect_pixel(flat: usize) -> i16 {
    let (x, y) = (flat % 24, flat / 24);
    (x as i16) * 7 - (y as i16) * 5
}

fn check_decoded(name: &str) {
    let mut f = open(name);
    let img = f.read_image(1).unwrap();
    assert_eq!(img.shape, vec![24, 16]);
    match img.decode() {
        ImageData::I16(v) => {
            assert_eq!(v.len(), 24 * 16);
            for (i, &got) in v.iter().enumerate() {
                assert_eq!(got, expect_pixel(i), "pixel {i} of {name}");
            }
        }
        other => panic!("expected I16, got {other:?}"),
    }
}

#[test]
fn decompresses_gzip_1_tiled_image() {
    check_decoded("comp_gzip_i16.fits");
}

#[test]
fn decompresses_rice_1_tiled_image() {
    check_decoded("comp_rice_i16.fits");
}

#[test]
fn decompresses_hcompress_1_tiled_image() {
    // Lossless HCOMPRESS (SCALE=0), single 24×16 tile.
    check_decoded("comp_hcomp_i16.fits");
}

/// Decode an i32 image and compare pixel-exact against astropy's reconstruction
/// stored as a plain-image reference.
fn check_i32_against_ref(compressed: &str, reference: &str) {
    let got = match open(compressed).read_image(1).unwrap().decode() {
        ImageData::I32(v) => v,
        other => panic!("expected I32, got {other:?}"),
    };
    let want = match open(reference).read_image(0).unwrap().decode() {
        ImageData::I32(v) => v,
        other => panic!("expected I32 reference, got {other:?}"),
    };
    assert_eq!(got, want, "{compressed} must match astropy {reference}");
}

#[test]
fn decompresses_hcompress_lossy() {
    // Lossy HCOMPRESS (SCALE=4, SMOOTH=0): exercises undigitize (×scale).
    check_i32_against_ref("comp_hcomp_lossy.fits", "comp_ref_hcomp_lossy.fits");
}

#[test]
fn decompresses_hcompress_smoothed() {
    // SMOOTH=1: the SMOOTH ZVAL triggers inverse-transform smoothing, which must
    // reproduce astropy's smoothed reconstruction bit-for-bit.
    check_i32_against_ref("comp_hcomp_smooth.fits", "comp_ref_hcomp_smooth.fits");
}

#[test]
fn decompresses_subtractive_dither_2() {
    // SUBTRACTIVE_DITHER_2 float: must match astropy's dithered reconstruction.
    check_float("comp_dither2_f32.fits", "comp_ref_dither2_f32.fits");
}

#[test]
fn decompresses_float_with_nan_nulls() {
    // SUBTRACTIVE_DITHER_1 with ZBLANK: null pixels decode to NaN, the rest match.
    let got = match open("comp_nan_f32.fits").read_image(1).unwrap().decode() {
        ImageData::F32(v) => v,
        other => panic!("expected F32, got {other:?}"),
    };
    let want = match open("comp_ref_nan_f32.fits")
        .read_image(0)
        .unwrap()
        .decode()
    {
        ImageData::F32(v) => v,
        other => panic!("expected F32 reference, got {other:?}"),
    };
    assert_eq!(got.len(), want.len());
    let mut nan_count = 0;
    for (i, (&g, &w)) in got.iter().zip(&want).enumerate() {
        if w.is_nan() {
            assert!(g.is_nan(), "pixel {i} should be NaN");
            nan_count += 1;
        } else {
            assert_eq!(g, w, "pixel {i}");
        }
    }
    assert_eq!(nan_count, 2, "expected 2 null pixels");
}

/// Emit compressed files written by this crate for external (astropy) validation.
/// Run with `cargo test --features compression -- --ignored emit_`.
#[test]
#[ignore]
fn emit_compressed_files_for_astropy() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::fs::File;

    let samples: Vec<i16> = (0..24 * 16)
        .map(|i| (i % 24) as i16 * 7 - (i / 24) as i16 * 5)
        .collect();
    let image = Image {
        shape: vec![24, 16],
        samples: ImageData::I16(samples),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    for (cmptype, tiles) in [
        ("GZIP_1", &[][..]),
        ("GZIP_2", &[]),
        ("RICE_1", &[]),
        ("HCOMPRESS_1", &[24, 16]),
    ] {
        let f = File::create(format!(".tmp/wr_{}.fits", cmptype.to_lowercase())).unwrap();
        let mut w = FitsWriter::new(f);
        w.write_compressed_image(&image, cmptype, &CompressOptions::tiled(tiles))
            .unwrap();
    }

    // PLIO needs a non-negative mask image.
    let mask: Vec<i32> = (0..24 * 16).map(|i| (i % 24 + i / 24) % 7).collect();
    let mask_image = Image {
        shape: vec![24, 16],
        samples: ImageData::I32(mask),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    let f = File::create(".tmp/wr_plio_1.fits").unwrap();
    let mut w = FitsWriter::new(f);
    w.write_compressed_image(&mask_image, "PLIO_1", &CompressOptions::default())
        .unwrap();

    // Quantized float (SUBTRACTIVE_DITHER_1) for astropy to reconstruct.
    let fimage = Image {
        shape: vec![24, 16],
        samples: ImageData::F32(float_field()),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    let f = File::create(".tmp/wr_ricef.fits").unwrap();
    let mut w = FitsWriter::new(f);
    w.write_compressed_image(&fimage, "RICE_1", &CompressOptions::tiled([24, 16]))
        .unwrap();
}

#[test]
fn compression_write_round_trips_through_decode() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;

    let samples: Vec<i16> = (0..24 * 16)
        .map(|i| (i % 24) as i16 * 7 - (i / 24) as i16 * 5)
        .collect();
    let image = Image {
        shape: vec![24, 16],
        samples: ImageData::I16(samples.clone()),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    // Row tiling for the byte codecs; HCOMPRESS needs a genuinely 2-D tile.
    for (cmptype, tiles) in [
        ("GZIP_1", &[][..]),
        ("GZIP_2", &[]),
        ("RICE_1", &[]),
        ("HCOMPRESS_1", &[24, 16]),
    ] {
        let mut w = FitsWriter::new(Cursor::new(Vec::new()));
        w.write_compressed_image(&image, cmptype, &CompressOptions::tiled(tiles))
            .unwrap();
        let mut r = FitsReader::open(Cursor::new(w.into_inner().into_inner())).unwrap();
        let back = r.read_image(1).unwrap();
        assert_eq!(back.shape, vec![24, 16], "{cmptype}");
        match back.decode() {
            ImageData::I16(v) => assert_eq!(v, samples, "{cmptype} round-trip"),
            other => panic!("{cmptype}: expected I16, got {other:?}"),
        }
    }
}

/// A 24×16 float field: a smooth ramp plus genuine high-frequency noise (a
/// splitmix64 hash, decorrelated neighbour-to-neighbour) so the 3rd-order MAD
/// estimate is realistic (≈ 1) and the tile genuinely quantizes.
fn float_field() -> Vec<f32> {
    let mix = |i: u64| {
        // splitmix64 finalizer — uncorrelated output for consecutive inputs.
        let mut z = i.wrapping_add(0x9E37_79B9_7F4A_7C15);
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
        z ^ (z >> 31)
    };
    (0..24 * 16)
        .map(|i| {
            let (x, y) = (i % 24, i / 24);
            let smooth = 100.0 + 3.0 * x as f32 - 2.0 * y as f32;
            let noise = (mix(i as u64) % 2000) as f32 / 1000.0 - 1.0; // ±1.0
            smooth + noise
        })
        .collect()
}

#[test]
fn float_quantize_write_round_trips_within_tolerance() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;

    let orig = float_field();
    let image = Image {
        shape: vec![24, 16],
        samples: ImageData::F32(orig.clone()),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    for cmptype in ["RICE_1", "GZIP_1", "GZIP_2"] {
        let mut w = FitsWriter::new(Cursor::new(Vec::new()));
        // Whole-image tile so the noise estimate sees the full field.
        w.write_compressed_image(&image, cmptype, &CompressOptions::tiled([24, 16]))
            .unwrap();
        let mut r = FitsReader::open(Cursor::new(w.into_inner().into_inner())).unwrap();
        let back = match r.read_image(1).unwrap().decode() {
            ImageData::F32(v) => v,
            other => panic!("{cmptype}: expected F32, got {other:?}"),
        };
        assert_eq!(back.len(), orig.len(), "{cmptype}");
        // Quantization error is bounded by ~0.5·delta; delta ≈ noise/4 ≈ 0.07 for
        // this field, so 0.2 is a safe ceiling. Also confirm it actually quantized.
        let max_err = orig
            .iter()
            .zip(&back)
            .map(|(a, b)| (a - b).abs())
            .fold(0.0f32, f32::max);
        assert!(max_err < 0.2, "{cmptype} max error {max_err} too large");
        assert!(
            orig.iter().zip(&back).any(|(a, b)| a != b),
            "{cmptype} stored losslessly — quantized path not exercised"
        );
    }
}

#[test]
fn dither_option_sets_zquantiz_and_round_trips() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;

    let image = Image {
        shape: vec![24, 16],
        samples: ImageData::F32(float_field()),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    // Each `DitherMethod` writes its `ZQUANTIZ` keyword and round-trips; the option is
    // honored rather than always emitting the hardcoded SUBTRACTIVE_DITHER_1.
    for (dither, zquantiz) in [
        (DitherMethod::None, "NO_DITHER"),
        (DitherMethod::Subtractive1, "SUBTRACTIVE_DITHER_1"),
        (DitherMethod::Subtractive2, "SUBTRACTIVE_DITHER_2"),
    ] {
        let mut w = FitsWriter::new(Cursor::new(Vec::new()));
        w.write_compressed_image(
            &image,
            "RICE_1",
            &CompressOptions {
                dither,
                ..CompressOptions::tiled([24, 16])
            },
        )
        .unwrap();
        let mut r = FitsReader::open(Cursor::new(w.into_inner().into_inner())).unwrap();
        assert_eq!(
            r.hdus()[1].header.get_text("ZQUANTIZ"),
            Some(zquantiz),
            "{dither:?} must write {zquantiz}"
        );
        match r.read_image(1).unwrap().decode() {
            ImageData::F32(v) => assert_eq!(v.len(), 24 * 16, "{dither:?}"),
            other => panic!("{dither:?}: expected F32, got {other:?}"),
        }
    }
}

#[test]
fn tile_shape_with_wrong_rank_is_rejected() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::error::FitsError;
    use crate::writer::FitsWriter;
    use std::io::Cursor;

    // A 2-D image with a 3-axis tile shape: a wrong-rank request errors rather than
    // silently row-tiling (an empty shape is still the row-tiling default).
    let image = Image {
        shape: vec![4, 3],
        samples: ImageData::I16((0..12).map(|i| i as i16).collect()),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    let mut w = FitsWriter::new(Cursor::new(Vec::new()));
    let err = w.write_compressed_image(&image, "RICE_1", &CompressOptions::tiled([2, 2, 2]));
    assert!(
        matches!(
            err,
            Err(FitsError::TileShapeRankMismatch {
                tile_rank: 3,
                image_rank: 2,
            })
        ),
        "got {err:?}"
    );
}

#[test]
fn float_write_preserves_nan_nulls() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;

    let mut orig = float_field();
    orig[5 + 3 * 24] = f32::NAN;
    orig[20 + 10 * 24] = f32::NAN;
    let image = Image {
        shape: vec![24, 16],
        samples: ImageData::F32(orig.clone()),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    let mut w = FitsWriter::new(Cursor::new(Vec::new()));
    w.write_compressed_image(&image, "RICE_1", &CompressOptions::tiled([24, 16]))
        .unwrap();
    let mut r = FitsReader::open(Cursor::new(w.into_inner().into_inner())).unwrap();
    let back = match r.read_image(1).unwrap().decode() {
        ImageData::F32(v) => v,
        other => panic!("expected F32, got {other:?}"),
    };
    for (i, (&o, &b)) in orig.iter().zip(&back).enumerate() {
        if o.is_nan() {
            assert!(b.is_nan(), "null pixel {i} must round-trip to NaN");
        } else {
            assert!((o - b).abs() < 0.2, "pixel {i}: {o} vs {b}");
        }
    }
}

#[test]
fn dither2_quantize_round_trips() {
    use super::DitherMethod;
    use super::quantize::{dequantize_into, quantize_tile};

    // 8×8 field with genuine noise and a scattering of exact zeros.
    let mut data: Vec<f64> = (0..64)
        .map(|i| {
            let mut z = (i as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
            z ^= z >> 31;
            10.0 + (z % 1000) as f64 / 100.0
        })
        .collect();
    for &k in &[0usize, 13, 27, 40, 63] {
        data[k] = 0.0;
    }
    let irow = 7;
    let q = quantize_tile(&data, 8, 8, 0.0, DitherMethod::Subtractive2, irow).unwrap();
    // Exact zeros must encode to the reserved ZERO_VALUE.
    for &k in &[0usize, 13, 27, 40, 63] {
        assert_eq!(q.idata[k], super::quantize::ZERO_VALUE, "zero pixel {k}");
    }
    let ints: Vec<i64> = q.idata.iter().map(|&v| v as i64).collect();
    let mut back = Vec::new();
    dequantize_into(
        &ints,
        q.bscale,
        q.bzero,
        DitherMethod::Subtractive2,
        irow,
        None,
        &mut back,
    );
    for (i, (&o, &b)) in data.iter().zip(&back).enumerate() {
        if o == 0.0 {
            assert_eq!(b, 0.0, "zero pixel {i} must decode to exactly 0.0");
        } else {
            assert!(
                (o - b).abs() <= 0.5 * q.bscale + 1e-9,
                "pixel {i}: {o} vs {b}"
            );
        }
    }
}

#[test]
fn hcompress_lossy_write_round_trips_within_scale() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;

    // A smooth-ish 32×32 i32 image; lossy HCOMPRESS with SCALE=4.
    let samples: Vec<i32> = (0..32 * 32)
        .map(|i| 100 + 5 * (i % 32) + 3 * (i / 32))
        .collect();
    let image = Image {
        shape: vec![32, 32],
        samples: ImageData::I32(samples.clone()),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    let mut w = FitsWriter::new(Cursor::new(Vec::new()));
    w.write_compressed_image(
        &image,
        "HCOMPRESS_1",
        &CompressOptions {
            hcompress_scale: 4,
            ..CompressOptions::tiled([32, 32])
        },
    )
    .unwrap();
    let mut r = FitsReader::open(Cursor::new(w.into_inner().into_inner())).unwrap();
    let back = match r.read_image(1).unwrap().decode() {
        ImageData::I32(v) => v,
        other => panic!("expected I32, got {other:?}"),
    };
    // Lossy: each pixel within the quantization scale (±scale), and not identical.
    let max_err = samples
        .iter()
        .zip(&back)
        .map(|(a, b)| (a - b).abs())
        .max()
        .unwrap();
    assert!(
        max_err <= 4,
        "HCOMPRESS lossy error {max_err} exceeds scale"
    );
}

#[test]
fn plio_write_round_trips_through_decode() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;

    // PLIO is a mask codec: non-negative i32 values. value(x, y) = (x + y) % 7,
    // with a few longer runs to exercise multi-word counts.
    let samples: Vec<i32> = (0..24 * 16).map(|i| (i % 24 + i / 24) % 7).collect();
    let image = Image {
        shape: vec![24, 16],
        samples: ImageData::I32(samples.clone()),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    let mut w = FitsWriter::new(Cursor::new(Vec::new()));
    w.write_compressed_image(&image, "PLIO_1", &CompressOptions::default())
        .unwrap();
    let mut r = FitsReader::open(Cursor::new(w.into_inner().into_inner())).unwrap();
    match r.read_image(1).unwrap().decode() {
        ImageData::I32(v) => assert_eq!(v, samples, "PLIO_1 round-trip"),
        other => panic!("PLIO_1: expected I32, got {other:?}"),
    }
}

#[test]
fn decompresses_gzip_2_tiled_image() {
    check_decoded("comp_gzip2_i16.fits");
}

#[test]
fn decompresses_plio_1_mask() {
    // PLIO fixture encodes value(x, y) = (x + y) % 7 as an i32 mask.
    let mut f = open("comp_plio_i32.fits");
    let img = f.read_image(1).unwrap();
    assert_eq!(img.shape, vec![24, 16]);
    match img.decode() {
        ImageData::I32(v) => {
            assert_eq!(v.len(), 24 * 16);
            for (i, &got) in v.iter().enumerate() {
                let (x, y) = (i % 24, i / 24);
                assert_eq!(got, ((x + y) % 7) as i32, "pixel {i}");
            }
        }
        other => panic!("expected I32, got {other:?}"),
    }
}

/// Compare a compressed-float decode against astropy's reconstructed reference.
fn check_float(compressed: &str, reference: &str) {
    let got = match open(compressed).read_image(1).unwrap().decode() {
        ImageData::F32(v) => v,
        other => panic!("expected F32, got {other:?}"),
    };
    let want = match open(reference).read_image(0).unwrap().decode() {
        ImageData::F32(v) => v,
        other => panic!("expected F32 reference, got {other:?}"),
    };
    assert_eq!(got.len(), 24 * 16);
    assert_eq!(got, want, "{compressed} must match astropy");
}

#[test]
fn decompresses_unquantized_float_via_gzip_fallback() {
    // Smooth data stored losslessly: ZSCALE=0, raw floats gzip'd in
    // GZIP_COMPRESSED_DATA (COMPRESSED_DATA empty).
    check_float("comp_ricef_nodither.fits", "comp_ref_f32.fits");
}

#[test]
fn decompresses_quantized_float_no_dither() {
    // Noisy data genuinely quantized: per-tile ZSCALE≠0, integers RICE-packed in
    // COMPRESSED_DATA, dequantized as ZSCALE·int + ZZERO.
    check_float("comp_ricef_quant.fits", "comp_ref_quant_f32.fits");
}

/// Build a fixed-width BINTABLE, write it, then round-trip it through table
/// compression with `algo`/`rows_per_tile` and assert the data is byte-identical.
#[test]
fn decompresses_nocompress_tile_verbatim() {
    use crate::data::ImageData;
    use crate::header::Header;
    use crate::table::BinTable;
    // A 2×2 i16 image as a single NOCOMPRESS tile: the COMPRESSED_DATA cell holds
    // the four pixels verbatim as big-endian i16.
    let mut h = Header::new();
    h.set("XTENSION", "BINTABLE")
        .set("BITPIX", 8)
        .set("NAXIS", 2)
        .set("NAXIS1", 8) // one 1P descriptor
        .set("NAXIS2", 1) // one tile
        .set("PCOUNT", 8) // heap = 8 raw bytes
        .set("GCOUNT", 1)
        .set("TFIELDS", 1)
        .set("TFORM1", "1PB(8)")
        .set("TTYPE1", "COMPRESSED_DATA")
        .set("ZIMAGE", true)
        .set("ZCMPTYPE", "NOCOMPRESS")
        .set("ZBITPIX", 16)
        .set("ZNAXIS", 2)
        .set("ZNAXIS1", 2)
        .set("ZNAXIS2", 2)
        .set("ZTILE1", 2)
        .set("ZTILE2", 2);
    let mut data = Vec::new();
    data.extend_from_slice(&8i32.to_be_bytes()); // descriptor nelem = 8 bytes
    data.extend_from_slice(&0i32.to_be_bytes()); // descriptor offset = 0
    for x in [1i16, 2, 3, 4] {
        data.extend_from_slice(&x.to_be_bytes());
    }
    let table = BinTable::from_data(&h, data).unwrap();
    let img = decompress_image(&h, &table).unwrap();
    assert_eq!(img.shape, vec![2, 2]);
    assert_eq!(img.samples, ImageData::I16(vec![1, 2, 3, 4]));
}

#[test]
fn zblank_column_overrides_keyword_per_tile() {
    use crate::data::ImageData;
    use crate::header::Header;
    use crate::table::BinTable;
    // A 2×1 float image, one NOCOMPRESS tile of quantized i32 [10, 99]. ZSCALE=2,
    // ZZERO=5 ⇒ pixel 0 = 25.0; pixel 1's quantized int equals the per-tile ZBLANK
    // *column* value (99), so it decodes to NaN — proving the column drives nulls.
    let mut h = Header::new();
    h.set("XTENSION", "BINTABLE")
        .set("BITPIX", 8)
        .set("NAXIS", 2)
        .set("NAXIS1", 28) // 1P(8) + 1D + 1D + 1J
        .set("NAXIS2", 1)
        .set("PCOUNT", 8)
        .set("GCOUNT", 1)
        .set("TFIELDS", 4)
        .set("TFORM1", "1PB(8)")
        .set("TTYPE1", "COMPRESSED_DATA")
        .set("TFORM2", "1D")
        .set("TTYPE2", "ZSCALE")
        .set("TFORM3", "1D")
        .set("TTYPE3", "ZZERO")
        .set("TFORM4", "1J")
        .set("TTYPE4", "ZBLANK")
        .set("ZIMAGE", true)
        .set("ZCMPTYPE", "NOCOMPRESS")
        .set("ZBITPIX", -32)
        .set("ZNAXIS", 2)
        .set("ZNAXIS1", 2)
        .set("ZNAXIS2", 1)
        .set("ZTILE1", 2)
        .set("ZTILE2", 1);
    let mut data = Vec::new();
    data.extend_from_slice(&8i32.to_be_bytes()); // descriptor nelem
    data.extend_from_slice(&0i32.to_be_bytes()); // descriptor offset
    data.extend_from_slice(&2.0f64.to_be_bytes()); // ZSCALE
    data.extend_from_slice(&5.0f64.to_be_bytes()); // ZZERO
    data.extend_from_slice(&99i32.to_be_bytes()); // ZBLANK column value
    data.extend_from_slice(&10i32.to_be_bytes()); // heap: quantized int 0
    data.extend_from_slice(&99i32.to_be_bytes()); // heap: quantized int 1 (== ZBLANK)
    let table = BinTable::from_data(&h, data).unwrap();
    let img = decompress_image(&h, &table).unwrap();
    let ImageData::F32(px) = img.samples else {
        panic!("expected F32")
    };
    assert_eq!(px[0], 25.0);
    assert!(px[1].is_nan());
}

fn check_table_roundtrip(algo: &str, rows_per_tile: usize) {
    use crate::table::ColumnData;
    use crate::writer::{FitsWriter, WriteColumn};
    use std::io::Cursor;

    let nrows = 10;
    let col = |name: &str, data, repeat| WriteColumn::fixed(name, data, repeat);
    let columns = vec![
        col(
            "SHORT",
            ColumnData::I16((0..nrows).map(|i| i as i16 * 7 - 30).collect()),
            1,
        ),
        col(
            "INT",
            ColumnData::I32((0..nrows).map(|i| (i as i32) * 100_000 - 5).collect()),
            1,
        ),
        col(
            "FLT",
            ColumnData::F32((0..nrows).map(|i| i as f32 * 1.5 - 3.25).collect()),
            1,
        ),
        col(
            "DBL",
            ColumnData::F64((0..nrows).map(|i| i as f64 * 0.1).collect()),
            1,
        ),
        col(
            "BYTE",
            ColumnData::Bytes((0..nrows).map(|i| (i * 3) as u8).collect()),
            1,
        ),
        // A multi-element (repeat=3) short column.
        col(
            "VEC",
            ColumnData::I16((0..nrows * 3).map(|i| (i * 2) as i16).collect()),
            3,
        ),
    ];

    // 1. Write an uncompressed table and read it back.
    let mut w = FitsWriter::new(Cursor::new(Vec::new()));
    w.write_table(nrows, &columns).unwrap();
    let bytes = w.into_inner().into_inner();
    let mut r = FitsReader::open(Cursor::new(bytes)).unwrap();
    let orig = r.read_table(1).unwrap();
    let orig_header = r.hdus[1].header.clone();

    // 2. Compress it, then read + uncompress.
    let mut cw = FitsWriter::new(Cursor::new(Vec::new()));
    cw.write_compressed_table(&orig_header, &orig, rows_per_tile, algo)
        .unwrap();
    let cbytes = cw.into_inner().into_inner();
    let mut cr = FitsReader::open(Cursor::new(cbytes)).unwrap();
    let restored = cr.read_compressed_table(1).unwrap();

    // 3. The uncompressed table must be byte-identical to the original.
    assert_eq!(restored.nrows, orig.nrows, "{algo}/{rows_per_tile} nrows");
    assert_eq!(
        restored.row_len, orig.row_len,
        "{algo}/{rows_per_tile} row width"
    );
    assert_eq!(
        restored.raw_rows(),
        orig.raw_rows(),
        "{algo}/{rows_per_tile} data mismatch"
    );
}

#[test]
fn table_compression_round_trips() {
    // One tile, several tiles, and a tile smaller than the table — across codecs.
    for &rpt in &[10usize, 4, 1] {
        check_table_roundtrip("GZIP_1", rpt);
        check_table_roundtrip("GZIP_2", rpt);
        check_table_roundtrip("RICE_1", rpt);
    }
}

/// Emit our `write_compressed_table` output for external (cfitsio) validation of
/// the *encode* direction. After running, verify with cfitsio:
///   `funpack -O .tmp/my_unpk.fits .tmp/my_ctable.fits`
/// then compare `.tmp/my_unpk.fits` against `comp_table_ref.fits` — they match.
/// Run with `cargo test --features compression -- --ignored emit_compressed_table`.
#[test]
#[ignore]
fn emit_compressed_table_for_funpack() {
    use crate::writer::FitsWriter;
    use std::fs::File;

    let src = std::fs::read("tests/data/fits/comp_table_ref.fits").unwrap();
    let mut r = FitsReader::open(std::io::Cursor::new(src)).unwrap();
    let table = r.read_table(1).unwrap();
    let header = r.hdus[1].header.clone();
    let mut w = FitsWriter::new(File::create(".tmp/my_ctable.fits").unwrap());
    w.write_compressed_table(&header, &table, 100, "RICE_1")
        .unwrap();
}

#[test]
fn decodes_a_cfitsio_compressed_table() {
    // Ground truth: `comp_table_cfitsio.fits` was produced by cfitsio's `fpack
    // -tableonly` from `comp_table_ref.fits` (500 rows, 6 fixed-width columns).
    // fpack chose a real mix of per-column codecs — GZIP_2 (short/float/double),
    // RICE_1 (the int32 columns), GZIP_1 (byte) — so this exercises every decode
    // path against an independent implementation. Our uncompressed output must be
    // byte-identical to the original table.
    let restored = open("comp_table_cfitsio.fits")
        .read_compressed_table(1)
        .unwrap();
    let original = open("comp_table_ref.fits").read_table(1).unwrap();

    assert_eq!(restored.nrows, 500);
    assert_eq!(restored.nrows, original.nrows);
    assert_eq!(restored.row_len, original.row_len);
    assert_eq!(restored.columns.len(), 6);
    assert_eq!(
        restored.raw_rows(),
        original.raw_rows(),
        "decoded cfitsio-compressed table must match the original bytes"
    );
    // Spot-check a decoded value against the known formula (INT = i·100000 − 5).
    match original.column_by_idx(1).unwrap().raw().unwrap() {
        ColumnData::I32(v) => assert_eq!(v[3], 3 * 100_000 - 5),
        other => panic!("expected I32, got {other:?}"),
    }
}

#[test]
fn compressed_table_with_vla_column_is_rejected_cleanly() {
    // `comp_table_vla.fits` (from `fpack -tableonly` of a table with a `PJ` VLA
    // column) — decoding VLA columns inside a compressed table is not yet
    // implemented, so it must error rather than misread or panic.
    let mut f = open("comp_table_vla.fits");
    assert!(matches!(
        f.read_compressed_table(1),
        Err(FitsError::UnsupportedCompression { .. })
    ));
}

#[test]
fn read_compressed_table_rejects_a_plain_bintable() {
    let mut f = open("DDTSUVDATA.fits");
    assert!(matches!(
        f.read_compressed_table(1),
        Err(FitsError::NotCompressedTable)
    ));
}

#[test]
fn reading_a_plain_bintable_as_an_image_is_rejected() {
    // DDTSUVDATA hdu 1 is an ordinary BINTABLE (no ZIMAGE).
    let mut f = open("DDTSUVDATA.fits");
    // Public path: `read_image` sees a non-ZIMAGE bintable and rejects it as a
    // non-image (it never reaches the decompressor).
    assert!(matches!(f.read_image(1), Err(FitsError::NotAnImage)));
    // The decompressor itself still guards its `ZIMAGE` precondition.
    let table = f.read_table(1).unwrap();
    assert!(matches!(
        decompress_image(&f.hdus[1].header, &table),
        Err(FitsError::NotCompressedImage)
    ));
}

#[test]
fn integer_image_compression_preserves_bscale_bzero_and_blank() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;
    // §10.2: the compressed tiles store *raw* stored integers, so BSCALE/BZERO and
    // the BLANK sentinel must survive in the rebuilt header (was dropped before).
    let samples: Vec<i16> = (0..24 * 16).map(|i| (i % 50) as i16 - 5).collect();
    let image = Image {
        shape: vec![24, 16],
        samples: ImageData::I16(samples.clone()),
        scaling: Scaling {
            bscale: 2.5,
            bzero: 100.0,
            blank: Some(-5),
        },
    };
    let mut w = FitsWriter::new(Cursor::new(Vec::new()));
    w.write_compressed_image(&image, "GZIP_1", &CompressOptions::default())
        .unwrap();
    let mut r = FitsReader::open(Cursor::new(w.into_inner().into_inner())).unwrap();
    let back = r.read_image(1).unwrap();

    match back.decode() {
        ImageData::I16(v) => assert_eq!(v, samples, "raw samples"),
        other => panic!("expected I16, got {other:?}"),
    }
    assert_eq!(back.scaling.bscale, 2.5);
    assert_eq!(back.scaling.bzero, 100.0);
    assert_eq!(back.scaling.blank, Some(-5));
}

#[test]
fn rice_rejects_64_bit_pixels() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;
    // Table 37 permits BYTEPIX 8, but the 64-bit RICE bitstream is unsupported; it
    // must error cleanly rather than panic / silently corrupt.
    let image = Image {
        shape: vec![4],
        samples: ImageData::I64(vec![1, 2, 3, 4]),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    let mut w = FitsWriter::new(Cursor::new(Vec::new()));
    assert!(matches!(
        w.write_compressed_image(&image, "RICE_1", &CompressOptions::default()),
        Err(FitsError::UnsupportedCompression { .. })
    ));
    // GZIP handles 64-bit fine — the rejection is RICE-specific.
    let mut w2 = FitsWriter::new(Cursor::new(Vec::new()));
    assert!(
        w2.write_compressed_image(&image, "GZIP_1", &CompressOptions::default())
            .is_ok()
    );
}

#[test]
fn nocompress_image_round_trips() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;
    // §10.4: tiles stored verbatim (uncompressed big-endian pixels) round-trip.
    let samples: Vec<i16> = (0..24 * 16)
        .map(|i| (i % 24) as i16 * 7 - (i / 24) as i16 * 5)
        .collect();
    let image = Image {
        shape: vec![24, 16],
        samples: ImageData::I16(samples.clone()),
        scaling: Scaling {
            bscale: 1.0,
            bzero: 0.0,
            blank: None,
        },
    };
    let mut w = FitsWriter::new(Cursor::new(Vec::new()));
    w.write_compressed_image(&image, "NOCOMPRESS", &CompressOptions::default())
        .unwrap();
    let mut r = FitsReader::open(Cursor::new(w.into_inner().into_inner())).unwrap();
    match r.read_image(1).unwrap().decode() {
        ImageData::I16(v) => assert_eq!(v, samples),
        other => panic!("expected I16, got {other:?}"),
    }
}

#[test]
fn empty_naxis0_image_round_trips() {
    use crate::data::{Image, ImageData, Scaling};
    use crate::writer::FitsWriter;
    use std::io::Cursor;
    // A `NAXIS = 0` image has no data array. The encoder must emit an empty
    // `NAXIS2 = 0` ZIMAGE table (not panic on a fabricated phantom tile), and the
    // decoder must restore the empty image. Exercise both the integer and float
    // encoder paths.
    let cases = [
        ImageData::I16(Vec::new()),
        ImageData::I32(Vec::new()),
        ImageData::F32(Vec::new()),
    ];
    for samples in cases {
        let image = Image {
            shape: Vec::new(),
            samples: samples.clone(),
            scaling: Scaling {
                bscale: 1.0,
                bzero: 0.0,
                blank: None,
            },
        };
        let mut w = FitsWriter::new(Cursor::new(Vec::new()));
        w.write_compressed_image(&image, "GZIP_1", &CompressOptions::default())
            .unwrap();
        let mut r = FitsReader::open(Cursor::new(w.into_inner().into_inner())).unwrap();
        let back = r.read_image(1).unwrap();
        assert!(back.shape.is_empty(), "shape for {samples:?}");
        // Same empty variant back out, no phantom pixel.
        match (back.decode(), &samples) {
            (ImageData::I16(v), ImageData::I16(_)) => assert!(v.is_empty()),
            (ImageData::I32(v), ImageData::I32(_)) => assert!(v.is_empty()),
            (ImageData::F32(v), ImageData::F32(_)) => assert!(v.is_empty()),
            (other, _) => panic!("variant mismatch: {other:?}"),
        }
    }
}

#[test]
fn compressed_image_descriptor_switches_to_q_for_large_offsets() {
    // §10.1.3: a heap offset beyond the 32-bit P range needs a 64-bit Q descriptor.
    let mut q = Vec::new();
    super::push_pq_descriptor(&mut q, true, 3, u32::MAX as u64 + 8);
    assert_eq!(q.len(), 16);
    assert_eq!(i64::from_be_bytes(q[0..8].try_into().unwrap()), 3);
    assert_eq!(
        i64::from_be_bytes(q[8..16].try_into().unwrap()),
        u32::MAX as i64 + 8
    );
    let mut p = Vec::new();
    super::push_pq_descriptor(&mut p, false, 3, 40);
    assert_eq!(p.len(), 8);
    assert_eq!(i32::from_be_bytes(p[4..8].try_into().unwrap()), 40);
}

#[test]
fn hcompress_tile_rejects_dimension_mismatch() {
    // Encode a valid 2×3 tile, then decode it claiming a different element count.
    // The decoder reads nx/ny from the stream and must cross-check them against the
    // tile size it was handed — rather than allocate/transform `nx*ny` blindly
    // (a wild-allocation / overflow / empty-buffer-panic guard, R2-4).
    let vals: Vec<i64> = vec![10, 20, 30, 40, 50, 60];
    let bytes = hcompress::hcompress_tile_encode(&vals, &[2, 3], 0).unwrap();
    // The correct element count round-trips losslessly (scale 0).
    let mut out = Vec::new();
    hcompress::hcompress_tile_into(&bytes, false, 6, &mut out).unwrap();
    assert_eq!(out, vals);
    // A mismatched element count is rejected, not decoded.
    assert!(hcompress::hcompress_tile_into(&bytes, false, 7, &mut out).is_err());
    assert!(hcompress::hcompress_tile_into(&bytes, false, 5, &mut out).is_err());
}

#[test]
fn decompress_image_rejects_overflowing_znaxis_product() {
    use crate::error::FitsError;
    use crate::header::Header;
    use crate::table::BinTable;
    // ZNAXIS1·ZNAXIS2 = 5e9·5e9 = 2.5e19 overflows usize; decode must reject the
    // header up front (before allocating the output plane), not wrap to a small
    // buffer and then scatter out of bounds (R2-2).
    let mut h = Header::new();
    h.set("XTENSION", "BINTABLE")
        .set("BITPIX", 8)
        .set("NAXIS", 2)
        .set("NAXIS1", 8)
        .set("NAXIS2", 1)
        .set("PCOUNT", 0)
        .set("GCOUNT", 1)
        .set("TFIELDS", 1)
        .set("TFORM1", "1PB(0)")
        .set("TTYPE1", "COMPRESSED_DATA")
        .set("ZIMAGE", true)
        .set("ZCMPTYPE", "GZIP_1")
        .set("ZBITPIX", 16)
        .set("ZNAXIS", 2)
        .set("ZNAXIS1", 5_000_000_000i64)
        .set("ZNAXIS2", 5_000_000_000i64);
    let mut data = Vec::new();
    data.extend_from_slice(&0i32.to_be_bytes()); // empty P descriptor: nelem
    data.extend_from_slice(&0i32.to_be_bytes()); // offset
    let table = BinTable::from_data(&h, data).unwrap();
    assert!(matches!(
        decompress_image(&h, &table),
        Err(FitsError::DataUnitOverflow)
    ));
}

#[test]
fn uncompress_table_rejects_overflowing_row_product() {
    use crate::error::FitsError;
    use crate::header::Header;
    use crate::table::BinTable;
    // ZNAXIS2·ZNAXIS1 = 3e18·8 = 2.4e19 overflows usize; uncompress must reject the
    // header before allocating the row buffer (R2-3).
    let mut h = Header::new();
    h.set("XTENSION", "BINTABLE")
        .set("BITPIX", 8)
        .set("NAXIS", 2)
        .set("NAXIS1", 16) // one 1QB descriptor row
        .set("NAXIS2", 1)
        .set("PCOUNT", 0)
        .set("GCOUNT", 1)
        .set("TFIELDS", 1)
        .set("TFORM1", "1QB")
        .set("TTYPE1", "C1")
        .set("ZTABLE", true)
        .set("ZTILELEN", 1)
        .set("ZNAXIS1", 8)
        .set("ZNAXIS2", 3_000_000_000_000_000_000i64)
        .set("ZFORM1", "1K");
    let mut data = Vec::new();
    data.extend_from_slice(&0i64.to_be_bytes()); // Q descriptor: nelem
    data.extend_from_slice(&0i64.to_be_bytes()); // offset
    let table = BinTable::from_data(&h, data).unwrap();
    assert!(matches!(
        uncompress_table(&h, &table),
        Err(FitsError::DataUnitOverflow)
    ));
}

#[test]
fn decompress_image_rejects_oversized_znaxis_product() {
    use crate::error::FitsError;
    use crate::header::Header;
    use crate::table::BinTable;
    // ZNAXIS1 = 2^60 does NOT overflow usize (so the overflow guard passes), but
    // allocating that many bytes would abort the process. The output plane is
    // allocated fallibly (`try_reserve`), so decode must return a recoverable error.
    let mut h = Header::new();
    h.set("XTENSION", "BINTABLE")
        .set("BITPIX", 8)
        .set("NAXIS", 2)
        .set("NAXIS1", 8)
        .set("NAXIS2", 1)
        .set("PCOUNT", 0)
        .set("GCOUNT", 1)
        .set("TFIELDS", 1)
        .set("TFORM1", "1PB(0)")
        .set("TTYPE1", "COMPRESSED_DATA")
        .set("ZIMAGE", true)
        .set("ZCMPTYPE", "GZIP_1")
        .set("ZBITPIX", 8)
        .set("ZNAXIS", 1)
        .set("ZNAXIS1", 1i64 << 60);
    let mut data = Vec::new();
    data.extend_from_slice(&0i32.to_be_bytes()); // empty P descriptor: nelem
    data.extend_from_slice(&0i32.to_be_bytes()); // offset
    let table = BinTable::from_data(&h, data).unwrap();
    assert!(matches!(
        decompress_image(&h, &table),
        Err(FitsError::DataUnitTooLarge { .. })
    ));
}

#[test]
fn gunzip_rejects_a_stream_larger_than_its_tile() {
    use crate::error::FitsError;
    // A decompression bomb: a tiny gzip stream inflating far past the tile's known
    // byte size must be rejected, not allocated unbounded.
    let big = vec![0u8; 100_000];
    let bomb = super::gzip::gzip_encode(&big, 9);
    assert!(bomb.len() < 1000, "an all-zero buffer compresses tiny");
    // Bounded at 1 KiB (a small tile): inflating to 100 KB overruns → error.
    assert!(matches!(
        super::gzip::gunzip(&bomb, 1024),
        Err(FitsError::UnsupportedCompression { .. })
    ));
    // One byte short of the true size still overruns (the +1 detection boundary).
    assert!(super::gzip::gunzip(&bomb, big.len() - 1).is_err());
    // Bounded at the true size: decodes to exactly the original bytes.
    assert_eq!(super::gzip::gunzip(&bomb, big.len()).unwrap(), big);
}

#[test]
fn i64_be_round_trip_and_buffer_reuse() {
    // I16: narrowing (`as i16`) then big-endian packing, hand-computed. -1 → 0xFFFF,
    // 258 → 0x0102, -2 → 0xFFFE.
    let vals = [0i64, 1, -1, 258, -2];
    let want = [0, 0, 0, 1, 0xFF, 0xFF, 0x01, 0x02, 0xFF, 0xFE];
    assert_eq!(i64_to_be(&vals, Bitpix::I16), want);
    // Decode is the exact inverse (sign-extending back to i64), into a reused buffer.
    let mut d = Vec::new();
    be_to_i64_into(&want, Bitpix::I16, &mut d);
    assert_eq!(d, vals);

    // I32 packs four bytes/elem; 0x00010203 = 66051, -1 → all-ones.
    let i32_vals = [66051i64, -1];
    assert_eq!(
        i64_to_be(&i32_vals, Bitpix::I32),
        [0x00, 0x01, 0x02, 0x03, 0xFF, 0xFF, 0xFF, 0xFF]
    );
    be_to_i64_into(&[0, 1, 2, 3, 0xFF, 0xFF, 0xFF, 0xFF], Bitpix::I32, &mut d);
    assert_eq!(d, i32_vals);

    // U8 and I64 ends of the range.
    assert_eq!(i64_to_be(&[255, 0], Bitpix::U8), [0xFF, 0x00]);
    be_to_i64_into(&[0xFF, 0x00], Bitpix::U8, &mut d);
    assert_eq!(d, [255, 0]);
    assert_eq!(i64_to_be(&[-1], Bitpix::I64), [0xFF; 8]);

    // `i64_to_be_into` clears + resizes its scratch: a long fill followed by a short
    // one must leave exactly the short result (no stale tail), matching the owning form.
    let mut buf = Vec::new();
    i64_to_be_into(&[7, 8, 9, 10], Bitpix::I16, &mut buf);
    i64_to_be_into(&vals, Bitpix::I16, &mut buf);
    assert_eq!(buf, want);
}

#[test]
fn tile_into_rows_match_layout() {
    // Expand the row representation (row_bases + row_len) back to the full flat-index
    // list it stands for, so the hand-computed expectations still describe pixels.
    fn flat(s: &TileScratch) -> Vec<usize> {
        let mut v = Vec::new();
        for &b in &s.row_bases {
            v.extend(b..b + s.row_len);
        }
        v
    }

    // 4×3 image (x fastest, strides [1, 4]) tiled 2×2 → 4 tiles; the top row of tiles
    // is edge-clipped to height 1. Same scratch reused across tiles (no stale rows).
    let geom = TileGeometry::new(&[4, 3], &[2, 2]);
    assert_eq!(geom.ntiles(), 4);
    let mut s = TileScratch::default();
    let expect = [
        (vec![0, 4], 2, vec![0, 1, 4, 5]), // origin (0,0), full 2×2: two rows of 2
        (vec![2, 6], 2, vec![2, 3, 6, 7]), // origin (2,0)
        (vec![8], 2, vec![8, 9]),          // origin (0,2), clipped to height 1
        (vec![10], 2, vec![10, 11]),       // origin (2,2), clipped to height 1
    ];
    for (t, (bases, row_len, pixels)) in expect.iter().enumerate() {
        geom.tile_into(t, &mut s);
        assert_eq!(&s.row_bases, bases, "tile {t} row bases");
        assert_eq!(s.row_len, *row_len, "tile {t} row len");
        assert_eq!(s.nelem(), pixels.len(), "tile {t} nelem");
        assert_eq!(&flat(&s), pixels, "tile {t} pixels");
    }

    // 4×2×2 image (strides [1, 4, 8]) tiled 2×2×2 → 2 tiles: exercises the odometer
    // carrying across the two higher axes; each tile is 4 rows of 2.
    let geom3 = TileGeometry::new(&[4, 2, 2], &[2, 2, 2]);
    assert_eq!(geom3.ntiles(), 2);
    geom3.tile_into(0, &mut s);
    assert_eq!(s.row_bases, vec![0, 4, 8, 12]);
    assert_eq!(flat(&s), vec![0, 1, 4, 5, 8, 9, 12, 13]);
    geom3.tile_into(1, &mut s);
    assert_eq!(flat(&s), vec![2, 3, 6, 7, 10, 11, 14, 15]);
}