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

.gitignore

@@ -566,3 +566,8 @@ cython_debug/
 # setuptools-scm generated version file
 src/nyx/python/nyxus/_version.py
 .planning/
+
+# Feature tables written by the test suites into the tree
+NyxusFeatures.arrow
+NyxusFeatures.parquet
+NyxusFeatures.csv

CMakeLists.txt

@@ -127,6 +127,10 @@ set(SOURCE
 	src/nyx/raw_tiff.h
 	# sources
 	src/nyx/3rdparty/dsyevj3.cpp
+	# native OME metadata parsing + unified format detection
+	src/nyx/ome/format_detect.cpp
+	src/nyx/ome/ome_tiff_meta.cpp
+	src/nyx/ome/ome_zarr_meta.cpp
 	src/nyx/features/basic_morphology.cpp
 	src/nyx/features/caliper_feret.cpp
 	src/nyx/features/caliper_martin.cpp
@@ -154,6 +158,7 @@ set(SOURCE
 	src/nyx/features/hexagonality_polygonality.cpp
 	src/nyx/features/image_matrix.cpp
 	src/nyx/features/image_matrix_nontriv.cpp
+	src/nyx/features/voxel_cloud_nontriv.cpp
 	src/nyx/features/2d_geomoments.cpp
 	src/nyx/features/2d_geomoments_basic.cpp
 	src/nyx/features/2d_geomoments_basic_nt.cpp
@@ -223,6 +228,7 @@ set(SOURCE
 	src/nyx/phase2_25d.cpp
 	src/nyx/phase2_3d.cpp
 	src/nyx/phase3.cpp
+	src/nyx/phase3_3d.cpp
 	src/nyx/pixel_feed.cpp
 	src/nyx/raw_image_loader.cpp
 	src/nyx/reduce_trivial_rois.cpp
@@ -275,7 +281,25 @@ if(USE_Z5)
 		include_directories (${ZLIB_INCLUDE_DIR})
 		add_definitions(-DHAVE_ZLIB)	# telling NIFTI code of zlib availability, otherwise it'll support only uncompressed NIFTI
 	endif()
-
+
+	# zstd: the DEFAULT codec of Zarr v3 (OME-Zarr 0.5). z5's ZstdCompressor is gated by
+	# WITH_ZSTD; without it, zstd-compressed v3 chunks fail to decode. Optional: if zstd
+	# isn't found we still read uncompressed / blosc / zlib v3 stores. Defining WITH_ZSTD
+	# makes a z5 header include <zstd.h>, so the header must be found as well as the
+	# library -- a build environment can carry libzstd while its headers sit outside the
+	# dependency prefix, and defining the macro on the strength of the library alone turns
+	# that into a fatal include error.
+	find_library(ZSTD_LIBRARY NAMES zstd libzstd)
+	find_path(ZSTD_INCLUDE_DIR NAMES zstd.h)
+	if(ZSTD_LIBRARY AND ZSTD_INCLUDE_DIR)
+		add_definitions(-DWITH_ZSTD)
+		list(APPEND Nyxus_LIBRARIES ${ZSTD_LIBRARY})
+		include_directories(${ZSTD_INCLUDE_DIR})
+		message(STATUS "zstd found (${ZSTD_LIBRARY}) -> Zarr v3 zstd codec enabled")
+	else()
+		message(STATUS "zstd not fully available (lib='${ZSTD_LIBRARY}' include='${ZSTD_INCLUDE_DIR}') -> Zarr v3 zstd-compressed chunks will be unsupported")
+	endif()
+
 	if(nlohmann_json_FOUND)
 		include_directories(${nlohmann_json_INCLUDE_DIR})
 	endif()

README.md

@@ -15,7 +15,7 @@ Nyxus is a feature-rich, optimized, Python/C++ application capable of analyzing
 
 Nyxus can be used via Python or command line and is available in containerized form for reproducible execution. Nyxus computes over 450 combined intensity, texture, and morphological features at the ROI or whole image level with more in development. Key features that make Nyxus unique among other image feature extraction applications is its ability to operate at any scale, its highly validated algorithms, and its modular nature that makes the addition of new features straightforward.
 
-Currently, Nyxus can read 2D image data from OME-TIFF, OME-Zarr, and DICOM 2D Grayscale images. Nyxus also reads compressed and uncompressed NIFTI 3D files. Nyxus Python API supports featurizing in-memory 2D image data represented by NumPy arrays. 
+Currently, Nyxus can read 2D image data from OME-TIFF, OME-Zarr, and DICOM 2D Grayscale images. Nyxus reads 3D volumes from compressed and uncompressed NIFTI files, and natively from OME-TIFF and OME-Zarr (NGFF 0.4 and 0.5), including 5-dimensional (X, Y, Z, C, T) datasets whose axis order is taken from the file's own metadata rather than assumed. Nyxus Python API supports featurizing in-memory 2D image data represented by NumPy arrays. 
 
 The docs can be found at [Read the Docs](https://nyxus.readthedocs.io/en/latest/).
 
@@ -34,7 +34,7 @@ conda install nyxus -c conda-forge
 
 ## Usage
 
-The library provides class `Nyxus` for 2-dimensional TIFF, OME.TIFF, OME.ZARR, and DICOM slides and intensity-mask slide pairs, and class `Nyxus3D` for NIFTI volumes and intensity-volume volume pairs. Additionally to a single-file representation, a volume can be represented by its z-slices residing in separate 2-dimensional TIFF or ONE.TIFF slides (so called "layout A"). Slides and volumes can be featurized as all the files in a  directory filtered with a file pattern passed specified. Alternatively, explicit file name pairs and pair lists can be featurized. Alternatively, 2D and 3D NumPy arrays can be featurized. 
+The library provides class `Nyxus` for 2-dimensional TIFF, OME.TIFF, OME.ZARR, and DICOM slides and intensity-mask slide pairs, and class `Nyxus3D` for NIFTI, OME-TIFF, and OME-Zarr volumes and intensity-volume volume pairs. Additionally to a single-file representation, a volume can be represented by its z-slices residing in separate 2-dimensional TIFF or ONE.TIFF slides (so called "layout A"). Slides and volumes can be featurized as all the files in a  directory filtered with a file pattern passed specified. Alternatively, explicit file name pairs and pair lists can be featurized. Alternatively, 2D and 3D NumPy arrays can be featurized. 
 
 ### 2D usage
 
@@ -190,6 +190,25 @@ features2 = nyx.featurize_directory (idir, mdir, file_pattern=".*\.nii\.gz")
 ```
 <u>Note on CT datasets:</u> voxel intensities recorded in the Hounsfield units are automatically read by Nyxus as 0-based values by adding the minimum value of an original file (typically -1024).
 
+Volumes are also read natively from OME-TIFF and OME-Zarr, including 5-dimensional (X, Y, Z, C, T) datasets. The axis order comes from the file's own metadata, so non-canonical layouts are addressed correctly rather than assumed to be `TCZYX`:
+
+```python
+import nyxus
+nyx = nyxus.Nyxus3D (["*3D_ALL*"])
+
+features = nyx.featurize_directory (idir, mdir, file_pattern=".*\.ome\.zarr")
+```
+
+Each channel and timeframe of a multi-channel or time-series volume is featurized as a volume of its own, and the resulting rows carry `c_index` and `t_index` columns identifying which one they came from. A single-channel mask is reused across all intensity channels.
+
+<u>Note on physical voxel spacing:</u> OME-TIFF and OME-Zarr files may declare a physical size per axis. Passing `use_physical_spacing=True` makes Nyxus use each volume's declared spacing as its voxel anisotropy, ratio-normalized so that the smallest axis is 1:
+
+```python
+nyx = nyxus.Nyxus3D (["*3D_ALL*"], use_physical_spacing=True)
+```
+
+It is off by default, so the default cubic-voxel behavior is unchanged, and explicitly specified `anisotropy_x/y/z` values take precedence over the declared spacing. The spacing actually applied is reported in the `phys_x`, `phys_y`, `phys_z`, and `phys_unit` output columns, canonicalized to micrometers so that files declaring their sizes in different units stay numerically comparable.
+
 Featurizing explicitly specified volume files is straightforward, too:
 
 ```python
@@ -817,6 +836,12 @@ Applying a 1.5 x 2.0 x 2.5 anisotropy correction to all the volumes in a compres
 --anisox=1.5  --anisoy=2 --filePattern="*\.ome\.tiff" --features=*ALL*  --resultFname=features1 --outputType=singlecsv --intDir=/data/plate123/int --segDir=/data/plate123/seg --outDir=/output/plate123
 ```
 
+Instead of stating the anisotropy explicitly, it can be derived from the physical voxel spacing that OME-TIFF and OME-Zarr volumes declare in their own metadata:
+
+```
+--use-physical-spacing=true --dim=3 --filePattern=".*\.ome\.zarr" --features=*3D_ALL* --resultFname=3d-features --outputType=singlecsv --intDir=/data/patient123/intensity --segDir=/data/patient123/masks --outDir=/output/patient123
+```
+
 ## Nested ROIs 
 
 Hierarchical ROI analysis in a form of finding ROIs nested geometrically as nested AABBs and aggregating features of child ROIs within corresponding parent is available as an optional extra step after the feature extraction of the whole image set is finished. To enable this step, all the command line options '--hsig', '--hpar', '--hchi', and '--hag' need to have non-blank valid values. 

src/nyx/abs_tile_loader.h

@@ -32,12 +32,16 @@ class AbstractTileLoader
   /// @param tile Tile to load
   /// @param indexRowGlobalTile Tile's row index in the file to load
   /// @param indexColGlobalTile Tile's col index in the file to load
-  /// @param indexLayerGlobalTile Tile's layer index in the file to load
+  /// @param indexLayerGlobalTile Tile's layer (Z) index in the file to load
+  /// @param indexChannel Channel (C) index to load; 0 for single-channel data
+  /// @param indexTimeframe Time (T) index to load; 0 for non-time-series data
   /// @param level Tile's pyramidal level in the file to load
   virtual void loadTileFromFile(std::shared_ptr<std::vector<DataType>> tile,
                                 size_t indexRowGlobalTile,
                                 size_t indexColGlobalTile,
                                 size_t indexLayerGlobalTile,
+                                size_t indexChannel,
+                                size_t indexTimeframe,
                                 size_t level) = 0;
 
   /// \brief Getter to full Height
@@ -70,6 +74,14 @@ class AbstractTileLoader
     return 1;
   }
 
+  /// \brief Physical voxel spacing along X/Y/Z (default 1.0 == uncalibrated).
+  /// OME loaders return the parsed PhysicalSize*; everything else stays 1.0.
+  [[nodiscard]] virtual double physicalSizeX() const { return 1.0; }
+  [[nodiscard]] virtual double physicalSizeY() const { return 1.0; }
+  [[nodiscard]] virtual double physicalSizeZ() const { return 1.0; }
+  /// \brief Physical-size unit string (e.g. "micrometer"); empty if unknown.
+  [[nodiscard]] virtual std::string physicalSizeUnit() const { return std::string(); }
+
   /// \brief Getter to tile Width
   /// @param level tile's level considered
   /// \return Tile Width

src/nyx/cli_option_constants.h

@@ -69,6 +69,7 @@
 #define clo_AGGREGATE "--aggr"				// -> raw_aggregate
 #define clo_ANNOTATE "--annot"				// -> raw_annotate
 #define clo_ANNOT_SEP "--annotsep"		// -> raw_anno_separator
+#define clo_USE_PHYSICAL_SPACING "--use-physical-spacing"	// -> rawUsePhysicalSpacing (true/false); opt-in OME PhysicalSize* voxel spacing
 
 // Valid values of 'OUTPUTTYPE'
 #define clo_OT_SEPCSV "separatecsv"

src/nyx/constants.h

@@ -1,4 +1,5 @@
 #pragma once
 
 constexpr double DEFAULT_T_INDEX = 0.0;
+constexpr double DEFAULT_C_INDEX = 0.0;		// FIX (IO): default channel index for single-channel / non-OME input
 constexpr int DEFAULT_NUM_HISTO_BINS = 24;

src/nyx/environment.cpp

@@ -473,6 +473,7 @@ bool Environment::parse_cmdline(int argc, char** argv)
 			|| find_string_argument(i, clo_FPIMAGE_MIN, fpimageOptions.raw_min_intensity)
 			|| find_string_argument(i, clo_FPIMAGE_MAX, fpimageOptions.raw_max_intensity)
 			|| find_string_argument(i, clo_PRESERVE_HU, fpimageOptions.raw_preserve_hu)		// CT/HU mode
+			|| find_string_argument(i, clo_USE_PHYSICAL_SPACING, rawUsePhysicalSpacing)		// FIX (IO): opt-in physical voxel spacing
 			|| find_string_argument(i, clo_RESULTFNAME, nyxus_result_fname)
 			|| find_string_argument(i, clo_CLI_DIM, raw_dim)
 
@@ -799,6 +800,14 @@ bool Environment::parse_cmdline(int argc, char** argv)
 		}
 	}
 
+	//==== Parse opt-in physical-spacing flag (FIX (IO); default false)
+	if (!rawUsePhysicalSpacing.empty())
+	{
+		std::string v = rawUsePhysicalSpacing;
+		std::transform(v.begin(), v.end(), v.begin(), ::tolower);
+		use_physical_spacing_ = (v == "true" || v == "t" || v == "1" || v == "on" || v == "yes");
+	}
+
 	//==== Parse exclusive-inclusive timing
 	#ifdef CHECKTIMING
 	if (!rawExclusiveTiming.empty())

src/nyx/environment.h

@@ -2,6 +2,7 @@
 
 #include <map>
 #include <string>
+#include <unordered_set>		// FIX: csv_paths_written registry (see below)
 #include <vector>
 
 #include "arrow_output_stream.h"
@@ -80,13 +81,33 @@ class Environment: public BasicEnvironment
 	std::string rawOutpType; //= ""; // Valid values: "separatecsv", "singlecsv", "arrow", "parquet"
 	bool separateCsv; //= true;
 
+	// FIX: CSV output paths already written during this run. A slide now yields one row per
+	// (channel, timeframe), and the CSV sinks are called once per plane -- but in separatecsv
+	// mode they opened the same path with "w" every time, so each plane truncated the previous
+	// and only the LAST one survived. The first write to a path truncates and renders the
+	// header; later writes append. Held per-Environment rather than in a function-static so a
+	// long-lived Python session starts each run clean instead of appending to a stale file.
+	std::unordered_set<std::string> csv_paths_written;
+
+	// True the first time `path` is written in this run (and records it).
+	bool csv_claim_first_write (const std::string& path) { return csv_paths_written.insert(path).second; }
+	void reset_csv_output_state() { csv_paths_written.clear(); }
+
 	Nyxus::SaveOption saveOption;
 
 	// x- and y- resolution in pixels per centimeter
 	std::string rawXYRes; //= "";
 	float xyRes; //= 0.0,
 	float pixelSizeUm; //= 0.0;
 
+	// FIX (IO): opt-in physical-voxel-spacing calibration. When set, the 3D pipeline uses
+	// each slide's OME PhysicalSize* (ratio-normalized) as voxel spacing (anisotropy).
+	// Off by default -> today's cube-voxel behavior is unchanged. Spacing is still emitted
+	// as phys_x/y/z + phys_unit output columns regardless of this flag.
+	std::string rawUsePhysicalSpacing; //= "";
+	bool use_physical_spacing_ = false;
+	bool use_physical_spacing() const { return use_physical_spacing_; }
+
 	int get_pixel_distance();
 	void set_pixel_distance(int pixelDistance);
 	size_t get_ram_limit();

src/nyx/feature_method.cpp

@@ -43,7 +43,10 @@ void FeatureMethod::add_dependencies(const std::initializer_list<Nyxus::Feature3
 		dependencies.push_back((int)f);
 }
 
-void FeatureMethod::osized_scan_whole_image (LR& r, const Fsettings& s, ImageLoader& ldr)
+// Default out-of-core dispatch: ignore the Dataset and run the Dataset-less osized_calculate,
+// so features that don't need slide props are unaffected; only features that do (intensity,
+// intensity-histogram) override this and consume the Dataset.
+void FeatureMethod::osized_scan_whole_image (LR& r, const Fsettings& s, const Dataset& ds, ImageLoader& ldr)
 {
 	this->osized_calculate (r, s, ldr);
 	this->save_value (r.fvals);

src/nyx/feature_method.h

@@ -5,6 +5,8 @@
 #include "image_loader.h"
 #include "features/image_matrix_nontriv.h"
 
+class Dataset;	// forward decl for the Dataset-aware osized_scan_whole_image overload
+
 /// @brief Abstract class encapsulating basic feature functionality e.g. dependency on other features, dependency on helper objects, state (calculated or pending), etc.
 class FeatureMethod
 {
@@ -25,7 +27,10 @@ class FeatureMethod
 	virtual void parallel_process(std::vector<int>& roi_labels, std::unordered_map <int, LR>& roiData, int n_threads)  {}
 
 	//=== Oversized ROI
-	virtual void osized_scan_whole_image (LR& roi, const Fsettings& s, ImageLoader& ldr);
+	// Out-of-core dispatch. Features that need slide props (intensity, intensity-histogram)
+	// override this to reach their Dataset-taking osized_calculate; the default ignores the
+	// Dataset and calls the Dataset-less osized_calculate, so every other feature is unaffected.
+	virtual void osized_scan_whole_image (LR& roi, const Fsettings& s, const Dataset& ds, ImageLoader& ldr);
 
 	virtual void osized_add_online_pixel (size_t x, size_t y, uint32_t intensity) = 0;	// Called each time the ROI pixel is being scanned in the raster order
 	virtual void osized_reduce() final {};	// Get rid of this method in all derived

src/nyx/features/3d_glcm.cpp

@@ -179,7 +179,14 @@ void D3_GLCM_feature::save_value(std::vector<std::vector<double>>& fvals)
 void D3_GLCM_feature::extract_texture_features_at_angle(int dx, int dy, int dz, const SimpleCube<PixIntens>& binned_greys, PixIntens min_val, PixIntens max_val, int n_greys, bool ibsi, double soft_nan)
 {
 	calculateCoocMatAtAngle (P_matrix, dx, dy, dz, binned_greys, min_val, max_val, n_greys, ibsi);
+	finalize_angle (soft_nan);
+}
 
+// Computes this angle's GLCM feature values from the already-built P_matrix + sum_p + I. Split out
+// of extract_texture_features_at_angle so the out-of-core path (which builds P_matrix by streaming a
+// 2-plane window instead of scanning the in-RAM cube) reuses the exact same feature math.
+void D3_GLCM_feature::finalize_angle (double soft_nan)
+{
 	// Blank cooc-matrix? -- no point to use it, assign each feature value '0' and return.
 	if (sum_p == 0)
 	{

src/nyx/features/3d_glcm.h

@@ -145,6 +145,8 @@ class D3_GLCM_feature : public FeatureMethod, public TextureFeature
 		bool normalize);
 
 	void extract_texture_features_at_angle (int dx, int dy, int dz, const SimpleCube<PixIntens> & grays, PixIntens min_val, PixIntens max_val, int n_greys, bool ibsi, double soft_nan);
+	// Compute this angle's features from the already-built P_matrix (shared by in-core + out-of-core)
+	void finalize_angle (double soft_nan);
 
 	void calculateCoocMatAtAngle(
 		// out