NIfTI-2, NRRD and NWB decoders (0.9.0)

CHANGELOG.md

@@ -1,5 +1,36 @@
 # Changelog
 
+## 0.9.0
+
+- NIfTI-2: the 540-byte header (`n+2` / `ni2` magic, 64-bit dims, double
+  pixdim) decodes through the same `nifti` decoder as NIfTI-1, so an axis past
+  32767 voxels reads correctly. A `niftiVersion` field (1 or 2) is added to
+  both. Verified against nibabel-written files, both byte orders.
+- NRRD: a decoder for `.nrrd` and detached `.nhdr` headers (teem, ITK, 3D
+  Slicer, the Allen CCFv3 atlases): sizes, dtype, encoding, endian, space,
+  space directions and derived per-axis spacing, origin, units, kinds, labels,
+  measurement frame, data file, `key:=value` pairs, and whether the header fit
+  the read. Verified on the real CCFv3 `annotation_25.nrrd` header and against
+  pynrrd.
+- NWB: `.nwb` files (Neurodata Without Borders, the DANDI archive format)
+  decode through the HDF5 reader as their own format: version, identifier,
+  session description, start time and id, institution, lab, experimenters,
+  keywords, publications, subject, devices, electrode groups and count,
+  imaging planes, acquisition series (type, shape, dtype, unit, rate),
+  processing modules and interfaces, unit count and columns, interval tables,
+  stimuli, embedded schema namespaces and extensions. Targeted lookups rather
+  than a walk. Verified against h5py on four DANDI files and pynwb output.
+- HDF5: a soft link in a v1 symbol-table group (an undefined object address
+  in the symbol node) no longer ends the listing of that node; the entries
+  after it are read. An NWB TimeSeries keeps its `electrode` link in front of
+  `starting_time`, so the sampling rate of every MIES series was unreadable.
+- HDF5: the global heap scan widens in steps (8 KiB, 256 KiB, 1 MiB, never
+  past the collection's own size) instead of a fixed 8 KiB window, so strings
+  on the large collections a session file builds are read.
+- HDF5: identical reads and group listings are memoized within one decode,
+  so repeated lookups through the same group no longer count against the
+  read cap; no byte is fetched that was not fetched before.
+
 ## 0.8.0
 
 - A3M: a decoder for multiple sequence alignments (`.a3m`, `.a3m.gz`), the

README.md

@@ -52,8 +52,16 @@ Decoders today:
 
 - **MRC / MRCS**: cryo-EM micrographs, movie stacks, tilt series, EMDB maps.
 - **NPY**: NumPy arrays. Cross-checked against numpy.
-- **NIfTI-1**: neuroimaging volumes, `.nii` and `.nii.gz`. Cross-checked against
-  nibabel. Big- and little-endian.
+- **NIfTI-1 and NIfTI-2**: neuroimaging volumes, `.nii` and `.nii.gz`. Both
+  header layouts (348 and 540 bytes), big- and little-endian; `niftiVersion`
+  says which. Cross-checked against nibabel.
+- **NRRD** (`.nrrd`, detached `.nhdr`): the volume format of teem, ITK, 3D
+  Slicer and the Allen Institute's mouse brain atlases (CCFv3). Version,
+  dimension, sizes, dtype, encoding, endianness, space, space directions and
+  the per-axis spacing derived from them (or `spacings`), origin, units, kinds,
+  labels, measurement frame, the data file of a detached header, and the
+  `key:=value` pairs Slicer writes. Verified on the real CCFv3 annotation
+  header and cross-checked against pynrrd.
 - **CryoSPARC `.cs`**: a structured-array dataset. Reports record count and the
   field schema.
 - **SerialEM `.mdoc`**: the acquisition record SerialEM (and Thermo Fisher
@@ -136,6 +144,18 @@ pointers, see "Pointer-based formats" below):
     x), dtype and compression. Two Range requests describe a 51 GB file.
   - **Luxendo `.lux.h5`**: the generic structure labelled, with the `Data`
     block shape. No public sample was available; shape only.
