Lightsheet pyramids: ingest from OME-Zarr or SmartSPIM, chunk materialization, cloud round trip, and napari viewing

src/ndi/+ndi/+fun/+doc/+lightsheet/README-lightsheet-zarr.md

@@ -0,0 +1,101 @@
+# lightsheet OME-Zarr documents (design + status)
+
+This package (`+ndi/+fun/+doc/+lightsheet/`) lets a session hold a
+lightsheet OME-Zarr (NGFF v0.4) multiscale pyramid as NDI documents,
+and hands one to napari for viewing through NDI-python.
+
+The pattern mirrors `+ndi/+fun/+doc/+gene/` (spatial transcriptomics
+Stereo-seq / GEF) so a reader of that package will recognise every
+piece here.
+
+## Document shape
+
+Two new document classes:
+
+- **`lightsheetZarrPyramid`** - one per REDUCTION variant (`mean`,
+  `max`, ...). Holds the frame (axes, level-0 shape, voxel size,
+  translation, channels, dtype). Depends on subject + optional element
+  + optional `fileReference` for the source store.
+- **`lightsheetZarrLevel`** - one per (pyramid, level). Owns the chunk
+  bytes as a `chunk.bin_#` file series and the per-level metadata
+  (shape, chunks, voxel size, translation, codec, ...). Depends on its
+  parent pyramid.
+
+NDI documents are immutable, so:
+- A new reduction is a new pyramid document + its own level ladder.
+- A deeper coarsening is a new level document; the parent's `n_levels`
+  is out of date until a new parent version is written for that count.
+
+## MATLAB API surface
+
+- `ndi.fun.doc.lightsheet.fromOMEZarr(session, zarrPath, ...)` -
+  probe an OME-Zarr store, write ONE pyramid document per multiscales
+  entry and ONE level document per level. Metadata-only; does not copy
+  chunk bytes. **Working.**
+- `ndi.fun.doc.lightsheet.makePyramid(session, pyramidEntry, ...)` -
+  the per-multiscales-entry worker `fromOMEZarr` calls. Reusable when
+  the caller already has an `ndr.format.omezarr.listPyramids` struct
+  in hand. **Working.**
+- `ndi.fun.doc.lightsheet.makeSourceFile(session, zarrPath)` -
+  fileReference for the store. **Working.**
+- `ndi.fun.doc.lightsheet.levelTable(session, pyramidDoc)` -
+  finest-first table of a pyramid's levels; `Properties.UserData`
+  carries the parent frame. **Working.**
+- `ndi.fun.doc.lightsheet.chooseLevel(levels, targetVoxelSize)` -
+  pick the coarsest level whose voxel size still meets the target.
+  **Working.**
+- `ndi.fun.doc.lightsheet.viewCommand(launcher, sessionPath, pyramidID, ...)`
+  - pure command builder. Same pattern as
+  `ndi.gui.app.GEFManager.viewCommand`; the launcher is
+  `/usr/local/bin/napariViewLightsheet`, a shell wrapper around the
+  `napariViewLightsheet` console script that ships in
+  `Waltham-Data-Science/NDI-python`. **Working.**
+
+## GUI
+
+`ndi.gui.app.LightsheetZarrManager(session)` mirrors `GEFManager`:
+lists pyramids in a session, and drives Add / View / Delete. `View`
+builds the napari command from `viewCommand` and shells out.
+
+## What still needs building
+
+1. **Chunk materialization** into `chunk.bin_#` on
+   `lightsheetZarrLevel`. Uses the same file-series upload path as
+   `ndi.cloud.sync.internal.uploadFilesForDatasetDocuments`, plus a
+   small `chunk_index` side-field so sparse pyramids can elide
+   fill-value chunks without breaking the 1-based `chunk.bin_#`
+   addressing. Land that side-field on the schema first, then a
+   `makePyramid` option (or a stand-alone `materializeChunks`
+   function) can write bytes.
+2. **NDI-python napari reader.** The MATLAB `viewCommand` targets
+   `/usr/local/bin/napariViewLightsheet`, a shell wrapper around
+   `ndi.gui.app.lightsheetZarr.cli:main` in
+   `Waltham-Data-Science/NDI-python` (installed as
+   `[project.scripts].napariViewLightsheet`). The CLI, viewer, and
+   metadata reader are landed; the per-chunk dask fetcher that
+   `levelArrays` calls is a scaffold today. See that repo's
+   `src/ndi/gui/app/lightsheetZarr/README-lightsheet-zarr.md`.
+3. **`readViewport` (MATLAB).** A MATLAB path for reading a rectangle
+   out of a pyramid at the appropriate level. Two backends are
+   plausible: read the source store via `ndr.format.omezarr.readArray`
+   (works today, but only when the store is on disk), or read the
+   materialized `chunk.bin_#` files via the cloud API (blocked on #1).
+   The scaffold from an earlier version of this PR was removed to
+   keep the package free of `notImplemented` bodies; add it back with
+   real code when #1 lands.
+4. **NDR helpers.** `ndr.format.omezarr.probe` is small; when it lands
+   in NDR-matlab / NDR-python, `fromOMEZarr` here can call it instead
+   of `listPyramids` for cheaper metadata reads. (Landing on the same
+   branch name in the NDR repos.)
+5. **Tests.** Once #1 lands, extend
+   `tests/+ndi/+unittest/+fun/+doc/+lightsheet/` with round-trip and
+   viewport tests analogous to `+gene/`'s.
+
+## Related
+
+- `ndi.gui.app.GEFManager` - the sibling manager for Stereo-seq GEF
+  spatial transcriptomics. Same launcher pattern (a shell wrapper
+  around an NDI-python console script). Read it first if this is your
+  first NDI file-format ingest.
+- `ndr.format.omezarr.*` - the metadata / read primitives this
+  package builds on.

src/ndi/+ndi/+fun/+doc/+lightsheet/chooseLevel.m

@@ -0,0 +1,64 @@
+function [levelIndex, levelRow] = chooseLevel(levels, targetVoxelSize)
+% NDI.FUN.DOC.LIGHTSHEET.CHOOSELEVEL - pick the coarsest level meeting a target
+%
+%   [LEVELINDEX, LEVELROW] = NDI.FUN.DOC.LIGHTSHEET.CHOOSELEVEL(LEVELS,
+%       TARGETVOXELSIZE)
+%
+%   Given LEVELS (a struct array or table with per-level `voxel_size`
+%   and `level` fields, finest-first, e.g. from LEVELTABLE), returns
+%   the coarsest level whose per-axis voxel size is still <=
+%   TARGETVOXELSIZE. If every level is coarser than the target, returns
+%   level 0. If every level is finer, returns the coarsest.
+%
+%   LEVELINDEX is the 0-based level index (matches
+%   lightsheetZarrLevel.level).
+%
+%   TARGETVOXELSIZE may be a scalar (applied to every spatial axis) or a
+%   vector aligned with the level's `voxel_size`.
+%
+%   See also: ndi.fun.doc.lightsheet.levelTable
+
+    arguments
+        levels
+        targetVoxelSize (1,:) double
+    end
+
+    if istable(levels)
+        vs = levels.voxel_size;
+        li = levels.level;
+        n = height(levels);
+    else
+        n = numel(levels);
+        vs = cell(n, 1);
+        li = zeros(n, 1);
+        for k = 1:n
+            vs{k} = reshape(double(levels(k).voxel_size), 1, []);
+            li(k) = levels(k).level;
+        end
+    end
+    if n == 0
+        error('NDI:lightsheet:chooseLevel:emptyLevels', ...
+            'levels is empty.');
+    end
+
+    bestK = 1;
+    for k = 1:n
+        v = vs{k};
+        if isscalar(targetVoxelSize)
+            ok = all(v <= targetVoxelSize);
+        else
+            m = min(numel(v), numel(targetVoxelSize));
+            ok = all(v(1:m) <= targetVoxelSize(1:m));
+        end
+        if ok
+            bestK = k;
+        end
+    end
+
+    levelIndex = li(bestK);
+    if istable(levels)
+        levelRow = levels(bestK, :);
+    else
+        levelRow = levels(bestK);
+    end
+end

src/ndi/+ndi/+fun/+doc/+lightsheet/fromOMEZarr.m

@@ -0,0 +1,185 @@
+function [pyramidDoc, levelDocs, info] = fromOMEZarr(session, zarrPath, options)
+% NDI.FUN.DOC.LIGHTSHEET.FROMOMEZARR - build lightsheetZarrPyramid docs from an OME-Zarr store
+%
+%   [PYRAMIDDOC, LEVELDOCS, INFO] = NDI.FUN.DOC.LIGHTSHEET.FROMOMEZARR(...
+%       SESSION, ZARRPATH)
+%
+%   Walks an OME-Zarr (NGFF v0.4) store at ZARRPATH, enumerates its
+%   multiscale pyramids with NDR.FORMAT.OMEZARR.LISTPYRAMIDS, and writes
+%   ONE lightsheetZarrPyramid document (the source volume) plus one
+%   lightsheetZarrLevel document per unique level to SESSION.
+%
+%   ONE DOCUMENT PER SOURCE VOLUME. The parent describes the source
+%   volume; reductions (mean, max, ...) live on the level children.
+%   A store with mean and max pyramids produces one parent, not two.
+%
+%   LEVEL 0 IS DEDUPED. When two NGFF multiscales entries reference
+%   the same underlying array (typically level 0: their
+%   `datasets[0].path` is the same string), one level document is
+%   written -- with reduction_function='none' -- and both readers use
+%   it. Reduced levels are per-reduction and get
+%   reduction_function='mean'/'max'/... .
+%
+%   READING IS METADATA-ONLY. No chunk bytes are read here; the
+%   source store's .zattrs and per-level .zarray are enough to fill
+%   the documents. Materializing chunk bytes into a level document's
+%   `chunk.bin_#` file series is follow-up work (see the package
+%   README).
+%
+%   Optional Name-Value Arguments:
+%   subjectID  - char/string, id of the subject depended on. If empty,
+%                the caller is refused: a lightsheetZarrPyramid without
+%                a subject is not a valid document.
+%   elementID  - char/string, id of the imaging element the volume was
+%                captured on. May be empty.
+%   pyramidNames - cellstr/string, restrict to these multiscales.name
+%                  values. Empty means all.
+%   reductionMap - containers.Map / struct mapping pyramid name (as
+%                  reported by listPyramids) or type to a reduction
+%                  label ('mean', 'max', ...). When a name is not in
+%                  the map, `pyramid_type` is inspected first, then a
+%                  fallback of 'mean' is used with a warning.
+%   pipelineVersion - char/string tag written into
+%                     lightsheetZarrPyramid.pipeline_version.
+%   sourceFileID - id of an existing fileReference document for
+%                  ZARRPATH. When empty, a new fileReference is created
+%                  via NDI.FUN.DOC.LIGHTSHEET.MAKESOURCEFILE.
+%   label      - human-readable label written into the parent pyramid.
+%
+%   INFO returns a struct with fields:
+%     zarrPath          - absolute path resolved
+%     pyramids          - the raw output of ndr.format.omezarr.listPyramids
+%     reductions        - cellstr, unique reductions present in the ladder
+%     sourceFileID      - id of the fileReference document used
+%     sharedLevel0      - logical, true if level 0 was shared across
+%                         reductions and deduped
+%
+%   Example:
+%     S = ndi.session.dir('mysession','/path/to/session');
+%     [P, L, info] = ndi.fun.doc.lightsheet.fromOMEZarr(S, ...
+%         '/cloudcache/subject42/volume.ome.zarr', ...
+%         'subjectID', subj.id());
+%
+%   See also: ndi.fun.doc.lightsheet.makePyramid,
+%             ndi.fun.doc.lightsheet.chooseLevel,
+%             ndi.fun.doc.lightsheet.viewCommand,
+%             ndr.format.omezarr.listPyramids,
+%             ndr.format.omezarr.probe
+
+    arguments
+        session (1,1)
+        zarrPath (1,:) char
+        options.subjectID char = ''
+        options.elementID char = ''
+        options.pyramidNames = {}
+        options.reductionMap = struct()
+        options.pipelineVersion char = ''
+        options.sourceFileID char = ''
+        options.label char = ''
+    end
+
+    if isempty(strtrim(options.subjectID))
+        error('NDI:lightsheet:fromOMEZarr:noSubject', ...
+            ['A lightsheetZarrPyramid depends_on a subject; pass ' ...
+             '''subjectID''.']);
+    end
+
+    if ~isfolder(zarrPath)
+        error('NDI:lightsheet:fromOMEZarr:noSuchStore', ...
+            'OME-Zarr store not found at %s.', zarrPath);
+    end
+
+    if ~ndr.format.omezarr.isOMEZarr(zarrPath)
+        error('NDI:lightsheet:fromOMEZarr:notOMEZarr', ...
+            '%s does not look like an OME-Zarr store (.zattrs missing multiscales).', ...
+            zarrPath);
+    end
+
+    pyramids = ndr.format.omezarr.listPyramids(zarrPath);
+    if isempty(pyramids)
+        error('NDI:lightsheet:fromOMEZarr:noPyramids', ...
+            '%s reports zero multiscales entries.', zarrPath);
+    end
+
+    keep = true(numel(pyramids), 1);
+    if ~isempty(options.pyramidNames)
+        wanted = cellstr(options.pyramidNames);
+        for i = 1:numel(pyramids)
+            keep(i) = any(strcmp(pyramids(i).name, wanted));
+        end
+        if ~any(keep)
+            error('NDI:lightsheet:fromOMEZarr:noNamedPyramid', ...
+                ['None of the requested pyramid names (%s) matched. ' ...
+                 'Available: %s.'], strjoin(wanted, ', '), ...
+                strjoin({pyramids.name}, ', '));
+        end
+    end
+    pyramids = pyramids(keep);
+
+    n = numel(pyramids);
+    perEntryReduction = cell(n, 1);
+    for i = 1:n
+        perEntryReduction{i} = pickReduction(pyramids(i), options.reductionMap);
+    end
+
+    sourceFileID = options.sourceFileID;
+    if isempty(sourceFileID)
+        srcDoc = ndi.fun.doc.lightsheet.makeSourceFile(session, zarrPath);
+        session.database_add(srcDoc);
+        sourceFileID = srcDoc.id();
+    end
+
+    [pyramidDoc, levelDocs, sharedLevel0] = ndi.fun.doc.lightsheet.makePyramid( ...
+        session, pyramids, perEntryReduction, ...
+        'subjectID', options.subjectID, ...
+        'elementID', options.elementID, ...
+        'sourceFileID', sourceFileID, ...
+        'pipelineVersion', options.pipelineVersion, ...
+        'label', options.label);
+
+    info = struct( ...
+        'zarrPath', char(zarrPath), ...
+        'pyramids', pyramids, ...
+        'reductions', {unique(perEntryReduction)}, ...
+        'sourceFileID', sourceFileID, ...
+        'sharedLevel0', sharedLevel0);
+end
+
+function r = pickReduction(pyramid, reductionMap)
+% pickReduction - choose 'mean' or 'max' for this pyramid entry
+%
+%   Preference: explicit map by name, then by type, then a heuristic:
+%   pyramid.type == 'max' -> 'max', anything else -> 'mean'. Warns when
+%   falling back because pyramid.type is unreliable.
+    r = '';
+    if ~isempty(reductionMap)
+        if isstruct(reductionMap)
+            if ~isempty(pyramid.name) && isfield(reductionMap, pyramid.name)
+                r = char(reductionMap.(pyramid.name));
+            elseif ~isempty(pyramid.type) && isfield(reductionMap, pyramid.type)
+                r = char(reductionMap.(pyramid.type));
+            end
+        elseif isa(reductionMap, 'containers.Map')
+            if ~isempty(pyramid.name) && isKey(reductionMap, pyramid.name)
+                r = char(reductionMap(pyramid.name));
+            elseif ~isempty(pyramid.type) && isKey(reductionMap, pyramid.type)
+                r = char(reductionMap(pyramid.type));
+            end
+        end
+    end
+    if isempty(r)
+        t = lower(strtrim(char(pyramid.type)));
+        n = lower(strtrim(char(pyramid.name)));
+        if strcmp(t, 'max') || contains(n, 'max')
+            r = 'max';
+        else
+            r = 'mean';
+            if ~isempty(t) && ~any(strcmp(t, {'mean','box','average'}))
+                warning('NDI:lightsheet:fromOMEZarr:reductionFallback', ...
+                    ['Pyramid %s has type ''%s''; falling back to ' ...
+                     '''mean''. Pass reductionMap to override.'], ...
+                    pyramid.name, pyramid.type);
+            end
+        end
+    end
+end