feat(io): Native 5D IO + True 3D Out-of-Core Featurization

[darkclad]

Status: feature-complete and validated. This document is the consolidated, self-contained
write-up of the feature — design, what shipped, test results, and usage — including the two
out-of-core design notes as appendices.

This work delivers two related-but-distinct pieces of work that grew out of the same effort:

1. Native 5D IO — Nyxus reads X,Y,Z,C,T directly from OME-TIFF and OME-Zarr metadata instead
   of assuming a single 2D plane / hard-coded axis order.
2. True 3D out-of-core (OOC) featurization — a disk-streaming path for every 3D feature
   family (previously 3D only worked if the whole volume fit in RAM), discovered to be necessary
   once native 5D reading exposed that the volumetric pipeline had no real streaming path.

---

1. Design

1.1 One descriptor, many producers

All axis knowledge funnels through a single format-agnostic struct, OmeAxes
(src/nyx/ome/ome_axes.h): a dimensionOrder string (e.g. "XYZCT") plus a storageAxes
vector carrying kind/label/size/physical-size/unit per axis, independent of on-disk order,
with a storageIndexOf(char label) lookup. Sizes are resolved by axis role, not by fixed
position — this is what makes TCZYX, TZCYX, CZYX, ZYX, and YX all correctly readable
through the same code path, replacing the old hard-coded "always TCZYX, always plane 0" reader.

OME-TIFF populates OmeAxes from the OME-XML DimensionOrder/TiffData block; OME-Zarr
populates it from the NGFF axes + multiscales metadata. Everything downstream — the loader
facades, the 3D consumer (phase2_3d.cpp), the output schema, SlideProps — only ever consumes
OmeAxes, never format-specific metadata directly.

1.2 Design principles

• Additive, not disruptive. Existing 2D/2.5D/NIfTI paths are untouched; new behavior only
  triggers when OME metadata is actually present. Every phase was independently shippable.
• Symmetry across the two loader stacks. Every interface change (e.g. axis-aware
  loadTileFromFile(...,channel,timeframe,level)) was applied to both AbstractTileLoader and
  RawFormatLoader, and to all three format-dispatch sites, to avoid the build fragmenting into
  two half-migrated stacks.
• Channel is a first-class but optional axis. Output defaults to per-channel
  featurization (one row per ROI × timepoint × channel) to preserve information rather than
  silently collapsing channels.
• Physical calibration flows from metadata, overridable by CLI. Precedence is: explicit
  CLI --anisox/y/z > OME PhysicalSize* (opt-in via --use-physical-spacing) > isotropic
  1.0. Reconciled in one place (resolve_slide_anisotropy).
• Zarr version isolation. Zarr 0.4 (v2) keeps the existing vendored z5 library
  unconditionally; Zarr 0.5 (v3 + sharding) support was added behind the same code path once it
  was discovered z5 3.0.0 already reads v3 natively — the original design called for a new
  tensorstore dependency for this, which turned out to be unnecessary.

1.3 Chunking convention

Every OME-Zarr fixture in this project follows one Z-plane per chunk. This is what makes
plane-by-plane out-of-core streaming (stream_volume_planes) possible: Nyxus reads exactly one
Z-plane per tile request. A store chunked with Z-depth > 1 per chunk is still supported for
in-RAM reads (fixed in a bug described in §3.1), but a store where a chunk spans multiple Z
AND cannot be read plane-by-plane is correctly refused for OOC streaming with an explicit
error ("the volume is not delivered plane-by-plane") — the same category of guard used for
genuinely unstreamable formats like NIfTI (see §1.5).

1.4 Anisotropic resampling

Nyxus's anisotropic path works by resampling (nearest-neighbor duplication along the
coarser axis), not by carrying absolute physical units through every calculation. So when
--use-physical-spacing is used, the OME PhysicalSizeX/Y/Z values are converted to a ratio
normalized so the minimum axis = 1 (e.g. Z spacing 4× coarser than X/Y becomes a resampling
factor of 4, not "4.0 micrometers"). Physical units across axes are canonicalized to
micrometers before the ratio is computed, so a file with Z calibrated in millimeters and X/Y in
micrometers is still handled correctly.

1.5 Out-of-core streaming design

This was not part of the original plan — it emerged mid-project once native 5D reading
surfaced that the 3D texture feature families had no disk-streaming path at all (large 3D ROIs
either fit in RAM or the run failed). The full original design reasoning is embedded as
Appendix A (intensity/
histogram/surface streaming foundation) and Appendix B
(the GLSZM/GLDZM connected-component design).

• OutOfRamVoxelCloud (voxel_cloud_nontriv.{h,cpp}) — a disk-backed 3D voxel source, the
  3D analog of the existing 2D OutOfRamPixelCloud, with a block cache (65536 records/block) so
  reads aren't done one fseek/fread per voxel.
