From 25729dd3d3883b317d8684a72a1caccc091dbe3a Mon Sep 17 00:00:00 2001 From: Brandon Date: Tue, 3 Mar 2026 10:44:34 -0600 Subject: [PATCH] Add feature maps (fmaps) mode for sliding-kernel feature extraction --- CMakeLists.txt | 2 + README.md | 75 ++- ci-utils/envs/conda_cpp.txt | 4 +- docs/source/References/Classes.rst | 13 +- docs/source/cmdline_and_examples.rst | 44 +- src/nyx/cli_option_constants.h | 5 + src/nyx/environment.cpp | 43 +- src/nyx/environment.h | 18 +- src/nyx/features/3d_glcm.cpp | 9 +- src/nyx/features/3d_ngtdm.cpp | 2 +- src/nyx/features/glcm.cpp | 9 +- src/nyx/features/gldm.cpp | 2 +- src/nyx/features/gldzm.cpp | 2 +- src/nyx/features/glrlm.cpp | 2 +- src/nyx/features/glszm.cpp | 2 +- src/nyx/features/ngtdm.cpp | 2 +- src/nyx/globals.h | 113 ++++ src/nyx/helpers/helpers.h | 9 +- src/nyx/main_nyxus.cpp | 50 +- src/nyx/output_2_buffer.cpp | 603 +++++------------ src/nyx/output_2_csv.cpp | 954 +++++++++++---------------- src/nyx/parallel.h | 5 + src/nyx/python/new_bindings_py.cpp | 258 ++++++-- src/nyx/python/nyxus/__init__.py | 1 + src/nyx/python/nyxus/fmap_io.py | 174 +++++ src/nyx/python/nyxus/nyxus.py | 222 +++++-- src/nyx/results_cache.h | 25 +- src/nyx/workflow_2d_fmaps.cpp | 397 +++++++++++ src/nyx/workflow_2d_segmented.cpp | 3 +- src/nyx/workflow_2d_whole.cpp | 2 + src/nyx/workflow_3d_fmaps.cpp | 378 +++++++++++ src/nyx/workflow_pythonapi.cpp | 208 +++++- tests/CMakeLists.txt | 2 +- tests/test_all.cc | 90 +++ tests/test_binning_origin.h | 145 ++++ tests/test_fmaps.h | 214 ++++++ tests/test_fmaps_3d.h | 216 ++++++ 37 files changed, 3127 insertions(+), 1176 deletions(-) create mode 100644 src/nyx/python/nyxus/fmap_io.py create mode 100644 src/nyx/workflow_2d_fmaps.cpp create mode 100644 src/nyx/workflow_3d_fmaps.cpp create mode 100644 tests/test_binning_origin.h create mode 100644 tests/test_fmaps.h create mode 100644 tests/test_fmaps_3d.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 4a10fa2c..7bbb76b7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -237,7 +237,9 @@ set(SOURCE src/nyx/slideprops.cpp src/nyx/strpat.cpp src/nyx/workflow_2d_segmented.cpp + src/nyx/workflow_2d_fmaps.cpp src/nyx/workflow_2d_whole.cpp + src/nyx/workflow_3d_fmaps.cpp src/nyx/workflow_3d_segmented.cpp src/nyx/workflow_3d_whole.cpp src/nyx/workflow_pythonapi.cpp diff --git a/README.md b/README.md index 44c13091..9e96d00a 100644 --- a/README.md +++ b/README.md @@ -234,6 +234,54 @@ f = nyx.featurize_directory (intensity_dir=dir, label_dir=dir) ``` +### Feature maps (sliding kernel) mode + +Feature maps mode computes features at every position of a sliding kernel across each ROI, producing spatial feature maps instead of a single feature vector per ROI. This is useful for generating spatially resolved feature representations for downstream analysis or machine learning. + +#### 2D feature maps + +```python +from nyxus import Nyxus, save_fmaps_to_tiff + +nyx = Nyxus(["*ALL_INTENSITY*"], fmaps=True, fmaps_radius=2) # 5x5 kernel +results = nyx.featurize_directory("/path/to/intensities", "/path/to/labels") + +# results is a list of dicts, one per parent ROI: +# [ +# { +# "parent_roi_label": 1, +# "intensity_image": "img1.tif", +# "mask_image": "seg1.tif", +# "origin_x": 10, "origin_y": 20, +# "features": { +# "MEAN": numpy.array(shape=(map_h, map_w)), +# "STDDEV": numpy.array(shape=(map_h, map_w)), +# ... +# } +# }, +# ... +# ] + +# Save feature maps as TIFF stacks (requires tifffile) +save_fmaps_to_tiff(results, "output/tiff/") +``` + +#### 3D feature maps + +```python +from nyxus import Nyxus3D, save_fmaps_to_nifti + +nyx = Nyxus3D(["*3D_ALL_INTENSITY*"], fmaps=True, fmaps_radius=1) # 3x3x3 kernel +results = nyx.featurize_directory("/path/to/volumes", "/path/to/masks") + +# Each feature map is a 3D numpy array shaped (map_d, map_h, map_w) + +# Save as NIfTI volumes (requires nibabel) +save_fmaps_to_nifti(results, "output/nifti/", voxel_size=(0.5, 0.5, 1.0)) +``` + +Note: Feature maps mode returns numpy arrays rather than DataFrames, since the output is inherently image-shaped. Arrow/Parquet output is not supported in this mode. + ## Further steps For more information on all of the available options and features, check out [the documentation](#). @@ -271,19 +319,20 @@ print(nyx.get_params()) will print the dictionary ```bash -{'coarse_gray_depth': 256, -'features': ['*ALL*'], -'gabor_f0': 0.1, -'gabor_freqs': [1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0], -'gabor_gamma': 0.1, -'gabor_kersize': 16, -'gabor_sig2lam': 0.8, -'gabor_theta': 45.0, -'gabor_thold': 0.025, -'ibsi': 0, -'n_loader_threads': 1, -'n_feature_calc_threads': 4, -'neighbor_distance': 5, +{'binning_origin': 'zero', +'coarse_gray_depth': 256, +'features': ['*ALL*'], +'gabor_f0': 0.1, +'gabor_freqs': [1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0], +'gabor_gamma': 0.1, +'gabor_kersize': 16, +'gabor_sig2lam': 0.8, +'gabor_theta': 45.0, +'gabor_thold': 0.025, +'ibsi': 0, +'n_loader_threads': 1, +'n_feature_calc_threads': 4, +'neighbor_distance': 5, 'pixels_per_micron': 1.0} ``` diff --git a/ci-utils/envs/conda_cpp.txt b/ci-utils/envs/conda_cpp.txt index 8a963274..2cb5b864 100644 --- a/ci-utils/envs/conda_cpp.txt +++ b/ci-utils/envs/conda_cpp.txt @@ -9,5 +9,5 @@ xsimd >=13,<14 cmake dcmtk >=3.6.9 fmjpeg2koj >=1.0.3 -libarrow -libparquet +libarrow <=23.0.0 +libparquet <=23.0.0 diff --git a/docs/source/References/Classes.rst b/docs/source/References/Classes.rst index 113ae302..1f334b0e 100644 --- a/docs/source/References/Classes.rst +++ b/docs/source/References/Classes.rst @@ -1,7 +1,16 @@ Nyxus Classes ================ .. autosummary:: - :toctree: stubs + :toctree: stubs nyxus.Nyxus - nyxus.Nested \ No newline at end of file + nyxus.Nyxus3D + nyxus.Nested + +Feature Map I/O +================ +.. autosummary:: + :toctree: stubs + + nyxus.save_fmaps_to_tiff + nyxus.save_fmaps_to_nifti \ No newline at end of file diff --git a/docs/source/cmdline_and_examples.rst b/docs/source/cmdline_and_examples.rst index 6783729d..a0d5a770 100644 --- a/docs/source/cmdline_and_examples.rst +++ b/docs/source/cmdline_and_examples.rst @@ -53,6 +53,11 @@ should adhere to columns "WIPP I/O role" and "WIPP type". - integer - input - integer + * - --binningOrigin + - (Optional) Origin of the intensity binning range for texture features. 'zero' bins span [0, max] (default). 'min' bins span [min, max], adapting to the actual data range (PyRadiomics-compatible). Default: '--binningOrigin=zero' + - string constant + - input + - enum * - --gaborfreqs - (Optional) Feature GABOR: custom denominators of :math:`\pi` as frequencies of Gabor filter's harmonic factor. Default: '--gaborfreqs=1,2,4,8,16,32,64' - list of integer constants @@ -135,9 +140,19 @@ should adhere to columns "WIPP I/O role" and "WIPP type". - path * - --arrowOutputType - (Optional) Type of Arrow file to write the feature results to. Current options are 'arrow' for Arrow IPC or 'parquet' for Parquet - - string + - string - output - enum + * - --fmaps + - (Optional) Enable feature maps mode. When enabled, a sliding kernel is moved across each ROI and features are computed at every position, producing spatial feature maps instead of a single feature vector per ROI. Acceptable values: true, false. Default: '--fmaps=false'. Not compatible with Arrow/Parquet output. + - string constant + - input + - enum + * - --fmapsRadius + - (Optional) Radius of the sliding kernel in feature maps mode. The kernel size is (2 * radius + 1). For example, '--fmapsRadius=2' produces a 5x5 kernel (2D) or 5x5x5 kernel (3D). Default: '--fmapsRadius=2' + - integer + - input + - integer Examples ======== @@ -378,19 +393,20 @@ will print the dictionary .. code-block:: bash - {'coarse_gray_depth': 256, - 'features': ['*ALL*'], - 'gabor_f0': 0.1, - 'gabor_freqs': [1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0], - 'gabor_gamma': 0.1, - 'gabor_kersize': 16, - 'gabor_sig2lam': 0.8, - 'gabor_theta': 45.0, - 'gabor_thold': 0.025, - 'ibsi': 0, - 'n_loader_threads': 1, - 'n_feature_calc_threads': 4, - 'neighbor_distance': 5, + {'binning_origin': 'zero', + 'coarse_gray_depth': 256, + 'features': ['*ALL*'], + 'gabor_f0': 0.1, + 'gabor_freqs': [1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0], + 'gabor_gamma': 0.1, + 'gabor_kersize': 16, + 'gabor_sig2lam': 0.8, + 'gabor_theta': 45.0, + 'gabor_thold': 0.025, + 'ibsi': 0, + 'n_loader_threads': 1, + 'n_feature_calc_threads': 4, + 'neighbor_distance': 5, 'pixels_per_micron': 1.0} There is also the option to pass arguments to this function to only receive a subset of parameter values. The arguments should be diff --git a/src/nyx/cli_option_constants.h b/src/nyx/cli_option_constants.h index a782d715..6b6272f6 100644 --- a/src/nyx/cli_option_constants.h +++ b/src/nyx/cli_option_constants.h @@ -16,6 +16,7 @@ #define clo_XYRESOLUTION "--pixelsPerCentimeter" // pixels per centimeter #define clo_PXLDIST "--pixelDistance" // used in neighbor features #define clo_COARSEGRAYDEPTH "--coarseGrayDepth" // Environment :: raw_coarse_grayscale_depth +#define clo_BINNINGORIGIN "--binningOrigin" // Environment :: "zero" (default) or "min" (PyRadiomics-style) #define clo_RAMLIMIT "--ramLimit" // Optional. Limit for treating ROIs as non-trivial and for setting the batch size of trivial ROIs. Default - amount of available system RAM #define clo_TEMPDIR "--tempDir" // Optional. Used in processing non-trivial features. Default - system temp directory #define clo_IBSICOMPLIANCE "--ibsi" // skip binning for grey level and grey tone features @@ -61,6 +62,10 @@ #define clo_ANISO_Y "--anisoy" #define clo_ANISO_Z "--anisoz" +// Feature maps +#define clo_FMAPS "--fmaps" // Enable feature maps mode. "true" or "false" +#define clo_FMAPS_RADIUS "--fmapsRadius" // Kernel radius for feature maps. Integer >= 1. Example: "2" (produces 5x5 kernel) + // Result options #define clo_NOVAL "--noval" // -> raw_noval #define clo_TINYVAL "--tinyval" // -> raw_tiny diff --git a/src/nyx/environment.cpp b/src/nyx/environment.cpp index 1e0712d2..6a543210 100644 --- a/src/nyx/environment.cpp +++ b/src/nyx/environment.cpp @@ -191,6 +191,8 @@ void Environment::show_cmdline_help() << "\t\t\tDefault: 5 \n" << "\t\t" << OPT << clo_COARSEGRAYDEPTH << "= \n" << "\t\t\tDefault: 64 \n" + << "\t\t" << OPT << clo_BINNINGORIGIN << "= Origin of intensity binning range \n" + << "\t\t\t'zero' bins from [0, max] (default), 'min' bins from [min, max] (PyRadiomics-compatible) \n" << "\t\t" << OPT << clo_GLCMANGLES << "= \n" << "\t\t\tDefault: 0,45,90,135 \n" << "\t\t" << OPT << clo_VERBOSITY << "= \n" @@ -441,6 +443,7 @@ bool Environment::parse_cmdline(int argc, char** argv) find_string_argument(i, clo_GLCMOFFSET, glcmOptions.rawOffs) || find_string_argument(i, clo_PXLDIST, pixel_distance) || find_string_argument(i, clo_COARSEGRAYDEPTH, raw_coarse_grayscale_depth) || + find_string_argument(i, clo_BINNINGORIGIN, raw_binning_origin) || find_string_argument(i, clo_VERBOSITY, rawVerbosity) || find_string_argument(i, clo_IBSICOMPLIANCE, raw_ibsi_compliance) || find_string_argument(i, clo_RAMLIMIT, rawRamLimit) || @@ -474,6 +477,9 @@ bool Environment::parse_cmdline(int argc, char** argv) || find_string_argument(i, clo_RESULTFNAME, nyxus_result_fname) || find_string_argument(i, clo_CLI_DIM, raw_dim) + || find_string_argument(i, clo_FMAPS, raw_fmaps) + || find_string_argument(i, clo_FMAPS_RADIUS, raw_fmaps_radius) + #ifdef CHECKTIMING || find_string_argument(i, clo_EXCLUSIVETIMING, rawExclusiveTiming) #endif @@ -663,9 +669,21 @@ bool Environment::parse_cmdline(int argc, char** argv) if (!raw_coarse_grayscale_depth.empty()) { // string -> integer - if (sscanf(raw_coarse_grayscale_depth.c_str(), "%d", &coarse_grayscale_depth) != 1) + if (sscanf(raw_coarse_grayscale_depth.c_str(), "%d", &coarse_grayscale_depth) != 1 || coarse_grayscale_depth < 1) + { + std::cerr << "Error: " << clo_COARSEGRAYDEPTH << "=" << raw_coarse_grayscale_depth << ": expecting a positive integer constant\n"; + return false; + } + } + + // parse BINNINGORIGIN -- negate coarse_grayscale_depth for "min" origin (PyRadiomics-style) + if (!raw_binning_origin.empty()) + { + if (raw_binning_origin == "min") + coarse_grayscale_depth = -std::abs(coarse_grayscale_depth); + else if (raw_binning_origin != "zero") { - std::cerr << "Error: " << clo_COARSEGRAYDEPTH << "=" << raw_coarse_grayscale_depth << ": expecting an integer constant\n"; + std::cerr << "Error: " << clo_BINNINGORIGIN << "=" << raw_binning_origin << ": expecting 'zero' or 'min'\n"; return false; } } @@ -900,6 +918,25 @@ bool Environment::parse_cmdline(int argc, char** argv) ibsi_compliance = false; } + //==== Parse feature maps options + if (!raw_fmaps.empty()) + { + std::string tmp = raw_fmaps; + std::transform(tmp.begin(), tmp.end(), tmp.begin(), ::tolower); + fmaps_mode = (tmp == "true" || tmp == "1" || tmp == "on"); + } + + if (!raw_fmaps_radius.empty()) + { + int r = 0; + if (sscanf(raw_fmaps_radius.c_str(), "%d", &r) != 1 || r < 1) + { + std::cerr << "Error: " << clo_FMAPS_RADIUS << "=" << raw_fmaps_radius << ": expecting an integer >= 1\n"; + return false; + } + fmaps_kernel_radius = r; + } + // Success return true; } @@ -915,7 +952,7 @@ int Environment::get_coarse_gray_depth() return coarse_grayscale_depth; } -void Environment::set_coarse_gray_depth(unsigned int new_depth) +void Environment::set_coarse_gray_depth(int new_depth) { coarse_grayscale_depth = new_depth; } diff --git a/src/nyx/environment.h b/src/nyx/environment.h index c089a898..28c0f8d0 100644 --- a/src/nyx/environment.h +++ b/src/nyx/environment.h @@ -112,7 +112,7 @@ class Environment: public BasicEnvironment int get_floating_point_precision(); int get_coarse_gray_depth(); - void set_coarse_gray_depth(unsigned int new_depth); + void set_coarse_gray_depth(int new_depth); // implementation of SKIPROI bool roi_is_blacklisted (const std::string& fname, int roi_label); @@ -148,6 +148,21 @@ class Environment: public BasicEnvironment ResultOptions resultOptions; std::tuple> parse_result_options_4cli (); + // feature maps options + bool fmaps_mode = false; + int fmaps_kernel_radius = 2; + std::string raw_fmaps; + std::string raw_fmaps_radius; + + /// @brief Returns the kernel side length: 2*radius+1 + int fmaps_kernel_size() const { return 2 * fmaps_kernel_radius + 1; } + + /// @brief Returns true if fmaps mode conflicts with the current save option. + bool fmaps_prevents_arrow() const + { + return fmaps_mode && (saveOption == Nyxus::SaveOption::saveArrowIPC || saveOption == Nyxus::SaveOption::saveParquet); + } + // feature settings Fsettings fsett_PixelIntensity, fsett_BasicMorphology, @@ -242,6 +257,7 @@ class Environment: public BasicEnvironment int coarse_grayscale_depth; //= 64; std::string raw_coarse_grayscale_depth; //= ""; + std::string raw_binning_origin; //= "" (default "zero"; alternative "min") // data members implementing RAMLIMIT std::string rawRamLimit; //= ""; diff --git a/src/nyx/features/3d_glcm.cpp b/src/nyx/features/3d_glcm.cpp index 505c4dce..7c34d1fa 100644 --- a/src/nyx/features/3d_glcm.cpp +++ b/src/nyx/features/3d_glcm.cpp @@ -48,8 +48,11 @@ void D3_GLCM_feature::calculate (LR& r, const Fsettings& s) SimpleCube D; D.allocate (w,h,d); - auto greyBinningInfo = STNGS_GLCM_GREYDEPTH(s); // former Nyxus::theEnvironment.get_coarse_gray_depth() - if (STNGS_IBSI(s)) // former Nyxus::theEnvironment.ibsi_compliance + // Use GLCM-specific grey depth if set via metaparams, otherwise fall back to global + auto greyBinningInfo = STNGS_GLCM_GREYDEPTH(s); + if (greyBinningInfo == 0) + greyBinningInfo = s[(int)NyxSetting::GREYDEPTH].ival; + if (STNGS_IBSI(s)) greyBinningInfo = 0; bin_intensities_3d (D, r.aux_image_cube, r.aux_min, r.aux_max, greyBinningInfo); @@ -64,7 +67,7 @@ void D3_GLCM_feature::calculate (LR& r, const Fsettings& s) D, r.aux_min, r.aux_max, - STNGS_GLCM_GREYDEPTH(s), + greyBinningInfo, STNGS_IBSI(s), STNGS_NAN(s)); } diff --git a/src/nyx/features/3d_ngtdm.cpp b/src/nyx/features/3d_ngtdm.cpp index cf6dce21..ff0290fa 100644 --- a/src/nyx/features/3d_ngtdm.cpp +++ b/src/nyx/features/3d_ngtdm.cpp @@ -299,7 +299,7 @@ void D3_NGTDM_feature::osized_calculate (LR& r, const Fsettings& s, ImageLoader& D.allocate_from_cloud(r.raw_pixels_NT, r.aabb, false); // Gather zones - unsigned int nGrays = STNGS_NGTDM_GREYDEPTH(s); // former theEnvironment.get_coarse_gray_depth() + int nGrays = STNGS_NGTDM_GREYDEPTH(s); // former theEnvironment.get_coarse_gray_depth() for (int row = 0; row < D.get_height(); row++) for (int col = 0; col < D.get_width(); col++) { diff --git a/src/nyx/features/glcm.cpp b/src/nyx/features/glcm.cpp index 11748ff4..bb12a521 100644 --- a/src/nyx/features/glcm.cpp +++ b/src/nyx/features/glcm.cpp @@ -18,10 +18,11 @@ void GLCMFeature::calculate (LR& r, const Fsettings& s) // Clear the feature values buffers clear_result_buffers(); - // Skip feature calculation in case of bad data - // (We need to smart-select the greyInfo rather than just theEnvironment.get_coarse_gray_depth()) - int nGreys = STNGS_GLCM_GREYDEPTH(s), - offset = STNGS_GLCM_OFFSET(s); + // Use GLCM-specific grey depth if set via metaparams, otherwise fall back to global + int nGreys = STNGS_GLCM_GREYDEPTH(s); + if (nGreys == 0) + nGreys = s[(int)NyxSetting::GREYDEPTH].ival; + int offset = STNGS_GLCM_OFFSET(s); double softNAN = s[(int)NyxSetting::SOFTNAN].rval; auto binnedMin = bin_pixel(r.aux_min, r.aux_min, r.aux_max, nGreys); diff --git a/src/nyx/features/gldm.cpp b/src/nyx/features/gldm.cpp index 2378ef82..493ba2bc 100644 --- a/src/nyx/features/gldm.cpp +++ b/src/nyx/features/gldm.cpp @@ -274,7 +274,7 @@ void GLDMFeature::osized_calculate (LR& r, const Fsettings& s, ImageLoader&) // Prepare ROI's intensity range for normalize_I() PixIntens piRange = r.aux_max - r.aux_min; - unsigned int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() + int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() size_t height = D.get_height(), width = D.get_width(); diff --git a/src/nyx/features/gldzm.cpp b/src/nyx/features/gldzm.cpp index 49410823..9b094ae2 100644 --- a/src/nyx/features/gldzm.cpp +++ b/src/nyx/features/gldzm.cpp @@ -546,7 +546,7 @@ void GLDZMFeature::osized_calculate (LR& r, const Fsettings& s, ImageLoader&) // -- Squeeze pixels' intensity range for getting more prominent zones PixIntens piRange = r.aux_max - 0; // reflecting the fact that the original image's pixel intensity range is [0-r.aux_max] where 0 represents off-ROI pixels - unsigned int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() + int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() for (size_t i = 0; i < D.size(); i++) { // raw intensity diff --git a/src/nyx/features/glrlm.cpp b/src/nyx/features/glrlm.cpp index 4d59f803..592a22f2 100644 --- a/src/nyx/features/glrlm.cpp +++ b/src/nyx/features/glrlm.cpp @@ -368,7 +368,7 @@ void GLRLMFeature::osized_calculate (LR& r, const Fsettings& s, ImageLoader&) D.copy(im); // Squeeze the intensity range - unsigned int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() + int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() for (size_t i = 0; i < D.size(); i++) D.set_at(i, Nyxus::to_grayscale(D[i], r.aux_min, piRange, nGrays, STNGS_IBSI(s))); diff --git a/src/nyx/features/glszm.cpp b/src/nyx/features/glszm.cpp index d37e5a87..2c59231d 100644 --- a/src/nyx/features/glszm.cpp +++ b/src/nyx/features/glszm.cpp @@ -73,7 +73,7 @@ void GLSZMFeature::osized_calculate (LR& r, const Fsettings& s, ImageLoader&) // Squeeze the intensity range PixIntens piRange = r.aux_max - r.aux_min; // Prepare ROI's intensity range - unsigned int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() + int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() for (size_t i = 0; i < D.size(); i++) D.set_at(i, Nyxus::to_grayscale(D[i], r.aux_min, piRange, nGrays, STNGS_IBSI(s))); diff --git a/src/nyx/features/ngtdm.cpp b/src/nyx/features/ngtdm.cpp index 63de8971..e179c74e 100644 --- a/src/nyx/features/ngtdm.cpp +++ b/src/nyx/features/ngtdm.cpp @@ -237,7 +237,7 @@ void NGTDMFeature::osized_calculate (LR& r, const Fsettings& s, ImageLoader&) D.allocate_from_cloud (r.raw_pixels_NT, r.aabb, false); // Gather zones - unsigned int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() + int nGrays = STNGS_NGREYS(s); // former theEnvironment.get_coarse_gray_depth() for (int row = 0; row < D.get_height(); row++) for (int col = 0; col < D.get_width(); col++) { diff --git a/src/nyx/globals.h b/src/nyx/globals.h index 02c4896f..4e015edf 100644 --- a/src/nyx/globals.h +++ b/src/nyx/globals.h @@ -49,6 +49,24 @@ namespace Nyxus const SaveOption saveOption, const std::string& outputPath); + // feature maps 2D workflow + int processDataset_2D_fmaps( + Environment& env, + const std::vector& intensFiles, + const std::vector& labelFiles, + int numReduceThreads, + const SaveOption saveOption, + const std::string& outputPath); + + // feature maps 3D workflow + int processDataset_3D_fmaps( + Environment& env, + const std::vector & intensFiles, + const std::vector & labelFiles, + int numReduceThreads, + const SaveOption saveOption, + const std::string& outputPath); + // single-segment 2D workflow int processDataset_2D_wholeslide( Environment & env, @@ -75,6 +93,10 @@ namespace Nyxus const SaveOption saveOption, const std::string& outputPath); + // Reset CSV header state between workflow invocations (fixes static flag persistence). + // Must be called before any worker threads begin processing (not thread-safe). + void reset_csv_header_state(); + std::string getPureFname(const std::string& fpath); bool gatherRoisMetrics(int slide_idx, const std::string& intens_fpath, const std::string& label_fpath, Environment & env, ImageLoader & L); bool gather_wholeslide_metrics(const std::string& intens_fpath, ImageLoader& L, LR& roi); @@ -90,7 +112,10 @@ namespace Nyxus bool scan_trivial_wholevolume (LR& vroi, const std::string& intens_fpath, ImageLoader& ldr); bool scan_trivial_wholevolume_anisotropic (LR& vroi, const std::string& intens_fpath, ImageLoader& ldr, double aniso_x, double aniso_y, double aniso_z); + bool scanTrivialRois (const std::vector& batch_labels, const std::string& intens_fpath, const std::string& label_fpath, Environment & env, ImageLoader & ldr); + bool scanTrivialRois_anisotropic (const std::vector& batch_labels, const std::string& intens_fpath, const std::string& label_fpath, Environment & env, ImageLoader & ldr, double sf_x, double sf_y); bool scanTrivialRois_3D (Environment& env, const std::vector& batch_labels, const std::string& intens_fpath, const std::string& label_fpath, size_t t_index); + bool scanTrivialRois_3D_anisotropic (Environment& env, const std::vector& batch_labels, const std::string& intens_fpath, const std::string& label_fpath, size_t t_index, double aniso_x, double aniso_y, double aniso_z); void dump_roi_metrics (const int dim, const std::string& output_dir, const size_t ram_limit, const std::string& seg_fpath, const Uniqueids& uniqueLabels, const Roidata& roiData); void dump_roi_pixels (const int dim, const std::string& output_dir, const std::vector& batch_labels, const std::string& seg_fpath, const Uniqueids& uniqueLabels, const Roidata& roiData); void dump_2d_image_with_halfcontour( @@ -135,7 +160,95 @@ namespace Nyxus extern const std::vector mandatory_output_columns; bool save_features_2_csv (Environment & env, const std::string & intFpath, const std::string & segFpath, const std::string & outputDir, size_t t_index, bool need_aggregation); bool save_features_2_csv_wholeslide (Environment & env, const LR & r, const std::string & ifpath, const std::string & mfpath, const std::string & outdir, size_t t_index); + + // Shared helpers for feature column names and value collection (eliminates duplication across output paths) + std::vector get_feature_column_names (Environment & env, const std::vector> & F); + std::vector collect_feature_values (const LR & r, const std::vector> & F); + + /// @brief RAII guard that temporarily replaces env's ROI data with child data and restores on destruction. + /// Ensures env is restored even if an exception is thrown during child ROI processing. + struct EnvRoiSwapGuard + { + Environment & env; + Uniqueids savedUniqueLabels; + Roidata savedRoiData; + + EnvRoiSwapGuard (Environment & e, std::unordered_set childLabels, std::unordered_map childRoiData) + : env(e) + { + savedUniqueLabels = std::move(env.uniqueLabels); + savedRoiData = std::move(env.roiData); + env.uniqueLabels = std::move(childLabels); + env.roiData = std::move(childRoiData); + } + + ~EnvRoiSwapGuard () + { + env.uniqueLabels = std::move(savedUniqueLabels); + env.roiData = std::move(savedRoiData); + } + + EnvRoiSwapGuard (const EnvRoiSwapGuard &) = delete; + EnvRoiSwapGuard & operator= (const EnvRoiSwapGuard &) = delete; + }; + + // Feature maps + struct FmapChildInfo + { + int parent_label; + int center_x; + int center_y; + int center_z = 0; // unused in 2D + }; + + int generateChildRois ( + const LR & parent, + int kernel_size, + std::unordered_set & childLabels, + std::unordered_map & childRoiData, + std::unordered_map & childToParentMap, + int64_t startLabel); + int generateChildRois_3D ( + const LR & parent, + int kernel_size, + std::unordered_set & childLabels, + std::unordered_map & childRoiData, + std::unordered_map & childToParentMap, + int64_t startLabel); + bool save_features_2_csv_fmaps ( + Environment & env, + const std::string & intFpath, + const std::string & segFpath, + const std::string & outputDir, + int slide_idx, + const std::unordered_set & childLabels, + const std::unordered_map & childRoiData, + const std::unordered_map & childToParentMap); bool save_features_2_buffer (ResultsCache &results_cache, Environment &env, size_t t_index); + void save_features_2_fmap_arrays ( + ResultsCache & rescache, + Environment & env, + const std::string & intens_name, + const std::string & seg_name, + int parent_label, + int parent_xmin, int parent_ymin, + int parent_w, int parent_h, + int kernel_size, + const std::unordered_set & childLabels, + const std::unordered_map & childRoiData, + const std::unordered_map & childToParentMap); + void save_features_2_fmap_arrays_3D ( + ResultsCache & rescache, + Environment & env, + const std::string & intens_name, + const std::string & seg_name, + int parent_label, + int parent_xmin, int parent_ymin, int parent_zmin, + int parent_w, int parent_h, int parent_d, + int kernel_size, + const std::unordered_set & childLabels, + const std::unordered_map & childRoiData, + const std::unordered_map & childToParentMap); bool save_features_2_buffer_wholeslide( ResultsCache& rescache, Environment& env, diff --git a/src/nyx/helpers/helpers.h b/src/nyx/helpers/helpers.h index 9b0f4f7f..03b4537d 100644 --- a/src/nyx/helpers/helpers.h +++ b/src/nyx/helpers/helpers.h @@ -334,12 +334,13 @@ namespace Nyxus /// @param min_i Minimum ROI's intensity /// @param i_range Precalculated ROI's intensity range (= max-min) /// @return Squeezed intensity within range [0,255] - inline unsigned int to_grayscale (unsigned int i, unsigned int min_i, unsigned int i_range, unsigned int n_levels, bool disable_binning=false) + inline unsigned int to_grayscale (unsigned int i, unsigned int min_i, unsigned int i_range, int n_levels, bool disable_binning=false) { - if (disable_binning) + if (disable_binning) return i; - - double pi = ((double(i-min_i) / double(i_range) * double(n_levels))); + + unsigned int abs_n = (unsigned int)std::abs(n_levels); + double pi = ((double(i-min_i) / double(i_range) * double(abs_n))); unsigned int new_pi = (unsigned int)pi; return new_pi; } diff --git a/src/nyx/main_nyxus.cpp b/src/nyx/main_nyxus.cpp index fb314d8d..9cd32b14 100644 --- a/src/nyx/main_nyxus.cpp +++ b/src/nyx/main_nyxus.cpp @@ -80,7 +80,22 @@ int main (int argc, char** argv) // Process the image data int errorCode = 0; - if (env.singleROI) + if (env.fmaps_prevents_arrow()) + { + std::cerr << "Error: Arrow/Parquet output is not supported in feature maps (fmaps) mode. Use CSV output instead.\n"; + return 1; + } + if (env.fmaps_mode) + { + errorCode = processDataset_2D_fmaps( + env, + intensFiles, + labelFiles, + env.n_reduce_threads, + env.saveOption, + env.output_dir); + } + else if (env.singleROI) { errorCode = processDataset_2D_wholeslide( env, @@ -181,13 +196,32 @@ int main (int argc, char** argv) return 1; } - int errorCode = processDataset_3D_segmented( - env, - intensFiles, - labelFiles, - env.n_reduce_threads, - env.saveOption, - env.output_dir); + int errorCode = 0; + if (env.fmaps_prevents_arrow()) + { + std::cerr << "Error: Arrow/Parquet output is not supported in feature maps (fmaps) mode. Use CSV output instead.\n"; + return 1; + } + if (env.fmaps_mode) + { + errorCode = processDataset_3D_fmaps( + env, + intensFiles, + labelFiles, + env.n_reduce_threads, + env.saveOption, + env.output_dir); + } + else + { + errorCode = processDataset_3D_segmented( + env, + intensFiles, + labelFiles, + env.n_reduce_threads, + env.saveOption, + env.output_dir); + } // Report feature extraction error, if any switch (errorCode) diff --git a/src/nyx/output_2_buffer.cpp b/src/nyx/output_2_buffer.cpp index 7c7f1348..26f2d760 100644 --- a/src/nyx/output_2_buffer.cpp +++ b/src/nyx/output_2_buffer.cpp @@ -23,6 +23,31 @@ namespace Nyxus { static std::mutex mx1; + /// @brief Writes feature header columns for all enabled features to a ResultsCache. + /// Uses get_feature_column_names() as the single source of truth. + void write_feature_header_buffer ( + ResultsCache & rescache, + Environment & env, + const std::vector> & F) + { + auto cols = get_feature_column_names(env, F); + for (const auto& c : cols) + rescache.add_to_header(c); + } + + /// @brief Writes all feature values for a single ROI to a ResultsCache. + /// Uses collect_feature_values() as the single source of truth. + void write_feature_values_buffer ( + ResultsCache & rescache, + const LR & r, + const std::vector> & F, + Environment & env) + { + auto vals = collect_feature_values(r, F); + for (auto v : vals) + rescache.add_numeric(Nyxus::force_finite_number(v, env.resultOptions.noval())); + } + bool save_features_2_buffer_wholeslide ( // out ResultsCache & rescache, @@ -45,119 +70,7 @@ namespace Nyxus if (fill_header) { rescache.add_to_header({ Nyxus::colname_intensity_image, Nyxus::colname_mask_image, Nyxus::colname_roi_label, Nyxus::colname_t_index }); - - for (auto& enabdF : F) - { - auto fn = std::get<0>(enabdF); // feature name - auto fc = std::get<1>(enabdF); // feature code - - // Handle missing feature name (which is a significant issue!) in order to at least be able to trace back to the feature code - if (fn.empty()) - fn = "feature" + std::to_string(fc); - - // Parameterized feature - // --GLCM family - bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); - bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glcmFeature && nonAngledGlcmFeature == false) - { - // Populate with angles - for (auto ang : env.glcmOptions.glcmAngles) - { - std::string col = fn + "_" + std::to_string(ang); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - // --GLRLM family - bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); - bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glrlmFeature && nonAngledGlrlmFeature == false) - { - // Populate with angles - for (auto ang : GLRLMFeature::rotAngles) - { - std::string col = fn + "_" + std::to_string(ang); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - // --Gabor - if (fc == (int)Feature2D::GABOR) - { - // Generate the feature value list - for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) - { - std::string col = fn + "_" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - if (fc == (int)Feature2D::FRAC_AT_D) - { - // Generate the feature value list - for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) - { - std::string col = fn + "_" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - if (fc == (int)Feature2D::MEAN_FRAC) - { - // Generate the feature value list - for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) - { - std::string col = fn + "_" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - if (fc == (int)Feature2D::RADIAL_CV) - { - // Generate the feature value list - for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) - { - std::string col = fn + "_" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - // --Zernike family - if (fc == (int)Feature2D::ZERNIKE2D) - { - // Populate with indices - for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) - { - std::string col = fn + "_Z" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - // Regular feature - rescache.add_to_header(fn); - } + write_feature_header_buffer(rescache, env, F); } // -- Values @@ -171,103 +84,7 @@ namespace Nyxus rescache.add_numeric (DEFAULT_T_INDEX); // - features - for (auto& enabdF : F) - { - auto fc = std::get<1>(enabdF); - auto fn = std::get<0>(enabdF); // debug - - auto vv = r.get_fvals(fc); - - // Parameterized feature - // --GLCM family - bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); - bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glcmFeature && nonAngledGlcmFeature == false) - { - // Mock angled values if they haven't been calculated for some error reason - if (vv.size() < GLCMFeature::angles.size()) - vv.resize(GLCMFeature::angles.size(), 0.0); - - // Populate with angles - int nAng = GLCMFeature::angles.size(); - for (int i = 0; i < nAng; i++) - { - rescache.add_numeric(Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - - // Proceed with other features - continue; - } - - // --GLRLM family - bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); - bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glrlmFeature && nonAngledGlrlmFeature == false) - { - // Populate with angles - int nAng = 4; // check GLRLMFeature::rotAngles - for (int i = 0; i < nAng; i++) - { - rescache.add_numeric(Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - - // --Gabor - if (fc == (int)Feature2D::GABOR) - { - for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) - { - rescache.add_numeric(Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - - // --Zernike family - if (fc == (int)Feature2D::ZERNIKE2D) - { - for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) - { - rescache.add_numeric(Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - - // --Radial distribution features - if (fc == (int)Feature2D::FRAC_AT_D) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) - { - rescache.add_numeric(Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - if (fc == (int)Feature2D::MEAN_FRAC) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) - { - rescache.add_numeric(Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - if (fc == (int)Feature2D::RADIAL_CV) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) - { - rescache.add_numeric(Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - - // Regular feature - rescache.add_numeric(Nyxus::force_finite_number(vv[0], env.resultOptions.noval())); - } + write_feature_values_buffer(rescache, r, F, env); return true; } @@ -287,119 +104,7 @@ namespace Nyxus if (fill_header) { rescache.add_to_header({ Nyxus::colname_intensity_image, Nyxus::colname_mask_image, Nyxus::colname_roi_label, Nyxus::colname_t_index }); - - for (auto& enabdF : F) - { - auto fn = std::get<0>(enabdF); // feature name - auto fc = std::get<1>(enabdF); // feature code - - // Handle missing feature name (which is a significant issue!) in order to at least be able to trace back to the feature code - if (fn.empty()) - fn = "feature" + std::to_string(fc); - - // Parameterized feature - // --GLCM family - bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); - bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glcmFeature && nonAngledGlcmFeature == false) - { - // Populate with angles - for (auto ang : env.glcmOptions.glcmAngles) - { - std::string col = fn + "_" + std::to_string(ang); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - // --GLRLM family - bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); - bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glrlmFeature && nonAngledGlrlmFeature == false) - { - // Populate with angles - for (auto ang : GLRLMFeature::rotAngles) - { - std::string col = fn + "_" + std::to_string(ang); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - // --Gabor - if (fc == (int) Feature2D::GABOR) - { - // Generate the feature value list - for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) - { - std::string col = fn + "_" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - if (fc == (int) Feature2D::FRAC_AT_D) - { - // Generate the feature value list - for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) - { - std::string col = fn + "_" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - if (fc == (int) Feature2D::MEAN_FRAC) - { - // Generate the feature value list - for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) - { - std::string col = fn + "_" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - if (fc == (int) Feature2D::RADIAL_CV) - { - // Generate the feature value list - for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) - { - std::string col = fn + "_" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - // --Zernike family - if (fc == (int) Feature2D::ZERNIKE2D) - { - // Populate with indices - for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) - { - std::string col = fn + "_Z" + std::to_string(i); - rescache.add_to_header(col); - } - - // Proceed with other features - continue; - } - - // Regular feature - rescache.add_to_header(fn); - } + write_feature_header_buffer(rescache, env, F); } // -- Values @@ -415,7 +120,7 @@ namespace Nyxus // Tear off pure file names from segment and intensity file paths const SlideProps & slide = env.dataset.dataset_props[r.slide_idx]; - fs::path pseg(slide.fname_seg), + fs::path pseg(slide.fname_seg), pint(slide.fname_int); std::string segfname = pseg.filename().string(), intfname = pint.filename().string(); @@ -425,106 +130,162 @@ namespace Nyxus rescache.add_numeric (l); rescache.add_numeric (t_index); - for (auto& enabdF : F) - { - auto fc = std::get<1>(enabdF); - auto fn = std::get<0>(enabdF); // debug - - auto vv = r.get_fvals(fc); - - // Parameterized feature - // --GLCM family - bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); - bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glcmFeature && nonAngledGlcmFeature == false) - { - // Mock angled values if they haven't been calculated for some error reason - if (vv.size() < GLCMFeature::angles.size()) - vv.resize(GLCMFeature::angles.size(), 0.0); - - // Populate with angles - int nAng = GLCMFeature::angles.size(); - for (int i = 0; i < nAng; i++) - { - rescache.add_numeric (Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - - // Proceed with other features - continue; - } - - // --GLRLM family - bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); - bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glrlmFeature && nonAngledGlrlmFeature == false) - { - // Populate with angles - int nAng = 4; // check GLRLMFeature::rotAngles - for (int i = 0; i < nAng; i++) - { - rescache.add_numeric (Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - - // --Gabor - if (fc == (int) Feature2D::GABOR) - { - for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) - { - rescache.add_numeric (Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - - // --Zernike family - if (fc == (int) Feature2D::ZERNIKE2D) - { - for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) - { - rescache.add_numeric (Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - - // --Radial distribution features - if (fc == (int) Feature2D::FRAC_AT_D) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) - { - rescache.add_numeric (Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - if (fc == (int) Feature2D::MEAN_FRAC) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) - { - rescache.add_numeric (Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - if (fc == (int) Feature2D::RADIAL_CV) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) - { - rescache.add_numeric (Nyxus::force_finite_number(vv[i], env.resultOptions.noval())); - } - // Proceed with other features - continue; - } - - // Regular feature - rescache.add_numeric (Nyxus::force_finite_number(vv[0], env.resultOptions.noval())); - } + write_feature_values_buffer(rescache, r, F, env); } return true; } + void save_features_2_fmap_arrays ( + ResultsCache & rescache, + Environment & env, + const std::string & intens_name, + const std::string & seg_name, + int parent_label, + int parent_xmin, int parent_ymin, + int parent_w, int parent_h, + int kernel_size, + const std::unordered_set & childLabels, + const std::unordered_map & childRoiData, + const std::unordered_map & childToParentMap) + { + int half = kernel_size / 2; + int map_w = parent_w - kernel_size + 1; + int map_h = parent_h - kernel_size + 1; + + if (map_w <= 0 || map_h <= 0) + return; + + std::vector> F = env.theFeatureSet.getEnabledFeatures(); + std::vector feature_names = get_feature_column_names(env, F); + size_t n_features = feature_names.size(); + size_t map_size = (size_t)map_w * map_h; + + // Initialize with NaN + std::vector feature_data(n_features * map_size, std::numeric_limits::quiet_NaN()); + + // Fill in values from child ROIs + for (auto l : childLabels) + { + auto it_roi = childRoiData.find(l); + if (it_roi == childRoiData.end()) + continue; + const LR& r = it_roi->second; + if (r.blacklisted) + continue; + + auto it_map = childToParentMap.find(l); + if (it_map == childToParentMap.end()) + continue; + + const FmapChildInfo& info = it_map->second; + + // Convert global center coords to map indices + int map_row = info.center_y - parent_ymin - half; + int map_col = info.center_x - parent_xmin - half; + + if (map_row < 0 || map_row >= map_h || map_col < 0 || map_col >= map_w) + continue; + + size_t pixel_idx = (size_t)map_row * map_w + map_col; + + // Collect feature values for this child + auto vals = collect_feature_values(r, F); + for (size_t fi = 0; fi < n_features && fi < vals.size(); fi++) + feature_data[fi * map_size + pixel_idx] = force_finite_number(vals[fi], env.resultOptions.noval()); + } + + // Store result + FmapArrayResult result; + result.parent_label = parent_label; + result.intens_name = intens_name; + result.seg_name = seg_name; + result.map_w = map_w; + result.map_h = map_h; + result.origin_x = parent_xmin + half; + result.origin_y = parent_ymin + half; + result.feature_names = std::move(feature_names); + result.feature_data = std::move(feature_data); + rescache.get_fmapArrayResults().push_back(std::move(result)); + } + + void save_features_2_fmap_arrays_3D ( + ResultsCache & rescache, + Environment & env, + const std::string & intens_name, + const std::string & seg_name, + int parent_label, + int parent_xmin, int parent_ymin, int parent_zmin, + int parent_w, int parent_h, int parent_d, + int kernel_size, + const std::unordered_set & childLabels, + const std::unordered_map & childRoiData, + const std::unordered_map & childToParentMap) + { + int half = kernel_size / 2; + int map_w = parent_w - kernel_size + 1; + int map_h = parent_h - kernel_size + 1; + int map_d = parent_d - kernel_size + 1; + + if (map_w <= 0 || map_h <= 0 || map_d <= 0) + return; + + std::vector> F = env.theFeatureSet.getEnabledFeatures(); + std::vector feature_names = get_feature_column_names(env, F); + size_t n_features = feature_names.size(); + size_t map_size = (size_t)map_d * map_h * map_w; + + // Initialize with NaN + std::vector feature_data(n_features * map_size, std::numeric_limits::quiet_NaN()); + + // Fill in values from child ROIs + for (auto l : childLabels) + { + auto it_roi = childRoiData.find(l); + if (it_roi == childRoiData.end()) + continue; + const LR& r = it_roi->second; + if (r.blacklisted) + continue; + + auto it_map = childToParentMap.find(l); + if (it_map == childToParentMap.end()) + continue; + + const FmapChildInfo& info = it_map->second; + + // Convert global center coords to map indices + int map_col = info.center_x - parent_xmin - half; + int map_row = info.center_y - parent_ymin - half; + int map_z = info.center_z - parent_zmin - half; + + if (map_col < 0 || map_col >= map_w || + map_row < 0 || map_row >= map_h || + map_z < 0 || map_z >= map_d) + continue; + + size_t voxel_idx = (size_t)map_z * map_h * map_w + (size_t)map_row * map_w + map_col; + + // Collect feature values for this child + auto vals = collect_feature_values(r, F); + for (size_t fi = 0; fi < n_features && fi < vals.size(); fi++) + feature_data[fi * map_size + voxel_idx] = force_finite_number(vals[fi], env.resultOptions.noval()); + } + + // Store result + FmapArrayResult result; + result.parent_label = parent_label; + result.intens_name = intens_name; + result.seg_name = seg_name; + result.map_w = map_w; + result.map_h = map_h; + result.map_d = map_d; + result.origin_x = parent_xmin + half; + result.origin_y = parent_ymin + half; + result.origin_z = parent_zmin + half; + result.feature_names = std::move(feature_names); + result.feature_data = std::move(feature_data); + rescache.get_fmapArrayResults().push_back(std::move(result)); + } + } diff --git a/src/nyx/output_2_csv.cpp b/src/nyx/output_2_csv.cpp index b64ed936..c6c47519 100644 --- a/src/nyx/output_2_csv.cpp +++ b/src/nyx/output_2_csv.cpp @@ -1,6 +1,7 @@ +#include #include #include -#include +#include #include #include #include @@ -26,6 +27,211 @@ namespace Nyxus #define fopen_s(pFile,filename,mode) ((*(pFile))=fopen((filename),(mode)))==NULL #endif + /// @brief Returns feature column names for all enabled features. + /// Single source of truth for header generation across CSV, buffer, and Arrow output paths. + std::vector get_feature_column_names ( + Environment & env, + const std::vector> & F) + { + std::vector cols; + + for (const auto& enabdF : F) + { + auto fn = std::get<0>(enabdF); + auto fc = std::get<1>(enabdF); + + if (fn.empty()) + fn = "feature" + std::to_string(fc); + + // GLCM family + bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); + bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); + if (glcmFeature && nonAngledGlcmFeature == false) + { + for (auto ang : env.glcmOptions.glcmAngles) + cols.push_back(fn + "_" + std::to_string(ang)); + continue; + } + + // GLRLM family + bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); + bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); + if (glrlmFeature && nonAngledGlrlmFeature == false) + { + for (auto ang : GLRLMFeature::rotAngles) + cols.push_back(fn + "_" + std::to_string(ang)); + continue; + } + + // Gabor + if (fc == (int)Feature2D::GABOR) + { + for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) + cols.push_back(fn + "_" + std::to_string(i)); + continue; + } + + if (fc == (int)Feature2D::FRAC_AT_D) + { + for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) + cols.push_back(fn + "_" + std::to_string(i)); + continue; + } + + if (fc == (int)Feature2D::MEAN_FRAC) + { + for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) + cols.push_back(fn + "_" + std::to_string(i)); + continue; + } + + if (fc == (int)Feature2D::RADIAL_CV) + { + for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) + cols.push_back(fn + "_" + std::to_string(i)); + continue; + } + + // Zernike + if (fc == (int)Feature2D::ZERNIKE2D) + { + for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) + cols.push_back(fn + "_Z" + std::to_string(i)); + continue; + } + + // Regular feature + cols.push_back(fn); + } + + return cols; + } + + /// @brief Collects raw feature values for a single ROI into a flat vector. + /// Single source of truth for value collection across CSV, buffer, and Python output paths. + /// Callers apply force_finite_number or other formatting as needed. + std::vector collect_feature_values ( + const LR & r, + const std::vector> & F) + { + std::vector vals; + + for (const auto& enabdF : F) + { + auto fc = std::get<1>(enabdF); + auto vv = r.get_fvals(fc); + + // GLCM family + bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); + bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); + if (glcmFeature && nonAngledGlcmFeature == false) + { + if (vv.size() < GLCMFeature::angles.size()) + vv.resize(GLCMFeature::angles.size(), 0.0); + int nAng = (int)GLCMFeature::angles.size(); + for (int i = 0; i < nAng; i++) + vals.push_back(vv[i]); + continue; + } + + // GLRLM family + bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); + bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); + if (glrlmFeature && nonAngledGlrlmFeature == false) + { + int nAng = (int)std::size(GLRLMFeature::rotAngles); + if ((int)vv.size() < nAng) + vv.resize(nAng, 0.0); + for (int i = 0; i < nAng; i++) + vals.push_back(vv[i]); + continue; + } + + // Gabor + if (fc == (int)Feature2D::GABOR) + { + auto nGabor = GaborFeature::f0_theta_pairs.size(); + if (vv.size() < nGabor) + vv.resize(nGabor, 0.0); + for (size_t i = 0; i < nGabor; i++) + vals.push_back(vv[i]); + continue; + } + + // Zernike + if (fc == (int)Feature2D::ZERNIKE2D) + { + if ((int)vv.size() < ZernikeFeature::NUM_FEATURE_VALS) + vv.resize(ZernikeFeature::NUM_FEATURE_VALS, 0.0); + for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) + vals.push_back(vv[i]); + continue; + } + + // Radial distribution + if (fc == (int)Feature2D::FRAC_AT_D) + { + if ((int)vv.size() < RadialDistributionFeature::num_features_FracAtD) + vv.resize(RadialDistributionFeature::num_features_FracAtD, 0.0); + for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) + vals.push_back(vv[i]); + continue; + } + if (fc == (int)Feature2D::MEAN_FRAC) + { + if ((int)vv.size() < RadialDistributionFeature::num_features_MeanFrac) + vv.resize(RadialDistributionFeature::num_features_MeanFrac, 0.0); + for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) + vals.push_back(vv[i]); + continue; + } + if (fc == (int)Feature2D::RADIAL_CV) + { + if ((int)vv.size() < RadialDistributionFeature::num_features_RadialCV) + vv.resize(RadialDistributionFeature::num_features_RadialCV, 0.0); + for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) + vals.push_back(vv[i]); + continue; + } + + // Regular feature + vals.push_back(vv.empty() ? 0.0 : vv[0]); + } + + return vals; + } + + /// @brief Writes all feature values for a single ROI to a stringstream in CSV format. + /// Uses collect_feature_values() as the single source of truth for value dispatch. + void write_feature_values_csv ( + std::stringstream & ssVals, + const LR & r, + const std::vector> & F, + Environment & env, + char* rvbuf, + int rvbuf_len, + const char* rvfmt) + { + auto vals = collect_feature_values(r, F); + +#ifdef DIAGNOSE_NYXUS_OUTPUT + auto names = get_feature_column_names(env, F); + for (size_t i = 0; i < vals.size(); i++) + { + double fv = Nyxus::force_finite_number(vals[i], env.resultOptions.noval()); + snprintf(rvbuf, rvbuf_len, rvfmt, fv); + ssVals << "," << names[i] << "-" << rvbuf; + } +#else + for (auto v : vals) + { + double fv = Nyxus::force_finite_number(v, env.resultOptions.noval()); + snprintf(rvbuf, rvbuf_len, rvfmt, fv); + ssVals << "," << rvbuf; + } +#endif + } + double auto_precision (Environment & env, std::stringstream& ss, double x) { if (std::abs(x) >= 1.0) @@ -80,107 +286,10 @@ namespace Nyxus // time head.push_back (Nyxus::colname_t_index); - // Optional columns - for (auto& enabdF : F) - { - auto fn = std::get<0>(enabdF); // feature name - auto fc = std::get<1>(enabdF); // feature code - - // Handle missing feature name (which is a significant issue!) in order to at least be able to trace back to the feature code - if (fn.empty()) - { - std::stringstream temp; - temp << "feature" << fc; - fn = temp.str(); - } - - // Parameterized feature - // --GLCM family - bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); - bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glcmFeature && nonAngledGlcmFeature == false) - { - // Populate with angles - for (auto ang : env.glcmOptions.glcmAngles) - { - // CSV separator - //if (ang != env.rotAngles[0]) - // ssHead << ","; - head.emplace_back(fn + "_" + std::to_string(ang)); - } - // Proceed with other features - continue; - } - - // --GLRLM family - bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); - bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glrlmFeature && nonAngledGlrlmFeature == false) - { - // Populate with angles - for (auto ang : GLRLMFeature::rotAngles) - { - head.emplace_back(fn + "_" + std::to_string(ang)); - } - // Proceed with other features - continue; - } - - // --Gabor - if (fc == (int) Feature2D::GABOR) - { - // Generate the feature value list - for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) - head.emplace_back(fn + "_" + std::to_string(i)); - - // Proceed with other features - continue; - } - - if (fc == (int) Feature2D::FRAC_AT_D) - { - // Generate the feature value list - for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) - head.emplace_back(fn + "_" + std::to_string(i)); - - // Proceed with other features - continue; - } - - if (fc == (int) Feature2D::MEAN_FRAC) - { - // Generate the feature value list - for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) - head.emplace_back(fn + "_" + std::to_string(i)); - - // Proceed with other features - continue; - } - - if (fc == (int) Feature2D::RADIAL_CV) - { - // Generate the feature value list - for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) - head.emplace_back(fn + "_" + std::to_string(i)); - - // Proceed with other features - continue; - } - - // --Zernike features header - if (fc == (int) Feature2D::ZERNIKE2D) - { - // Populate with indices - for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) // i < ZernikeFeature::num_feature_values_calculated - head.emplace_back(fn + "_Z" + std::to_string(i)); - - // Proceed with other features - continue; - } - - // Regular feature - head.emplace_back(fn); - } + // Feature columns (single source of truth) + auto feature_cols = get_feature_column_names(env, F); + for (const auto& c : feature_cols) + head.push_back(c); return head; } @@ -197,6 +306,19 @@ namespace Nyxus static std::mutex mutex1; + // Namespace-level header flags (resettable between workflow invocations). + // std::atomic ensures safe reset-before-use even without mutex protection. + static std::atomic mustRenderHeader_wholeslide{true}; + static std::atomic mustRenderHeader_segmented{true}; + static std::atomic mustRenderHeader_fmaps{true}; + + void reset_csv_header_state() + { + mustRenderHeader_wholeslide.store(true); + mustRenderHeader_segmented.store(true); + mustRenderHeader_fmaps.store(true); + } + bool save_features_2_csv_wholeslide ( Environment & env, const LR & r, @@ -212,8 +334,6 @@ namespace Nyxus char rvbuf[VAL_BUF_LEN]; // real value buffer large enough to fit a float64 value in range (2.2250738585072014E-308 to 1.79769313486231570e+308) const char rvfmt[] = "%g"; // instead of "%20.12f" which produces a too massive output - static bool mustRenderHeader = true; // In 'singlecsv' scenario this flag flips from 'T' to 'F' when necessary (after the header is rendered) - // Make the file name and write mode std::string fullPath = get_feature_output_fname (env, ifpath, mfpath); VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\t--> " << fullPath << "\n"); @@ -221,7 +341,7 @@ namespace Nyxus // Single CSV: create or continue? const char* mode = "w"; if (!env.separateCsv) - mode = mustRenderHeader ? "w" : "a"; + mode = mustRenderHeader_wholeslide ? "w" : "a"; // Open it FILE* fp = nullptr; @@ -235,7 +355,7 @@ namespace Nyxus } // configure buffered write - if (std::setvbuf(fp, nullptr, _IOFBF, 32768) != 0) + if (std::setvbuf(fp, nullptr, _IOFBF, 32768) != 0) std::cerr << "setvbuf failed \n"; // Learn what features need to be displayed @@ -243,7 +363,7 @@ namespace Nyxus // ********** Header - if (mustRenderHeader) + if (mustRenderHeader_wholeslide) { std::stringstream ssHead; @@ -263,7 +383,7 @@ namespace Nyxus // Prevent rendering the header again for another image's portion of labels if (env.separateCsv == false) - mustRenderHeader = false; + mustRenderHeader_wholeslide = false; } // ********** Values @@ -282,117 +402,7 @@ namespace Nyxus ssVals << "," << r.label << "," << t_index; // features - for (auto& enabdF : F) - { - auto fc = std::get<1>(enabdF); - auto fn = std::get<0>(enabdF); // feature name, for debugging - auto vv = r.get_fvals(fc); - - // Parameterized feature - // --GLCM family - bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); - bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glcmFeature && nonAngledGlcmFeature == false) - { - // Mock angled values if they haven't been calculated for some error reason - if (vv.size() < GLCMFeature::angles.size()) - vv.resize(GLCMFeature::angles.size(), 0.0); - // Output the sub-values - int nAng = (int)GLCMFeature::angles.size(); - for (int i = 0; i < nAng; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - ssVals << "," << rvbuf; - } - // Proceed with other features - continue; - } - - // --GLRLM family - bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); - bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glrlmFeature && nonAngledGlrlmFeature == false) - { - // Populate with angles - int nAng = 4; - for (int i = 0; i < nAng; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - ssVals << "," << rvbuf; - } - // Proceed with other features - continue; - } - - // --Gabor - if (fc == (int)Feature2D::GABOR) - { - for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - ssVals << "," << rvbuf; - } - - // Proceed with other features - continue; - } - - // --Zernike feature values - if (fc == (int)Feature2D::ZERNIKE2D) - { - for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - ssVals << "," << rvbuf; - } - - // Proceed with other features - continue; - } - - // --Radial distribution features - if (fc == (int)Feature2D::FRAC_AT_D) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - ssVals << "," << rvbuf; - } - // Proceed with other features - continue; - } - if (fc == (int)Feature2D::MEAN_FRAC) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - ssVals << "," << rvbuf; - } - // Proceed with other features - continue; - } - if (fc == (int)Feature2D::RADIAL_CV) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf (rvbuf, VAL_BUF_LEN, rvfmt, fv); - ssVals << "," << rvbuf; - } - // Proceed with other features - continue; - } - - // Regular feature - snprintf (rvbuf, VAL_BUF_LEN, rvfmt, Nyxus::force_finite_number(vv[0], env.resultOptions.noval())); - ssVals << "," << rvbuf; // Alternatively: auto_precision(ssVals, vv[0]); - } + write_feature_values_csv(ssVals, r, F, env, rvbuf, VAL_BUF_LEN, rvfmt); fprintf(fp, "%s\n", ssVals.str().c_str()); } @@ -421,8 +431,6 @@ namespace Nyxus std::vector L{ env.uniqueLabels.begin(), env.uniqueLabels.end() }; std::sort(L.begin(), L.end()); - static bool mustRenderHeader = true; // In 'singlecsv' scenario this flag flips from 'T' to 'F' when necessary (after the header is rendered) - // Make the file name and write mode std::string fullPath = get_feature_output_fname (env, intFpath, segFpath); VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\t--> " << fullPath << "\n"); @@ -430,7 +438,7 @@ namespace Nyxus // Single CSV: create or continue? const char* mode = "w"; if (!env.separateCsv) - mode = mustRenderHeader ? "w" : "a"; + mode = mustRenderHeader_segmented ? "w" : "a"; // Open it FILE* fp = nullptr; @@ -453,7 +461,7 @@ namespace Nyxus std::vector> F = env.theFeatureSet.getEnabledFeatures(); // -- Header - if (mustRenderHeader) + if (mustRenderHeader_segmented) { std::stringstream ssHead; @@ -473,7 +481,7 @@ namespace Nyxus // Prevent rendering the header again for another image's portion of labels if (env.separateCsv == false) - mustRenderHeader = false; + mustRenderHeader_segmented = false; } if (need_aggregation) @@ -560,157 +568,7 @@ namespace Nyxus // ROI label and time ssVals << "," << l << "," << t_index; - for (auto& enabdF : F) - { - auto fc = std::get<1>(enabdF); - auto fn = std::get<0>(enabdF); // debug - auto vv = r.get_fvals(fc); - - // Parameterized feature - // --GLCM family - bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); - bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glcmFeature && nonAngledGlcmFeature == false) - { - // Mock angled values if they haven't been calculated for some error reason - if (vv.size() < GLCMFeature::angles.size()) - vv.resize(GLCMFeature::angles.size(), 0.0); - // Output the sub-values - int nAng = (int) GLCMFeature::angles.size(); - for (int i = 0; i < nAng; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - #ifndef DIAGNOSE_NYXUS_OUTPUT - ssVals << "," << rvbuf; - #else - //--diagnoze misalignment-- - ssVals << "," << fn << "-" << rvbuf; - #endif - } - // Proceed with other features - continue; - } - - // --GLRLM family - bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); - bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glrlmFeature && nonAngledGlrlmFeature == false) - { - // Populate with angles - int nAng = 4; - for (int i = 0; i < nAng; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - #ifndef DIAGNOSE_NYXUS_OUTPUT - ssVals << "," << rvbuf; - #else - //--diagnoze misalignment-- - ssVals << "," << fn << "-" << rvbuf; - #endif - } - // Proceed with other features - continue; - } - - // --Gabor - if (fc == (int) Feature2D::GABOR) - { - for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - #ifndef DIAGNOSE_NYXUS_OUTPUT - ssVals << "," << rvbuf; - #else - //--diagnoze misalignment-- - ssVals << "," << fn << "-" << rvbuf; - #endif - } - - // Proceed with other features - continue; - } - - // --Zernike feature values - if (fc == (int) Feature2D::ZERNIKE2D) - { - for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - #ifndef DIAGNOSE_NYXUS_OUTPUT - ssVals << "," << rvbuf; - #else - //--diagnoze misalignment-- - ssVals << "," << fn << "-" << rvbuf; - #endif - } - - // Proceed with other features - continue; - } - - // --Radial distribution features - if (fc == (int) Feature2D::FRAC_AT_D) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - #ifndef DIAGNOSE_NYXUS_OUTPUT - ssVals << "," << rvbuf; - #else - //--diagnoze misalignment-- - ssVals << "," << fn << "-" << rvbuf; - #endif - } - // Proceed with other features - continue; - } - if (fc == (int) Feature2D::MEAN_FRAC) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - #ifndef DIAGNOSE_NYXUS_OUTPUT - ssVals << "," << rvbuf; - #else - //--diagnoze misalignment-- - ssVals << "," << fn << "-" << rvbuf; - #endif - } - // Proceed with other features - continue; - } - if (fc == (int) Feature2D::RADIAL_CV) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) - { - double fv = Nyxus::force_finite_number(vv[i], env.resultOptions.noval()); // safe feature value (no NAN, no inf) - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, fv); - #ifndef DIAGNOSE_NYXUS_OUTPUT - ssVals << "," << rvbuf; - #else - //--diagnoze misalignment-- - ssVals << "," << fn << "-" << rvbuf; - #endif - } - // Proceed with other features - continue; - } - - // Regular feature - snprintf(rvbuf, VAL_BUF_LEN, rvfmt, Nyxus::force_finite_number(vv[0], env.resultOptions.noval())); - #ifndef DIAGNOSE_NYXUS_OUTPUT - ssVals << "," << rvbuf; // Alternatively: auto_precision(ssVals, vv[0]); - #else - //--diagnoze misalignment-- - ssVals << "," << fn << "-" << rvbuf; - #endif - } + write_feature_values_csv(ssVals, r, F, env, rvbuf, VAL_BUF_LEN, rvfmt); fprintf(fp, "%s\n", ssVals.str().c_str()); } @@ -730,112 +588,13 @@ namespace Nyxus { std::vector features; - // user's feature selection std::vector> F = fset.getEnabledFeatures(); + auto fvals = collect_feature_values(r, F); - // numeric columns - std::vector fvals; - - for (auto& enabdF : F) - { - auto fc = std::get<1>(enabdF); - auto vv = r.get_fvals(std::get<1>(enabdF)); - - // Parameterized feature - // --GLCM family - bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); - bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glcmFeature && nonAngledGlcmFeature == false) - { - // Mock angled values if they haven't been calculated for some error reason - if (vv.size() < GLCMFeature::angles.size()) - vv.resize(GLCMFeature::angles.size(), 0.0); - // Output the sub-values - int nAng = (int)GLCMFeature::angles.size(); - for (int i = 0; i < nAng; i++) - { - fvals.push_back(vv[i]); - } - // Proceed with other features - continue; - } - - // --GLRLM family - bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); - bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glrlmFeature && nonAngledGlrlmFeature == false) - { - // Polulate with angles - int nAng = 4; - for (int i = 0; i < nAng; i++) - { - fvals.push_back(vv[i]); - } - // Proceed with other features - continue; - } - - // --Gabor - if (fc == (int)Feature2D::GABOR) - { - for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) - { - fvals.push_back(vv[i]); - } - - // Proceed with other features - continue; - } - - // --Zernike feature values - if (fc == (int)Feature2D::ZERNIKE2D) - { - for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) - { - fvals.push_back(vv[i]); - } - - // Proceed with other features - continue; - } - - // --Radial distribution features - if (fc == (int)Feature2D::FRAC_AT_D) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) - { - fvals.push_back(vv[i]); - } - // Proceed with other features - continue; - } - if (fc == (int)Feature2D::MEAN_FRAC) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) - { - fvals.push_back(vv[i]); - } - // Proceed with other features - continue; - } - if (fc == (int)Feature2D::RADIAL_CV) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) - { - fvals.push_back(vv[i]); - } - // Proceed with other features - continue; - } - - fvals.push_back(vv[0]); - } - - // other columns std::vector textcols; textcols.push_back ((fs::path(ifpath)).filename().string()); textcols.push_back (""); - int roilabl = r.label; // whole-slide roi # + int roilabl = r.label; features.push_back (std::make_tuple(textcols, roilabl, -999.88, fvals)); @@ -862,114 +621,19 @@ namespace Nyxus { const LR& r = roiData.at (l); - std::vector feature_values; - // Skip blacklisted ROI if (r.blacklisted) continue; // Tear off pure file names from segment and intensity file paths const SlideProps& sli = dataset.dataset_props [r.slide_idx]; - fs::path pseg (sli.fname_seg), + fs::path pseg (sli.fname_seg), pint (sli.fname_int); std::vector filenames; filenames.push_back(pint.filename().string()); filenames.push_back(pseg.filename().string()); - for (auto& enabdF : F) - { - auto fc = std::get<1>(enabdF); - auto vv = r.get_fvals(std::get<1>(enabdF)); - - // Parameterized feature - // --GLCM family - bool glcmFeature = std::find(GLCMFeature::featureset.begin(), GLCMFeature::featureset.end(), (Feature2D)fc) != GLCMFeature::featureset.end(); - bool nonAngledGlcmFeature = std::find(GLCMFeature::nonAngledFeatures.begin(), GLCMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLCMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glcmFeature && nonAngledGlcmFeature == false) - { - // Mock angled values if they haven't been calculated for some error reason - if (vv.size() < GLCMFeature::angles.size()) - vv.resize(GLCMFeature::angles.size(), 0.0); - // Output the sub-values - int nAng = (int) GLCMFeature::angles.size(); - for (int i = 0; i < nAng; i++) - { - feature_values.push_back(vv[i]); - } - // Proceed with other features - continue; - } - - // --GLRLM family - bool glrlmFeature = std::find(GLRLMFeature::featureset.begin(), GLRLMFeature::featureset.end(), (Feature2D)fc) != GLRLMFeature::featureset.end(); - bool nonAngledGlrlmFeature = std::find(GLRLMFeature::nonAngledFeatures.begin(), GLRLMFeature::nonAngledFeatures.end(), (Feature2D)fc) != GLRLMFeature::nonAngledFeatures.end(); // prevent output of a non-angled feature in an angled way - if (glrlmFeature && nonAngledGlrlmFeature == false) - { - // Polulate with angles - int nAng = 4; - for (int i = 0; i < nAng; i++) - { - feature_values.push_back(vv[i]); - } - // Proceed with other features - continue; - } - - // --Gabor - if (fc == (int) Feature2D::GABOR) - { - for (auto i = 0; i < GaborFeature::f0_theta_pairs.size(); i++) - { - feature_values.push_back(vv[i]); - } - - // Proceed with other features - continue; - } - - // --Zernike feature values - if (fc == (int) Feature2D::ZERNIKE2D) - { - for (int i = 0; i < ZernikeFeature::NUM_FEATURE_VALS; i++) - { - feature_values.push_back(vv[i]); - } - - // Proceed with other features - continue; - } - - // --Radial distribution features - if (fc == (int) Feature2D::FRAC_AT_D) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_FracAtD; i++) - { - feature_values.push_back(vv[i]); - } - // Proceed with other features - continue; - } - if (fc == (int) Feature2D::MEAN_FRAC) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_MeanFrac; i++) - { - feature_values.push_back(vv[i]); - } - // Proceed with other features - continue; - } - if (fc == (int) Feature2D::RADIAL_CV) - { - for (auto i = 0; i < RadialDistributionFeature::num_features_RadialCV; i++) - { - feature_values.push_back(vv[i]); - } - // Proceed with other features - continue; - } - - feature_values.push_back(vv[0]); - } + auto feature_values = collect_feature_values(r, F); features.push_back (std::make_tuple(filenames, l, DEFAULT_T_INDEX, feature_values)); } @@ -1017,4 +681,148 @@ namespace Nyxus } + bool save_features_2_csv_fmaps ( + Environment & env, + const std::string & intFpath, + const std::string & segFpath, + const std::string & outputDir, + int slide_idx, + const std::unordered_set & childLabels, + const std::unordered_map & childRoiData, + const std::unordered_map & childToParentMap) + { + std::lock_guard lg (mutex1); + + constexpr int VAL_BUF_LEN = 450; + char rvbuf[VAL_BUF_LEN]; + const char rvfmt[] = "%g"; + + // Sort child labels + std::vector L(childLabels.begin(), childLabels.end()); + std::sort(L.begin(), L.end()); + + if (L.empty()) + return true; + + // Output file — fmaps results are written to a separate CSV with "_fmaps.csv" suffix + // to keep them distinct from standard per-ROI feature results + std::string fullPath; + if (env.separateCsv) + fullPath = outputDir + "/_INT_" + getPureFname(intFpath) + "_SEG_" + getPureFname(segFpath) + "_fmaps.csv"; + else + fullPath = outputDir + "/" + env.nyxus_result_fname + "_fmaps.csv"; + + const char* mode = "w"; + if (!env.separateCsv) + mode = mustRenderHeader_fmaps ? "w" : "a"; + + FILE* fp = nullptr; + auto eno = fopen_s(&fp, fullPath.c_str(), mode); + if (!fp) + { + std::cerr << "Cannot open file " << fullPath << " for writing. Errno=" << eno << "\n"; + return false; + } + + if (std::setvbuf(fp, nullptr, _IOFBF, 32768) != 0) + std::cerr << "setvbuf failed \n"; + + std::vector> F = env.theFeatureSet.getEnabledFeatures(); + + // Header + if (mustRenderHeader_fmaps) + { + std::stringstream ssHead; + + auto head_vector = Nyxus::get_header(env); + + // Build fmaps header: insert fmaps columns after mask_image + annotations, before ROI_label + // get_header returns: [intensity_image, mask_image, , ROI_label, t_index, ] + // We want: [intensity_image, mask_image, , parent_roi_label, kernel_center_x, kernel_center_y, ROI_label, t_index, ] + std::vector fmaps_head; + fmaps_head.push_back(head_vector[0]); // intensity_image + fmaps_head.push_back(head_vector[1]); // mask_image + + // Annotation columns (if any) + int n_annot_cols = 0; + if (env.resultOptions.need_annotation()) + { + auto slp = env.dataset.dataset_props[0]; + n_annot_cols = (int) slp.annots.size(); + } + int roi_label_pos = 2 + n_annot_cols; // position of ROI_label in head_vector + for (int i = 2; i < roi_label_pos; i++) + fmaps_head.push_back(head_vector[i]); + + // Fmaps-specific columns + fmaps_head.push_back("parent_roi_label"); + fmaps_head.push_back("kernel_center_x"); + fmaps_head.push_back("kernel_center_y"); + + // Remaining standard columns (ROI_label, t_index, features...) + for (size_t i = roi_label_pos; i < head_vector.size(); i++) + fmaps_head.push_back(head_vector[i]); + + for (int i = 0; i < (int)fmaps_head.size(); i++) + { + if (i) + ssHead << ','; + ssHead << '\"' << fmaps_head[i] << '\"'; + } + + fprintf(fp, "%s\n", ssHead.str().c_str()); + + if (env.separateCsv == false) + mustRenderHeader_fmaps = false; + } + + // Values + for (auto l : L) + { + const LR& r = childRoiData.at(l); + + if (r.blacklisted) + continue; + + auto it = childToParentMap.find(l); + if (it == childToParentMap.end()) + continue; + + const FmapChildInfo& info = it->second; + + std::stringstream ssVals; + ssVals << std::fixed; + + // File names + fs::path pint(intFpath), pseg(segFpath); + ssVals << pint.filename().string() << "," << pseg.filename().string(); + + // Annotation columns + if (env.resultOptions.need_annotation()) + { + const SlideProps& sli = env.dataset.dataset_props[slide_idx]; + for (const auto& a : sli.annots) + ssVals << "," << a; + } + + // Feature maps columns + ssVals << "," << info.parent_label; + ssVals << "," << info.center_x; + ssVals << "," << info.center_y; + + // ROI label (child) and time + ssVals << "," << l << "," << DEFAULT_T_INDEX; + + // Feature values + write_feature_values_csv(ssVals, r, F, env, rvbuf, VAL_BUF_LEN, rvfmt); + + fprintf(fp, "%s\n", ssVals.str().c_str()); + } + + std::fflush(fp); + std::fclose(fp); + + return true; + } + } diff --git a/src/nyx/parallel.h b/src/nyx/parallel.h index 80dc73f3..5662c192 100644 --- a/src/nyx/parallel.h +++ b/src/nyx/parallel.h @@ -39,6 +39,11 @@ namespace Nyxus idxE = datasetSize; // include the tail T.push_back(std::async(std::launch::async, f, idxS, idxE, ptrLabels, ptrLabelData, std::cref(f_settings), std::cref(dataset))); } + // Wait for all threads to complete before returning. + // Critical for fmaps mode where EnvRoiSwapGuard will move + // env.roiData on scope exit — threads must finish first. + for (auto& fut : T) + fut.get(); } void parallelReduceConvHull (size_t start, size_t end, std::vector* ptrLabels, std::unordered_map * ptrLabelData, const Fsettings & fst, const Dataset & ds); diff --git a/src/nyx/python/new_bindings_py.cpp b/src/nyx/python/new_bindings_py.cpp index 8a377ba3..f4f13136 100644 --- a/src/nyx/python/new_bindings_py.cpp +++ b/src/nyx/python/new_bindings_py.cpp @@ -36,7 +36,7 @@ namespace Nyxus { }; -using ParameterTypes = std::variant, std::vector>; +using ParameterTypes = std::variant, std::vector>; // Defined in nested.cpp bool mine_segment_relations ( @@ -52,13 +52,60 @@ bool mine_segment_relations ( template inline py::array_t as_pyarray(Sequence &&seq) { - auto size = seq.size(); - auto data = seq.data(); - std::unique_ptr seq_ptr = std::make_unique(std::move(seq)); - auto capsule = py::capsule(seq_ptr.get(), [](void *p) - { std::unique_ptr(reinterpret_cast(p)); }); - seq_ptr.release(); - return py::array(size, data, capsule); + // Copy data into a numpy-owned array rather than wrapping the C++ memory + // with a pybind11 capsule. A capsule's destructor is compiled code inside + // this .so — if Python's GC frees the array after the .so is unloaded + // (e.g. during interpreter shutdown), the destructor call jumps to + // unmapped memory and segfaults. Copying avoids that by letting numpy + // manage the memory entirely with no reference back to this module. + using T = typename Sequence::value_type; + py::array_t arr(seq.size()); + std::copy(seq.begin(), seq.end(), arr.mutable_data()); + return arr; +} + +/// @brief Converts accumulated FmapArrayResult entries into a Python list of dicts. +/// Each dict contains: parent_roi_label, intensity_image, mask_image, origin_x, origin_y, +/// and a 'features' dict mapping feature names to numpy arrays. +/// For 2D: arrays are (map_h, map_w). For 3D: arrays are (map_d, map_h, map_w). +py::list fmap_results_to_python(ResultsCache & rescache) +{ + py::list result; + for (auto & fr : rescache.get_fmapArrayResults()) + { + py::dict d; + d["parent_roi_label"] = fr.parent_label; + d["intensity_image"] = fr.intens_name; + d["mask_image"] = fr.seg_name; + d["origin_x"] = fr.origin_x; + d["origin_y"] = fr.origin_y; + + bool is_3d = fr.map_d > 1; + if (is_3d) + d["origin_z"] = fr.origin_z; + + size_t map_size = (size_t)fr.map_d * fr.map_h * fr.map_w; + size_t n_features = fr.feature_names.size(); + + py::dict features; + for (size_t fi = 0; fi < n_features; fi++) + { + py::array_t arr; + if (is_3d) + arr = py::array_t({fr.map_d, fr.map_h, fr.map_w}); + else + arr = py::array_t({fr.map_h, fr.map_w}); + auto ptr = arr.mutable_data(); + std::copy( + fr.feature_data.begin() + fi * map_size, + fr.feature_data.begin() + (fi + 1) * map_size, + ptr); + features[py::str(fr.feature_names[fi])] = arr; + } + d["features"] = features; + result.append(d); + } + return result; } void initialize_environment( @@ -67,7 +114,7 @@ void initialize_environment( const std::vector &features, int neighbor_distance, float pixels_per_micron, - uint32_t coarse_gray_depth, + int coarse_gray_depth, uint32_t n_reduce_threads, int using_gpu, bool ibsi, @@ -79,7 +126,9 @@ void initialize_environment( int verb_lvl, float aniso_x, float aniso_y, - float aniso_z) + float aniso_z, + bool fmaps = false, + int fmaps_radius = 2) { Environment & theEnvironment = Nyxus::findenv (instid); @@ -93,6 +142,11 @@ void initialize_environment( theEnvironment.n_reduce_threads = n_reduce_threads; theEnvironment.ibsi_compliance = ibsi; + // feature maps + theEnvironment.fmaps_mode = fmaps; + if (fmaps_radius >= 1) + theEnvironment.fmaps_kernel_radius = fmaps_radius; + // Throws exception if invalid feature is passed theEnvironment.expand_featuregroups(); if (!theEnvironment.theFeatureMgr.compile()) @@ -108,7 +162,7 @@ void initialize_environment( theEnvironment.fpimageOptions.set_max_intensity(max_intensity); // RAM limit controlling trivial-nontrivial featurization - if (ram_limit_mb >= 0) + if (ram_limit_mb >= 0) theEnvironment.set_ram_limit(ram_limit_mb); // anisotropy @@ -118,17 +172,17 @@ void initialize_environment( // GPU related #ifdef USE_GPU - if(using_gpu == -1) + if(using_gpu == -1) { theEnvironment.set_using_gpu(false); - } - else + } + else { theEnvironment.set_gpu_device_id(using_gpu); theEnvironment.set_using_gpu(true); } - #else - if (using_gpu != -1) + #else + if (using_gpu != -1) { throw std::runtime_error ("this Nyxus backend was built without the GPU support"); } @@ -142,12 +196,29 @@ void set_if_ibsi_imp (uint64_t instid, bool ibsi) env.set_ibsi_compliance (ibsi); } +void set_fmaps_imp (uint64_t instid, int set_mode, int radius) +{ + Environment & env = Nyxus::findenv (instid); + // Always validate radius when explicitly provided (not sentinel -1), + // even if fmaps_mode is false, to prevent stale invalid values + if (radius != -1) + { + if (radius >= 1) + env.fmaps_kernel_radius = radius; + else + throw std::invalid_argument("fmaps_radius must be an integer >= 1"); + } + // set_mode: 0=disable, 1=enable, -1=leave unchanged (radius-only update) + if (set_mode >= 0) + env.fmaps_mode = (set_mode != 0); +} + void set_environment_params_imp ( uint64_t instid, const std::vector &features = {}, int neighbor_distance = -1, float pixels_per_micron = -1, - uint32_t coarse_gray_depth = 0, + int coarse_gray_depth = 0, uint32_t n_reduce_threads = 0, int using_gpu = -2, float dynamic_range = -1, @@ -267,7 +338,18 @@ py::tuple featurize_directory_imp( // run the workflow int ercode = 0; - if (env.singleROI) + if (env.fmaps_prevents_arrow()) + throw std::runtime_error("Arrow/Parquet output is not supported in feature maps (fmaps) mode. Use CSV or buffer output instead."); + + if (env.fmaps_mode) + ercode = processDataset_2D_fmaps( + env, + intensFiles, + labelFiles, + env.n_reduce_threads, + env.saveOption, + output_path); + else if (env.singleROI) ercode = processDataset_2D_wholeslide( env, intensFiles, @@ -287,10 +369,15 @@ py::tuple featurize_directory_imp( if (ercode) throw std::runtime_error("Error: " + std::to_string(ercode)); - // shut down the output structures + // Fmaps mode returns spatial arrays, not a DataFrame + if (env.fmaps_mode && env.saveOption == Nyxus::SaveOption::saveBuffer) + { + auto fmaps = fmap_results_to_python(env.theResultsCache); + return py::make_tuple(fmaps); + } // shape the resulting data frame - + if (env.saveOption == Nyxus::SaveOption::saveBuffer) { // has the backend produced any result ? @@ -511,19 +598,39 @@ py::tuple featurize_directory_3D_imp( if (mayBerror.has_value()) throw std::runtime_error ("Error traversing dataset: " + *mayBerror); - int errorCode = processDataset_3D_segmented( - theEnvironment, - intensFiles, - labelFiles, - theEnvironment.n_reduce_threads, - theEnvironment.saveOption, - output_path); + int errorCode = 0; + if (theEnvironment.fmaps_mode) + { + errorCode = processDataset_3D_fmaps( + theEnvironment, + intensFiles, + labelFiles, + theEnvironment.n_reduce_threads, + theEnvironment.saveOption, + output_path); + } + else + { + errorCode = processDataset_3D_segmented( + theEnvironment, + intensFiles, + labelFiles, + theEnvironment.n_reduce_threads, + theEnvironment.saveOption, + output_path); + } if (errorCode) throw std::runtime_error ("Error " + std::to_string(errorCode) + " occurred during dataset processing"); } // Output the result + if (theEnvironment.fmaps_mode && theEnvironment.saveOption == Nyxus::SaveOption::saveBuffer) + { + py::list fmaps = fmap_results_to_python(theEnvironment.theResultsCache); + return py::make_tuple(fmaps); + } + if (theEnvironment.saveOption == Nyxus::SaveOption::saveBuffer) { auto pyHeader = py::array(py::cast(theEnvironment.theResultsCache.get_headerBuf())); @@ -608,6 +715,12 @@ py::tuple featurize_montage_imp ( if (mayBerror.has_value()) throw std::runtime_error ("Error occurred during dataset processing: " + *mayBerror); + if (theEnvironment.fmaps_mode) + { + auto fmaps = fmap_results_to_python(theEnvironment.theResultsCache); + return py::make_tuple(fmaps, error_message); + } + if (theEnvironment.saveOption == Nyxus::SaveOption::saveBuffer) { auto pyHeader = py::array(py::cast(theEnvironment.theResultsCache.get_headerBuf())); @@ -619,7 +732,7 @@ py::tuple featurize_montage_imp ( pyNumData = pyNumData.reshape({ nRows, pyNumData.size() / nRows }); return py::make_tuple(pyHeader, pyStrData, pyNumData, error_message); - } + } return py::make_tuple(error_message); } @@ -681,17 +794,35 @@ py::tuple featurize_fname_lists_imp (uint64_t instid, const py::list& int_fnames } else {return SaveOption::saveBuffer;} }(); - errorCode = processDataset_2D_segmented ( - theEnvironment, - intensFiles, - labelFiles, - theEnvironment.n_reduce_threads, - theEnvironment.saveOption, - output_path); + if (theEnvironment.fmaps_prevents_arrow()) + throw std::runtime_error("Arrow/Parquet output is not supported in feature maps (fmaps) mode. Use CSV or buffer output instead."); + + if (theEnvironment.fmaps_mode) + errorCode = processDataset_2D_fmaps ( + theEnvironment, + intensFiles, + labelFiles, + theEnvironment.n_reduce_threads, + theEnvironment.saveOption, + output_path); + else + errorCode = processDataset_2D_segmented ( + theEnvironment, + intensFiles, + labelFiles, + theEnvironment.n_reduce_threads, + theEnvironment.saveOption, + output_path); if (errorCode) throw std::runtime_error("Error occurred during dataset processing."); + if (theEnvironment.fmaps_mode && theEnvironment.saveOption == Nyxus::SaveOption::saveBuffer) + { + auto fmaps = fmap_results_to_python(theEnvironment.theResultsCache); + return py::make_tuple(fmaps); + } + if (theEnvironment.saveOption == Nyxus::SaveOption::saveBuffer) { auto pyHeader = py::array(py::cast(theEnvironment.theResultsCache.get_headerBuf())); @@ -703,7 +834,7 @@ py::tuple featurize_fname_lists_imp (uint64_t instid, const py::list& int_fnames pyNumData = pyNumData.reshape({ nRows, pyNumData.size() / nRows }); return py::make_tuple(pyHeader, pyStrData, pyNumData); - } + } // Return "nothing" when output will be an Arrow format return py::make_tuple(); @@ -784,18 +915,37 @@ py::tuple featurize_fname_lists_3D_imp ( } }(); - errorCode = processDataset_3D_segmented( - theEnvironment, - ifiles, - mfiles, - theEnvironment.n_reduce_threads, - theEnvironment.saveOption, - output_path); + if (theEnvironment.fmaps_mode) + { + errorCode = processDataset_3D_fmaps( + theEnvironment, + ifiles, + mfiles, + theEnvironment.n_reduce_threads, + theEnvironment.saveOption, + output_path); + } + else + { + errorCode = processDataset_3D_segmented( + theEnvironment, + ifiles, + mfiles, + theEnvironment.n_reduce_threads, + theEnvironment.saveOption, + output_path); + } if (errorCode) throw std::runtime_error("Error occurred during dataset processing"); // save the result + if (theEnvironment.fmaps_mode && theEnvironment.saveOption == Nyxus::SaveOption::saveBuffer) + { + py::list fmaps = fmap_results_to_python(theEnvironment.theResultsCache); + return py::make_tuple(fmaps); + } + if (theEnvironment.saveOption == Nyxus::SaveOption::saveBuffer) { auto pyHeader = py::array(py::cast(theEnvironment.theResultsCache.get_headerBuf())); @@ -946,7 +1096,9 @@ std::map get_params_imp (uint64_t instid, const std params["features"] = theEnvironment.recognizedFeatureNames; params["neighbor_distance"] = theEnvironment.n_pixel_distance; params["pixels_per_micron"] = theEnvironment.xyRes; - params["coarse_gray_depth"] = theEnvironment.get_coarse_gray_depth(); + int cgd = theEnvironment.get_coarse_gray_depth(); + params["coarse_gray_depth"] = std::abs(cgd); + params["binning_origin"] = std::string(cgd < 0 ? "min" : "zero"); params["n_feature_calc_threads"] = theEnvironment.n_reduce_threads; params["ibsi"] = theEnvironment.ibsi_compliance; @@ -970,7 +1122,10 @@ std::map get_params_imp (uint64_t instid, const std params["max_intensity"] = theEnvironment.fpimageOptions.max_intensity(); params["ram_limit"] = (int)(theEnvironment.get_ram_limit()/1048576); // convert from bytes to megabytes - if (vars.size() == 0) + params["fmaps"] = theEnvironment.fmaps_mode; + params["fmaps_radius"] = theEnvironment.fmaps_kernel_radius; + + if (vars.size() == 0) return params; std::map params_subset; @@ -1048,6 +1203,18 @@ py::tuple get_metaparam_imp (uint64_t instid, const std::string p_name) PYBIND11_MODULE(backend, m) { m.doc() = "Nyxus"; + + // Register an atexit handler to destroy all Environment objects while + // the Python interpreter is still alive. Without this, the global + // pynyxus_cache is destroyed during C++ static destruction — after + // Python has already torn down modules — causing segfaults from + // stale references in LR/ImageMatrix/ResultsCache destructors. + auto atexit = py::module_::import("atexit"); + atexit.attr("register")(py::cpp_function([]() { + Nyxus::pynyxus_cache.clear(); + Nyxus::unique_pynyxus_ids.clear(); + })); + m.def("initialize_environment", &initialize_environment, "Environment initialization"); m.def("featurize_directory_imp", &featurize_directory_imp, "Calculate features of images defined by intensity and mask image collection directories"); m.def("featurize_directory_3D_imp", &featurize_directory_3D_imp, "Calculate 3D features of images defined by intensity and mask image collection directories"); @@ -1063,6 +1230,7 @@ PYBIND11_MODULE(backend, m) m.def("roi_blacklist_get_summary_imp", &roi_blacklist_get_summary_imp, "Returns a summary of the ROI blacklist"); m.def("customize_gabor_feature_imp", &customize_gabor_feature_imp, "Sets custom GABOR feature's parameters"); m.def("set_if_ibsi_imp", &set_if_ibsi_imp, "Set if the features will be ibsi compliant"); + m.def("set_fmaps_imp", &set_fmaps_imp, "Enable/disable feature maps mode and set kernel radius"); m.def("set_environment_params_imp", &set_environment_params_imp, "Set the environment variables of Nyxus"); m.def("get_params_imp", &get_params_imp, "Get parameters of Nyxus"); m.def("arrow_is_enabled_imp", &arrow_is_enabled_imp, "Check if arrow is enabled."); diff --git a/src/nyx/python/nyxus/__init__.py b/src/nyx/python/nyxus/__init__.py index 5dc5a8d1..22f44e59 100644 --- a/src/nyx/python/nyxus/__init__.py +++ b/src/nyx/python/nyxus/__init__.py @@ -3,6 +3,7 @@ from .nyxus import Nested from .nyxus import ImageQuality from .functions import gpu_is_available, get_gpu_properties +from .fmap_io import save_fmaps_to_tiff, save_fmaps_to_nifti from . import _version __version__ = _version.get_versions()['version'] diff --git a/src/nyx/python/nyxus/fmap_io.py b/src/nyx/python/nyxus/fmap_io.py new file mode 100644 index 00000000..2d1f22f6 --- /dev/null +++ b/src/nyx/python/nyxus/fmap_io.py @@ -0,0 +1,174 @@ +"""Utilities for saving feature-map (fmaps) results to image formats. + +Feature-map results from Nyxus/Nyxus3D are returned as a list of dicts, +one per parent ROI, each containing numpy arrays of spatial feature maps. +These utilities write those arrays to TIFF stacks or NIfTI volumes so +they can be consumed by downstream imaging and ML pipelines. +""" + +import os +import numpy as np +from typing import List, Dict, Optional + + +def save_fmaps_to_tiff( + fmaps: List[Dict], + output_dir: str, + prefix: str = "fmap", + features: Optional[List[str]] = None, +): + """Save 2D feature maps to multi-page TIFF files. + + For each parent ROI, creates one TIFF per feature. If the feature map + is 3D (D, H, W), each z-slice becomes a page in the TIFF stack. + + Requires ``tifffile`` (``pip install tifffile``). + + Parameters + ---------- + fmaps : list[dict] + Feature-map results from ``Nyxus.featurize_directory()`` or similar, + obtained with ``fmaps=True``. + output_dir : str + Directory to write TIFF files into. Created if it does not exist. + prefix : str, optional + Filename prefix. Default ``"fmap"``. + features : list[str], optional + Subset of feature names to save. Default saves all features. + + Returns + ------- + list[str] + Paths of all TIFF files written. + """ + try: + import tifffile + except ImportError: + raise ImportError( + "tifffile is required for TIFF export. Install it with: pip install tifffile" + ) + + os.makedirs(output_dir, exist_ok=True) + written = [] + + for roi_result in fmaps: + roi_label = roi_result["parent_roi_label"] + feat_dict = roi_result["features"] + + names = features if features else list(feat_dict.keys()) + + for feat_name in names: + if feat_name not in feat_dict: + continue + arr = feat_dict[feat_name] + + # Ensure float32 for compatibility + arr = np.asarray(arr, dtype=np.float32) + + fname = f"{prefix}_roi{roi_label}_{feat_name}.tif" + fpath = os.path.join(output_dir, fname) + + if arr.ndim == 2: + tifffile.imwrite(fpath, arr) + elif arr.ndim == 3: + # (D, H, W) -> multi-page TIFF + tifffile.imwrite(fpath, arr) + else: + raise ValueError( + f"Unexpected array dimensions {arr.ndim} for feature '{feat_name}'" + ) + + written.append(fpath) + + return written + + +def save_fmaps_to_nifti( + fmaps: List[Dict], + output_dir: str, + prefix: str = "fmap", + features: Optional[List[str]] = None, + voxel_size: tuple = (1.0, 1.0, 1.0), +): + """Save 3D feature maps to NIfTI-1 (.nii.gz) volumes. + + For each parent ROI, creates one NIfTI file per feature. + The affine is set so the origin matches the ROI's global coordinates. + + Requires ``nibabel`` (``pip install nibabel``). + + Parameters + ---------- + fmaps : list[dict] + Feature-map results from ``Nyxus3D.featurize_directory()`` or similar, + obtained with ``fmaps=True``. + output_dir : str + Directory to write NIfTI files into. Created if it does not exist. + prefix : str, optional + Filename prefix. Default ``"fmap"``. + features : list[str], optional + Subset of feature names to save. Default saves all features. + voxel_size : tuple of float, optional + Voxel dimensions (x, y, z) in physical units. Default ``(1.0, 1.0, 1.0)``. + + Returns + ------- + list[str] + Paths of all NIfTI files written. + """ + try: + import nibabel as nib + except ImportError: + raise ImportError( + "nibabel is required for NIfTI export. Install it with: pip install nibabel" + ) + + os.makedirs(output_dir, exist_ok=True) + written = [] + + for roi_result in fmaps: + roi_label = roi_result["parent_roi_label"] + feat_dict = roi_result["features"] + origin_x = roi_result.get("origin_x", 0) + origin_y = roi_result.get("origin_y", 0) + origin_z = roi_result.get("origin_z", 0) + + names = features if features else list(feat_dict.keys()) + + for feat_name in names: + if feat_name not in feat_dict: + continue + arr = feat_dict[feat_name] + arr = np.asarray(arr, dtype=np.float32) + + # For 2D maps, add a singleton z-dimension + if arr.ndim == 2: + arr = arr[np.newaxis, :, :] + + if arr.ndim != 3: + raise ValueError( + f"Unexpected array dimensions {arr.ndim} for feature '{feat_name}'" + ) + + # NIfTI expects (x, y, z) axis order; our arrays are (D, H, W) + # Transpose to (W, H, D) = (x, y, z) for NIfTI convention + arr_nifti = np.transpose(arr, (2, 1, 0)) + + # Build affine: diagonal scaling + translation to ROI origin + affine = np.eye(4) + affine[0, 0] = voxel_size[0] + affine[1, 1] = voxel_size[1] + affine[2, 2] = voxel_size[2] + affine[0, 3] = origin_x * voxel_size[0] + affine[1, 3] = origin_y * voxel_size[1] + affine[2, 3] = origin_z * voxel_size[2] + + img = nib.Nifti1Image(arr_nifti, affine) + + fname = f"{prefix}_roi{roi_label}_{feat_name}.nii.gz" + fpath = os.path.join(output_dir, fname) + nib.save(img, fpath) + + written.append(fpath) + + return written diff --git a/src/nyx/python/nyxus/nyxus.py b/src/nyx/python/nyxus/nyxus.py index a2725e07..94e09d51 100644 --- a/src/nyx/python/nyxus/nyxus.py +++ b/src/nyx/python/nyxus/nyxus.py @@ -13,10 +13,11 @@ roi_blacklist_get_summary_imp, customize_gabor_feature_imp, set_if_ibsi_imp, + set_fmaps_imp, set_environment_params_imp, get_params_imp, arrow_is_enabled_imp, - get_arrow_file_imp, + get_arrow_file_imp, get_parquet_file_imp, set_metaparam_imp, get_metaparam_imp) @@ -120,6 +121,21 @@ class Nyxus: X-dimension scale factor anisotropy_y: float (optional, default 1.0) Y-dimension scale factor + fmaps: bool (optional, default False) + Enable feature maps mode. Instead of computing a single feature vector per ROI, + a sliding kernel is moved across each ROI and features are computed at every + position, producing spatial feature maps as numpy arrays. + When enabled, featurize methods return a list of dicts (one per parent ROI) + instead of a DataFrame. Each dict contains numpy arrays shaped (map_h, map_w). + fmaps_radius: int (optional, default 2) + Radius of the sliding kernel used in feature maps mode. + The kernel size is (2 * fmaps_radius + 1). For example, fmaps_radius=2 + produces a 5x5 kernel. Must be >= 1. + binning_origin: str (optional, default "zero") + Origin of the intensity binning range for texture features. + "zero" - bins span [0, max] (default Nyxus/MATLAB behavior). + "min" - bins span [min, max], adapting to the actual data range + (PyRadiomics-compatible behavior). """ def __init__( @@ -131,10 +147,12 @@ def __init__( 'neighbor_distance', 'pixels_per_micron', 'coarse_gray_depth', 'n_feature_calc_threads', 'use_gpu_device', 'ibsi', 'gabor_kersize', 'gabor_gamma', 'gabor_sig2lam', 'gabor_f0', - 'gabor_thold', 'gabor_thetas', 'gabor_freqs', 'channel_signature', + 'gabor_thold', 'gabor_thetas', 'gabor_freqs', 'channel_signature', 'parent_channel', 'child_channel', 'aggregate', 'dynamic_range', 'min_intensity', 'max_intensity', 'ram_limit', 'verbose', - 'anisotropy_x', 'anisotropy_y' + 'anisotropy_x', 'anisotropy_y', + 'fmaps', 'fmaps_radius', + 'binning_origin' } # Check for unexpected keyword arguments @@ -164,15 +182,25 @@ def __init__( verb_lvl = kwargs.get('verbose', 0) aniso_x = kwargs.get('anisotropy_x', 1.0) aniso_y = kwargs.get('anisotropy_y', 1.0) - + fmaps = kwargs.get('fmaps', False) + fmaps_radius = kwargs.get('fmaps_radius', 2) + binning_origin = kwargs.get('binning_origin', 'zero') + if neighbor_distance <= 0: raise ValueError("Neighbor distance must be greater than zero.") if pixels_per_micron <= 0: raise ValueError("Pixels per micron must be greater than zero.") - if coarse_gray_depth <= 0: - raise ValueError("Custom number of grayscale levels (parameter coarse_gray_depth, default=64) must be non-negative.") + if coarse_gray_depth < 1: + raise ValueError("coarse_gray_depth must be >= 1.") + + if binning_origin not in ('zero', 'min'): + raise ValueError("binning_origin must be 'zero' or 'min'.") + + # Negative coarse_gray_depth signals data-adaptive (min-based) binning to the C++ backend + if binning_origin == 'min': + coarse_gray_depth = -coarse_gray_depth if n_feature_calc_threads < 1: raise ValueError("There must be at least one feature calculation thread.") @@ -189,6 +217,11 @@ def __init__( if aniso_y <= 0: raise ValueError ("anisotropy_y must be positive") + if fmaps and fmaps_radius < 1: + raise ValueError("fmaps_radius must be an integer >= 1") + + self._fmaps = fmaps + aniso_z = 1.0 # not used in 2D initialize_environment( @@ -197,7 +230,7 @@ def __init__( features, neighbor_distance, pixels_per_micron, - coarse_gray_depth, + coarse_gray_depth, n_feature_calc_threads, use_gpu_device, ibsi, @@ -209,7 +242,9 @@ def __init__( verb_lvl, aniso_x, aniso_y, - aniso_z) + aniso_z, + fmaps, + fmaps_radius) self.set_gabor_feature_params( kersize = gabor_kersize, @@ -320,9 +355,13 @@ def featurize_directory( if (output_type not in self._valid_output_types): raise ValueError(f'Invalid output type {output_type}. Valid output types are {self._valid_output_types}.') - + + if self._fmaps and output_type == 'pandas': + result = featurize_directory_imp(id(self), intensity_dir, label_dir, file_pattern, output_type, "") + return result[0] # list of dicts with numpy array feature maps + if (output_type == 'pandas'): - + header, string_data, numeric_data = featurize_directory_imp (id(self), intensity_dir, label_dir, file_pattern, output_type, "") df = pd.concat( @@ -338,11 +377,11 @@ def featurize_directory( df.ROI_label = df.ROI_label.astype(np.uint32) return df - + else: - + featurize_directory_imp (id(self), intensity_dir, label_dir, file_pattern, output_type, output_path) - + return get_arrow_file_imp (id(self)) # return path to file @@ -453,6 +492,13 @@ def featurize( M = label_images.astype (np.uint32) # featurize + if self._fmaps and output_type == 'pandas': + fmaps_list, error_message = featurize_montage_imp(id(self), I, M, intensity_names, label_names, output_type, "") + self.error_message = error_message + if error_message != '': + print(error_message) + return fmaps_list # list of dicts with numpy array feature maps + if (output_type == 'pandas'): header, string_data, numeric_data, error_message = featurize_montage_imp (id(self), I, M, intensity_names, label_names, output_type, "") self.error_message = error_message @@ -472,13 +518,13 @@ def featurize( df.ROI_label = df.ROI_label.astype(np.uint32) return df - + else: ret = featurize_montage_imp (id(self), I, M, intensity_names, label_names, output_type, output_path) self.error_message = ret[0] if(self.error_message != ''): raise RuntimeError('Error calculating features: ' + error_message[0]) - + return get_arrow_file_imp (id(self)) # return path to file @@ -531,8 +577,12 @@ def featurize_files ( if (output_type not in self._valid_output_types): raise ValueError(f'Invalid output type {output_type}. Valid output types are {self._valid_output_types}') + if self._fmaps and output_type == 'pandas': + result = featurize_fname_lists_imp(id(self), intensity_files, mask_files, single_roi, output_type, "") + return result[0] # list of dicts with numpy array feature maps + if (output_type == 'pandas'): - + header, string_data, numeric_data = featurize_fname_lists_imp (id(self), intensity_files, mask_files, single_roi, output_type, "") df = pd.concat( @@ -548,9 +598,9 @@ def featurize_files ( df.ROI_label = df.ROI_label.astype(np.uint32) return df - + else: - + featurize_fname_lists_imp (id(self), intensity_files, mask_files, single_roi, output_type, output_path) return get_arrow_file_imp (id(self)) @@ -775,19 +825,46 @@ def set_params(self, **params): for key, value in params.items(): if key.startswith("gabor_"): gabor_params[key[len("gabor_"):]] = value - + elif (key == "ibsi"): set_if_ibsi_imp (id(self), value) - + + elif key in ("fmaps", "fmaps_radius"): + pass # handled after loop to avoid double-processing + + elif key == "binning_origin": + if value not in ('zero', 'min'): + raise ValueError("binning_origin must be 'zero' or 'min'.") + # Get current depth magnitude + cur = get_params_imp(id(self), ["coarse_gray_depth"]) + depth = params.get("coarse_gray_depth", cur["coarse_gray_depth"]) + depth = abs(int(depth)) + environment_params["coarse_gray_depth"] = -depth if value == 'min' else depth + else: if (key not in available_environment_params): raise ValueError ("Invalid parameter: ", key) else: environment_params[key] = value - + + # If coarse_gray_depth was set without binning_origin, apply current binning_origin sign + if "coarse_gray_depth" in params and "binning_origin" not in params and "coarse_gray_depth" not in environment_params: + cur = get_params_imp(id(self), ["binning_origin"]) + depth = abs(int(params["coarse_gray_depth"])) + environment_params["coarse_gray_depth"] = -depth if cur.get("binning_origin") == 'min' else depth + + # Update fmaps settings if either param was provided. + # -1 sentinel means "leave unchanged" for that param. + if "fmaps" in params or "fmaps_radius" in params: + mode = int(bool(params["fmaps"])) if "fmaps" in params else -1 + radius = params.get("fmaps_radius", -1) + set_fmaps_imp(id(self), mode, radius) + if "fmaps" in params: + self._fmaps = bool(params["fmaps"]) + if (len(gabor_params) > 0): self.set_gabor_feature_params(**gabor_params) - + if (len(environment_params) > 0): self.set_environment_params(**environment_params) @@ -944,6 +1021,21 @@ class Nyxus3D: Y-dimension scale factor anisotropy_z: float (optional, default 1.0) Z-dimension scale factor + fmaps: bool (optional, default False) + Enable feature maps mode. Instead of computing a single feature vector per ROI, + a sliding kernel is moved across each ROI and features are computed at every + voxel position, producing spatial feature maps as numpy arrays. + When enabled, featurize methods return a list of dicts (one per parent ROI) + instead of a DataFrame. Each dict contains numpy arrays shaped (map_d, map_h, map_w). + fmaps_radius: int (optional, default 2) + Radius of the sliding kernel used in feature maps mode. + The kernel size is (2 * fmaps_radius + 1). For example, fmaps_radius=2 + produces a 5x5x5 kernel. Must be >= 1. + binning_origin: str (optional, default "zero") + Origin of the intensity binning range for texture features. + "zero" - bins span [0, max] (default Nyxus/MATLAB behavior). + "min" - bins span [min, max], adapting to the actual data range + (PyRadiomics-compatible behavior). """ def __init__( @@ -953,20 +1045,22 @@ def __init__( ): valid_keys = { 'neighbor_distance', 'pixels_per_micron', 'coarse_gray_depth', - 'n_feature_calc_threads', 'use_gpu_device', 'ibsi', 'channel_signature', - 'parent_channel', 'child_channel', 'aggregate', + 'n_feature_calc_threads', 'use_gpu_device', 'ibsi', 'channel_signature', + 'parent_channel', 'child_channel', 'aggregate', 'dynamic_range', 'min_intensity', 'max_intensity', 'ram_limit', 'verbose', - 'anisotropy_x', + 'anisotropy_x', 'anisotropy_y', - 'anisotropy_z' + 'anisotropy_z', + 'fmaps', 'fmaps_radius', + 'binning_origin' } # Check for unexpected keyword arguments invalid_keys = set(kwargs.keys()) - valid_keys if invalid_keys: print(f"Warning: unexpected keyword argument(s): {', '.join(invalid_keys)}") - + # parse kwargs features = features neighbor_distance = kwargs.get('neighbor_distance', 5) @@ -982,19 +1076,29 @@ def __init__( aniso_x = kwargs.get('anisotropy_x', 1.0) aniso_y = kwargs.get('anisotropy_y', 1.0) aniso_z = kwargs.get('anisotropy_z', 1.0) - + fmaps = kwargs.get('fmaps', False) + fmaps_radius = kwargs.get('fmaps_radius', 2) + binning_origin = kwargs.get('binning_origin', 'zero') + if neighbor_distance <= 0: raise ValueError("Neighbor distance must be greater than zero.") if pixels_per_micron <= 0: raise ValueError("Pixels per micron must be greater than zero.") - if coarse_gray_depth <= 0: - raise ValueError("Custom number of grayscale levels (parameter coarse_gray_depth, default=64) must be non-negative.") + if coarse_gray_depth < 1: + raise ValueError("coarse_gray_depth must be >= 1.") + + if binning_origin not in ('zero', 'min'): + raise ValueError("binning_origin must be 'zero' or 'min'.") + + # Negative coarse_gray_depth signals data-adaptive (min-based) binning to the C++ backend + if binning_origin == 'min': + coarse_gray_depth = -coarse_gray_depth if n_feature_calc_threads < 1: raise ValueError("There must be at least one feature calculation thread.") - + if(use_gpu_device > -1 and n_feature_calc_threads != 1): print("Gpu features only support a single thread. Defaulting to one thread.") n_feature_calc_threads = 1 @@ -1017,13 +1121,18 @@ def __init__( if aniso_z <= 0: raise ValueError ("anisotropy_z must be positive") + if fmaps and fmaps_radius < 1: + raise ValueError("fmaps_radius must be an integer >= 1") + + self._fmaps = fmaps + initialize_environment( id(self), 3, # 3D features, neighbor_distance, pixels_per_micron, - coarse_gray_depth, + coarse_gray_depth, n_feature_calc_threads, use_gpu_device, ibsi, @@ -1035,8 +1144,10 @@ def __init__( verb_lvl, aniso_x, aniso_y, - aniso_z) - + aniso_z, + fmaps, + fmaps_radius) + # list of valid outputs that are used throughout featurize functions self._valid_output_types = ['pandas', 'arrowipc', 'parquet'] @@ -1135,7 +1246,11 @@ def featurize_directory( if (output_type not in self._valid_output_types): raise ValueError(f'Invalid output type {output_type}. Valid output types are {self._valid_output_types}.') - + + if self._fmaps and output_type == 'pandas': + result = featurize_directory_3D_imp(id(self), intensity_dir, label_dir, file_pattern, output_type, "") + return result[0] # list of dicts with numpy array feature maps + if (output_type == 'pandas'): header, string_data, numeric_data = featurize_directory_3D_imp (id(self), intensity_dir, label_dir, file_pattern, output_type, "") df = pd.concat( @@ -1200,8 +1315,12 @@ def featurize_files ( if (output_type not in self._valid_output_types): raise ValueError(f'Invalid output type {output_type}. Valid output types are {self._valid_output_types}') + if self._fmaps and output_type == 'pandas': + result = featurize_fname_lists_3D_imp(id(self), intensity_files, mask_files, single_roi, output_type, "") + return result[0] # list of dicts with numpy array feature maps + if (output_type == 'pandas'): - + header, string_data, numeric_data = featurize_fname_lists_3D_imp (id(self), intensity_files, mask_files, single_roi, output_type, "") df = pd.concat( @@ -1217,9 +1336,9 @@ def featurize_files ( df.ROI_label = df.ROI_label.astype(np.uint32) return df - + else: - + featurize_fname_lists_3D_imp (id(self), intensity_files, mask_files, single_roi, output_type, output_path) return get_arrow_file_imp (id(self)) @@ -1326,17 +1445,30 @@ def set_params(self, **params): for key, value in params.items(): - + if (key == "ibsi"): set_if_ibsi_imp (id(self), value) - + + elif key == "binning_origin": + if value not in ('zero', 'min'): + raise ValueError("binning_origin must be 'zero' or 'min'.") + cur = get_params_imp(id(self), ["coarse_gray_depth"]) + depth = params.get("coarse_gray_depth", cur["coarse_gray_depth"]) + depth = abs(int(depth)) + environment_params["coarse_gray_depth"] = -depth if value == 'min' else depth + else: if (key not in available_environment_params): raise ValueError(f"Invalid parameter {key}.") else: environment_params[key] = value - - + + # If coarse_gray_depth was set without binning_origin, apply current binning_origin sign + if "coarse_gray_depth" in params and "binning_origin" not in params and "coarse_gray_depth" not in environment_params: + cur = get_params_imp(id(self), ["binning_origin"]) + depth = abs(int(params["coarse_gray_depth"])) + environment_params["coarse_gray_depth"] = -depth if cur.get("binning_origin") == 'min' else depth + if (len(environment_params) > 0): self.set_environment_params(**environment_params) @@ -1562,8 +1694,10 @@ def __init__( verb_lvl, aniso_x, aniso_y, - aniso_z) - + aniso_z, + False, + 2) + # list of valid outputs that are used throughout featurize functions self._valid_output_types = ['pandas', 'arrowipc', 'parquet'] diff --git a/src/nyx/results_cache.h b/src/nyx/results_cache.h index 75973105..f3a7c0bf 100644 --- a/src/nyx/results_cache.h +++ b/src/nyx/results_cache.h @@ -1,7 +1,26 @@ #pragma once +#include +#include #include #include +/// @brief Holds a spatial feature map for one parent ROI: one 2D/3D array per feature. +struct FmapArrayResult +{ + int parent_label; + std::string intens_name; + std::string seg_name; + int map_w, map_h; // dimensions of the feature map arrays + int map_d = 1; // depth dimension (1 for 2D feature maps) + int origin_x, origin_y; // global image coords of the map's top-left pixel + int origin_z = 0; // global z-coord of the map's first slice (unused in 2D) + std::vector feature_names; + // Flat storage: n_features contiguous blocks of (map_d * map_h * map_w) doubles each. + // 2D access: feature_data[f * map_h * map_w + row * map_w + col] + // 3D access: feature_data[f * map_d * map_h * map_w + z * map_h * map_w + row * map_w + col] + std::vector feature_data; +}; + class ResultsCache { public: @@ -14,6 +33,7 @@ class ResultsCache stringColBuf_.clear(); calcResultBuf_.clear(); totalNumLabels_ = 0; + fmapArrayResults_.clear(); } std::vector& get_headerBuf() { return headerBuf_; } @@ -25,7 +45,7 @@ class ResultsCache for (auto c : cols) add_to_header(c); } - void add_to_header(std::string& col) + void add_to_header(const std::string& col) { headerBuf_.push_back(col); } @@ -35,9 +55,12 @@ class ResultsCache void inc_num_rows() { totalNumLabels_++; } size_t get_num_rows() { return totalNumLabels_; } + std::vector& get_fmapArrayResults() { return fmapArrayResults_; } + private: std::vector calcResultBuf_; size_t totalNumLabels_ = 0; std::vector stringColBuf_, headerBuf_; + std::vector fmapArrayResults_; }; diff --git a/src/nyx/workflow_2d_fmaps.cpp b/src/nyx/workflow_2d_fmaps.cpp new file mode 100644 index 00000000..5e1095fb --- /dev/null +++ b/src/nyx/workflow_2d_fmaps.cpp @@ -0,0 +1,397 @@ +#include "helpers/fsystem.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef WITH_PYTHON_H + #include + #include + #include + namespace py = pybind11; +#endif + +#include "constants.h" +#include "dirs_and_files.h" +#include "environment.h" +#include "globals.h" +#include "helpers/helpers.h" +#include "helpers/system_resource.h" +#include "helpers/timing.h" +#include "raw_image_loader.h" + +namespace Nyxus +{ + + /// @brief Generates kernel-sized child ROIs by sliding a kernel across a parent ROI's image matrix. + /// Only creates child ROIs where the kernel center pixel is within the mask (non-zero intensity in the image matrix). + /// @param parent The parent ROI (must have aux_image_matrix populated) + /// @param kernel_size The kernel size (odd integer >= 3) + /// @param childLabels [out] Set of child ROI labels + /// @param childRoiData [out] Map of child label -> child LR + /// @param childToParentMap [out] Map of child label -> (parent_label, center_x, center_y) + /// @param startLabel Starting label for child ROIs to avoid collisions across parents + /// @return Number of child ROIs generated + // + // Disable optimizations for this function: at -O3 the loop vectorizer + // miscompiles inlined code (LR ctor, push_back, init_aabb, etc.) within + // the sliding-kernel loops, causing heap corruption. The function is + // dominated by memory allocation (LR construction, hash-map insertion) + // so vectorization provides no benefit here. +#pragma clang optimize off + int generateChildRois ( + const LR & parent, + int kernel_size, + // out: + std::unordered_set & childLabels, + std::unordered_map & childRoiData, + std::unordered_map & childToParentMap, + int64_t startLabel) + { + int half = kernel_size / 2; + int childCount = 0; + + int parentW = parent.aabb.get_width(); + int parentH = parent.aabb.get_height(); + int parentXmin = parent.aabb.get_xmin(); + int parentYmin = parent.aabb.get_ymin(); + + const auto& pixelMatrix = parent.aux_image_matrix.ReadablePixels(); + const auto* pixelData = pixelMatrix.data(); + int matW = parentW; // image matrix stride (width) + + // Slide kernel across the parent ROI's bounding box + for (int cy = half; cy < parentH - half; cy++) + { + for (int cx = half; cx < parentW - half; cx++) + { + // Check if the kernel center pixel is within the mask + auto centerVal = pixelData[matW * cy + cx]; + if (centerVal == 0) + continue; + + int64_t childLabel64 = startLabel + childCount; + if (childLabel64 > std::numeric_limits::max()) + throw std::overflow_error("Child ROI label overflow — too many child ROIs generated across all parents"); + int childLabel = (int)childLabel64; + + // Create child LR + LR child(childLabel); + + // Set up child's AABB in image coordinates + int childXmin = parentXmin + cx - half; + int childYmin = parentYmin + cy - half; + int childXmax = parentXmin + cx + half; + int childYmax = parentYmin + cy + half; + + child.init_aabb(childXmin, childYmin); + child.update_aabb(childXmax, childYmax); + child.make_nonanisotropic_aabb(); + + // Extract pixels from the parent's image matrix + double cmin = std::numeric_limits::max(); + double cmax = std::numeric_limits::lowest(); + for (int ky = 0; ky < kernel_size; ky++) + { + for (int kx = 0; kx < kernel_size; kx++) + { + int localY = cy - half + ky; + int localX = cx - half + kx; + auto val = pixelData[matW * localY + localX]; + if (val != 0) + { + int imgX = parentXmin + localX; + int imgY = parentYmin + localY; + child.raw_pixels.push_back(Pixel2(imgX, imgY, val)); + child.aux_area++; + if (val < cmin) cmin = val; + if (val > cmax) cmax = val; + } + } + } + if (child.aux_area > 0) + { + // Use parent ROI's min/max so all kernel patches share + // the same binning range (consistent with pyradiomics + // fixed-width global binning). Local cmin/cmax would + // cause each tiny patch to get its own adaptive range, + // collapsing most pixels to the same bin. + child.aux_min = parent.aux_min; + child.aux_max = parent.aux_max; + } + + // Allocate and populate image matrix + child.aux_image_matrix.allocate(kernel_size, kernel_size); + child.aux_image_matrix.calculate_from_pixelcloud(child.raw_pixels, child.aabb); + + // Initialize feature value buffers + child.initialize_fvals(); + + // Store + childLabels.insert(childLabel); + childRoiData[childLabel] = std::move(child); + childToParentMap[childLabel] = { parent.label, parentXmin + cx, parentYmin + cy }; + + childCount++; + } + } + + return childCount; + } +#pragma clang optimize on + + /// @brief Processes a single intensity-segmentation image pair in feature maps mode. + /// Loads parent ROIs, generates child ROIs per kernel, computes features, outputs results. + /// @param globalChildLabel [in/out] Running child label counter, persists across file pairs to avoid label collisions. + static bool processIntSegImagePair_fmaps ( + Environment & env, + const std::string & intens_fpath, + const std::string & label_fpath, + int filepair_index, + int tot_num_filepairs, + int64_t & globalChildLabel) + { + // Display progress + if (tot_num_filepairs > 0) + { + int digits = (tot_num_filepairs >= 100) ? (int)std::log10(tot_num_filepairs/100.) + 1 : 1, + k = (int)std::pow(10.f, digits); + float perCent = float(filepair_index + 1) * 100.f / float(tot_num_filepairs); + perCent = std::round(perCent * k) / k; + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "[ " << filepair_index+1 << " = " << std::setw(digits + 2) << perCent << "% ]\t" << intens_fpath << "\n") + } + + // Phase 1: gather ROI metrics (unchanged) + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "Gathering ROI metrics\n"); + bool okGather = gatherRoisMetrics (filepair_index, intens_fpath, label_fpath, env, env.theImLoader); + if (!okGather) + { + std::string msg = "Error gathering ROI metrics from " + intens_fpath + " / " + label_fpath + "\n"; + std::cerr << msg; + throw (std::runtime_error(msg)); + } + + if (env.uniqueLabels.size() == 0) + return true; + + // Collect parent ROI labels (all trivial for fmaps — kernel-sized children are always small) + std::vector parentLabels; + for (auto lab : env.uniqueLabels) + { + LR& r = env.roiData[lab]; + if (env.roi_is_blacklisted("", lab)) + { + r.blacklisted = true; + continue; + } + parentLabels.push_back(lab); + } + + // Phase 2 (fmaps): Process each parent ROI + for (auto parentLab : parentLabels) + { + LR& parentROI = env.roiData[parentLab]; + + // Check that the parent ROI is large enough for the kernel + int parentW = parentROI.aabb.get_width(); + int parentH = parentROI.aabb.get_height(); + if (parentW < env.fmaps_kernel_size() || parentH < env.fmaps_kernel_size()) + { + VERBOSLVL2 (env.get_verbosity_level(), + std::cout << "Skipping ROI " << parentLab + << " (too small for kernel: " << parentW << "x" << parentH + << " < " << env.fmaps_kernel_size() << ")\n"); + continue; + } + + // Step a: Load parent ROI pixels + std::vector singleParent = { parentLab }; + + if (env.anisoOptions.customized() == false) + scanTrivialRois (singleParent, intens_fpath, label_fpath, env, env.theImLoader); + else + { + double ax = env.anisoOptions.get_aniso_x(), + ay = env.anisoOptions.get_aniso_y(); + scanTrivialRois_anisotropic (singleParent, intens_fpath, label_fpath, env, env.theImLoader, ax, ay); + } + + // Step b: Allocate image matrix for the parent + allocateTrivialRoisBuffers (singleParent, env.roiData, env.hostCache); + + // Step c: Generate child ROIs + std::unordered_set childLabels; + std::unordered_map childRoiData; + std::unordered_map childToParentMap; + + int nChildren = generateChildRois ( + parentROI, + env.fmaps_kernel_size(), + childLabels, + childRoiData, + childToParentMap, + globalChildLabel); + + VERBOSLVL2 (env.get_verbosity_level(), + std::cout << "ROI " << parentLab << ": generated " << nChildren << " child ROIs\n"); + + globalChildLabel += nChildren; + + if (nChildren > 0) + { + // Capture parent geometry before the swap invalidates the reference + int parentXmin = parentROI.aabb.get_xmin(); + int parentYmin = parentROI.aabb.get_ymin(); + + // Step d: Compute features on child ROIs + // RAII guard swaps env's ROI data with child data and restores on scope exit (even on exception) + EnvRoiSwapGuard guard (env, std::move(childLabels), std::move(childRoiData)); + + std::vector childLabelVec(env.uniqueLabels.begin(), env.uniqueLabels.end()); + std::sort(childLabelVec.begin(), childLabelVec.end()); + + reduce_trivial_rois_manual (childLabelVec, env); + + // Step e: Output child ROI features + if (env.saveOption == SaveOption::saveBuffer) + { + save_features_2_fmap_arrays ( + env.theResultsCache, + env, + intens_fpath, + label_fpath, + parentLab, + parentXmin, parentYmin, + parentW, parentH, + env.fmaps_kernel_size(), + env.uniqueLabels, + env.roiData, + childToParentMap); + } + else + { + save_features_2_csv_fmaps ( + env, + intens_fpath, + label_fpath, + env.output_dir, + filepair_index, + env.uniqueLabels, + env.roiData, + childToParentMap); + } + } + + // Free parent ROI buffers + freeTrivialRoisBuffers (singleParent, env.roiData); + +#ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + { + sureprint("\nAborting per user input\n"); + throw pybind11::error_already_set(); + } +#endif + } + + return true; + } + + int processDataset_2D_fmaps ( + Environment & env, + const std::vector & intensFiles, + const std::vector & labelFiles, + int numReduceThreads, + const SaveOption saveOption, + const std::string & outputPath) + { + reset_csv_header_state(); + +#ifdef CHECKTIMING + if (Stopwatch::inclusive()) + Stopwatch::reset(); +#endif + + //********************** prescan (unchanged from segmented workflow) *********************** + + size_t nf = intensFiles.size(); + + { STOPWATCH("prescan/p0/P/#ccbbaa", "\t="); + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "phase 0 (prescanning)\n"); + + env.dataset.reset_dataset_props(); + + for (size_t i = 0; i < nf; i++) + { + SlideProps& p = env.dataset.dataset_props.emplace_back (intensFiles[i], labelFiles[i]); + + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "prescanning " << p.fname_int); + if (! scan_slide_props(p, 2, env.anisoOptions, env.resultOptions.need_annotation())) + { + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "error prescanning pair " << p.fname_int << " and " << p.fname_seg << std::endl); + return 1; + } + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "\t " << p.slide_w << " W x " << p.slide_h << " H\tmax ROI " << p.max_roi_w << " x " << p.max_roi_h << "\tmin-max I " << Nyxus::virguler_real(p.min_preroi_inten) << "-" << Nyxus::virguler_real(p.max_preroi_inten) << "\t" << p.lolvl_slide_descr << "\n"); + } + + env.dataset.update_dataset_props_extrema(); + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "\t finished prescanning \n"); + } // prescan timing + + // One-time initialization + init_slide_rois (env.uniqueLabels, env.roiData); + + bool ok = true; + int64_t globalChildLabel = 1; // Persists across file pairs to avoid label collisions + + // Iterate intensity-segmentation pairs + for (int i = 0; i < nf; i++) + { +#ifdef CHECKTIMING + if (Stopwatch::exclusive()) + Stopwatch::reset(); +#endif + + clear_slide_rois (env.uniqueLabels, env.roiData); + + auto& ifp = intensFiles[i], + & lfp = labelFiles[i]; + + SlideProps& p = env.dataset.dataset_props[i]; + ok = env.theImLoader.open (p, env.fpimageOptions); + if (ok == false) + { + std::cerr << "Terminating\n"; + return 1; + } + + // Feature maps processing + if (! processIntSegImagePair_fmaps(env, ifp, lfp, i, nf, globalChildLabel)) + { + std::cerr << "Error in feature maps processing for " << ifp << " @ " << __FILE__ << ":" << __LINE__ << "\n"; + return 1; + } + + env.theImLoader.close(); + +#ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + { + sureprint("\nAborting per user input\n"); + throw pybind11::error_already_set(); + } +#endif + } + + return 0; // success + } + +} // namespace Nyxus diff --git a/src/nyx/workflow_2d_segmented.cpp b/src/nyx/workflow_2d_segmented.cpp index ad747ec0..45b27aa8 100644 --- a/src/nyx/workflow_2d_segmented.cpp +++ b/src/nyx/workflow_2d_segmented.cpp @@ -166,11 +166,12 @@ namespace Nyxus const SaveOption saveOption, const std::string& outputPath) { + reset_csv_header_state(); #ifdef CHECKTIMING if (Stopwatch::inclusive()) Stopwatch::reset(); -#endif +#endif //********************** prescan *********************** diff --git a/src/nyx/workflow_2d_whole.cpp b/src/nyx/workflow_2d_whole.cpp index 7db03583..3542d4b1 100644 --- a/src/nyx/workflow_2d_whole.cpp +++ b/src/nyx/workflow_2d_whole.cpp @@ -207,6 +207,8 @@ namespace Nyxus const SaveOption saveOption, const std::string& outputPath) { + reset_csv_header_state(); + //**** prescan all slides size_t nf = intensFiles.size(); diff --git a/src/nyx/workflow_3d_fmaps.cpp b/src/nyx/workflow_3d_fmaps.cpp new file mode 100644 index 00000000..b36ba03b --- /dev/null +++ b/src/nyx/workflow_3d_fmaps.cpp @@ -0,0 +1,378 @@ +#include "helpers/fsystem.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef WITH_PYTHON_H + #include + #include + #include + namespace py = pybind11; +#endif + +#include "constants.h" +#include "dirs_and_files.h" +#include "environment.h" +#include "globals.h" +#include "helpers/helpers.h" +#include "helpers/system_resource.h" +#include "helpers/timing.h" +#include "raw_image_loader.h" + +namespace Nyxus +{ + + /// @brief Generates kernel-sized child ROIs by sliding a 3D kernel across a parent ROI's image cube. + /// Only creates child ROIs where the kernel center voxel is non-zero. + /// @param parent The parent ROI (must have aux_image_cube populated) + /// @param kernel_size The kernel size (odd integer >= 3) + /// @param childLabels [out] Set of child ROI labels + /// @param childRoiData [out] Map of child label -> child LR + /// @param childToParentMap [out] Map of child label -> (parent_label, center_x, center_y, center_z) + /// @param startLabel Starting label for child ROIs to avoid collisions across parents + /// @return Number of child ROIs generated +#pragma clang optimize off + int generateChildRois_3D ( + const LR & parent, + int kernel_size, + // out: + std::unordered_set & childLabels, + std::unordered_map & childRoiData, + std::unordered_map & childToParentMap, + int64_t startLabel) + { + int half = kernel_size / 2; + int childCount = 0; + + int parentW = parent.aabb.get_width(); + int parentH = parent.aabb.get_height(); + int parentD = parent.aabb.get_z_depth(); + int parentXmin = parent.aabb.get_xmin(); + int parentYmin = parent.aabb.get_ymin(); + int parentZmin = parent.aabb.get_zmin(); + + const auto& cube = parent.aux_image_cube; + + // Slide kernel across the parent ROI's bounding box in 3D + for (int cz = half; cz < parentD - half; cz++) + { + for (int cy = half; cy < parentH - half; cy++) + { + for (int cx = half; cx < parentW - half; cx++) + { + // Check if the kernel center voxel is within the mask + auto centerVal = cube.xyz(cx, cy, cz); + if (centerVal == 0) + continue; + + int64_t childLabel64 = startLabel + childCount; + if (childLabel64 > std::numeric_limits::max()) + throw std::overflow_error("Child ROI label overflow — too many child ROIs generated across all parents"); + int childLabel = (int)childLabel64; + + // Create child LR + LR child(childLabel); + + // Set up child's AABB in image coordinates + int childXmin = parentXmin + cx - half; + int childYmin = parentYmin + cy - half; + int childZmin = parentZmin + cz - half; + int childXmax = parentXmin + cx + half; + int childYmax = parentYmin + cy + half; + int childZmax = parentZmin + cz + half; + + child.init_aabb_3D(childXmin, childYmin, childZmin); + child.update_aabb_3D(childXmax, childYmax, childZmax); + child.make_nonanisotropic_aabb(); + + // Extract voxels from the parent's image cube + double cmin = std::numeric_limits::max(); + double cmax = std::numeric_limits::lowest(); + for (int kz = 0; kz < kernel_size; kz++) + { + for (int ky = 0; ky < kernel_size; ky++) + { + for (int kx = 0; kx < kernel_size; kx++) + { + int localX = cx - half + kx; + int localY = cy - half + ky; + int localZ = cz - half + kz; + auto val = cube.xyz(localX, localY, localZ); + if (val != 0) + { + int imgX = parentXmin + localX; + int imgY = parentYmin + localY; + int imgZ = parentZmin + localZ; + child.raw_pixels_3D.push_back(Pixel3(imgX, imgY, imgZ, val)); + child.aux_area++; + if (val < cmin) cmin = val; + if (val > cmax) cmax = val; + } + } + } + } + if (child.aux_area > 0) + { + // Use parent ROI's min/max so all kernel patches share + // the same binning range (consistent with pyradiomics + // fixed-width global binning). + child.aux_min = parent.aux_min; + child.aux_max = parent.aux_max; + } + + // Allocate and populate image cube + child.aux_image_cube.allocate(kernel_size, kernel_size, kernel_size); + child.aux_image_cube.calculate_from_pixelcloud(child.raw_pixels_3D, child.aabb); + + // Initialize feature value buffers + child.initialize_fvals(); + + // Store + childLabels.insert(childLabel); + childRoiData[childLabel] = std::move(child); + childToParentMap[childLabel] = { parent.label, parentXmin + cx, parentYmin + cy, parentZmin + cz }; + + childCount++; + } + } + } + + return childCount; + } +#pragma clang optimize on + + /// @brief Processes a single intensity-segmentation image pair in 3D feature maps mode. + static bool processIntSegImagePair_3D_fmaps ( + Environment & env, + const std::string & intens_fpath, + const std::string & label_fpath, + size_t filepair_index, + size_t t_index, + const std::vector& z_indices, + int64_t & globalChildLabel) + { + // Phase 1: gather ROI metrics + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "Gathering ROI metrics (3D fmaps)\n"); + bool okGather = false; + if (z_indices.size()) + okGather = gatherRoisMetrics_25D (env, filepair_index, intens_fpath, label_fpath, z_indices); + else + okGather = gatherRoisMetrics_3D (env, filepair_index, intens_fpath, label_fpath, t_index); + if (!okGather) + { + std::string msg = "Error gathering ROI metrics from " + intens_fpath + " / " + label_fpath + "\n"; + std::cerr << msg; + throw (std::runtime_error(msg)); + } + + if (env.uniqueLabels.size() == 0) + return true; + + // Collect parent ROI labels + std::vector parentLabels; + for (auto lab : env.uniqueLabels) + { + LR& r = env.roiData[lab]; + r.initialize_fvals(); + if (env.roi_is_blacklisted("", lab)) + { + r.blacklisted = true; + continue; + } + parentLabels.push_back(lab); + } + + // Phase 2 (fmaps): Process each parent ROI + for (auto parentLab : parentLabels) + { + LR& parentROI = env.roiData[parentLab]; + + // Check that the parent ROI is large enough for the kernel in all dimensions + int parentW = parentROI.aabb.get_width(); + int parentH = parentROI.aabb.get_height(); + int parentD = parentROI.aabb.get_z_depth(); + if (parentW < env.fmaps_kernel_size() || parentH < env.fmaps_kernel_size() || parentD < env.fmaps_kernel_size()) + { + VERBOSLVL2 (env.get_verbosity_level(), + std::cout << "Skipping ROI " << parentLab + << " (too small for kernel: " << parentW << "x" << parentH << "x" << parentD + << " < " << env.fmaps_kernel_size() << ")\n"); + continue; + } + + // Step a: Load parent ROI voxels and allocate image cube + std::vector singleParent = { parentLab }; + + if (env.anisoOptions.customized() == false) + scanTrivialRois_3D (env, singleParent, intens_fpath, label_fpath, t_index); + else + { + double ax = env.anisoOptions.get_aniso_x(), + ay = env.anisoOptions.get_aniso_y(), + az = env.anisoOptions.get_aniso_z(); + scanTrivialRois_3D_anisotropic (env, singleParent, intens_fpath, label_fpath, t_index, ax, ay, az); + } + + allocateTrivialRoisBuffers_3D (singleParent, env.roiData, env.hostCache); + + // Step b: Generate child ROIs + std::unordered_set childLabels; + std::unordered_map childRoiData; + std::unordered_map childToParentMap; + + int nChildren = generateChildRois_3D ( + parentROI, + env.fmaps_kernel_size(), + childLabels, + childRoiData, + childToParentMap, + globalChildLabel); + + VERBOSLVL2 (env.get_verbosity_level(), + std::cout << "ROI " << parentLab << ": generated " << nChildren << " child ROIs (3D)\n"); + + globalChildLabel += nChildren; + + if (nChildren > 0) + { + // Capture parent geometry before the swap invalidates the reference + int parentXmin = parentROI.aabb.get_xmin(); + int parentYmin = parentROI.aabb.get_ymin(); + int parentZmin = parentROI.aabb.get_zmin(); + + // Step c: Compute features on child ROIs + EnvRoiSwapGuard guard (env, std::move(childLabels), std::move(childRoiData)); + + std::vector childLabelVec(env.uniqueLabels.begin(), env.uniqueLabels.end()); + std::sort(childLabelVec.begin(), childLabelVec.end()); + + reduce_trivial_rois_manual (childLabelVec, env); + + // Step d: Output child ROI features + if (env.saveOption == SaveOption::saveBuffer) + { + save_features_2_fmap_arrays_3D ( + env.theResultsCache, + env, + intens_fpath, + label_fpath, + parentLab, + parentXmin, parentYmin, parentZmin, + parentW, parentH, parentD, + env.fmaps_kernel_size(), + env.uniqueLabels, + env.roiData, + childToParentMap); + } + else + { + save_features_2_csv_fmaps ( + env, + intens_fpath, + label_fpath, + env.output_dir, + filepair_index, + env.uniqueLabels, + env.roiData, + childToParentMap); + } + } + + // Free parent ROI buffers + freeTrivialRoisBuffers_3D (singleParent, env.roiData); + +#ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + { + sureprint("\nAborting per user input\n"); + throw pybind11::error_already_set(); + } +#endif + } + + return true; + } + + int processDataset_3D_fmaps ( + Environment & env, + const std::vector & intensFiles, + const std::vector & labelFiles, + int numReduceThreads, + const SaveOption saveOption, + const std::string & outputPath) + { + reset_csv_header_state(); + + //********************** prescan *********************** + + size_t nf = intensFiles.size(); + + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "phase 0 (3D fmaps prescanning)\n"); + env.dataset.reset_dataset_props(); + + for (size_t i = 0; i < nf; i++) + { + SlideProps& p = env.dataset.dataset_props.emplace_back (intensFiles[i].fdir + intensFiles[i].fname, labelFiles[i].fdir + labelFiles[i].fname); + + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "prescanning " << p.fname_int); + if (! scan_slide_props(p, 3, env.anisoOptions, env.resultOptions.need_annotation())) + { + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "error prescanning pair " << p.fname_int << " and " << p.fname_seg << std::endl); + return 1; + } + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "\t " << p.slide_w << " W x " << p.slide_h << " H x " << p.volume_d << " D\n"); + } + + env.dataset.update_dataset_props_extrema(); + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "\t finished prescanning \n"); + + // One-time initialization + init_slide_rois (env.uniqueLabels, env.roiData); + + bool ok = true; + int64_t globalChildLabel = 1; + + // Iterate intensity-mask pairs + for (size_t i = 0; i < nf; i++) + { + // Iterate time frames + for (size_t t = 0; t < env.dataset.dataset_props[i].inten_time; t++) + { + clear_slide_rois (env.uniqueLabels, env.roiData); + + auto& ifile = intensFiles[i], + & mfile = labelFiles[i]; + + int digits = 2, k = (int)std::pow(10.f, digits); + float perCent = float(i * 100 * k / nf) / float(k); + VERBOSLVL1 (env.get_verbosity_level(), std::cout << "[ " << std::setw(digits + 2) << perCent << "% ]\t" << " INT: " << ifile.fname << " SEG: " << mfile.fname << " T:" << t << "\n") + + if (! processIntSegImagePair_3D_fmaps(env, ifile.fdir+ifile.fname, mfile.fdir+mfile.fname, i, t, intensFiles[i].z_indices, globalChildLabel)) + { + std::cerr << "Error in 3D feature maps processing for " << ifile.fname << " @ " << __FILE__ << ":" << __LINE__ << "\n"; + return 1; + } + +#ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + { + sureprint("\nAborting per user input\n"); + throw pybind11::error_already_set(); + } +#endif + } //- time frames + } //- inten-mask pairs + + return 0; // success + } + +} // namespace Nyxus diff --git a/src/nyx/workflow_pythonapi.cpp b/src/nyx/workflow_pythonapi.cpp index 2f38a515..3499a7fa 100644 --- a/src/nyx/workflow_pythonapi.cpp +++ b/src/nyx/workflow_pythonapi.cpp @@ -82,6 +82,140 @@ namespace Nyxus return true; } + /// @brief In-memory feature maps processing for a single image pair. + /// Gathers parent ROIs, generates child ROIs per kernel, computes features, saves to buffer. + /// @param globalChildLabel [in/out] Running child label counter, persists across file pairs to avoid label collisions. + bool processIntSegImagePairInMemory_fmaps ( + Environment & env, + const py::array_t& intens, + const py::array_t& label, + int pair_index, + const std::string& intens_name, + const std::string& seg_name, + int64_t & globalChildLabel) + { + // Phase 1: gather parent ROI metrics + if (! gatherRoisMetricsInMemory (env, intens, label, pair_index)) + return false; + + if (env.uniqueLabels.size() == 0) + return true; + + // Set up parent ROIs + for (auto lab : env.uniqueLabels) + { + LR& r = env.roiData[lab]; + r.make_nonanisotropic_aabb(); + r.initialize_fvals(); + } + + // Collect parent labels, skipping blacklisted and oversized ROIs + std::vector parentLabels; + for (auto lab : env.uniqueLabels) + { + LR& r = env.roiData[lab]; + if (env.roi_is_blacklisted("", lab)) + { + r.blacklisted = true; + continue; + } + + // Skip oversized ROIs that exceed RAM limit (matching non-fmaps behavior) + size_t roiFootprint = r.get_ram_footprint_estimate(env.uniqueLabels.size()); + size_t ramLim = env.get_ram_limit(); + if (roiFootprint >= ramLim) + { + VERBOSLVL1 (env.get_verbosity_level(), + std::cout << "Skipping oversized ROI " << lab + << " (estimated " << roiFootprint << " bytes >= RAM limit " << ramLim << " bytes)\n"); + continue; + } + + parentLabels.push_back(lab); + } + + // Phase 2: process each parent ROI + for (auto parentLab : parentLabels) + { + LR& parentROI = env.roiData[parentLab]; + + int parentW = parentROI.aabb.get_width(); + int parentH = parentROI.aabb.get_height(); + if (parentW < env.fmaps_kernel_size() || parentH < env.fmaps_kernel_size()) + { + VERBOSLVL2 (env.get_verbosity_level(), + std::cout << "Skipping ROI " << parentLab + << " (too small for kernel: " << parentW << "x" << parentH + << " < " << env.fmaps_kernel_size() << ")\n"); + continue; + } + + // Scan parent ROI pixels from in-memory arrays + std::vector singleParent = { parentLab }; + scanTrivialRoisInMemory (singleParent, intens, label, pair_index, env); + + // Allocate image matrix for parent + allocateTrivialRoisBuffers (singleParent, env.roiData, env.hostCache); + + // Generate child ROIs + std::unordered_set childLabels; + std::unordered_map childRoiData; + std::unordered_map childToParentMap; + + int nChildren = generateChildRois ( + parentROI, + env.fmaps_kernel_size(), + childLabels, + childRoiData, + childToParentMap, + globalChildLabel); + + VERBOSLVL2 (env.get_verbosity_level(), + std::cout << "ROI " << parentLab << ": generated " << nChildren << " child ROIs\n"); + + globalChildLabel += nChildren; + + if (nChildren > 0) + { + // Capture parent geometry before the swap invalidates the reference + int parentXmin = parentROI.aabb.get_xmin(); + int parentYmin = parentROI.aabb.get_ymin(); + + // RAII guard swaps env's ROI data with child data and restores on scope exit (even on exception) + EnvRoiSwapGuard guard (env, std::move(childLabels), std::move(childRoiData)); + + std::vector childLabelVec(env.uniqueLabels.begin(), env.uniqueLabels.end()); + std::sort(childLabelVec.begin(), childLabelVec.end()); + + // Compute features on child ROIs + reduce_trivial_rois_manual (childLabelVec, env); + + // Save as spatial feature map arrays + save_features_2_fmap_arrays ( + env.theResultsCache, + env, + intens_name, + seg_name, + parentLab, + parentXmin, parentYmin, + parentW, parentH, + env.fmaps_kernel_size(), + env.uniqueLabels, + env.roiData, + childToParentMap); + } + + // Free parent ROI buffers + freeTrivialRoisBuffers (singleParent, env.roiData); + + // Allow keyboard interrupt + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); + } + + return true; + } + std::optional processMontage( Environment & env, const py::array_t& intensity_images, @@ -91,16 +225,19 @@ namespace Nyxus const std::vector& seg_names, const SaveOption saveOption, const std::string& outputPath) - { + { // prepare the output bool write_apache = (saveOption == SaveOption::saveArrowIPC || saveOption == SaveOption::saveParquet); - if (write_apache) + if (env.fmaps_prevents_arrow()) + return { "Arrow/Parquet output is not supported in feature maps (fmaps) mode. Use CSV or buffer output instead." }; + + if (write_apache) { env.arrow_stream = ArrowOutputStream(); auto [status, msg] = env.arrow_stream.create_arrow_file (saveOption, get_arrow_filename(outputPath, env.nyxus_result_fname, saveOption), Nyxus::get_header(env)); - if (!status) + if (!status) return { "error creating Arrow file: " + msg.value() }; } @@ -129,44 +266,55 @@ namespace Nyxus VERBOSLVL1(env.get_verbosity_level(), std::cout << "\t finished prescanning \n"); //****** extract features - + + int64_t globalChildLabel = 1; // Persists across file pairs to avoid label collisions in fmaps mode + for (int i_pair = 0; i_pair < n_pairs; i_pair++) { VERBOSLVL4 (env.get_verbosity_level(), std::cout << "processMontage() pair " << i_pair << "/" << n_pairs << "\n"); clear_slide_rois (env.uniqueLabels, env.roiData); // Clear ROI label list, ROI data, etc. - std::vector unprocessed_rois; - - if (! processIntSegImagePairInMemory (env, intensity_images, label_images, i_pair, intensity_names[i_pair], seg_names[i_pair], unprocessed_rois)) - return { "error processing a slide pair" }; - - if (write_apache) + if (env.fmaps_mode) { - auto [status, msg] = env.arrow_stream.write_arrow_file (Nyxus::get_feature_values(env.theFeatureSet, env.uniqueLabels, env.roiData, env.dataset)); - if (!status) - return { "error writing Arrow file: " + msg.value() }; - } - else - { - if (! save_features_2_buffer(env.theResultsCache, env, DEFAULT_T_INDEX)) - return { "error saving results to a buffer" }; + // Feature maps mode: generate child ROIs and compute features + if (! processIntSegImagePairInMemory_fmaps (env, intensity_images, label_images, i_pair, intensity_names[i_pair], seg_names[i_pair], globalChildLabel)) + return { "error processing fmaps for a slide pair" }; } - - if (unprocessed_rois.size() > 0) + else { - std::string erm = "the following ROIS are oversized and cannot be processed: "; - for (const auto& roi: unprocessed_rois) + std::vector unprocessed_rois; + + if (! processIntSegImagePairInMemory (env, intensity_images, label_images, i_pair, intensity_names[i_pair], seg_names[i_pair], unprocessed_rois)) + return { "error processing a slide pair" }; + + if (write_apache) { - erm += std::to_string(roi); - erm += ", "; + auto [status, msg] = env.arrow_stream.write_arrow_file (Nyxus::get_feature_values(env.theFeatureSet, env.uniqueLabels, env.roiData, env.dataset)); + if (!status) + return { "error writing Arrow file: " + msg.value() }; + } + else + { + if (! save_features_2_buffer(env.theResultsCache, env, DEFAULT_T_INDEX)) + return { "error saving results to a buffer" }; } - - // remove trailing space and comma - erm.pop_back(); - erm.pop_back(); - return { erm }; + if (unprocessed_rois.size() > 0) + { + std::string erm = "the following ROIS are oversized and cannot be processed: "; + for (const auto& roi: unprocessed_rois) + { + erm += std::to_string(roi); + erm += ", "; + } + + // remove trailing space and comma + erm.pop_back(); + erm.pop_back(); + + return { erm }; + } } // allow keyboard interrupt @@ -174,7 +322,7 @@ namespace Nyxus throw pybind11::error_already_set(); } - if (write_apache) + if (write_apache) { // close arrow file after use auto [status, msg] = env.arrow_stream.close_arrow_file(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 32acafab..5b0d3555 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 3.20) -set (CMAKE_CXX_STANDARD 17) +set (CMAKE_CXX_STANDARD 20) message(STATUS "Build tests") #==== Compiler Options set(CMAKE_CXX_STANDARD 17) diff --git a/tests/test_all.cc b/tests/test_all.cc index f40d39f9..03407f4d 100644 --- a/tests/test_all.cc +++ b/tests/test_all.cc @@ -31,12 +31,102 @@ #include "test_compat_3d_ngtdm.h" #include "test_compat_3d_glrlm.h" #include "test_compat_3d_glszm.h" +#include "test_fmaps.h" +#include "test_fmaps_3d.h" +#include "test_binning_origin.h" #ifdef USE_ARROW #include "test_arrow.h" #include "test_arrow_file_name.h" #endif +//***** Feature maps ***** + +TEST(TEST_NYXUS, TEST_FMAPS_CHILD_ROI_COUNT) { + ASSERT_NO_THROW(Nyxus::test_fmaps_child_roi_count()); +} + +TEST(TEST_NYXUS, TEST_FMAPS_CHILD_ROI_DIMENSIONS) { + ASSERT_NO_THROW(Nyxus::test_fmaps_child_roi_dimensions()); +} + +TEST(TEST_NYXUS, TEST_FMAPS_PARENT_MAPPING) { + ASSERT_NO_THROW(Nyxus::test_fmaps_parent_mapping()); +} + +TEST(TEST_NYXUS, TEST_FMAPS_CHILD_PIXEL_VALUES) { + ASSERT_NO_THROW(Nyxus::test_fmaps_child_pixel_values()); +} + +TEST(TEST_NYXUS, TEST_FMAPS_NONORIGIN_PARENT) { + ASSERT_NO_THROW(Nyxus::test_fmaps_nonorigin_parent()); +} + +TEST(TEST_NYXUS, TEST_FMAPS_SPARSE_MASK) { + ASSERT_NO_THROW(Nyxus::test_fmaps_sparse_mask()); +} + +TEST(TEST_NYXUS, TEST_FMAPS_PARENT_TOO_SMALL) { + ASSERT_NO_THROW(Nyxus::test_fmaps_parent_too_small()); +} + +TEST(TEST_NYXUS, TEST_FMAPS_START_LABEL_OFFSET) { + ASSERT_NO_THROW(Nyxus::test_fmaps_start_label_offset()); +} + +//***** 3D Feature maps ***** + +TEST(TEST_NYXUS, TEST_FMAPS3D_CHILD_ROI_COUNT) { + ASSERT_NO_THROW(Nyxus::test_fmaps3d_child_roi_count()); +} + +TEST(TEST_NYXUS, TEST_FMAPS3D_CHILD_ROI_DIMENSIONS) { + ASSERT_NO_THROW(Nyxus::test_fmaps3d_child_roi_dimensions()); +} + +TEST(TEST_NYXUS, TEST_FMAPS3D_PARENT_MAPPING) { + ASSERT_NO_THROW(Nyxus::test_fmaps3d_parent_mapping()); +} + +TEST(TEST_NYXUS, TEST_FMAPS3D_CHILD_VOXEL_VALUES) { + ASSERT_NO_THROW(Nyxus::test_fmaps3d_child_voxel_values()); +} + +TEST(TEST_NYXUS, TEST_FMAPS3D_NONORIGIN_PARENT) { + ASSERT_NO_THROW(Nyxus::test_fmaps3d_nonorigin_parent()); +} + +TEST(TEST_NYXUS, TEST_FMAPS3D_SPARSE_MASK) { + ASSERT_NO_THROW(Nyxus::test_fmaps3d_sparse_mask()); +} + +TEST(TEST_NYXUS, TEST_FMAPS3D_PARENT_TOO_SMALL) { + ASSERT_NO_THROW(Nyxus::test_fmaps3d_parent_too_small()); +} + +TEST(TEST_NYXUS, TEST_FMAPS3D_START_LABEL_OFFSET) { + ASSERT_NO_THROW(Nyxus::test_fmaps3d_start_label_offset()); +} + + +//***** Binning origin ***** + +TEST(TEST_NYXUS, TEST_BIN_PIXEL_ZERO_ORIGIN) { + ASSERT_NO_THROW(Nyxus::test_bin_pixel_zero_origin()); +} + +TEST(TEST_NYXUS, TEST_BIN_PIXEL_MIN_ORIGIN) { + ASSERT_NO_THROW(Nyxus::test_bin_pixel_min_origin()); +} + +TEST(TEST_NYXUS, TEST_GLCM_BINNING_ORIGIN_DIVERGENCE) { + ASSERT_NO_THROW(Nyxus::test_glcm_binning_origin_divergence()); +} + +TEST(TEST_NYXUS, TEST_NGTDM_BINNING_ORIGIN_DIVERGENCE) { + ASSERT_NO_THROW(Nyxus::test_ngtdm_binning_origin_divergence()); +} + //***** 2D contour and multicontour ***** TEST(TEST_NYXUS, TEST_CONTOUR_MULTI_1) { diff --git a/tests/test_binning_origin.h b/tests/test_binning_origin.h new file mode 100644 index 00000000..bb57cc14 --- /dev/null +++ b/tests/test_binning_origin.h @@ -0,0 +1,145 @@ +#pragma once + +#include +#include "../src/nyx/features/texture_feature.h" +#include "../src/nyx/features/glcm.h" +#include "../src/nyx/features/ngtdm.h" +#include "../src/nyx/roi_cache.h" +#include "test_data.h" +#include "test_main_nyxus.h" + +namespace Nyxus +{ + // "zero" origin (positive greybin_info), MATLAB-style: + // slope = n_levels / max, intercept = 1 + // bin = floor(slope * x + intercept), clipped to [1, n_levels] + static void test_bin_pixel_zero_origin() + { + int n_bins = 10; + PixIntens min_I = 50, max_I = 200, x = 100; + + // slope = 10 / 200 = 0.05, intercept = 1 - 0.05*0 = 1 + // bin = floor(0.05 * 100 + 1) = floor(6) = 6 + auto result = TextureFeature::bin_pixel(x, min_I, max_I, n_bins); + ASSERT_EQ(result, 6); + } + + // "min" origin (negative greybin_info), PyRadiomics-style: + // binWidth = (max - min) / binCount + // bin = floor((x - min) / binWidth) + 1 + static void test_bin_pixel_min_origin() + { + int n_bins = 10; + PixIntens min_I = 50, max_I = 200, x = 100; + + // binWidth = (200 - 50) / 10 = 15 + // bin = floor((100 - 50) / 15) + 1 = floor(3.33) + 1 = 4 + auto result = TextureFeature::bin_pixel(x, min_I, max_I, -n_bins); + ASSERT_EQ(result, 4); + } + + // End-to-end: GLCM_ASM computed with "zero" vs "min" binning origin + // must differ, and each must match its known ground truth. + static void test_glcm_binning_origin_divergence() + { + Fsettings s; + s.resize((int)NyxSetting::__COUNT__); + s[(int)NyxSetting::SOFTNAN].rval = 0.0; + s[(int)NyxSetting::TINY].rval = 0.0; + s[(int)NyxSetting::SINGLEROI].bval = false; + s[(int)NyxSetting::PIXELSIZEUM].rval = 100; + s[(int)NyxSetting::PIXELDISTANCE].ival = 5; + s[(int)NyxSetting::USEGPU].bval = false; + s[(int)NyxSetting::VERBOSLVL].ival = 0; + s[(int)NyxSetting::IBSI].bval = false; + s[(int)NyxSetting::GLCM_OFFSET].ival = 1; + GLCMFeature::symmetric_glcm = false; + GLCMFeature::angles = { 0, 45, 90, 135 }; + + int feature = int(Feature2D::GLCM_ASM); + + // --- "zero" origin (MATLAB) --- + s[(int)NyxSetting::GREYDEPTH].ival = 100; + s[(int)NyxSetting::GLCM_GREYDEPTH].ival = 100; + + LR roi_zero; + GLCMFeature f_zero; + load_masked_test_roi_data(roi_zero, ibsi_phantom_z1_intensity, ibsi_phantom_z1_mask, + sizeof(ibsi_phantom_z1_mask) / sizeof(NyxusPixel)); + ASSERT_NO_THROW(f_zero.calculate(roi_zero, s)); + roi_zero.initialize_fvals(); + f_zero.save_value(roi_zero.fvals); + double val_zero = roi_zero.fvals[feature][0]; + + // --- "min" origin (PyRadiomics) --- + s[(int)NyxSetting::GREYDEPTH].ival = -100; + s[(int)NyxSetting::GLCM_GREYDEPTH].ival = -100; + + LR roi_min; + GLCMFeature f_min; + load_masked_test_roi_data(roi_min, ibsi_phantom_z1_intensity, ibsi_phantom_z1_mask, + sizeof(ibsi_phantom_z1_mask) / sizeof(NyxusPixel)); + ASSERT_NO_THROW(f_min.calculate(roi_min, s)); + roi_min.initialize_fvals(); + f_min.save_value(roi_min.fvals); + double val_min = roi_min.fvals[feature][0]; + + // The two binning origins must produce different GLCM values + ASSERT_NE(val_zero, val_min) + << "zero_origin=" << val_zero << " min_origin=" << val_min; + + // Ground truth for GLCM_ASM angle-0 on ibsi_phantom_z1 at 100 bins + ASSERT_TRUE(agrees_gt(val_zero, 0.148438, 100.)) + << "zero_origin GLCM_ASM=" << val_zero; + ASSERT_TRUE(agrees_gt(val_min, 0.140625, 100.)) + << "min_origin GLCM_ASM=" << val_min; + } + + // End-to-end: NGTDM_COARSENESS computed with "zero" vs "min" binning origin + // must differ — exercises the non-GLCM texture feature path (bin_intensities). + static void test_ngtdm_binning_origin_divergence() + { + Fsettings s; + s.resize((int)NyxSetting::__COUNT__); + s[(int)NyxSetting::SOFTNAN].rval = 0.0; + s[(int)NyxSetting::TINY].rval = 0.0; + s[(int)NyxSetting::SINGLEROI].bval = false; + s[(int)NyxSetting::PIXELSIZEUM].rval = 100; + s[(int)NyxSetting::PIXELDISTANCE].ival = 5; + s[(int)NyxSetting::USEGPU].bval = false; + s[(int)NyxSetting::VERBOSLVL].ival = 0; + s[(int)NyxSetting::IBSI].bval = false; + + int feature = int(Feature2D::NGTDM_COARSENESS); + + // --- "zero" origin (MATLAB) --- + s[(int)NyxSetting::GREYDEPTH].ival = 100; + NGTDMFeature::n_levels = 0; // let it use GREYDEPTH + + LR roi_zero; + NGTDMFeature f_zero; + load_masked_test_roi_data(roi_zero, ibsi_phantom_z1_intensity, ibsi_phantom_z1_mask, + sizeof(ibsi_phantom_z1_mask) / sizeof(NyxusPixel)); + ASSERT_NO_THROW(f_zero.calculate(roi_zero, s)); + roi_zero.initialize_fvals(); + f_zero.save_value(roi_zero.fvals); + double val_zero = roi_zero.fvals[feature][0]; + + // --- "min" origin (PyRadiomics) --- + s[(int)NyxSetting::GREYDEPTH].ival = -100; + NGTDMFeature::n_levels = 0; + + LR roi_min; + NGTDMFeature f_min; + load_masked_test_roi_data(roi_min, ibsi_phantom_z1_intensity, ibsi_phantom_z1_mask, + sizeof(ibsi_phantom_z1_mask) / sizeof(NyxusPixel)); + ASSERT_NO_THROW(f_min.calculate(roi_min, s)); + roi_min.initialize_fvals(); + f_min.save_value(roi_min.fvals); + double val_min = roi_min.fvals[feature][0]; + + // The two binning origins must produce different NGTDM values + ASSERT_NE(val_zero, val_min) + << "zero_origin=" << val_zero << " min_origin=" << val_min; + } +} diff --git a/tests/test_fmaps.h b/tests/test_fmaps.h new file mode 100644 index 00000000..6d57f1da --- /dev/null +++ b/tests/test_fmaps.h @@ -0,0 +1,214 @@ +#pragma once + +#include "../src/nyx/roi_cache.h" +#include "../src/nyx/globals.h" +#include "../src/nyx/environment.h" +#include "test_main_nyxus.h" + +namespace Nyxus +{ + /// Helper: populates a test parent ROI (output param to avoid copy). + static void make_test_parent ( + LR& parent, int W, int H, + int xoff = 0, int yoff = 0, + PixIntens (*intensity_fn)(int, int) = nullptr) + { + for (int y = 0; y < H; y++) + for (int x = 0; x < W; x++) + { + int gx = xoff + x, gy = yoff + y; + PixIntens val = intensity_fn ? intensity_fn(x, y) : 100; + if (val == 0) + continue; + if (parent.aux_area == 0) + init_label_record_3(parent, gx, gy, val); + else + update_label_record_3(parent, gx, gy, val); + parent.raw_pixels.push_back(Pixel2(gx, gy, val)); + } + parent.make_nonanisotropic_aabb(); + parent.aux_image_matrix.allocate(W, H); + parent.aux_image_matrix.calculate_from_pixelcloud(parent.raw_pixels, parent.aabb); + } + + /// Helper: run generateChildRois and return the count. + struct FmapTestResult + { + int nChildren; + std::unordered_set childLabels; + std::unordered_map childRoiData; + std::unordered_map childToParentMap; + }; + + static FmapTestResult run_generate (const LR& parent, int kernel_size, int64_t startLabel = 1) + { + FmapTestResult r; + r.nChildren = generateChildRois(parent, kernel_size, r.childLabels, r.childRoiData, r.childToParentMap, startLabel); + return r; + } + + // Intensity functions for tests that need non-uniform values + static PixIntens intensity_varying (int x, int y) { return 10 * (y + 1) + (x + 1); } + static PixIntens intensity_offset (int x, int y) { return (PixIntens)(x + y + 1); } + static PixIntens intensity_checkerboard (int x, int y) { return ((x + y) % 2 == 0) ? 100 : 0; } + + /// Test correct child count and output container sizes + static void test_fmaps_child_roi_count() + { + LR parent(1); + make_test_parent(parent, 8, 8); + auto res = run_generate(parent, 3); + + int expected = (8 - 2) * (8 - 2); // 36 + if (res.nChildren != expected) + throw std::runtime_error("Expected " + std::to_string(expected) + " children, got " + std::to_string(res.nChildren)); + if ((int)res.childLabels.size() != res.nChildren) + throw std::runtime_error("childLabels size mismatch"); + if ((int)res.childRoiData.size() != res.nChildren) + throw std::runtime_error("childRoiData size mismatch"); + if ((int)res.childToParentMap.size() != res.nChildren) + throw std::runtime_error("childToParentMap size mismatch"); + } + + /// Test that each child AABB is kernel_size x kernel_size + static void test_fmaps_child_roi_dimensions() + { + LR parent(1); + make_test_parent(parent, 6, 6); + auto res = run_generate(parent, 3); + + for (auto& [label, child] : res.childRoiData) + { + int cw = child.aabb.get_width(), ch = child.aabb.get_height(); + if (cw != 3 || ch != 3) + throw std::runtime_error("Child " + std::to_string(label) + + " dimensions " + std::to_string(cw) + "x" + std::to_string(ch) + ", expected 3x3"); + } + } + + /// Test that all children map back to the correct parent label + static void test_fmaps_parent_mapping() + { + LR parent(42); + make_test_parent(parent, 5, 5); + auto res = run_generate(parent, 3); + + for (auto& [label, info] : res.childToParentMap) + { + if (info.parent_label != 42) + throw std::runtime_error("Child " + std::to_string(label) + + " has parent_label " + std::to_string(info.parent_label) + ", expected 42"); + } + } + + /// Test that child pixel intensities match the parent's image matrix + static void test_fmaps_child_pixel_values() + { + LR parent(1); + make_test_parent(parent, 5, 5, 0, 0, intensity_varying); + auto res = run_generate(parent, 3); + + // Find the child centered at (2,2) + for (auto& [label, info] : res.childToParentMap) + { + if (info.center_x == 2 && info.center_y == 2) + { + const LR& child = res.childRoiData.at(label); + if (child.aux_area != 9) + throw std::runtime_error("Center child should have 9 pixels, got " + std::to_string(child.aux_area)); + + for (const auto& px : child.raw_pixels) + { + PixIntens expected = 10 * (px.y + 1) + (px.x + 1); + if (px.inten != expected) + throw std::runtime_error("Pixel (" + std::to_string(px.x) + "," + std::to_string(px.y) + + ") intensity " + std::to_string(px.inten) + ", expected " + std::to_string(expected)); + } + return; + } + } + throw std::runtime_error("Could not find child ROI at center (2,2)"); + } + + /// Test coordinate handling with an offset parent (not at origin) + static void test_fmaps_nonorigin_parent() + { + const int xoff = 50, yoff = 100; + LR parent(1); + make_test_parent(parent, 6, 6, xoff, yoff, intensity_offset); + auto res = run_generate(parent, 3); + + int expected = (6 - 2) * (6 - 2); + if (res.nChildren != expected) + throw std::runtime_error("Nonorigin parent: expected " + std::to_string(expected) + + " children, got " + std::to_string(res.nChildren)); + + for (auto& [label, child] : res.childRoiData) + { + if (child.aabb.get_xmin() < xoff || child.aabb.get_ymin() < yoff) + throw std::runtime_error("Child AABB not offset correctly"); + } + + for (auto& [label, info] : res.childToParentMap) + { + if (info.center_x < xoff || info.center_y < yoff) + throw std::runtime_error("Child center not in global coordinates"); + } + } + + /// Test sparse mask: checkerboard pattern exercises the zero-center skip path + static void test_fmaps_sparse_mask() + { + LR parent(1); + make_test_parent(parent, 7, 7, 0, 0, intensity_checkerboard); + auto res = run_generate(parent, 3); + + // Count expected: kernel centers at (1..5, 1..5), only non-zero centers produce children + int expectedCount = 0; + for (int cy = 1; cy <= 5; cy++) + for (int cx = 1; cx <= 5; cx++) + if ((cx + cy) % 2 == 0) + expectedCount++; + + if (res.nChildren != expectedCount) + throw std::runtime_error("Sparse mask: expected " + std::to_string(expectedCount) + + " children, got " + std::to_string(res.nChildren)); + + for (auto& [label, child] : res.childRoiData) + { + if (child.aux_area > 9) + throw std::runtime_error("Sparse mask: child has too many pixels"); + if (child.aux_area == 0) + throw std::runtime_error("Sparse mask: child has zero pixels"); + } + } + + /// Test that parent smaller than kernel produces zero children + static void test_fmaps_parent_too_small() + { + LR parent(1); + make_test_parent(parent, 3, 3); + auto res = run_generate(parent, 5); + + if (res.nChildren != 0) + throw std::runtime_error("Parent too small: expected 0 children, got " + std::to_string(res.nChildren)); + } + + /// Test that startLabel offsets prevent label collisions across batches + static void test_fmaps_start_label_offset() + { + LR parent(1); + make_test_parent(parent, 5, 5); + auto r1 = run_generate(parent, 3, 1); + auto r2 = run_generate(parent, 3, 1 + r1.nChildren); + + for (auto lab : r2.childLabels) + if (r1.childLabels.count(lab) > 0) + throw std::runtime_error("Label collision: label " + std::to_string(lab) + " in both batches"); + + int minLabel2 = *std::min_element(r2.childLabels.begin(), r2.childLabels.end()); + if (minLabel2 != 1 + r1.nChildren) + throw std::runtime_error("Second batch min label should be " + std::to_string(1 + r1.nChildren) + + ", got " + std::to_string(minLabel2)); + } +} diff --git a/tests/test_fmaps_3d.h b/tests/test_fmaps_3d.h new file mode 100644 index 00000000..db0d3ca1 --- /dev/null +++ b/tests/test_fmaps_3d.h @@ -0,0 +1,216 @@ +#pragma once + +#include "../src/nyx/roi_cache.h" +#include "../src/nyx/globals.h" +#include "../src/nyx/environment.h" +#include "test_main_nyxus.h" + +namespace Nyxus +{ + /// Helper: populates a test 3D parent ROI (output param to avoid copy). + static void make_test_parent_3D ( + LR& parent, int W, int H, int D, + int xoff = 0, int yoff = 0, int zoff = 0, + PixIntens (*intensity_fn)(int, int, int) = nullptr) + { + for (int z = 0; z < D; z++) + for (int y = 0; y < H; y++) + for (int x = 0; x < W; x++) + { + int gx = xoff + x, gy = yoff + y, gz = zoff + z; + PixIntens val = intensity_fn ? intensity_fn(x, y, z) : 100; + if (val == 0) + continue; + if (parent.aux_area == 0) + init_label_record_3D(parent, gx, gy, gz, parent.label, val); + else + update_label_record_3D(parent, gx, gy, gz, parent.label, val); + parent.raw_pixels_3D.push_back(Pixel3(gx, gy, gz, val)); + } + parent.make_nonanisotropic_aabb(); + parent.aux_image_cube.allocate(W, H, D); + parent.aux_image_cube.calculate_from_pixelcloud(parent.raw_pixels_3D, parent.aabb); + } + + /// Helper: run generateChildRois_3D and return the count. + struct FmapTestResult3D + { + int nChildren; + std::unordered_set childLabels; + std::unordered_map childRoiData; + std::unordered_map childToParentMap; + }; + + static FmapTestResult3D run_generate_3D (const LR& parent, int kernel_size, int64_t startLabel = 1) + { + FmapTestResult3D r; + r.nChildren = generateChildRois_3D(parent, kernel_size, r.childLabels, r.childRoiData, r.childToParentMap, startLabel); + return r; + } + + // Intensity functions for 3D tests + static PixIntens intensity_varying_3D (int x, int y, int z) { return 100 * (z + 1) + 10 * (y + 1) + (x + 1); } + static PixIntens intensity_offset_3D (int x, int y, int z) { return (PixIntens)(x + y + z + 1); } + static PixIntens intensity_checkerboard_3D (int x, int y, int z) { return ((x + y + z) % 2 == 0) ? 100 : 0; } + + /// Test correct child count and output container sizes for 3D + static void test_fmaps3d_child_roi_count() + { + LR parent(1); + make_test_parent_3D(parent, 8, 8, 8); + auto res = run_generate_3D(parent, 3); + + int expected = (8 - 2) * (8 - 2) * (8 - 2); // 216 + if (res.nChildren != expected) + throw std::runtime_error("Expected " + std::to_string(expected) + " children, got " + std::to_string(res.nChildren)); + if ((int)res.childLabels.size() != res.nChildren) + throw std::runtime_error("childLabels size mismatch"); + if ((int)res.childRoiData.size() != res.nChildren) + throw std::runtime_error("childRoiData size mismatch"); + if ((int)res.childToParentMap.size() != res.nChildren) + throw std::runtime_error("childToParentMap size mismatch"); + } + + /// Test that each child AABB is kernel_size x kernel_size x kernel_size + static void test_fmaps3d_child_roi_dimensions() + { + LR parent(1); + make_test_parent_3D(parent, 6, 6, 6); + auto res = run_generate_3D(parent, 3); + + for (auto& [label, child] : res.childRoiData) + { + int cw = child.aabb.get_width(), ch = child.aabb.get_height(), cd = child.aabb.get_z_depth(); + if (cw != 3 || ch != 3 || cd != 3) + throw std::runtime_error("Child " + std::to_string(label) + + " dimensions " + std::to_string(cw) + "x" + std::to_string(ch) + "x" + std::to_string(cd) + ", expected 3x3x3"); + } + } + + /// Test that all children map back to the correct parent label + static void test_fmaps3d_parent_mapping() + { + LR parent(42); + make_test_parent_3D(parent, 5, 5, 5); + auto res = run_generate_3D(parent, 3); + + for (auto& [label, info] : res.childToParentMap) + { + if (info.parent_label != 42) + throw std::runtime_error("Child " + std::to_string(label) + + " has parent_label " + std::to_string(info.parent_label) + ", expected 42"); + } + } + + /// Test that child voxel intensities match the parent's image cube + static void test_fmaps3d_child_voxel_values() + { + LR parent(1); + make_test_parent_3D(parent, 5, 5, 5, 0, 0, 0, intensity_varying_3D); + auto res = run_generate_3D(parent, 3); + + // Find the child centered at (2,2,2) + for (auto& [label, info] : res.childToParentMap) + { + if (info.center_x == 2 && info.center_y == 2 && info.center_z == 2) + { + const LR& child = res.childRoiData.at(label); + if (child.aux_area != 27) + throw std::runtime_error("Center child should have 27 voxels, got " + std::to_string(child.aux_area)); + + for (const auto& px : child.raw_pixels_3D) + { + PixIntens expected = 100 * (px.z + 1) + 10 * (px.y + 1) + (px.x + 1); + if (px.inten != expected) + throw std::runtime_error("Voxel (" + std::to_string(px.x) + "," + std::to_string(px.y) + "," + std::to_string(px.z) + + ") intensity " + std::to_string(px.inten) + ", expected " + std::to_string(expected)); + } + return; + } + } + throw std::runtime_error("Could not find child ROI at center (2,2,2)"); + } + + /// Test coordinate handling with an offset parent (not at origin) + static void test_fmaps3d_nonorigin_parent() + { + const int xoff = 50, yoff = 100, zoff = 25; + LR parent(1); + make_test_parent_3D(parent, 6, 6, 6, xoff, yoff, zoff, intensity_offset_3D); + auto res = run_generate_3D(parent, 3); + + int expected = (6 - 2) * (6 - 2) * (6 - 2); // 64 + if (res.nChildren != expected) + throw std::runtime_error("Nonorigin parent: expected " + std::to_string(expected) + + " children, got " + std::to_string(res.nChildren)); + + for (auto& [label, child] : res.childRoiData) + { + if (child.aabb.get_xmin() < xoff || child.aabb.get_ymin() < yoff || child.aabb.get_zmin() < zoff) + throw std::runtime_error("Child AABB not offset correctly"); + } + + for (auto& [label, info] : res.childToParentMap) + { + if (info.center_x < xoff || info.center_y < yoff || info.center_z < zoff) + throw std::runtime_error("Child center not in global coordinates"); + } + } + + /// Test sparse mask: 3D checkerboard pattern exercises the zero-center skip path + static void test_fmaps3d_sparse_mask() + { + LR parent(1); + make_test_parent_3D(parent, 7, 7, 7, 0, 0, 0, intensity_checkerboard_3D); + auto res = run_generate_3D(parent, 3); + + // Count expected: kernel centers at (1..5, 1..5, 1..5), only non-zero centers produce children + int expectedCount = 0; + for (int cz = 1; cz <= 5; cz++) + for (int cy = 1; cy <= 5; cy++) + for (int cx = 1; cx <= 5; cx++) + if ((cx + cy + cz) % 2 == 0) + expectedCount++; + + if (res.nChildren != expectedCount) + throw std::runtime_error("Sparse mask: expected " + std::to_string(expectedCount) + + " children, got " + std::to_string(res.nChildren)); + + for (auto& [label, child] : res.childRoiData) + { + if (child.aux_area > 27) + throw std::runtime_error("Sparse mask: child has too many voxels"); + if (child.aux_area == 0) + throw std::runtime_error("Sparse mask: child has zero voxels"); + } + } + + /// Test that parent smaller than kernel produces zero children + static void test_fmaps3d_parent_too_small() + { + LR parent(1); + make_test_parent_3D(parent, 3, 3, 3); + auto res = run_generate_3D(parent, 5); + + if (res.nChildren != 0) + throw std::runtime_error("Parent too small: expected 0 children, got " + std::to_string(res.nChildren)); + } + + /// Test that startLabel offsets prevent label collisions across batches + static void test_fmaps3d_start_label_offset() + { + LR parent(1); + make_test_parent_3D(parent, 5, 5, 5); + auto r1 = run_generate_3D(parent, 3, 1); + auto r2 = run_generate_3D(parent, 3, 1 + r1.nChildren); + + for (auto lab : r2.childLabels) + if (r1.childLabels.count(lab) > 0) + throw std::runtime_error("Label collision: label " + std::to_string(lab) + " in both batches"); + + int minLabel2 = *std::min_element(r2.childLabels.begin(), r2.childLabels.end()); + if (minLabel2 != 1 + r1.nChildren) + throw std::runtime_error("Second batch min label should be " + std::to_string(1 + r1.nChildren) + + ", got " + std::to_string(minLabel2)); + } +}