[dirkboye]

Add the Terse HDF5 filter (registered filter ID 32029) from
TERSE/PROLIX v2.0 (TRPX, https://github.com/senikm/trpx2.0), embedded
as a git subtree in lib/trpx2.0 (MIT license).

TERSE/PROLIX is a fast, lossless compression algorithm for integer
greyscale images with a high dynamic range, such as X-ray and
electron diffraction patterns. It encodes the data in small blocks of
values (12 by default): each block is stripped of its unused most
significant bits and is preceded by a short header giving the number
of bits per value, which is reused from the previous block when
unchanged. Version 2.0 compresses in parallel and adds compression
modes for unsigned data with overloads and for very weak data.
See S. Matinyan et al., Ultramicroscopy 281 (2026) 114298,
https://doi.org/10.1016/j.ultramic.2025.114298.

The filter supports signed and unsigned 8-, 16- and 32-bit integer
datasets. It takes no user options: its set_local callback stores the
dataset's data type code in cd_values. hdf5plugin.Terse provides the
compression arguments for h5py.

The filter is written in C++20, so the build now also handles the C++20
compile flag for filter extensions (previously only for static
libraries).

[kif]

Hi, I tried the plugin and found that recompressing uint16 data (this dataset: https://www.silx.org/pub/pyFAI/kevlar/) was randomly wrong ... Same for uin32. Signed datasets looks OK. I believe this is a killer for a lossless compression plugin.

Performances are also 50% worse the bslz4 when it comes to compression ratio and 100% worse for the read-speed, but a laptop is not a decent benchmarking tool.

Claude/Opus provided a minimum reproducing example:
terse_bug_repro.py

[kif]

Here is also the report (AI-generated):

Terse filter is not lossless for uint16/uint32 chunks larger than 2¹⁸ elements

Testing this PR (commit f1e3e45, branch add-trpx-filter) on Eiger 4M data from ESRF-ID13, I found
that reading back a uint16 chunk does not always return the data that was written. A few pixels
come back with wrong values, and the number of wrong pixels changes from one read to the next.

Environment: Linux x86_64 (Debian 13, 4 cores), Python 3.13.5, h5py 3.16.0, HDF5 2.0.0,
numpy 2.5.3. hdf5plugin was built from this branch with pip install ..

Reproducer

import numpy, h5py, hdf5plugin

S = 1 << 18  # TERSE_DEFAULT_CHUNK_SIZE in lib/trpx2.0/tersecodec/src/filter.cpp
rng = numpy.random.default_rng(0)
frame = rng.poisson(5.0, size=(2167, 2070))  # Eiger 4M sized frame, weak Poisson signal

with h5py.File("terse_bug.h5", "w") as f:
    for dtype in ("uint8", "int8", "uint16", "int16", "uint32", "int32"):
        data = frame.clip(0, 127).astype(dtype)
        ds = f.create_dataset(dtype, data=data, chunks=data.shape, compression=hdf5plugin.Terse())
        nerr = [int(numpy.count_nonzero(ds[()] != data)) for _ in range(5)]  # read the same chunk 5 times
        bad = numpy.flatnonzero(ds[()] != data)
        print(f"{dtype:>6}: mismatching pixels over 5 reads {nerr}; flat index mod 2**18 - 2**18: "
              f"{sorted(set((bad % S - S).tolist()))}")

Output:

 uint8: mismatching pixels over 5 reads [0, 0, 0, 0, 0]; flat index mod 2**18 - 2**18: []
  int8: mismatching pixels over 5 reads [0, 0, 0, 0, 0]; flat index mod 2**18 - 2**18: []
uint16: mismatching pixels over 5 reads [7, 7, 7, 6, 6]; flat index mod 2**18 - 2**18: [-4, -3, -2, -1]
 int16: mismatching pixels over 5 reads [0, 0, 0, 0, 0]; flat index mod 2**18 - 2**18: []
uint32: mismatching pixels over 5 reads [7, 6, 6, 6, 6]; flat index mod 2**18 - 2**18: [-4, -3, -2, -1]
 int32: mismatching pixels over 5 reads [0, 0, 0, 0, 0]; flat index mod 2**18 - 2**18: []

The values are not clipped, since Poisson(5) never reaches 127, so every dtype receives the same
integers.

Observations

• Only unsigned 16- and 32-bit data are affected. On the same values, int8, int16 and
  int32 round-trip exactly, and so does uint8.

• Compression is deterministic; decompression is not. Five writes of the same frame give
  byte-identical chunks (checked with read_direct_chunk). Reading one stored chunk several times
  gives 6 or 7 wrong pixels depending on the read.

• The wrong values are always the last partial block of a sub-frame. filter.cpp splits each
  chunk into Terse frames of TERSE_DEFAULT_CHUNK_SIZE = 1 << 18 values. Because
  2¹⁸ = 262144 ≡ 4 (mod 12), with the default block of 12 values each sub-frame ends with a partial
  block of 4 values. All wrong pixels are at flat indices k·2¹⁸ − 4 … k·2¹⁸ − 1. Setting those last
  4 values of every sub-frame to 0 makes the errors disappear.

• The chunk size matters. In 1D tests with Poisson(5) uint16 data:

  Number of elements	Sub-frames	Result
  up to 3·2²⁰ + 12	up to 12	always correct
  4·2²⁰ + 4	16	reproducibly wrong

  A full Eiger 4M frame (4.49 M pixels, 17 sub-frames) is affected.

• Constant data is not affected. All-zero, all-one and constant-7 data round-trip correctly.

• Real detector data is affected badly. On real Eiger 4M uint16 frames, which contain many
  65535 gap/masked pixels, up to 1.8 M pixels out of 4.5 M are wrong after a round-trip. 18 of the
  20 frames tested were corrupted. The same frames converted to int16 (bitwise view) or int32
  were lossless on 200 frames.

Possible cause

This is only a guess. The read-time non-determinism, together with the fact that only multi-frame
chunks fail, points to Terse<Concurrent>::prolix(container), which decodes sub-frames in parallel
into one shared buffer. The unsigned decoders (f_prolix_unsigned, f_prolix_small_unsigned*)
could read or write past to on the last partial block, for example by working on whole blocks of
d_block values. Neighbouring sub-frames would then race on those few values.

Two quick checks could confirm this:

• set TERSE_DEFAULT_CHUNK_SIZE to a multiple of 12;
• decode the sub-frames sequentially.

Suggestions for the PR

• Add a round-trip test on a large, realistic unsigned array to src/hdf5plugin/test.py, for
  example a Poisson-distributed uint16 array of at least 4·2²⁰ elements with a few 65535 values.
  The current test uses numpy.arange(256**2), which is a single sub-frame, so it cannot catch
  this.
• Report the bug upstream in senikm/trpx2.0 if it is in Terse.hpp rather than in filter.cpp.

[dirkboye]

Thank you for the review!

I've created 2 merge requests into https://github.com/senikm/trpx2.0/ which should fix the issues:

senikm/trpx2.0#3
senikm/trpx2.0#4