• Every 3D feature family streams differently, chosen per the family's actual data-dependency
  shape:
  • Intensity / histogram — stream the voxel cloud directly.
  • Surface — per-plane contours via read_slab(z); surface area reformulated as
    6N − 2·(adjacent voxel pairs) over a 2-plane occupancy-bitmap window; PCA axis lengths
    accumulated online; hull volume via det4.
  • GLCM — a ring of (offset+1) dense grey-binned Z-planes.
  • GLDM — 26-neighborhood matrix over a 3-plane sliding window.
  • NGLDM — per-interior-voxel 24-neighborhood.
  • NGTDM — a (2·radius+1)-plane streaming window.
  • GLRLM — all 13 direction shifts have dz ∈ {0,1}; dz=0 reuses the existing
    same-plane run-gathering code, dz=1 uses a "carry" array tracking in-progress runs across
    a 2-plane window.
  • GLSZM / GLDZM — both reduce to the same mechanism: streaming two-pass connected-
    component labeling via a growable union-find (features/streaming_ccl.h). The key
    insight, elaborated in Appendix B, is that only causally-earlier neighbors (13 for
    26-connectivity, 3 for 6-connectivity) need to be checked in scan order, which is what makes
    connected-component labeling streamable at all.
• Whole-volume OOC (no mask, entire image treated as one ROI) shares the same underlying
  primitives as the segmented (masked ROI) OOC path via populate_3d_voxel_cloud /
  run_3d_ooc_features.
• Guard: any 3D feature not yet ported to a streaming implementation raises a clear
  "not yet supported out-of-core" error (per-feature), rather than silently producing wrong
  values. The only unconditional loud-fail is for loader formats that cannot deliver data
  plane-by-plane at all (currently: NIfTI, which loads its whole 4D blob in one read).

1.6 Memory-limit / RAM budgeting design

--ramLimit (CLI) / ram_limit (Python) sets the threshold above which an ROI or whole volume
is treated as "oversized" (streamed OOC) rather than trivial (loaded whole). Default, if not
specified, is 50% of available system RAM (getAvailPhysMemory()); requesting a limit above
actual available RAM is rejected rather than silently overcommitting.

---

2. Features delivered

• Native 5D reading (X,Y,Z,C,T) for both OME-TIFF (strip and tiled, OME-XML
  DimensionOrder-driven) and OME-Zarr (both v2/"0.4" and v3/"0.5", both via the existing
  z5 library — no new dependency), replacing the old assumption of a single (c=0,t=0) plane.
• Per-channel, per-timepoint output: one row per (ROI, t, c). New c_index output column
  (alongside the pre-existing t_index) in every sink — CSV, in-memory buffer, and Arrow.
• Physical calibration, opt-in via --use-physical-spacing / use_physical_spacing=True:
  new phys_x, phys_y, phys_z (numeric) and phys_unit (string) output columns, populated
  from OME PhysicalSize*/coordinateTransformations metadata and canonicalized to
  micrometers across mismatched per-axis units. Manual override via --anisox/y/z still works
  and takes precedence, unchanged from before this project.
• True out-of-core 3D featurization, both workflows (segmented ROI and whole-volume, no
  mask), across all 7 texture feature families plus intensity and histogram — previously
  3D OOC did not exist at all; oversized 3D input either loaded fully into RAM or failed.
• Anisotropic resampling now works correctly for both whole-volume and segmented 3D
  workflows (a hang and 3 related correctness bugs in the whole-volume path were found and
  fixed — see §3.1).
• OME-Zarr v3 (0.5) support, including sharded stores, via the same z5-based reader used
  for v2 — the version is auto-detected (zarr.json vs .zattrs).

Not delivered / explicitly out of scope (see §7).

---

3. What was done

Work proceeded in two overlapping tracks. Numbers below are feature/bug counts; full narrative
detail is captured below and in the appendices.

3.1 Native 5D IO track (delivered across eight incremental phases)

1. OmeAxes descriptor + format detection (detect_input_format) for both containers.
2. Axis-aware loader interface (loadTileFromFile(..., channel, timeframe, level)) applied
   symmetrically across both loader stacks (10 concrete loader overrides).
3. OME-Zarr native 5D via axis roles instead of fixed shape indices.
4. OME-TIFF native 5D via OME-XML DimensionOrder (strip loaders, then tiled loaders).
5. Output schema: c_index column, per-channel default policy.
6. Physical calibration: --use-physical-spacing, phys_*/phys_unit columns, unit
   canonicalization.
7. OME-Zarr 0.5 (v3/sharding) support — discovered z5 3.0.0 reads v3 natively, so the originally
   planned tensorstore dependency was never added.
8. HCS / multi-file / labels-group masks — scoped, then explicitly declined (§7).

Along the way, this track found and fixed real bugs, independent of the OOC track: 3D whole-
volume featurization was completely non-functional at first (extract() called a stub that
always threw); a 3-layer prescan bug read out-of-bounds garbage and exposed a use-after-free in
a TIFF strip loader's buffer lifecycle; a mask read used the intensity channel index, dropping
ROIs for channel > 0; a Python binding call site was left on the old (pre-5D) argument count;
an OME-Zarr loadTile() always read exactly one Z-plane regardless of the chunk's actual
Z-depth (any store chunked with Z>1 silently returned zero past the first plane of each chunk);
COVERED_IMAGE_INTENSITY_RANGE was degenerate (always exactly 1.0) for any segmented 2D or
3D image, because the "whole slide" intensity baseline was computed through a voxel-visitor
that (by contract) only visits in-mask voxels — fixed for both 2D and 3D; physical-size units
were stored verbatim with no cross-unit conversion; the whole-volume anisotropic-resampling path
had a hang (a loop control variable was clobbered mid-body) plus three follow-on correctness
bugs surfaced while fixing it — a wrong row stride, a stale bounding box left over from before
resampling, and a stale voxel count that made 3MEAN come out exactly 4x too large — all four
fixed, and the fix was extended to the equivalent segmented-ROI code path once the same gap was
found there; and a size_t underflow in the ROI memory-footprint estimator (triggered whenever
the running count was zero) had been silently defeating ROI batching for essentially all
segmented 2D/2.5D/3D runs — not a crash or wrong output, but a performance/memory-management
regression, fixed narrowly without changing batching semantics further.

