wsi2zarr

Made with the assistance of Claude AI

Sorry for the vibe coding, but it helped me getting back to C. I started with implementations of the tile extraction within openslide code, but I had to implement it for each formats, working directly on the file bytes, not benefiting much from the openslide library. I've used the help of Claude AI to reorganize the code into this library. Further work will focus on building more tests to validate the conversion.

I'm at the stage of saying this is a prototype. If I can integrate this to other conversion library, this would be great.

Description

NOTE: Hopefully most information is correct about the formats. Go to openslide.org if you are seeking to understand them better. Read through if you want to understand about how this tool extracts the jpeg.

Extracts already-compressed tile/strip streams from TIFF-based WSI formats (Aperio SVS, Hamamatsu NDPI, generic tiled TIFF, ...) and writes them into a Zarr v2 layout without decoding/re-encoding pixels.:

• JPEG tiles: TIFF stores these in "abbreviated" form (quant/Huffman tables live once in the page's JPEGTables tag, tags 347). We splice the tables into each tile to get a standalone JPEG, then run a jpegtran- equivalent coefficient-level transcode via libjpeg-turbo (jpeg_read_coefficients / jpeg_write_coefficients — no iDCT, no color conversion) to optionally re-optimize Huffman tables, switch to arithmetic coding, or make the stream progressive. No new generation loss.
• JPEG2000 tiles (Aperio-style compression tags 33003/33005, or generic 34712): pure passthrough, byte-for-byte, no decode step at all.
• Zarr layout: each chunk file on disk is the standalone JPEG/J2K stream for that tile (openable directly by any JPEG/J2K decoder). The .zarray metadata records shape/chunks in pixel terms and a compressor id (imagecodecs_jpeg / imagecodecs_jpeg2k) so any zarr+imagecodecs reader can also decode the array as an ordinary array.

Stage 1 (generic tiled-TIFF dispatch for JPEG/JPEG2000) is done. Stage 2, the NDPI giant-strip case, is now in as well: a page stored as one enormous JPEG-compressed TIFF strip (RST0-RST7 restart markers standing in for tile boundaries) gets its own header/DRI parser, a single-pass restart-marker scan, batched reads (one read per output-tile-row, not per segment), and marker-renumbering reassembly into standalone per-tile JPEGs. JPEG-XL tiles and the VMS directory-of-JPEG case are still follow-ups.

How MCU offsets are found (updated after seeing OpenSlide's actual openslide-vendor-hamamatsu.c): NDPI ships precomputed-but-"unreliable" restart-marker byte offsets in private TIFF tags NDPI_MCU_STARTS_LOW (65426) / NDPI_MCU_STARTS_HIGH (65432). wsi_ndpi_open_page() reads these as a fast path (O(1) per tile, no scanning), validates every hint (each must really be preceded by a 0xFF <RST0-RST7> pair — that's the "unreliable" part, per OpenSlide's own naming), and only falls back to a single O(strip size) linear scan if the tags are absent or any hint fails validation. page->used_precomputed_offsets records which path ran.

Known limitation of the NDPI path: assumes the page is a single TIFF strip. Multi-strip giant-JPEG concatenation isn't implemented yet. (The earlier note about "sparse/holed segment grids" was a misreading on my part before seeing OpenSlide's actual compute_mcu_start — it asserts every tile index has real data, so that case doesn't exist.)

A property worth knowing about, not a bug: splitting one continuous JPEG bitstream into standalone per-tile JPEGs at the coefficient level is genuinely lossless -- no DCT coefficient is touched. But if chroma is subsampled (the common 4:2:0 case) and a decoder uses "fancy" (edge-aware) upsampling, pixels reconstructed near a newly-introduced tile boundary can differ slightly from decoding the same region as part of the original whole image, because the upsampler no longer sees the real neighboring block on that side. This is inherent to DCT-domain tiling of subsampled JPEG, not specific to this tool -- it's the same reason real tiled WSI viewers either accept it or decode with simple (non-fancy) upsampling. test_ndpi_reassemble.c disables fancy upsampling (do_fancy_upsampling = FALSE) precisely to get a bit-exact comparison against a whole-image decode; downstream consumers who care about this at tile seams should do the same.

Zarr chunk layout

Chunks are written as <root>/<level>/<t>/<z>/<row>/<col>/0:

• <level> is the bare pyramid level index (0, 1, 2, ... -- level 0 is full resolution), not level_0.
• <t> and <z> are the time point and focal-plane/Z-stack index. Every format currently extracted only ever has a single T and a single Z (both always 0), but the folder levels exist now so NDPI's own multiple focal-plane pages (read by OpenSlide's backend but not yet exposed this way here) have somewhere to go later without another layout migration.
• <row>/<col>/0 is .zarray's own dimension_separator: "/" nested- folder convention for a 3-D [height, width, samples] array (the trailing 0 is the always-singleton chunk index along the sample axis).

NDPI output tile size

--tile-size WxH (default 2048x2048) sets the requested output tile size for the NDPI giant-strip path (wsi_ndpi_extract rounds up to whole restart-interval segments, so the actual tile size may come out somewhat larger). Every other format currently passes its tiles through at their native size rather than re-tiling them, so this option is silently ignored there for now.

Code organization

include/            public headers
  wsi_common.h       shared status/codec enums
  wsi_debug.h         debug logging macros (WSI2ZARR_DEBUG=1 env var)
  tiff_source.h       libtiff wrapper: page introspection + raw tile read
  jpeg_tile.h         JPEGTables splice + libjpeg-turbo coefficient transcode
  jp2_tile.h          JPEG2000 stream classification/passthrough
  zarr_writer.h       Zarr v2 group/array/chunk writer
  ndpi_reassemble.h   NDPI giant-strip header/DRI parse + segment reassembly
  wsi_extract.h       ties the above together per-page / per-file
src/                 implementations of the above, plus main.c (CLI)
tests/
  test_unit_edge_cases.c    pure-logic tests (JP2 classify, splice edge cases)
  test_extract_pipeline.c   end-to-end: builds a synthetic tiled-JPEG TIFF
                            fixture with libtiff itself, runs extraction,
                            verifies every output chunk decodes to
                            byte-identical pixels vs. libtiff's own decode
                            of the original tile. No network / sample files
                            needed. Covers default, --progressive and
                            --arithmetic transform variants.
  test_ndpi_reassemble.c    end-to-end: encodes one complete JPEG (with a
                            restart interval, default 4:2:0 subsampling) via
                            raw libjpeg, drops it into a single TIFF strip
                            with TIFFWriteRawStrip (the actual NDPI shape),
                            runs the reassembly, and checks every output
                            tile against a crop of a direct decode of the
                            original whole-image JPEG. Dimensions are
                            deliberately non-multiples of the MCU/tile grid
                            to exercise real-dimension edge cropping.

Build

Requires libtiff-dev and libjpeg-turbo8-dev (or any libjpeg providing jpeg_mem_src/jpeg_mem_dest/jpeg_read_coefficients, i.e. libjpeg >= 8 or libjpeg-turbo).

sudo apt-get install -y libtiff-dev libjpeg-turbo8-dev cmake build-essential
cmake -B build -DCMAKE_BUILD_TYPE=RelWithDebInfo
cmake --build build -j
ctest --test-dir build --output-on-failure

Run

Usage: ./build/wsi2zarr <input.svs|.ndpi|...> <output_zarr_dir> [options]
Options:
  --level N          only extract TIFF directory N (default: all)
  --tile-size WxH     preferred output tile size in pixels (default: 2048x2048).
                      Currently only honored by the NDPI giant-strip path --
                      other formats pass their native tiles through as-is
                      and ignore this option, for now.
  --no-optimize      disable Huffman optimize_coding on JPEG tiles
  --arithmetic       re-encode JPEG tiles with arithmetic coding
  --progressive      re-encode JPEG tiles as progressive JPEG
  --debug            verbose debug logging (same as WSI2ZARR_DEBUG=1)

./build/wsi2zarr input.svs output_zarr_dir/
./build/wsi2zarr input.svs output_zarr_dir/ --level 0 --arithmetic --debug

Debug logging can also be enabled via WSI2ZARR_DEBUG=1 in the environment instead of --debug.

Getting real sample files

openslide.org publishes real SVS/NDPI/etc. samples for exactly this kind of testing, e.g.:

• https://openslide.org/formats/aperio/ (Aperio SVS samples)
• https://openslide.org/formats/hamamatsu/ (NDPI samples — the giant-strip path is now implemented; see the known limitations above for what's still not handled, in particular multi-strip files and sparse/holed segment grids)

Testing the ouput with zarr

Go to the next section if you want the description of additional work I did to get my test Zarr.

Currently, jpeg and jpeg-2000 are not in the default codecs accepted by zarr. To enable the codecs, see https://github.com/d-v-b/zarr-jpeg (for zarr v2 only?).

You can use my hugging-face bucket to test first. But for performance these need larger tiles or sharding (jpeg-2000 extracted tiles are 256x256).

https://huggingface.co/buckets/TomTBT/wsi_zarr_jpeg

I couldn't build & install the library zarr-jpeg, so I just took what I needed, and registered the codecs with the name I have in my zarr, imagecodecs_jpeg & imagecodecs_jpeg2k

from numcodecs.abc import Codec
from numcodecs.compat import ensure_ndarray, ensure_contiguous_ndarray, ndarray_copy
from numcodecs.registry import register_codec
from imagecodecs import jpeg_decode, jpeg2k_decode
import numpy as np

class jpeg(Codec):
    """Codec providing jpeg compression via imagecodecs.
    """

    codec_id = "imagecodecs_jpeg"

    def __init__(self):
        super().__init__()

    def encode(self, buf):
        return None

    def decode(self, buf, out=None):
        buf = ensure_contiguous_ndarray(buf)

        if out is not None:
            out = ensure_contiguous_ndarray(out)

        tiled = jpeg_decode(buf)
        
        return ndarray_copy(tiled, out)


register_codec(jpeg)

class jpeg2k(Codec):
    """Codec providing jpeg2k compression via imagecodecs.
    """

    codec_id = "imagecodecs_jpeg2k"

    def __init__(self):
        super().__init__()

    def encode(self, buf):
        return None

    def decode(self, buf, out=None):
        buf = ensure_contiguous_ndarray(buf)

        if out is not None:
            out = ensure_contiguous_ndarray(out)

        tiled = jpeg2k_decode(buf)
        
        return ndarray_copy(tiled, out)


register_codec(jpeg2k)

import zarr
import matplotlib.pyplot as plt
z_jpg = zarr.open('https://huggingface.co/buckets/TomTBT/wsi_zarr_jpeg/resolve/ndpi_CMU-2_arithmetic.zarr/0', mode='r')
print(z_jpg["0"])
plt.figure()
plt.imshow(z_jpg["0"][0, 0, 5000:8000, 5000:8000, :])

Finalizing the Zarr from the output

Thing's I'll implement very soon because it was fun for test, but explaining it just sounds crazy.

The Zarr coming out of this library is not OME-Zarr. I used the NGFF converter to produce a valid OME-Zarr from the same slide (better skip the conversion of the WSI tile to lossless compression if you don't want to blow up your disk with x10 the file size). Then add (or replace if you did the mistake) the jpeg-zarr serie to your OME-Zarr. You need then to modify the .zattrs "coordinate transformations" to what the existing levels were. (NGFF converter creates more) Also, the C dimension must be changed to last. See this example https://huggingface.co/buckets/TomTBT/wsi_zarr_jpeg/tree/ndpi_CMU-2_arithmetic.zarr/0/.zattrs

Moreover, the tiles produced are not padded to the chunk shape, so Zarr will fail reading the tiles if they are of the wrong shape (could be fixed by https://github.com/zarr-developers/zarr-extensions/tree/main/chunk-grids/rectilinear, but it's not yet accepted)

Padding with jpegtran in bash:

• assuming all tiles have the same dimension across levels...)
• parallel processing with &
• better to use the same options -optimize or -arithmetic that were use for tile generation

ZARR_DIR=/home/tom/jpeg/ndpi_CMU-2-optimize.zarr
SERIE_IDX=0
END_L=$(ls $ZARR_DIR/$SERIE_IDX | sort -V | tail -n 1)

TILE_W=2048
TILE_H=2048
for l in $(seq 0 $END_L); do
  END_Y=$(ls $ZARR_DIR/$SERIE_IDX/$l/0/0 | sort -V | tail -n 1);
  END_X=$(ls $ZARR_DIR/$SERIE_IDX/$l/0/0/0 | sort -V | tail -n 1); 
  for y in $(seq 0 $END_Y); do
    jpegtran -crop $TILE_W"x"$TILE_H"+0+0" -optimize -outfile $ZARR_DIR/$SERIE_IDX/$l/0/0/$y/$END_X/0 $ZARR_DIR/$SERIE_IDX/$l/0/0/$y/$END_X/0 &
  done;
  for x in $(seq 0 $END_X); do
    jpegtran -crop $TILE_W"x"$TILE_H"+0+0" -optimize -outfile $ZARR_DIR/$SERIE_IDX/$l/0/0/$END_Y/$x/0 $ZARR_DIR/$SERIE_IDX/$l/0/0/$END_Y/$x/0 &
   done;
done