+  - **NWB `.nwb`** (Neurodata Without Borders, the format of the DANDI
+    archive): version, identifier, session description and start time,
+    session id, institution, lab, experimenters, keywords, publications, the
+    subject (id, species, sex, age, genotype, strain), devices, electrode
+    groups and electrode count, imaging planes, every acquisition series with
+    its type, shape, dtype, unit and rate, processing modules and their
+    interfaces, unit count and columns, interval tables with row counts,
+    stimuli, and the schema namespaces the file embeds, with any extension
+    (`ndx-*`) named. Targeted lookups, not a walk, so a session file with
+    thousands of objects still reads in a few hundred bounded reads.
+    Verified against h5py on four DANDI files from four labs and on
+    pynwb-written files.
   - **generic `.h5` / `.hdf5` / `.he5`**: root attributes (fixed and
     variable-length strings, numbers), groups, and datasets with shape, dtype,
     layout, chunks and filters, bounded and with the truncation stated.

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "scigantic-headers"
-version = "0.8.0"
+version = "0.9.0"
 description = "Decode scientific file headers (cryo-EM, light microscopy and light-sheet imaging incl. TIFF/OME-TIFF, CZI, HDF5/Imaris/BigDataViewer, LIF, ND2, Zarr/N5; tabular, flow cytometry, mass spec, sequencing, structural biology) into typed fields. Pure, zero-dependency, with bounded readers that follow internal pointers without downloading the file."
 readme = "README.md"
 requires-python = ">=3.9"

src/scigantic_headers/__init__.py

@@ -21,6 +21,7 @@
 from .lif import decode_lif
 from .mdoc import decode_mdoc
 from .nd2 import decode_nd2
+from .nrrd import decode_nrrd
 from .obf import decode_obf
 from .ome_xml import parse_ome_xml
 from .tiff import decode_tiff
@@ -88,7 +89,7 @@
     "decode_parquet", "decode_fcs", "decode_mzml",
     "decode_a3m", "decode_fastq", "decode_illumina_run", "decode_vcf", "decode_sam", "decode_pdb",
     "decode_mmcif", "decode_genbank", "decode_gff", "decode_bed", "decode_dicom",
-    "decode_czi", "decode_hdf5", "decode_ics", "decode_klb", "decode_lif", "decode_nd2", "decode_obf",
+    "decode_czi", "decode_hdf5", "decode_ics", "decode_klb", "decode_lif", "decode_nd2", "decode_nrrd", "decode_obf",
     "decode_tiff", "parse_ome_xml",
     "decode_bdv_xml", "decode_ome_xml", "decode_xml_sidecar",
     "decode_n5_attributes", "decode_zarr_json", "decode_zarray", "decode_zattrs", "decode_zgroup",

src/scigantic_headers/decoders.py

@@ -89,6 +89,13 @@ def __init__(self, budget: int):
         self._budget = budget
         self._reads = 0
 
+    @property
+    def reads_remaining(self) -> int:
+        """Reads left before `SEEK_MAX_READS`. A decoder with optional detail
+        (NWB's per-series descriptions) checks this so the detail never starves
+        the fields that come after it."""
+        return max(0, SEEK_MAX_READS - self._reads)
+
     def _charge(self, n: int) -> int:
         """Bytes allowed for a read of `n`: 0 once the read count or byte budget
         is spent, else min(n, remaining budget). Counts the read."""
@@ -467,48 +474,77 @@ def decode_npy_header(data: bytes) -> Optional[DecodedHeader]:
 register_decoder("npy", decode_npy_header)
 
 
-# ── NIfTI-1 (neuroimaging volumes, MRI/fMRI; not cryo-EM) ──────────────────
+# ── NIfTI-1 and NIfTI-2 (neuroimaging volumes, MRI/fMRI; not cryo-EM) ──────
 #
-# 348-byte fixed binary header. Endianness is detected by reading sizeof_hdr:
-# it must equal 348 in the file's byte order. Validated by the magic at byte 344
-# ('n+1\0' single-file, 'ni1\0' header/image pair). Only raw .nii is handled;
-# .nii.gz is gzip-wrapped and, like BAM, needs a decompress step first.
+# Two fixed binary headers share one decoder. Endianness is detected by reading
+# sizeof_hdr: it must equal 348 (NIfTI-1) or 540 (NIfTI-2) in the file's byte
+# order, and the magic confirms which. `.nii.gz` is handled by the gzip step in
+# the reader, which hands the inflated leading block to this function.
 #