3.2 True 3D out-of-core track (not part of the original plan)

Delivered incrementally, one feature family (or family group) per increment:

1. Intensity/histogram streaming (the disk-backed voxel cloud and plane-by-plane streaming
   primitive, detailed in Appendix A).
2. Surface (area / PCA axis lengths / hull volume).
3. GLCM, GLDM, NGLDM, NGTDM.
4. GLRLM (plus two latent correctness bugs found and fixed: the GLCM/GLDM grey-level lookup
   table was built from mask-only voxels instead of including background voxels the way the
   in-core calculation does; GLRLM conflated two different IBSI-related flags between lookup-
   table construction and row lookup).
5. GLSZM + GLDZM (detailed in Appendix B; both passed their out-of-core-vs-in-RAM equivalence
   tests on the first try, attributed to scoping the streaming design in a written note before
   writing code).
6. A dedicated negative-path testing pass, prompted by asking "do we have negative tests for
   out-of-core features?" — the answer was no, so tests were added for degenerate/constant-
   intensity ROIs, unsupported-format loud-failure, and the unsupported-feature guard. One of
   these caught a real bug: EXCESS_KURTOSIS used an algebraic shortcut that diverges from the
   correct value specifically on zero-variance data.
7. Whole-volume out-of-core (previously only the segmented/masked-ROI path streamed;
   whole-volume input just failed loudly on an oversized volume).

26 real bugs were found and fixed across both tracks over the course of this project; none
were caught by CI before being found here — partway through, this codebase was resynced onto
16 commits of upstream churn, which required 2 additional semantic (non-textual) fixes to
reconcile an upstream PCA axis-length fix and a new upstream test's changed function signature
with the OOC/channel additions described here.

3.3 Final validation phase

The final planned validation phase — a 2D/3D/4D/5D scale/stress test on both Windows and a
memory-constrained Linux VM — had never been executed before the final session of this project.
It is now done; see §4 for results. This was the first time multi-channel or multi-timepoint
data (C>1 or T>1) had been exercised at real scale on either platform in this project.

---

4. Tests

4.1 Automated suites (Windows, current baseline)

Suite	Result
build-test (gtest, no z5)	799 / 799
build-test-z5 (gtest, with z5/OME-Zarr)	852 / 852
pytest tests/python/	96 passed, 1 skipped
Phase-0 oracle (io_oracle/, independent tifffile/ome-zarr-py cross-check)	PASS 21 / FAIL 0

New gtest coverage added by this project includes: whole-volume + segmented anisotropic
resampling correctness (TEST_3D_WHOLEVOLUME_ANISOTROPIC_REDUCE_MATCHES_ISOTROPIC,
TEST_3D_SEGMENTED_ANISOTROPIC_AUX_AREA_MATCHES_VOXELCLOUD), the n_rois==0 footprint-
estimate underflow guard, degenerate/constant-intensity OOC equivalence, unsupported-format
loud-failure, and the unsupported-feature OOC guard.

4.2 Manual memory/scale validation (Linux VM, 3.8 GB RAM / 2 cores)

Used throughout the OOC track to prove memory-boundedness under real memory pressure
(not just a --ramLimit argument), including one case where the in-RAM path was OOM-killed
by the kernel while the OOC path on the same fixture succeeded:

Fixture	OOC peak RSS	In-RAM peak RSS	Values
Single oversized ROI, 46.1M voxels (~2.21 GB footprint) — intensity/histogram	886 MB	2936 MB	byte-identical
Same ROI, surface features only	656 MB	OOM-killed (3530 MB)	3AREA exact match to analytical formula
Same ROI, all 4 GLCM/GLDM/NGLDM/NGTDM families	669 MB	3095 MB	byte-identical
GLRLM (separate ellipsoid-mask fixture, real background)	849 MB	3275 MB	byte-identical
Same ellipsoid fixture, all 7 texture families together	940 MB	3275 MB	byte-identical
Whole-volume, no mask (2.21 GB)	365 MB	2641 MB	byte-identical
4D-T / 4D-C / 5D (~0.75 GB per (c,t) frame)	200–232 MB	910–941 MB	byte-identical;3MEAN exact match to a hand-derived closed form for all frames

The 4D/5D result confirms peak RSS does not grow with frame count (3 vs. 4 frames gave
nearly identical peaks in both OOC and in-RAM modes) — no cross-frame memory accumulation in
either path.

4.3 Windows format-parity validation

A 5D (T=2, C=2, Z=60, Y=512, X=512) fixture generated as both OME-TIFF and OME-Zarr:
forced-OOC (ram_limit=200) vs. generous in-RAM (ram_limit=8000) — 4/4 rows, exact match to
the same hand-derived formula, np.allclose on every numeric column, for both formats.

