diff --git a/crates/clawhdf5-format/src/filters_szip.rs b/crates/clawhdf5-format/src/filters_szip.rs index b05a5e4..4e7fd71 100644 --- a/crates/clawhdf5-format/src/filters_szip.rs +++ b/crates/clawhdf5-format/src/filters_szip.rs @@ -1,19 +1,39 @@ //! SZIP (libaec Adaptive Entropy Coding) decompression. //! //! Gated by the `szip` feature which links against the system libaec library. +//! +//! libhdf5's SZIP filter (`H5Zszip.c`) prefixes each chunk with its +//! uncompressed size and hands the rest to szlib's `SZ_BufftoBuffDecompress`. +//! libaec implements that call (`sz_compat.c`) on top of `aec_buffer_decode` +//! with some reshaping — 32/64-bit samples are coded as byte planes of 8-bit +//! samples, and scanlines that are not a whole number of blocks are padded — +//! which [`szip_decompress`] reproduces so its output matches libhdf5's. #[cfg(not(feature = "std"))] use alloc::vec::Vec; use crate::error::FormatError; -/// Decompress SZIP-compressed data using libaec. +/// `SZ_MSB_OPTION_MASK`: samples are big-endian. +#[cfg(feature = "szip")] +const SZ_MSB_OPTION_MASK: u32 = 16; +/// `SZ_NN_OPTION_MASK`: nearest-neighbour preprocessing. +#[cfg(feature = "szip")] +const SZ_NN_OPTION_MASK: u32 = 32; + +/// Decompress one SZIP-filtered chunk. /// -/// `cd` is the HDF5 SZIP filter client data (matches `H5Z_SZIP_PARM_*` indices): -/// cd[0] = options mask (`H5_SZIP_NN_OPTION_MASK = 0x20` enables NN preprocessing) -/// cd[1] = pixels per block (H5Z_SZIP_PARM_PPB; 8, 10, 16, or 32) -/// cd[2] = bits per sample (H5Z_SZIP_PARM_BPP; element bit width) -/// cd[3] = pixels per scan line (H5Z_SZIP_PARM_PPS; informational only) +/// `cd` is the HDF5 SZIP filter client data (`H5Z_SZIP_PARM_*` indices): +/// cd[0] = options mask (`SZ_*_OPTION_MASK`: 16 = MSB byte order, +/// 32 = nearest-neighbour preprocessing; K13/EC/LSB/RAW bits carry +/// no decoding information for libaec) +/// cd[1] = pixels per block +/// cd[2] = bits per pixel (sample precision, rounded up to 32 or 64 above +/// 24 by libhdf5) +/// cd[3] = pixels per scanline +/// +/// The chunk is a 4-byte little-endian uncompressed size followed by the +/// szlib stream. pub(crate) fn szip_decompress( _data: &[u8], _cd: &[u32], @@ -33,62 +53,174 @@ pub(crate) fn szip_decompress( #[cfg(feature = "szip")] fn szip_decode_impl(data: &[u8], cd: &[u32], chunk_size: usize) -> Result, FormatError> { - if cd.len() < 3 { - return Err(FormatError::ChunkedReadError( - "szip: missing client data".into(), - )); + let err = |m: &str| FormatError::ChunkedReadError(format!("szip: {m}")); + if cd.len() < 4 { + return Err(err("missing client data")); } let options = cd[0]; - let pixels_per_block = cd[1]; - let bits_per_sample = cd[2]; // H5Z_SZIP_PARM_BPP - if bits_per_sample == 0 || bits_per_sample > 32 { - return Err(FormatError::ChunkedReadError( - "szip: invalid bits per sample".into(), - )); + let pixels_per_block = cd[1] as usize; + let bits_per_pixel = cd[2]; + let pixels_per_scanline = cd[3] as usize; + if !(1..=32).contains(&bits_per_pixel) && bits_per_pixel != 64 { + return Err(err("invalid bits per sample")); } - if chunk_size == 0 { - return Err(FormatError::ChunkedReadError( - "szip: unknown output size".into(), - )); + if pixels_per_block == 0 || pixels_per_scanline == 0 { + return Err(err("invalid block or scanline size")); } - if data.is_empty() { - return Err(FormatError::ChunkedReadError("szip: empty input".into())); + if data.len() < 4 { + return Err(err("chunk too short")); } + // H5Zszip.c: UINT32DECODE of the uncompressed size, then the stream. + let dest_len = u32::from_le_bytes([data[0], data[1], data[2], data[3]]) as usize; + let limit = if chunk_size != 0 { + chunk_size + } else { + crate::filters::MAX_DECOMPRESS_SIZE + }; + if dest_len > limit { + return Err(err("declared size exceeds chunk size")); + } + let stream = &data[4..]; - // Map HDF5 option mask to libaec flags. - // HDF5 always stores SZIP data in MSB order, so AEC_DATA_MSB is unconditional. - // H5_SZIP_NN_OPTION_MASK (0x20): NN differential preprocessing. - let mut flags: u32 = libaec_sys::AEC_DATA_MSB; - if options & 0x20 != 0 { + // --- libaec sz_compat.c: SZ_BufftoBuffDecompress --- + let rsi = pixels_per_scanline.div_ceil(pixels_per_block); + let mut flags = 0; + if options & SZ_MSB_OPTION_MASK != 0 { + flags |= libaec_sys::AEC_DATA_MSB; + } + if options & SZ_NN_OPTION_MASK != 0 { flags |= libaec_sys::AEC_DATA_PREPROCESS; } + let pad_scanline = !pixels_per_scanline.is_multiple_of(pixels_per_block); + let deinterleave = bits_per_pixel == 32 || bits_per_pixel == 64; + let bits_per_sample = if deinterleave { 8 } else { bits_per_pixel }; + let pixel_size = match bits_per_sample { + 17.. => 4, + 9.. => 2, + _ => 1, + }; + let scanlines = (dest_len / pixel_size).div_ceil(pixels_per_scanline); + let buf_size = if pad_scanline { + rsi.checked_mul(pixels_per_block) + .and_then(|n| n.checked_mul(pixel_size)) + .and_then(|n| n.checked_mul(scanlines)) + .filter(|&n| n <= crate::filters::MAX_DECOMPRESS_SIZE.max(limit)) + .ok_or_else(|| err("scanline padding too large"))? + } else { + dest_len + }; - let mut out = vec![0u8; chunk_size]; + let mut buf = vec![0u8; buf_size]; let mut strm = libaec_sys::AecStream::zeroed(); - strm.next_in = data.as_ptr(); - strm.avail_in = data.len(); - strm.next_out = out.as_mut_ptr(); - strm.avail_out = chunk_size; + strm.next_in = stream.as_ptr(); + strm.avail_in = stream.len(); + strm.next_out = buf.as_mut_ptr(); + strm.avail_out = buf_size; strm.bits_per_sample = bits_per_sample; - strm.block_size = pixels_per_block; - strm.rsi = 128; // HDF5 default: 128 blocks per reference sample interval + strm.block_size = pixels_per_block as u32; + strm.rsi = rsi as u32; strm.flags = flags; - + // SAFETY: next_in/avail_in and next_out/avail_out describe live buffers + // (`stream` and `buf`) that outlive the call. let result = unsafe { libaec_sys::aec_buffer_decode(&mut strm) }; if result != 0 { return Err(FormatError::DecompressionError(format!( "szip: libaec error {result}" ))); } - let decoded_len = chunk_size - strm.avail_out; - out.truncate(decoded_len); - Ok(out) + let mut total_out = strm.total_out; + if pad_scanline { + let line = pixels_per_scanline * pixel_size; + let padded_line = rsi * pixels_per_block * pixel_size; + // remove_padding: compact each padded line down to `line` bytes. + let mut i = line; + let mut j = padded_line; + while j < total_out { + let end = (j + line).min(buf.len()); + buf.copy_within(j..end, i); + i += line; + j += padded_line; + } + total_out = scanlines * line; + } + if total_out < dest_len { + return Err(err("stream decoded to fewer bytes than declared")); + } + buf.truncate(dest_len); + if deinterleave { + // deinterleave_buffer: byte planes back into words. + let w = (bits_per_pixel / 8) as usize; + let n = dest_len / w; + let mut out = vec![0u8; dest_len]; + for i in 0..n { + for j in 0..w { + out[i * w + j] = buf[j * n + i]; + } + } + Ok(out) + } else { + Ok(buf) + } } #[cfg(test)] mod tests { use super::*; + #[cfg(feature = "szip")] + fn unhex(s: &str) -> Vec { + (0..s.len()) + .step_by(2) + .map(|i| u8::from_str_radix(&s[i..i + 2], 16).unwrap()) + .collect() + } + + /// SZIP chunks written by libhdf5, decoded exactly as libhdf5 decodes + /// them. Each case: fixture, chunk byte offset and size (from h5py's + /// `get_chunk_info`), the filter's cd_values, and the chunk's values as + /// h5py reads them (file byte order, hex). Before the fix every one of + /// these came back as garbage or zeros (or "invalid bits per sample" for + /// 64-bit): the 4-byte size prefix was fed to libaec, 32/64-bit samples + /// were not de-interleaved from byte planes, the reference sample + /// interval was fixed at 128 instead of derived from the scanline, padded + /// scanlines were not unpadded, and LE data was decoded as MSB. + #[cfg(feature = "szip")] + #[test] + fn szip_decodes_libhdf5_chunks_exactly() { + /// (name, file, chunk offset, chunk size, cd_values, decoded hex) + type Case<'a> = (&'a str, &'a [u8], usize, usize, [u32; 4], &'a str); + let noencoder: &[u8] = include_bytes!("../tests/fixtures/filters/noencoder.h5"); + let le_data: &[u8] = include_bytes!("../tests/fixtures/filters/le_data.h5"); + let h5py: &[u8] = include_bytes!("../tests/fixtures/filters/szip_h5py.h5"); + #[rustfmt::skip] + let cases: &[Case] = &[ + // f4, MSB + NN. + ("le_data /Szip_float_data_be", le_data, 55396, 48, [177, 4, 32, 12], + "3eaaaaab3f2aaaab3f8000003f2aaaab3f8000003faaaaab3f8000003faaaaab3fd555553faaaaab3fd5555540000000"), + //