+# NIfTI-1 (348 bytes; magic at the end):
 #   0    sizeof_hdr   int32   == 348
 #   40   dim[8]       int16   dim[0]=ndim, dim[1..]=sizes
 #   70   datatype     int16   NIfTI datatype code
 #   72   bitpix       int16   bits per voxel
 #   76   pixdim[8]    float32 pixdim[1..3]=voxel sizes in mm
+#   108  vox_offset   float32
 #   344  magic        char[4] 'n+1\0' | 'ni1\0'
+#
+# NIfTI-2 (540 bytes; magic right after sizeof_hdr, 64-bit dims and doubles,
+# so a volume may exceed 32767 voxels on an axis):
+#   0    sizeof_hdr   int32   == 540
+#   4    magic        char[8] 'n+2\0\r\n\032\n' | 'ni2\0\r\n\032\n'
+#   12   datatype     int16
+#   14   bitpix       int16
+#   16   dim[8]       int64
+#   104  pixdim[8]    float64
+#   168  vox_offset   int64
 
 _NIFTI_DTYPE = {
     2: "uint8", 4: "int16", 8: "int32", 16: "float32", 32: "complex64",
     64: "float64", 256: "int8", 512: "uint16", 768: "uint32",
     1024: "int64", 1280: "uint64", 1792: "complex128",
 }
 _NIFTI_MAGIC = (b"n+1\x00", b"ni1\x00")
+_NIFTI2_MAGIC = (b"n+2\x00\r\n\x1a\n", b"ni2\x00\r\n\x1a\n")
+NIFTI1_HEADER_BYTES = 348
+NIFTI2_HEADER_BYTES = 540
 
 
 def decode_nifti_header(data: bytes) -> Optional[DecodedHeader]:
-    if len(data) < 348:
+    if len(data) < NIFTI1_HEADER_BYTES:
         return None
-    # Detect byte order: sizeof_hdr reads as 348 only in the correct one.
-    order = None
+    # Detect version and byte order together: sizeof_hdr reads as 348 or 540
+    # only in the correct one.
+    version = order = None
     for bo in ("<", ">"):
-        if struct.unpack_from(bo + "i", data, 0)[0] == 348:
-            order = bo
+        size = struct.unpack_from(bo + "i", data, 0)[0]
+        if size == NIFTI1_HEADER_BYTES:
+            version, order = 1, bo
+            break
+        if size == NIFTI2_HEADER_BYTES:
+            version, order = 2, bo
             break
-    if order is None or data[344:348] not in _NIFTI_MAGIC:
+    if order is None:
         return None
 
-    dim = struct.unpack_from(order + "8h", data, 40)
-    (datatype,) = struct.unpack_from(order + "h", data, 70)
-    pixdim = struct.unpack_from(order + "8f", data, 76)
+    if version == 1:
+        if data[344:348] not in _NIFTI_MAGIC:
+            return None
+        dim = struct.unpack_from(order + "8h", data, 40)
+        (datatype,) = struct.unpack_from(order + "h", data, 70)
+        pixdim = struct.unpack_from(order + "8f", data, 76)
+    else:
+        if len(data) < NIFTI2_HEADER_BYTES or data[4:12] not in _NIFTI2_MAGIC:
+            return None
+        (datatype,) = struct.unpack_from(order + "h", data, 12)
+        dim = struct.unpack_from(order + "8q", data, 16)
+        pixdim = struct.unpack_from(order + "8d", data, 104)
     ndim = dim[0]
     if not 1 <= ndim <= 7:
         return None
 
-    shape = list(dim[1:1 + ndim])
+    shape = [int(d) for d in dim[1:1 + ndim]]
     # _finite guards against uninitialized / garbage pixdim decoding to NaN/inf.
     voxel = [_finite(round(pixdim[i], 4)) for i in range(1, min(1 + ndim, 4))]
     dtype = _NIFTI_DTYPE.get(datatype, f"code{datatype}")
@@ -519,10 +555,11 @@ def decode_nifti_header(data: bytes) -> Optional[DecodedHeader]:
     tail = f", {mm} mm" if mm else ""
     return DecodedHeader(
         format="nifti",
-        summary=f"NIfTI-1 volume {dims}, {dtype}{tail}",
+        summary=f"NIfTI-{version} volume {dims}, {dtype}{tail}",
         fields={
             "ndim": ndim, "shape": shape, "dtype": dtype,
             "voxelSizesMm": voxel, "datatypeCode": datatype,
+            "niftiVersion": version,
         },
     )