---

5. How to use, with examples

5.1 CLI

Standard flags (--intDir, --segDir, --outDir, --filePattern, --outputType,
--features, --resultFname) are unchanged. New/relevant flags for this feature:

Flag	Meaning
--dim=3	3D (volumetric) mode
--ramLimit=<MB>	Threshold above which an ROI/volume streams out-of-core instead of loading whole. Default: 50% of available RAM.
--anisox, --anisoy, --anisoz	Manual per-axis resampling factors (unchanged from before this project).
--use-physical-spacing	Opt-in: derive resampling factors from the file's own OMEPhysicalSize* metadata instead of manual --aniso*.

Example — native 3D OME-TIFF, in-RAM, per-channel/per-timepoint output:

nyxus --dim=3 --intDir=/data/vol_intensity --segDir=/data/vol_masks --outDir=/output \
      --filePattern=".*\.ome\.tif" --features=*3D_ALL* --outputType=singlecsv \
      --resultFname=3d-features

Output rows now include one row per (ROI, t_index, c_index) for every channel/timepoint
found in the file's OME metadata — previously only (c=0,t=0) was ever read.

Example — same, but forced out-of-core (large volumes, constrained memory):

nyxus --dim=3 --intDir=/data/big_vol --segDir=/data/big_masks --outDir=/output \
      --filePattern=".*\.ome\.zarr" --features=*3D_ALL* --ramLimit=500 \
      --resultFname=3d-features-ooc

Any ROI (or the whole volume, in whole-volume mode) whose estimated footprint exceeds 500 MB
streams from disk instead of loading fully into RAM. Output is numerically identical to the
in-RAM path.

Example — physical-spacing-driven anisotropic resampling (opt-in):

nyxus --dim=3 --intDir=/data/calibrated --segDir=/data/calibrated_masks --outDir=/output \
      --filePattern=".*\.ome\.tif" --features=*3D_ALL* --use-physical-spacing \
      --resultFname=3d-features-calibrated

Resampling factors are derived from the file's own PhysicalSizeX/Y/Z (converted to a common
unit, then ratio-normalized) instead of requiring manual --anisox/y/z. Output gains phys_x,
phys_y, phys_z, phys_unit columns.

Example — manual anisotropy override (pre-existing, still supported, e.g. for NIfTI):

nyxus --anisox=1.5 --anisoy=2 --anisoz=2.5 --dim=3 --filePattern="*\.nii\.gz" \
      --features=*3D_ALL* --resultFname=3d-features --outputType=singlecsv \
      --intDir=/data/patient123/intensity --segDir=/data/patient123/masks \
      --outDir=/output/patient123

5.2 Python API (Nyxus3D)

import nyxus

# Directory-of-files, out-of-core, per-channel/per-timepoint output
nyx = nyxus.Nyxus3D(["*3D_ALL*"], ram_limit=500)
df = nyx.featurize_directory(
    intensity_dir="/data/big_vol",
    label_dir="/data/big_masks",
    file_pattern=".*\\.ome\\.zarr",
)

# Whole-volume (no mask), explicit file lists, single-ROI mode
nyx = nyxus.Nyxus3D(["*3D_ALL*"], ram_limit=500)
df = nyx.featurize_files(
    intensity_files=["/data/vol.ome.tif"],
    mask_files=[],          # ignored when single_roi=True
    single_roi=True,
)

# Physical-spacing-driven anisotropic resampling
nyx = nyxus.Nyxus3D(["*3D_ALL*"], use_physical_spacing=True)
df = nyx.featurize_directory("/data/calibrated", "/data/calibrated_masks", ".*\\.ome\\.tif")

# Manual anisotropy override (equivalent to --anisox/y/z)
nyx = nyxus.Nyxus3D(["*3D_ALL*"], anisotropy_x=1.5, anisotropy_y=2.0, anisotropy_z=2.5)

Relevant Nyxus3D constructor kwargs: ram_limit (int, MB; default -1 = no limit/whole-RAM),
use_physical_spacing (bool, default False), anisotropy_x/y/z (float, default 1.0).

Gotcha: ram_limit, use_physical_spacing, and the anisotropy factors are constructor-only
— set_environment_params(**params) does not accept them and always resets ram_limit to
"no limit" internally. To change these, construct a new Nyxus3D instance.

---

6. Changes to current usage

• Input format coverage expanded. Nyxus can now natively read 3D/4D/5D volumes directly
  from OME-TIFF and OME-Zarr (both v2/"0.4" and v3/"0.5"/sharded), not just from NIfTI. The
  in-repo README.md ("Nyxus also reads compressed and uncompressed NIFTI 3D files") and
  CLAUDE.md ("3D from NIFTI...") summaries predate this project and have not been updated
  — they describe the prior state, not current capability. Anyone updating those files
  should add native 3D OME-TIFF/OME-Zarr reading, --use-physical-spacing, out-of-core 3D
  featurization, and the c_index/phys_* output columns.
• New output columns, present in all sinks (CSV, buffer, Arrow), additive/non-breaking to
  existing consumers that read by column name: c_index (always present for 3D output now,
  even single-channel files get c_index=0), phys_x, phys_y, phys_z, phys_unit
  (populated when --use-physical-spacing/use_physical_spacing=True, otherwise default
  values).
