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<Vec<u8>, 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);
+    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<u8> {
+        (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] = &[
+            // <i4, 10 px/scanline over 4 px/block: padded scanlines + byte planes.
+            ("noencoder /noencoder_szip_dset.h5", noencoder, 6040, 16, [168, 4, 32, 10],
+             "00000000010000000200000003000000040000000500000006000000070000000800000009000000"),
+            // <f4, LSB + NN.
+            ("le_data /Szip_float_data_le", le_data, 55224, 48, [169, 4, 32, 12],
+             "abaaaa3eabaa2a3f0000803fabaa2a3f0000803fabaaaa3f0000803fabaaaa3f5555d53fabaaaa3f5555d53f00000040"),
+            // >f4, MSB + NN.
+            ("le_data /Szip_float_data_be", le_data, 55396, 48, [177, 4, 32, 12],
+             "3eaaaaab3f2aaaab3f8000003f2aaaab3f8000003faaaaab3f8000003faaaaab3fd555553faaaaab3fd5555540000000"),
+            // <f8 (64-bit), NN.
+            ("szip_h5py /f8", h5py, 4016, 100, [169, 8, 64, 10],
+             "00000000000008c000000000000008c000000000000008c000000000000008c000000000000004c000000000000004c000000000000004c000000000000004c000000000000000c000000000000000c000000000000000c000000000000000c0000000000000f8bf000000000000f8bf000000000000f8bf000000000000f8bf000000000000f0bf000000000000f0bf000000000000f0bf000000000000f0bf000000000000e0bf000000000000e0bf000000000000e0bf000000000000e0bf0000000000000000000000000000000000000000000000000000000000000000000000000000e03f000000000000e03f000000000000e03f000000000000e03f000000000000f03f000000000000f03f000000000000f03f000000000000f03f000000000000f83f000000000000f83f000000000000f83f000000000000f83f"),
+            // <i8 (64-bit), entropy coding without NN.
+            ("szip_h5py /i8", h5py, 4188, 53, [141, 4, 64, 10],
+             "000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000300000000000000030000000000000003000000000000000300000000000000030000000000000003000000000000000300000000000000030000000000000006000000000000000600000000000000060000000000000006000000000000000600000000000000060000000000000006000000000000000600000000000000090000000000000009000000000000000900000000000000090000000000000009000000000000000900000000000000090000000000000009000000000000000c000000000000000c000000000000000c000000000000000c000000000000000c000000000000000c000000000000000c000000000000000c00000000000000"),
+            // <u2, 35 px/scanline over 8 px/block: padded scanlines, 16-bit samples.
+            ("szip_h5py /u2", h5py, 4308, 43, [169, 8, 16, 35],
+             "00000000000000006100610061006100c200c200c200c20023012301230123018401840184018401e501e501e501e5014602460246024602a702a702a702a702080308030803"),
+        ];
+        for (name, file, off, len, cd, want) in cases {
+            let want = unhex(want);
+            let got = szip_decompress(&file[*off..off + len], cd, want.len())
+                .unwrap_or_else(|e| panic!("{name}: {e:?}"));
+            assert_eq!(got, want, "{name}");
+        }
+    }
+
     #[test]
     fn szip_disabled_returns_unsupported() {
         #[cfg(not(feature = "szip"))]
@@ -132,6 +264,8 @@ mod tests {
         assert_eq!(rc, 0, "aec_buffer_encode failed: {rc}");
         let enc_len = encoded.len() - enc.avail_out;
         encoded.truncate(enc_len);
+        // H5Zszip.c prefixes the stream with the uncompressed size.
+        encoded.splice(0..0, (original.len() as u32).to_le_bytes());
 
         // Decode through our public interface.
         // cd[0]=0 (no NN bit 0x20), cd[1]=8 (ppb), cd[2]=8 (bpp), cd[3]=1024 (pps).
@@ -163,6 +297,8 @@ mod tests {
         assert_eq!(rc, 0, "aec_buffer_encode with NN failed: {rc}");
         let enc_len = encoded.len() - enc.avail_out;
         encoded.truncate(enc_len);
+        // H5Zszip.c prefixes the stream with the uncompressed size.
+        encoded.splice(0..0, (original.len() as u32).to_le_bytes());
 
         // cd[0] = 0x20 (H5_SZIP_NN_OPTION_MASK) → decoder must set AEC_DATA_PREPROCESS.
         let cd = [0x20u32, 8, 8, 1024];

crates/clawhdf5-format/tests/fixtures/filters/README.md

@@ -8,3 +8,6 @@ the tests come from h5py's `DatasetID.get_chunk_info`.
 | File | Origin | Licence |
 |------|--------|---------|
 | `h5ex_d_lz4.h5` | HDF Group `HDF5Examples/C/H5FLT/tfiles/h5ex_d_lz4.h5` (hdf5 repository) | HDF5 licence (BSD-3-Clause style) |
+| `noencoder.h5` | HDF Group `test/testfiles/noencoder.h5` (hdf5 repository) | HDF5 licence (BSD-3-Clause style) |
+| `le_data.h5` | HDF Group `test/testfiles/le_data.h5` (hdf5 repository) | HDF5 licence (BSD-3-Clause style) |
+| `szip_h5py.h5` | Written for these tests with h5py 3 / libhdf5 2.0.0 (libaec szip): `f8` (8x10, chunks 4x10, `('nn', 8)`), `i8` (8x10, chunks 4x10, `('ec', 4)`), `u2` (70, chunks 35, `('nn', 8)`) | Same as this repository |