From 61064d701fea46f4090f589a12c5f657d63c8ea1 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Thu, 12 Mar 2026 15:19:43 -0400 Subject: [PATCH 01/10] Fix: unify wholeslide logic and honor single_roi=True in featurize_files (#327) Make theEnvironment.singleROI the single source of truth for wholeslide mode across all entry points and scan sites. - dirs_and_files.cpp: set singleROI=true when wholeslide detected from directory structure, so downstream scan sites see the flag - new_bindings_py.cpp: fill labelFiles with empty strings when single_roi=True; fix validation to allow empty seg list in that case; only check seg file existence when path is non-empty - nyxus.py: allow mask_files=None when single_roi=True; pass [] instead of None to the C++ binding - slideprops.cpp, phase1.cpp, phase2.cpp: replace independent spL==nullptr / fname_seg.empty() checks with theEnvironment.singleROI (3 sites in phase2.cpp) Co-Authored-By: Claude Sonnet 4.6 --- src/nyx/dirs_and_files.cpp | 1 + src/nyx/phase1.cpp | 802 +++++++------- src/nyx/phase2_2d.cpp | 1524 +++++++++++++------------- src/nyx/phase2_3d.cpp | 1640 ++++++++++++++-------------- src/nyx/python/new_bindings_py.cpp | 34 +- src/nyx/python/nyxus/nyxus.py | 14 +- src/nyx/slideprops.cpp | 987 +++++++++-------- 7 files changed, 2501 insertions(+), 2501 deletions(-) diff --git a/src/nyx/dirs_and_files.cpp b/src/nyx/dirs_and_files.cpp index f5b5cce1..a9145833 100644 --- a/src/nyx/dirs_and_files.cpp +++ b/src/nyx/dirs_and_files.cpp @@ -78,6 +78,7 @@ namespace Nyxus if (wholeslide) { + theEnvironment.singleROI = true; // populate with empty mask file names labelFiles.insert (labelFiles.begin(), intensFiles.size(), ""); } diff --git a/src/nyx/phase1.cpp b/src/nyx/phase1.cpp index 1bc94940..31bc8c71 100644 --- a/src/nyx/phase1.cpp +++ b/src/nyx/phase1.cpp @@ -1,402 +1,402 @@ -#include -#include -#include -#include -#include - -#ifdef WITH_PYTHON_H - #include - #include - #include - namespace py = pybind11; -#endif - -#include "environment.h" -#include "globals.h" -#include "helpers/timing.h" - -namespace Nyxus -{ - // - // segmented 2D case - // - bool gatherRoisMetrics (int sidx, const std::string & intens_fpath, const std::string & label_fpath, Environment & env, ImageLoader & L) - { - // Reset per-image counters and extrema - // -- disabling this due to new prescan functionality--> LR::reset_global_stats(); - - int lvl = 0, // Pyramid level - lyr = 0; // Layer - - // Read the tiff. The image loader is put in the open state in processDataset() - size_t nth = L.get_num_tiles_hor(), - ntv = L.get_num_tiles_vert(), - fw = L.get_tile_width(), - th = L.get_tile_height(), - tw = L.get_tile_width(), - tileSize = L.get_tile_size(), - fullwidth = L.get_full_width(), - fullheight = L.get_full_height(); - - int cnt = 1; - for (unsigned int row = 0; row < nth; row++) - for (unsigned int col = 0; col < ntv; col++) - { - // Fetch the tile - bool ok = L.load_tile(row, col); - if (!ok) - { - std::string erm = "Error fetching tile row:" + std::to_string(row) + " col:" + std::to_string(col) + " from I:" + intens_fpath + " M:" + label_fpath; - #ifdef WITH_PYTHON_H - throw erm; - #endif - std::cerr << erm << "\n"; - return false; - } - - // Get ahold of tile's pixel buffer - const std::vector& dataI = L.get_int_tile_buffer(); - const std::shared_ptr>& spL = L.get_seg_tile_sptr(); - bool wholeslide = spL == nullptr; - - // Iterate pixels - for (size_t i = 0; i < tileSize; i++) - { - // mask label if not in the wholeslide mode - PixIntens label = 1; - if (!wholeslide) - label = (*spL)[i]; - - // Skip non-mask pixels - if (! label) - continue; - - int y = row * th + i / tw, - x = col * tw + i % tw; - - // Skip tile buffer pixels beyond the image's bounds - if (x >= fullwidth || y >= fullheight) - continue; - - // Update pixel's ROI metrics - feed_pixel_2_metrics (env.uniqueLabels, env.roiData, x, y, dataI[i], label, sidx); - } - -#ifdef WITH_PYTHON_H - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); -#endif - - // Show progress info - VERBOSLVL2 (env.get_verbosity_level(), - if (cnt++ % 4 == 0) - std::cout << "\t" << int((row * nth + col) * 100 / float(nth * ntv) * 100) / 100. << "%\t" << env.uniqueLabels.size() << " ROIs" << "\n"; - ); - } - - // fix ROIs' AABBs with respect to anisotropy - if (env.anisoOptions.customized() == false) - { - for (auto& rd : env.roiData) - { - LR& r = rd.second; - r.make_nonanisotropic_aabb (); - } - } - else - { - double ax = env.anisoOptions.get_aniso_x(), - ay = env.anisoOptions.get_aniso_y(); - - for (auto& rd : env.roiData) - { - LR& r = rd.second; - r.make_anisotropic_aabb (ax, ay); - } - } - - return true; - } - - // - // segmented 2.5D case (aka layoutA, collections of 2D slice images e.g. blah_z1_blah.ome.tif, blah_z2_blah.ome.tif, ..., blah_z500_blah.ome.tif) - // prerequisite: 'env.theImLoader' needs to be pre-opened ! - // - bool gatherRoisMetrics_25D (Environment & env, size_t sidx, const std::string& intens_fpath, const std::string& mask_fpath, const std::vector& z_indices) - { - for (size_t z=0; z= fullwidth || y >= fullheight) - continue; - - // Collapse all the labels to one if single-ROI mde is requested - if (env.singleROI) - label = 1; - - // Update pixel's ROI metrics - feed_pixel_2_metrics_3D (env.uniqueLabels, env.roiData, x, y, z, dataI[i], label, sidx); // Updates 'uniqueLabels' and 'roiData' - } - - #ifdef WITH_PYTHON_H - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); - #endif - - // Show stayalive progress info - VERBOSLVL2 (env.get_verbosity_level(), - if (cnt++ % 4 == 0) - std::cout << "\t" << int((row * nth + col) * 100 / float(nth * ntv) * 100) / 100. << "%\t" << env.uniqueLabels.size() << " ROIs" << "\n"; - ); - } - - env.theImLoader.close(); - } - - // fix ROIs' AABBs with respect to anisotropy - if (env.anisoOptions.customized() == false) - { - for (auto& rd : env.roiData) - { - LR& r = rd.second; - r.make_nonanisotropic_aabb(); - } - } - else - { - double ax = env.anisoOptions.get_aniso_x(), - ay = env.anisoOptions.get_aniso_y(), - az = env.anisoOptions.get_aniso_z(); - - for (auto& rd : env.roiData) - { - LR& r = rd.second; - r.make_anisotropic_aabb (ax, ay, az); - } - } - - return true; - } - - // - // segmented 3D case (true volumetric images e.g. .nii, .nii.gz, .dcm, etc) - // prerequisite: 'env.theImLoader' needs to be pre-opened ! - // - bool gatherRoisMetrics_3D (Environment& env, size_t sidx, const std::string& intens_fpath, const std::string& mask_fpath, size_t t_index) - { - SlideProps & sprp = env.dataset.dataset_props [sidx]; - if (! env.theImLoader.open(sprp, env.fpimageOptions)) - { - std::cerr << "Error opening a file pair with ImageLoader. Terminating\n"; - return false; - } - - // Read the tiff. The image loader is put in the open state in processDataset() - size_t - w = env.theImLoader.get_full_width(), - h = env.theImLoader.get_full_height(), - d = env.theImLoader.get_full_depth(), - sliceSize = w * h, - timeFrameSize = sliceSize * d, - timeI = env.theImLoader.get_inten_time(), - timeM = env.theImLoader.get_mask_time(), // may be ==0 for WSI - nVoxI = timeFrameSize * timeI, - nVoxM = timeFrameSize * timeM; - - // is this intensity-mask pair's shape supported? - if (nVoxI < nVoxM) - { - std::string erm = "Error: unsupported shape - intensity file: " + std::to_string(nVoxI) + ", mask file: " + std::to_string(nVoxM); -#ifdef WITH_PYTHON_H - throw erm; -#endif - std::cerr << erm << "\n"; - return false; - } - - int cnt = 1; - - // Fetch a tile - bool ok = env.theImLoader.load_tile (0/*row*/ , 0/*col*/); - if (!ok) - { - std::string erm = "Error fetching tile (0,0) from I:" + intens_fpath + " M:" + mask_fpath; - #ifdef WITH_PYTHON_H - throw erm; - #endif - std::cerr << erm << "\n"; - return false; - } - - // Get ahold of tile's pixel buffer - auto dataI = env.theImLoader.get_int_tile_buffer(), - dataL = env.theImLoader.get_seg_tile_buffer(); - - size_t baseI, baseM; - if (nVoxI == nVoxM) - { - baseM = baseI = t_index * timeFrameSize; - } - else // nVoxI > nVoxM - { - baseM = 0; - baseI = t_index * timeFrameSize; - } - - // Iterate voxels - for (size_t i=0; i= w || y >= h || z >= d) - continue; - - // Collapse all the labels to one if single-ROI mde is requested - if (env.singleROI) - label = 1; - - // Update pixel's ROI metrics - feed_pixel_2_metrics_3D (env.uniqueLabels, env.roiData, x, y, z, dataI[j], label, sidx); - } - -#ifdef WITH_PYTHON_H - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); -#endif - - env.theImLoader.close(); - - // fix ROIs' AABBs with respect to anisotropy - if (env.anisoOptions.customized() == false) - { - for (auto& rd : env.roiData) - { - LR& r = rd.second; - r.make_nonanisotropic_aabb(); - } - } - else - { - double ax = env.anisoOptions.get_aniso_x(), - ay = env.anisoOptions.get_aniso_y(), - az = env.anisoOptions.get_aniso_z(); - - for (auto& rd : env.roiData) - { - LR& r = rd.second; - r.make_anisotropic_aabb(ax, ay, az); - } - } - - return true; - } - -#ifdef WITH_PYTHON_H - - // - // segmented 2D case - // - bool gatherRoisMetricsInMemory (Environment & env, const py::array_t& intens_images, const py::array_t& label_images, int pair_index) - { - VERBOSLVL4 (env.get_verbosity_level(), std::cout << "gatherRoisMetricsInMemory (pair_index=" << pair_index << ") \n"); - - auto rI = intens_images.unchecked<3>(); - auto rL = label_images.unchecked<3>(); - - size_t w = rI.shape(2); - size_t h = rI.shape(1); - - for (size_t col = 0; col < w; col++) - for (size_t row = 0; row < h; row++) - { - // Skip non-mask pixels - auto label = rL (pair_index, row, col); - if (!label) - continue; - - // Collapse all the labels to one if single-ROI mde is requested - if (env.singleROI) - label = 1; - - // Update pixel's ROI metrics - auto inten = rI (pair_index, row, col); - feed_pixel_2_metrics (env.uniqueLabels, env.roiData, col, row, inten, label, pair_index); // Updates 'uniqueLabels' and 'roiData' - - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); - } - - return true; - } - -#endif +#include +#include +#include +#include +#include + +#ifdef WITH_PYTHON_H + #include + #include + #include + namespace py = pybind11; +#endif + +#include "environment.h" +#include "globals.h" +#include "helpers/timing.h" + +namespace Nyxus +{ + // + // segmented 2D case + // + bool gatherRoisMetrics (int sidx, const std::string & intens_fpath, const std::string & label_fpath, Environment & env, ImageLoader & L) + { + // Reset per-image counters and extrema + // -- disabling this due to new prescan functionality--> LR::reset_global_stats(); + + int lvl = 0, // Pyramid level + lyr = 0; // Layer + + // Read the tiff. The image loader is put in the open state in processDataset() + size_t nth = L.get_num_tiles_hor(), + ntv = L.get_num_tiles_vert(), + fw = L.get_tile_width(), + th = L.get_tile_height(), + tw = L.get_tile_width(), + tileSize = L.get_tile_size(), + fullwidth = L.get_full_width(), + fullheight = L.get_full_height(); + + int cnt = 1; + for (unsigned int row = 0; row < nth; row++) + for (unsigned int col = 0; col < ntv; col++) + { + // Fetch the tile + bool ok = L.load_tile(row, col); + if (!ok) + { + std::string erm = "Error fetching tile row:" + std::to_string(row) + " col:" + std::to_string(col) + " from I:" + intens_fpath + " M:" + label_fpath; + #ifdef WITH_PYTHON_H + throw erm; + #endif + std::cerr << erm << "\n"; + return false; + } + + // Get ahold of tile's pixel buffer + const std::vector& dataI = L.get_int_tile_buffer(); + const std::shared_ptr>& spL = L.get_seg_tile_sptr(); + bool wholeslide = env.singleROI; + + // Iterate pixels + for (size_t i = 0; i < tileSize; i++) + { + // mask label if not in the wholeslide mode + PixIntens label = 1; + if (!wholeslide) + label = (*spL)[i]; + + // Skip non-mask pixels + if (! label) + continue; + + int y = row * th + i / tw, + x = col * tw + i % tw; + + // Skip tile buffer pixels beyond the image's bounds + if (x >= fullwidth || y >= fullheight) + continue; + + // Update pixel's ROI metrics + feed_pixel_2_metrics (env.uniqueLabels, env.roiData, x, y, dataI[i], label, sidx); + } + +#ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); +#endif + + // Show progress info + VERBOSLVL2 (env.get_verbosity_level(), + if (cnt++ % 4 == 0) + std::cout << "\t" << int((row * nth + col) * 100 / float(nth * ntv) * 100) / 100. << "%\t" << env.uniqueLabels.size() << " ROIs" << "\n"; + ); + } + + // fix ROIs' AABBs with respect to anisotropy + if (env.anisoOptions.customized() == false) + { + for (auto& rd : env.roiData) + { + LR& r = rd.second; + r.make_nonanisotropic_aabb (); + } + } + else + { + double ax = env.anisoOptions.get_aniso_x(), + ay = env.anisoOptions.get_aniso_y(); + + for (auto& rd : env.roiData) + { + LR& r = rd.second; + r.make_anisotropic_aabb (ax, ay); + } + } + + return true; + } + + // + // segmented 2.5D case (aka layoutA, collections of 2D slice images e.g. blah_z1_blah.ome.tif, blah_z2_blah.ome.tif, ..., blah_z500_blah.ome.tif) + // prerequisite: 'env.theImLoader' needs to be pre-opened ! + // + bool gatherRoisMetrics_25D (Environment & env, size_t sidx, const std::string& intens_fpath, const std::string& mask_fpath, const std::vector& z_indices) + { + for (size_t z=0; z= fullwidth || y >= fullheight) + continue; + + // Collapse all the labels to one if single-ROI mde is requested + if (env.singleROI) + label = 1; + + // Update pixel's ROI metrics + feed_pixel_2_metrics_3D (env.uniqueLabels, env.roiData, x, y, z, dataI[i], label, sidx); // Updates 'uniqueLabels' and 'roiData' + } + + #ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); + #endif + + // Show stayalive progress info + VERBOSLVL2 (env.get_verbosity_level(), + if (cnt++ % 4 == 0) + std::cout << "\t" << int((row * nth + col) * 100 / float(nth * ntv) * 100) / 100. << "%\t" << env.uniqueLabels.size() << " ROIs" << "\n"; + ); + } + + env.theImLoader.close(); + } + + // fix ROIs' AABBs with respect to anisotropy + if (env.anisoOptions.customized() == false) + { + for (auto& rd : env.roiData) + { + LR& r = rd.second; + r.make_nonanisotropic_aabb(); + } + } + else + { + double ax = env.anisoOptions.get_aniso_x(), + ay = env.anisoOptions.get_aniso_y(), + az = env.anisoOptions.get_aniso_z(); + + for (auto& rd : env.roiData) + { + LR& r = rd.second; + r.make_anisotropic_aabb (ax, ay, az); + } + } + + return true; + } + + // + // segmented 3D case (true volumetric images e.g. .nii, .nii.gz, .dcm, etc) + // prerequisite: 'env.theImLoader' needs to be pre-opened ! + // + bool gatherRoisMetrics_3D (Environment& env, size_t sidx, const std::string& intens_fpath, const std::string& mask_fpath, size_t t_index) + { + SlideProps & sprp = env.dataset.dataset_props [sidx]; + if (! env.theImLoader.open(sprp, env.fpimageOptions)) + { + std::cerr << "Error opening a file pair with ImageLoader. Terminating\n"; + return false; + } + + // Read the tiff. The image loader is put in the open state in processDataset() + size_t + w = env.theImLoader.get_full_width(), + h = env.theImLoader.get_full_height(), + d = env.theImLoader.get_full_depth(), + sliceSize = w * h, + timeFrameSize = sliceSize * d, + timeI = env.theImLoader.get_inten_time(), + timeM = env.theImLoader.get_mask_time(), // may be ==0 for WSI + nVoxI = timeFrameSize * timeI, + nVoxM = timeFrameSize * timeM; + + // is this intensity-mask pair's shape supported? + if (nVoxI < nVoxM) + { + std::string erm = "Error: unsupported shape - intensity file: " + std::to_string(nVoxI) + ", mask file: " + std::to_string(nVoxM); +#ifdef WITH_PYTHON_H + throw erm; +#endif + std::cerr << erm << "\n"; + return false; + } + + int cnt = 1; + + // Fetch a tile + bool ok = env.theImLoader.load_tile (0/*row*/ , 0/*col*/); + if (!ok) + { + std::string erm = "Error fetching tile (0,0) from I:" + intens_fpath + " M:" + mask_fpath; + #ifdef WITH_PYTHON_H + throw erm; + #endif + std::cerr << erm << "\n"; + return false; + } + + // Get ahold of tile's pixel buffer + auto dataI = env.theImLoader.get_int_tile_buffer(), + dataL = env.theImLoader.get_seg_tile_buffer(); + + size_t baseI, baseM; + if (nVoxI == nVoxM) + { + baseM = baseI = t_index * timeFrameSize; + } + else // nVoxI > nVoxM + { + baseM = 0; + baseI = t_index * timeFrameSize; + } + + // Iterate voxels + for (size_t i=0; i= w || y >= h || z >= d) + continue; + + // Collapse all the labels to one if single-ROI mde is requested + if (env.singleROI) + label = 1; + + // Update pixel's ROI metrics + feed_pixel_2_metrics_3D (env.uniqueLabels, env.roiData, x, y, z, dataI[j], label, sidx); + } + +#ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); +#endif + + env.theImLoader.close(); + + // fix ROIs' AABBs with respect to anisotropy + if (env.anisoOptions.customized() == false) + { + for (auto& rd : env.roiData) + { + LR& r = rd.second; + r.make_nonanisotropic_aabb(); + } + } + else + { + double ax = env.anisoOptions.get_aniso_x(), + ay = env.anisoOptions.get_aniso_y(), + az = env.anisoOptions.get_aniso_z(); + + for (auto& rd : env.roiData) + { + LR& r = rd.second; + r.make_anisotropic_aabb(ax, ay, az); + } + } + + return true; + } + +#ifdef WITH_PYTHON_H + + // + // segmented 2D case + // + bool gatherRoisMetricsInMemory (Environment & env, const py::array_t& intens_images, const py::array_t& label_images, int pair_index) + { + VERBOSLVL4 (env.get_verbosity_level(), std::cout << "gatherRoisMetricsInMemory (pair_index=" << pair_index << ") \n"); + + auto rI = intens_images.unchecked<3>(); + auto rL = label_images.unchecked<3>(); + + size_t w = rI.shape(2); + size_t h = rI.shape(1); + + for (size_t col = 0; col < w; col++) + for (size_t row = 0; row < h; row++) + { + // Skip non-mask pixels + auto label = rL (pair_index, row, col); + if (!label) + continue; + + // Collapse all the labels to one if single-ROI mde is requested + if (env.singleROI) + label = 1; + + // Update pixel's ROI metrics + auto inten = rI (pair_index, row, col); + feed_pixel_2_metrics (env.uniqueLabels, env.roiData, col, row, inten, label, pair_index); // Updates 'uniqueLabels' and 'roiData' + + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); + } + + return true; + } + +#endif } \ No newline at end of file diff --git a/src/nyx/phase2_2d.cpp b/src/nyx/phase2_2d.cpp index 8f7f7620..6cded7c0 100644 --- a/src/nyx/phase2_2d.cpp +++ b/src/nyx/phase2_2d.cpp @@ -1,762 +1,762 @@ -#include -#include -#include -#include -#include -#include -#include - -#ifdef WITH_PYTHON_H - #include -#endif - -#include "environment.h" -#include "globals.h" -#include "helpers/timing.h" - -namespace Nyxus -{ - - #define disable_DUMP_ALL_ROI - #ifdef DUMP_ALL_ROI - void dump_all_roi() - { - std::string fpath = theEnvironment.output_dir + "/all_roi.txt"; - std::cout << "Dumping all the ROIs to " << fpath << " ...\n"; - - std::ofstream f(fpath); - - for (auto lab : uniqueLabels) - { - auto& r = roiData[lab]; - std::cout << "Dumping ROI " << lab << "\n"; - - r.aux_image_matrix.print(f); - - f << "ROI " << lab << ": \n" - << "xmin = " << r.aabb.get_xmin() << "; \n" - << "width=" << r.aabb.get_width() << "; \n" - << "ymin=" << r.aabb.get_ymin() << "; \n" - << "height=" << r.aabb.get_height() << "; \n" - << "area=" << r.aux_area << "; \n"; - - // C++ constant: - f << "// C:\n" - << "struct NyxusPixel {\n" - << "\tsize_t x, y; \n" - << "\tunsigned int intensity; \n" - << "}; \n" - << "NyxusPixel testData[] = {\n"; - for (auto i=0; i 0 && i % 4 == 0) - f << "\n"; - } - f << "}; \n"; - - // Matlab constant: - f << "// MATLAB:\n" - << "%==== begin \n"; - f << "pixelCloud = [ \n"; - for (auto i = 0; i < r.raw_pixels.size(); i++) - { - auto& px = r.raw_pixels[i]; - f << px.inten << "; % [" << i << "] \n"; - } - f << "]; \n"; - - f << "testData = zeros(" << r.aabb.get_height() << "," << r.aabb.get_width() << ");\n"; - for (auto i = 0; i < r.raw_pixels.size(); i++) - { - auto& px = r.raw_pixels[i]; - f << "testData(" << (px.y - r.aabb.get_ymin() + 1) << "," << (px.x - r.aabb.get_xmin() + 1) << ")=" << px.inten << "; "; // +1 due to 1-based nature of Matlab - if (i > 0 && i % 4 == 0) - f << "\n"; - } - f << "\n"; - f << "testVecZ = reshape(testData, 1, []); \n"; - f << "testVecNZ = nonzeros(testData); \n"; - f << "[mean(testVecNZ) mean(testVecZ) mean2(testData)] \n"; - f << "%==== end \n"; - } - - f.flush(); - } - #endif - - bool scanTrivialRois ( - const std::vector& batch_labels, - const std::string& intens_fpath, - const std::string& label_fpath, - Environment & env, - ImageLoader & ldr) - { - // Sort the batch's labels to enable binary searching in it - std::vector whiteList = batch_labels; - std::sort (whiteList.begin(), whiteList.end()); - - int lvl = 0, // Pyramid level - lyr = 0; // Layer - - // Read the tiffs - size_t nth = ldr.get_num_tiles_hor(), - ntv = ldr.get_num_tiles_vert(), - fw = ldr.get_tile_width(), - th = ldr.get_tile_height(), - tw = ldr.get_tile_width(), - tileSize = ldr.get_tile_size(), - fullwidth = ldr.get_full_width(), - fullheight = ldr.get_full_height(); - - int cnt = 1; - for (unsigned int row = 0; row < nth; row++) - for (unsigned int col = 0; col < ntv; col++) - { - // Fetch the tile - bool ok = ldr.load_tile(row, col); - if (!ok) - { - std::stringstream ss; - ss << "Error fetching tile row=" << row << " col=" << col; - #ifdef WITH_PYTHON_H - throw ss.str(); - #endif - std::cerr << ss.str() << "\n"; - return false; - } - - // Get ahold of tile's pixel buffer - const std::vector& dataI = ldr.get_int_tile_buffer(); - const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); - bool wholeslide = spL == nullptr; // alternatively, theEnvironment.singleROI - - // Iterate pixels - for (unsigned long i = 0; i < tileSize; i++) - { - // mask label if not in the wholeslide mode - PixIntens label = 1; - if (!wholeslide) - label = (*spL)[i]; - - // Skip this ROI if the label isn't in the pending set of a multi-ROI mode - if (! env.singleROI && ! std::binary_search(whiteList.begin(), whiteList.end(), label)) - continue; - - auto inten = dataI[i]; - int y = row * th + i / tw, - x = col * tw + i % tw; - - // Skip tile buffer pixels beyond the image's bounds - if (x >= fullwidth || y >= fullheight) - continue; - - // Cache this pixel - LR& r = env.roiData [label]; - feed_pixel_2_cache_LR (x, y, dataI[i], r); - } - - VERBOSLVL2 (env.get_verbosity_level(), - // Show stayalive progress info - if (cnt++ % 4 == 0) - { - static int prevIntPc = 0; - float pc = int((row * nth + col) * 100 / float(nth * ntv) * 100) / 100.; - if (int(pc) != prevIntPc) - { - std::cout << "\t scan trivial " << int(pc) << " %\n"; - prevIntPc = int(pc); - } - } - ); - } - - return true; - } - - 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) - { - // Sort the batch's labels to enable binary searching in it - std::vector whiteList = batch_labels; - std::sort(whiteList.begin(), whiteList.end()); - - int lvl = 0, // pyramid level - lyr = 0; // layer - - // physical slide properties - size_t nth = ldr.get_num_tiles_hor(), - ntv = ldr.get_num_tiles_vert(), - fw = ldr.get_tile_width(), - th = ldr.get_tile_height(), - tw = ldr.get_tile_width(), - tileSize = ldr.get_tile_size(), - fullwidth = ldr.get_full_width(), - fullheight = ldr.get_full_height(); - - // virtual slide properties - size_t vh = (size_t) (double(fullheight) * sf_y), - vw = (size_t) (double(fullwidth) * sf_x), - vth = (size_t)(double(th) * sf_y), - vtw = (size_t)(double(tw) * sf_x); - - // current tile to skip tile reloads - size_t curt_x = 999, curt_y = 999; - - for (size_t vr = 0; vr < vh; vr++) - { - for (size_t vc = 0; vc < vw; vc++) - { - // tile position - size_t tidx_y = size_t(vr / vth), - tidx_x = size_t(vc / vtw); - - // load it - if (tidx_y != curt_y || tidx_x != curt_x) - { - bool ok = ldr.load_tile(tidx_y, tidx_x); - if (!ok) - { - std::string s = "Error fetching tile row=" + std::to_string(tidx_y) + " col=" + std::to_string(tidx_x); -#ifdef WITH_PYTHON_H - throw s; -#endif - std::cerr << s << "\n"; - return false; - } - - // cache tile position to avoid reloading - curt_y = tidx_y; - curt_x = tidx_x; - } - - // within-tile virtual pixel position - size_t vx = vc - tidx_x * vtw, - vy = vr - tidx_y * vth; - - // within-tile physical pixel position - size_t ph_x = size_t (double(vx) / sf_x), - ph_y = size_t (double(vy) / sf_y), - i = ph_y * tw + ph_x; - - // read buffered physical pixel - const std::vector& dataI = ldr.get_int_tile_buffer(); - const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); - bool wholeslide = spL == nullptr; // alternatively, theEnvironment.singleROI - - PixIntens label = 1; - if (!wholeslide) - label = (*spL)[i]; - - // not a ROI ? - if (!label) - continue; - - // skip this ROI if the label isn't in the to-do list 'whiteList' that's only possible in multi-ROI mode - if (wholeslide==false && !std::binary_search(whiteList.begin(), whiteList.end(), label)) - continue; - - auto inten = dataI[i]; - - // cache this pixel - // (ROI 'label' is known to the cache by means of gatherRoisMetrics() called previously.) - LR& r = env.roiData [label]; - feed_pixel_2_cache_LR (vc, vr, inten, r); - } - } - - return true; - } - - // - // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' - // - bool scan_trivial_wholeslide ( - LR & vroi, - const std::string& intens_fpath, - ImageLoader& ldr) - { - int lvl = 0, // Pyramid level - lyr = 0; // Layer - - // Read the tiffs - size_t nth = ldr.get_num_tiles_hor(), - ntv = ldr.get_num_tiles_vert(), - fw = ldr.get_tile_width(), - th = ldr.get_tile_height(), - tw = ldr.get_tile_width(), - tileSize = ldr.get_tile_size(), - fullwidth = ldr.get_full_width(), - fullheight = ldr.get_full_height(); - - int cnt = 1; - for (unsigned int row = 0; row < nth; row++) - for (unsigned int col = 0; col < ntv; col++) - { - // Fetch the tile - bool ok = ldr.load_tile(row, col); - if (!ok) - { - std::stringstream ss; - ss << "Error fetching tile row=" << row << " col=" << col; -#ifdef WITH_PYTHON_H - throw ss.str(); -#endif - std::cerr << ss.str() << "\n"; - return false; - } - - // Get ahold of tile's pixel buffer - const std::vector& dataI = ldr.get_int_tile_buffer(); - - // Iterate pixels - for (unsigned long i = 0; i < tileSize; i++) - { - auto inten = dataI[i]; - int y = row * th + i / tw, - x = col * tw + i % tw; - - // Skip tile buffer pixels beyond the image's bounds - if (x >= fullwidth || y >= fullheight) - continue; - - // Cache this pixel - feed_pixel_2_cache_LR (x, y, dataI[i], vroi); - } - } - - return true; - } - - // - // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' - // performing anisotropy correction - // - bool scan_trivial_wholeslide_anisotropic ( - LR& vroi, - const std::string& intens_fpath, - ImageLoader& ldr, - double aniso_x, - double aniso_y) - { - int lvl = 0, // Pyramid level - lyr = 0; // Layer - - // physical slide properties - size_t nth = ldr.get_num_tiles_hor(), - ntv = ldr.get_num_tiles_vert(), - fw = ldr.get_tile_width(), - th = ldr.get_tile_height(), - tw = ldr.get_tile_width(), - tileSize = ldr.get_tile_size(), - fullwidth = ldr.get_full_width(), - fullheight = ldr.get_full_height(); - - // virtual slide properties - size_t vh = (size_t)(double(fullheight) * aniso_y), - vw = (size_t)(double(fullwidth) * aniso_x), - vth = (size_t)(double(th) * aniso_y), - vtw = (size_t)(double(tw) * aniso_x); - - // current tile to skip tile reloads - size_t curt_x = 999, curt_y = 999; - - for (size_t vr = 0; vr < vh; vr++) - { - for (size_t vc = 0; vc < vw; vc++) - { - // tile position for virtual pixel (vc, vr) - size_t tidx_y = size_t(vr / vth), - tidx_x = size_t(vc / vtw); - - // load it - if (tidx_y != curt_y || tidx_x != curt_x) - { - bool ok = ldr.load_tile(tidx_y, tidx_x); - if (!ok) - { - std::string s = "Error fetching tile row=" + std::to_string(tidx_y) + " col=" + std::to_string(tidx_x); -#ifdef WITH_PYTHON_H - throw s; -#endif - std::cerr << s << "\n"; - return false; - } - - // cache tile position to avoid reloading - curt_y = tidx_y; - curt_x = tidx_x; - } - - // within-tile virtual pixel position - size_t vx = vc - tidx_x * vtw, - vy = vr - tidx_y * vth; - - // within-tile physical pixel position - size_t ph_x = size_t(double(vx) / aniso_x), - ph_y = size_t(double(vy) / aniso_y), - i = ph_y * tw + ph_x; - - // read buffered physical pixel - const std::vector& dataI = ldr.get_int_tile_buffer(); - const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); - bool wholeslide = spL == nullptr; // alternatively, theEnvironment.singleROI - - // Cache this pixel - feed_pixel_2_cache_LR (vc, vr, dataI[i], vroi); - } - } - - return true; - } - - void allocateTrivialRoisBuffers (const std::vector& Pending, Roidata& roiData, CpusideCache& cache) - { - // Calculate the total memory demand (in # of items) of all segments' image matrices - cache.imageMatrixBufferLen = 0; - for (auto lab : Pending) - { - LR& r = roiData[lab]; - - size_t w = r.aabb.get_width(), - h = r.aabb.get_height(), - imatrSize = w * h; - cache.imageMatrixBufferLen += imatrSize; - - cache.largest_roi_imatr_buf_len = cache.largest_roi_imatr_buf_len == 0 ? imatrSize : std::max (cache.largest_roi_imatr_buf_len, imatrSize); - } - - // Lagest ROI - cache.largest_roi_imatr_buf_len = 0; - for (auto lab : Pending) - { - LR& r = roiData[lab]; - cache.largest_roi_imatr_buf_len = cache.largest_roi_imatr_buf_len ? std::max(cache.largest_roi_imatr_buf_len, r.raw_pixels.size()) : r.raw_pixels.size(); - } - - // Allocate image matrices and remember each ROI's image matrix offset in 'ImageMatrixBuffer' - size_t baseIdx = 0; - for (auto lab : Pending) - { - LR& r = roiData[lab]; - - // matrix data - size_t imatrSize = r.aabb.get_width() * r.aabb.get_height(); - r.aux_image_matrix.allocate (r.aabb.get_width(), r.aabb.get_height()); - baseIdx += imatrSize; - - // Calculate the image matrix or cube - r.aux_image_matrix.calculate_from_pixelcloud (r.raw_pixels, r.aabb); - } - } - - void freeTrivialRoisBuffers (const std::vector& roi_labels, Roidata& roiData) - { - // Dispose memory of ROIs having their feature calculation finished - // in order to give memory ROIs of the next ROI batch. - // (Vector 'Pending' is the set of ROIs of the finished batch.) - for (auto lab : roi_labels) - { - LR& r = roiData[lab]; - std::vector().swap(r.raw_pixels); - std::vector().swap(r.aux_image_matrix._pix_plane); - std::vector().swap(r.convHull_CH); // convex hull is not a large object but there's no point to keep it beyond the batch, unlike contour - } - } - - void freeTrivialRoisBuffers_3D (const std::vector& roi_labels, Roidata& roiData) - { - // Dispose memory of ROIs having their feature calculation finished - // in order to give memory ROIs of the next ROI batch. - // (Vector 'Pending' is the set of ROIs of the finished batch.) - for (auto lab : roi_labels) - { - LR& r = roiData[lab]; - std::vector().swap(r.raw_pixels_3D); - - // - // Deallocate image matrices, cubes, and convex shells here (in the future). - // - } - - // Dispose the buffer of batches' ROIs' image matrices. (We allocate the - // image matrix buffer externally to minimize host-GPU transfers.) -// delete ImageMatrixBuffer; - } - - bool processTrivialRois (Environment & env, const std::vector& trivRoiLabels, const std::string& intens_fpath, const std::string& label_fpath, size_t memory_limit) - { - std::vector Pending; - size_t batchDemand = 0; - int roiBatchNo = 1; - - for (auto lab : trivRoiLabels) - { - LR& r = env.roiData[lab]; - - size_t itemFootprint = r.get_ram_footprint_estimate (trivRoiLabels.size()); - - // Check if we are good to accumulate this ROI in the current batch or should close the batch and reduce it - if (batchDemand + itemFootprint < memory_limit) - { - // There is room in the ROI batch. Insert another ROI in it - Pending.push_back(lab); - batchDemand += itemFootprint; - } - else - { - // The ROI batch is full. Let's process it - std::sort(Pending.begin(), Pending.end()); - VERBOSLVL2 (env.get_verbosity_level(), std::cout << ">>> Scanning batch #" << roiBatchNo << " of " << Pending.size() << " pending ROIs of total " << env.uniqueLabels.size() << " ROIs\n"); - VERBOSLVL2(env.get_verbosity_level(), - if (Pending.size() == 1) - std::cout << ">>> (single ROI label " << Pending[0] << ")\n"; - else - std::cout << ">>> (ROI labels " << Pending[0] << " ... " << Pending[Pending.size() - 1] << ")\n"; - ); - - if (env.anisoOptions.customized() == false) - scanTrivialRois (Pending, intens_fpath, label_fpath, env, env.theImLoader); - else - { - double ax = env.anisoOptions.get_aniso_x(), - ay = env.anisoOptions.get_aniso_y(); - scanTrivialRois_anisotropic (Pending, intens_fpath, label_fpath, env, env.theImLoader, ax, ay); - } - - // Allocate memory - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tallocating ROI buffers\n"); - allocateTrivialRoisBuffers (Pending, env.roiData, env.hostCache); - - // Reduce them - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing ROIs\n"); - // reduce_trivial_rois(Pending); - reduce_trivial_rois_manual (Pending, env); - - // Free memory - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tfreeing ROI buffers\n"); - freeTrivialRoisBuffers (Pending, env.roiData); // frees what's allocated by feed_pixel_2_cache() and allocateTrivialRoisBuffers() - - // Reset the RAM footprint accumulator - batchDemand = 0; - - // Clear the freshly processed ROIs from pending list - Pending.clear(); - - // Start a new pending set by adding the batch-overflowing ROI - Pending.push_back(lab); - - // Advance the batch counter - roiBatchNo++; - } - - // Allow keyboard interrupt - #ifdef WITH_PYTHON_H - if (PyErr_CheckSignals() != 0) - { - sureprint("\nAborting per user input\n"); - throw pybind11::error_already_set(); - } - #endif - } - - // Process what's remaining pending - if (Pending.size() > 0) - { - // Scan pixels of pending trivial ROIs - std::sort (Pending.begin(), Pending.end()); - VERBOSLVL2 (env.get_verbosity_level(), std::cout << ">>> Scanning batch #" << roiBatchNo << " of " << Pending.size() << " pending ROIs of " << env.uniqueLabels.size() << " all ROIs\n"); - VERBOSLVL2 (env.get_verbosity_level(), - if (Pending.size() == 1) - std::cout << ">>> (single ROI " << Pending[0] << ")\n"; - else - std::cout << ">>> (ROIs " << Pending[0] << " ... " << Pending[Pending.size() - 1] << ")\n"; - ); - - if (env.anisoOptions.customized() == false) - { - scanTrivialRois (Pending, intens_fpath, label_fpath, env, env.theImLoader); - } - else - { - double ax = env.anisoOptions.get_aniso_x(), - ay = env.anisoOptions.get_aniso_y(); - scanTrivialRois_anisotropic (Pending, intens_fpath, label_fpath, env, env.theImLoader, ax, ay); - } - - // Allocate memory - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tallocating ROI buffers\n"); - allocateTrivialRoisBuffers (Pending, env.roiData, env.hostCache); - - // Dump ROIs for use in unit testing -#ifdef DUMP_ALL_ROI - dump_all_roi(); -#endif - - // Reduce them - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing ROIs\n"); - //reduce_trivial_rois(Pending): - reduce_trivial_rois_manual (Pending, env); - - // Free memory - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tfreeing ROI buffers\n"); - freeTrivialRoisBuffers (Pending, env.roiData); - - #ifdef WITH_PYTHON_H - // Allow keyboard interrupt - if (PyErr_CheckSignals() != 0) - { - sureprint("\nAborting per user input\n"); - throw pybind11::error_already_set(); - } - #endif - } - - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing neighbor features and their depends for all ROIs\n"); - reduce_neighbors_and_dependencies_manual (env); - - return true; - } - -#ifdef WITH_PYTHON_H - - bool scanTrivialRoisInMemory ( - const std::vector& batch_labels, - const py::array_t& intens_images, - const py::array_t& label_images, - int pair_idx, - Environment & env) - { - // Sort the batch's labels to enable binary searching in it - std::vector whiteList = batch_labels; - std::sort(whiteList.begin(), whiteList.end()); - - auto rI = intens_images.unchecked<3>(); - auto rL = label_images.unchecked<3>(); - size_t width = rI.shape(2); - size_t height = rI.shape(1); - - int cnt = 1; - for (size_t col = 0; col < width; col++) - for (size_t row = 0; row < height; row++) - { - // Skip non-mask pixels - auto label = rL (pair_idx, row, col); - if (!label) - continue; - - // Skip this ROI if the label isn't in the pending set of a multi-ROI mode - if (! env.singleROI && !std::binary_search(whiteList.begin(), whiteList.end(), label)) - continue; - - auto inten = rI (pair_idx, row, col); - - // Collapse all the labels to one if single-ROI mde is requested - if (env.singleROI) - label = 1; - - // Cache this pixel - LR& r = env.roiData [label]; - feed_pixel_2_cache_LR (col, row, inten, r); - } - - return true; - } - - bool processTrivialRoisInMemory (Environment& env, const std::vector& trivRoiLabels, const py::array_t& intens, const py::array_t& label, int pair_idx, size_t memory_limit) - { - VERBOSLVL4(env.get_verbosity_level(), std::cout << "processTrivialRoisInMemory (pair_idx=" << pair_idx << ") \n"); - - std::vector Pending; - size_t batchDemand = 0; - int roiBatchNo = 1; - - for (auto lab : trivRoiLabels) - { - LR& r = env.roiData[lab]; - - size_t itemFootprint = r.get_ram_footprint_estimate(env.uniqueLabels.size()); - - // Check if we are good to accumulate this ROI in the current batch or should close the batch and reduce it - if (batchDemand + itemFootprint < memory_limit) - { - Pending.push_back(lab); - batchDemand += itemFootprint; - } - else - { - // Scan pixels of pending trivial ROIs - std::sort(Pending.begin(), Pending.end()); - scanTrivialRoisInMemory (Pending, intens, label, pair_idx, env); - - // Allocate memory - allocateTrivialRoisBuffers (Pending, env.roiData, env.hostCache); - - // reduce_trivial_rois(Pending); - reduce_trivial_rois_manual (Pending, env); - - // Free memory - freeTrivialRoisBuffers (Pending, env.roiData); // frees what's allocated by feed_pixel_2_cache() and allocateTrivialRoisBuffers() - - // Reset the RAM footprint accumulator - batchDemand = 0; - - // Clear the freshly processed ROIs from pending list - Pending.clear(); - - // Start a new pending set by adding the stopper ROI - Pending.push_back(lab); - - // Advance the batch counter - roiBatchNo++; - } - - // Allow keyboard interrupt - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); - - } - - // Process what's remaining pending - if (Pending.size() > 0) - { - // Scan pixels of pending trivial ROIs - std::sort(Pending.begin(), Pending.end()); - scanTrivialRoisInMemory(Pending, intens, label, pair_idx, env); - - // Allocate memory - allocateTrivialRoisBuffers (Pending, env.roiData, env.hostCache); - - // Dump ROIs for use in unit testing -#ifdef DUMP_ALL_ROI - dump_all_roi(); -#endif - - // Reduce them - reduce_trivial_rois_manual (Pending, env); - - // Free memory - freeTrivialRoisBuffers (Pending, env.roiData); - - // Allow keyboard interrupt - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); - } - - reduce_neighbors_and_dependencies_manual (env); - - return true; - - } - -#endif // python api - -} +#include +#include +#include +#include +#include +#include +#include + +#ifdef WITH_PYTHON_H + #include +#endif + +#include "environment.h" +#include "globals.h" +#include "helpers/timing.h" + +namespace Nyxus +{ + + #define disable_DUMP_ALL_ROI + #ifdef DUMP_ALL_ROI + void dump_all_roi() + { + std::string fpath = theEnvironment.output_dir + "/all_roi.txt"; + std::cout << "Dumping all the ROIs to " << fpath << " ...\n"; + + std::ofstream f(fpath); + + for (auto lab : uniqueLabels) + { + auto& r = roiData[lab]; + std::cout << "Dumping ROI " << lab << "\n"; + + r.aux_image_matrix.print(f); + + f << "ROI " << lab << ": \n" + << "xmin = " << r.aabb.get_xmin() << "; \n" + << "width=" << r.aabb.get_width() << "; \n" + << "ymin=" << r.aabb.get_ymin() << "; \n" + << "height=" << r.aabb.get_height() << "; \n" + << "area=" << r.aux_area << "; \n"; + + // C++ constant: + f << "// C:\n" + << "struct NyxusPixel {\n" + << "\tsize_t x, y; \n" + << "\tunsigned int intensity; \n" + << "}; \n" + << "NyxusPixel testData[] = {\n"; + for (auto i=0; i 0 && i % 4 == 0) + f << "\n"; + } + f << "}; \n"; + + // Matlab constant: + f << "// MATLAB:\n" + << "%==== begin \n"; + f << "pixelCloud = [ \n"; + for (auto i = 0; i < r.raw_pixels.size(); i++) + { + auto& px = r.raw_pixels[i]; + f << px.inten << "; % [" << i << "] \n"; + } + f << "]; \n"; + + f << "testData = zeros(" << r.aabb.get_height() << "," << r.aabb.get_width() << ");\n"; + for (auto i = 0; i < r.raw_pixels.size(); i++) + { + auto& px = r.raw_pixels[i]; + f << "testData(" << (px.y - r.aabb.get_ymin() + 1) << "," << (px.x - r.aabb.get_xmin() + 1) << ")=" << px.inten << "; "; // +1 due to 1-based nature of Matlab + if (i > 0 && i % 4 == 0) + f << "\n"; + } + f << "\n"; + f << "testVecZ = reshape(testData, 1, []); \n"; + f << "testVecNZ = nonzeros(testData); \n"; + f << "[mean(testVecNZ) mean(testVecZ) mean2(testData)] \n"; + f << "%==== end \n"; + } + + f.flush(); + } + #endif + + bool scanTrivialRois ( + const std::vector& batch_labels, + const std::string& intens_fpath, + const std::string& label_fpath, + Environment & env, + ImageLoader & ldr) + { + // Sort the batch's labels to enable binary searching in it + std::vector whiteList = batch_labels; + std::sort (whiteList.begin(), whiteList.end()); + + int lvl = 0, // Pyramid level + lyr = 0; // Layer + + // Read the tiffs + size_t nth = ldr.get_num_tiles_hor(), + ntv = ldr.get_num_tiles_vert(), + fw = ldr.get_tile_width(), + th = ldr.get_tile_height(), + tw = ldr.get_tile_width(), + tileSize = ldr.get_tile_size(), + fullwidth = ldr.get_full_width(), + fullheight = ldr.get_full_height(); + + int cnt = 1; + for (unsigned int row = 0; row < nth; row++) + for (unsigned int col = 0; col < ntv; col++) + { + // Fetch the tile + bool ok = ldr.load_tile(row, col); + if (!ok) + { + std::stringstream ss; + ss << "Error fetching tile row=" << row << " col=" << col; + #ifdef WITH_PYTHON_H + throw ss.str(); + #endif + std::cerr << ss.str() << "\n"; + return false; + } + + // Get ahold of tile's pixel buffer + const std::vector& dataI = ldr.get_int_tile_buffer(); + const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); + bool wholeslide = env.singleROI; + + // Iterate pixels + for (unsigned long i = 0; i < tileSize; i++) + { + // mask label if not in the wholeslide mode + PixIntens label = 1; + if (!wholeslide) + label = (*spL)[i]; + + // Skip this ROI if the label isn't in the pending set of a multi-ROI mode + if (! env.singleROI && ! std::binary_search(whiteList.begin(), whiteList.end(), label)) + continue; + + auto inten = dataI[i]; + int y = row * th + i / tw, + x = col * tw + i % tw; + + // Skip tile buffer pixels beyond the image's bounds + if (x >= fullwidth || y >= fullheight) + continue; + + // Cache this pixel + LR& r = env.roiData [label]; + feed_pixel_2_cache_LR (x, y, dataI[i], r); + } + + VERBOSLVL2 (env.get_verbosity_level(), + // Show stayalive progress info + if (cnt++ % 4 == 0) + { + static int prevIntPc = 0; + float pc = int((row * nth + col) * 100 / float(nth * ntv) * 100) / 100.; + if (int(pc) != prevIntPc) + { + std::cout << "\t scan trivial " << int(pc) << " %\n"; + prevIntPc = int(pc); + } + } + ); + } + + return true; + } + + 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) + { + // Sort the batch's labels to enable binary searching in it + std::vector whiteList = batch_labels; + std::sort(whiteList.begin(), whiteList.end()); + + int lvl = 0, // pyramid level + lyr = 0; // layer + + // physical slide properties + size_t nth = ldr.get_num_tiles_hor(), + ntv = ldr.get_num_tiles_vert(), + fw = ldr.get_tile_width(), + th = ldr.get_tile_height(), + tw = ldr.get_tile_width(), + tileSize = ldr.get_tile_size(), + fullwidth = ldr.get_full_width(), + fullheight = ldr.get_full_height(); + + // virtual slide properties + size_t vh = (size_t) (double(fullheight) * sf_y), + vw = (size_t) (double(fullwidth) * sf_x), + vth = (size_t)(double(th) * sf_y), + vtw = (size_t)(double(tw) * sf_x); + + // current tile to skip tile reloads + size_t curt_x = 999, curt_y = 999; + + for (size_t vr = 0; vr < vh; vr++) + { + for (size_t vc = 0; vc < vw; vc++) + { + // tile position + size_t tidx_y = size_t(vr / vth), + tidx_x = size_t(vc / vtw); + + // load it + if (tidx_y != curt_y || tidx_x != curt_x) + { + bool ok = ldr.load_tile(tidx_y, tidx_x); + if (!ok) + { + std::string s = "Error fetching tile row=" + std::to_string(tidx_y) + " col=" + std::to_string(tidx_x); +#ifdef WITH_PYTHON_H + throw s; +#endif + std::cerr << s << "\n"; + return false; + } + + // cache tile position to avoid reloading + curt_y = tidx_y; + curt_x = tidx_x; + } + + // within-tile virtual pixel position + size_t vx = vc - tidx_x * vtw, + vy = vr - tidx_y * vth; + + // within-tile physical pixel position + size_t ph_x = size_t (double(vx) / sf_x), + ph_y = size_t (double(vy) / sf_y), + i = ph_y * tw + ph_x; + + // read buffered physical pixel + const std::vector& dataI = ldr.get_int_tile_buffer(); + const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); + bool wholeslide = env.singleROI; + + PixIntens label = 1; + if (!wholeslide) + label = (*spL)[i]; + + // not a ROI ? + if (!label) + continue; + + // skip this ROI if the label isn't in the to-do list 'whiteList' that's only possible in multi-ROI mode + if (wholeslide==false && !std::binary_search(whiteList.begin(), whiteList.end(), label)) + continue; + + auto inten = dataI[i]; + + // cache this pixel + // (ROI 'label' is known to the cache by means of gatherRoisMetrics() called previously.) + LR& r = env.roiData [label]; + feed_pixel_2_cache_LR (vc, vr, inten, r); + } + } + + return true; + } + + // + // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' + // + bool scan_trivial_wholeslide ( + LR & vroi, + const std::string& intens_fpath, + ImageLoader& ldr) + { + int lvl = 0, // Pyramid level + lyr = 0; // Layer + + // Read the tiffs + size_t nth = ldr.get_num_tiles_hor(), + ntv = ldr.get_num_tiles_vert(), + fw = ldr.get_tile_width(), + th = ldr.get_tile_height(), + tw = ldr.get_tile_width(), + tileSize = ldr.get_tile_size(), + fullwidth = ldr.get_full_width(), + fullheight = ldr.get_full_height(); + + int cnt = 1; + for (unsigned int row = 0; row < nth; row++) + for (unsigned int col = 0; col < ntv; col++) + { + // Fetch the tile + bool ok = ldr.load_tile(row, col); + if (!ok) + { + std::stringstream ss; + ss << "Error fetching tile row=" << row << " col=" << col; +#ifdef WITH_PYTHON_H + throw ss.str(); +#endif + std::cerr << ss.str() << "\n"; + return false; + } + + // Get ahold of tile's pixel buffer + const std::vector& dataI = ldr.get_int_tile_buffer(); + + // Iterate pixels + for (unsigned long i = 0; i < tileSize; i++) + { + auto inten = dataI[i]; + int y = row * th + i / tw, + x = col * tw + i % tw; + + // Skip tile buffer pixels beyond the image's bounds + if (x >= fullwidth || y >= fullheight) + continue; + + // Cache this pixel + feed_pixel_2_cache_LR (x, y, dataI[i], vroi); + } + } + + return true; + } + + // + // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' + // performing anisotropy correction + // + bool scan_trivial_wholeslide_anisotropic ( + LR& vroi, + const std::string& intens_fpath, + ImageLoader& ldr, + double aniso_x, + double aniso_y) + { + int lvl = 0, // Pyramid level + lyr = 0; // Layer + + // physical slide properties + size_t nth = ldr.get_num_tiles_hor(), + ntv = ldr.get_num_tiles_vert(), + fw = ldr.get_tile_width(), + th = ldr.get_tile_height(), + tw = ldr.get_tile_width(), + tileSize = ldr.get_tile_size(), + fullwidth = ldr.get_full_width(), + fullheight = ldr.get_full_height(); + + // virtual slide properties + size_t vh = (size_t)(double(fullheight) * aniso_y), + vw = (size_t)(double(fullwidth) * aniso_x), + vth = (size_t)(double(th) * aniso_y), + vtw = (size_t)(double(tw) * aniso_x); + + // current tile to skip tile reloads + size_t curt_x = 999, curt_y = 999; + + for (size_t vr = 0; vr < vh; vr++) + { + for (size_t vc = 0; vc < vw; vc++) + { + // tile position for virtual pixel (vc, vr) + size_t tidx_y = size_t(vr / vth), + tidx_x = size_t(vc / vtw); + + // load it + if (tidx_y != curt_y || tidx_x != curt_x) + { + bool ok = ldr.load_tile(tidx_y, tidx_x); + if (!ok) + { + std::string s = "Error fetching tile row=" + std::to_string(tidx_y) + " col=" + std::to_string(tidx_x); +#ifdef WITH_PYTHON_H + throw s; +#endif + std::cerr << s << "\n"; + return false; + } + + // cache tile position to avoid reloading + curt_y = tidx_y; + curt_x = tidx_x; + } + + // within-tile virtual pixel position + size_t vx = vc - tidx_x * vtw, + vy = vr - tidx_y * vth; + + // within-tile physical pixel position + size_t ph_x = size_t(double(vx) / aniso_x), + ph_y = size_t(double(vy) / aniso_y), + i = ph_y * tw + ph_x; + + // read buffered physical pixel + const std::vector& dataI = ldr.get_int_tile_buffer(); + const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); + bool wholeslide = env.singleROI; + + // Cache this pixel + feed_pixel_2_cache_LR (vc, vr, dataI[i], vroi); + } + } + + return true; + } + + void allocateTrivialRoisBuffers (const std::vector& Pending, Roidata& roiData, CpusideCache& cache) + { + // Calculate the total memory demand (in # of items) of all segments' image matrices + cache.imageMatrixBufferLen = 0; + for (auto lab : Pending) + { + LR& r = roiData[lab]; + + size_t w = r.aabb.get_width(), + h = r.aabb.get_height(), + imatrSize = w * h; + cache.imageMatrixBufferLen += imatrSize; + + cache.largest_roi_imatr_buf_len = cache.largest_roi_imatr_buf_len == 0 ? imatrSize : std::max (cache.largest_roi_imatr_buf_len, imatrSize); + } + + // Lagest ROI + cache.largest_roi_imatr_buf_len = 0; + for (auto lab : Pending) + { + LR& r = roiData[lab]; + cache.largest_roi_imatr_buf_len = cache.largest_roi_imatr_buf_len ? std::max(cache.largest_roi_imatr_buf_len, r.raw_pixels.size()) : r.raw_pixels.size(); + } + + // Allocate image matrices and remember each ROI's image matrix offset in 'ImageMatrixBuffer' + size_t baseIdx = 0; + for (auto lab : Pending) + { + LR& r = roiData[lab]; + + // matrix data + size_t imatrSize = r.aabb.get_width() * r.aabb.get_height(); + r.aux_image_matrix.allocate (r.aabb.get_width(), r.aabb.get_height()); + baseIdx += imatrSize; + + // Calculate the image matrix or cube + r.aux_image_matrix.calculate_from_pixelcloud (r.raw_pixels, r.aabb); + } + } + + void freeTrivialRoisBuffers (const std::vector& roi_labels, Roidata& roiData) + { + // Dispose memory of ROIs having their feature calculation finished + // in order to give memory ROIs of the next ROI batch. + // (Vector 'Pending' is the set of ROIs of the finished batch.) + for (auto lab : roi_labels) + { + LR& r = roiData[lab]; + std::vector().swap(r.raw_pixels); + std::vector().swap(r.aux_image_matrix._pix_plane); + std::vector().swap(r.convHull_CH); // convex hull is not a large object but there's no point to keep it beyond the batch, unlike contour + } + } + + void freeTrivialRoisBuffers_3D (const std::vector& roi_labels, Roidata& roiData) + { + // Dispose memory of ROIs having their feature calculation finished + // in order to give memory ROIs of the next ROI batch. + // (Vector 'Pending' is the set of ROIs of the finished batch.) + for (auto lab : roi_labels) + { + LR& r = roiData[lab]; + std::vector().swap(r.raw_pixels_3D); + + // + // Deallocate image matrices, cubes, and convex shells here (in the future). + // + } + + // Dispose the buffer of batches' ROIs' image matrices. (We allocate the + // image matrix buffer externally to minimize host-GPU transfers.) +// delete ImageMatrixBuffer; + } + + bool processTrivialRois (Environment & env, const std::vector& trivRoiLabels, const std::string& intens_fpath, const std::string& label_fpath, size_t memory_limit) + { + std::vector Pending; + size_t batchDemand = 0; + int roiBatchNo = 1; + + for (auto lab : trivRoiLabels) + { + LR& r = env.roiData[lab]; + + size_t itemFootprint = r.get_ram_footprint_estimate (trivRoiLabels.size()); + + // Check if we are good to accumulate this ROI in the current batch or should close the batch and reduce it + if (batchDemand + itemFootprint < memory_limit) + { + // There is room in the ROI batch. Insert another ROI in it + Pending.push_back(lab); + batchDemand += itemFootprint; + } + else + { + // The ROI batch is full. Let's process it + std::sort(Pending.begin(), Pending.end()); + VERBOSLVL2 (env.get_verbosity_level(), std::cout << ">>> Scanning batch #" << roiBatchNo << " of " << Pending.size() << " pending ROIs of total " << env.uniqueLabels.size() << " ROIs\n"); + VERBOSLVL2(env.get_verbosity_level(), + if (Pending.size() == 1) + std::cout << ">>> (single ROI label " << Pending[0] << ")\n"; + else + std::cout << ">>> (ROI labels " << Pending[0] << " ... " << Pending[Pending.size() - 1] << ")\n"; + ); + + if (env.anisoOptions.customized() == false) + scanTrivialRois (Pending, intens_fpath, label_fpath, env, env.theImLoader); + else + { + double ax = env.anisoOptions.get_aniso_x(), + ay = env.anisoOptions.get_aniso_y(); + scanTrivialRois_anisotropic (Pending, intens_fpath, label_fpath, env, env.theImLoader, ax, ay); + } + + // Allocate memory + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tallocating ROI buffers\n"); + allocateTrivialRoisBuffers (Pending, env.roiData, env.hostCache); + + // Reduce them + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing ROIs\n"); + // reduce_trivial_rois(Pending); + reduce_trivial_rois_manual (Pending, env); + + // Free memory + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tfreeing ROI buffers\n"); + freeTrivialRoisBuffers (Pending, env.roiData); // frees what's allocated by feed_pixel_2_cache() and allocateTrivialRoisBuffers() + + // Reset the RAM footprint accumulator + batchDemand = 0; + + // Clear the freshly processed ROIs from pending list + Pending.clear(); + + // Start a new pending set by adding the batch-overflowing ROI + Pending.push_back(lab); + + // Advance the batch counter + roiBatchNo++; + } + + // Allow keyboard interrupt + #ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + { + sureprint("\nAborting per user input\n"); + throw pybind11::error_already_set(); + } + #endif + } + + // Process what's remaining pending + if (Pending.size() > 0) + { + // Scan pixels of pending trivial ROIs + std::sort (Pending.begin(), Pending.end()); + VERBOSLVL2 (env.get_verbosity_level(), std::cout << ">>> Scanning batch #" << roiBatchNo << " of " << Pending.size() << " pending ROIs of " << env.uniqueLabels.size() << " all ROIs\n"); + VERBOSLVL2 (env.get_verbosity_level(), + if (Pending.size() == 1) + std::cout << ">>> (single ROI " << Pending[0] << ")\n"; + else + std::cout << ">>> (ROIs " << Pending[0] << " ... " << Pending[Pending.size() - 1] << ")\n"; + ); + + if (env.anisoOptions.customized() == false) + { + scanTrivialRois (Pending, intens_fpath, label_fpath, env, env.theImLoader); + } + else + { + double ax = env.anisoOptions.get_aniso_x(), + ay = env.anisoOptions.get_aniso_y(); + scanTrivialRois_anisotropic (Pending, intens_fpath, label_fpath, env, env.theImLoader, ax, ay); + } + + // Allocate memory + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tallocating ROI buffers\n"); + allocateTrivialRoisBuffers (Pending, env.roiData, env.hostCache); + + // Dump ROIs for use in unit testing +#ifdef DUMP_ALL_ROI + dump_all_roi(); +#endif + + // Reduce them + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing ROIs\n"); + //reduce_trivial_rois(Pending): + reduce_trivial_rois_manual (Pending, env); + + // Free memory + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tfreeing ROI buffers\n"); + freeTrivialRoisBuffers (Pending, env.roiData); + + #ifdef WITH_PYTHON_H + // Allow keyboard interrupt + if (PyErr_CheckSignals() != 0) + { + sureprint("\nAborting per user input\n"); + throw pybind11::error_already_set(); + } + #endif + } + + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing neighbor features and their depends for all ROIs\n"); + reduce_neighbors_and_dependencies_manual (env); + + return true; + } + +#ifdef WITH_PYTHON_H + + bool scanTrivialRoisInMemory ( + const std::vector& batch_labels, + const py::array_t& intens_images, + const py::array_t& label_images, + int pair_idx, + Environment & env) + { + // Sort the batch's labels to enable binary searching in it + std::vector whiteList = batch_labels; + std::sort(whiteList.begin(), whiteList.end()); + + auto rI = intens_images.unchecked<3>(); + auto rL = label_images.unchecked<3>(); + size_t width = rI.shape(2); + size_t height = rI.shape(1); + + int cnt = 1; + for (size_t col = 0; col < width; col++) + for (size_t row = 0; row < height; row++) + { + // Skip non-mask pixels + auto label = rL (pair_idx, row, col); + if (!label) + continue; + + // Skip this ROI if the label isn't in the pending set of a multi-ROI mode + if (! env.singleROI && !std::binary_search(whiteList.begin(), whiteList.end(), label)) + continue; + + auto inten = rI (pair_idx, row, col); + + // Collapse all the labels to one if single-ROI mde is requested + if (env.singleROI) + label = 1; + + // Cache this pixel + LR& r = env.roiData [label]; + feed_pixel_2_cache_LR (col, row, inten, r); + } + + return true; + } + + bool processTrivialRoisInMemory (Environment& env, const std::vector& trivRoiLabels, const py::array_t& intens, const py::array_t& label, int pair_idx, size_t memory_limit) + { + VERBOSLVL4(env.get_verbosity_level(), std::cout << "processTrivialRoisInMemory (pair_idx=" << pair_idx << ") \n"); + + std::vector Pending; + size_t batchDemand = 0; + int roiBatchNo = 1; + + for (auto lab : trivRoiLabels) + { + LR& r = env.roiData[lab]; + + size_t itemFootprint = r.get_ram_footprint_estimate(env.uniqueLabels.size()); + + // Check if we are good to accumulate this ROI in the current batch or should close the batch and reduce it + if (batchDemand + itemFootprint < memory_limit) + { + Pending.push_back(lab); + batchDemand += itemFootprint; + } + else + { + // Scan pixels of pending trivial ROIs + std::sort(Pending.begin(), Pending.end()); + scanTrivialRoisInMemory (Pending, intens, label, pair_idx, env); + + // Allocate memory + allocateTrivialRoisBuffers (Pending, env.roiData, env.hostCache); + + // reduce_trivial_rois(Pending); + reduce_trivial_rois_manual (Pending, env); + + // Free memory + freeTrivialRoisBuffers (Pending, env.roiData); // frees what's allocated by feed_pixel_2_cache() and allocateTrivialRoisBuffers() + + // Reset the RAM footprint accumulator + batchDemand = 0; + + // Clear the freshly processed ROIs from pending list + Pending.clear(); + + // Start a new pending set by adding the stopper ROI + Pending.push_back(lab); + + // Advance the batch counter + roiBatchNo++; + } + + // Allow keyboard interrupt + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); + + } + + // Process what's remaining pending + if (Pending.size() > 0) + { + // Scan pixels of pending trivial ROIs + std::sort(Pending.begin(), Pending.end()); + scanTrivialRoisInMemory(Pending, intens, label, pair_idx, env); + + // Allocate memory + allocateTrivialRoisBuffers (Pending, env.roiData, env.hostCache); + + // Dump ROIs for use in unit testing +#ifdef DUMP_ALL_ROI + dump_all_roi(); +#endif + + // Reduce them + reduce_trivial_rois_manual (Pending, env); + + // Free memory + freeTrivialRoisBuffers (Pending, env.roiData); + + // Allow keyboard interrupt + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); + } + + reduce_neighbors_and_dependencies_manual (env); + + return true; + + } + +#endif // python api + +} diff --git a/src/nyx/phase2_3d.cpp b/src/nyx/phase2_3d.cpp index 661b902b..d4122070 100644 --- a/src/nyx/phase2_3d.cpp +++ b/src/nyx/phase2_3d.cpp @@ -1,821 +1,821 @@ -#include -#include -#include -#include -#include -#include -#include -#include - -#ifdef WITH_PYTHON_H - #include -#endif - -#include "environment.h" -#include "globals.h" -#include "helpers/fsystem.h" -#include "helpers/timing.h" - -namespace Nyxus -{ - // - // Loads ROI voxels into voxel clouds - // - bool scanTrivialRois_3D (Environment & env, const std::vector& batch_labels, const std::string& intens_fpath, const std::string& label_fpath, size_t t_index) - { - // Sort the batch's labels to enable binary searching in it - std::vector whiteList = batch_labels; - std::sort(whiteList.begin(), whiteList.end()); - - // Scan this Z intensity-mask pair - SlideProps p (intens_fpath, label_fpath); - if (! env.theImLoader.open(p, env.fpimageOptions)) - { - std::cerr << "Error opening a file pair with ImageLoader. Terminating\n"; - return false; - } - - // thanks to ImageLoader::open() we are guaranteed that the mask and intensity volumes' - // width, height, and depth match. Mask and intensity may only differ in the number of - // time frames: 1:1, 1:N, and N:1 cases are permitted. - size_t - /* - nth = env.theImLoader.get_num_tiles_hor(), - ntv = env.theImLoader.get_num_tiles_vert(), - fw = env.theImLoader.get_tile_width(), - th = env.theImLoader.get_tile_height(), - tw = env.theImLoader.get_tile_width(), - tileSize = env.theImLoader.get_tile_size(), - */ - w = env.theImLoader.get_full_width(), - h = env.theImLoader.get_full_height(), - d = env.theImLoader.get_full_depth(), - sliceSize = w * h, - timeFrameSize = sliceSize * d, - timeI = env.theImLoader.get_inten_time(), - timeM = env.theImLoader.get_mask_time(), - nVoxI = timeFrameSize * timeI, - nVoxM = timeFrameSize * timeM; - - // is this intensity-mask pair's shape supported? - if (nVoxI < nVoxM) - { - std::string erm = "Error: unsupported shape - intensity file: " + std::to_string(nVoxI) + ", mask file: " + std::to_string(nVoxM); -#ifdef WITH_PYTHON_H - throw erm; -#endif - std::cerr << erm << "\n"; - return false; - } - - int cnt = 1; - - // fetch 3D data - bool ok = env.theImLoader.load_tile (0/*row*/, 0/*col*/); - if (!ok) - { - std::string erm = "Error fetching segmented data from " + intens_fpath + "(I) " + label_fpath + "(M)"; -#ifdef WITH_PYTHON_H - throw erm; -#endif - std::cerr << erm << "\n"; - return false; - } - - // Get ahold of tile's pixel buffer - auto dataI = env.theImLoader.get_int_tile_buffer(), - dataL = env.theImLoader.get_seg_tile_buffer(); - - size_t baseI, baseM; - if (timeI == timeM) - { - baseM = baseI = t_index * timeFrameSize; - } - else // nVoxI > nVoxM - { - baseM = 0; - baseI = t_index * timeFrameSize; - } - - // Iterate voxels - for (size_t i=0; i= w || y >= h || z >= d) - continue; - - // Collapse all the labels to one if single-ROI mde is requested - if (env.singleROI) - label = 1; - - // Cache this pixel - LR& r = env.roiData[label]; - feed_pixel_2_cache_3D_LR (x, y, z, dataI[j], r); - } - -#ifdef WITH_PYTHON_H - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); -#endif - - // Close the image pair - env.theImLoader.close(); - - // Dump ROI pixel clouds to the output directory - VERBOSLVL5 (env.get_verbosity_level(), dump_roi_pixels(env.dim(), Nyxus::get_temp_dir_path(), batch_labels, label_fpath, env.uniqueLabels, env.roiData)); - - return true; - } - - bool scanTrivialRois_3D_anisotropic__BEFORE4D( - Environment & env, - const std::vector& batch_labels, - const std::string& intens_fpath, - const std::string& label_fpath, - double aniso_x, - double aniso_y, - double aniso_z) - { - // Sort the batch's labels to enable binary searching in it - std::vector whiteList = batch_labels; - std::sort(whiteList.begin(), whiteList.end()); - - int lvl = 0, // pyramid level - lyr = 0; // layer - - // Scan this Z intensity-mask pair - SlideProps p (intens_fpath, label_fpath); - if (! env.theImLoader.open(p, env.fpimageOptions)) - { - std::cerr << "Error opening a file pair with ImageLoader. Terminating\n"; - return false; - } - - size_t nth = env.theImLoader.get_num_tiles_hor(), - ntv = env.theImLoader.get_num_tiles_vert(), - fw = env.theImLoader.get_tile_width(), - th = env.theImLoader.get_tile_height(), - tw = env.theImLoader.get_tile_width(), - tileSize = env.theImLoader.get_tile_size(), - fullW = env.theImLoader.get_full_width(), - fullH = env.theImLoader.get_full_height(), - fullD = env.theImLoader.get_full_depth(); - - size_t vD = (size_t)(double(fullD) * aniso_z); // virtual depth - - for (size_t vz = 0; vz < vD; vz++) - { - size_t z = size_t(double(vz) / aniso_z); // physical z - - // virtual slide properties - size_t vh = (size_t)(double(fullH) * aniso_y), - vw = (size_t)(double(fullW) * aniso_x), - vth = (size_t)(double(th) * aniso_y), - vtw = (size_t)(double(tw) * aniso_x); - - // current tile to skip tile reloads - size_t curt_x = 999, curt_y = 999; - - for (size_t vr = 0; vr < vh; vr++) - { - for (size_t vc = 0; vc < vw; vc++) - { - // tile position - size_t tidx_y = size_t(vr / vth), - tidx_x = size_t(vc / vtw); - - // load it - if (tidx_y != curt_y || tidx_x != curt_x) - { - bool ok = env.theImLoader.load_tile(tidx_y, tidx_x); - if (!ok) - { - std::string s = "Error fetching tile row=" + std::to_string(tidx_y) + " col=" + std::to_string(tidx_x); -#ifdef WITH_PYTHON_H - throw s; -#endif - std::cerr << s << "\n"; - return false; - } - - // cache tile position to avoid reloading - curt_y = tidx_y; - curt_x = tidx_x; - } - - // within-tile virtual pixel position - size_t vx = vc - tidx_x * vtw, - vy = vr - tidx_y * vth; - - // within-tile physical pixel position - size_t ph_x = size_t(double(vx) / aniso_x), - ph_y = size_t(double(vy) / aniso_y), - i = ph_y * tw + ph_x; - - // read buffered physical pixel - auto dataI = env.theImLoader.get_int_tile_buffer(), - dataL = env.theImLoader.get_seg_tile_buffer(); - - // skip non-mask pixels - auto label = dataL[i]; - if (!label) - continue; - - // skip this ROI if the label isn't in the pending set of a multi-ROI mode - if (!env.singleROI && !std::binary_search(whiteList.begin(), whiteList.end(), label)) - continue; - - // skip tile buffer pixels beyond the image's bounds - if (vc >= fullW || vr >= fullH) - continue; - - // collapse all the labels to one if single-ROI mde is requested - if (env.singleROI) - label = 1; - - // cache this voxel - auto inten = dataI[i]; - LR& r = env.roiData[label]; - feed_pixel_2_cache_3D_LR (vc, vr, vz, dataI[i], r); - } - } - - // Close the image pair - env.theImLoader.close(); - } - - // Dump ROI pixel clouds to the output directory - VERBOSLVL5 (env.get_verbosity_level(), dump_roi_pixels(env.dim(), Nyxus::get_temp_dir_path(), batch_labels, label_fpath, env.uniqueLabels, env.roiData)); - - return true; - } - - // - // Loads ROI voxels into voxel clouds - // - 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) - { - // Sort the batch's labels to enable binary searching in it - std::vector whiteList = batch_labels; - std::sort (whiteList.begin(), whiteList.end()); - - // temp slideprops instance to pass some into to ImageLoader - SlideProps p (intens_fpath, label_fpath); - if (!env.theImLoader.open(p, env.fpimageOptions)) - { - std::cerr << "Error opening a file pair with ImageLoader. Terminating\n"; - return false; - } - - // thanks to ImageLoader::open() we are guaranteed that the mask's and intensity's - // W, H, and D match. Mask and intensity may only differ in the number of - // time frames: 1:1, 1:N, and N:1 cases are permitted - size_t - w = env.theImLoader.get_full_width(), - h = env.theImLoader.get_full_height(), - d = env.theImLoader.get_full_depth(), - slice = w * h, - timeFrameSize = slice * d, - timeI = env.theImLoader.get_inten_time(), - timeM = env.theImLoader.get_mask_time(), - nVoxI = timeFrameSize * timeI, - nVoxM = timeFrameSize * timeM; - - // is this intensity-mask pair's shape supported? - if (nVoxI < nVoxM) - { - std::string erm = "Error: unsupported shape - intensity file: " + std::to_string(nVoxI) + ", mask file: " + std::to_string(nVoxM); - #ifdef WITH_PYTHON_H - throw erm; - #endif - std::cerr << erm << "\n"; - return false; - } - - int cnt = 1; - - // fetch 3D data - if (!env.theImLoader.load_tile (0/*row*/, 0/*col*/)) - { - std::string erm = "Error fetching data from file pair " + intens_fpath + "(I) " + label_fpath + "(M)"; - #ifdef WITH_PYTHON_H - throw erm; - #endif - std::cerr << erm << "\n"; - return false; - } - - // get ahold of voxel buffers - auto dataI = env.theImLoader.get_int_tile_buffer(), - dataL = env.theImLoader.get_seg_tile_buffer(); - - // align time frame's mask and intensity volumes - size_t baseI, baseM; - if (timeI == timeM) - { - // trivial N mask : N intensity - baseM = - baseI = t_index * timeFrameSize; - } - else - { - // nontrivial 1 mask : N intensity - baseM = 0; - baseI = t_index * timeFrameSize; - } - - // virtual dimensions - size_t virt_h = h * aniso_y, - virt_w = w * aniso_x, - virt_d = d * aniso_z; - size_t vSliceLen = virt_h * virt_w, - virt_v = vSliceLen * virt_d; - - // iterate virtual voxels and fill them with corresponding physical intensities - for (size_t vIdx = 0; vIdx < virt_v; vIdx++) - { - // virtual Cartesian position - size_t vZ = vIdx / vSliceLen, - vLastSliceLen = vIdx % vSliceLen, - vY = vLastSliceLen / virt_w, - vX = vLastSliceLen % virt_w; - - // physical Cartesian position - size_t pZ = vZ / aniso_z + 0.5, - pY = vY / aniso_y + 0.5, - pX = vX / aniso_x + 0.5; - - // skip a position outside the bounds - // (since we are casting coorinates from virtual to physical, - // we may get positions outside the physical bounds) - if (pX >= w || pY >= h || pZ >= d) - continue; - - // physical offset - size_t i = pZ * slice + pY * w + pX; - - // - // interpret the mask intensity - // - - // skip non-mask pixels - auto lbl = dataL[baseM + i]; - if (!lbl) - continue; - - // skip this ROI if the label isn't in the pending set of a multi-ROI mode - if (!env.singleROI && !std::binary_search(whiteList.begin(), whiteList.end(), lbl)) - continue; - - // collapse all the labels to one if single-ROI mde is requested - if (env.singleROI) - lbl = 1; - -#if !defined(NDEBUG) - if (vZ >= virt_d) - std::cout << "vZ=" << vZ << " < virt_d =" << virt_d << "\n"; - assert(vZ < virt_d); - if (vY >= virt_h) - std::cout << "vY=" << vY << " < virt_h =" << virt_h << "\n"; - assert(vY < virt_h); - if (vX >= virt_w) - std::cout << "vX=" << vX << " < virt_w =" << virt_w << "\n"; - assert(vX < virt_w); -#endif - - // cache this voxel - auto inten = dataI[baseI + i]; - LR& r = env.roiData[lbl]; - feed_pixel_2_cache_3D_LR (vX, vY, vZ, inten, r); - } - - #ifdef WITH_PYTHON_H - // allow keyboard interrupt - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); - #endif - - // Close the image pair - env.theImLoader.close(); - return true; - } - - // - // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' - // - bool scan_trivial_wholevolume ( - LR& vroi, - const std::string& intens_fpath, - ImageLoader& ilo) - { - int lvl = 0, // Pyramid level - lyr = 0; // Layer - - // Read the tiffs - - size_t fullwidth = ilo.get_full_width(), - fullheight = ilo.get_full_height(), - fullD = ilo.get_full_depth(), - sliceSize = fullwidth * fullheight, - nVox = sliceSize * fullD; - - // in the 3D case tiling is a formality, so fetch the only tile in the file - if (!ilo.load_tile(0, 0)) - { -#ifdef WITH_PYTHON_H - throw "Error fetching tile"; -#endif - std::cerr << "Error fetching tile\n"; - return false; - } - - // Get ahold of tile's pixel buffer - const std::vector& dataI = ilo.get_int_tile_buffer(); - - // iterate abstract tiles (in a tiled slide /e.g. tiled tiff/ they correspond to physical tiles, in a nontiled slide /e.g. scanline tiff or strip tiff/ they correspond to ) - int cnt = 1; - - // iterate voxels - for (size_t i = 0; i < nVox; i++) - { - int z = i / sliceSize, - y = (i - z * sliceSize) / fullwidth, - x = (i - z * sliceSize) % fullwidth; - - // Skip tile buffer pixels beyond the image's bounds - if (x >= fullwidth || y >= fullheight || z >= fullD) - continue; - - // dynamic range within- and off-ROI - auto inten = dataI[i]; - - // Cache this pixel - feed_pixel_2_cache_3D_LR (x, y, z, inten, vroi); - - } //- all voxels - - return true; - } - - // - // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' - // - bool scan_trivial_wholevolume_anisotropic ( - LR& vroi, - const std::string& intens_fpath, - ImageLoader& ilo, - double aniso_x, - double aniso_y, - double aniso_z) - { - int lvl = 0, // Pyramid level - lyr = 0; // Layer - - // Read the tiffs - - size_t fullW = ilo.get_full_width(), - fullH = ilo.get_full_height(), - fullD = ilo.get_full_depth(), - sliceSize = fullW * fullH; - - size_t vh = (size_t) (double(fullH) * aniso_y), - vw = (size_t) (double(fullW) * aniso_x), - vd = (size_t) (double(fullD) * aniso_z); - - // in the 3D case tiling is a formality, so fetch the only tile in the file - if (! ilo.load_tile(0, 0)) - { -#ifdef WITH_PYTHON_H - throw "Error loading volume data"; -#endif - std::cerr << "Error loading volume data\n"; - return false; - } - - // Get ahold of tile's pixel buffer - const std::vector& dataI = ilo.get_int_tile_buffer(); - - // iterate virtual voxels - size_t vSliceSize = vh * vw, - nVox = vh * vw * vd; - for (size_t i = 0; i < nVox; i++) - { - // virtual voxel position - int z = i / vSliceSize, - y = (i - z * vSliceSize) / vw, - x = (i - z * vSliceSize) % vw; - - // physical voxel position - size_t ph_x = (size_t) (double(x) / aniso_x), - ph_y = (size_t) (double(y) / aniso_y), - ph_z = (size_t) (double(z) / aniso_z); - i = ph_z * sliceSize + ph_y * fullH + ph_x; - - // Cache this pixel - feed_pixel_2_cache_3D_LR (x, y, z, dataI[i], vroi); - - } - - return true; - } - - // - // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' - // performing anisotropy correction - // - bool scan_trivial_wholevolume_anisotropic__OLD( - LR& vroi, - const std::string& intens_fpath, - ImageLoader& ldr, - double aniso_x, - double aniso_y) - { - int lvl = 0, // Pyramid level - lyr = 0; // Layer - - // physical slide properties - size_t nth = ldr.get_num_tiles_hor(), - ntv = ldr.get_num_tiles_vert(), - fw = ldr.get_tile_width(), - th = ldr.get_tile_height(), - tw = ldr.get_tile_width(), - tileSize = ldr.get_tile_size(), - fullwidth = ldr.get_full_width(), - fullheight = ldr.get_full_height(); - - // virtual slide properties - size_t vh = (size_t)(double(fullheight) * aniso_y), - vw = (size_t)(double(fullwidth) * aniso_x), - vth = (size_t)(double(th) * aniso_y), - vtw = (size_t)(double(tw) * aniso_x); - - // current tile to skip tile reloads - size_t curt_x = 999, curt_y = 999; - - for (size_t vr = 0; vr < vh; vr++) - { - for (size_t vc = 0; vc < vw; vc++) - { - // tile position for virtual pixel (vc, vr) - size_t tidx_y = size_t(vr / vth), - tidx_x = size_t(vc / vtw); - - // load it - if (tidx_y != curt_y || tidx_x != curt_x) - { - bool ok = ldr.load_tile(tidx_y, tidx_x); - if (!ok) - { - std::string s = "Error fetching tile row=" + std::to_string(tidx_y) + " col=" + std::to_string(tidx_x); -#ifdef WITH_PYTHON_H - throw s; -#endif - std::cerr << s << "\n"; - return false; - } - - // cache tile position to avoid reloading - curt_y = tidx_y; - curt_x = tidx_x; - } - - // within-tile virtual pixel position - size_t vx = vc - tidx_x * vtw, - vy = vr - tidx_y * vth; - - // within-tile physical pixel position - size_t ph_x = size_t(double(vx) / aniso_x), - ph_y = size_t(double(vy) / aniso_y), - i = ph_y * tw + ph_x; - - // read buffered physical pixel - const std::vector& dataI = ldr.get_int_tile_buffer(); - const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); - bool wholeslide = spL == nullptr; // alternatively, theEnvironment.singleROI - - // Cache this pixel - feed_pixel_2_cache_LR(vc, vr, dataI[i], vroi); - } - } - - return true; - } - - - bool processTrivialRois_3D (Environment & env, size_t sidx, size_t t_index, const std::vector& trivRoiLabels, const std::string& intens_fpath, const std::string& label_fpath, size_t memory_limit) - { - std::vector Pending; - size_t batchDemand = 0; - int roiBatchNo = 1; - - for (auto lab : trivRoiLabels) - { - LR& r = env.roiData[lab]; - - size_t itemFootprint = r.get_ram_footprint_estimate (Pending.size()); - - // Check if we are good to accumulate this ROI in the current batch or should close the batch and reduce it - if (batchDemand + itemFootprint < memory_limit) - { - Pending.push_back(lab); - batchDemand += itemFootprint; - } - else - { - // Scan pixels of pending trivial ROIs - std::sort (Pending.begin(), Pending.end()); - - VERBOSLVL2 (env.get_verbosity_level(), std::cout << ">>> Scanning batch #" << roiBatchNo << " of " << Pending.size() << " pending ROIs of total " << env.uniqueLabels.size() << " ROIs\n"); - VERBOSLVL2 (env.get_verbosity_level(), - if (Pending.size() == 1) - std::cout << ">>> (single ROI label " << Pending[0] << ")\n"; - else - std::cout << ">>> (ROI labels " << Pending[0] << " ... " << Pending[Pending.size() - 1] << ")\n"; - ); - - if (env.anisoOptions.customized() == false) - { - scanTrivialRois_3D (env, Pending, 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, Pending, intens_fpath, label_fpath, t_index, ax, ay, az); - - // rescan and update ROI's AABB - for (auto lbl : Pending) - { - LR& r = env.roiData[lbl]; - r.aabb.update_from_voxelcloud (r.raw_pixels_3D); - } - } - - // Allocate memory - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tallocating ROI buffers\n";) - allocateTrivialRoisBuffers_3D (Pending, env.roiData, env.hostCache); - - // Reduce them - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing ROIs\n";) - // reduce_trivial_rois(Pending); - reduce_trivial_rois_manual (Pending, env); - - // Free memory - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tfreeing ROI buffers\n";) - freeTrivialRoisBuffers_3D (Pending, env.roiData); // frees what's allocated by feed_pixel_2_cache() and allocateTrivialRoisBuffers() - - // Reset the RAM footprint accumulator - batchDemand = 0; - - // Clear the freshly processed ROIs from pending list - Pending.clear(); - - // Start a new pending set by adding the stopper ROI - Pending.push_back(lab); - - // Advance the batch counter - roiBatchNo++; - } - - // Allow keyboard interrupt -#ifdef WITH_PYTHON_H - if (PyErr_CheckSignals() != 0) - { - sureprint("\nAborting per user input\n"); - throw pybind11::error_already_set(); - } -#endif - } - - // Process what's remaining pending - if (Pending.size() > 0) - { - // Read raw pixels of pending trivial ROIs - std::sort(Pending.begin(), Pending.end()); - - VERBOSLVL2(env.get_verbosity_level(), - std::cout << ">>> Scanning batch #" << roiBatchNo << " of " << Pending.size() << "(" << env.uniqueLabels.size() << ") ROIs\n"; - std::cout << ">>> (labels " << Pending[0] << " ... " << Pending[Pending.size() - 1] << ")\n"; - ); - - if (env.anisoOptions.customized() == false) - { - scanTrivialRois_3D (env, Pending, 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, Pending, intens_fpath, label_fpath, t_index, ax, ay, az); - - // rescan and update ROI's AABB - for (auto lbl : Pending) - { - LR& r = env.roiData[lbl]; - r.aabb.update_from_voxelcloud(r.raw_pixels_3D); - } - } - - for (auto lab : Pending) - { - LR& r = env.roiData[lab]; - for (Pixel3& vox : r.raw_pixels_3D) - { - assert (vox.x >= r.aabb.get_xmin()); - assert (vox.x <= r.aabb.get_xmax()); - assert (vox.y >= r.aabb.get_ymin()); - assert (vox.y <= r.aabb.get_ymax()); - assert (vox.z >= r.aabb.get_zmin()); - assert (vox.z <= r.aabb.get_zmax()); - } - } - - // Allocate memory - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tallocating ROI buffers\n"); - allocateTrivialRoisBuffers_3D (Pending, env.roiData, env.hostCache); - - // Dump ROIs for use in unit testing -#ifdef DUMP_ALL_ROI - dump_all_roi(); -#endif - - // Reduce them - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing ROIs\n"); - //reduce_trivial_rois(Pending); - reduce_trivial_rois_manual (Pending, env); - - // Free memory - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tfreeing ROI buffers\n"); - freeTrivialRoisBuffers_3D (Pending, env.roiData); - -#ifdef WITH_PYTHON_H - // Allow keyboard interrupt - if (PyErr_CheckSignals() != 0) - { - sureprint("\nAborting per user input\n"); - throw pybind11::error_already_set(); - } -#endif - } - - VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing neighbor features and their depends for all ROIs\n"); - reduce_neighbors_and_dependencies_manual (env); - - return true; - } - - void allocateTrivialRoisBuffers_3D (const std::vector& roi_labels, Roidata& roiData, CpusideCache & cache) - { - // Calculate the total memory demand (in # of items) of all segments' image matrices - cache.imageMatrixBufferLen = 0; - for (auto lab : roi_labels) - { - LR& r = roiData[lab]; - size_t w = r.aabb.get_width(), - h = r.aabb.get_height(), - d = r.aabb.get_z_depth(), - v = w * h * d; - cache.imageMatrixBufferLen += v; - - cache.largest_roi_imatr_buf_len = cache.largest_roi_imatr_buf_len == 0 ? v : std::max (cache.largest_roi_imatr_buf_len, v); - } - - // - // Preallocate image matrices and cubes here (in the future). - // - for (auto lab : roi_labels) - { - LR& r = roiData[lab]; - r.aux_image_cube.calculate_from_pixelcloud(r.raw_pixels_3D, r.aabb); - } - } - +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef WITH_PYTHON_H + #include +#endif + +#include "environment.h" +#include "globals.h" +#include "helpers/fsystem.h" +#include "helpers/timing.h" + +namespace Nyxus +{ + // + // Loads ROI voxels into voxel clouds + // + bool scanTrivialRois_3D (Environment & env, const std::vector& batch_labels, const std::string& intens_fpath, const std::string& label_fpath, size_t t_index) + { + // Sort the batch's labels to enable binary searching in it + std::vector whiteList = batch_labels; + std::sort(whiteList.begin(), whiteList.end()); + + // Scan this Z intensity-mask pair + SlideProps p (intens_fpath, label_fpath); + if (! env.theImLoader.open(p, env.fpimageOptions)) + { + std::cerr << "Error opening a file pair with ImageLoader. Terminating\n"; + return false; + } + + // thanks to ImageLoader::open() we are guaranteed that the mask and intensity volumes' + // width, height, and depth match. Mask and intensity may only differ in the number of + // time frames: 1:1, 1:N, and N:1 cases are permitted. + size_t + /* + nth = env.theImLoader.get_num_tiles_hor(), + ntv = env.theImLoader.get_num_tiles_vert(), + fw = env.theImLoader.get_tile_width(), + th = env.theImLoader.get_tile_height(), + tw = env.theImLoader.get_tile_width(), + tileSize = env.theImLoader.get_tile_size(), + */ + w = env.theImLoader.get_full_width(), + h = env.theImLoader.get_full_height(), + d = env.theImLoader.get_full_depth(), + sliceSize = w * h, + timeFrameSize = sliceSize * d, + timeI = env.theImLoader.get_inten_time(), + timeM = env.theImLoader.get_mask_time(), + nVoxI = timeFrameSize * timeI, + nVoxM = timeFrameSize * timeM; + + // is this intensity-mask pair's shape supported? + if (nVoxI < nVoxM) + { + std::string erm = "Error: unsupported shape - intensity file: " + std::to_string(nVoxI) + ", mask file: " + std::to_string(nVoxM); +#ifdef WITH_PYTHON_H + throw erm; +#endif + std::cerr << erm << "\n"; + return false; + } + + int cnt = 1; + + // fetch 3D data + bool ok = env.theImLoader.load_tile (0/*row*/, 0/*col*/); + if (!ok) + { + std::string erm = "Error fetching segmented data from " + intens_fpath + "(I) " + label_fpath + "(M)"; +#ifdef WITH_PYTHON_H + throw erm; +#endif + std::cerr << erm << "\n"; + return false; + } + + // Get ahold of tile's pixel buffer + auto dataI = env.theImLoader.get_int_tile_buffer(), + dataL = env.theImLoader.get_seg_tile_buffer(); + + size_t baseI, baseM; + if (timeI == timeM) + { + baseM = baseI = t_index * timeFrameSize; + } + else // nVoxI > nVoxM + { + baseM = 0; + baseI = t_index * timeFrameSize; + } + + // Iterate voxels + for (size_t i=0; i= w || y >= h || z >= d) + continue; + + // Collapse all the labels to one if single-ROI mde is requested + if (env.singleROI) + label = 1; + + // Cache this pixel + LR& r = env.roiData[label]; + feed_pixel_2_cache_3D_LR (x, y, z, dataI[j], r); + } + +#ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); +#endif + + // Close the image pair + env.theImLoader.close(); + + // Dump ROI pixel clouds to the output directory + VERBOSLVL5 (env.get_verbosity_level(), dump_roi_pixels(env.dim(), Nyxus::get_temp_dir_path(), batch_labels, label_fpath, env.uniqueLabels, env.roiData)); + + return true; + } + + bool scanTrivialRois_3D_anisotropic__BEFORE4D( + Environment & env, + const std::vector& batch_labels, + const std::string& intens_fpath, + const std::string& label_fpath, + double aniso_x, + double aniso_y, + double aniso_z) + { + // Sort the batch's labels to enable binary searching in it + std::vector whiteList = batch_labels; + std::sort(whiteList.begin(), whiteList.end()); + + int lvl = 0, // pyramid level + lyr = 0; // layer + + // Scan this Z intensity-mask pair + SlideProps p (intens_fpath, label_fpath); + if (! env.theImLoader.open(p, env.fpimageOptions)) + { + std::cerr << "Error opening a file pair with ImageLoader. Terminating\n"; + return false; + } + + size_t nth = env.theImLoader.get_num_tiles_hor(), + ntv = env.theImLoader.get_num_tiles_vert(), + fw = env.theImLoader.get_tile_width(), + th = env.theImLoader.get_tile_height(), + tw = env.theImLoader.get_tile_width(), + tileSize = env.theImLoader.get_tile_size(), + fullW = env.theImLoader.get_full_width(), + fullH = env.theImLoader.get_full_height(), + fullD = env.theImLoader.get_full_depth(); + + size_t vD = (size_t)(double(fullD) * aniso_z); // virtual depth + + for (size_t vz = 0; vz < vD; vz++) + { + size_t z = size_t(double(vz) / aniso_z); // physical z + + // virtual slide properties + size_t vh = (size_t)(double(fullH) * aniso_y), + vw = (size_t)(double(fullW) * aniso_x), + vth = (size_t)(double(th) * aniso_y), + vtw = (size_t)(double(tw) * aniso_x); + + // current tile to skip tile reloads + size_t curt_x = 999, curt_y = 999; + + for (size_t vr = 0; vr < vh; vr++) + { + for (size_t vc = 0; vc < vw; vc++) + { + // tile position + size_t tidx_y = size_t(vr / vth), + tidx_x = size_t(vc / vtw); + + // load it + if (tidx_y != curt_y || tidx_x != curt_x) + { + bool ok = env.theImLoader.load_tile(tidx_y, tidx_x); + if (!ok) + { + std::string s = "Error fetching tile row=" + std::to_string(tidx_y) + " col=" + std::to_string(tidx_x); +#ifdef WITH_PYTHON_H + throw s; +#endif + std::cerr << s << "\n"; + return false; + } + + // cache tile position to avoid reloading + curt_y = tidx_y; + curt_x = tidx_x; + } + + // within-tile virtual pixel position + size_t vx = vc - tidx_x * vtw, + vy = vr - tidx_y * vth; + + // within-tile physical pixel position + size_t ph_x = size_t(double(vx) / aniso_x), + ph_y = size_t(double(vy) / aniso_y), + i = ph_y * tw + ph_x; + + // read buffered physical pixel + auto dataI = env.theImLoader.get_int_tile_buffer(), + dataL = env.theImLoader.get_seg_tile_buffer(); + + // skip non-mask pixels + auto label = dataL[i]; + if (!label) + continue; + + // skip this ROI if the label isn't in the pending set of a multi-ROI mode + if (!env.singleROI && !std::binary_search(whiteList.begin(), whiteList.end(), label)) + continue; + + // skip tile buffer pixels beyond the image's bounds + if (vc >= fullW || vr >= fullH) + continue; + + // collapse all the labels to one if single-ROI mde is requested + if (env.singleROI) + label = 1; + + // cache this voxel + auto inten = dataI[i]; + LR& r = env.roiData[label]; + feed_pixel_2_cache_3D_LR (vc, vr, vz, dataI[i], r); + } + } + + // Close the image pair + env.theImLoader.close(); + } + + // Dump ROI pixel clouds to the output directory + VERBOSLVL5 (env.get_verbosity_level(), dump_roi_pixels(env.dim(), Nyxus::get_temp_dir_path(), batch_labels, label_fpath, env.uniqueLabels, env.roiData)); + + return true; + } + + // + // Loads ROI voxels into voxel clouds + // + 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) + { + // Sort the batch's labels to enable binary searching in it + std::vector whiteList = batch_labels; + std::sort (whiteList.begin(), whiteList.end()); + + // temp slideprops instance to pass some into to ImageLoader + SlideProps p (intens_fpath, label_fpath); + if (!env.theImLoader.open(p, env.fpimageOptions)) + { + std::cerr << "Error opening a file pair with ImageLoader. Terminating\n"; + return false; + } + + // thanks to ImageLoader::open() we are guaranteed that the mask's and intensity's + // W, H, and D match. Mask and intensity may only differ in the number of + // time frames: 1:1, 1:N, and N:1 cases are permitted + size_t + w = env.theImLoader.get_full_width(), + h = env.theImLoader.get_full_height(), + d = env.theImLoader.get_full_depth(), + slice = w * h, + timeFrameSize = slice * d, + timeI = env.theImLoader.get_inten_time(), + timeM = env.theImLoader.get_mask_time(), + nVoxI = timeFrameSize * timeI, + nVoxM = timeFrameSize * timeM; + + // is this intensity-mask pair's shape supported? + if (nVoxI < nVoxM) + { + std::string erm = "Error: unsupported shape - intensity file: " + std::to_string(nVoxI) + ", mask file: " + std::to_string(nVoxM); + #ifdef WITH_PYTHON_H + throw erm; + #endif + std::cerr << erm << "\n"; + return false; + } + + int cnt = 1; + + // fetch 3D data + if (!env.theImLoader.load_tile (0/*row*/, 0/*col*/)) + { + std::string erm = "Error fetching data from file pair " + intens_fpath + "(I) " + label_fpath + "(M)"; + #ifdef WITH_PYTHON_H + throw erm; + #endif + std::cerr << erm << "\n"; + return false; + } + + // get ahold of voxel buffers + auto dataI = env.theImLoader.get_int_tile_buffer(), + dataL = env.theImLoader.get_seg_tile_buffer(); + + // align time frame's mask and intensity volumes + size_t baseI, baseM; + if (timeI == timeM) + { + // trivial N mask : N intensity + baseM = + baseI = t_index * timeFrameSize; + } + else + { + // nontrivial 1 mask : N intensity + baseM = 0; + baseI = t_index * timeFrameSize; + } + + // virtual dimensions + size_t virt_h = h * aniso_y, + virt_w = w * aniso_x, + virt_d = d * aniso_z; + size_t vSliceLen = virt_h * virt_w, + virt_v = vSliceLen * virt_d; + + // iterate virtual voxels and fill them with corresponding physical intensities + for (size_t vIdx = 0; vIdx < virt_v; vIdx++) + { + // virtual Cartesian position + size_t vZ = vIdx / vSliceLen, + vLastSliceLen = vIdx % vSliceLen, + vY = vLastSliceLen / virt_w, + vX = vLastSliceLen % virt_w; + + // physical Cartesian position + size_t pZ = vZ / aniso_z + 0.5, + pY = vY / aniso_y + 0.5, + pX = vX / aniso_x + 0.5; + + // skip a position outside the bounds + // (since we are casting coorinates from virtual to physical, + // we may get positions outside the physical bounds) + if (pX >= w || pY >= h || pZ >= d) + continue; + + // physical offset + size_t i = pZ * slice + pY * w + pX; + + // + // interpret the mask intensity + // + + // skip non-mask pixels + auto lbl = dataL[baseM + i]; + if (!lbl) + continue; + + // skip this ROI if the label isn't in the pending set of a multi-ROI mode + if (!env.singleROI && !std::binary_search(whiteList.begin(), whiteList.end(), lbl)) + continue; + + // collapse all the labels to one if single-ROI mde is requested + if (env.singleROI) + lbl = 1; + +#if !defined(NDEBUG) + if (vZ >= virt_d) + std::cout << "vZ=" << vZ << " < virt_d =" << virt_d << "\n"; + assert(vZ < virt_d); + if (vY >= virt_h) + std::cout << "vY=" << vY << " < virt_h =" << virt_h << "\n"; + assert(vY < virt_h); + if (vX >= virt_w) + std::cout << "vX=" << vX << " < virt_w =" << virt_w << "\n"; + assert(vX < virt_w); +#endif + + // cache this voxel + auto inten = dataI[baseI + i]; + LR& r = env.roiData[lbl]; + feed_pixel_2_cache_3D_LR (vX, vY, vZ, inten, r); + } + + #ifdef WITH_PYTHON_H + // allow keyboard interrupt + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); + #endif + + // Close the image pair + env.theImLoader.close(); + return true; + } + + // + // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' + // + bool scan_trivial_wholevolume ( + LR& vroi, + const std::string& intens_fpath, + ImageLoader& ilo) + { + int lvl = 0, // Pyramid level + lyr = 0; // Layer + + // Read the tiffs + + size_t fullwidth = ilo.get_full_width(), + fullheight = ilo.get_full_height(), + fullD = ilo.get_full_depth(), + sliceSize = fullwidth * fullheight, + nVox = sliceSize * fullD; + + // in the 3D case tiling is a formality, so fetch the only tile in the file + if (!ilo.load_tile(0, 0)) + { +#ifdef WITH_PYTHON_H + throw "Error fetching tile"; +#endif + std::cerr << "Error fetching tile\n"; + return false; + } + + // Get ahold of tile's pixel buffer + const std::vector& dataI = ilo.get_int_tile_buffer(); + + // iterate abstract tiles (in a tiled slide /e.g. tiled tiff/ they correspond to physical tiles, in a nontiled slide /e.g. scanline tiff or strip tiff/ they correspond to ) + int cnt = 1; + + // iterate voxels + for (size_t i = 0; i < nVox; i++) + { + int z = i / sliceSize, + y = (i - z * sliceSize) / fullwidth, + x = (i - z * sliceSize) % fullwidth; + + // Skip tile buffer pixels beyond the image's bounds + if (x >= fullwidth || y >= fullheight || z >= fullD) + continue; + + // dynamic range within- and off-ROI + auto inten = dataI[i]; + + // Cache this pixel + feed_pixel_2_cache_3D_LR (x, y, z, inten, vroi); + + } //- all voxels + + return true; + } + + // + // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' + // + bool scan_trivial_wholevolume_anisotropic ( + LR& vroi, + const std::string& intens_fpath, + ImageLoader& ilo, + double aniso_x, + double aniso_y, + double aniso_z) + { + int lvl = 0, // Pyramid level + lyr = 0; // Layer + + // Read the tiffs + + size_t fullW = ilo.get_full_width(), + fullH = ilo.get_full_height(), + fullD = ilo.get_full_depth(), + sliceSize = fullW * fullH; + + size_t vh = (size_t) (double(fullH) * aniso_y), + vw = (size_t) (double(fullW) * aniso_x), + vd = (size_t) (double(fullD) * aniso_z); + + // in the 3D case tiling is a formality, so fetch the only tile in the file + if (! ilo.load_tile(0, 0)) + { +#ifdef WITH_PYTHON_H + throw "Error loading volume data"; +#endif + std::cerr << "Error loading volume data\n"; + return false; + } + + // Get ahold of tile's pixel buffer + const std::vector& dataI = ilo.get_int_tile_buffer(); + + // iterate virtual voxels + size_t vSliceSize = vh * vw, + nVox = vh * vw * vd; + for (size_t i = 0; i < nVox; i++) + { + // virtual voxel position + int z = i / vSliceSize, + y = (i - z * vSliceSize) / vw, + x = (i - z * vSliceSize) % vw; + + // physical voxel position + size_t ph_x = (size_t) (double(x) / aniso_x), + ph_y = (size_t) (double(y) / aniso_y), + ph_z = (size_t) (double(z) / aniso_z); + i = ph_z * sliceSize + ph_y * fullH + ph_x; + + // Cache this pixel + feed_pixel_2_cache_3D_LR (x, y, z, dataI[i], vroi); + + } + + return true; + } + + // + // Reads pixels of whole slide 'intens_fpath' into virtual ROI 'vroi' + // performing anisotropy correction + // + bool scan_trivial_wholevolume_anisotropic__OLD( + LR& vroi, + const std::string& intens_fpath, + ImageLoader& ldr, + double aniso_x, + double aniso_y) + { + int lvl = 0, // Pyramid level + lyr = 0; // Layer + + // physical slide properties + size_t nth = ldr.get_num_tiles_hor(), + ntv = ldr.get_num_tiles_vert(), + fw = ldr.get_tile_width(), + th = ldr.get_tile_height(), + tw = ldr.get_tile_width(), + tileSize = ldr.get_tile_size(), + fullwidth = ldr.get_full_width(), + fullheight = ldr.get_full_height(); + + // virtual slide properties + size_t vh = (size_t)(double(fullheight) * aniso_y), + vw = (size_t)(double(fullwidth) * aniso_x), + vth = (size_t)(double(th) * aniso_y), + vtw = (size_t)(double(tw) * aniso_x); + + // current tile to skip tile reloads + size_t curt_x = 999, curt_y = 999; + + for (size_t vr = 0; vr < vh; vr++) + { + for (size_t vc = 0; vc < vw; vc++) + { + // tile position for virtual pixel (vc, vr) + size_t tidx_y = size_t(vr / vth), + tidx_x = size_t(vc / vtw); + + // load it + if (tidx_y != curt_y || tidx_x != curt_x) + { + bool ok = ldr.load_tile(tidx_y, tidx_x); + if (!ok) + { + std::string s = "Error fetching tile row=" + std::to_string(tidx_y) + " col=" + std::to_string(tidx_x); +#ifdef WITH_PYTHON_H + throw s; +#endif + std::cerr << s << "\n"; + return false; + } + + // cache tile position to avoid reloading + curt_y = tidx_y; + curt_x = tidx_x; + } + + // within-tile virtual pixel position + size_t vx = vc - tidx_x * vtw, + vy = vr - tidx_y * vth; + + // within-tile physical pixel position + size_t ph_x = size_t(double(vx) / aniso_x), + ph_y = size_t(double(vy) / aniso_y), + i = ph_y * tw + ph_x; + + // read buffered physical pixel + const std::vector& dataI = ldr.get_int_tile_buffer(); + const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); + bool wholeslide = env.singleROI; + + // Cache this pixel + feed_pixel_2_cache_LR(vc, vr, dataI[i], vroi); + } + } + + return true; + } + + + bool processTrivialRois_3D (Environment & env, size_t sidx, size_t t_index, const std::vector& trivRoiLabels, const std::string& intens_fpath, const std::string& label_fpath, size_t memory_limit) + { + std::vector Pending; + size_t batchDemand = 0; + int roiBatchNo = 1; + + for (auto lab : trivRoiLabels) + { + LR& r = env.roiData[lab]; + + size_t itemFootprint = r.get_ram_footprint_estimate (Pending.size()); + + // Check if we are good to accumulate this ROI in the current batch or should close the batch and reduce it + if (batchDemand + itemFootprint < memory_limit) + { + Pending.push_back(lab); + batchDemand += itemFootprint; + } + else + { + // Scan pixels of pending trivial ROIs + std::sort (Pending.begin(), Pending.end()); + + VERBOSLVL2 (env.get_verbosity_level(), std::cout << ">>> Scanning batch #" << roiBatchNo << " of " << Pending.size() << " pending ROIs of total " << env.uniqueLabels.size() << " ROIs\n"); + VERBOSLVL2 (env.get_verbosity_level(), + if (Pending.size() == 1) + std::cout << ">>> (single ROI label " << Pending[0] << ")\n"; + else + std::cout << ">>> (ROI labels " << Pending[0] << " ... " << Pending[Pending.size() - 1] << ")\n"; + ); + + if (env.anisoOptions.customized() == false) + { + scanTrivialRois_3D (env, Pending, 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, Pending, intens_fpath, label_fpath, t_index, ax, ay, az); + + // rescan and update ROI's AABB + for (auto lbl : Pending) + { + LR& r = env.roiData[lbl]; + r.aabb.update_from_voxelcloud (r.raw_pixels_3D); + } + } + + // Allocate memory + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tallocating ROI buffers\n";) + allocateTrivialRoisBuffers_3D (Pending, env.roiData, env.hostCache); + + // Reduce them + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing ROIs\n";) + // reduce_trivial_rois(Pending); + reduce_trivial_rois_manual (Pending, env); + + // Free memory + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tfreeing ROI buffers\n";) + freeTrivialRoisBuffers_3D (Pending, env.roiData); // frees what's allocated by feed_pixel_2_cache() and allocateTrivialRoisBuffers() + + // Reset the RAM footprint accumulator + batchDemand = 0; + + // Clear the freshly processed ROIs from pending list + Pending.clear(); + + // Start a new pending set by adding the stopper ROI + Pending.push_back(lab); + + // Advance the batch counter + roiBatchNo++; + } + + // Allow keyboard interrupt +#ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + { + sureprint("\nAborting per user input\n"); + throw pybind11::error_already_set(); + } +#endif + } + + // Process what's remaining pending + if (Pending.size() > 0) + { + // Read raw pixels of pending trivial ROIs + std::sort(Pending.begin(), Pending.end()); + + VERBOSLVL2(env.get_verbosity_level(), + std::cout << ">>> Scanning batch #" << roiBatchNo << " of " << Pending.size() << "(" << env.uniqueLabels.size() << ") ROIs\n"; + std::cout << ">>> (labels " << Pending[0] << " ... " << Pending[Pending.size() - 1] << ")\n"; + ); + + if (env.anisoOptions.customized() == false) + { + scanTrivialRois_3D (env, Pending, 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, Pending, intens_fpath, label_fpath, t_index, ax, ay, az); + + // rescan and update ROI's AABB + for (auto lbl : Pending) + { + LR& r = env.roiData[lbl]; + r.aabb.update_from_voxelcloud(r.raw_pixels_3D); + } + } + + for (auto lab : Pending) + { + LR& r = env.roiData[lab]; + for (Pixel3& vox : r.raw_pixels_3D) + { + assert (vox.x >= r.aabb.get_xmin()); + assert (vox.x <= r.aabb.get_xmax()); + assert (vox.y >= r.aabb.get_ymin()); + assert (vox.y <= r.aabb.get_ymax()); + assert (vox.z >= r.aabb.get_zmin()); + assert (vox.z <= r.aabb.get_zmax()); + } + } + + // Allocate memory + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tallocating ROI buffers\n"); + allocateTrivialRoisBuffers_3D (Pending, env.roiData, env.hostCache); + + // Dump ROIs for use in unit testing +#ifdef DUMP_ALL_ROI + dump_all_roi(); +#endif + + // Reduce them + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing ROIs\n"); + //reduce_trivial_rois(Pending); + reduce_trivial_rois_manual (Pending, env); + + // Free memory + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\tfreeing ROI buffers\n"); + freeTrivialRoisBuffers_3D (Pending, env.roiData); + +#ifdef WITH_PYTHON_H + // Allow keyboard interrupt + if (PyErr_CheckSignals() != 0) + { + sureprint("\nAborting per user input\n"); + throw pybind11::error_already_set(); + } +#endif + } + + VERBOSLVL2 (env.get_verbosity_level(), std::cout << "\treducing neighbor features and their depends for all ROIs\n"); + reduce_neighbors_and_dependencies_manual (env); + + return true; + } + + void allocateTrivialRoisBuffers_3D (const std::vector& roi_labels, Roidata& roiData, CpusideCache & cache) + { + // Calculate the total memory demand (in # of items) of all segments' image matrices + cache.imageMatrixBufferLen = 0; + for (auto lab : roi_labels) + { + LR& r = roiData[lab]; + size_t w = r.aabb.get_width(), + h = r.aabb.get_height(), + d = r.aabb.get_z_depth(), + v = w * h * d; + cache.imageMatrixBufferLen += v; + + cache.largest_roi_imatr_buf_len = cache.largest_roi_imatr_buf_len == 0 ? v : std::max (cache.largest_roi_imatr_buf_len, v); + } + + // + // Preallocate image matrices and cubes here (in the future). + // + for (auto lab : roi_labels) + { + LR& r = roiData[lab]; + r.aux_image_cube.calculate_from_pixelcloud(r.raw_pixels_3D, r.aabb); + } + } + } \ No newline at end of file diff --git a/src/nyx/python/new_bindings_py.cpp b/src/nyx/python/new_bindings_py.cpp index 8a377ba3..e26cb243 100644 --- a/src/nyx/python/new_bindings_py.cpp +++ b/src/nyx/python/new_bindings_py.cpp @@ -643,28 +643,26 @@ py::tuple featurize_fname_lists_imp (uint64_t instid, const py::list& int_fnames labelFiles.push_back(fn); } - // Check the file names + // Check the file names if (intensFiles.size() == 0) throw std::runtime_error("Intensity file list is blank"); - if (labelFiles.size() == 0) - throw std::runtime_error("Segmentation mask file list is blank"); - if (intensFiles.size() != labelFiles.size()) - throw std::runtime_error("Imbalanced intensity and segmentation mask file lists"); - for (auto i = 0; i < intensFiles.size(); i++) + + if (single_roi) + labelFiles.assign(intensFiles.size(), ""); // wholeslide — ignore any provided seg files + else { - const std::string& i_fname = intensFiles[i]; - const std::string& s_fname = labelFiles[i]; + if (labelFiles.size() == 0) + throw std::runtime_error("Segmentation mask file list is blank"); + if (intensFiles.size() != labelFiles.size()) + throw std::runtime_error("Imbalanced intensity and segmentation mask file lists"); + } - if (!existsOnFilesystem(i_fname)) - { - auto msg = "File does not exist: " + i_fname; - throw std::runtime_error(msg); - } - if (!existsOnFilesystem(s_fname)) - { - auto msg = "File does not exist: " + s_fname; - throw std::runtime_error(msg); - } + for (size_t i = 0; i < intensFiles.size(); i++) + { + if (!existsOnFilesystem(intensFiles[i])) + throw std::runtime_error("File does not exist: " + intensFiles[i]); + if (!labelFiles[i].empty() && !existsOnFilesystem(labelFiles[i])) + throw std::runtime_error("File does not exist: " + labelFiles[i]); } theEnvironment.theResultsCache.clear(); diff --git a/src/nyx/python/nyxus/nyxus.py b/src/nyx/python/nyxus/nyxus.py index a2725e07..bc3a8224 100644 --- a/src/nyx/python/nyxus/nyxus.py +++ b/src/nyx/python/nyxus/nyxus.py @@ -525,15 +525,17 @@ def featurize_files ( if intensity_files is None: raise IOError ("The list of intensity file paths is empty") - if mask_files is None: - raise IOError ("The list of segment file paths is empty") - + if mask_files is None and not single_roi: + raise IOError ("The list of segment file paths is empty. Supply mask images or set single_roi to True") + 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}') + mask_files_arg = mask_files if mask_files is not None else [] + if (output_type == 'pandas'): - header, string_data, numeric_data = featurize_fname_lists_imp (id(self), intensity_files, mask_files, single_roi, output_type, "") + header, string_data, numeric_data = featurize_fname_lists_imp (id(self), intensity_files, mask_files_arg, single_roi, output_type, "") df = pd.concat( [ @@ -548,10 +550,10 @@ 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) + featurize_fname_lists_imp (id(self), intensity_files, mask_files_arg, single_roi, output_type, output_path) return get_arrow_file_imp (id(self)) diff --git a/src/nyx/slideprops.cpp b/src/nyx/slideprops.cpp index 0881970b..41c97c59 100644 --- a/src/nyx/slideprops.cpp +++ b/src/nyx/slideprops.cpp @@ -1,495 +1,494 @@ -#include -#include -#include -#include "globals.h" -#include "helpers/fsystem.h" -#include "helpers/timing.h" -#include "raw_image_loader.h" - -namespace Nyxus -{ - bool gatherRoisMetrics_2_slideprops_2D_montage ( - // in - const AnisotropyOptions& aniso, - // out - SlideProps& p) - { - // low-level slide properties (intensity and mask, if available) - p.lolvl_slide_descr = "from montage"; - p.fp_phys_pivoxels = false; - - // time series - p.inten_time = 0; - p.mask_time = 0; - - // scan intensity slide's data - - bool wholeslide = false; - - double slide_I_max = (std::numeric_limits::lowest)(), - slide_I_min = (std::numeric_limits::max)(); - - std::unordered_set U; // unique ROI mask labels - std::unordered_map R; // ROI data - - //****** fix ROIs' AABBs with respect to anisotropy - - if (!aniso.customized()) - { - for (auto& pair : R) - { - LR& r = pair.second; - r.make_nonanisotropic_aabb(); - } - } - else - { - for (auto& pair : R) - { - LR& r = pair.second; - r.make_anisotropic_aabb(aniso.get_aniso_x(), aniso.get_aniso_y()); - } - } - - //****** Analysis - - // slide-wide (max ROI area) x (number of ROIs) - size_t maxArea = 0; - size_t max_w = 0, max_h = 0; - for (const auto& pair : R) - { - const LR& r = pair.second; - maxArea = maxArea > r.aux_area ? maxArea : r.aux_area; //std::max (maxArea, r.aux_area); - const AABB& bb = r.aabb; - auto w = bb.get_width(); - auto h = bb.get_height(); - max_w = max_w > w ? max_w : w; - max_h = max_h > h ? max_h : h; - } - - p.slide_w = 2; - p.slide_h = 2; - - p.max_preroi_inten = slide_I_max; // in case fp_phys_pivoxels==true, max/min _preroi_inten - p.min_preroi_inten = slide_I_min; // needs adjusting (grey-binning) before using in wsi scenarios (assigning ROI's min and max) - - p.max_roi_area = maxArea; - p.n_rois = R.size(); - p.max_roi_w = max_w; - p.max_roi_h = max_h; - - return true; - } - - bool gatherRoisMetrics_2_slideprops_2D( - // in - RawImageLoader& ilo, - const AnisotropyOptions& aniso, - // out - SlideProps& p) - { - // low-level slide properties (intensity and mask, if available) - p.lolvl_slide_descr = ilo.get_slide_descr(); - p.fp_phys_pivoxels = ilo.get_fp_phys_pixvoxels(); - - // time series - p.inten_time = ilo.get_inten_time(); - p.mask_time = ilo.get_mask_time(); - - // scan intensity slide's data - - bool wholeslide = p.fname_seg.empty(); - - double slide_I_max = (std::numeric_limits::lowest)(), - slide_I_min = (std::numeric_limits::max)(); - - std::unordered_set U; // unique ROI mask labels - std::unordered_map R; // ROI data - - int lvl = 0, // pyramid level - lyr = 0; // layer - - // Read the image/volume. The image loader is put in the open state in processDataset_XX_YY () - size_t nth = ilo.get_num_tiles_hor(), - ntv = ilo.get_num_tiles_vert(), - fw = ilo.get_tile_width(), - th = ilo.get_tile_height(), - tw = ilo.get_tile_width(), - tileSize = ilo.get_tile_size(), - fullwidth = ilo.get_full_width(), - fullheight = ilo.get_full_height(); - - // iterate abstract tiles (in a tiled slide /e.g. tiled tiff/ they correspond to physical tiles, in a nontiled slide /e.g. scanline tiff or strip tiff/ they correspond to ) - int cnt = 1; - for (unsigned int row = 0; row < nth; row++) - for (unsigned int col = 0; col < ntv; col++) - { - // Fetch the tile - if (!ilo.load_tile(row, col)) - { -#ifdef WITH_PYTHON_H - throw "Error fetching tile"; -#endif - std::cerr << "Error fetching tile\n"; - return false; - } - - // Get ahold of tile's pixel buffer - auto tidx = row * nth + col; - - // Iterate pixels - for (size_t i = 0; i < tileSize; i++) - { - // Mask - uint32_t msk = 1; // wholeslide by default - if (!wholeslide) - msk = ilo.get_cur_tile_seg_pixel(i); - - // Skip non-mask pixels - if (!msk) - continue; - - int y = row * th + i / tw, - x = col * tw + i % tw; - - // Skip tile buffer pixels beyond the image's bounds - if (x >= fullwidth || y >= fullheight) - continue; - - // dynamic range within- and off-ROI - double dxequiv_I = ilo.get_cur_tile_dpequiv_pixel(i); - slide_I_max = (std::max)(slide_I_max, dxequiv_I); - slide_I_min = (std::min)(slide_I_min, dxequiv_I); - - // Update pixel's ROI metrics - // - the following block mocks feed_pixel_2_metrics (x, y, dataI[i], msk, tidx) - if (U.find(msk) == U.end()) - { - // Remember this label - U.insert(msk); - - // Initialize the ROI label record - LR r(msk); - - // - mocking init_label_record_3 (roi, theSegFname, theIntFname, x, y, label, intensity, tile_index) - // Initialize basic counters - r.aux_area = 1; - r.aux_min = r.aux_max = 0; //we don't have uint-cast intensities at this moment - r.init_aabb(x, y); - - // - not storing file names (r.segFname = segFile, r.intFname = intFile) but will do so in the future - - // Attach - R[msk] = r; - } - else - { - // Update basic ROI info (info that doesn't require costly calculations) - LR& r = R[msk]; - - // - mocking update_label_record_2 (r, x, y, label, intensity, tile_index) - - // Per-ROI - r.aux_area++; - - // save - r.update_aabb(x, y); - } - } // scan tile - -#ifdef WITH_PYTHON_H - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); -#endif - - } // foreach tile - - //****** fix ROIs' AABBs with respect to anisotropy - - if (!aniso.customized()) - { - for (auto& pair : R) - { - LR& r = pair.second; - r.make_nonanisotropic_aabb(); - } - } - else - { - for (auto& pair : R) - { - LR& r = pair.second; - r.make_anisotropic_aabb(aniso.get_aniso_x(), aniso.get_aniso_y()); - } - } - - //****** Analysis - - // slide-wide (max ROI area) x (number of ROIs) - size_t maxArea = 0; - size_t max_w = 0, max_h = 0; - for (const auto& pair : R) - { - const LR& r = pair.second; - maxArea = maxArea > r.aux_area ? maxArea : r.aux_area; //std::max (maxArea, r.aux_area); - const AABB& bb = r.aabb; - auto w = bb.get_width(); - auto h = bb.get_height(); - max_w = max_w > w ? max_w : w; - max_h = max_h > h ? max_h : h; - } - - p.slide_w = fullwidth; - p.slide_h = fullheight; - - p.max_preroi_inten = slide_I_max; // in case fp_phys_pivoxels==true, max/min _preroi_inten - p.min_preroi_inten = slide_I_min; // needs adjusting (grey-binning) before using in wsi scenarios (assigning ROI's min and max) - - p.max_roi_area = maxArea; - p.n_rois = R.size(); - p.max_roi_w = max_w; - p.max_roi_h = max_h; - - return true; - } - - bool gatherRoisMetrics_2_slideprops_3D( - // in - RawImageLoader& ilo, - const AnisotropyOptions& aniso, - // out - SlideProps& p) - { - // low-level slide properties (intensity and mask, if available) - p.lolvl_slide_descr = ilo.get_slide_descr(); - p.fp_phys_pivoxels = ilo.get_fp_phys_pixvoxels(); - - // time series - p.inten_time = ilo.get_inten_time(); - p.mask_time = ilo.get_mask_time(); - - // scan intensity slide's data - - bool wholeslide = p.fname_seg.empty(); - - double slide_I_max = (std::numeric_limits::lowest)(), - slide_I_min = (std::numeric_limits::max)(); - - std::unordered_set U; // unique ROI mask labels - std::unordered_map R; // ROI data - - // Read the volume. The image loader is in the open state by previously called processDataset_XX_YY () - size_t fullW = ilo.get_full_width(), - fullH = ilo.get_full_height(), - fullD = ilo.get_full_depth(), - sliceSize = fullW * fullH, - nVox = sliceSize * fullD; - - // in the 3D case tiling is a formality, so fetch the only tile in the file - if (!ilo.load_tile(0, 0)) - { -#ifdef WITH_PYTHON_H - throw "Error fetching tile"; -#endif - std::cerr << "Error fetching tile\n"; - return false; - } - - // iterate abstract tiles (in a tiled slide /e.g. tiled tiff/ they correspond to physical tiles, in a nontiled slide /e.g. scanline tiff or strip tiff/ they correspond to ) - int cnt = 1; - - // iterate voxels - for (size_t i = 0; i < nVox; i++) - { - // Mask - uint32_t msk = 1; // wholeslide by default - if (!wholeslide) - msk = ilo.get_cur_tile_seg_pixel(i); - - // Skip non-mask pixels - if (!msk) - continue; - - int z = i / sliceSize, - y = (i - z * sliceSize) / fullW, - x = (i - z * sliceSize) % fullW; - - // Skip tile buffer pixels beyond the image's bounds - if (x >= fullW || y >= fullH || z >= fullD) - continue; - - // dynamic range within- and off-ROI - double dxequiv_I = ilo.get_cur_tile_dpequiv_pixel(i); - slide_I_max = (std::max)(slide_I_max, dxequiv_I); - slide_I_min = (std::min)(slide_I_min, dxequiv_I); - - // Update pixel's ROI metrics - // - the following block mocks feed_pixel_2_metrics (x, y, dataI[i], msk, tidx) - if (U.find(msk) == U.end()) - { - // Remember this label - U.insert(msk); - - // Initialize the ROI label record - LR r(msk); - - // - mocking init_label_record_3 (newData, theSegFname, theIntFname, x, y, label, intensity, tile_index) - // Initialize basic counters - r.aux_area = 1; - r.aux_min = r.aux_max = 0; //we don't have uint-cast intensities at this moment - r.init_aabb_3D(x, y, z); - - // - not storing file names (r.segFname = segFile, r.intFname = intFile) but will do so in the future - - // Attach - R[msk] = r; - } - else - { - // Update basic ROI info (info that doesn't require costly calculations) - LR& r = R[msk]; - - // - mocking update_label_record_2 (r, x, y, label, intensity, tile_index) - - // Per-ROI - r.aux_area++; - - // save - r.update_aabb_3D(x, y, z); - } - -#ifdef WITH_PYTHON_H - // keyboard interrupt - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); -#endif - - } //- all voxels - - //****** fix ROIs' AABBs with respect to anisotropy - - if (!aniso.customized()) - { - for (auto& pair : R) - { - LR& r = pair.second; - r.make_nonanisotropic_aabb(); - } - } - else - { - for (auto& pair : R) - { - LR& r = pair.second; - r.make_anisotropic_aabb(aniso.get_aniso_x(), aniso.get_aniso_y(), aniso.get_aniso_z()); - } - } - - //****** Analysis - - // slide-wide (max ROI area) x (number of ROIs) - size_t maxArea = 0; - size_t max_w = 0, max_h = 0, max_d = 0; - for (const auto& pair : R) - { - const LR& r = pair.second; - maxArea = (std::max)(maxArea, (size_t)r.aux_area); - const AABB& bb = r.aabb; - auto w = bb.get_width(); - auto h = bb.get_height(); - auto d = bb.get_z_depth(); - max_w = (std::max)(max_w, (size_t)w); - max_h = (std::max)(max_h, (size_t)h); - max_d = (std::max)(max_d, (size_t)d); - } - - p.slide_w = fullW; - p.slide_h = fullH; - p.volume_d = fullD; - - p.max_preroi_inten = slide_I_max; // in case fp_phys_pivoxels==true, max/min _preroi_inten - p.min_preroi_inten = slide_I_min; // needs adjusting (grey-binning) before using in wsi scenarios (assigning ROI's min and max) - - p.max_roi_area = maxArea; - p.n_rois = R.size(); - p.max_roi_w = max_w; - p.max_roi_h = max_h; - p.max_roi_d = max_d; - - return true; - } - - std::pair split_alnum(const std::string& annot) - { - std::string A = annot; // a string that we can edit - std::string al; - for (auto c : A) - { - if (!std::isdigit(c)) - al += c; - else - { - A.erase(0, al.size()); - break; - } - } - - return { al, A }; - } - - bool scan_slide_props_montage (SlideProps & p, int dim, const AnisotropyOptions & aniso) - { - if (dim != 2) - return false; - - gatherRoisMetrics_2_slideprops_2D_montage (aniso, p); - - return true; - } - - // - // prerequisite: initialized fields fname_int and fname_seg - // - bool scan_slide_props (SlideProps & p, int dim, const AnisotropyOptions & aniso, bool need_annot) - { - RawImageLoader ilo; - if (! ilo.open(p.fname_int, p.fname_seg)) - { - std::cerr << "error opening an ImageLoader for " << p.fname_int << " | " << p.fname_seg << "\n"; - return false; - } - - bool ok = dim==2 ? gatherRoisMetrics_2_slideprops_2D(ilo, aniso, p) : gatherRoisMetrics_2_slideprops_3D(ilo, aniso, p); - if (!ok) - { - std::cerr << "error gathering ROI metrics to slide/volume props \n"; - return false; - } - - ilo.close(); - - // annotations - if (need_annot) - { - // throw away the directory part - fs::path pth(p.fname_seg); - auto purefn = pth.filename().string(); - - // LHS part till the 1st dot is the annotation info that we need to parse - std::vectortoks1; - Nyxus::parse_delimited_string (purefn, ".", toks1); - - // the result tokens is the annotation info that we need - p.annots.clear(); - Nyxus::parse_delimited_string (toks1[0], "_", p.annots); - - // prune blank annotations (usually caused by multiple separator e.g. '__' in 'slide_blah1_bla_2__something3.ome.tiff') - p.annots.erase ( - std::remove_if (p.annots.begin(), p.annots.end(), [](const std::string & s) { return s.empty(); }), - p.annots.end()); - } - - - return true; - } +#include +#include +#include +#include "globals.h" +#include "helpers/fsystem.h" +#include "helpers/timing.h" +#include "raw_image_loader.h" + +namespace Nyxus +{ + bool gatherRoisMetrics_2_slideprops_2D_montage ( + // in + const AnisotropyOptions& aniso, + // out + SlideProps& p) + { + // low-level slide properties (intensity and mask, if available) + p.lolvl_slide_descr = "from montage"; + p.fp_phys_pivoxels = false; + + // time series + p.inten_time = 0; + p.mask_time = 0; + + // scan intensity slide's data + + bool wholeslide = false; + + double slide_I_max = (std::numeric_limits::lowest)(), + slide_I_min = (std::numeric_limits::max)(); + + std::unordered_set U; // unique ROI mask labels + std::unordered_map R; // ROI data + + //****** fix ROIs' AABBs with respect to anisotropy + + if (!aniso.customized()) + { + for (auto& pair : R) + { + LR& r = pair.second; + r.make_nonanisotropic_aabb(); + } + } + else + { + for (auto& pair : R) + { + LR& r = pair.second; + r.make_anisotropic_aabb(aniso.get_aniso_x(), aniso.get_aniso_y()); + } + } + + //****** Analysis + + // slide-wide (max ROI area) x (number of ROIs) + size_t maxArea = 0; + size_t max_w = 0, max_h = 0; + for (const auto& pair : R) + { + const LR& r = pair.second; + maxArea = maxArea > r.aux_area ? maxArea : r.aux_area; //std::max (maxArea, r.aux_area); + const AABB& bb = r.aabb; + auto w = bb.get_width(); + auto h = bb.get_height(); + max_w = max_w > w ? max_w : w; + max_h = max_h > h ? max_h : h; + } + + p.slide_w = 2; + p.slide_h = 2; + + p.max_preroi_inten = slide_I_max; // in case fp_phys_pivoxels==true, max/min _preroi_inten + p.min_preroi_inten = slide_I_min; // needs adjusting (grey-binning) before using in wsi scenarios (assigning ROI's min and max) + + p.max_roi_area = maxArea; + p.n_rois = R.size(); + p.max_roi_w = max_w; + p.max_roi_h = max_h; + + return true; + } + + bool gatherRoisMetrics_2_slideprops_2D( + // in + RawImageLoader& ilo, + const AnisotropyOptions& aniso, + // out + SlideProps& p) + { + // low-level slide properties (intensity and mask, if available) + p.lolvl_slide_descr = ilo.get_slide_descr(); + p.fp_phys_pivoxels = ilo.get_fp_phys_pixvoxels(); + + // time series + p.inten_time = ilo.get_inten_time(); + p.mask_time = ilo.get_mask_time(); + + // scan intensity slide's data + + bool wholeslide = env.singleROI; + double slide_I_max = (std::numeric_limits::lowest)(), + slide_I_min = (std::numeric_limits::max)(); + + std::unordered_set U; // unique ROI mask labels + std::unordered_map R; // ROI data + + int lvl = 0, // pyramid level + lyr = 0; // layer + + // Read the image/volume. The image loader is put in the open state in processDataset_XX_YY () + size_t nth = ilo.get_num_tiles_hor(), + ntv = ilo.get_num_tiles_vert(), + fw = ilo.get_tile_width(), + th = ilo.get_tile_height(), + tw = ilo.get_tile_width(), + tileSize = ilo.get_tile_size(), + fullwidth = ilo.get_full_width(), + fullheight = ilo.get_full_height(); + + // iterate abstract tiles (in a tiled slide /e.g. tiled tiff/ they correspond to physical tiles, in a nontiled slide /e.g. scanline tiff or strip tiff/ they correspond to ) + int cnt = 1; + for (unsigned int row = 0; row < nth; row++) + for (unsigned int col = 0; col < ntv; col++) + { + // Fetch the tile + if (!ilo.load_tile(row, col)) + { +#ifdef WITH_PYTHON_H + throw "Error fetching tile"; +#endif + std::cerr << "Error fetching tile\n"; + return false; + } + + // Get ahold of tile's pixel buffer + auto tidx = row * nth + col; + + // Iterate pixels + for (size_t i = 0; i < tileSize; i++) + { + // Mask + uint32_t msk = 1; // wholeslide by default + if (!wholeslide) + msk = ilo.get_cur_tile_seg_pixel(i); + + // Skip non-mask pixels + if (!msk) + continue; + + int y = row * th + i / tw, + x = col * tw + i % tw; + + // Skip tile buffer pixels beyond the image's bounds + if (x >= fullwidth || y >= fullheight) + continue; + + // dynamic range within- and off-ROI + double dxequiv_I = ilo.get_cur_tile_dpequiv_pixel(i); + slide_I_max = (std::max)(slide_I_max, dxequiv_I); + slide_I_min = (std::min)(slide_I_min, dxequiv_I); + + // Update pixel's ROI metrics + // - the following block mocks feed_pixel_2_metrics (x, y, dataI[i], msk, tidx) + if (U.find(msk) == U.end()) + { + // Remember this label + U.insert(msk); + + // Initialize the ROI label record + LR r(msk); + + // - mocking init_label_record_3 (roi, theSegFname, theIntFname, x, y, label, intensity, tile_index) + // Initialize basic counters + r.aux_area = 1; + r.aux_min = r.aux_max = 0; //we don't have uint-cast intensities at this moment + r.init_aabb(x, y); + + // - not storing file names (r.segFname = segFile, r.intFname = intFile) but will do so in the future + + // Attach + R[msk] = r; + } + else + { + // Update basic ROI info (info that doesn't require costly calculations) + LR& r = R[msk]; + + // - mocking update_label_record_2 (r, x, y, label, intensity, tile_index) + + // Per-ROI + r.aux_area++; + + // save + r.update_aabb(x, y); + } + } // scan tile + +#ifdef WITH_PYTHON_H + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); +#endif + + } // foreach tile + + //****** fix ROIs' AABBs with respect to anisotropy + + if (!aniso.customized()) + { + for (auto& pair : R) + { + LR& r = pair.second; + r.make_nonanisotropic_aabb(); + } + } + else + { + for (auto& pair : R) + { + LR& r = pair.second; + r.make_anisotropic_aabb(aniso.get_aniso_x(), aniso.get_aniso_y()); + } + } + + //****** Analysis + + // slide-wide (max ROI area) x (number of ROIs) + size_t maxArea = 0; + size_t max_w = 0, max_h = 0; + for (const auto& pair : R) + { + const LR& r = pair.second; + maxArea = maxArea > r.aux_area ? maxArea : r.aux_area; //std::max (maxArea, r.aux_area); + const AABB& bb = r.aabb; + auto w = bb.get_width(); + auto h = bb.get_height(); + max_w = max_w > w ? max_w : w; + max_h = max_h > h ? max_h : h; + } + + p.slide_w = fullwidth; + p.slide_h = fullheight; + + p.max_preroi_inten = slide_I_max; // in case fp_phys_pivoxels==true, max/min _preroi_inten + p.min_preroi_inten = slide_I_min; // needs adjusting (grey-binning) before using in wsi scenarios (assigning ROI's min and max) + + p.max_roi_area = maxArea; + p.n_rois = R.size(); + p.max_roi_w = max_w; + p.max_roi_h = max_h; + + return true; + } + + bool gatherRoisMetrics_2_slideprops_3D( + // in + RawImageLoader& ilo, + const AnisotropyOptions& aniso, + // out + SlideProps& p) + { + // low-level slide properties (intensity and mask, if available) + p.lolvl_slide_descr = ilo.get_slide_descr(); + p.fp_phys_pivoxels = ilo.get_fp_phys_pixvoxels(); + + // time series + p.inten_time = ilo.get_inten_time(); + p.mask_time = ilo.get_mask_time(); + + // scan intensity slide's data + + bool wholeslide = env.singleROI; + + double slide_I_max = (std::numeric_limits::lowest)(), + slide_I_min = (std::numeric_limits::max)(); + + std::unordered_set U; // unique ROI mask labels + std::unordered_map R; // ROI data + + // Read the volume. The image loader is in the open state by previously called processDataset_XX_YY () + size_t fullW = ilo.get_full_width(), + fullH = ilo.get_full_height(), + fullD = ilo.get_full_depth(), + sliceSize = fullW * fullH, + nVox = sliceSize * fullD; + + // in the 3D case tiling is a formality, so fetch the only tile in the file + if (!ilo.load_tile(0, 0)) + { +#ifdef WITH_PYTHON_H + throw "Error fetching tile"; +#endif + std::cerr << "Error fetching tile\n"; + return false; + } + + // iterate abstract tiles (in a tiled slide /e.g. tiled tiff/ they correspond to physical tiles, in a nontiled slide /e.g. scanline tiff or strip tiff/ they correspond to ) + int cnt = 1; + + // iterate voxels + for (size_t i = 0; i < nVox; i++) + { + // Mask + uint32_t msk = 1; // wholeslide by default + if (!wholeslide) + msk = ilo.get_cur_tile_seg_pixel(i); + + // Skip non-mask pixels + if (!msk) + continue; + + int z = i / sliceSize, + y = (i - z * sliceSize) / fullW, + x = (i - z * sliceSize) % fullW; + + // Skip tile buffer pixels beyond the image's bounds + if (x >= fullW || y >= fullH || z >= fullD) + continue; + + // dynamic range within- and off-ROI + double dxequiv_I = ilo.get_cur_tile_dpequiv_pixel(i); + slide_I_max = (std::max)(slide_I_max, dxequiv_I); + slide_I_min = (std::min)(slide_I_min, dxequiv_I); + + // Update pixel's ROI metrics + // - the following block mocks feed_pixel_2_metrics (x, y, dataI[i], msk, tidx) + if (U.find(msk) == U.end()) + { + // Remember this label + U.insert(msk); + + // Initialize the ROI label record + LR r(msk); + + // - mocking init_label_record_3 (newData, theSegFname, theIntFname, x, y, label, intensity, tile_index) + // Initialize basic counters + r.aux_area = 1; + r.aux_min = r.aux_max = 0; //we don't have uint-cast intensities at this moment + r.init_aabb_3D(x, y, z); + + // - not storing file names (r.segFname = segFile, r.intFname = intFile) but will do so in the future + + // Attach + R[msk] = r; + } + else + { + // Update basic ROI info (info that doesn't require costly calculations) + LR& r = R[msk]; + + // - mocking update_label_record_2 (r, x, y, label, intensity, tile_index) + + // Per-ROI + r.aux_area++; + + // save + r.update_aabb_3D(x, y, z); + } + +#ifdef WITH_PYTHON_H + // keyboard interrupt + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); +#endif + + } //- all voxels + + //****** fix ROIs' AABBs with respect to anisotropy + + if (!aniso.customized()) + { + for (auto& pair : R) + { + LR& r = pair.second; + r.make_nonanisotropic_aabb(); + } + } + else + { + for (auto& pair : R) + { + LR& r = pair.second; + r.make_anisotropic_aabb(aniso.get_aniso_x(), aniso.get_aniso_y(), aniso.get_aniso_z()); + } + } + + //****** Analysis + + // slide-wide (max ROI area) x (number of ROIs) + size_t maxArea = 0; + size_t max_w = 0, max_h = 0, max_d = 0; + for (const auto& pair : R) + { + const LR& r = pair.second; + maxArea = (std::max)(maxArea, (size_t)r.aux_area); + const AABB& bb = r.aabb; + auto w = bb.get_width(); + auto h = bb.get_height(); + auto d = bb.get_z_depth(); + max_w = (std::max)(max_w, (size_t)w); + max_h = (std::max)(max_h, (size_t)h); + max_d = (std::max)(max_d, (size_t)d); + } + + p.slide_w = fullW; + p.slide_h = fullH; + p.volume_d = fullD; + + p.max_preroi_inten = slide_I_max; // in case fp_phys_pivoxels==true, max/min _preroi_inten + p.min_preroi_inten = slide_I_min; // needs adjusting (grey-binning) before using in wsi scenarios (assigning ROI's min and max) + + p.max_roi_area = maxArea; + p.n_rois = R.size(); + p.max_roi_w = max_w; + p.max_roi_h = max_h; + p.max_roi_d = max_d; + + return true; + } + + std::pair split_alnum(const std::string& annot) + { + std::string A = annot; // a string that we can edit + std::string al; + for (auto c : A) + { + if (!std::isdigit(c)) + al += c; + else + { + A.erase(0, al.size()); + break; + } + } + + return { al, A }; + } + + bool scan_slide_props_montage (SlideProps & p, int dim, const AnisotropyOptions & aniso) + { + if (dim != 2) + return false; + + gatherRoisMetrics_2_slideprops_2D_montage (aniso, p); + + return true; + } + + // + // prerequisite: initialized fields fname_int and fname_seg + // + bool scan_slide_props (SlideProps & p, int dim, const AnisotropyOptions & aniso, bool need_annot) + { + RawImageLoader ilo; + if (! ilo.open(p.fname_int, p.fname_seg)) + { + std::cerr << "error opening an ImageLoader for " << p.fname_int << " | " << p.fname_seg << "\n"; + return false; + } + + bool ok = dim==2 ? gatherRoisMetrics_2_slideprops_2D(ilo, aniso, p) : gatherRoisMetrics_2_slideprops_3D(ilo, aniso, p); + if (!ok) + { + std::cerr << "error gathering ROI metrics to slide/volume props \n"; + return false; + } + + ilo.close(); + + // annotations + if (need_annot) + { + // throw away the directory part + fs::path pth(p.fname_seg); + auto purefn = pth.filename().string(); + + // LHS part till the 1st dot is the annotation info that we need to parse + std::vectortoks1; + Nyxus::parse_delimited_string (purefn, ".", toks1); + + // the result tokens is the annotation info that we need + p.annots.clear(); + Nyxus::parse_delimited_string (toks1[0], "_", p.annots); + + // prune blank annotations (usually caused by multiple separator e.g. '__' in 'slide_blah1_bla_2__something3.ome.tiff') + p.annots.erase ( + std::remove_if (p.annots.begin(), p.annots.end(), [](const std::string & s) { return s.empty(); }), + p.annots.end()); + } + + + return true; + } } \ No newline at end of file From 114212098846688a589f9a26e61888e59d9024bf Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Thu, 12 Mar 2026 15:24:50 -0400 Subject: [PATCH 02/10] tests: add TestSingleRoi regression tests for issue #327 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Five pytest cases covering the fixed single_roi=True path in featurize_files: - no mask files → succeeds, one row per image with label=1 - mask files supplied but single_roi=True → mask is ignored, results identical to no-mask case - single_roi=False (segmented) → multiple rows per image (unchanged behavior) - single_roi=False, no mask → raises IOError (unchanged behavior) - single_roi=True results consistent with featurize_directory when no seg dir given Co-Authored-By: Claude Sonnet 4.6 --- tests/python/test_nyxus.py | 70 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/tests/python/test_nyxus.py b/tests/python/test_nyxus.py index d5e578a1..16ffddf9 100644 --- a/tests/python/test_nyxus.py +++ b/tests/python/test_nyxus.py @@ -1095,3 +1095,73 @@ def test_3d_glszm_compatibility (self): assert np.isclose (f.at[0, "3GLSZM_ZV"], radiomics_gt["Case-1_original_glszm_ZoneVariance"], rtol=1.e-1, atol=1.e-2) + +class TestSingleRoi(): + """Regression tests for issue #327: featurize_files with single_roi=True.""" + + def setup_method(self): + path = str(pathlib.Path(__file__).parent.resolve()) + self.data_path = path + '/data/' + self.int_files = [ + self.data_path + 'int/p0_y1_r1_c0.ome.tif', + self.data_path + 'int/p0_y1_r1_c1.ome.tif', + ] + self.seg_files = [ + self.data_path + 'seg/p0_y1_r1_c0.ome.tif', + self.data_path + 'seg/p0_y1_r1_c1.ome.tif', + ] + self.nyx = nyxus.Nyxus(["*ALL_INTENSITY*"]) + + def test_single_roi_no_mask_succeeds(self): + """featurize_files(int_files, None, single_roi=True) must succeed and return one row per image.""" + df = self.nyx.featurize_files(self.int_files, None, single_roi=True) + assert df is not None + # one row per intensity file, all with label == 1 + assert len(df) == len(self.int_files) + assert (df["label"] == 1).all() + + def test_single_roi_mask_ignored(self): + """featurize_files with single_roi=True should produce the same result whether or not mask files are supplied.""" + df_no_mask = self.nyx.featurize_files(self.int_files, None, single_roi=True) + df_with_mask = self.nyx.featurize_files(self.int_files, self.seg_files, single_roi=True) + + df_no_mask.replace([np.inf, -np.inf, np.nan], 0, inplace=True) + df_with_mask.replace([np.inf, -np.inf, np.nan], 0, inplace=True) + + assert list(df_no_mask.columns) == list(df_with_mask.columns) + for col in df_no_mask.columns: + for v1, v2 in zip(df_no_mask[col].tolist(), df_with_mask[col].tolist()): + assert v1 == pytest.approx(v2, rel=1e-5, abs=1e-5) + + def test_single_roi_false_multi_row(self): + """featurize_files with single_roi=False (normal segmented mode) must return multiple rows per image.""" + df = self.nyx.featurize_files(self.int_files, self.seg_files, single_roi=False) + assert df is not None + # segmented mode produces more than one row per image + assert len(df) > len(self.int_files) + + def test_single_roi_false_no_mask_raises(self): + """featurize_files(int_files, None, single_roi=False) must raise IOError.""" + with pytest.raises(IOError): + self.nyx.featurize_files(self.int_files, None, single_roi=False) + + def test_single_roi_consistent_with_featurize_directory(self): + """single_roi=True via featurize_files should match featurize_directory when no seg dir is given.""" + df_dir = self.nyx.featurize_directory(self.data_path + 'int/', '') + df_files = self.nyx.featurize_files(self.int_files, None, single_roi=True) + + df_dir.replace([np.inf, -np.inf, np.nan], 0, inplace=True) + df_files.replace([np.inf, -np.inf, np.nan], 0, inplace=True) + + # Sort both by filename so row order is comparable + df_dir = df_dir.sort_values("intensity_image").reset_index(drop=True) + df_files = df_files.sort_values("intensity_image").reset_index(drop=True) + + shared_cols = [c for c in df_dir.columns if c in df_files.columns] + not_equal = [] + for col in shared_cols: + for v1, v2 in zip(df_dir[col].tolist(), df_files[col].tolist()): + if v1 != pytest.approx(v2, rel=1e-5, abs=1e-5): + not_equal.append(col) + break + assert len(not_equal) == 0, f"Columns differ: {not_equal}" From fb85e4b317ae5d00d453e81952f7d0ff4623dbd0 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Fri, 13 Mar 2026 14:37:51 -0400 Subject: [PATCH 03/10] Fix: remove theEnvironment.singleROI from dirs_and_files.cpp The polus version of read_2D_dataset() has no Environment parameter, so theEnvironment is undeclared there. The callers (featurize_directory_imp et al.) already set env.singleROI before calling read_2D_dataset, making this assignment redundant anyway. Co-Authored-By: Claude Sonnet 4.6 --- src/nyx/dirs_and_files.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/nyx/dirs_and_files.cpp b/src/nyx/dirs_and_files.cpp index a9145833..f5b5cce1 100644 --- a/src/nyx/dirs_and_files.cpp +++ b/src/nyx/dirs_and_files.cpp @@ -78,7 +78,6 @@ namespace Nyxus if (wholeslide) { - theEnvironment.singleROI = true; // populate with empty mask file names labelFiles.insert (labelFiles.begin(), intensFiles.size(), ""); } From 57e03970fd011bc3fbdfacca4f4a201686ce28b7 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Fri, 13 Mar 2026 17:33:34 -0400 Subject: [PATCH 04/10] Fix: remove unused wholeslide variable from scan_trivial_wholeslide_anisotropic MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This function has no Environment& parameter, so env.singleROI would fail to compile. The variable was also never used — the function is wholeslide- only by definition and feeds pixels unconditionally. Co-Authored-By: Claude Sonnet 4.6 --- src/nyx/phase2_2d.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/nyx/phase2_2d.cpp b/src/nyx/phase2_2d.cpp index 6cded7c0..8b46f996 100644 --- a/src/nyx/phase2_2d.cpp +++ b/src/nyx/phase2_2d.cpp @@ -404,12 +404,10 @@ namespace Nyxus ph_y = size_t(double(vy) / aniso_y), i = ph_y * tw + ph_x; - // read buffered physical pixel + // read buffered physical pixel const std::vector& dataI = ldr.get_int_tile_buffer(); - const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); - bool wholeslide = env.singleROI; - // Cache this pixel + // Cache this pixel feed_pixel_2_cache_LR (vc, vr, dataI[i], vroi); } } From d0663289a6a74582ce2cbab7af7df335d70f2b97 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Fri, 13 Mar 2026 17:37:08 -0400 Subject: [PATCH 05/10] Fix: revert slideprops.cpp wholeslide to p.fname_seg.empty() gatherRoisMetrics_2_slideprops_2D and _3D have no Environment parameter so env.singleROI doesn't compile there. p.fname_seg.empty() is the correct local signal: in the wholeslide path labelFiles are empty strings, so fname_seg is empty iff singleROI is true. Co-Authored-By: Claude Sonnet 4.6 --- src/nyx/slideprops.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nyx/slideprops.cpp b/src/nyx/slideprops.cpp index 41c97c59..d88cc349 100644 --- a/src/nyx/slideprops.cpp +++ b/src/nyx/slideprops.cpp @@ -98,7 +98,7 @@ namespace Nyxus // scan intensity slide's data - bool wholeslide = env.singleROI; + bool wholeslide = p.fname_seg.empty(); double slide_I_max = (std::numeric_limits::lowest)(), slide_I_min = (std::numeric_limits::max)(); @@ -269,7 +269,7 @@ namespace Nyxus // scan intensity slide's data - bool wholeslide = env.singleROI; + bool wholeslide = p.fname_seg.empty(); double slide_I_max = (std::numeric_limits::lowest)(), slide_I_min = (std::numeric_limits::max)(); From 7e9b897d6e432a8643fcd197c743633f49c00579 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Fri, 13 Mar 2026 18:09:30 -0400 Subject: [PATCH 06/10] Fix: remove unused wholeslide variable from scan_trivial_wholevolume_anisotropic__OLD MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Same issue as scan_trivial_wholeslide_anisotropic in phase2_2d.cpp: no Environment parameter so env.singleROI doesn't compile, and the variable was never used anyway — this is a wholevolume-only function. Co-Authored-By: Claude Sonnet 4.6 --- src/nyx/phase2_3d.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/nyx/phase2_3d.cpp b/src/nyx/phase2_3d.cpp index d4122070..f88f9650 100644 --- a/src/nyx/phase2_3d.cpp +++ b/src/nyx/phase2_3d.cpp @@ -614,12 +614,10 @@ namespace Nyxus ph_y = size_t(double(vy) / aniso_y), i = ph_y * tw + ph_x; - // read buffered physical pixel + // read buffered physical pixel const std::vector& dataI = ldr.get_int_tile_buffer(); - const std::shared_ptr>& spL = ldr.get_seg_tile_sptr(); - bool wholeslide = env.singleROI; - // Cache this pixel + // Cache this pixel feed_pixel_2_cache_LR(vc, vr, dataI[i], vroi); } } From 05363ccc05d2c3121e032a282c8bbf553bbcbad3 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Fri, 13 Mar 2026 18:26:05 -0400 Subject: [PATCH 07/10] tests: fix TestSingleRoi for polus codebase - Use ROI_label instead of label (polus column name) - featurize_directory wholeslide: pass int_dir as label_dir (same-dir detection) instead of '' which polus rejects with OSError Co-Authored-By: Claude Sonnet 4.6 --- tests/python/test_nyxus.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/python/test_nyxus.py b/tests/python/test_nyxus.py index 16ffddf9..97381f8b 100644 --- a/tests/python/test_nyxus.py +++ b/tests/python/test_nyxus.py @@ -1118,7 +1118,7 @@ def test_single_roi_no_mask_succeeds(self): assert df is not None # one row per intensity file, all with label == 1 assert len(df) == len(self.int_files) - assert (df["label"] == 1).all() + assert (df["ROI_label"] == 1).all() def test_single_roi_mask_ignored(self): """featurize_files with single_roi=True should produce the same result whether or not mask files are supplied.""" @@ -1147,7 +1147,8 @@ def test_single_roi_false_no_mask_raises(self): def test_single_roi_consistent_with_featurize_directory(self): """single_roi=True via featurize_files should match featurize_directory when no seg dir is given.""" - df_dir = self.nyx.featurize_directory(self.data_path + 'int/', '') + # Pass the same dir as label dir — polus detects same-dir as wholeslide + df_dir = self.nyx.featurize_directory(self.data_path + 'int/', self.data_path + 'int/') df_files = self.nyx.featurize_files(self.int_files, None, single_roi=True) df_dir.replace([np.inf, -np.inf, np.nan], 0, inplace=True) From 8d0f1d56d71baaff49b72ad817e5d8aaa382f8d7 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Wed, 18 Mar 2026 11:31:57 -0400 Subject: [PATCH 08/10] tests: exclude path columns from featurize_directory/featurize_files comparison intensity_image legitimately differs in format (full path vs basename) between the two code paths, so skip string columns in numeric comparison. Co-Authored-By: Claude Sonnet 4.6 --- tests/python/test_nyxus.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/python/test_nyxus.py b/tests/python/test_nyxus.py index 97381f8b..bedf95ae 100644 --- a/tests/python/test_nyxus.py +++ b/tests/python/test_nyxus.py @@ -1158,7 +1158,9 @@ def test_single_roi_consistent_with_featurize_directory(self): df_dir = df_dir.sort_values("intensity_image").reset_index(drop=True) df_files = df_files.sort_values("intensity_image").reset_index(drop=True) - shared_cols = [c for c in df_dir.columns if c in df_files.columns] + # Exclude path/string columns — they legitimately differ in format (full path vs basename) + string_cols = {"intensity_image", "mask_image"} + shared_cols = [c for c in df_dir.columns if c in df_files.columns and c not in string_cols] not_equal = [] for col in shared_cols: for v1, v2 in zip(df_dir[col].tolist(), df_files[col].tolist()): From 08d2692b376f5a0ac920e8c0b1a57cca2e4921cd Mon Sep 17 00:00:00 2001 From: Sameeul B Samee Date: Wed, 18 Mar 2026 13:38:40 -0400 Subject: [PATCH 09/10] fix: restrict libarrow and libparquet so that C++20 stuff (span etc.) does not get it. --- ci-utils/envs/conda_cpp.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ci-utils/envs/conda_cpp.txt b/ci-utils/envs/conda_cpp.txt index 8a963274..2a5c719b 100644 --- a/ci-utils/envs/conda_cpp.txt +++ b/ci-utils/envs/conda_cpp.txt @@ -9,5 +9,5 @@ xsimd >=13,<14 cmake dcmtk >=3.6.9 fmjpeg2koj >=1.0.3 -libarrow -libparquet +libarrow <=23.0.0 +libparquet <=23.0.0 \ No newline at end of file From e5425814abd1640d9f5a704d63b9f1f2fd3a4dbd Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Wed, 18 Mar 2026 16:11:07 -0400 Subject: [PATCH 10/10] fix: python 3.12 fix for async() --- src/nyx/phase1.cpp | 5 ++--- src/nyx/python/nyxus/nyxus.py | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/nyx/phase1.cpp b/src/nyx/phase1.cpp index 31bc8c71..1d521675 100644 --- a/src/nyx/phase1.cpp +++ b/src/nyx/phase1.cpp @@ -390,10 +390,9 @@ namespace Nyxus // Update pixel's ROI metrics auto inten = rI (pair_index, row, col); feed_pixel_2_metrics (env.uniqueLabels, env.roiData, col, row, inten, label, pair_index); // Updates 'uniqueLabels' and 'roiData' - - if (PyErr_CheckSignals() != 0) - throw pybind11::error_already_set(); } + if (PyErr_CheckSignals() != 0) + throw pybind11::error_already_set(); return true; } diff --git a/src/nyx/python/nyxus/nyxus.py b/src/nyx/python/nyxus/nyxus.py index bc3a8224..646d458e 100644 --- a/src/nyx/python/nyxus/nyxus.py +++ b/src/nyx/python/nyxus/nyxus.py @@ -446,8 +446,8 @@ def featurize( min_raw_I = np.min (intensity_images) if (min_raw_I < 0): I -= min_raw_I - if (not isinstance(I.flat[0], np.uint32)): - I = I.astype (np.uint32) + if I.dtype != np.uint32: + I = I.astype(np.uint32) # cast mask data to unsigned integer, too M = label_images.astype (np.uint32)