• Multi-channel/multi-timepoint files now produce multiple rows per ROI (one per (t,c))
  where previously only the first plane's data was read — existing pipelines that assumed
  exactly one row per ROI for 3D input will see more rows for any file with C>1 or T>1. This is
  the "per-channel" default policy; there is no config to collapse channels back to one row —
  a select=/reduce= channel-policy mode was designed early on but never implemented, only
  perChannel shipped.
• 3D no longer requires the whole volume to fit in RAM. Previously, an oversized 3D
  ROI/volume either loaded fully into RAM (risking OOM) or, for whole-volume mode, failed
  outright. Now --ramLimit/ram_limit governs a genuine disk-streaming path with identical
  output to in-RAM, for all feature families. This is a strict capability addition — behavior
  for inputs that already fit under the RAM limit is unchanged.
• A pre-existing 2D-and-3D bug was fixed: COVERED_IMAGE_INTENSITY_RANGE previously always
  computed to exactly 1.0 for any segmented (non-whole-slide) image, in both 2D and 3D. It now
  reflects the true ratio of the ROI's intensity range to the whole slide's. Any downstream
  analysis that depended on the old (always-1.0) value will see different numbers after
  updating past this fix — this is a correctness fix, not a new feature, but is a behavior
  change on data that hasn't changed.
• A pre-existing performance bug was fixed: LR::get_ram_footprint_estimate(_3D)
  underflowed at n_rois==0, which defeated ROI batching for essentially all segmented 2D/2.5D/
  3D runs (every ROI was processed as a batch of one instead of being grouped). Fixed; segmented
  runs may now batch ROIs together where they previously didn't, which can change wall-clock
  time but not output values.
• No breaking changes to existing 2D/2.5D/NIfTI-3D usage. All changes were additive by
  design (§1.2); the full existing gtest/pytest suites pass unchanged.

---

7. Explicitly out of scope (declined, not deferred by oversight)

Decided 2026-07-22, after explicit scoping discussion.

• Multi-file OME-TIFF (<TiffData>/<UUID> cross-file plane references,
  .companion.ome). Detected (multiFileTiff flag) and rejected, never resolved — real
  support would need coordinating multiple open loader instances per logical dataset, a
  structural change rather than a bolt-on.
• HCS (plate/well) layouts. Never attempted — a fundamentally different discovery/
  addressing model, by far the largest of the three declined items.
• OME-Zarr labels/-group masks. Declined despite being comparatively small, to keep the
  format-coverage story consistent; Nyxus's existing two-file (separate intensity/mask)
  convention remains the only supported pairing.

None of the three are required by this project's actual goal (native 5D single-image reading +
out-of-core featurization, already delivered). Revisit only if a real dataset needs one of
these.

Other known gaps, carried over rather than newly introduced:

• GLDZM's dist2border remains 2D-only (a pre-existing limitation, replicated not fixed in
  the new streaming implementation).
• OME-TIFF explicit TiffData plane maps and SubIFD pyramids are not parsed — reads fall
  back to default rasterization.
• Whether real-world production data actually needs Zarr v3/sharding urgently was never
  resolved either way — it was implemented regardless since it came at no extra cost once z5's
  existing v3 support was discovered.

---

8. Open backlog (not blocking, not yet actioned)

• Batching-margin semantics: processTrivialRois/_25D/_3D's in-progress-batch-size
  callers pass Pending.size() (starts at 0 for every batch) to the footprint estimator, while
  other callers correctly pass the total slide ROI count (env.uniqueLabels.size()). The
  n==0 underflow this exposed is fixed (§3.1), but whether Pending.size() is the
  semantically right quantity for the neighbors-margin term is an open design question, not
  filed as a bug.

---

Appendix A: Design note — true out-of-core for 3D/4D/5D (foundation)

