libaec / szip support #420

mkitti · 2023-02-21T04:04:28Z

NASA Earth Observing System (EOS) uses szip compression:
https://www.earthdata.nasa.gov/esdis/esco/standards-and-practices/hdf-eos5

It is an integrated compression codec in HDF5:
https://portal.hdfgroup.org/display/HDF5/Szip+Compression+in+HDF+Products

Deutschen Klimarechenzentrum (DKRZ) has a freely available implementation available here under a 2-Clause BSD License:
https://gitlab.dkrz.de/k202009/libaec
Github Mirror
https://github.com/MathisRosenhauer/libaec

Conda-forge has libaec binaries here:
https://anaconda.org/conda-forge/libaec/files

They are produced by the following feedstock:
https://github.com/conda-forge/libaec-feedstock

Based on the outstanding use of SZIP and the availability of the libaec implementation, I recommend that numcodecs support SZIP / AEC compression.

mkitti · 2023-02-21T04:09:05Z

Also szip support via libaec was added to conda-forge's build of HDF5 in conda-forge/hdf5-feedstock#179

martindurant · 2023-02-21T14:16:07Z

I don't know about the merits of if sz versus any other compressor, but if there is data out there using it, then at least kerchunk would like numcodecs to be able to deal with it. So we need a little cython code to link with the available binary? We will probably not want numcodecs to explicitly depend on libaec.

martindurant · 2023-03-05T18:04:35Z

Although I can see szlib.h is included with libaec via conda-forge, there are no instructions on how to call the function that seems to be the one we want, SZ_BufftoBuffDecompress .

mkitti · 2023-03-05T18:31:10Z

Note that there are two shared libraries

libaec.so
libszip.so

I think you want the second one?

$ objdump -T libaec.so

libaec.so:     file format elf64-x86-64

DYNAMIC SYMBOL TABLE:
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) free
0000000000000000  w   D  *UND*	0000000000000000  Base        _ITM_deregisterTMCloneTable
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) memset
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) calloc
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) memcpy
0000000000000000  w   D  *UND*	0000000000000000  Base        __gmon_start__
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) malloc
0000000000000000  w   D  *UND*	0000000000000000  Base        _ITM_registerTMCloneTable
0000000000000000  w   DF *UND*	0000000000000000 (GLIBC_2.2.5) __cxa_finalize
0000000000002ad0 g    DF .text	000000000000008d  Base        aec_encode
0000000000005d40 g    DF .text	0000000000000072  Base        aec_decode
0000000000002bc0 g    DF .text	0000000000000099  Base        aec_buffer_encode
0000000000005980 g    DF .text	00000000000003bc  Base        aec_decode_init
0000000000005e2c g    DF .fini	0000000000000000  Base        _fini
0000000000005df0 g    DF .text	000000000000003c  Base        aec_buffer_decode
0000000000001000 g    DF .init	0000000000000000  Base        _init
0000000000002b60 g    DF .text	0000000000000057  Base        aec_encode_end
0000000000005dc0 g    DF .text	0000000000000026  Base        aec_decode_end
0000000000002790 g    DF .text	0000000000000338  Base        aec_encode_init


$ objdump -T libsz.so

libsz.so:     file format elf64-x86-64

