From 5935e13866e30b9a8c362a7733f0fc90f916711a Mon Sep 17 00:00:00 2001 From: osobh Date: Fri, 25 Sep 2026 21:11:24 -0500 Subject: [PATCH] fix(format): decode SZIP chunks the way libhdf5 does SZIP-filtered datasets from libhdf5 came back as garbage or zeros with no error (ref_szip.h5, h5repack_szip.h5, noencoder.h5, le_data/be_data Szip_float_data_*), and 64-bit ones failed with "invalid bits per sample" (h5wasm compressed.h5). The decoder called aec_buffer_decode directly, but libhdf5 goes through szlib's SZ_BufftoBuffDecompress (H5Zszip.c), which libaec implements with reshaping (sz_compat.c). Differences, all fixed: - H5Zszip.c prefixes the stream with the 4-byte LE uncompressed size; it was fed to libaec as data. - 32- and 64-bit samples are coded as byte planes of 8-bit samples and must be de-interleaved. - The reference sample interval is ceil(pixels_per_scanline / pixels_per_block), not a fixed 128. - Scanlines that are not a whole number of blocks are padded and must be unpadded. - Byte order comes from the MSB option bit; LE data was decoded as MSB. Test: szip_decodes_libhdf5_chunks_exactly compares chunks from HDF Group test files (noencoder.h5, le_data.h5) and an h5py-written file (64-bit, 16-bit, padded scanlines, NN and EC) byte for byte with h5py's values; it failed before on the first case. Co-Authored-By: Claude Opus 5.5 (1M context) --- crates/clawhdf5-format/src/filters_szip.rs | 212 ++++++++++++++---- .../tests/fixtures/filters/README.md | 3 + .../tests/fixtures/filters/le_data.h5 | Bin 0 -> 72368 bytes .../tests/fixtures/filters/noencoder.h5 | Bin 0 -> 8088 bytes .../tests/fixtures/filters/szip_h5py.h5 | Bin 0 -> 11320 bytes 5 files changed, 177 insertions(+), 38 deletions(-) create mode 100644 crates/clawhdf5-format/tests/fixtures/filters/le_data.h5 create mode 100644 crates/clawhdf5-format/tests/fixtures/filters/noencoder.h5 create mode 100644 crates/clawhdf5-format/tests/fixtures/filters/szip_h5py.h5 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"), + // %-xy%vXYjAc^Eb-%m&P4+hT z?j4sW>dbfH-oD-MXTQ7q$-8%-@4k2M9HJ*VW4jQ+{G}|xcd7&qjKv$Ie`g03&k9#Ruou$r$M`kb0Y^-(v3tjG6%Ur%=27v1XHaWy@+ z<0>B5wX2hpq90dS2V9LDSNku6GGmwaxI#X{>Nv>%eKw8H5mJ~<4Y8@)#XfwpsP|bx zFaY-By(MuljpH+iF0pO`+bqffHI5?DuAFUI-PF2v)4FwA)^3GnBEcbD!e7{4*2*&a z$2zZknEzlwqxuAfAwd^nI_CQtKCnZ!ok%~Em7m(M;Dw>b)Q2Y%pqNu*aqvhoK;oQe zi==tWc;=L-fssB<)JP}VUZWZ&&rhqS=OE{0;i}Vd2GOM@>vO{Kc0s`Nl+_>Ncz=Z! z2fOua;y5$}cyoMNaN_|F$>J|B93Mxiz`^+*@VVrXgl#u+aOAYg&+!732UK~`aXNX~ z;x9EhPqfn-&eWD|S6E<&96moV z=T(;T_N}|!aPX;@>2=S*N0-abA&dL6xo?@roE9q+9_mi?O`2NkII&0#e?|G340EW1 zGUN%{#+*B4?>zU<4!sO^sJAM#ocj@F$hq-tr}9^BxOhqL@l3CQd&=*(RC_TY?nmS; zPAn51>iSMBoY%hl_QNeZw>9_M_ZJlv6_5K?lKp0R_ZttsX-zL?t5cmgpicpmesGGNoH;Ka<6i8D_zS1H^#BwUImwXbKlx!_EHKT%Uh-i40bJErW+pmzV zKzf1<6%y234+PoD*8?Gqc@@v+)dxBHk*!!S$Tu!H7L<;>;(`wq&MPiBq~diGhk74} z{hkrmV;q1uEF+$BVt^PR28aP-fEXYKhyi~Xz__4VU7yeYJ^_#WCir~&YR zKY6Uf>-?F0P?;DY28aP-fEXYKhyh|Ce+GQh+g|r6E4>Z(?eWxM-lq(8FpIJ3ZCsBZ zgWf(gThf_TZ^OFqD2+;Y%JCbqv)^GT?~r;gX|MdGq@?td(0vuVCzc)_zvru!l~=I8 zUq5AoVXj}}b1HYJ4#L-D z@2Odt_nS9{7|rTb`59Ay8@Os!-ce*2zGOwW`)oX>&z?5o3Y@kXSL`=%ZgI-1E;!TQ z2Oe>XI9PBnloJEQ05L!e5Cg;jF+dFX#(-~e3cn5@HCyI5C8ZYX^*9CU#u}##UkAuP zPKmCTY~po*jB!f7`g+LF!}y0Y^>AgWw4h#h)7RIy;7otFDIE5P{t!p{8kwP-7$63S z0b+m{AO?s5e;M#iUx(BKC0SSuP`?fU+HQzHupJDjudyy)eGTo?uLDT`&_^ywyAIH| zQ8oaVAU`1=_4`-(>g#oe9tPi#Q0K}LV^ISq27yK`?{LAH{tkNRYd8N|SF}qE5Cg;j zF+dCu1H=F^kOu?4>1)1zEj3$MjKb8{TwaS6Y+YYN-B|VY$oh5O`kHN%yut8Y^|gdB zHn**<>E4Fcwdv(m!{vJT{mcgW1DXx%S~oRqZCKN^wQ0D8n;ku*+1G5kYZdpVJ8-ot z2Wj@z%}s4-tr>H%WB2LyEzO(SvW;J}qk}a6>o*#$jsA!2*xInSZrXTfgV+pHkmld8 zgO{6I*KS?iytd8cK2ts5=dE4W+O&183C~n723ddK`pu)j{9rKaaYpsccdc8;|AywF zkfB~QS-5Kb)=^;Y->gUf7-QmnSa@8TdA{-CidoVVUROcCW<1ws;N0T3pSj>ne=jH; z_J{ruhh?-=P7DwO!~iis3=jjv05RYX0~jBMmC(HS-ytXy*FoSpLi|5Zme}g?A=Jb5 z7fAZ^P_YggpsD+SJJM-?0;B-F-{cIar-_-o$9pm-t0Z#?>ijYkR7)} z`Rji&JFZ#b<~)=g_eHnn8?VHT_yXsNaRZi58929i;$s(_>2IQXST6cUe~6REK3XX! z28aP-fEXYKhyh}N81Rb$-{KX1JxgjfOAO77S7fx38q(txtb+z)jaNX%!)pxjO65)% zh`cNZGBjhnlCRFLR(jX#y1@4g9O_~*GHvq%7aZzmgx*a$8>yh27$63S0b+m{AO?tm zOUgii%So}+tk#<^{<=V=Xg@hTZTP;isUbFXyPAcrw?Q)=?=6XgNxy##>L#$T@Dotu zC}8sY)BK(raUaku^*khP9lx;?9F6E=AxX`3e8#FJD#xvxjyB}#C~qg!_aIE1se6MH zWmCkt>FDcRaHhW&g~I`%Kg3}f?UWM(!~iis3=jjv05L!e_``s2I(m9;Iyzp`dmKit z>u9LULr3q8NeVzJf@~C$@$=Hrzc6&K9NOsbF~4Ep+;ncyBc+CSO~RKKFY`0Z7HqTt`p$ zp!;zpcugGYS!}*h>E5uaaMQVuDcm#zXZq_V9gWmdP7DwO!~iis3=jjvz$InCHy!Qs zdUOstdd{CEQXTy@%l33DdIdWmiABMV972m9Q?#cQU zi%w7e3BZ#*H9a>i+gHYJTK4?1Jx=25tT>r^A5EMuwBcVd7uRjxi1s{~Vh)ORxm1U==kLHC7~IMJ$yh__#^l-)W7(G7BJt9-#Oq_z$E?pqpQAoDea!k?^f~Hd)5omO zMeM6F>toZ$tj|TC;|Ps{ z7!x21kbE;f7N2!v!K{%D6QH^F(lW+3Rin(GXWx8ol!yUhfEdV$0pHfq#(X_G`#Re1 z_-tTJD9m5RI@+;5UdI!dK~0FIqtV}!3UiZp%6A=Yjatw0T1R_I;k?$-UQsx&b+n|y z;mmN|4A(S}Y>k5kJhy~ad`8#=19P=q_ z*P&OOl@xgB=Zkw@SN+08^8)55jXUByV z?$G<$aUq4<+?yTew>95*r0g9z5?IvK6UGHtXBjxRcww^(&h*!*aOfZXAr8xEr<@od z28aP-fEXYKhyh~29|nAjN5;Iaf6?(s=XptiuklE}I{I}(_u^j%^(>azDY=Jw+fC=r zbitYa8Wj%xqd&x98SRu41H=F^KnxHA!~iis4EV!oKEs`e_};$Lc&8Fr&;eZ6^-YLEc~Rm}=_ zGMF8AOtoKpi;E*wy;tSuYO~|~wz|b5cwD&k_qd7+OdQ4q*sj~a;qicW%ul%BOn+x7 z9zn_|CkBWCVt^PR28aP-;F2=nTRbv7w|Jx|uqz4^)Z-DT%O@TI#hxJ^X)lvfKCgHr zUmacYF3gGzKkrSJWI%-bkY8=!+;sFUE;v(1uU9zi1N|Wm%V?*Z7$63S0b+m{AO?s5 JV!$5;{ts@6(^>!k literal 0 HcmV?d00001 diff --git a/crates/clawhdf5-format/tests/fixtures/filters/noencoder.h5 b/crates/clawhdf5-format/tests/fixtures/filters/noencoder.h5 new file mode 100644 index 0000000000000000000000000000000000000000..84f87526ffe4943f03299168c5ff17fdbec94706 GIT binary patch literal 8088 zcmeHKF;Buk6n@uMQy?lD6OD^EFm!Uzn6QZn7+BP(0~;n#QjGy4xHz#mxH2(1@)L}^ zt8sR6)L)_ZdiS815E2px>3hlD>wEY0uJ_&b`t_lhKUhd^BmvEJ!GSc`VW6vTW4+DM zm>`c)FeM|tSlB-fL!pyJ)7CMbT%+kXs?C~zQ)%5@Tvuu>zr7Xg0{qIig4=ri%CD%> zt_BKCK+`=uQHTG;pPdv600M|-^7Qy0j=qwVU>_?61XyZ!07x_IW%Yhy+%X&e+~mm# zaKm<80n8b^WrHwVDTF!(M>*Tz(k`eo08NIqfOYMK9&j~`V6b0P;Dtk(W_2%)?*P8V z2<@wa{Yh2uka1Oi*pqEp-n{b4FinpED++;zn6J;596%PssqOUeD32WKUIE@FoLAC} zR9>8V$eFwm&ncl+h0Z7P=`pkIt&KM54@_jCN1Hif2ABb6fEi#0n1LB!;IDNk=?T&y zDJX}UQ5bzUNr$2$L_K4ECel;r)UBjAHF%tUB0a59 z!@AG9B7j0 literal 0 HcmV?d00001 diff --git a/crates/clawhdf5-format/tests/fixtures/filters/szip_h5py.h5 b/crates/clawhdf5-format/tests/fixtures/filters/szip_h5py.h5 new file mode 100644 index 0000000000000000000000000000000000000000..f197d796871fe81bf8eaea6e29ab3e46b40009fb GIT binary patch literal 11320 zcmeHMO=uHA6n^`YW@)sW;!jJ7?N%cQS}0PXhhkTnlz8xm+Rz>f=1_zl1VN;+muy8e zR)kvhDm{3p0X>il$&sU zVeT4QJQo$`a7D#AEjn3-Zsa3kEDD|!DGG*}pxvqF()kN14UFI!(5FM|XyjkrCgrY1 z8j{~l;;q9!*ZWtbc;9KD*W7orO%_8Y1D?acvukK^6aEz5Dxk6HXW^?Ci=LR)F_{n6~B zwe-a7(u=cQ%iy?+x@pC02pTP|R}K`m_dPR{k5_K5JSji@(4YKtHCrr}-a>Zz{;Sgb z{M2ixYmUruu?iPDpckAr!~X^XtuT#%y(zQ+Rda*hNQWb0@t)|SUymMyKHo!aE{p&p zzz8q`jKCj4;Men!aljBFE|Q3tN0gaFfL8O64}>M=lw46mc}Smfc_^;RLk8(v!P?o< zD35lk$uu*R9~&b*eAs&BA?MNR!NNR7fDvE>7y(9r5nu#%9Rjs+w2WFc`tXj4AuNg- zB8e2l&4Q-~5*Zk_8aJ!&_h@}Y)F0pPbt@V84cKyuV4RcO-`qBeqmMV3_anu{10=__ fRF_I9ZdUTDg;lF_pAldL7y(9r5nu#%0RrCvf