Original scoping note for the first out-of-core streaming increment: the disk-backed 3D voxel
source and streaming population, plus the first feature wins (intensity, histogram, surface).
Later families (GLCM/GLDM/NGLDM/NGTDM/GLRLM/GLSZM/GLDZM) reuse the same primitive and were
designed in follow-on notes (GLSZM/GLDZM's is Appendix B). All of this was implemented as
described; kept here as-written for the design rationale, not as a status report — see §3–4
above for what actually shipped.

A.1 The problem

A 3D segmented ROI whose footprint exceeds the RAM limit was routed to the same nontrivial-ROI
processing function used by the 2D path. That function scanned with a 2D-only tile loop —
no Z loop, and its pixel record had no Z coordinate. So for a 3D ROI the on-disk cloud collapsed
to a single 2D plane (Z was discarded and only plane 0 was ever read). Every downstream 3D
feature that consumed it was therefore wrong or fell back to requiring the whole volume in RAM:

• Intensity / histogram: read the (2D, plane-0) disk cloud for mean/energy/centroid/moments
  and the full in-RAM voxel vector for variance/MAD/median-AD/robust-mean. It only completed
  when the whole cube fit — i.e. it was not out-of-core at all.
• Texture (GLCM/GLRLM/GLDM/GLDZM/GLSZM/NGLDM): the out-of-core entry point was a stub
  falling through to the whole-cube trivial calculation. Surface: the out-of-core entry
  point was empty.
• Disk primitives were all 2D/Z-less: the existing out-of-core pixel cloud and its
  read/write matrix helpers had no Z dimension. No disk-backed 3D voxel source existed.

The full in-RAM voxel vector was the single object forcing the whole cube into RAM. Goal:
eliminate that dependency for intensity, histogram, and surface by streaming voxels off disk,
with bounded resident memory.

A.2 The foundation — a disk-backed 3D voxel cloud

A Z-slab-paged, disk-backed 3D voxel source — the direct 3D analog of the existing 2D
out-of-core pixel cloud — storing (x, y, z, intensity) records instead of (x, y, intensity).
Implemented as a new, separate primitive so the existing 2D path stays untouched.

On-disk layout — exploit Z-ascending population. The streaming population visits planes in
ascending Z, so records naturally arrive Z-sorted; they are simply appended, with an in-RAM
index remembering where each slab begins:

file: <one temp file per ROI per thread>
  record  = { int64 x, int64 y, int64 z, uint32 intensity }
  records written strictly in ascending z (then arbitrary y,x within a plane)

in-RAM index (built during population, negligible — one entry per Z-plane):
  slab_first[z] = index of first record with that z
  slab_count[z] = number of records in plane z  (0 for empty planes)

Depth is thousands at most, so the index is a few tens of KB; everything else stays on disk.

API shape:

class OutOfRamVoxelCloud
{
public:
    void   init (unsigned int roi_label, const std::string& name);
    void   close();
    void   clear();

    // Population (called in ascending-z order by the streaming scan)
    void   begin_slab (size_t z);
    void   add_voxel (const Pixel3& v);
    void   end_slab (size_t z);

    // Whole-cloud linear access (intensity / histogram — z-agnostic scalars)
    size_t size() const;
    Pixel3 get_at (size_t idx) const;
    struct iterator { ... };  iterator begin() const; iterator end() const;

    // Slab access (surface: 2-slab window; texture later: k-slab window)
    size_t depth() const;
    void   read_slab (size_t z, std::vector<Pixel3>& out) const;
    void   read_slab_window (size_t z0, size_t z1, std::vector<Pixel3>& out) const;
};

read_slab's resident set is one plane's ROI voxels, never the cube.

A.3 Streaming population

The existing nontrivial-ROI processing function is shared between 2D and 3D and could not
simply gain a Z-loop without risking the byte-identical 2D behavior already validated. So a
sibling 3D-specific population function was added instead, called from the 3D segmented
workflow's oversized branch in place of the 2D one — the 2D path is left exactly as-is.

The scan mirrors the loader's own volume-assembly Z-loop, but writes to disk per-plane instead
of into one giant in-RAM buffer:

for lz in 0 .. numberTileDepth-1:                    # z-slab of the loader
    cloud.begin_slab(gz)
    for (tr, tc) over the plane's tile grid:
        intFL->loadTileFromFile(ptrI, tr, tc, lz, channel, timeframe, level)
        segFL->loadTileFromFile(ptrL, tr, tc, lz, mask_channel, mask_timeframe, level)
        for each voxel in tile:
            if segL == r.label:
                cloud.add_voxel( Pixel3(globalX, globalY, gz, intI) )
    cloud.end_slab(gz)

Resident set per iteration is the plane's tile buffers (already allocated by the loader) plus
the current slab's ROI voxels being written — no whole-volume buffer is ever allocated,
which is the crux of the memory win. Channel/timeframe come from the workflow; mask
channel/timeframe use the same clamp the whole-volume loader already applies (mask channel
falls back to 0 if the mask has fewer channels than the intensity image).

An open item flagged for implementation time: how the loader's Z-slab argument maps for each
loader implementation at a given pyramid level, and confirming a loader that delivers its whole
volume in one read (rather than plane-by-plane) still produces identical voxels when the loop
degenerates to a single read.

A.4 Feature streaming rewrites (first wins)

Intensity + histogram — low effort. All intensity stats are per-voxel scalars
(Z-agnostic), so a linear pass over the voxel cloud suffices. Two-pass: pass 1 accumulates
mean/energy/min/max/moments; pass 2 (needs mean & median) accumulates variance/MAD/median-AD/
robust-mean. Both are sequential disk reads, resident memory ~ one record. Histogram gets a
cloud-based initializer mirroring the existing one. The math itself (formulas, robust-mean
window, hyperskew/flatness) is unchanged — only the data source changes, so the existing in-RAM
values are the oracle: forced-out-of-core must equal in-RAM exactly.

Surface (marching cubes) — moderate effort (previously unimplemented for out-of-core).
Marching cubes needs only two adjacent Z-slices at a time:

read_slab(0) -> prev
for z in 1 .. depth-1:
    read_slab(z) -> cur
    run marching-cubes / area+volume kernels on the (prev, cur) slice pair; accumulate
    prev = std::move(cur)

Peak resident memory is 2 planes of ROI voxels, reusing the in-core surface kernel on a 2-plane
dense window bounded by the plane bounding box, not the cube.

What this increment does not touch. Texture families stay loud-failing (or explicitly
skipped with a logged message) until streamed in later increments — never silently returning
zeros, matching the same principle already used for the whole-volume loud-fail case.

A.5 Validation approach

• Unit tests: a forced-out-of-core test per feature — populate the voxel cloud from a small
  fixture volume, run the streaming calculation, and assert equality with the trivial
  whole-cube calculation over the same volume (exact for integer stats, tight relative
  tolerance for float ones).
• Memory-constrained VM (small RAM, swap off): a segmented 3D volume whose ROI footprint
  exceeds the machine's RAM must featurize out-of-core with bounded peak memory (far below the
  full cube) and produce feature values identical to an in-RAM run on the same data (using a
  generous RAM limit for that baseline).