DYNAMIC SYMBOL TABLE:
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) free
0000000000000000  w   D  *UND*	0000000000000000  Base        _ITM_deregisterTMCloneTable
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.4)  __stack_chk_fail
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) memset
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) calloc
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) memcpy
0000000000000000  w   D  *UND*	0000000000000000  Base        __gmon_start__
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) malloc
0000000000000000      DF *UND*	0000000000000000 (GLIBC_2.2.5) memmove
0000000000000000  w   D  *UND*	0000000000000000  Base        _ITM_registerTMCloneTable
0000000000000000  w   DF *UND*	0000000000000000 (GLIBC_2.2.5) __cxa_finalize
00000000000014f0 g    DF .text	0000000000000384  Base        SZ_BufftoBuffDecompress
00000000000065dc g    DF .fini	0000000000000000  Base        _fini
00000000000064f0 g    DF .text	0000000000000072  Base        aec_decode
0000000000001000 g    DF .init	0000000000000000  Base        _init
0000000000003370 g    DF .text	0000000000000099  Base        aec_buffer_encode
0000000000006130 g    DF .text	00000000000003bc  Base        aec_decode_init
0000000000003280 g    DF .text	000000000000008d  Base        aec_encode
0000000000006570 g    DF .text	0000000000000026  Base        aec_decode_end
00000000000010e0 g    DF .text	0000000000000405  Base        SZ_BufftoBuffCompress
0000000000001880 g    DF .text	0000000000000006  Base        SZ_encoder_enabled
0000000000002f40 g    DF .text	0000000000000338  Base        aec_encode_init
0000000000003310 g    DF .text	0000000000000057  Base        aec_encode_end
00000000000065a0 g    DF .text	000000000000003c  Base        aec_buffer_decode

Are you looking for documentation of the function?

mkitti · 2023-03-05T18:39:28Z

Here's the link to the HDF5 SZIP filter showing usage of SZ_BufftoBuffDecompress:

https://github.com/HDFGroup/hdf5/blob/7b833f04b5146bdad339ff10d42aadc416fb2f00/src/H5Zszip.c#L294-L309

mkitti · 2023-03-05T19:56:35Z

Does this help?

#include<stdio.h>
#include<stdint.h>
#include<szlib.h>


int main(int argc, char *args[]) {
    printf("SZIP demo test\n______________\n\n");

    // Pixel type
    typedef uint32_t pixel_t;
    // Status, check for errors
    int status;

    // INPUT DATA

    const size_t inbuf_len = SZ_MAX_PIXELS_PER_SCANLINE*4;
    const size_t inbuf_nbytes = SZ_MAX_PIXELS_PER_SCANLINE*4*sizeof(pixel_t);
    pixel_t inbuf[inbuf_len];

    for(int i = 0; i < SZ_MAX_PIXELS_PER_SCANLINE*4; ++i) {
        inbuf[i] = i % 16;
    }

    // SZIP PARAMETERS

    SZ_com_t sz_params;
    //little endian
    sz_params.options_mask = SZ_LSB_OPTION_MASK;
    sz_params.bits_per_pixel = sizeof(pixel_t)*8;
    sz_params.pixels_per_block = SZ_MAX_PIXELS_PER_BLOCK;
    sz_params.pixels_per_scanline = SZ_MAX_PIXELS_PER_SCANLINE;

    // COMPRESSION

    size_t compressed_buf_nbytes = SZ_MAX_PIXELS_PER_SCANLINE*4*sizeof(pixel_t);
    char compressed_buf[compressed_buf_nbytes];


    status = SZ_BufftoBuffCompress(
        compressed_buf, &compressed_buf_nbytes,
        inbuf, inbuf_nbytes,
        &sz_params
    );

    printf("SZ_BufftoBuffCompress status: %d\n", status);
    printf("inbuf_nbytes: %zd\n", inbuf_nbytes);
    printf("compressed_buf_nbytes: %zd\n\n", compressed_buf_nbytes);

    if (status != SZ_OK) {
        return status;
    }

    // DECOMPRESSION
    
    size_t decompressed_buf_nbytes = inbuf_nbytes;
    char decompressed_buf[decompressed_buf_nbytes];

    status = SZ_BufftoBuffDecompress(
        decompressed_buf, &decompressed_buf_nbytes,
        compressed_buf, compressed_buf_nbytes,
        &sz_params
    );

    printf("status: %d\n", status);
    printf("compressed_buf_nbytes: %zd\n", compressed_buf_nbytes);
    printf("decompressed_buf_nbytes: %zd\n", decompressed_buf_nbytes);
    
    if (status != SZ_OK) {
        return status;
    }

    return 0;
}

To compile and run:

gcc test_szip.c -lsz -Wall -std=c2x -o test_szip && ./test_szip

