diff --git a/.gitignore b/.gitignore index a9895a92..11763884 100644 --- a/.gitignore +++ b/.gitignore @@ -562,3 +562,13 @@ dmypy.json # Cython debug symbols cython_debug/ .DS_Store + +# CMake build artifacts (when cmake runs in the repo root) +CMakeCache.txt +CMakeFiles/ +cmake_install.cmake +Makefile + +# Test output files +*.arrow +*.parquet 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..9dbaa93a 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,19 @@ 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, +{'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..901bdeb9 100644 --- a/ci-utils/envs/conda_cpp.txt +++ b/ci-utils/envs/conda_cpp.txt @@ -9,5 +9,7 @@ xsimd >=13,<14 cmake dcmtk >=3.6.9 fmjpeg2koj >=1.0.3 -libarrow -libparquet +# Pin to <=23.0.0 — newer versions introduce breaking API changes +# that cause build failures. +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..3f4b354b 100644 --- a/docs/source/cmdline_and_examples.rst +++ b/docs/source/cmdline_and_examples.rst @@ -135,9 +135,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 +388,19 @@ 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, + {'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..8605dda6 100644 --- a/src/nyx/cli_option_constants.h +++ b/src/nyx/cli_option_constants.h @@ -61,6 +61,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..df5ed28d 100644 --- a/src/nyx/environment.cpp +++ b/src/nyx/environment.cpp @@ -474,6 +474,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 @@ -900,6 +903,27 @@ bool Environment::parse_cmdline(int argc, char** argv) ibsi_compliance = false; } + //==== Parse feature maps options + // --fmaps: enable/disable sliding-kernel feature extraction + 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"); + } + + // --fmaps_radius: kernel half-width (full kernel size = 2*radius+1) + 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; } diff --git a/src/nyx/environment.h b/src/nyx/environment.h index c089a898..cb239f27 100644 --- a/src/nyx/environment.h +++ b/src/nyx/environment.h @@ -148,6 +148,21 @@ class Environment: public BasicEnvironment ResultOptions resultOptions; std::tuple> parse_result_options_4cli (); + // Feature maps (fmaps) options — sliding-kernel feature extraction mode + bool fmaps_mode = false; ///< When true, features are computed per kernel position (spatial maps) + int fmaps_kernel_radius = 2; ///< Half-width of the kernel; full kernel side = 2*radius+1 + std::string raw_fmaps; ///< Raw CLI string for --fmaps flag (parsed in process_input) + std::string raw_fmaps_radius; ///< Raw CLI string for --fmaps_radius (parsed in process_input) + + /// @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, diff --git a/src/nyx/globals.h b/src/nyx/globals.h index 02c4896f..b1968478 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,92 @@ 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 — per-child metadata linking a child ROI back to its parent + /// @brief Tracks which parent ROI a child belongs to and the global image + /// coordinates of the kernel center that produced it. + struct FmapChildInfo + { + int parent_label; ///< Label of the parent ROI this child was generated from + int center_x; ///< Global x-coordinate of the kernel center + int center_y; ///< Global y-coordinate of the kernel center + int center_z = 0; ///< Global z-coordinate of the kernel center (unused in 2D) + }; + + /// @brief Shared helper: swaps child ROI data into env, reduces features, and outputs fmap arrays. + /// Used by both 2D and 3D fmaps workflows. + void reduceAndOutputChildRois_fmaps ( + Environment & env, + std::unordered_set childLabels, + std::unordered_map childRoiData, + const std::unordered_map & childToParentMap, + const std::string & intens_fpath, + const std::string & label_fpath, + int parentLab, + int parentXmin, int parentYmin, int parentZmin, + int parentW, int parentH, int parentD); + + /// @brief Generates 2D child ROIs by sliding a kernel across a parent ROI. + int generateChildRois ( + const LR & parent, + int kernel_size, + std::unordered_set & childLabels, + std::unordered_map & childRoiData, + std::unordered_map & childToParentMap, + int64_t startLabel); + /// @brief Generates 3D child ROIs by sliding a cubic kernel across a parent ROI. + 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_buffer (ResultsCache &results_cache, Environment &env, size_t t_index); + /// @brief Writes fmaps child ROI features into flat arrays for the Python buffer output path. + 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_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..698c2a1f 100644 --- a/src/nyx/helpers/helpers.h +++ b/src/nyx/helpers/helpers.h @@ -336,9 +336,9 @@ namespace Nyxus /// @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) { - if (disable_binning) + if (disable_binning) return i; - + double pi = ((double(i-min_i) / double(i_range) * double(n_levels))); 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..5b094b61 100644 --- a/src/nyx/main_nyxus.cpp +++ b/src/nyx/main_nyxus.cpp @@ -80,6 +80,13 @@ int main (int argc, char** argv) // Process the image data int errorCode = 0; + // Feature maps mode outputs spatial arrays and is only supported via the Python API + if (env.fmaps_mode) + { + std::cerr << "Error: feature maps (fmaps) mode is only supported via the Python API.\n"; + return 1; + } + if (env.singleROI) { errorCode = processDataset_2D_wholeslide( @@ -181,7 +188,15 @@ int main (int argc, char** argv) return 1; } - int errorCode = processDataset_3D_segmented( + int errorCode = 0; + // Feature maps mode outputs spatial arrays and is only supported via the Python API + if (env.fmaps_mode) + { + std::cerr << "Error: feature maps (fmaps) mode is only supported via the Python API.\n"; + return 1; + } + + errorCode = processDataset_3D_segmented( env, intensFiles, labelFiles, diff --git a/src/nyx/output_2_buffer.cpp b/src/nyx/output_2_buffer.cpp index 7c7f1348..dc4c8e32 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,94 @@ 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; } + /// @brief Assembles child ROI features into spatial arrays for one parent ROI. + /// Produces one FmapArrayResult with n_features arrays of shape (map_d, map_h, map_w). + /// Map dimensions are parent_dim - kernel_size + 1 (the valid convolution output size). + /// Positions with no valid child ROI remain NaN. + 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_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; + // Map dimensions = valid convolution output size + 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..2cca7bdc 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,20 @@ namespace Nyxus static std::mutex mutex1; + // Moved from function-local statics to namespace scope so they can be + // reset between workflow invocations (e.g. when processing multiple + // datasets in one Python session). Function-local statics persisted + // across calls, causing headers to be missing on subsequent runs. + // std::atomic ensures safe reset-before-use even without mutex protection. + static std::atomic mustRenderHeader_wholeslide{true}; + static std::atomic mustRenderHeader_segmented{true}; + + void reset_csv_header_state() + { + mustRenderHeader_wholeslide.store(true); + mustRenderHeader_segmented.store(true); + } + bool save_features_2_csv_wholeslide ( Environment & env, const LR & r, @@ -212,8 +335,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 +342,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 +356,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 +364,7 @@ namespace Nyxus // ********** Header - if (mustRenderHeader) + if (mustRenderHeader_wholeslide) { std::stringstream ssHead; @@ -263,7 +384,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 +403,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 +432,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 +439,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 +462,7 @@ namespace Nyxus std::vector> F = env.theFeatureSet.getEnabledFeatures(); // -- Header - if (mustRenderHeader) + if (mustRenderHeader_segmented) { std::stringstream ssHead; @@ -473,7 +482,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 +569,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 +589,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 +622,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)); } diff --git a/src/nyx/parallel.h b/src/nyx/parallel.h index 80dc73f3..76025ff5 100644 --- a/src/nyx/parallel.h +++ b/src/nyx/parallel.h @@ -39,6 +39,14 @@ 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. This is + // necessary for correctness in all callers, not just fmaps mode: + // without explicit .get(), futures are only awaited by their + // destructors, which can race with the caller's scope exit and + // cause use-after-free when referenced data (labels, LR map, + // settings) is destroyed while threads are still running. + 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..84cc308a 100644 --- a/src/nyx/python/new_bindings_py.cpp +++ b/src/nyx/python/new_bindings_py.cpp @@ -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( @@ -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()) @@ -142,6 +196,23 @@ 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 = {}, @@ -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."); + + 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."); + + 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())); @@ -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())); @@ -970,7 +1120,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 +1201,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 +1228,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..a7cf078f 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,16 @@ 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. """ def __init__( @@ -131,10 +142,11 @@ 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', } # Check for unexpected keyword arguments @@ -164,7 +176,9 @@ 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) + if neighbor_distance <= 0: raise ValueError("Neighbor distance must be greater than zero.") @@ -189,6 +203,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 +216,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 +228,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 +341,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 +363,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 +478,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 +504,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 +563,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 +584,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)) @@ -698,11 +734,11 @@ def set_environment_params (self, **params): 'max_intensity', 'ram_limit', ] - + for key in params: if key not in valid_params: raise ValueError(f'Invalid environment parameter {key}. Value parameters are {params}') - + features = params.get('features', []) neighbor_distance = params.get ('neighbor_distance', -1) pixels_per_micron = params.get ('pixels_per_micron', -1) @@ -714,10 +750,10 @@ def set_environment_params (self, **params): min_intensity = params.get('min_intensity', -1) max_intensity = params.get('max_intensity', -1) ram_limit = params.get('ram_limit', -1) - + set_environment_params_imp (id(self), - features, - neighbor_distance, + features, + neighbor_distance, pixels_per_micron, coarse_gray_depth, n_reduce_threads, @@ -775,19 +811,31 @@ 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 + else: if (key not in available_environment_params): raise ValueError ("Invalid parameter: ", key) else: environment_params[key] = value - + + # 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 +992,16 @@ 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. """ def __init__( @@ -953,20 +1011,21 @@ 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', } # 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,7 +1041,9 @@ 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) + if neighbor_distance <= 0: raise ValueError("Neighbor distance must be greater than zero.") @@ -994,7 +1055,7 @@ def __init__( 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 +1078,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 +1101,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 +1203,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 +1272,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 +1293,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)) @@ -1257,13 +1333,13 @@ def set_environment_params (self, **params): 'verbose', 'dynamic_range', 'min_intensity', - 'max_intensity' + 'max_intensity', ] - + for key in params: if key not in valid_params: raise ValueError(f'Invalid environment parameter {key}. Value parameters are {params}') - + features = params.get('features', []) neighbor_distance = params.get ('neighbor_distance', -1) pixels_per_micron = params.get ('pixels_per_micron', -1) @@ -1275,10 +1351,10 @@ def set_environment_params (self, **params): min_intensity = params.get('min_intensity', -1) max_intensity = params.get('max_intensity', -1) ram_limit = -1 # no limit - + set_environment_params_imp (id(self), - features, - neighbor_distance, + features, + neighbor_distance, pixels_per_micron, coarse_gray_depth, n_reduce_threads, @@ -1288,12 +1364,12 @@ def set_environment_params (self, **params): max_intensity, ram_limit, verb_lvl) - + def set_params(self, **params): """Sets parameters of the Nyxus class Keyword args: - + * features: List[str], * neighbor_distance * pixels_per_micron @@ -1326,17 +1402,16 @@ def set_params(self, **params): for key, value in params.items(): - + if (key == "ibsi"): set_if_ibsi_imp (id(self), value) - + else: if (key not in available_environment_params): raise ValueError(f"Invalid parameter {key}.") else: environment_params[key] = value - - + if (len(environment_params) > 0): self.set_environment_params(**environment_params) @@ -1550,7 +1625,7 @@ def __init__( features, neighbor_distance, pixels_per_micron, - coarse_gray_depth, + coarse_gray_depth, n_feature_calc_threads, using_gpu, ibsi, @@ -1562,8 +1637,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'] @@ -1991,13 +2068,13 @@ def set_environment_params (self, **params): 'verbose', 'dynamic_range', 'min_intensity', - 'max_intensity' + 'max_intensity', ] - + for key in params: if key not in valid_params: raise ValueError(f'Invalid environment parameter {key}. Value parameters are {params}') - + features = params.get('features', []) neighbor_distance = params.get ('neighbor_distance', -1) pixels_per_micron = params.get ('pixels_per_micron', -1) @@ -2009,10 +2086,10 @@ def set_environment_params (self, **params): min_intensity = params.get('min_intensity', -1) max_intensity = params.get('max_intensity', -1) ram_limit = -1 # no limit - + set_environment_params_imp (id(self), features, - neighbor_distance, + neighbor_distance, pixels_per_micron, coarse_gray_depth, n_reduce_threads, @@ -2022,7 +2099,7 @@ def set_environment_params (self, **params): max_intensity, ram_limit, verb_lvl) - + def set_params(self, **params): """Sets parameters of the Nyxus class diff --git a/src/nyx/results_cache.h b/src/nyx/results_cache.h index 75973105..e78a7e81 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,8 @@ class ResultsCache for (auto c : cols) add_to_header(c); } - void add_to_header(std::string& col) + // const ref: callers may pass temporaries or string expressions + void add_to_header(const std::string& col) { headerBuf_.push_back(col); } @@ -35,9 +56,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..57cc02ff --- /dev/null +++ b/src/nyx/workflow_2d_fmaps.cpp @@ -0,0 +1,415 @@ +/// @file workflow_2d_fmaps.cpp +/// @brief 2D feature maps (fmaps) workflow: slides a kernel across each parent ROI to +/// produce spatially-resolved feature maps. For every valid kernel position a +/// "child" ROI is created, its features are computed, and the results are stored +/// either as flat arrays (Python buffer mode) or as a CSV with per-pixel rows. + +#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 Reduces child ROIs to feature values and outputs results as spatial arrays. + /// Temporarily swaps child ROI data into the environment (via RAII guard), + /// computes features for all children, and writes results to fmap arrays. + void reduceAndOutputChildRois_fmaps ( + Environment & env, + std::unordered_set childLabels, + std::unordered_map childRoiData, + const std::unordered_map & childToParentMap, + const std::string & intens_fpath, + const std::string & label_fpath, + int parentLab, + int parentXmin, int parentYmin, int parentZmin, + int parentW, int parentH, int parentD) + { + // RAII guard swaps env's ROI data with child data and restores on scope exit + 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); + + // Output child ROI features as spatial arrays + save_features_2_fmap_arrays ( + env.theResultsCache, + env, + intens_fpath, + label_fpath, + parentLab, + parentXmin, parentYmin, parentZmin, + parentW, parentH, parentD, + env.fmaps_kernel_size(), + env.uniqueLabels, + env.roiData, + childToParentMap); + } + + /// @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 and output results + reduceAndOutputChildRois_fmaps ( + env, + std::move(childLabels), std::move(childRoiData), childToParentMap, + intens_fpath, label_fpath, + parentLab, + parentXmin, parentYmin, /*parentZmin=*/ 0, + parentW, parentH, /*parentD=*/ 1); + } + + // 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; + } + + /// @brief Entry point for the 2D feature maps workflow. + /// Prescans all image pairs, then processes each pair by generating child ROIs + /// per parent and computing features at each kernel position. + /// @return 0 on success, nonzero on error. + 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..0996d95d --- /dev/null +++ b/src/nyx/workflow_3d_fmaps.cpp @@ -0,0 +1,362 @@ +/// @file workflow_3d_fmaps.cpp +/// @brief 3D feature maps (fmaps) workflow: slides a cubic kernel across each parent ROI +/// to produce volumetric feature maps. Analogous to workflow_2d_fmaps.cpp but +/// operates on 3D image cubes and supports time-frame iteration. + +#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. + /// Gathers parent ROI metrics, loads voxels, generates 3D child ROIs per kernel, + /// and outputs features for each child. + /// @param globalChildLabel [in/out] Running child label counter, persists across file pairs. + 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 and output results + reduceAndOutputChildRois_fmaps ( + env, + std::move(childLabels), std::move(childRoiData), childToParentMap, + intens_fpath, label_fpath, + parentLab, + parentXmin, parentYmin, parentZmin, + parentW, parentH, parentD); + } + + // 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; + } + + /// @brief Entry point for the 3D feature maps workflow. + /// Prescans all image pairs, then iterates time frames and file pairs, + /// generating 3D child ROIs and computing volumetric features at each kernel position. + /// @return 0 on success, nonzero on error. + 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..c65b3814 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, /*parent_zmin=*/ 0, + parentW, parentH, /*parent_d=*/ 1, + 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." }; + + 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..b88f2094 100644 --- a/tests/test_all.cc +++ b/tests/test_all.cc @@ -31,12 +31,83 @@ #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" #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()); +} + + //***** 2D contour and multicontour ***** TEST(TEST_NYXUS, TEST_CONTOUR_MULTI_1) { 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..088ecae0 --- /dev/null +++ b/tests/test_fmaps_3d.h @@ -0,0 +1,209 @@ +#pragma once + +#include "../src/nyx/roi_cache.h" +#include "../src/nyx/globals.h" +#include "../src/nyx/environment.h" +#include "test_main_nyxus.h" +#include "test_fmaps.h" // reuse FmapTestResult + +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. + static FmapTestResult run_generate_3D (const LR& parent, int kernel_size, int64_t startLabel = 1) + { + FmapTestResult 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)); + } +}