• Full regression suite (unit tests across build configurations, plus the Python test suite)
  green before any push.

A.6 Later families (reuse this same primitive, designed separately)

feature(s)	streaming approach	window
GLCM	co-occurrence over 13 directions; stream (z, z+dz) slab pairs; counts additive	2-slab
GLDM, NGLDM, NGTDM	3×3×3 neighborhood	3-slab sliding
GLRLM	runs along 13 directions; lines cross slabs, needs a halo/carry	3-slab + carry
GLSZM, GLDZM	connected grey-level zones (flood-fill)	whole-volume — hardest; needs an on-disk union-find / multi-pass design (Appendix B)

---

Appendix B: Design note — out-of-core GLSZM/GLDZM

Companion to Appendix A, which flagged GLSZM/GLDZM as "hardest — connected grey-level zones
(flood-fill), needs an on-disk union-find or multi-pass design." This note works out that
design concretely: it turns out to be one shared streaming mechanism — two-pass connected-
component labeling via a growable union-find — not two unrelated problems. As with Appendix A,
this was implemented as designed; both families passed their out-of-core-vs-in-RAM equivalence
tests on the first try, attributed to scoping this design before writing code.

B.1 Why these are different from the other texture families

Every feature streamed before this (GLCM/GLDM/NGLDM/NGTDM/GLRLM) only ever needed a bounded
sliding window (1–3 Z-planes, or a 2-plane carry for GLRLM), because each feature's basic
unit — a co-occurrence pair, a dependency count, a run — only involves voxels within a fixed,
small neighborhood.

GLSZM and GLDZM are different: their basic unit is a connected component ("zone") of
same-intensity voxels under 26- or 6-connectivity, and a zone's extent is a genuinely global
property — it can span the entire volume (e.g. one solid uniform-intensity region). A
component's membership only ever grows one voxel at a time via a local adjacency step, so no
single computation needs far-apart data at once — but knowing a zone's final size (or GLDZM's
minimum border distance) requires having fully explored every voxel that belongs to it,
wherever in the volume that turns out to be.

The classic answer for exactly this shape of problem is two-pass (or single-pass-with-
deferred-merges) connected-component labeling with union-find — a well-established streaming
technique, not something novel needed here. It fits a bounded window plus a small side table.

B.2 The shared mechanism: streaming connected-component labeling via union-find

Causal-neighbor trick. Process Z-planes in ascending order; within each plane, process
voxels in raster (row, column) order. For a voxel being visited, only check its causal
neighbors — those already fully labeled by scan order:

• Same plane: West, North, and (for 26-connectivity) North-West, North-East.
• Previous plane (z-1): already fully labeled by the time the current plane is reached, so
  every neighbor there counts as causal (not just the causal half).

The symmetric "future" neighbors (East, South, next-Z, …) are not missed — they get caught from
the other voxel's perspective when the scan reaches it (checking its own causal West/North/
previous-Z back to here). This is the standard raster-scan connected-component-labeling
neighbor set; for 26-connectivity it is 4 same-plane + 9 previous-plane = 13 causal neighbors
(exactly half of 26, mirroring the same canonical-half trick used for GLRLM's 13 direction
shifts, though here it decides what to check rather than how to walk a run). For GLDZM's
6-connectivity it reduces to 2 same-plane (West, North) + 1 previous-plane (straight down) = 3.

Union-find over labels, not voxels.

• A voxel with no matching causal neighbor gets a brand-new label (appended to a growable
  table).
• A voxel with one or more matching causal neighbors reuses (unions into) their label(s).
  If two different labels both match (the zone was discovered from two separate arms not yet
  known to be connected), union them.
• Each label's table entry carries the accumulated zone size (GLSZM) or running minimum
  distance (GLDZM), combined on union (sum for GLSZM, min for GLDZM) via standard
  union-by-size with path compression. No second pass over voxels is needed — after the whole
  volume streams once, each root label directly holds its zone's final size/min-distance.
• Final step: walk the (small) label table once, one (intensity, metric) pair per root, and
  feed that into the unchanged matrix-fill and feature-calculation functions — exactly the
  pattern used for the earlier texture families: the streaming path only replaces zone
  discovery, never the feature math.

Memory bound. The 2-plane window of dense grey-binned planes is O(area), same as every
other streamed feature. The union-find table is bounded by the number of distinct labels
created — typically far smaller than the voxel count (real zones are usually large), but worst
case (checkerboard-like data) approaches O(volume) in label count. This is the same class of
bound already accepted for GLRLM's per-row run-length histograms, and each label entry is a
handful of integers versus the far larger per-voxel footprint of holding the actual voxel cube.

B.3 GLSZM specifics