martindurant · 2023-03-05T20:23:04Z

Thank you, that is helpful.

Are those sz_params values guaranteed to be like that in some HDF file, or is that something stored in the HDF filter metadata? Some question for knowing how big the decompressed buffer should be.

mkitti · 2023-03-06T00:49:48Z

In HDF5, H5Z__set_local_szip configures some of the parameters (bits_per_pixel, pixels_per_scanline, options_mask) based on the the properties of the dataset:
https://github.com/HDFGroup/hdf5/blob/7b833f04b5146bdad339ff10d42aadc416fb2f00/src/H5Zszip.c#L106-L244

A HDF5 chunk stores the number of bytes of a decompressed buffer in the first four bytes of the chunk as a little endian uint32:
https://github.com/HDFGroup/hdf5/blob/7b833f04b5146bdad339ff10d42aadc416fb2f00/src/H5Zszip.c#L298-L300

The macro UINT32DECODE is defined here and assumes the uint32 value is little endian:
https://github.com/HDFGroup/hdf5/blob/7b833f04b5146bdad339ff10d42aadc416fb2f00/src/H5Fprivate.h#L199-L209

Note that the practice of storing extra metadata at the beginning of a chunk is common in HDF5. For example, the LZ4 filter uses the first 8 bytes to store the original size of the buffer followed by 4 bytes to encode the block size:
https://github.com/HDFGroup/hdf5_plugins/blob/5908c83a2e20c2d1ae5ce5fdf35953238619b719/LZ4/src/H5Zlz4.c#L111-L118

mkitti · 2023-03-06T00:55:34Z

Some of the parameters are also configurable by the user. See the documentation for H5Pset_szip:
https://docs.hdfgroup.org/hdf5/develop/group___d_c_p_l.html#ga37de4b6071a94574cfab5cd6de9c3fc6

martindurant · 2023-03-06T02:13:37Z

Some of the parameters are also configurable by the user.

I mean at read time, when these must be fixed.

Perhaps easiest would be for me to get my hands on a sz-containing hdf and see what info is available.

mkitti · 2023-03-06T06:18:24Z

It's pretty easy to create such a hdf5 file with conda-forge's h5py now.

conda environment:

mamba create -n h5py_hdf5_1_14 -c conda-forge hdf5=1.14 h5py ipython

In [1]: import h5py

In [2]: with h5py.File("test.h5", "w") as h5f:
   ...:     ds = h5f.create_dataset("test", (1024, 1024), chunks=(16,16), dtype="uint16", compression="szip", compression_opts=('ec',32))
   ...:     ds[:] = 1

In [3]: with h5py.File("test.h5", "r") as h5f:
    ...:     ds = h5f["test"]
    ...:     print(ds.compression)
    ...:     print(ds.compression_opts)
    ...:     dcpl = ds.id.get_create_plist()
    ...:     print(dcpl.get_filter(0))
    ...: 
szip
('ec', 32)
(4, 1, (141, 32, 16, 256), b'szip')

4 is h5py.h5z.FILTER_SZIP
1 is h5py.h5z.FLAG_OPTIONAL
141 is 1 | 4 | 8 | 128 or SZ_ALLOW_K13_OPTION_MASK | SZ_EC_OPTION_MASK | SZ_LSB_OPTION_MASK | SZ_RAW_OPTION_MASK
32 is pixels per block
16 is bits per pixel (from uint16)
256 is pixels per scanline. The chunk is only 256 elements.

You can use the h5ls utility to get the information as well.

$ h5ls -v test.h5 
Opened "test.h5" with sec2 driver.
test                     Dataset {1024/1024, 1024/1024}
    Location:  1:800
    Links:     1
    Chunks:    {16, 16} 512 bytes
    Storage:   2097152 logical bytes, 294912 allocated bytes, 711.11% utilization
    Filter-0:  szip-4 OPT {141, 32, 16, 256}
    Type:      native unsigned short