• Connectivity: full 26 (all 8 in-plane, 8 upper-Z-diagonal, 8 lower-Z-diagonal, plus
  straight up/down).
• Per-zone metric: size (voxel count).
• Background/skip convention: matlab-aware, differing from GLRLM's unconditional
  convention — must not reuse GLRLM's constant here.
• Grey-level lookup table: construction branches on the numeric IBSI-grey-binning flag
  (dense linear range vs. sparse unique-value set); row lookup separately branches on the
  boolean IBSI-compliance flag. This is the exact same two-flag split already found (and
  fixed) in GLRLM — the design must use the two flags for their respective purposes, not
  conflate them.
• Matrix fill is straightforward once zones are known: column = size − 1, row = grey-level
  index.

B.4 GLDZM specifics

• Connectivity: only 6 (no diagonals at all) — a simpler causal-neighbor set (3, not 13)
  than GLSZM.
• Per-zone metric: minimum, over all zone voxels, of each voxel's distance to the ROI
  border — and critically, this border-distance calculation is 2D-only: it scans
  left/right/up/down within the same Z-plane to the nearest background voxel or the plane
  edge, takes the minimum of the four, 1-based, clamped to at least 1. It does not consider
  the Z axis at all. This is an existing convention of the underlying calculation, not something
  to "fix" here — it must be replicated exactly. It also means the border-distance calculation
  for a voxel only needs the current plane's dense buffer (already resident in the 2-plane
  window) — no extra window cost versus GLSZM.
• Per-label running minimum (not sum), updated the same way as GLSZM's size but combined
  with min on union.
• Grey-level lookup table: uses the numeric IBSI-grey-binning flag for construction (already
  the correct convention) and a plain linear search (no separate boolean-flag branch) for row
  lookup — simpler than GLSZM, no two-flag trap here.
• Column = distance − 1 (0-based), row = grey-level index.

B.5 What must be replicated exactly (lessons from the earlier texture families)

Every texture feature streamed before this had at least one subtle convention that had to match
the whole-cube calculation exactly (background inclusion in the lookup table, the two-flag IBSI
split, a feature-specific background/shift quirk). Before implementing GLSZM/GLDZM, each of the
following needed confirming:

1. The exact causal-neighbor set (13 for GLSZM's 26-connectivity, 3 for GLDZM's
   6-connectivity) — getting this wrong causes zones to split or merge incorrectly, a much
   harder bug class to spot than a scalar value being off.
2. The background-skip convention for each family separately: GLSZM is matlab-aware; GLDZM
   skips voxels only in IBSI-binning mode when intensity is exactly zero, otherwise processes
   every value including background — different again from both GLRLM and GLSZM, and needed its
   own careful replication.
3. The two-flag split (grey-binning flag for lookup-table construction vs. IBSI-compliance flag
   for row lookup), confirmed present in GLSZM, absent in GLDZM.
4. Background inclusion in the lookup table when the ROI bounding box contains background
   voxels (the class of bug found and fixed for GLCM/GLDM earlier in this project) — checked
   for both.
5. Validated from the start with a non-cuboid-mask fixture (not as an afterthought) — it's what
   caught the background-lookup-table gap the first time this class of bug appeared.

B.6 Implementation plan (as executed)

1. Shared union-find helper: a growable parent[]/metric[] table with find (path
   compression) and a caller-supplied combine function (sum for GLSZM, min for GLDZM), so one
   helper serves both features.
2. GLSZM streaming calculation: 2-plane window, 13 causal neighbors, size-tracking
   union-find, final root walk → matrix fill → reuse existing feature-calculation functions.
3. GLDZM streaming calculation: 2-plane window, 3 causal neighbors, per-voxel border
   distance from the current plane, minimum-tracking union-find, final root walk → matrix fill
   → reuse existing feature-calculation functions.
4. Guard + tests: allow both families through the existing out-of-core feature guard;
   equivalence tests against the in-RAM calculation, including the partial/non-cuboid-mask
   fixture.
5. Memory-constrained VM validation: same methodology as the earlier texture families —
   in-RAM-vs-out-of-core value diff plus peak-memory comparison.

Both families completed all 7 texture families' out-of-core coverage — no 3D feature remained
loud-failed after this increment.

[vjaganat90]

Split or stage the merge
51 commits / 518 files / +10k lines is hard to review and hard to revert. I insist on stacked PRs or merge order: (A) OmeAxes + loaders + tests, (B) load_volume + c_index/(c,t) loop, (C) OOC voxel cloud + intensity/surface, (D) texture OOC, (E) physical spacing + docs. Even if history stays one branch, review in those layers.

Update user-facing docs in this PR (or a same-day follow-up)
README/CLAUDE still say “3D = NIfTI.” That is a merge blocker for usability: native OME-TIFF/Zarr 3D/5D, OOC, --use-physical-spacing, c_index/phys_*. Without docs, “feature complete” is only true for people who read the PR body.

CI must exercise the new paths
Local Windows 799/852 + one 3.8 GB Linux VM is strong evidence, but merge confidence needs CI for: USE_Z5 5D fixtures, OOC vs in-RAM equivalence tests, and ideally at least one forced-ramLimit 3D job. If OOC tests only run on the your machine, regressions will land on main unnoticed.

Even with full help of AI tools, I am having hard time right scoping and fully critiquing this PR.