The order of the parameters come from here:
https://github.com/HDFGroup/hdf5/blob/7b833f04b5146bdad339ff10d42aadc416fb2f00/src/H5Zpublic.h#L165-L180

martindurant · 2023-03-06T14:50:29Z

I have installed h5py from conda-forge or from pip, but

h.create_dataset("sz", (100, 100), "i", compression="szip")

says

ValueError: Compression filter "szip" is unavailable

I have h5py 3.8.0 and libaec 1.0.6 (and I do see libsz.dylib in my lib/)

martindurant · 2023-03-06T14:53:31Z

(nvm, had to update hdf too)

cgohlke · 2023-03-13T04:54:50Z

FWIW, the next version of imagecodecs will include a SZIP codec. It already has an AEC codec.

mkitti · 2023-03-13T06:17:15Z

Thanks @cgohlke . Does that mean you have and will be working on a Cython based version?

martindurant · 2023-03-13T13:34:41Z

If imagecodecs is getting SZIP, I don't think we mind how it is implemented :)
That probably means it'll get released sooner than a new numcodecs would - any ETA @cgohlke ? It certainly solves everything for us, so I'll happily close this PR.

cgohlke · 2023-03-14T04:20:51Z

Does that mean you have and will be working on a Cython based version?

Yes. Imagecodecs is Cython based and includes numcodecs compatible codecs.

any ETA

It's ready for release but depends on libjpeg-turbo 3, which is currently in beta...

cgohlke · 2023-03-17T03:48:41Z

any ETA

Imagecodecs v2023.3.16 is available on PyPI (conda-forge will have to wait for libjpeg-turbo 3) and includes a numcodecs compatible SZIP codec based on an implementation in Cython

martindurant · 2023-03-17T13:14:45Z

Thank you, @cgohlke

martindurant mentioned this issue Mar 6, 2023

Add sz using ctypes #422

Closed

7 tasks

mkitti mentioned this issue Mar 16, 2023

Adding zarr-like "filters" to the n5 api saalfeldlab/n5#87

Open

martindurant closed this as completed Mar 17, 2023

mkitti mentioned this issue Jul 12, 2023

Add Jenkin's lookup3 as a 32-bit checksum for HDF5 #445

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

libaec / szip support #420

libaec / szip support #420

mkitti commented Feb 21, 2023 •

edited

Loading

mkitti commented Feb 21, 2023

martindurant commented Feb 21, 2023

martindurant commented Mar 5, 2023

mkitti commented Mar 5, 2023

mkitti commented Mar 5, 2023

mkitti commented Mar 5, 2023

martindurant commented Mar 5, 2023

mkitti commented Mar 6, 2023

mkitti commented Mar 6, 2023

martindurant commented Mar 6, 2023

mkitti commented Mar 6, 2023

martindurant commented Mar 6, 2023

martindurant commented Mar 6, 2023

cgohlke commented Mar 13, 2023

mkitti commented Mar 13, 2023

martindurant commented Mar 13, 2023

cgohlke commented Mar 14, 2023

cgohlke commented Mar 17, 2023

martindurant commented Mar 17, 2023

libaec / szip support #420

libaec / szip support #420

Comments

mkitti commented Feb 21, 2023 • edited Loading

mkitti commented Feb 21, 2023

martindurant commented Feb 21, 2023

martindurant commented Mar 5, 2023

mkitti commented Mar 5, 2023

mkitti commented Mar 5, 2023

mkitti commented Mar 5, 2023

martindurant commented Mar 5, 2023

mkitti commented Mar 6, 2023

mkitti commented Mar 6, 2023

martindurant commented Mar 6, 2023

mkitti commented Mar 6, 2023

martindurant commented Mar 6, 2023

martindurant commented Mar 6, 2023

cgohlke commented Mar 13, 2023

mkitti commented Mar 13, 2023

martindurant commented Mar 13, 2023

cgohlke commented Mar 14, 2023

cgohlke commented Mar 17, 2023

martindurant commented Mar 17, 2023

mkitti commented Feb 21, 2023 •

edited

Loading