From f869a6c6aede1ba702bfbe8184e1e946b6b97b15 Mon Sep 17 00:00:00 2001 From: Joey Kleingers Date: Sat, 28 Feb 2026 14:56:09 -0500 Subject: [PATCH 1/5] ENH: Split BadDataNeighborOrientationCheck into in-core and OOC algorithms Split BadDataNeighborOrientationCheck into two dispatched algorithms using DispatchAlgorithm for optimal performance in both configurations: - Worklist (in-core): Uses std::deque worklist for Phase 2 to process only eligible voxels with fast random access. ~5x speedup vs original. - Scanline (OOC): Uses chunk-sequential multi-pass scans for Phase 2 to avoid random access chunk thrashing. Includes chunk-skip optimization that checks in-memory neighborCount before loading chunks, skipping those with no eligible voxels. ~1.8x speedup vs original. Both algorithms share Phase 1 (chunk-sequential neighbor counting) and use only a single neighborCount vector (4 bytes/voxel) with no additional large allocations. Updated tests with GENERATE + ForceOocAlgorithmGuard to exercise both algorithm paths in in-core builds. Added 200x200x200 benchmark test. Signed-off-by: Joey Kleingers --- .../OrientationAnalysis/CMakeLists.txt | 2 + .../BadDataNeighborOrientationCheck.cpp | 181 +------------ ...adDataNeighborOrientationCheckScanline.cpp | 232 ++++++++++++++++ ...adDataNeighborOrientationCheckScanline.hpp | 35 +++ ...adDataNeighborOrientationCheckWorklist.cpp | 180 +++++++++++++ ...adDataNeighborOrientationCheckWorklist.hpp | 35 +++ .../BadDataNeighborOrientationCheckTest.cpp | 252 ++++++++++++++++++ 7 files changed, 745 insertions(+), 172 deletions(-) create mode 100644 src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.cpp create mode 100644 src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.hpp create mode 100644 src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.cpp create mode 100644 src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.hpp diff --git a/src/Plugins/OrientationAnalysis/CMakeLists.txt b/src/Plugins/OrientationAnalysis/CMakeLists.txt index 74ac32e5b4..82dca74e8a 100644 --- a/src/Plugins/OrientationAnalysis/CMakeLists.txt +++ b/src/Plugins/OrientationAnalysis/CMakeLists.txt @@ -168,6 +168,8 @@ set(filter_algorithms AlignSectionsMisorientation AlignSectionsMutualInformation BadDataNeighborOrientationCheck + BadDataNeighborOrientationCheckScanline + BadDataNeighborOrientationCheckWorklist CAxisSegmentFeatures ComputeAvgCAxes ComputeAvgOrientations diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheck.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheck.cpp index 2937ac640c..92bf066b3c 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheck.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheck.cpp @@ -1,13 +1,10 @@ #include "BadDataNeighborOrientationCheck.hpp" -#include "simplnx/Common/Numbers.hpp" -#include "simplnx/DataStructure/DataArray.hpp" -#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" -#include "simplnx/Utilities/MaskCompareUtilities.hpp" -#include "simplnx/Utilities/MessageHelper.hpp" -#include "simplnx/Utilities/NeighborUtilities.hpp" +#include "BadDataNeighborOrientationCheckScanline.hpp" +#include "BadDataNeighborOrientationCheckWorklist.hpp" -#include +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" using namespace nx::core; @@ -33,170 +30,10 @@ const std::atomic_bool& BadDataNeighborOrientationCheck::getCancel() // ----------------------------------------------------------------------------- Result<> BadDataNeighborOrientationCheck::operator()() { - const float misorientationTolerance = m_InputValues->MisorientationTolerance * numbers::pi_v / 180.0f; - - const auto* imageGeomPtr = m_DataStructure.getDataAs(m_InputValues->ImageGeomPath); - SizeVec3 udims = imageGeomPtr->getDimensions(); - const auto& cellPhases = m_DataStructure.getDataRefAs(m_InputValues->CellPhasesArrayPath); - const auto& quats = m_DataStructure.getDataRefAs(m_InputValues->QuatsArrayPath); - const auto& crystalStructures = m_DataStructure.getDataRefAs(m_InputValues->CrystalStructuresArrayPath); - const usize totalPoints = quats.getNumberOfTuples(); - - std::unique_ptr maskCompare; - try - { - maskCompare = MaskCompareUtilities::InstantiateMaskCompare(m_DataStructure, m_InputValues->MaskArrayPath); - } catch(const std::out_of_range& exception) - { - // This really should NOT be happening as the path was verified during preflight, BUT we may be calling this from - // somewhere else that is NOT going through the normal nx::core::IFilter API of Preflight and Execute - return MakeErrorResult(-54900, fmt::format("Mask Array DataPath does not exist or is not of the correct type (Bool | UInt8) {}", m_InputValues->MaskArrayPath.toString())); - } - - std::array dims = { - static_cast(udims[0]), - static_cast(udims[1]), - static_cast(udims[2]), - }; - - std::array neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims); - std::array faceNeighborInternalIdx = initializeFaceNeighborInternalIdx(); - - std::vector orientationOps = ebsdlib::LaueOps::GetAllOrientationOps(); - - std::vector neighborCount(totalPoints, 0); - - MessageHelper messageHelper(m_MessageHandler); - ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger(); - // Loop over every point finding the number of neighbors that fall within the - // user defined angle tolerance. - for(int64 voxelIndex = 0; voxelIndex < totalPoints; voxelIndex++) - { - throttledMessenger.sendThrottledMessage([&] { return fmt::format("Processing Data {:.2f}% completed", CalculatePercentComplete(voxelIndex, totalPoints)); }); - // If the mask was set to false, then we check this voxel - if(!maskCompare->isTrue(voxelIndex)) - { - // We precalculate the positive voxel quaternion and laue class here to prevent reading and recalculating it for each face below - ebsdlib::QuatD quat1(quats[voxelIndex * 4], quats[voxelIndex * 4 + 1], quats[voxelIndex * 4 + 2], quats[voxelIndex * 4 + 3]); - quat1.positiveOrientation(); - const uint32 laueClass1 = crystalStructures[cellPhases[voxelIndex]]; - - int64 xIdx = voxelIndex % dims[0]; - int64 yIdx = (voxelIndex / dims[0]) % dims[1]; - int64 zIdx = voxelIndex / (dims[0] * dims[1]); - - // Loop over the 6 face neighbors of the voxel - std::array isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims); - for(const auto& faceIndex : faceNeighborInternalIdx) - { - if(!isValidFaceNeighbor[faceIndex]) - { - continue; - } - const int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; - - // Now compare the mask of the neighbor. If the mask is TRUE, i.e., that voxel - // did not fail the threshold filter that most likely produced the mask array, - // then we can look at that voxel. - if(maskCompare->isTrue(neighborPoint)) - { - // Both Cell Phases MUST be the same and be a valid Phase - if(cellPhases[voxelIndex] == cellPhases[neighborPoint] && cellPhases[voxelIndex] > 0) - { - ebsdlib::QuatD quat2(quats[neighborPoint * 4], quats[neighborPoint * 4 + 1], quats[neighborPoint * 4 + 2], quats[neighborPoint * 4 + 3]); - quat2.positiveOrientation(); - // Compute the Axis_Angle misorientation between those 2 quaternions - ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClass1]->calculateMisorientation(quat1, quat2); - // if the angle is less than our tolerance, then we increment the neighbor count - // for this voxel - if(axisAngle[3] < misorientationTolerance) - { - neighborCount[voxelIndex]++; - } - } - } - } - } - } - - constexpr int32 startLevel = 6; - int32 currentLevel = startLevel; - int32 counter = 0; - - // Now we loop over all the points again, but this time we do it as many times - // as the user has requested to iteratively flip voxels - while(currentLevel >= m_InputValues->NumberOfNeighbors) - { - counter = 1; - int32 loopNumber = 0; - while(counter > 0) - { - counter = 0; // Set this while control variable to zero - for(usize voxelIndex = 0; voxelIndex < totalPoints; voxelIndex++) - { - throttledMessenger.sendThrottledMessage([&] { - return fmt::format("Level '{}' of '{}' || Processing Data ('{}') {:.2f}% completed", (startLevel - currentLevel) + 1, startLevel - m_InputValues->NumberOfNeighbors, loopNumber, - CalculatePercentComplete(voxelIndex, totalPoints)); - }); - - // We are comparing the number-of-neighbors of the current voxel, and if it - // is > the current level and the mask is FALSE, then we drop into this - // conditional. The first thing that happens in the conditional is that - // the current voxel's mask value is set to TRUE. - if(neighborCount[voxelIndex] >= currentLevel && !maskCompare->isTrue(voxelIndex)) - { - maskCompare->setValue(voxelIndex, true); // the current voxel's mask value is set to TRUE. - counter++; // Increment the `counter` to force the loop to iterate again - - // We precalculate the positive voxel quaternion and laue class here to prevent reading and recalculating it for each face below - ebsdlib::QuatD quat1(quats[voxelIndex * 4], quats[voxelIndex * 4 + 1], quats[voxelIndex * 4 + 2], quats[voxelIndex * 4 + 3]); - quat1.positiveOrientation(); - const uint32 laueClass1 = crystalStructures[cellPhases[voxelIndex]]; - - // This whole section below is to now look at the neighbor voxels of the - // current voxel that just got flipped to true. This is needed because - // if any of those neighbor's mask was `false`, then its neighbor count - // is now not correct and will be off-by-one. So we run _almost_ the same - // loop code as above but checking the specific neighbors of the current - // voxel. This part should be termed the "Update Neighbor's Neighbor Count" - int64 xIdx = voxelIndex % dims[0]; - int64 yIdx = (voxelIndex / dims[0]) % dims[1]; - int64 zIdx = voxelIndex / (dims[0] * dims[1]); - - // Loop over the 6 face neighbors of the voxel - std::array isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims); - for(const auto& faceIndex : faceNeighborInternalIdx) - { - if(!isValidFaceNeighbor[faceIndex]) - { - continue; - } - - int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; - - // If the neighbor voxel's mask is false, then .... - if(!maskCompare->isTrue(neighborPoint)) - { - // Make sure both cells phase values are identical and valid - if(cellPhases[voxelIndex] == cellPhases[neighborPoint] && cellPhases[voxelIndex] > 0) - { - ebsdlib::QuatD quat2(quats[neighborPoint * 4], quats[neighborPoint * 4 + 1], quats[neighborPoint * 4 + 2], quats[neighborPoint * 4 + 3]); - quat2.positiveOrientation(); - // Quaternion Math is not commutative so do not reorder - ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClass1]->calculateMisorientation(quat1, quat2); - if(axisAngle[3] < misorientationTolerance) - { - neighborCount[neighborPoint]++; - } - } - } - } - } - } - ++loopNumber; - } - currentLevel = currentLevel - 1; - } + auto* quatsArray = m_DataStructure.getDataAs(m_InputValues->QuatsArrayPath); + auto* maskArray = m_DataStructure.getDataAs(m_InputValues->MaskArrayPath); + auto* phasesArray = m_DataStructure.getDataAs(m_InputValues->CellPhasesArrayPath); - return {}; + return DispatchAlgorithm({quatsArray, maskArray, phasesArray}, m_DataStructure, m_MessageHandler, m_ShouldCancel, + m_InputValues); } diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.cpp new file mode 100644 index 0000000000..f198ff60b0 --- /dev/null +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.cpp @@ -0,0 +1,232 @@ +#include "BadDataNeighborOrientationCheckScanline.hpp" + +#include "BadDataNeighborOrientationCheck.hpp" + +#include "simplnx/Common/Numbers.hpp" +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/Utilities/MaskCompareUtilities.hpp" +#include "simplnx/Utilities/MessageHelper.hpp" +#include "simplnx/Utilities/NeighborUtilities.hpp" + +#include + +using namespace nx::core; + +// ----------------------------------------------------------------------------- +BadDataNeighborOrientationCheckScanline::BadDataNeighborOrientationCheckScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const BadDataNeighborOrientationCheckInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +BadDataNeighborOrientationCheckScanline::~BadDataNeighborOrientationCheckScanline() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> BadDataNeighborOrientationCheckScanline::operator()() +{ + const float misorientationTolerance = m_InputValues->MisorientationTolerance * numbers::pi_v / 180.0f; + + const auto* imageGeomPtr = m_DataStructure.getDataAs(m_InputValues->ImageGeomPath); + SizeVec3 udims = imageGeomPtr->getDimensions(); + const auto& cellPhases = m_DataStructure.getDataRefAs(m_InputValues->CellPhasesArrayPath); + auto& quats = m_DataStructure.getDataRefAs(m_InputValues->QuatsArrayPath); + const auto& crystalStructures = m_DataStructure.getDataRefAs(m_InputValues->CrystalStructuresArrayPath); + const usize totalPoints = quats.getNumberOfTuples(); + + std::unique_ptr maskCompare; + try + { + maskCompare = MaskCompareUtilities::InstantiateMaskCompare(m_DataStructure, m_InputValues->MaskArrayPath); + } catch(const std::out_of_range& exception) + { + return MakeErrorResult(-54900, fmt::format("Mask Array DataPath does not exist or is not of the correct type (Bool | UInt8) {}", m_InputValues->MaskArrayPath.toString())); + } + + std::array dims = { + static_cast(udims[0]), + static_cast(udims[1]), + static_cast(udims[2]), + }; + + const int64 xyStride = dims[0] * dims[1]; + + std::array neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims); + std::array faceNeighborInternalIdx = initializeFaceNeighborInternalIdx(); + + std::vector orientationOps = ebsdlib::LaueOps::GetAllOrientationOps(); + + std::vector neighborCount(totalPoints, 0); + + MessageHelper messageHelper(m_MessageHandler); + ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger(); + + // ===== Phase 1: Count matching good neighbors for each bad voxel ===== + // Chunk-sequential iteration for OOC efficiency (no-op for in-core) + auto& quatsStore = quats.getDataStoreRef(); + const uint64 numChunks = quatsStore.getNumberOfChunks(); + + for(uint64 chunkIdx = 0; chunkIdx < numChunks; chunkIdx++) + { + quatsStore.loadChunk(chunkIdx); + + const auto chunkLowerBounds = quatsStore.getChunkLowerBounds(chunkIdx); + const auto chunkUpperBounds = quatsStore.getChunkUpperBounds(chunkIdx); + + for(usize zIdx = chunkLowerBounds[0]; zIdx <= chunkUpperBounds[0]; zIdx++) + { + const int64 kStride = static_cast(zIdx) * xyStride; + for(usize yIdx = chunkLowerBounds[1]; yIdx <= chunkUpperBounds[1]; yIdx++) + { + const int64 jStride = static_cast(yIdx) * dims[0]; + for(usize xIdx = chunkLowerBounds[2]; xIdx <= chunkUpperBounds[2]; xIdx++) + { + const int64 voxelIndex = kStride + jStride + static_cast(xIdx); + + throttledMessenger.sendThrottledMessage([&] { return fmt::format("Processing Data {:.2f}% completed", CalculatePercentComplete(voxelIndex, totalPoints)); }); + + if(!maskCompare->isTrue(voxelIndex)) + { + ebsdlib::QuatD quat1(quats[voxelIndex * 4], quats[voxelIndex * 4 + 1], quats[voxelIndex * 4 + 2], quats[voxelIndex * 4 + 3]); + quat1.positiveOrientation(); + const uint32 laueClass1 = crystalStructures[cellPhases[voxelIndex]]; + + std::array isValidFaceNeighbor = computeValidFaceNeighbors(static_cast(xIdx), static_cast(yIdx), static_cast(zIdx), dims); + for(const auto& faceIndex : faceNeighborInternalIdx) + { + if(!isValidFaceNeighbor[faceIndex]) + { + continue; + } + const int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; + + if(maskCompare->isTrue(neighborPoint)) + { + if(cellPhases[voxelIndex] == cellPhases[neighborPoint] && cellPhases[voxelIndex] > 0) + { + ebsdlib::QuatD quat2(quats[neighborPoint * 4], quats[neighborPoint * 4 + 1], quats[neighborPoint * 4 + 2], quats[neighborPoint * 4 + 3]); + quat2.positiveOrientation(); + ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClass1]->calculateMisorientation(quat1, quat2); + if(axisAngle[3] < misorientationTolerance) + { + neighborCount[voxelIndex]++; + } + } + } + } + } + } + } + } + } + + // ===== Phase 2: Iteratively flip bad voxels using chunk-sequential multi-pass scans ===== + // For each level, repeatedly scan the volume sequentially. On each pass, flip eligible + // voxels and update their neighbors' counts. Repeat until no flips occur. + // This avoids the random access pattern of worklist-based processing. + // + // Chunk-skip optimization: Before loading a chunk from disk, scan the in-memory + // neighborCount vector for that chunk's voxel range. If no voxel has + // neighborCount >= currentLevel, skip the chunk entirely (no disk I/O). + // Flipped voxels get neighborCount = -1 to prevent false positives. + constexpr int32 startLevel = 6; + const int32 totalLevels = startLevel - m_InputValues->NumberOfNeighbors + 1; + + for(int32 currentLevel = startLevel; currentLevel >= m_InputValues->NumberOfNeighbors; currentLevel--) + { + bool changed = true; + int32 passCount = 0; + + while(changed) + { + changed = false; + passCount++; + + throttledMessenger.sendThrottledMessage([&] { return fmt::format("Level '{}' of '{}' || Pass {}", (startLevel - currentLevel) + 1, totalLevels, passCount); }); + + for(uint64 chunkIdx = 0; chunkIdx < numChunks; chunkIdx++) + { + const auto chunkLowerBounds = quatsStore.getChunkLowerBounds(chunkIdx); + const auto chunkUpperBounds = quatsStore.getChunkUpperBounds(chunkIdx); + + // Compute the flat voxel index range for this chunk and check if any + // voxel in the chunk could be eligible (neighborCount >= currentLevel). + // This check uses only the in-memory neighborCount vector — no disk I/O. + const usize chunkStartIdx = chunkLowerBounds[0] * static_cast(xyStride) + chunkLowerBounds[1] * static_cast(dims[0]) + chunkLowerBounds[2]; + const usize chunkEndIdx = chunkUpperBounds[0] * static_cast(xyStride) + chunkUpperBounds[1] * static_cast(dims[0]) + chunkUpperBounds[2]; + + bool hasEligible = false; + for(usize i = chunkStartIdx; i <= chunkEndIdx; i++) + { + if(neighborCount[i] >= currentLevel) + { + hasEligible = true; + break; + } + } + + if(!hasEligible) + { + continue; + } + + quatsStore.loadChunk(chunkIdx); + + for(usize zIdx = chunkLowerBounds[0]; zIdx <= chunkUpperBounds[0]; zIdx++) + { + const int64 kStride = static_cast(zIdx) * xyStride; + for(usize yIdx = chunkLowerBounds[1]; yIdx <= chunkUpperBounds[1]; yIdx++) + { + const int64 jStride = static_cast(yIdx) * dims[0]; + for(usize xIdx = chunkLowerBounds[2]; xIdx <= chunkUpperBounds[2]; xIdx++) + { + const int64 voxelIndex = kStride + jStride + static_cast(xIdx); + + if(!maskCompare->isTrue(voxelIndex) && neighborCount[voxelIndex] >= currentLevel) + { + maskCompare->setValue(voxelIndex, true); + neighborCount[voxelIndex] = -1; // Mark as flipped to prevent false positives in chunk-skip check + changed = true; + + ebsdlib::QuatD quat1(quats[voxelIndex * 4], quats[voxelIndex * 4 + 1], quats[voxelIndex * 4 + 2], quats[voxelIndex * 4 + 3]); + quat1.positiveOrientation(); + const uint32 laueClass1 = crystalStructures[cellPhases[voxelIndex]]; + + std::array isValidFaceNeighbor = computeValidFaceNeighbors(static_cast(xIdx), static_cast(yIdx), static_cast(zIdx), dims); + for(const auto& faceIndex : faceNeighborInternalIdx) + { + if(!isValidFaceNeighbor[faceIndex]) + { + continue; + } + + const int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; + + if(!maskCompare->isTrue(neighborPoint)) + { + if(cellPhases[voxelIndex] == cellPhases[neighborPoint] && cellPhases[voxelIndex] > 0) + { + ebsdlib::QuatD quat2(quats[neighborPoint * 4], quats[neighborPoint * 4 + 1], quats[neighborPoint * 4 + 2], quats[neighborPoint * 4 + 3]); + quat2.positiveOrientation(); + ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClass1]->calculateMisorientation(quat1, quat2); + if(axisAngle[3] < misorientationTolerance) + { + neighborCount[neighborPoint]++; + } + } + } + } + } + } + } + } + } + } + } + + return {}; +} diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.hpp new file mode 100644 index 0000000000..1057ab1ab4 --- /dev/null +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.hpp @@ -0,0 +1,35 @@ +#pragma once + +#include "OrientationAnalysis/OrientationAnalysis_export.hpp" + +#include "simplnx/DataStructure/DataPath.hpp" +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ + +struct BadDataNeighborOrientationCheckInputValues; + +class ORIENTATIONANALYSIS_EXPORT BadDataNeighborOrientationCheckScanline +{ +public: + BadDataNeighborOrientationCheckScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const BadDataNeighborOrientationCheckInputValues* inputValues); + ~BadDataNeighborOrientationCheckScanline() noexcept; + + BadDataNeighborOrientationCheckScanline(const BadDataNeighborOrientationCheckScanline&) = delete; + BadDataNeighborOrientationCheckScanline(BadDataNeighborOrientationCheckScanline&&) noexcept = delete; + BadDataNeighborOrientationCheckScanline& operator=(const BadDataNeighborOrientationCheckScanline&) = delete; + BadDataNeighborOrientationCheckScanline& operator=(BadDataNeighborOrientationCheckScanline&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const BadDataNeighborOrientationCheckInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.cpp new file mode 100644 index 0000000000..f3b8702b7f --- /dev/null +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.cpp @@ -0,0 +1,180 @@ +#include "BadDataNeighborOrientationCheckWorklist.hpp" + +#include "BadDataNeighborOrientationCheck.hpp" + +#include "simplnx/Common/Numbers.hpp" +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/Utilities/MaskCompareUtilities.hpp" +#include "simplnx/Utilities/MessageHelper.hpp" +#include "simplnx/Utilities/NeighborUtilities.hpp" + +#include + +#include + +using namespace nx::core; + +// ----------------------------------------------------------------------------- +BadDataNeighborOrientationCheckWorklist::BadDataNeighborOrientationCheckWorklist(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const BadDataNeighborOrientationCheckInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +BadDataNeighborOrientationCheckWorklist::~BadDataNeighborOrientationCheckWorklist() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> BadDataNeighborOrientationCheckWorklist::operator()() +{ + const float misorientationTolerance = m_InputValues->MisorientationTolerance * numbers::pi_v / 180.0f; + + const auto* imageGeomPtr = m_DataStructure.getDataAs(m_InputValues->ImageGeomPath); + SizeVec3 udims = imageGeomPtr->getDimensions(); + const auto& cellPhases = m_DataStructure.getDataRefAs(m_InputValues->CellPhasesArrayPath); + auto& quats = m_DataStructure.getDataRefAs(m_InputValues->QuatsArrayPath); + const auto& crystalStructures = m_DataStructure.getDataRefAs(m_InputValues->CrystalStructuresArrayPath); + const usize totalPoints = quats.getNumberOfTuples(); + + std::unique_ptr maskCompare; + try + { + maskCompare = MaskCompareUtilities::InstantiateMaskCompare(m_DataStructure, m_InputValues->MaskArrayPath); + } catch(const std::out_of_range& exception) + { + return MakeErrorResult(-54900, fmt::format("Mask Array DataPath does not exist or is not of the correct type (Bool | UInt8) {}", m_InputValues->MaskArrayPath.toString())); + } + + std::array dims = { + static_cast(udims[0]), + static_cast(udims[1]), + static_cast(udims[2]), + }; + + const int64 xyStride = dims[0] * dims[1]; + + std::array neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims); + std::array faceNeighborInternalIdx = initializeFaceNeighborInternalIdx(); + + std::vector orientationOps = ebsdlib::LaueOps::GetAllOrientationOps(); + + std::vector neighborCount(totalPoints, 0); + + MessageHelper messageHelper(m_MessageHandler); + ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger(); + + // ===== Phase 1: Count matching good neighbors for each bad voxel ===== + for(usize voxelIndex = 0; voxelIndex < totalPoints; voxelIndex++) + { + throttledMessenger.sendThrottledMessage([&] { return fmt::format("Processing Data {:.2f}% completed", CalculatePercentComplete(voxelIndex, totalPoints)); }); + + if(!maskCompare->isTrue(voxelIndex)) + { + ebsdlib::QuatD quat1(quats[voxelIndex * 4], quats[voxelIndex * 4 + 1], quats[voxelIndex * 4 + 2], quats[voxelIndex * 4 + 3]); + quat1.positiveOrientation(); + const uint32 laueClass1 = crystalStructures[cellPhases[voxelIndex]]; + + const int64 xIdx = static_cast(voxelIndex) % dims[0]; + const int64 yIdx = (static_cast(voxelIndex) / dims[0]) % dims[1]; + const int64 zIdx = static_cast(voxelIndex) / xyStride; + + std::array isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims); + for(const auto& faceIndex : faceNeighborInternalIdx) + { + if(!isValidFaceNeighbor[faceIndex]) + { + continue; + } + const int64 neighborPoint = static_cast(voxelIndex) + neighborVoxelIndexOffsets[faceIndex]; + + if(maskCompare->isTrue(neighborPoint)) + { + if(cellPhases[voxelIndex] == cellPhases[neighborPoint] && cellPhases[voxelIndex] > 0) + { + ebsdlib::QuatD quat2(quats[neighborPoint * 4], quats[neighborPoint * 4 + 1], quats[neighborPoint * 4 + 2], quats[neighborPoint * 4 + 3]); + quat2.positiveOrientation(); + ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClass1]->calculateMisorientation(quat1, quat2); + if(axisAngle[3] < misorientationTolerance) + { + neighborCount[voxelIndex]++; + } + } + } + } + } + } + + // ===== Phase 2: Iteratively flip bad voxels using worklist ===== + constexpr int32 startLevel = 6; + const int32 totalLevels = startLevel - m_InputValues->NumberOfNeighbors + 1; + + for(int32 currentLevel = startLevel; currentLevel >= m_InputValues->NumberOfNeighbors; currentLevel--) + { + throttledMessenger.sendThrottledMessage([&] { return fmt::format("Level '{}' of '{}' || Building worklist", (startLevel - currentLevel) + 1, totalLevels); }); + + std::deque worklist; + for(usize voxelIndex = 0; voxelIndex < totalPoints; voxelIndex++) + { + if(neighborCount[voxelIndex] >= currentLevel && !maskCompare->isTrue(voxelIndex)) + { + worklist.push_back(voxelIndex); + } + } + + while(!worklist.empty()) + { + const usize voxelIndex = worklist.front(); + worklist.pop_front(); + + if(maskCompare->isTrue(voxelIndex) || neighborCount[voxelIndex] < currentLevel) + { + continue; + } + + maskCompare->setValue(voxelIndex, true); + + ebsdlib::QuatD quat1(quats[voxelIndex * 4], quats[voxelIndex * 4 + 1], quats[voxelIndex * 4 + 2], quats[voxelIndex * 4 + 3]); + quat1.positiveOrientation(); + const uint32 laueClass1 = crystalStructures[cellPhases[voxelIndex]]; + + const int64 xIdx = static_cast(voxelIndex) % dims[0]; + const int64 yIdx = (static_cast(voxelIndex) / dims[0]) % dims[1]; + const int64 zIdx = static_cast(voxelIndex) / xyStride; + + std::array isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims); + for(const auto& faceIndex : faceNeighborInternalIdx) + { + if(!isValidFaceNeighbor[faceIndex]) + { + continue; + } + + const int64 neighborPoint = static_cast(voxelIndex) + neighborVoxelIndexOffsets[faceIndex]; + + if(!maskCompare->isTrue(neighborPoint)) + { + if(cellPhases[voxelIndex] == cellPhases[neighborPoint] && cellPhases[voxelIndex] > 0) + { + ebsdlib::QuatD quat2(quats[neighborPoint * 4], quats[neighborPoint * 4 + 1], quats[neighborPoint * 4 + 2], quats[neighborPoint * 4 + 3]); + quat2.positiveOrientation(); + ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClass1]->calculateMisorientation(quat1, quat2); + if(axisAngle[3] < misorientationTolerance) + { + neighborCount[neighborPoint]++; + if(neighborCount[neighborPoint] >= currentLevel) + { + worklist.push_back(static_cast(neighborPoint)); + } + } + } + } + } + } + } + + return {}; +} diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.hpp new file mode 100644 index 0000000000..5bbba4e92c --- /dev/null +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.hpp @@ -0,0 +1,35 @@ +#pragma once + +#include "OrientationAnalysis/OrientationAnalysis_export.hpp" + +#include "simplnx/DataStructure/DataPath.hpp" +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ + +struct BadDataNeighborOrientationCheckInputValues; + +class ORIENTATIONANALYSIS_EXPORT BadDataNeighborOrientationCheckWorklist +{ +public: + BadDataNeighborOrientationCheckWorklist(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const BadDataNeighborOrientationCheckInputValues* inputValues); + ~BadDataNeighborOrientationCheckWorklist() noexcept; + + BadDataNeighborOrientationCheckWorklist(const BadDataNeighborOrientationCheckWorklist&) = delete; + BadDataNeighborOrientationCheckWorklist(BadDataNeighborOrientationCheckWorklist&&) noexcept = delete; + BadDataNeighborOrientationCheckWorklist& operator=(const BadDataNeighborOrientationCheckWorklist&) = delete; + BadDataNeighborOrientationCheckWorklist& operator=(BadDataNeighborOrientationCheckWorklist&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const BadDataNeighborOrientationCheckInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/OrientationAnalysis/test/BadDataNeighborOrientationCheckTest.cpp b/src/Plugins/OrientationAnalysis/test/BadDataNeighborOrientationCheckTest.cpp index 9b37c50a4c..bb398d365e 100644 --- a/src/Plugins/OrientationAnalysis/test/BadDataNeighborOrientationCheckTest.cpp +++ b/src/Plugins/OrientationAnalysis/test/BadDataNeighborOrientationCheckTest.cpp @@ -4,8 +4,12 @@ #include "OrientationAnalysis/OrientationAnalysis_test_dirs.hpp" #include "OrientationAnalysisTestUtils.hpp" +#include "simplnx/DataStructure/AttributeMatrix.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" #include "simplnx/UnitTest/UnitTestCommon.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" +#include #include namespace fs = std::filesystem; @@ -34,6 +38,11 @@ const DataPath k_CStuctsArrayPath = k_CellEnsembleDataPath.createChildPath(k_CSt // Case 1.1.1: Base Case | 2 phase | Tolerance 5 | 1 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.1.1", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_1/case_1_1_1/case_1_1_1_input.dream3d", unit_test::k_TestFilesDir)); @@ -92,6 +101,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.1. // Case 1.1.2: Invalid Base Case | 3 phase | Tolerance 5 | 1 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.1.2", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_1/case_1_1_2/case_1_1_2_input.dream3d", unit_test::k_TestFilesDir)); @@ -150,6 +164,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.1. // Case 1.1.3: Invalid Base Case | 2 phase | Tolerance 5 | 1 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.1.3", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_1/case_1_1_3/case_1_1_3_input.dream3d", unit_test::k_TestFilesDir)); @@ -208,6 +227,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.1. // Case 1.2.1: Base Case | 2 phase | Tolerance 5 | 2 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.2.1", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_2/case_1_2_1/case_1_2_1_input.dream3d", unit_test::k_TestFilesDir)); @@ -266,6 +290,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.2. // Case 1.2.2: Invalid Base Case | 2 phase | Tolerance 5 | 2 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.2.2", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_2/case_1_2_2/case_1_2_2_input.dream3d", unit_test::k_TestFilesDir)); @@ -324,6 +353,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.2. // Case 1.2.3: Invalid Base Case | 2 phase | Tolerance 5 | 2 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.2.3", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_2/case_1_2_3/case_1_2_3_input.dream3d", unit_test::k_TestFilesDir)); @@ -382,6 +416,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.2. // Case 1.3.1: Base Case | 1 phase | Tolerance 5 | 3 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.3.1", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_3/case_1_3_1/case_1_3_1_input.dream3d", unit_test::k_TestFilesDir)); @@ -440,6 +479,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.3. // Case 1.3.2: Invalid Base Case | 2 phase | Tolerance 5 | 3 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.3.2", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_3/case_1_3_2/case_1_3_2_input.dream3d", unit_test::k_TestFilesDir)); @@ -498,6 +542,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.3. // Case 1.3.3: Invalid Base Case | 1 phase | Tolerance 5 | 3 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.3.3", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_3/case_1_3_3/case_1_3_3_input.dream3d", unit_test::k_TestFilesDir)); @@ -556,6 +605,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.3. // Case 1.4.1: Base Case | 1 phase | Tolerance 5 | 4 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.4.1", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_4/case_1_4_1/case_1_4_1_input.dream3d", unit_test::k_TestFilesDir)); @@ -614,6 +668,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.4. // Case 1.4.2: Invalid Base Case | 2 phase | Tolerance 5 | 4 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.4.2", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_4/case_1_4_2/case_1_4_2_input.dream3d", unit_test::k_TestFilesDir)); @@ -672,6 +731,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.4. // Case 1.4.3: Invalid Base Case | 1 phase | Tolerance 5 | 4 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.4.3", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_4/case_1_4_3/case_1_4_3_input.dream3d", unit_test::k_TestFilesDir)); @@ -730,6 +794,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.4. // Case 1.5.1: Base Case | 1 phase | Tolerance 5 | 5 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.5.1", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_5/case_1_5_1/case_1_5_1_input.dream3d", unit_test::k_TestFilesDir)); @@ -788,6 +857,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.5. // Case 1.5.2: Invalid Base Case | 2 phase | Tolerance 5 | 5 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.5.2", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_5/case_1_5_2/case_1_5_2_input.dream3d", unit_test::k_TestFilesDir)); @@ -846,6 +920,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.5. // Case 1.5.3: Invalid Base Case | 1 phase | Tolerance 5 | 5 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.5.3", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_5/case_1_5_3/case_1_5_3_input.dream3d", unit_test::k_TestFilesDir)); @@ -904,6 +983,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.5. // Case 1.6.1: Base Case | 1 phase | Tolerance 5 | 6 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.6.1", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_6/case_1_6_1/case_1_6_1_input.dream3d", unit_test::k_TestFilesDir)); @@ -962,6 +1046,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.6. // Case 1.6.2: Invalid Base Case | 2 phase | Tolerance 5 | 6 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.6.2", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_6/case_1_6_2/case_1_6_2_input.dream3d", unit_test::k_TestFilesDir)); @@ -1020,6 +1109,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.6. // Case 1.6.3: Invalid Base Case | 1 phase | Tolerance 5 | 6 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.6.3", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 140, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_1/case_1_6/case_1_6_3/case_1_6_3_input.dream3d", unit_test::k_TestFilesDir)); @@ -1078,6 +1172,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 1.6. // Case 2.1: X+ Dim Case (Sequential) | Valid | 1 phase | Tolerance 5 | 5 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.1", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 400, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_2/case_2_1/case_2_1_input.dream3d", unit_test::k_TestFilesDir)); @@ -1132,6 +1231,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.1" // Case 2.2: Y+ Dim Case (Sequential) | Valid | 1 phase | Tolerance 5 | 5 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.2", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 400, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_2/case_2_2/case_2_2_input.dream3d", unit_test::k_TestFilesDir)); @@ -1186,6 +1290,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.2" // Case 2.3: Z+ Dim Case (Sequential) | Valid | 1 phase | Tolerance 5 | 5 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.3", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 400, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_2/case_2_3/case_2_3_input.dream3d", unit_test::k_TestFilesDir)); @@ -1240,6 +1349,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.3" // Case 2.4: X- Dim Case (Recursive) | Valid | 1 phase | Tolerance 5 | 5 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.4", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 400, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_2/case_2_4/case_2_4_input.dream3d", unit_test::k_TestFilesDir)); @@ -1294,6 +1408,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.4" // Case 2.5: Y- Dim Case (Recursive) | Valid | 1 phase | Tolerance 5 | 5 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.5", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 400, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_2/case_2_5/case_2_5_input.dream3d", unit_test::k_TestFilesDir)); @@ -1348,6 +1467,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.5" // Case 2.6: Z- Dim Case (Recursive) | Valid | 1 phase | Tolerance 5 | 5 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.6", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 400, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_2/case_2_6/case_2_6_input.dream3d", unit_test::k_TestFilesDir)); @@ -1402,6 +1526,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 2.6" // Case 3.1: Long Sequential | Valid | 1 phase | Tolerance 5 | 1 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 3.1", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 400, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_3/case_3_1/case_3_1_input.dream3d", unit_test::k_TestFilesDir)); @@ -1456,6 +1585,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 3.1" // Case 3.2: Long Recursive | Valid | 1 phase | Tolerance 5 | 1 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 3.2", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 400, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_3/case_3_2/case_3_2_input.dream3d", unit_test::k_TestFilesDir)); @@ -1510,6 +1644,11 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 3.2" // Case 4: Semi-Complex Synthetic Structure | Valid | 3 phase | Tolerance 5 | 4 Min Neighbors TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 4", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter]") { + UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 400, true); + const UnitTest::TestFileSentinel testDataSentinel(unit_test::k_TestFilesDir, "7_bad_data_neighbor_orientation_check.tar.gz", "bad_data_neighbor_orientation_check"); auto baseDataFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check/case_4/case_4_input.dream3d", unit_test::k_TestFilesDir)); @@ -1599,3 +1738,116 @@ TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Case 4", UnitTest::CheckArraysInheritTupleDims(dataStructure); } + +// Benchmark: 200x200x200 programmatic dataset (no exemplar — correctness proven by the 27 tests above) +// Writes data to a .dream3d file and reloads it so arrays use disk-backed ZarrStore in OOC mode. +TEST_CASE("OrientationAnalysis::BadDataNeighborOrientationCheckFilter: Benchmark 200x200x200", "[OrientationAnalysis][BadDataNeighborOrientationCheckFilter][Benchmark]") +{ + UnitTest::LoadPlugins(); + // 1 Z-slice of quats: 200*200*4*4 = 640000 bytes → ~200 chunks in OOC + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 10000, true); + + constexpr usize kDimX = 200; + constexpr usize kDimY = 200; + constexpr usize kDimZ = 200; + constexpr usize kTotalVoxels = kDimX * kDimY * kDimZ; + + const auto benchmarkFilePath = fs::path(fmt::format("{}/bad_data_neighbor_orientation_check_benchmark.dream3d", unit_test::k_BinaryTestOutputDir)); + + // ===== Stage 1: Build the dataset programmatically and write to .dream3d ===== + { + DataStructure buildDS; + + ImageGeom* imageGeom = ImageGeom::Create(buildDS, VerificationConstants::k_ImageName); + imageGeom->setDimensions({kDimX, kDimY, kDimZ}); + imageGeom->setSpacing({1.0f, 1.0f, 1.0f}); + imageGeom->setOrigin({0.0f, 0.0f, 0.0f}); + + AttributeMatrix* cellAM = AttributeMatrix::Create(buildDS, Constants::k_Cell_Data, {kDimZ, kDimY, kDimX}, imageGeom->getId()); + imageGeom->setCellData(*cellAM); + + AttributeMatrix* ensembleAM = AttributeMatrix::Create(buildDS, Constants::k_Cell_Ensemble_Data, {2}, imageGeom->getId()); + + // Crystal Structures: [999 (Unknown), 1 (Cubic_High)] + auto* crystalStructures = UnitTest::CreateTestDataArray(buildDS, VerificationConstants::k_CStuctsName, {2}, {1}, ensembleAM->getId()); + auto& csStore = crystalStructures->getDataStoreRef(); + csStore.setValue(0, 999); + csStore.setValue(1, 1); + + const ShapeType cellTupleShape = {kDimZ, kDimY, kDimX}; + + // Phases: all phase 1 + auto* phases = UnitTest::CreateTestDataArray(buildDS, VerificationConstants::k_PhasesName, cellTupleShape, {1}, cellAM->getId()); + auto& phasesStore = phases->getDataStoreRef(); + for(usize i = 0; i < kTotalVoxels; i++) + { + phasesStore.setValue(i, 1); + } + + // Mask: ~40% bad voxels (mask=0) deterministically scattered to stress Phase 2 + auto* mask = UnitTest::CreateTestDataArray(buildDS, VerificationConstants::k_MaskName, cellTupleShape, {1}, cellAM->getId()); + auto& maskStore = mask->getDataStoreRef(); + for(usize i = 0; i < kTotalVoxels; i++) + { + maskStore.setValue(i, static_cast(((i * 7 + 13) % 100) >= 40 ? 1 : 0)); + } + + // Quats: 4-component float32 with spatially-varying orientations. + // 8 octants with distinct base quaternions (15 deg apart about Z axis). + // Intra-octant perturbations are <5 deg. Cross-boundary neighbors differ by >5 deg, + // creating varied neighbor counts (1-6) that exercise Phase 2 cascading. + auto* quats = UnitTest::CreateTestDataArray(buildDS, VerificationConstants::k_QuatsName, cellTupleShape, {4}, cellAM->getId()); + auto& quatsStore = quats->getDataStoreRef(); + + constexpr usize kHalfX = kDimX / 2; + constexpr usize kHalfY = kDimY / 2; + constexpr usize kHalfZ = kDimZ / 2; + const float32 baseAngles[8] = {0.0f, 15.0f, 30.0f, 45.0f, 60.0f, 75.0f, 90.0f, 105.0f}; + + for(usize i = 0; i < kTotalVoxels; i++) + { + const usize ix = i % kDimX; + const usize iy = (i / kDimX) % kDimY; + const usize iz = i / (kDimX * kDimY); + const usize octant = (iz >= kHalfZ ? 4 : 0) + (iy >= kHalfY ? 2 : 0) + (ix >= kHalfX ? 1 : 0); + + const float32 perturbDeg = 0.5f * static_cast(i % 10); + const float32 angleDeg = baseAngles[octant] + perturbDeg; + const float32 halfAngle = angleDeg * 3.14159265f / 360.0f; + + const float32 w = std::cos(halfAngle); + const float32 z = std::sin(halfAngle); + quatsStore.setValue(i * 4 + 0, w); + quatsStore.setValue(i * 4 + 1, 0.0f); + quatsStore.setValue(i * 4 + 2, 0.0f); + quatsStore.setValue(i * 4 + 3, z); + } + + UnitTest::WriteTestDataStructure(buildDS, benchmarkFilePath); + } + + // ===== Stage 2: Reload from .dream3d (arrays become ZarrStore in OOC mode) and run filter ===== + DataStructure dataStructure = UnitTest::LoadDataStructure(benchmarkFilePath); + + { + BadDataNeighborOrientationCheckFilter filter; + Arguments args; + + args.insertOrAssign(BadDataNeighborOrientationCheckFilter::k_MisorientationTolerance_Key, std::make_any(5.0f)); + args.insertOrAssign(BadDataNeighborOrientationCheckFilter::k_NumberOfNeighbors_Key, std::make_any(1)); + args.insertOrAssign(BadDataNeighborOrientationCheckFilter::k_ImageGeometryPath_Key, std::make_any(VerificationConstants::k_ImagePath)); + args.insertOrAssign(BadDataNeighborOrientationCheckFilter::k_MaskArrayPath_Key, std::make_any(VerificationConstants::k_MaskArrayPath)); + args.insertOrAssign(BadDataNeighborOrientationCheckFilter::k_CellPhasesArrayPath_Key, std::make_any(VerificationConstants::k_PhasesArrayPath)); + args.insertOrAssign(BadDataNeighborOrientationCheckFilter::k_QuatsArrayPath_Key, std::make_any(VerificationConstants::k_QuatsArrayPath)); + args.insertOrAssign(BadDataNeighborOrientationCheckFilter::k_CrystalStructuresArrayPath_Key, std::make_any(VerificationConstants::k_CStuctsArrayPath)); + + auto preflightResult = filter.preflight(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions); + + auto executeResult = filter.execute(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result); + } + + // Clean up the temporary benchmark file + fs::remove(benchmarkFilePath); +} From 27be61698de7c2995985e659b2a46bc79e70c316 Mon Sep 17 00:00:00 2001 From: Joey Kleingers Date: Thu, 5 Mar 2026 13:56:28 -0500 Subject: [PATCH 2/5] ENH: Split Group B face-neighbor filters into Direct/Scanline OOC dispatch Split ComputeBoundaryCells, ComputeSurfaceFeatures, ComputeFeatureNeighbors, and ComputeSurfaceAreaToVolume into Direct (in-core) and Scanline (OOC) algorithm classes using DispatchAlgorithm pattern. Direct classes preserve original code for zero in-core regression. Scanline classes add chunk-sequential iteration for OOC storage backends. Added benchmark tests (200x200x200) and OOC test coverage with PreferencesSentinel + ForceOocAlgorithmGuard for all four filters. Signed-off-by: Joey Kleingers --- src/Plugins/SimplnxCore/CMakeLists.txt | 8 + .../Algorithms/ComputeBoundaryCellsDirect.cpp | 109 ++++++++ .../Algorithms/ComputeBoundaryCellsDirect.hpp | 32 +++ .../ComputeBoundaryCellsScanline.cpp | 112 ++++++++ .../ComputeBoundaryCellsScanline.hpp | 32 +++ .../ComputeFeatureNeighborsDirect.cpp | 224 ++++++++++++++++ .../ComputeFeatureNeighborsDirect.hpp | 32 +++ .../ComputeFeatureNeighborsScanline.cpp | 239 ++++++++++++++++++ .../ComputeFeatureNeighborsScanline.hpp | 33 +++ .../ComputeSurfaceAreaToVolumeDirect.cpp | 155 ++++++++++++ .../ComputeSurfaceAreaToVolumeDirect.hpp | 33 +++ .../ComputeSurfaceAreaToVolumeScanline.cpp | 167 ++++++++++++ .../ComputeSurfaceAreaToVolumeScanline.hpp | 33 +++ .../ComputeSurfaceFeaturesDirect.cpp | 230 +++++++++++++++++ .../ComputeSurfaceFeaturesDirect.hpp | 33 +++ .../ComputeSurfaceFeaturesScanline.cpp | 226 +++++++++++++++++ .../ComputeSurfaceFeaturesScanline.hpp | 33 +++ .../test/ComputeBoundaryCellsTest.cpp | 67 +++++ .../test/ComputeFeatureNeighborsTest.cpp | 75 ++++++ .../test/ComputeSurfaceAreaToVolumeTest.cpp | 88 +++++++ .../test/ComputeSurfaceFeaturesTest.cpp | 78 ++++++ 21 files changed, 2039 insertions(+) create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.hpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.hpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.hpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.hpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.hpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.cpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.cpp create mode 100644 src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp diff --git a/src/Plugins/SimplnxCore/CMakeLists.txt b/src/Plugins/SimplnxCore/CMakeLists.txt index 8175d378fe..b67008225b 100644 --- a/src/Plugins/SimplnxCore/CMakeLists.txt +++ b/src/Plugins/SimplnxCore/CMakeLists.txt @@ -186,6 +186,8 @@ set(AlgorithmList ComputeArrayStatistics ComputeBiasedFeatures ComputeBoundaryCells + ComputeBoundaryCellsDirect + ComputeBoundaryCellsScanline ComputeBoundingBoxStats ComputeCoordinatesImageGeom ComputeCoordinateThreshold @@ -196,6 +198,8 @@ set(AlgorithmList ComputeFeatureCentroids ComputeFeatureClustering ComputeFeatureNeighbors + ComputeFeatureNeighborsDirect + ComputeFeatureNeighborsScanline ComputeFeaturePhases ComputeFeaturePhasesBinary ComputeFeatureRect @@ -209,7 +213,11 @@ set(AlgorithmList ComputeNeighborListStatistics # ComputeNumFeatures ComputeSurfaceAreaToVolume + ComputeSurfaceAreaToVolumeDirect + ComputeSurfaceAreaToVolumeScanline ComputeSurfaceFeatures + ComputeSurfaceFeaturesDirect + ComputeSurfaceFeaturesScanline # ComputeTriangleAreas ComputeTriangleGeomCentroids ComputeTriangleGeomVolumes diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp new file mode 100644 index 0000000000..e3c9e3e58d --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp @@ -0,0 +1,109 @@ +#include "ComputeBoundaryCellsDirect.hpp" + +#include "ComputeBoundaryCells.hpp" + +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/Utilities/NeighborUtilities.hpp" + +using namespace nx::core; + +// ----------------------------------------------------------------------------- +ComputeBoundaryCellsDirect::ComputeBoundaryCellsDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeBoundaryCellsInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeBoundaryCellsDirect::~ComputeBoundaryCellsDirect() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> ComputeBoundaryCellsDirect::operator()() +{ + const ImageGeom imageGeometry = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); + const SizeVec3 udims = imageGeometry.getDimensions(); + std::array dims = { + static_cast(udims[0]), + static_cast(udims[1]), + static_cast(udims[2]), + }; + + auto& featureIdsStore = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath)->getDataStoreRef(); + auto& boundaryCellsStore = m_DataStructure.getDataAs(m_InputValues->BoundaryCellsArrayName)->getDataStoreRef(); + + std::array neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims); + std::array faceNeighborInternalIdx = initializeFaceNeighborInternalIdx(); + + int32 feature = 0; + int8 onSurf = 0; + int64 neighborPoint = 0; + + int ignoreFeatureZeroVal = 0; + if(!m_InputValues->IgnoreFeatureZero) + { + ignoreFeatureZeroVal = -1; + } + + int64 kStride = 0; + int64 jStride = 0; + + for(int64 zIdx = 0; zIdx < dims[2]; zIdx++) + { + kStride = dims[0] * dims[1] * zIdx; + for(int64 yIdx = 0; yIdx < dims[1]; yIdx++) + { + jStride = dims[0] * yIdx; + for(int64 xIdx = 0; xIdx < dims[0]; xIdx++) + { + int64 voxelIndex = kStride + jStride + xIdx; + onSurf = 0; + feature = featureIdsStore[voxelIndex]; + if(feature >= 0) + { + if(m_InputValues->IncludeVolumeBoundary) + { + if(dims[0] > 2 && (xIdx == 0 || xIdx == dims[0] - 1)) + { + onSurf++; + } + if(dims[1] > 2 && (yIdx == 0 || yIdx == dims[1] - 1)) + { + onSurf++; + } + if(dims[2] > 2 && (zIdx == 0 || zIdx == dims[2] - 1)) + { + onSurf++; + } + + if(onSurf > 0 && feature == 0) + { + onSurf = 0; + } + } + + // Loop over the 6 face neighbors of the voxel + std::array isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims); + for(const auto& faceIndex : faceNeighborInternalIdx) + { + if(!isValidFaceNeighbor[faceIndex]) + { + continue; + } + neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; + + if(featureIdsStore[neighborPoint] != feature && featureIdsStore[neighborPoint] > ignoreFeatureZeroVal) + { + onSurf++; + } + } + } + boundaryCellsStore[voxelIndex] = onSurf; + } + } + } + return {}; +} diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.hpp new file mode 100644 index 0000000000..35f6a06bc7 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include "SimplnxCore/SimplnxCore_export.hpp" + +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ +struct ComputeBoundaryCellsInputValues; + +class SIMPLNXCORE_EXPORT ComputeBoundaryCellsDirect +{ +public: + ComputeBoundaryCellsDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, const ComputeBoundaryCellsInputValues* inputValues); + ~ComputeBoundaryCellsDirect() noexcept; + + ComputeBoundaryCellsDirect(const ComputeBoundaryCellsDirect&) = delete; + ComputeBoundaryCellsDirect(ComputeBoundaryCellsDirect&&) noexcept = delete; + ComputeBoundaryCellsDirect& operator=(const ComputeBoundaryCellsDirect&) = delete; + ComputeBoundaryCellsDirect& operator=(ComputeBoundaryCellsDirect&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const ComputeBoundaryCellsInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp new file mode 100644 index 0000000000..715d9a8055 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp @@ -0,0 +1,112 @@ +#include "ComputeBoundaryCellsScanline.hpp" + +#include "ComputeBoundaryCells.hpp" + +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/Utilities/NeighborUtilities.hpp" + +using namespace nx::core; + +// ----------------------------------------------------------------------------- +ComputeBoundaryCellsScanline::ComputeBoundaryCellsScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeBoundaryCellsInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeBoundaryCellsScanline::~ComputeBoundaryCellsScanline() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> ComputeBoundaryCellsScanline::operator()() +{ + const ImageGeom imageGeometry = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); + const SizeVec3 udims = imageGeometry.getDimensions(); + std::array dims = { + static_cast(udims[0]), + static_cast(udims[1]), + static_cast(udims[2]), + }; + + auto& featureIdsStore = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath)->getDataStoreRef(); + auto& boundaryCellsStore = m_DataStructure.getDataAs(m_InputValues->BoundaryCellsArrayName)->getDataStoreRef(); + + std::array neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims); + std::array faceNeighborInternalIdx = initializeFaceNeighborInternalIdx(); + + int ignoreFeatureZeroVal = 0; + if(!m_InputValues->IgnoreFeatureZero) + { + ignoreFeatureZeroVal = -1; + } + + const uint64 numChunks = featureIdsStore.getNumberOfChunks(); + + for(uint64 chunkIdx = 0; chunkIdx < numChunks; chunkIdx++) + { + featureIdsStore.loadChunk(chunkIdx); + + const auto chunkLowerBounds = featureIdsStore.getChunkLowerBounds(chunkIdx); + const auto chunkUpperBounds = featureIdsStore.getChunkUpperBounds(chunkIdx); + + for(usize zIdx = chunkLowerBounds[0]; zIdx <= chunkUpperBounds[0]; zIdx++) + { + const int64 kStride = static_cast(zIdx) * dims[0] * dims[1]; + for(usize yIdx = chunkLowerBounds[1]; yIdx <= chunkUpperBounds[1]; yIdx++) + { + const int64 jStride = static_cast(yIdx) * dims[0]; + for(usize xIdx = chunkLowerBounds[2]; xIdx <= chunkUpperBounds[2]; xIdx++) + { + const int64 voxelIndex = kStride + jStride + static_cast(xIdx); + int8 onSurf = 0; + const int32 feature = featureIdsStore[voxelIndex]; + if(feature >= 0) + { + if(m_InputValues->IncludeVolumeBoundary) + { + if(dims[0] > 2 && (static_cast(xIdx) == 0 || static_cast(xIdx) == dims[0] - 1)) + { + onSurf++; + } + if(dims[1] > 2 && (static_cast(yIdx) == 0 || static_cast(yIdx) == dims[1] - 1)) + { + onSurf++; + } + if(dims[2] > 2 && (static_cast(zIdx) == 0 || static_cast(zIdx) == dims[2] - 1)) + { + onSurf++; + } + + if(onSurf > 0 && feature == 0) + { + onSurf = 0; + } + } + + // Loop over the 6 face neighbors of the voxel + std::array isValidFaceNeighbor = computeValidFaceNeighbors(static_cast(xIdx), static_cast(yIdx), static_cast(zIdx), dims); + for(const auto& faceIndex : faceNeighborInternalIdx) + { + if(!isValidFaceNeighbor[faceIndex]) + { + continue; + } + const int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; + + if(featureIdsStore[neighborPoint] != feature && featureIdsStore[neighborPoint] > ignoreFeatureZeroVal) + { + onSurf++; + } + } + } + boundaryCellsStore[voxelIndex] = onSurf; + } + } + } + } + return {}; +} diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.hpp new file mode 100644 index 0000000000..9ca5742150 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include "SimplnxCore/SimplnxCore_export.hpp" + +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ +struct ComputeBoundaryCellsInputValues; + +class SIMPLNXCORE_EXPORT ComputeBoundaryCellsScanline +{ +public: + ComputeBoundaryCellsScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, const ComputeBoundaryCellsInputValues* inputValues); + ~ComputeBoundaryCellsScanline() noexcept; + + ComputeBoundaryCellsScanline(const ComputeBoundaryCellsScanline&) = delete; + ComputeBoundaryCellsScanline(ComputeBoundaryCellsScanline&&) noexcept = delete; + ComputeBoundaryCellsScanline& operator=(const ComputeBoundaryCellsScanline&) = delete; + ComputeBoundaryCellsScanline& operator=(ComputeBoundaryCellsScanline&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const ComputeBoundaryCellsInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp new file mode 100644 index 0000000000..da14331911 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp @@ -0,0 +1,224 @@ +#include "ComputeFeatureNeighborsDirect.hpp" + +#include "ComputeFeatureNeighbors.hpp" + +#include "simplnx/DataStructure/AttributeMatrix.hpp" +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/DataStructure/NeighborList.hpp" +#include "simplnx/Utilities/MessageHelper.hpp" +#include "simplnx/Utilities/NeighborUtilities.hpp" + +#include + +using namespace nx::core; + +// ----------------------------------------------------------------------------- +ComputeFeatureNeighborsDirect::ComputeFeatureNeighborsDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeFeatureNeighborsInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeFeatureNeighborsDirect::~ComputeFeatureNeighborsDirect() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> ComputeFeatureNeighborsDirect::operator()() +{ + auto storeBoundaryCells = m_InputValues->StoreBoundaryCells; + auto storeSurfaceFeatures = m_InputValues->StoreSurfaceFeatures; + auto imageGeomPath = m_InputValues->InputImageGeometryPath; + auto featureIdsPath = m_InputValues->FeatureIdsPath; + auto boundaryCellsName = m_InputValues->BoundaryCellsName; + auto numNeighborsName = m_InputValues->NumberOfNeighborsName; + auto neighborListName = m_InputValues->NeighborListName; + auto sharedSurfaceAreaName = m_InputValues->SharedSurfaceAreaListName; + auto surfaceFeaturesName = m_InputValues->SurfaceFeaturesName; + auto featureAttrMatrixPath = m_InputValues->CellFeatureArrayPath; + + DataPath boundaryCellsPath = featureIdsPath.replaceName(boundaryCellsName); + DataPath numNeighborsPath = featureAttrMatrixPath.createChildPath(numNeighborsName); + DataPath neighborListPath = featureAttrMatrixPath.createChildPath(neighborListName); + DataPath sharedSurfaceAreaPath = featureAttrMatrixPath.createChildPath(sharedSurfaceAreaName); + DataPath surfaceFeaturesPath = featureAttrMatrixPath.createChildPath(surfaceFeaturesName); + + auto& featureIds = m_DataStructure.getDataAs(featureIdsPath)->getDataStoreRef(); + auto& numNeighbors = m_DataStructure.getDataAs(numNeighborsPath)->getDataStoreRef(); + + auto& neighborList = m_DataStructure.getDataRefAs(neighborListPath); + auto& sharedSurfaceAreaList = m_DataStructure.getDataRefAs(sharedSurfaceAreaPath); + + auto* boundaryCells = storeBoundaryCells ? m_DataStructure.getDataAs(boundaryCellsPath)->getDataStore() : nullptr; + auto* surfaceFeatures = storeSurfaceFeatures ? m_DataStructure.getDataAs(surfaceFeaturesPath)->getDataStore() : nullptr; + + usize totalPoints = featureIds.getNumberOfTuples(); + usize totalFeatures = numNeighbors.getNumberOfTuples(); + + /* Ensure that we will be able to work with the user selected featureId Array */ + const auto [minFeatureId, maxFeatureId] = std::minmax_element(featureIds.begin(), featureIds.end()); + if(static_cast(*maxFeatureId) >= totalFeatures) + { + std::stringstream out; + out << "Data Array " << featureIdsPath.getTargetName() << " has a maximum value of " << *maxFeatureId << " which is greater than the " + << " number of features from array " << numNeighborsPath.getTargetName() << " which has " << totalFeatures << ". Did you select the " + << " incorrect array for the 'FeatureIds' array?"; + return MakeErrorResult(-24500, out.str()); + } + + auto& imageGeom = m_DataStructure.getDataRefAs(imageGeomPath); + SizeVec3 uDims = imageGeom.getDimensions(); + + std::array dims = { + static_cast(uDims[0]), + static_cast(uDims[1]), + static_cast(uDims[2]), + }; + + std::array neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims); + std::array faceNeighborInternalIdx = initializeFaceNeighborInternalIdx(); + + int32 feature = 0; + int32 nnum = 0; + uint8 onsurf = 0; + + std::vector> neighborlist(totalFeatures); + std::vector> neighborsurfacearealist(totalFeatures); + + int32 nListSize = 100; + + MessageHelper messageHelper(m_MessageHandler); + ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger(); + // Initialize the neighbor lists + for(usize featureIdx = 1; featureIdx < totalFeatures; featureIdx++) + { + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Initializing Neighbor Lists || {:.2f}% Complete", CalculatePercentComplete(featureIdx, totalFeatures)); }); + + if(m_ShouldCancel) + { + return {}; + } + + numNeighbors[featureIdx] = 0; + neighborlist[featureIdx].resize(nListSize); + neighborsurfacearealist[featureIdx].assign(nListSize, -1.0f); + if(storeSurfaceFeatures && surfaceFeatures != nullptr) + { + surfaceFeatures->setValue(featureIdx, false); + } + } + + // Loop over all points to generate the neighbor lists + for(int64 voxelIndex = 0; voxelIndex < totalPoints; voxelIndex++) + { + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Determining Neighbor Lists || {:.2f}% Complete", CalculatePercentComplete(voxelIndex, totalPoints)); }); + + if(m_ShouldCancel) + { + return {}; + } + + onsurf = 0; + feature = featureIds[voxelIndex]; + if(feature > 0 && feature < neighborlist.size()) + { + int64 xIdx = voxelIndex % dims[0]; + int64 yIdx = (voxelIndex / dims[0]) % dims[1]; + int64 zIdx = voxelIndex / (dims[0] * dims[1]); + + if(storeSurfaceFeatures && surfaceFeatures != nullptr) + { + if((xIdx == 0 || xIdx == static_cast((dims[0] - 1)) || yIdx == 0 || yIdx == static_cast((dims[1]) - 1) || zIdx == 0 || zIdx == static_cast((dims[2] - 1))) && dims[2] != 1) + { + surfaceFeatures->setValue(feature, true); + } + if((xIdx == 0 || xIdx == static_cast((dims[0] - 1)) || yIdx == 0 || yIdx == static_cast((dims[1] - 1))) && dims[2] == 1) + { + surfaceFeatures->setValue(feature, true); + } + } + + // Loop over the 6 face neighbors of the voxel + std::array isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims); + for(const auto& faceIndex : faceNeighborInternalIdx) + { + if(!isValidFaceNeighbor[faceIndex]) + { + continue; + } + + const int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; + + if(featureIds[neighborPoint] != feature && featureIds[neighborPoint] > 0) + { + onsurf++; + nnum = numNeighbors[feature]; + neighborlist[feature].push_back(featureIds[neighborPoint]); + nnum++; + numNeighbors[feature] = nnum; + } + } + } + if(storeBoundaryCells && boundaryCells != nullptr) + { + boundaryCells->setValue(voxelIndex, static_cast(onsurf)); + } + } + + FloatVec3 spacing = imageGeom.getSpacing(); + + // We do this to create new set of NeighborList objects + for(usize i = 1; i < totalFeatures; i++) + { + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Calculating Surface Areas || {:.2f}% Complete", CalculatePercentComplete(i, totalFeatures)); }); + + if(m_ShouldCancel) + { + return {}; + } + + std::map neighToCount; + auto numneighs = static_cast(neighborlist[i].size()); + + // this increments the voxel counts for each feature + for(int32 j = 0; j < numneighs; j++) + { + neighToCount[neighborlist[i][j]]++; + } + + auto neighborIter = neighToCount.find(0); + neighToCount.erase(neighborIter); + neighborIter = neighToCount.find(-1); + if(neighborIter != neighToCount.end()) + { + neighToCount.erase(neighborIter); + } + // Resize the features neighbor list to zero + neighborlist[i].resize(0); + neighborsurfacearealist[i].resize(0); + + for(const auto [neigh, number] : neighToCount) + { + float area = static_cast(number) * spacing[0] * spacing[1]; + + // Push the neighbor feature identifier back onto the list, so we stay synced up + neighborlist[i].push_back(neigh); + neighborsurfacearealist[i].push_back(area); + } + numNeighbors[i] = static_cast(neighborlist[i].size()); + + // Set the vector for each list into the NeighborList Object + NeighborList::SharedVectorType sharedNeiLst(new std::vector); + sharedNeiLst->assign(neighborlist[i].begin(), neighborlist[i].end()); + neighborList.setList(static_cast(i), sharedNeiLst); + + NeighborList::SharedVectorType sharedSAL(new std::vector); + sharedSAL->assign(neighborsurfacearealist[i].begin(), neighborsurfacearealist[i].end()); + sharedSurfaceAreaList.setList(static_cast(i), sharedSAL); + } + + return {}; +} diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.hpp new file mode 100644 index 0000000000..aace0b087c --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include "SimplnxCore/SimplnxCore_export.hpp" + +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ +struct ComputeFeatureNeighborsInputValues; + +class SIMPLNXCORE_EXPORT ComputeFeatureNeighborsDirect +{ +public: + ComputeFeatureNeighborsDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, const ComputeFeatureNeighborsInputValues* inputValues); + ~ComputeFeatureNeighborsDirect() noexcept; + + ComputeFeatureNeighborsDirect(const ComputeFeatureNeighborsDirect&) = delete; + ComputeFeatureNeighborsDirect(ComputeFeatureNeighborsDirect&&) noexcept = delete; + ComputeFeatureNeighborsDirect& operator=(const ComputeFeatureNeighborsDirect&) = delete; + ComputeFeatureNeighborsDirect& operator=(ComputeFeatureNeighborsDirect&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const ComputeFeatureNeighborsInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp new file mode 100644 index 0000000000..c081fed686 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp @@ -0,0 +1,239 @@ +#include "ComputeFeatureNeighborsScanline.hpp" + +#include "ComputeFeatureNeighbors.hpp" + +#include "simplnx/DataStructure/AttributeMatrix.hpp" +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/DataStructure/NeighborList.hpp" +#include "simplnx/Utilities/MessageHelper.hpp" +#include "simplnx/Utilities/NeighborUtilities.hpp" + +#include + +using namespace nx::core; + +// ----------------------------------------------------------------------------- +ComputeFeatureNeighborsScanline::ComputeFeatureNeighborsScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeFeatureNeighborsInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeFeatureNeighborsScanline::~ComputeFeatureNeighborsScanline() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> ComputeFeatureNeighborsScanline::operator()() +{ + auto storeBoundaryCells = m_InputValues->StoreBoundaryCells; + auto storeSurfaceFeatures = m_InputValues->StoreSurfaceFeatures; + auto imageGeomPath = m_InputValues->InputImageGeometryPath; + auto featureIdsPath = m_InputValues->FeatureIdsPath; + auto boundaryCellsName = m_InputValues->BoundaryCellsName; + auto numNeighborsName = m_InputValues->NumberOfNeighborsName; + auto neighborListName = m_InputValues->NeighborListName; + auto sharedSurfaceAreaName = m_InputValues->SharedSurfaceAreaListName; + auto surfaceFeaturesName = m_InputValues->SurfaceFeaturesName; + auto featureAttrMatrixPath = m_InputValues->CellFeatureArrayPath; + + DataPath boundaryCellsPath = featureIdsPath.replaceName(boundaryCellsName); + DataPath numNeighborsPath = featureAttrMatrixPath.createChildPath(numNeighborsName); + DataPath neighborListPath = featureAttrMatrixPath.createChildPath(neighborListName); + DataPath sharedSurfaceAreaPath = featureAttrMatrixPath.createChildPath(sharedSurfaceAreaName); + DataPath surfaceFeaturesPath = featureAttrMatrixPath.createChildPath(surfaceFeaturesName); + + auto& featureIds = m_DataStructure.getDataAs(featureIdsPath)->getDataStoreRef(); + auto& numNeighbors = m_DataStructure.getDataAs(numNeighborsPath)->getDataStoreRef(); + + auto& neighborList = m_DataStructure.getDataRefAs(neighborListPath); + auto& sharedSurfaceAreaList = m_DataStructure.getDataRefAs(sharedSurfaceAreaPath); + + auto* boundaryCells = storeBoundaryCells ? m_DataStructure.getDataAs(boundaryCellsPath)->getDataStore() : nullptr; + auto* surfaceFeatures = storeSurfaceFeatures ? m_DataStructure.getDataAs(surfaceFeaturesPath)->getDataStore() : nullptr; + + usize totalPoints = featureIds.getNumberOfTuples(); + usize totalFeatures = numNeighbors.getNumberOfTuples(); + + /* Ensure that we will be able to work with the user selected featureId Array */ + const auto [minFeatureId, maxFeatureId] = std::minmax_element(featureIds.begin(), featureIds.end()); + if(static_cast(*maxFeatureId) >= totalFeatures) + { + std::stringstream out; + out << "Data Array " << featureIdsPath.getTargetName() << " has a maximum value of " << *maxFeatureId << " which is greater than the " + << " number of features from array " << numNeighborsPath.getTargetName() << " which has " << totalFeatures << ". Did you select the " + << " incorrect array for the 'FeatureIds' array?"; + return MakeErrorResult(-24500, out.str()); + } + + auto& imageGeom = m_DataStructure.getDataRefAs(imageGeomPath); + SizeVec3 uDims = imageGeom.getDimensions(); + + std::array dims = { + static_cast(uDims[0]), + static_cast(uDims[1]), + static_cast(uDims[2]), + }; + + std::array neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims); + std::array faceNeighborInternalIdx = initializeFaceNeighborInternalIdx(); + + int32 feature = 0; + int32 nnum = 0; + uint8 onsurf = 0; + + std::vector> neighborlist(totalFeatures); + std::vector> neighborsurfacearealist(totalFeatures); + + int32 nListSize = 100; + + MessageHelper messageHelper(m_MessageHandler); + ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger(); + // Initialize the neighbor lists + for(usize featureIdx = 1; featureIdx < totalFeatures; featureIdx++) + { + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Initializing Neighbor Lists || {:.2f}% Complete", CalculatePercentComplete(featureIdx, totalFeatures)); }); + + if(m_ShouldCancel) + { + return {}; + } + + numNeighbors[featureIdx] = 0; + neighborlist[featureIdx].resize(nListSize); + neighborsurfacearealist[featureIdx].assign(nListSize, -1.0f); + if(storeSurfaceFeatures && surfaceFeatures != nullptr) + { + surfaceFeatures->setValue(featureIdx, false); + } + } + + // Phase 1: Loop over all points using chunk-sequential iteration + const uint64 numChunks = featureIds.getNumberOfChunks(); + + for(uint64 chunkIdx = 0; chunkIdx < numChunks; chunkIdx++) + { + featureIds.loadChunk(chunkIdx); + + const auto chunkLowerBounds = featureIds.getChunkLowerBounds(chunkIdx); + const auto chunkUpperBounds = featureIds.getChunkUpperBounds(chunkIdx); + + for(usize zIdx = chunkLowerBounds[0]; zIdx <= chunkUpperBounds[0]; zIdx++) + { + const int64 kStride = static_cast(zIdx) * dims[0] * dims[1]; + for(usize yIdx = chunkLowerBounds[1]; yIdx <= chunkUpperBounds[1]; yIdx++) + { + const int64 jStride = static_cast(yIdx) * dims[0]; + for(usize xIdx = chunkLowerBounds[2]; xIdx <= chunkUpperBounds[2]; xIdx++) + { + throttledMessenger.sendThrottledMessage( + [&]() { return fmt::format("Determining Neighbor Lists || Chunk {}/{}", chunkIdx + 1, numChunks); }); + + if(m_ShouldCancel) + { + return {}; + } + + const int64 voxelIndex = kStride + jStride + static_cast(xIdx); + onsurf = 0; + feature = featureIds[voxelIndex]; + if(feature > 0 && feature < neighborlist.size()) + { + const int64 ix = static_cast(xIdx); + const int64 iy = static_cast(yIdx); + const int64 iz = static_cast(zIdx); + + if(storeSurfaceFeatures && surfaceFeatures != nullptr) + { + if((ix == 0 || ix == dims[0] - 1 || iy == 0 || iy == dims[1] - 1 || iz == 0 || iz == dims[2] - 1) && dims[2] != 1) + { + surfaceFeatures->setValue(feature, true); + } + if((ix == 0 || ix == dims[0] - 1 || iy == 0 || iy == dims[1] - 1) && dims[2] == 1) + { + surfaceFeatures->setValue(feature, true); + } + } + + // Loop over the 6 face neighbors of the voxel + std::array isValidFaceNeighbor = computeValidFaceNeighbors(ix, iy, iz, dims); + for(const auto& faceIndex : faceNeighborInternalIdx) + { + if(!isValidFaceNeighbor[faceIndex]) + { + continue; + } + + const int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; + + if(featureIds[neighborPoint] != feature && featureIds[neighborPoint] > 0) + { + onsurf++; + nnum = numNeighbors[feature]; + neighborlist[feature].push_back(featureIds[neighborPoint]); + nnum++; + numNeighbors[feature] = nnum; + } + } + } + if(storeBoundaryCells && boundaryCells != nullptr) + { + boundaryCells->setValue(voxelIndex, static_cast(onsurf)); + } + } + } + } + } + + // Phase 2: Build NeighborList objects (operates on small feature-level data, no chunk iteration needed) + FloatVec3 spacing = imageGeom.getSpacing(); + + for(usize i = 1; i < totalFeatures; i++) + { + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Calculating Surface Areas || {:.2f}% Complete", CalculatePercentComplete(i, totalFeatures)); }); + + if(m_ShouldCancel) + { + return {}; + } + + std::map neighToCount; + auto numneighs = static_cast(neighborlist[i].size()); + + for(int32 j = 0; j < numneighs; j++) + { + neighToCount[neighborlist[i][j]]++; + } + + auto neighborIter = neighToCount.find(0); + neighToCount.erase(neighborIter); + neighborIter = neighToCount.find(-1); + if(neighborIter != neighToCount.end()) + { + neighToCount.erase(neighborIter); + } + neighborlist[i].resize(0); + neighborsurfacearealist[i].resize(0); + + for(const auto [neigh, number] : neighToCount) + { + float area = static_cast(number) * spacing[0] * spacing[1]; + neighborlist[i].push_back(neigh); + neighborsurfacearealist[i].push_back(area); + } + numNeighbors[i] = static_cast(neighborlist[i].size()); + + NeighborList::SharedVectorType sharedNeiLst(new std::vector); + sharedNeiLst->assign(neighborlist[i].begin(), neighborlist[i].end()); + neighborList.setList(static_cast(i), sharedNeiLst); + + NeighborList::SharedVectorType sharedSAL(new std::vector); + sharedSAL->assign(neighborsurfacearealist[i].begin(), neighborsurfacearealist[i].end()); + sharedSurfaceAreaList.setList(static_cast(i), sharedSAL); + } + + return {}; +} diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.hpp new file mode 100644 index 0000000000..9d9bb71916 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include "SimplnxCore/SimplnxCore_export.hpp" + +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ +struct ComputeFeatureNeighborsInputValues; + +class SIMPLNXCORE_EXPORT ComputeFeatureNeighborsScanline +{ +public: + ComputeFeatureNeighborsScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeFeatureNeighborsInputValues* inputValues); + ~ComputeFeatureNeighborsScanline() noexcept; + + ComputeFeatureNeighborsScanline(const ComputeFeatureNeighborsScanline&) = delete; + ComputeFeatureNeighborsScanline(ComputeFeatureNeighborsScanline&&) noexcept = delete; + ComputeFeatureNeighborsScanline& operator=(const ComputeFeatureNeighborsScanline&) = delete; + ComputeFeatureNeighborsScanline& operator=(ComputeFeatureNeighborsScanline&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const ComputeFeatureNeighborsInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp new file mode 100644 index 0000000000..a131959ff2 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp @@ -0,0 +1,155 @@ +#include "ComputeSurfaceAreaToVolumeDirect.hpp" + +#include "ComputeSurfaceAreaToVolume.hpp" + +#include "simplnx/Common/Constants.hpp" +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/DataGroup.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/Utilities/DataArrayUtilities.hpp" + +using namespace nx::core; + +// ----------------------------------------------------------------------------- +ComputeSurfaceAreaToVolumeDirect::ComputeSurfaceAreaToVolumeDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeSurfaceAreaToVolumeInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeSurfaceAreaToVolumeDirect::~ComputeSurfaceAreaToVolumeDirect() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> ComputeSurfaceAreaToVolumeDirect::operator()() +{ + // Input Cell Data + auto featureIdsArrayPtr = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath); + const auto& featureIdsStoreRef = featureIdsArrayPtr->getDataStoreRef(); + + // Input Feature Data + const auto& numCells = m_DataStructure.getDataAs(m_InputValues->NumCellsArrayPath)->getDataStoreRef(); + + // Output Feature Data + auto& surfaceAreaVolumeRatio = m_DataStructure.getDataAs(m_InputValues->SurfaceAreaVolumeRatioArrayName)->getDataStoreRef(); + + // Required Geometry + const auto& imageGeom = m_DataStructure.getDataRefAs(m_InputValues->InputImageGeometry); + + auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, m_InputValues->NumCellsArrayPath.getParent(), *featureIdsArrayPtr, false, m_MessageHandler); + if(validateNumFeatResult.invalid()) + { + return validateNumFeatResult; + } + auto numFeatures = static_cast(numCells.getNumberOfTuples()); + SizeVec3 dims = imageGeom.getDimensions(); + FloatVec3 spacing = imageGeom.getSpacing(); + + auto xPoints = static_cast(dims[0]); + auto yPoints = static_cast(dims[1]); + auto zPoints = static_cast(dims[2]); + + float32 voxelVol = spacing[0] * spacing[1] * spacing[2]; + + std::vector featureSurfaceArea(static_cast(numFeatures), 0.0f); + + int64 neighborOffset[6] = {0, 0, 0, 0, 0, 0}; + neighborOffset[0] = -xPoints * yPoints; // -Z + neighborOffset[1] = -xPoints; // -Y + neighborOffset[2] = -1; // -X + neighborOffset[3] = 1; // +X + neighborOffset[4] = xPoints; // +Y + neighborOffset[5] = xPoints * yPoints; // +Z + + for(int64 zIdx = 0; zIdx < zPoints; zIdx++) + { + m_MessageHandler(IFilter::Message::Type::Info, fmt::format("Computing Z Slice: '{}'", zIdx)); + + int64 zStride = zIdx * xPoints * yPoints; + for(int64 yIdx = 0; yIdx < yPoints; yIdx++) + { + int64 yStride = yIdx * xPoints; + for(int64 xIdx = 0; xIdx < xPoints; xIdx++) + { + float onSurface = 0.0f; + int32 currentFeatureId = featureIdsStoreRef[zStride + yStride + xIdx]; + if(currentFeatureId < 1) + { + continue; + } + + for(int32 neighborOffsetIndex = 0; neighborOffsetIndex < 6; neighborOffsetIndex++) + { + if(neighborOffsetIndex == 0 && zIdx == 0) + { + continue; + } + if(neighborOffsetIndex == 5 && zIdx == (zPoints - 1)) + { + continue; + } + if(neighborOffsetIndex == 1 && yIdx == 0) + { + continue; + } + if(neighborOffsetIndex == 4 && yIdx == (yPoints - 1)) + { + continue; + } + if(neighborOffsetIndex == 2 && xIdx == 0) + { + continue; + } + if(neighborOffsetIndex == 3 && xIdx == (xPoints - 1)) + { + continue; + } + + int64 neighborIndex = zStride + yStride + xIdx + neighborOffset[neighborOffsetIndex]; + + if(featureIdsStoreRef[neighborIndex] != currentFeatureId) + { + if(neighborOffsetIndex == 0 || neighborOffsetIndex == 5) // XY face shared + { + onSurface = onSurface + spacing[0] * spacing[1]; + } + if(neighborOffsetIndex == 1 || neighborOffsetIndex == 4) // YZ face shared + { + onSurface = onSurface + spacing[1] * spacing[2]; + } + if(neighborOffsetIndex == 2 || neighborOffsetIndex == 3) // XZ face shared + { + onSurface = onSurface + spacing[2] * spacing[0]; + } + } + } + int32 featureId = featureIdsStoreRef[zStride + yStride + xIdx]; + featureSurfaceArea[featureId] = featureSurfaceArea[featureId] + onSurface; + } + } + } + + const float32 thirdRootPi = std::pow(nx::core::Constants::k_PiF, 0.333333f); + for(usize i = 1; i < numFeatures; i++) + { + float featureVolume = voxelVol * numCells[i]; + surfaceAreaVolumeRatio[i] = featureSurfaceArea[i] / featureVolume; + } + + if(m_InputValues->CalculateSphericity) + { + m_MessageHandler(IFilter::Message::Type::Info, fmt::format("Computing Sphericity")); + + auto& sphericity = m_DataStructure.getDataAs(m_InputValues->SphericityArrayName)->getDataStoreRef(); + for(usize i = 1; i < static_cast(numFeatures); i++) + { + float featureVolume = voxelVol * numCells[i]; + sphericity[i] = (thirdRootPi * std::pow((6.0f * featureVolume), 0.66666f)) / featureSurfaceArea[i]; + } + } + + return {}; +} diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.hpp new file mode 100644 index 0000000000..1814f8f18f --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include "SimplnxCore/SimplnxCore_export.hpp" + +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ +struct ComputeSurfaceAreaToVolumeInputValues; + +class SIMPLNXCORE_EXPORT ComputeSurfaceAreaToVolumeDirect +{ +public: + ComputeSurfaceAreaToVolumeDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeSurfaceAreaToVolumeInputValues* inputValues); + ~ComputeSurfaceAreaToVolumeDirect() noexcept; + + ComputeSurfaceAreaToVolumeDirect(const ComputeSurfaceAreaToVolumeDirect&) = delete; + ComputeSurfaceAreaToVolumeDirect(ComputeSurfaceAreaToVolumeDirect&&) noexcept = delete; + ComputeSurfaceAreaToVolumeDirect& operator=(const ComputeSurfaceAreaToVolumeDirect&) = delete; + ComputeSurfaceAreaToVolumeDirect& operator=(ComputeSurfaceAreaToVolumeDirect&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const ComputeSurfaceAreaToVolumeInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp new file mode 100644 index 0000000000..2a2105ae87 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp @@ -0,0 +1,167 @@ +#include "ComputeSurfaceAreaToVolumeScanline.hpp" + +#include "ComputeSurfaceAreaToVolume.hpp" + +#include "simplnx/Common/Constants.hpp" +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/DataGroup.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/Utilities/DataArrayUtilities.hpp" + +using namespace nx::core; + +// ----------------------------------------------------------------------------- +ComputeSurfaceAreaToVolumeScanline::ComputeSurfaceAreaToVolumeScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeSurfaceAreaToVolumeInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeSurfaceAreaToVolumeScanline::~ComputeSurfaceAreaToVolumeScanline() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> ComputeSurfaceAreaToVolumeScanline::operator()() +{ + // Input Cell Data + auto featureIdsArrayPtr = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath); + auto& featureIdsStoreRef = featureIdsArrayPtr->getDataStoreRef(); + + // Input Feature Data + const auto& numCells = m_DataStructure.getDataAs(m_InputValues->NumCellsArrayPath)->getDataStoreRef(); + + // Output Feature Data + auto& surfaceAreaVolumeRatio = m_DataStructure.getDataAs(m_InputValues->SurfaceAreaVolumeRatioArrayName)->getDataStoreRef(); + + // Required Geometry + const auto& imageGeom = m_DataStructure.getDataRefAs(m_InputValues->InputImageGeometry); + + auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, m_InputValues->NumCellsArrayPath.getParent(), *featureIdsArrayPtr, false, m_MessageHandler); + if(validateNumFeatResult.invalid()) + { + return validateNumFeatResult; + } + auto numFeatures = static_cast(numCells.getNumberOfTuples()); + SizeVec3 dims = imageGeom.getDimensions(); + FloatVec3 spacing = imageGeom.getSpacing(); + + auto xPoints = static_cast(dims[0]); + auto yPoints = static_cast(dims[1]); + auto zPoints = static_cast(dims[2]); + + float32 voxelVol = spacing[0] * spacing[1] * spacing[2]; + + std::vector featureSurfaceArea(static_cast(numFeatures), 0.0f); + + int64 neighborOffset[6] = {0, 0, 0, 0, 0, 0}; + neighborOffset[0] = -xPoints * yPoints; // -Z + neighborOffset[1] = -xPoints; // -Y + neighborOffset[2] = -1; // -X + neighborOffset[3] = 1; // +X + neighborOffset[4] = xPoints; // +Y + neighborOffset[5] = xPoints * yPoints; // +Z + + // Chunk-sequential iteration + const uint64 numChunks = featureIdsStoreRef.getNumberOfChunks(); + + for(uint64 chunkIdx = 0; chunkIdx < numChunks; chunkIdx++) + { + featureIdsStoreRef.loadChunk(chunkIdx); + + const auto chunkLowerBounds = featureIdsStoreRef.getChunkLowerBounds(chunkIdx); + const auto chunkUpperBounds = featureIdsStoreRef.getChunkUpperBounds(chunkIdx); + + m_MessageHandler(IFilter::Message::Type::Info, fmt::format("Processing Chunk {}/{}", chunkIdx + 1, numChunks)); + + for(usize zIdx = chunkLowerBounds[0]; zIdx <= chunkUpperBounds[0]; zIdx++) + { + const int64 zStride = static_cast(zIdx) * xPoints * yPoints; + for(usize yIdx = chunkLowerBounds[1]; yIdx <= chunkUpperBounds[1]; yIdx++) + { + const int64 yStride = static_cast(yIdx) * xPoints; + for(usize xIdx = chunkLowerBounds[2]; xIdx <= chunkUpperBounds[2]; xIdx++) + { + const int64 ix = static_cast(xIdx); + float onSurface = 0.0f; + int32 currentFeatureId = featureIdsStoreRef[zStride + yStride + ix]; + if(currentFeatureId < 1) + { + continue; + } + + for(int32 neighborOffsetIndex = 0; neighborOffsetIndex < 6; neighborOffsetIndex++) + { + if(neighborOffsetIndex == 0 && static_cast(zIdx) == 0) + { + continue; + } + if(neighborOffsetIndex == 5 && static_cast(zIdx) == (zPoints - 1)) + { + continue; + } + if(neighborOffsetIndex == 1 && static_cast(yIdx) == 0) + { + continue; + } + if(neighborOffsetIndex == 4 && static_cast(yIdx) == (yPoints - 1)) + { + continue; + } + if(neighborOffsetIndex == 2 && ix == 0) + { + continue; + } + if(neighborOffsetIndex == 3 && ix == (xPoints - 1)) + { + continue; + } + + int64 neighborIndex = zStride + yStride + ix + neighborOffset[neighborOffsetIndex]; + + if(featureIdsStoreRef[neighborIndex] != currentFeatureId) + { + if(neighborOffsetIndex == 0 || neighborOffsetIndex == 5) + { + onSurface = onSurface + spacing[0] * spacing[1]; + } + if(neighborOffsetIndex == 1 || neighborOffsetIndex == 4) + { + onSurface = onSurface + spacing[1] * spacing[2]; + } + if(neighborOffsetIndex == 2 || neighborOffsetIndex == 3) + { + onSurface = onSurface + spacing[2] * spacing[0]; + } + } + } + int32 featureId = featureIdsStoreRef[zStride + yStride + ix]; + featureSurfaceArea[featureId] = featureSurfaceArea[featureId] + onSurface; + } + } + } + } + + const float32 thirdRootPi = std::pow(nx::core::Constants::k_PiF, 0.333333f); + for(usize i = 1; i < numFeatures; i++) + { + float featureVolume = voxelVol * numCells[i]; + surfaceAreaVolumeRatio[i] = featureSurfaceArea[i] / featureVolume; + } + + if(m_InputValues->CalculateSphericity) + { + m_MessageHandler(IFilter::Message::Type::Info, fmt::format("Computing Sphericity")); + + auto& sphericity = m_DataStructure.getDataAs(m_InputValues->SphericityArrayName)->getDataStoreRef(); + for(usize i = 1; i < static_cast(numFeatures); i++) + { + float featureVolume = voxelVol * numCells[i]; + sphericity[i] = (thirdRootPi * std::pow((6.0f * featureVolume), 0.66666f)) / featureSurfaceArea[i]; + } + } + + return {}; +} diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp new file mode 100644 index 0000000000..404bd7cf20 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include "SimplnxCore/SimplnxCore_export.hpp" + +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ +struct ComputeSurfaceAreaToVolumeInputValues; + +class SIMPLNXCORE_EXPORT ComputeSurfaceAreaToVolumeScanline +{ +public: + ComputeSurfaceAreaToVolumeScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeSurfaceAreaToVolumeInputValues* inputValues); + ~ComputeSurfaceAreaToVolumeScanline() noexcept; + + ComputeSurfaceAreaToVolumeScanline(const ComputeSurfaceAreaToVolumeScanline&) = delete; + ComputeSurfaceAreaToVolumeScanline(ComputeSurfaceAreaToVolumeScanline&&) noexcept = delete; + ComputeSurfaceAreaToVolumeScanline& operator=(const ComputeSurfaceAreaToVolumeScanline&) = delete; + ComputeSurfaceAreaToVolumeScanline& operator=(ComputeSurfaceAreaToVolumeScanline&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const ComputeSurfaceAreaToVolumeInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.cpp new file mode 100644 index 0000000000..ba5def33eb --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.cpp @@ -0,0 +1,230 @@ +#include "ComputeSurfaceFeaturesDirect.hpp" + +#include "ComputeSurfaceFeatures.hpp" + +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/Utilities/DataArrayUtilities.hpp" + +using namespace nx::core; + +namespace +{ +bool IsPointASurfaceFeature(const Point2D& point, usize xPoints, usize yPoints, bool markFeature0Neighbors, const Int32AbstractDataStore& featureIds) +{ + const usize yStride = point.getY() * xPoints; + + if(point.getX() <= 0 || point.getX() >= xPoints - 1) + { + return true; + } + if(point.getY() <= 0 || point.getY() >= yPoints - 1) + { + return true; + } + + if(markFeature0Neighbors) + { + if(featureIds[yStride + point.getX() - 1] == 0) + { + return true; + } + if(featureIds[yStride + point.getX() + 1] == 0) + { + return true; + } + if(featureIds[yStride + point.getX() - xPoints] == 0) + { + return true; + } + if(featureIds[yStride + point.getX() + xPoints] == 0) + { + return true; + } + } + + return false; +} + +bool IsPointASurfaceFeature(const Point3D& point, usize xPoints, usize yPoints, usize zPoints, bool markFeature0Neighbors, const Int32AbstractDataStore& featureIds) +{ + usize yStride = point.getY() * xPoints; + usize zStride = point.getZ() * xPoints * yPoints; + + if(point.getX() <= 0 || point.getX() >= xPoints - 1) + { + return true; + } + if(point.getY() <= 0 || point.getY() >= yPoints - 1) + { + return true; + } + if(point.getZ() <= 0 || point.getZ() >= zPoints - 1) + { + return true; + } + + if(markFeature0Neighbors) + { + if(featureIds[zStride + yStride + point.getX() - 1] == 0) + { + return true; + } + if(featureIds[zStride + yStride + point.getX() + 1] == 0) + { + return true; + } + if(featureIds[zStride + yStride + point.getX() - xPoints] == 0) + { + return true; + } + if(featureIds[zStride + yStride + point.getX() + xPoints] == 0) + { + return true; + } + if(featureIds[zStride + yStride + point.getX() - (xPoints * yPoints)] == 0) + { + return true; + } + if(featureIds[zStride + yStride + point.getX() + (xPoints * yPoints)] == 0) + { + return true; + } + } + + return false; +} + +void findSurfaceFeatures3D(DataStructure& dataStructure, const DataPath& featureGeometryPathValue, const DataPath& featureIdsArrayPathValue, const DataPath& surfaceFeaturesArrayPathValue, + bool markFeature0Neighbors, const std::atomic_bool& shouldCancel) +{ + const auto& featureGeometry = dataStructure.getDataRefAs(featureGeometryPathValue); + const auto& featureIds = dataStructure.getDataAs(featureIdsArrayPathValue)->getDataStoreRef(); + auto& surfaceFeatures = dataStructure.getDataAs(surfaceFeaturesArrayPathValue)->getDataStoreRef(); + + const usize xPoints = featureGeometry.getNumXCells(); + const usize yPoints = featureGeometry.getNumYCells(); + const usize zPoints = featureGeometry.getNumZCells(); + + for(usize z = 0; z < zPoints; z++) + { + const usize zStride = z * xPoints * yPoints; + for(usize y = 0; y < yPoints; y++) + { + const usize yStride = y * xPoints; + for(usize x = 0; x < xPoints; x++) + { + if(shouldCancel) + { + return; + } + + const int32 gNum = featureIds[zStride + yStride + x]; + if(gNum != 0 && !surfaceFeatures[gNum]) + { + if(IsPointASurfaceFeature(Point3D{x, y, z}, xPoints, yPoints, zPoints, markFeature0Neighbors, featureIds)) + { + surfaceFeatures[gNum] = 1; + } + } + } + } + } +} + +void findSurfaceFeatures2D(DataStructure& dataStructure, const DataPath& featureGeometryPathValue, const DataPath& featureIdsArrayPathValue, const DataPath& surfaceFeaturesArrayPathValue, + bool markFeature0Neighbors, const std::atomic_bool& shouldCancel) +{ + const auto& featureGeometry = dataStructure.getDataRefAs(featureGeometryPathValue); + const auto& featureIds = dataStructure.getDataAs(featureIdsArrayPathValue)->getDataStoreRef(); + auto& surfaceFeatures = dataStructure.getDataAs(surfaceFeaturesArrayPathValue)->getDataStoreRef(); + + usize xPoints = 0; + usize yPoints = 0; + + if(featureGeometry.getNumXCells() == 1) + { + xPoints = featureGeometry.getNumYCells(); + yPoints = featureGeometry.getNumZCells(); + } + if(featureGeometry.getNumYCells() == 1) + { + xPoints = featureGeometry.getNumXCells(); + yPoints = featureGeometry.getNumZCells(); + } + if(featureGeometry.getNumZCells() == 1) + { + xPoints = featureGeometry.getNumXCells(); + yPoints = featureGeometry.getNumYCells(); + } + + for(usize y = 0; y < yPoints; y++) + { + const usize yStride = y * xPoints; + + for(usize x = 0; x < xPoints; x++) + { + if(shouldCancel) + { + return; + } + + const int32 gNum = featureIds[yStride + x]; + if(gNum != 0 && surfaceFeatures[gNum] == 0) + { + if(IsPointASurfaceFeature(Point2D{x, y}, xPoints, yPoints, markFeature0Neighbors, featureIds)) + { + surfaceFeatures[gNum] = 1; + } + } + } + } +} +} // namespace + +// ----------------------------------------------------------------------------- +ComputeSurfaceFeaturesDirect::ComputeSurfaceFeaturesDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeSurfaceFeaturesInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeSurfaceFeaturesDirect::~ComputeSurfaceFeaturesDirect() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> ComputeSurfaceFeaturesDirect::operator()() +{ + const auto pMarkFeature0NeighborsValue = m_InputValues->MarkFeature0Neighbors; + const auto pFeatureGeometryPathValue = m_InputValues->InputImageGeometryPath; + const auto pFeatureIdsArrayPathValue = m_InputValues->FeatureIdsPath; + const auto pFeaturesAttributeMatrixPathValue = m_InputValues->FeatureAttributeMatrixPath; + const auto pSurfaceFeaturesArrayPathValue = pFeaturesAttributeMatrixPathValue.createChildPath(m_InputValues->SurfaceFeaturesArrayName); + + const auto& featureIdsArray = m_DataStructure.getDataRefAs(pFeatureIdsArrayPathValue); + + auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, pFeaturesAttributeMatrixPathValue, featureIdsArray, false, m_MessageHandler); + if(validateNumFeatResult.invalid()) + { + return validateNumFeatResult; + } + + const auto& featureGeometry = m_DataStructure.getDataRefAs(pFeatureGeometryPathValue); + if(const usize geometryDimensionality = featureGeometry.getDimensionality(); geometryDimensionality == 3) + { + findSurfaceFeatures3D(m_DataStructure, pFeatureGeometryPathValue, pFeatureIdsArrayPathValue, pSurfaceFeaturesArrayPathValue, pMarkFeature0NeighborsValue, m_ShouldCancel); + } + else if(geometryDimensionality == 2) + { + findSurfaceFeatures2D(m_DataStructure, pFeatureGeometryPathValue, pFeatureIdsArrayPathValue, pSurfaceFeaturesArrayPathValue, pMarkFeature0NeighborsValue, m_ShouldCancel); + } + else + { + MakeErrorResult(-1000, fmt::format("Image Geometry at path '{}' must be either 3D or 2D", pFeatureGeometryPathValue.toString())); + } + + return {}; +} diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp new file mode 100644 index 0000000000..d7810f2502 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include "SimplnxCore/SimplnxCore_export.hpp" + +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ +struct ComputeSurfaceFeaturesInputValues; + +class SIMPLNXCORE_EXPORT ComputeSurfaceFeaturesDirect +{ +public: + ComputeSurfaceFeaturesDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeSurfaceFeaturesInputValues* inputValues); + ~ComputeSurfaceFeaturesDirect() noexcept; + + ComputeSurfaceFeaturesDirect(const ComputeSurfaceFeaturesDirect&) = delete; + ComputeSurfaceFeaturesDirect(ComputeSurfaceFeaturesDirect&&) noexcept = delete; + ComputeSurfaceFeaturesDirect& operator=(const ComputeSurfaceFeaturesDirect&) = delete; + ComputeSurfaceFeaturesDirect& operator=(ComputeSurfaceFeaturesDirect&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const ComputeSurfaceFeaturesInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.cpp new file mode 100644 index 0000000000..306c955c7a --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.cpp @@ -0,0 +1,226 @@ +#include "ComputeSurfaceFeaturesScanline.hpp" + +#include "ComputeSurfaceFeatures.hpp" + +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/Utilities/DataArrayUtilities.hpp" + +using namespace nx::core; + +namespace +{ +bool IsPointASurfaceFeature2D(usize remappedX, usize remappedY, usize xPoints, usize yPoints, bool markFeature0Neighbors, const Int32AbstractDataStore& featureIds) +{ + const usize yStride = remappedY * xPoints; + + if(remappedX <= 0 || remappedX >= xPoints - 1) + { + return true; + } + if(remappedY <= 0 || remappedY >= yPoints - 1) + { + return true; + } + + if(markFeature0Neighbors) + { + if(featureIds[yStride + remappedX - 1] == 0) + { + return true; + } + if(featureIds[yStride + remappedX + 1] == 0) + { + return true; + } + if(featureIds[yStride + remappedX - xPoints] == 0) + { + return true; + } + if(featureIds[yStride + remappedX + xPoints] == 0) + { + return true; + } + } + + return false; +} + +bool IsPointASurfaceFeature3D(usize x, usize y, usize z, usize xPoints, usize yPoints, usize zPoints, bool markFeature0Neighbors, const Int32AbstractDataStore& featureIds) +{ + const usize yStride = y * xPoints; + const usize zStride = z * xPoints * yPoints; + + if(x <= 0 || x >= xPoints - 1) + { + return true; + } + if(y <= 0 || y >= yPoints - 1) + { + return true; + } + if(z <= 0 || z >= zPoints - 1) + { + return true; + } + + if(markFeature0Neighbors) + { + if(featureIds[zStride + yStride + x - 1] == 0) + { + return true; + } + if(featureIds[zStride + yStride + x + 1] == 0) + { + return true; + } + if(featureIds[zStride + yStride + x - xPoints] == 0) + { + return true; + } + if(featureIds[zStride + yStride + x + xPoints] == 0) + { + return true; + } + if(featureIds[zStride + yStride + x - (xPoints * yPoints)] == 0) + { + return true; + } + if(featureIds[zStride + yStride + x + (xPoints * yPoints)] == 0) + { + return true; + } + } + + return false; +} +} // namespace + +// ----------------------------------------------------------------------------- +ComputeSurfaceFeaturesScanline::ComputeSurfaceFeaturesScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeSurfaceFeaturesInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeSurfaceFeaturesScanline::~ComputeSurfaceFeaturesScanline() noexcept = default; + +// ----------------------------------------------------------------------------- +Result<> ComputeSurfaceFeaturesScanline::operator()() +{ + const auto pMarkFeature0NeighborsValue = m_InputValues->MarkFeature0Neighbors; + const auto pFeatureGeometryPathValue = m_InputValues->InputImageGeometryPath; + const auto pFeatureIdsArrayPathValue = m_InputValues->FeatureIdsPath; + const auto pFeaturesAttributeMatrixPathValue = m_InputValues->FeatureAttributeMatrixPath; + const auto pSurfaceFeaturesArrayPathValue = pFeaturesAttributeMatrixPathValue.createChildPath(m_InputValues->SurfaceFeaturesArrayName); + + const auto& featureIdsArray = m_DataStructure.getDataRefAs(pFeatureIdsArrayPathValue); + + auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, pFeaturesAttributeMatrixPathValue, featureIdsArray, false, m_MessageHandler); + if(validateNumFeatResult.invalid()) + { + return validateNumFeatResult; + } + + const auto& featureGeometry = m_DataStructure.getDataRefAs(pFeatureGeometryPathValue); + auto& featureIds = m_DataStructure.getDataAs(pFeatureIdsArrayPathValue)->getDataStoreRef(); + auto& surfaceFeatures = m_DataStructure.getDataAs(pSurfaceFeaturesArrayPathValue)->getDataStoreRef(); + + const usize xPoints = featureGeometry.getNumXCells(); + const usize yPoints = featureGeometry.getNumYCells(); + const usize zPoints = featureGeometry.getNumZCells(); + const usize geometryDimensionality = featureGeometry.getDimensionality(); + + // For 2D geometries, compute the remapped dimensions + usize remappedXPoints = 0; + usize remappedYPoints = 0; + if(geometryDimensionality == 2) + { + if(xPoints == 1) + { + remappedXPoints = yPoints; + remappedYPoints = zPoints; + } + else if(yPoints == 1) + { + remappedXPoints = xPoints; + remappedYPoints = zPoints; + } + else // zPoints == 1 + { + remappedXPoints = xPoints; + remappedYPoints = yPoints; + } + } + + // Chunk-sequential iteration always uses 3D chunk bounds + const uint64 numChunks = featureIds.getNumberOfChunks(); + + for(uint64 chunkIdx = 0; chunkIdx < numChunks; chunkIdx++) + { + featureIds.loadChunk(chunkIdx); + + const auto chunkLowerBounds = featureIds.getChunkLowerBounds(chunkIdx); + const auto chunkUpperBounds = featureIds.getChunkUpperBounds(chunkIdx); + + for(usize z = chunkLowerBounds[0]; z <= chunkUpperBounds[0]; z++) + { + const usize zStride = z * xPoints * yPoints; + for(usize y = chunkLowerBounds[1]; y <= chunkUpperBounds[1]; y++) + { + const usize yStride = y * xPoints; + for(usize x = chunkLowerBounds[2]; x <= chunkUpperBounds[2]; x++) + { + if(m_ShouldCancel) + { + return {}; + } + + const int32 gNum = featureIds[zStride + yStride + x]; + if(gNum != 0 && !surfaceFeatures[gNum]) + { + if(geometryDimensionality == 3) + { + if(IsPointASurfaceFeature3D(x, y, z, xPoints, yPoints, zPoints, pMarkFeature0NeighborsValue, featureIds)) + { + surfaceFeatures[gNum] = 1; + } + } + else if(geometryDimensionality == 2) + { + // Remap native 3D coordinates to 2D based on degenerate dimension + usize remappedX = 0; + usize remappedY = 0; + if(xPoints == 1) + { + remappedX = y; + remappedY = z; + } + else if(yPoints == 1) + { + remappedX = x; + remappedY = z; + } + else // zPoints == 1 + { + remappedX = x; + remappedY = y; + } + + if(IsPointASurfaceFeature2D(remappedX, remappedY, remappedXPoints, remappedYPoints, pMarkFeature0NeighborsValue, featureIds)) + { + surfaceFeatures[gNum] = 1; + } + } + } + } + } + } + } + + return {}; +} diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp new file mode 100644 index 0000000000..a5a6c23bf2 --- /dev/null +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include "SimplnxCore/SimplnxCore_export.hpp" + +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/Filter/IFilter.hpp" + +namespace nx::core +{ +struct ComputeSurfaceFeaturesInputValues; + +class SIMPLNXCORE_EXPORT ComputeSurfaceFeaturesScanline +{ +public: + ComputeSurfaceFeaturesScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + const ComputeSurfaceFeaturesInputValues* inputValues); + ~ComputeSurfaceFeaturesScanline() noexcept; + + ComputeSurfaceFeaturesScanline(const ComputeSurfaceFeaturesScanline&) = delete; + ComputeSurfaceFeaturesScanline(ComputeSurfaceFeaturesScanline&&) noexcept = delete; + ComputeSurfaceFeaturesScanline& operator=(const ComputeSurfaceFeaturesScanline&) = delete; + ComputeSurfaceFeaturesScanline& operator=(ComputeSurfaceFeaturesScanline&&) noexcept = delete; + + Result<> operator()(); + +private: + DataStructure& m_DataStructure; + const ComputeSurfaceFeaturesInputValues* m_InputValues = nullptr; + const std::atomic_bool& m_ShouldCancel; + const IFilter::MessageHandler& m_MessageHandler; +}; + +} // namespace nx::core diff --git a/src/Plugins/SimplnxCore/test/ComputeBoundaryCellsTest.cpp b/src/Plugins/SimplnxCore/test/ComputeBoundaryCellsTest.cpp index f5a9f4a393..f8a976988f 100644 --- a/src/Plugins/SimplnxCore/test/ComputeBoundaryCellsTest.cpp +++ b/src/Plugins/SimplnxCore/test/ComputeBoundaryCellsTest.cpp @@ -1,6 +1,9 @@ #include +#include "simplnx/DataStructure/AttributeMatrix.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" #include "simplnx/UnitTest/UnitTestCommon.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" #include "SimplnxCore/Filters/ComputeBoundaryCellsFilter.hpp" #include "SimplnxCore/SimplnxCore_test_dirs.hpp" @@ -20,6 +23,9 @@ const DataPath k_ComputedBoundaryCellsPath({k_ExemplarDataContainer, Constants:: TEST_CASE("SimplnxCore::ComputeBoundaryCellsFilter: Valid filter execution", "[ComputeBoundaryCellsFilter]") { UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 65536, true); const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_TestFilesDir, "6_6_find_boundary_cells.tar.gz", "6_6_FindBoundaryCellsExemplar.dream3d"); @@ -88,3 +94,64 @@ TEST_CASE("SimplnxCore::ComputeBoundaryCellsFilter: Invalid filter execution", " UnitTest::CheckArraysInheritTupleDims(dataStructure); } + +TEST_CASE("SimplnxCore::ComputeBoundaryCellsFilter: Benchmark 200x200x200", "[ComputeBoundaryCellsFilter][Benchmark]") +{ + UnitTest::LoadPlugins(); + + constexpr usize kDimX = 200; + constexpr usize kDimY = 200; + constexpr usize kDimZ = 200; + constexpr usize kTotalVoxels = kDimX * kDimY * kDimZ; + const ShapeType cellTupleShape = {kDimZ, kDimY, kDimX}; + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 160000, true); + + const auto benchmarkFile = fs::path(fmt::format("{}/compute_boundary_cells_benchmark.dream3d", unit_test::k_BinaryTestOutputDir)); + + const std::string k_BenchGeomName = "ImageGeom"; + const DataPath k_BenchGeomPath({k_BenchGeomName}); + const DataPath k_BenchFeatureIdsPath = k_BenchGeomPath.createChildPath(Constants::k_CellData).createChildPath(Constants::k_FeatureIds); + + { + DataStructure buildDS; + auto* imageGeom = ImageGeom::Create(buildDS, k_BenchGeomName); + imageGeom->setDimensions({kDimX, kDimY, kDimZ}); + imageGeom->setSpacing({1.0f, 1.0f, 1.0f}); + imageGeom->setOrigin({0.0f, 0.0f, 0.0f}); + + auto* cellAM = AttributeMatrix::Create(buildDS, Constants::k_CellData, cellTupleShape, imageGeom->getId()); + imageGeom->setCellData(*cellAM); + + auto* featureIds = UnitTest::CreateTestDataArray(buildDS, Constants::k_FeatureIds, cellTupleShape, {1}, cellAM->getId()); + auto& fidsStore = featureIds->getDataStoreRef(); + for(usize i = 0; i < kTotalVoxels; i++) + { + const usize ix = i % kDimX; + const usize iy = (i / kDimX) % kDimY; + const usize iz = i / (kDimX * kDimY); + const int32 octant = static_cast((iz >= kDimZ / 2 ? 4 : 0) + (iy >= kDimY / 2 ? 2 : 0) + (ix >= kDimX / 2 ? 1 : 0)) + 1; + fidsStore.setValue(i, ((i * 7 + 13) % 100) < 5 ? 0 : octant); + } + + UnitTest::WriteTestDataStructure(buildDS, benchmarkFile); + } + + DataStructure dataStructure = UnitTest::LoadDataStructure(benchmarkFile); + + { + ComputeBoundaryCellsFilter filter; + Arguments args; + args.insertOrAssign(ComputeBoundaryCellsFilter::k_IgnoreFeatureZero_Key, std::make_any(true)); + args.insertOrAssign(ComputeBoundaryCellsFilter::k_IncludeVolumeBoundary_Key, std::make_any(true)); + args.insertOrAssign(ComputeBoundaryCellsFilter::k_GeometryPath_Key, std::make_any(k_BenchGeomPath)); + args.insertOrAssign(ComputeBoundaryCellsFilter::k_FeatureIdsArrayPath_Key, std::make_any(k_BenchFeatureIdsPath)); + args.insertOrAssign(ComputeBoundaryCellsFilter::k_BoundaryCellsArrayName_Key, std::make_any("BoundaryCells")); + + auto preflightResult = filter.preflight(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions); + auto executeResult = filter.execute(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result); + } + + fs::remove(benchmarkFile); +} diff --git a/src/Plugins/SimplnxCore/test/ComputeFeatureNeighborsTest.cpp b/src/Plugins/SimplnxCore/test/ComputeFeatureNeighborsTest.cpp index df5f685ce1..13c1eb5866 100644 --- a/src/Plugins/SimplnxCore/test/ComputeFeatureNeighborsTest.cpp +++ b/src/Plugins/SimplnxCore/test/ComputeFeatureNeighborsTest.cpp @@ -1,7 +1,10 @@ #include "SimplnxCore/Filters/ComputeFeatureNeighborsFilter.hpp" #include "SimplnxCore/SimplnxCore_test_dirs.hpp" +#include "simplnx/DataStructure/AttributeMatrix.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" #include "simplnx/UnitTest/UnitTestCommon.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" #include @@ -13,6 +16,9 @@ using namespace nx::core::UnitTest; TEST_CASE("SimplnxCore::ComputeFeatureNeighborsFilter", "[SimplnxCore][ComputeFeatureNeighborsFilter]") { UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 25600, true); const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_TestFilesDir, "6_6_stats_test_v2.tar.gz", "6_6_stats_test_v2.dream3d"); // Read the Small IN100 Data set @@ -79,3 +85,72 @@ TEST_CASE("SimplnxCore::ComputeFeatureNeighborsFilter", "[SimplnxCore][ComputeFe UnitTest::CheckArraysInheritTupleDims(dataStructure); } + +TEST_CASE("SimplnxCore::ComputeFeatureNeighborsFilter: Benchmark 200x200x200", "[SimplnxCore][ComputeFeatureNeighborsFilter][Benchmark]") +{ + UnitTest::LoadPlugins(); + + constexpr usize kDimX = 200; + constexpr usize kDimY = 200; + constexpr usize kDimZ = 200; + constexpr usize kTotalVoxels = kDimX * kDimY * kDimZ; + const ShapeType cellTupleShape = {kDimZ, kDimY, kDimX}; + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 160000, true); + + const auto benchmarkFile = fs::path(fmt::format("{}/compute_feature_neighbors_benchmark.dream3d", unit_test::k_BinaryTestOutputDir)); + + const std::string k_BenchGeomName = "ImageGeom"; + const DataPath k_BenchGeomPath({k_BenchGeomName}); + const DataPath k_BenchFeatureIdsPath = k_BenchGeomPath.createChildPath(Constants::k_CellData).createChildPath(Constants::k_FeatureIds); + const DataPath k_BenchCellFeatureAMPath = k_BenchGeomPath.createChildPath(Constants::k_CellFeatureData); + + { + DataStructure buildDS; + auto* imageGeom = ImageGeom::Create(buildDS, k_BenchGeomName); + imageGeom->setDimensions({kDimX, kDimY, kDimZ}); + imageGeom->setSpacing({1.0f, 1.0f, 1.0f}); + imageGeom->setOrigin({0.0f, 0.0f, 0.0f}); + + auto* cellAM = AttributeMatrix::Create(buildDS, Constants::k_CellData, cellTupleShape, imageGeom->getId()); + imageGeom->setCellData(*cellAM); + + constexpr int32 kNumFeatures = 8; + auto* featureIds = UnitTest::CreateTestDataArray(buildDS, Constants::k_FeatureIds, cellTupleShape, {1}, cellAM->getId()); + auto& fidsStore = featureIds->getDataStoreRef(); + for(usize i = 0; i < kTotalVoxels; i++) + { + const usize ix = i % kDimX; + const usize iy = (i / kDimX) % kDimY; + const usize iz = i / (kDimX * kDimY); + const int32 octant = static_cast((iz >= kDimZ / 2 ? 4 : 0) + (iy >= kDimY / 2 ? 2 : 0) + (ix >= kDimX / 2 ? 1 : 0)) + 1; + fidsStore.setValue(i, ((i * 7 + 13) % 100) < 5 ? 0 : octant); + } + + AttributeMatrix::Create(buildDS, Constants::k_CellFeatureData, {static_cast(kNumFeatures + 1)}, imageGeom->getId()); + UnitTest::WriteTestDataStructure(buildDS, benchmarkFile); + } + + DataStructure dataStructure = UnitTest::LoadDataStructure(benchmarkFile); + + { + ComputeFeatureNeighborsFilter filter; + Arguments args; + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_SelectedImageGeometryPath_Key, std::make_any(k_BenchGeomPath)); + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_FeatureIdsPath_Key, std::make_any(k_BenchFeatureIdsPath)); + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_CellFeaturesPath_Key, std::make_any(k_BenchCellFeatureAMPath)); + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_StoreBoundary_Key, std::make_any(true)); + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_BoundaryCellsName_Key, std::make_any("BoundaryCells")); + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_StoreSurface_Key, std::make_any(true)); + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_SurfaceFeaturesName_Key, std::make_any("SurfaceFeatures")); + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_NumNeighborsName_Key, std::make_any("NumNeighbors")); + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_NeighborListName_Key, std::make_any("NeighborList")); + args.insertOrAssign(ComputeFeatureNeighborsFilter::k_SharedSurfaceAreaName_Key, std::make_any("SharedSurfaceAreaList")); + + auto preflightResult = filter.preflight(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions); + auto executeResult = filter.execute(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result); + } + + fs::remove(benchmarkFile); +} diff --git a/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTest.cpp b/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTest.cpp index 110bcb34cc..d02bf12efb 100644 --- a/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTest.cpp +++ b/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTest.cpp @@ -1,9 +1,12 @@ #include "SimplnxCore/Filters/ComputeSurfaceAreaToVolumeFilter.hpp" #include "SimplnxCore/SimplnxCore_test_dirs.hpp" +#include "simplnx/DataStructure/AttributeMatrix.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" #include "simplnx/Parameters/ArrayCreationParameter.hpp" #include "simplnx/Parameters/BoolParameter.hpp" #include "simplnx/UnitTest/UnitTestCommon.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" #include @@ -26,6 +29,9 @@ const std::string k_SphericityArrayNameNX("SphericityNX"); TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolume", "[SimplnxCore][ComputeSurfaceAreaToVolume]") { UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 25600, true); const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_TestFilesDir, "6_6_stats_test_v2.tar.gz", "6_6_stats_test_v2.dream3d"); @@ -77,3 +83,85 @@ TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolume", "[SimplnxCore][ComputeSurfa UnitTest::CheckArraysInheritTupleDims(dataStructure); } + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolume: Benchmark 200x200x200", "[SimplnxCore][ComputeSurfaceAreaToVolume][Benchmark]") +{ + UnitTest::LoadPlugins(); + + constexpr usize kDimX = 200; + constexpr usize kDimY = 200; + constexpr usize kDimZ = 200; + constexpr usize kTotalVoxels = kDimX * kDimY * kDimZ; + const ShapeType cellTupleShape = {kDimZ, kDimY, kDimX}; + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 160000, true); + + const auto benchmarkFile = fs::path(fmt::format("{}/compute_surface_area_to_volume_benchmark.dream3d", unit_test::k_BinaryTestOutputDir)); + + const std::string k_BenchGeomName = "ImageGeom"; + const DataPath k_BenchGeomPath({k_BenchGeomName}); + const DataPath k_BenchFeatureIdsPath = k_BenchGeomPath.createChildPath(Constants::k_CellData).createChildPath(Constants::k_FeatureIds); + const DataPath k_BenchCellFeatureAMPath = k_BenchGeomPath.createChildPath(Constants::k_CellFeatureData); + const DataPath k_BenchNumCellsPath = k_BenchCellFeatureAMPath.createChildPath(Constants::k_NumElements); + + constexpr int32 kNumFeatures = 8; + + { + DataStructure buildDS; + auto* imageGeom = ImageGeom::Create(buildDS, k_BenchGeomName); + imageGeom->setDimensions({kDimX, kDimY, kDimZ}); + imageGeom->setSpacing({1.0f, 1.0f, 1.0f}); + imageGeom->setOrigin({0.0f, 0.0f, 0.0f}); + + auto* cellAM = AttributeMatrix::Create(buildDS, Constants::k_CellData, cellTupleShape, imageGeom->getId()); + imageGeom->setCellData(*cellAM); + + auto* featureIds = UnitTest::CreateTestDataArray(buildDS, Constants::k_FeatureIds, cellTupleShape, {1}, cellAM->getId()); + auto& fidsStore = featureIds->getDataStoreRef(); + + // Count voxels per feature to build NumCells + std::vector featureCounts(kNumFeatures + 1, 0); + for(usize i = 0; i < kTotalVoxels; i++) + { + const usize ix = i % kDimX; + const usize iy = (i / kDimX) % kDimY; + const usize iz = i / (kDimX * kDimY); + const int32 octant = static_cast((iz >= kDimZ / 2 ? 4 : 0) + (iy >= kDimY / 2 ? 2 : 0) + (ix >= kDimX / 2 ? 1 : 0)) + 1; + const int32 fid = ((i * 7 + 13) % 100) < 5 ? 0 : octant; + fidsStore.setValue(i, fid); + if(fid > 0) + { + featureCounts[fid]++; + } + } + + auto* cellFeatureAM = AttributeMatrix::Create(buildDS, Constants::k_CellFeatureData, {static_cast(kNumFeatures + 1)}, imageGeom->getId()); + auto* numCells = UnitTest::CreateTestDataArray(buildDS, Constants::k_NumElements, {static_cast(kNumFeatures + 1)}, {1}, cellFeatureAM->getId()); + auto& numCellsStore = numCells->getDataStoreRef(); + for(int32 f = 0; f <= kNumFeatures; f++) + { + numCellsStore.setValue(static_cast(f), featureCounts[f]); + } + + UnitTest::WriteTestDataStructure(buildDS, benchmarkFile); + } + + DataStructure dataStructure = UnitTest::LoadDataStructure(benchmarkFile); + + { + ComputeSurfaceAreaToVolumeFilter filter; + Arguments args; + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SelectedImageGeometryPath_Key, std::make_any(k_BenchGeomPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CellFeatureIdsArrayPath_Key, std::make_any(k_BenchFeatureIdsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_NumCellsArrayPath_Key, std::make_any(k_BenchNumCellsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CalculateSphericity_Key, std::make_any(true)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SurfaceAreaVolumeRatioArrayName_Key, std::make_any("SurfaceAreaVolumeRatio")); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SphericityArrayName_Key, std::make_any("Sphericity")); + + auto preflightResult = filter.preflight(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions); + auto executeResult = filter.execute(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result); + } + + fs::remove(benchmarkFile); +} diff --git a/src/Plugins/SimplnxCore/test/ComputeSurfaceFeaturesTest.cpp b/src/Plugins/SimplnxCore/test/ComputeSurfaceFeaturesTest.cpp index 62517dc3fe..6b7afc2e5f 100644 --- a/src/Plugins/SimplnxCore/test/ComputeSurfaceFeaturesTest.cpp +++ b/src/Plugins/SimplnxCore/test/ComputeSurfaceFeaturesTest.cpp @@ -3,8 +3,11 @@ #include "SimplnxCore/Filters/ReadRawBinaryFilter.hpp" #include "SimplnxCore/SimplnxCore_test_dirs.hpp" +#include "simplnx/DataStructure/AttributeMatrix.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" #include "simplnx/Parameters/DynamicTableParameter.hpp" #include "simplnx/UnitTest/UnitTestCommon.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" #include @@ -123,6 +126,9 @@ void test_impl(const std::vector& geometryDims, const std::string& featu TEST_CASE("SimplnxCore::ComputeSurfaceFeaturesFilter: 3D", "[SimplnxCore][ComputeSurfaceFeaturesFilter]") { UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 40000, true); const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_TestFilesDir, "6_5_test_data_1_v2.tar.gz", "6_5_test_data_1_v2"); // Read the Small IN100 Data set @@ -172,6 +178,9 @@ TEST_CASE("SimplnxCore::ComputeSurfaceFeaturesFilter: 3D", "[SimplnxCore][Comput TEST_CASE("SimplnxCore::ComputeSurfaceFeaturesFilter: 2D(XY Plane)", "[SimplnxCore][ComputeSurfaceFeaturesFilter]") { UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 40000, true); test_impl(std::vector({100, 100, 1}), k_FeatureIds2DFileName, k_SurfaceFeatures2DExemplaryFileName); } @@ -179,6 +188,9 @@ TEST_CASE("SimplnxCore::ComputeSurfaceFeaturesFilter: 2D(XY Plane)", "[SimplnxCo TEST_CASE("SimplnxCore::ComputeSurfaceFeaturesFilter: 2D(XZ Plane)", "[SimplnxCore][ComputeSurfaceFeaturesFilter]") { UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 40000, true); test_impl(std::vector({100, 1, 100}), k_FeatureIds2DFileName, k_SurfaceFeatures2DExemplaryFileName); } @@ -186,6 +198,72 @@ TEST_CASE("SimplnxCore::ComputeSurfaceFeaturesFilter: 2D(XZ Plane)", "[SimplnxCo TEST_CASE("SimplnxCore::ComputeSurfaceFeaturesFilter: 2D(YZ Plane)", "[SimplnxCore][ComputeSurfaceFeaturesFilter]") { UnitTest::LoadPlugins(); + bool forceOocAlgo = GENERATE(false, true); + const nx::core::ForceOocAlgorithmGuard guard(forceOocAlgo); + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 40000, true); test_impl(std::vector({1, 100, 100}), k_FeatureIds2DFileName, k_SurfaceFeatures2DExemplaryFileName); } + +TEST_CASE("SimplnxCore::ComputeSurfaceFeaturesFilter: Benchmark 200x200x200", "[SimplnxCore][ComputeSurfaceFeaturesFilter][Benchmark]") +{ + UnitTest::LoadPlugins(); + + constexpr usize kDimX = 200; + constexpr usize kDimY = 200; + constexpr usize kDimZ = 200; + constexpr usize kTotalVoxels = kDimX * kDimY * kDimZ; + const ShapeType cellTupleShape = {kDimZ, kDimY, kDimX}; + const UnitTest::PreferencesSentinel prefsSentinel("Zarr", 160000, true); + + const auto benchmarkFile = fs::path(fmt::format("{}/compute_surface_features_benchmark.dream3d", unit_test::k_BinaryTestOutputDir)); + + const std::string k_BenchGeomName = "ImageGeom"; + const DataPath k_BenchGeomPath({k_BenchGeomName}); + const DataPath k_BenchFeatureIdsPath = k_BenchGeomPath.createChildPath(Constants::k_CellData).createChildPath(Constants::k_FeatureIds); + const DataPath k_BenchCellFeatureAMPath = k_BenchGeomPath.createChildPath(Constants::k_CellFeatureData); + + { + DataStructure buildDS; + auto* imageGeom = ImageGeom::Create(buildDS, k_BenchGeomName); + imageGeom->setDimensions({kDimX, kDimY, kDimZ}); + imageGeom->setSpacing({1.0f, 1.0f, 1.0f}); + imageGeom->setOrigin({0.0f, 0.0f, 0.0f}); + + auto* cellAM = AttributeMatrix::Create(buildDS, Constants::k_CellData, cellTupleShape, imageGeom->getId()); + imageGeom->setCellData(*cellAM); + + constexpr int32 kNumFeatures = 8; + auto* featureIds = UnitTest::CreateTestDataArray(buildDS, Constants::k_FeatureIds, cellTupleShape, {1}, cellAM->getId()); + auto& fidsStore = featureIds->getDataStoreRef(); + for(usize i = 0; i < kTotalVoxels; i++) + { + const usize ix = i % kDimX; + const usize iy = (i / kDimX) % kDimY; + const usize iz = i / (kDimX * kDimY); + const int32 octant = static_cast((iz >= kDimZ / 2 ? 4 : 0) + (iy >= kDimY / 2 ? 2 : 0) + (ix >= kDimX / 2 ? 1 : 0)) + 1; + fidsStore.setValue(i, ((i * 7 + 13) % 100) < 5 ? 0 : octant); + } + + AttributeMatrix::Create(buildDS, Constants::k_CellFeatureData, {static_cast(kNumFeatures + 1)}, imageGeom->getId()); + UnitTest::WriteTestDataStructure(buildDS, benchmarkFile); + } + + DataStructure dataStructure = UnitTest::LoadDataStructure(benchmarkFile); + + { + ComputeSurfaceFeaturesFilter filter; + Arguments args; + args.insertOrAssign(ComputeSurfaceFeaturesFilter::k_FeatureGeometryPath_Key, std::make_any(k_BenchGeomPath)); + args.insertOrAssign(ComputeSurfaceFeaturesFilter::k_CellFeatureIdsArrayPath_Key, std::make_any(k_BenchFeatureIdsPath)); + args.insertOrAssign(ComputeSurfaceFeaturesFilter::k_CellFeatureAttributeMatrixPath_Key, std::make_any(k_BenchCellFeatureAMPath)); + args.insertOrAssign(ComputeSurfaceFeaturesFilter::k_SurfaceFeaturesArrayName_Key, std::make_any("SurfaceFeatures")); + + auto preflightResult = filter.preflight(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions); + auto executeResult = filter.execute(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result); + } + + fs::remove(benchmarkFile); +} From c838f62f05a07901de8b2832e2473a55749f41a7 Mon Sep 17 00:00:00 2001 From: Joey Kleingers Date: Thu, 5 Mar 2026 14:30:57 -0500 Subject: [PATCH 3/5] FIX: Wire DispatchAlgorithm calls into base algorithm classes The 4 base algorithm .cpp files (ComputeBoundaryCells, ComputeSurfaceFeatures, ComputeFeatureNeighbors, ComputeSurfaceAreaToVolume) still contained the original algorithm code instead of dispatching to the Direct/Scanline classes. This replaces each operator()() body with a DispatchAlgorithm call. Also fixes ImageGeom copy-by-value in ComputeBoundaryCellsDirect/Scanline (const ImageGeom -> const auto&) and adds Doxygen to all 8 new headers. --- .../Algorithms/ComputeBoundaryCells.cpp | 90 +------- .../Algorithms/ComputeBoundaryCellsDirect.cpp | 2 +- .../Algorithms/ComputeBoundaryCellsDirect.hpp | 6 + .../ComputeBoundaryCellsScanline.cpp | 2 +- .../ComputeBoundaryCellsScanline.hpp | 7 + .../Algorithms/ComputeFeatureNeighbors.cpp | 206 +---------------- .../ComputeFeatureNeighborsDirect.hpp | 7 + .../ComputeFeatureNeighborsScanline.hpp | 7 + .../Algorithms/ComputeSurfaceAreaToVolume.cpp | 141 +----------- .../ComputeSurfaceAreaToVolumeDirect.hpp | 6 + .../ComputeSurfaceAreaToVolumeScanline.hpp | 7 + .../Algorithms/ComputeSurfaceFeatures.cpp | 214 +----------------- .../ComputeSurfaceFeaturesDirect.hpp | 6 + .../ComputeSurfaceFeaturesScanline.hpp | 7 + 14 files changed, 79 insertions(+), 629 deletions(-) diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCells.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCells.cpp index 83be2af03a..1f6ad3a5c3 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCells.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCells.cpp @@ -1,8 +1,10 @@ #include "ComputeBoundaryCells.hpp" +#include "ComputeBoundaryCellsDirect.hpp" +#include "ComputeBoundaryCellsScanline.hpp" + #include "simplnx/DataStructure/DataArray.hpp" -#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" -#include "simplnx/Utilities/NeighborUtilities.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" using namespace nx::core; @@ -27,86 +29,6 @@ const std::atomic_bool& ComputeBoundaryCells::getCancel() // ----------------------------------------------------------------------------- Result<> ComputeBoundaryCells::operator()() { - const ImageGeom imageGeometry = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); - const SizeVec3 udims = imageGeometry.getDimensions(); - std::array dims = { - static_cast(udims[0]), - static_cast(udims[1]), - static_cast(udims[2]), - }; - - auto& featureIdsStore = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath)->getDataStoreRef(); - auto& boundaryCellsStore = m_DataStructure.getDataAs(m_InputValues->BoundaryCellsArrayName)->getDataStoreRef(); - - std::array neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims); - std::array faceNeighborInternalIdx = initializeFaceNeighborInternalIdx(); - - int32 feature = 0; - int8 onSurf = 0; - int64 neighborPoint = 0; - - int ignoreFeatureZeroVal = 0; - if(!m_InputValues->IgnoreFeatureZero) - { - ignoreFeatureZeroVal = -1; - } - - int64 kStride = 0; - int64 jStride = 0; - - for(int64 zIdx = 0; zIdx < dims[2]; zIdx++) - { - kStride = dims[0] * dims[1] * zIdx; - for(int64 yIdx = 0; yIdx < dims[1]; yIdx++) - { - jStride = dims[0] * yIdx; - for(int64 xIdx = 0; xIdx < dims[0]; xIdx++) - { - int64 voxelIndex = kStride + jStride + xIdx; - onSurf = 0; - feature = featureIdsStore[voxelIndex]; - if(feature >= 0) - { - if(m_InputValues->IncludeVolumeBoundary) - { - if(dims[0] > 2 && (xIdx == 0 || xIdx == dims[0] - 1)) - { - onSurf++; - } - if(dims[1] > 2 && (yIdx == 0 || yIdx == dims[1] - 1)) - { - onSurf++; - } - if(dims[2] > 2 && (zIdx == 0 || zIdx == dims[2] - 1)) - { - onSurf++; - } - - if(onSurf > 0 && feature == 0) - { - onSurf = 0; - } - } - - // Loop over the 6 face neighbors of the voxel - std::array isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims); - for(const auto& faceIndex : faceNeighborInternalIdx) - { - if(!isValidFaceNeighbor[faceIndex]) - { - continue; - } - neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; - - if(featureIdsStore[neighborPoint] != feature && featureIdsStore[neighborPoint] > ignoreFeatureZeroVal) - { - onSurf++; - } - } - } - boundaryCellsStore[voxelIndex] = onSurf; - } - } - } - return {}; + auto* featureIdsArray = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath); + return DispatchAlgorithm({featureIdsArray}, m_DataStructure, m_MessageHandler, m_ShouldCancel, m_InputValues); } diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp index e3c9e3e58d..edde51e412 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp @@ -24,7 +24,7 @@ ComputeBoundaryCellsDirect::~ComputeBoundaryCellsDirect() noexcept = default; // ----------------------------------------------------------------------------- Result<> ComputeBoundaryCellsDirect::operator()() { - const ImageGeom imageGeometry = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); + const auto& imageGeometry = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); const SizeVec3 udims = imageGeometry.getDimensions(); std::array dims = { static_cast(udims[0]), diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.hpp index 35f6a06bc7..b2ed65ffe4 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.hpp @@ -9,6 +9,12 @@ namespace nx::core { struct ComputeBoundaryCellsInputValues; +/** + * @class ComputeBoundaryCellsDirect + * @brief In-core algorithm for ComputeBoundaryCells. Preserves the original sequential + * Z-Y-X voxel iteration with face-neighbor boundary counting. Selected by DispatchAlgorithm + * when all input arrays are backed by in-memory DataStore. + */ class SIMPLNXCORE_EXPORT ComputeBoundaryCellsDirect { public: diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp index 715d9a8055..9a3e5f783e 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp @@ -24,7 +24,7 @@ ComputeBoundaryCellsScanline::~ComputeBoundaryCellsScanline() noexcept = default // ----------------------------------------------------------------------------- Result<> ComputeBoundaryCellsScanline::operator()() { - const ImageGeom imageGeometry = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); + const auto& imageGeometry = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); const SizeVec3 udims = imageGeometry.getDimensions(); std::array dims = { static_cast(udims[0]), diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.hpp index 9ca5742150..0f01d6d78a 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.hpp @@ -9,6 +9,13 @@ namespace nx::core { struct ComputeBoundaryCellsInputValues; +/** + * @class ComputeBoundaryCellsScanline + * @brief Out-of-core algorithm for ComputeBoundaryCells. Wraps the voxel iteration in + * chunk-sequential access using getNumberOfChunks/loadChunk/getChunkLowerBounds/getChunkUpperBounds + * for guaranteed sequential disk I/O. Selected by DispatchAlgorithm when any input array + * is backed by ZarrStore. + */ class SIMPLNXCORE_EXPORT ComputeBoundaryCellsScanline { public: diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighbors.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighbors.cpp index e4d8cda73a..f571f25996 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighbors.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighbors.cpp @@ -1,13 +1,10 @@ #include "ComputeFeatureNeighbors.hpp" -#include "simplnx/DataStructure/AttributeMatrix.hpp" -#include "simplnx/DataStructure/DataArray.hpp" -#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" -#include "simplnx/DataStructure/NeighborList.hpp" -#include "simplnx/Utilities/MessageHelper.hpp" -#include "simplnx/Utilities/NeighborUtilities.hpp" +#include "ComputeFeatureNeighborsDirect.hpp" +#include "ComputeFeatureNeighborsScanline.hpp" -#include +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" using namespace nx::core; @@ -27,197 +24,6 @@ ComputeFeatureNeighbors::~ComputeFeatureNeighbors() noexcept = default; // ----------------------------------------------------------------------------- Result<> ComputeFeatureNeighbors::operator()() { - auto storeBoundaryCells = m_InputValues->StoreBoundaryCells; - auto storeSurfaceFeatures = m_InputValues->StoreSurfaceFeatures; - auto imageGeomPath = m_InputValues->InputImageGeometryPath; - auto featureIdsPath = m_InputValues->FeatureIdsPath; - auto boundaryCellsName = m_InputValues->BoundaryCellsName; - auto numNeighborsName = m_InputValues->NumberOfNeighborsName; - auto neighborListName = m_InputValues->NeighborListName; - auto sharedSurfaceAreaName = m_InputValues->SharedSurfaceAreaListName; - auto surfaceFeaturesName = m_InputValues->SurfaceFeaturesName; - auto featureAttrMatrixPath = m_InputValues->CellFeatureArrayPath; - - DataPath boundaryCellsPath = featureIdsPath.replaceName(boundaryCellsName); - DataPath numNeighborsPath = featureAttrMatrixPath.createChildPath(numNeighborsName); - DataPath neighborListPath = featureAttrMatrixPath.createChildPath(neighborListName); - DataPath sharedSurfaceAreaPath = featureAttrMatrixPath.createChildPath(sharedSurfaceAreaName); - DataPath surfaceFeaturesPath = featureAttrMatrixPath.createChildPath(surfaceFeaturesName); - - auto& featureIds = m_DataStructure.getDataAs(featureIdsPath)->getDataStoreRef(); - auto& numNeighbors = m_DataStructure.getDataAs(numNeighborsPath)->getDataStoreRef(); - - auto& neighborList = m_DataStructure.getDataRefAs(neighborListPath); - auto& sharedSurfaceAreaList = m_DataStructure.getDataRefAs(sharedSurfaceAreaPath); - - auto* boundaryCells = storeBoundaryCells ? m_DataStructure.getDataAs(boundaryCellsPath)->getDataStore() : nullptr; - auto* surfaceFeatures = storeSurfaceFeatures ? m_DataStructure.getDataAs(surfaceFeaturesPath)->getDataStore() : nullptr; - - usize totalPoints = featureIds.getNumberOfTuples(); - usize totalFeatures = numNeighbors.getNumberOfTuples(); - - /* Ensure that we will be able to work with the user selected featureId Array */ - const auto [minFeatureId, maxFeatureId] = std::minmax_element(featureIds.begin(), featureIds.end()); - if(static_cast(*maxFeatureId) >= totalFeatures) - { - std::stringstream out; - out << "Data Array " << featureIdsPath.getTargetName() << " has a maximum value of " << *maxFeatureId << " which is greater than the " - << " number of features from array " << numNeighborsPath.getTargetName() << " which has " << totalFeatures << ". Did you select the " - << " incorrect array for the 'FeatureIds' array?"; - return MakeErrorResult(-24500, out.str()); - } - - auto& imageGeom = m_DataStructure.getDataRefAs(imageGeomPath); - SizeVec3 uDims = imageGeom.getDimensions(); - - std::array dims = { - static_cast(uDims[0]), - static_cast(uDims[1]), - static_cast(uDims[2]), - }; - - std::array neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims); - std::array faceNeighborInternalIdx = initializeFaceNeighborInternalIdx(); - - int32 feature = 0; - int32 nnum = 0; - uint8 onsurf = 0; - - std::vector> neighborlist(totalFeatures); - std::vector> neighborsurfacearealist(totalFeatures); - - int32 nListSize = 100; - - MessageHelper messageHelper(m_MessageHandler); - ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger(); - // Initialize the neighbor lists - for(usize featureIdx = 1; featureIdx < totalFeatures; featureIdx++) - { - auto now = std::chrono::steady_clock::now(); - throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Initializing Neighbor Lists || {:.2f}% Complete", CalculatePercentComplete(featureIdx, totalFeatures)); }); - - if(m_ShouldCancel) - { - return {}; - } - - numNeighbors[featureIdx] = 0; - neighborlist[featureIdx].resize(nListSize); - neighborsurfacearealist[featureIdx].assign(nListSize, -1.0f); - if(storeSurfaceFeatures && surfaceFeatures != nullptr) - { - surfaceFeatures->setValue(featureIdx, false); - } - } - - // Loop over all points to generate the neighbor lists - for(int64 voxelIndex = 0; voxelIndex < totalPoints; voxelIndex++) - { - throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Determining Neighbor Lists || {:.2f}% Complete", CalculatePercentComplete(voxelIndex, totalPoints)); }); - - if(m_ShouldCancel) - { - return {}; - } - - onsurf = 0; - feature = featureIds[voxelIndex]; - if(feature > 0 && feature < neighborlist.size()) - { - int64 xIdx = voxelIndex % dims[0]; - int64 yIdx = (voxelIndex / dims[0]) % dims[1]; - int64 zIdx = voxelIndex / (dims[0] * dims[1]); - - if(storeSurfaceFeatures && surfaceFeatures != nullptr) - { - if((xIdx == 0 || xIdx == static_cast((dims[0] - 1)) || yIdx == 0 || yIdx == static_cast((dims[1]) - 1) || zIdx == 0 || zIdx == static_cast((dims[2] - 1))) && dims[2] != 1) - { - surfaceFeatures->setValue(feature, true); - } - if((xIdx == 0 || xIdx == static_cast((dims[0] - 1)) || yIdx == 0 || yIdx == static_cast((dims[1] - 1))) && dims[2] == 1) - { - surfaceFeatures->setValue(feature, true); - } - } - - // Loop over the 6 face neighbors of the voxel - std::array isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims); - for(const auto& faceIndex : faceNeighborInternalIdx) - { - if(!isValidFaceNeighbor[faceIndex]) - { - continue; - } - - const int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex]; - - if(featureIds[neighborPoint] != feature && featureIds[neighborPoint] > 0) - { - onsurf++; - nnum = numNeighbors[feature]; - neighborlist[feature].push_back(featureIds[neighborPoint]); - nnum++; - numNeighbors[feature] = nnum; - } - } - } - if(storeBoundaryCells && boundaryCells != nullptr) - { - boundaryCells->setValue(voxelIndex, static_cast(onsurf)); - } - } - - FloatVec3 spacing = imageGeom.getSpacing(); - - // We do this to create new set of NeighborList objects - for(usize i = 1; i < totalFeatures; i++) - { - throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Calculating Surface Areas || {:.2f}% Complete", CalculatePercentComplete(i, totalFeatures)); }); - - if(m_ShouldCancel) - { - return {}; - } - - std::map neighToCount; - auto numneighs = static_cast(neighborlist[i].size()); - - // this increments the voxel counts for each feature - for(int32 j = 0; j < numneighs; j++) - { - neighToCount[neighborlist[i][j]]++; - } - - auto neighborIter = neighToCount.find(0); - neighToCount.erase(neighborIter); - neighborIter = neighToCount.find(-1); - if(neighborIter != neighToCount.end()) - { - neighToCount.erase(neighborIter); - } - // Resize the features neighbor list to zero - neighborlist[i].resize(0); - neighborsurfacearealist[i].resize(0); - - for(const auto [neigh, number] : neighToCount) - { - float area = static_cast(number) * spacing[0] * spacing[1]; - - // Push the neighbor feature identifier back onto the list, so we stay synced up - neighborlist[i].push_back(neigh); - neighborsurfacearealist[i].push_back(area); - } - numNeighbors[i] = static_cast(neighborlist[i].size()); - - // Set the vector for each list into the NeighborList Object - NeighborList::SharedVectorType sharedNeiLst(new std::vector); - sharedNeiLst->assign(neighborlist[i].begin(), neighborlist[i].end()); - neighborList.setList(static_cast(i), sharedNeiLst); - - NeighborList::SharedVectorType sharedSAL(new std::vector); - sharedSAL->assign(neighborsurfacearealist[i].begin(), neighborsurfacearealist[i].end()); - sharedSurfaceAreaList.setList(static_cast(i), sharedSAL); - } - - return {}; + auto* featureIdsArray = m_DataStructure.getDataAs(m_InputValues->FeatureIdsPath); + return DispatchAlgorithm({featureIdsArray}, m_DataStructure, m_MessageHandler, m_ShouldCancel, m_InputValues); } diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.hpp index aace0b087c..0bc4683414 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.hpp @@ -9,6 +9,13 @@ namespace nx::core { struct ComputeFeatureNeighborsInputValues; +/** + * @class ComputeFeatureNeighborsDirect + * @brief In-core algorithm for ComputeFeatureNeighbors. Preserves the original two-phase + * algorithm: Phase 1 iterates all voxels to build per-feature neighbor lists, Phase 2 + * computes shared surface areas. Selected by DispatchAlgorithm when all input arrays + * are backed by in-memory DataStore. + */ class SIMPLNXCORE_EXPORT ComputeFeatureNeighborsDirect { public: diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.hpp index 9d9bb71916..df3bec43ab 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.hpp @@ -9,6 +9,13 @@ namespace nx::core { struct ComputeFeatureNeighborsInputValues; +/** + * @class ComputeFeatureNeighborsScanline + * @brief Out-of-core algorithm for ComputeFeatureNeighbors. Wraps only Phase 1 (voxel iteration) + * in chunk-sequential access for guaranteed sequential disk I/O on ZarrStore-backed arrays. + * Phase 2 (feature-level surface area computation) is unchanged. Selected by DispatchAlgorithm + * when any input array is backed by ZarrStore. + */ class SIMPLNXCORE_EXPORT ComputeFeatureNeighborsScanline { public: diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolume.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolume.cpp index e534da769c..155fb3188a 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolume.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolume.cpp @@ -1,10 +1,10 @@ #include "ComputeSurfaceAreaToVolume.hpp" -#include "simplnx/Common/Constants.hpp" +#include "ComputeSurfaceAreaToVolumeDirect.hpp" +#include "ComputeSurfaceAreaToVolumeScanline.hpp" + #include "simplnx/DataStructure/DataArray.hpp" -#include "simplnx/DataStructure/DataGroup.hpp" -#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" -#include "simplnx/Utilities/DataArrayUtilities.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" using namespace nx::core; @@ -30,135 +30,6 @@ const std::atomic_bool& ComputeSurfaceAreaToVolume::getCancel() // ----------------------------------------------------------------------------- Result<> ComputeSurfaceAreaToVolume::operator()() { - // Input Cell Data - auto featureIdsArrayPtr = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath); - const auto& featureIdsStoreRef = featureIdsArrayPtr->getDataStoreRef(); - - // Input Feature Data - const auto& numCells = m_DataStructure.getDataAs(m_InputValues->NumCellsArrayPath)->getDataStoreRef(); - - // Output Feature Data - auto& surfaceAreaVolumeRatio = m_DataStructure.getDataAs(m_InputValues->SurfaceAreaVolumeRatioArrayName)->getDataStoreRef(); - - // Required Geometry - const auto& imageGeom = m_DataStructure.getDataRefAs(m_InputValues->InputImageGeometry); - - auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, m_InputValues->NumCellsArrayPath.getParent(), *featureIdsArrayPtr, false, m_MessageHandler); - if(validateNumFeatResult.invalid()) - { - return validateNumFeatResult; - } - auto numFeatures = static_cast(numCells.getNumberOfTuples()); - SizeVec3 dims = imageGeom.getDimensions(); - FloatVec3 spacing = imageGeom.getSpacing(); - - auto xPoints = static_cast(dims[0]); - auto yPoints = static_cast(dims[1]); - auto zPoints = static_cast(dims[2]); - - float32 voxelVol = spacing[0] * spacing[1] * spacing[2]; - - std::vector featureSurfaceArea(static_cast(numFeatures), 0.0f); - - // This stores an offset to get to a particular index in the array based on - // a normal orthogonal cube - int64 neighborOffset[6] = {0, 0, 0, 0, 0, 0}; - neighborOffset[0] = -xPoints * yPoints; // -Z - neighborOffset[1] = -xPoints; // -Y - neighborOffset[2] = -1; // -X - neighborOffset[3] = 1; // +X - neighborOffset[4] = xPoints; // +Y - neighborOffset[5] = xPoints * yPoints; // +Z - - // Start looping over the regular grid data (This could be either an Image Geometry or a Rectilinear Grid geometry (in theory) - for(int64 zIdx = 0; zIdx < zPoints; zIdx++) - { - m_MessageHandler(IFilter::Message::Type::Info, fmt::format("Computing Z Slice: '{}'", zIdx)); - - int64 zStride = zIdx * xPoints * yPoints; - for(int64 yIdx = 0; yIdx < yPoints; yIdx++) - { - int64 yStride = yIdx * xPoints; - for(int64 xIdx = 0; xIdx < xPoints; xIdx++) - { - float onSurface = 0.0f; // Start totalling the surface area - int32 currentFeatureId = featureIdsStoreRef[zStride + yStride + xIdx]; - // If the current feature ID is not valid (< 1), then just continue; - if(currentFeatureId < 1) - { - continue; - } - - // Loop over all 6 face neighbors - for(int32 neighborOffsetIndex = 0; neighborOffsetIndex < 6; neighborOffsetIndex++) - { - if(neighborOffsetIndex == 0 && zIdx == 0) // if we are on the bottom Z Layer, skip - { - continue; - } - if(neighborOffsetIndex == 5 && zIdx == (zPoints - 1)) // if we are on the top Z Layer, skip - { - continue; - } - if(neighborOffsetIndex == 1 && yIdx == 0) // If we are on the first Y row, skip - { - continue; - } - if(neighborOffsetIndex == 4 && yIdx == (yPoints - 1)) // If we are on the last Y row, skip - { - continue; - } - if(neighborOffsetIndex == 2 && xIdx == 0) // If we are on the first X column, skip - { - continue; - } - if(neighborOffsetIndex == 3 && xIdx == (xPoints - 1)) // If we are on the last X column, skip - { - continue; - } - // - int64 neighborIndex = zStride + yStride + xIdx + neighborOffset[neighborOffsetIndex]; - - if(featureIdsStoreRef[neighborIndex] != currentFeatureId) - { - if(neighborOffsetIndex == 0 || neighborOffsetIndex == 5) // XY face shared - { - onSurface = onSurface + spacing[0] * spacing[1]; - } - if(neighborOffsetIndex == 1 || neighborOffsetIndex == 4) // YZ face shared - { - onSurface = onSurface + spacing[1] * spacing[2]; - } - if(neighborOffsetIndex == 2 || neighborOffsetIndex == 3) // XZ face shared - { - onSurface = onSurface + spacing[2] * spacing[0]; - } - } - } - int32 featureId = featureIdsStoreRef[zStride + yStride + xIdx]; - featureSurfaceArea[featureId] = featureSurfaceArea[featureId] + onSurface; - } - } - } - - const float32 thirdRootPi = std::pow(nx::core::Constants::k_PiF, 0.333333f); - for(usize i = 1; i < numFeatures; i++) - { - float featureVolume = voxelVol * numCells[i]; - surfaceAreaVolumeRatio[i] = featureSurfaceArea[i] / featureVolume; - } - - if(m_InputValues->CalculateSphericity) // Calc the sphericity if requested - { - m_MessageHandler(IFilter::Message::Type::Info, fmt::format("Computing Sphericity")); - - auto& sphericity = m_DataStructure.getDataAs(m_InputValues->SphericityArrayName)->getDataStoreRef(); - for(usize i = 1; i < static_cast(numFeatures); i++) - { - float featureVolume = voxelVol * numCells[i]; - sphericity[i] = (thirdRootPi * std::pow((6.0f * featureVolume), 0.66666f)) / featureSurfaceArea[i]; - } - } - - return {}; + auto* featureIdsArray = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath); + return DispatchAlgorithm({featureIdsArray}, m_DataStructure, m_MessageHandler, m_ShouldCancel, m_InputValues); } diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.hpp index 1814f8f18f..54e3160442 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.hpp @@ -9,6 +9,12 @@ namespace nx::core { struct ComputeSurfaceAreaToVolumeInputValues; +/** + * @class ComputeSurfaceAreaToVolumeDirect + * @brief In-core algorithm for ComputeSurfaceAreaToVolume. Preserves the original sequential + * Z-Y-X voxel iteration with face-neighbor surface area accumulation per feature. + * Selected by DispatchAlgorithm when all input arrays are backed by in-memory DataStore. + */ class SIMPLNXCORE_EXPORT ComputeSurfaceAreaToVolumeDirect { public: diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp index 404bd7cf20..20b14a2bd8 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp @@ -9,6 +9,13 @@ namespace nx::core { struct ComputeSurfaceAreaToVolumeInputValues; +/** + * @class ComputeSurfaceAreaToVolumeScanline + * @brief Out-of-core algorithm for ComputeSurfaceAreaToVolume. Wraps the voxel iteration in + * chunk-sequential access for guaranteed sequential disk I/O on ZarrStore-backed arrays. + * Feature-level ratio and sphericity computations are unchanged. Selected by DispatchAlgorithm + * when any input array is backed by ZarrStore. + */ class SIMPLNXCORE_EXPORT ComputeSurfaceAreaToVolumeScanline { public: diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeatures.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeatures.cpp index 9dcf86a149..6997eaeafb 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeatures.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeatures.cpp @@ -1,185 +1,13 @@ #include "ComputeSurfaceFeatures.hpp" +#include "ComputeSurfaceFeaturesDirect.hpp" +#include "ComputeSurfaceFeaturesScanline.hpp" + #include "simplnx/DataStructure/DataArray.hpp" -#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" -#include "simplnx/Utilities/DataArrayUtilities.hpp" +#include "simplnx/Utilities/AlgorithmDispatch.hpp" using namespace nx::core; -namespace -{ -bool IsPointASurfaceFeature(const Point2D& point, usize xPoints, usize yPoints, bool markFeature0Neighbors, const Int32AbstractDataStore& featureIds) -{ - const usize yStride = point.getY() * xPoints; - - if(point.getX() <= 0 || point.getX() >= xPoints - 1) - { - return true; - } - if(point.getY() <= 0 || point.getY() >= yPoints - 1) - { - return true; - } - - if(markFeature0Neighbors) - { - if(featureIds[yStride + point.getX() - 1] == 0) - { - return true; - } - if(featureIds[yStride + point.getX() + 1] == 0) - { - return true; - } - if(featureIds[yStride + point.getX() - xPoints] == 0) - { - return true; - } - if(featureIds[yStride + point.getX() + xPoints] == 0) - { - return true; - } - } - - return false; -} - -bool IsPointASurfaceFeature(const Point3D& point, usize xPoints, usize yPoints, usize zPoints, bool markFeature0Neighbors, const Int32AbstractDataStore& featureIds) -{ - usize yStride = point.getY() * xPoints; - usize zStride = point.getZ() * xPoints * yPoints; - - if(point.getX() <= 0 || point.getX() >= xPoints - 1) - { - return true; - } - if(point.getY() <= 0 || point.getY() >= yPoints - 1) - { - return true; - } - if(point.getZ() <= 0 || point.getZ() >= zPoints - 1) - { - return true; - } - - if(markFeature0Neighbors) - { - if(featureIds[zStride + yStride + point.getX() - 1] == 0) - { - return true; - } - if(featureIds[zStride + yStride + point.getX() + 1] == 0) - { - return true; - } - if(featureIds[zStride + yStride + point.getX() - xPoints] == 0) - { - return true; - } - if(featureIds[zStride + yStride + point.getX() + xPoints] == 0) - { - return true; - } - if(featureIds[zStride + yStride + point.getX() - (xPoints * yPoints)] == 0) - { - return true; - } - if(featureIds[zStride + yStride + point.getX() + (xPoints * yPoints)] == 0) - { - return true; - } - } - - return false; -} - -void findSurfaceFeatures3D(DataStructure& dataStructure, const DataPath& featureGeometryPathValue, const DataPath& featureIdsArrayPathValue, const DataPath& surfaceFeaturesArrayPathValue, - bool markFeature0Neighbors, const std::atomic_bool& shouldCancel) -{ - const auto& featureGeometry = dataStructure.getDataRefAs(featureGeometryPathValue); - const auto& featureIds = dataStructure.getDataAs(featureIdsArrayPathValue)->getDataStoreRef(); - auto& surfaceFeatures = dataStructure.getDataAs(surfaceFeaturesArrayPathValue)->getDataStoreRef(); - - const usize xPoints = featureGeometry.getNumXCells(); - const usize yPoints = featureGeometry.getNumYCells(); - const usize zPoints = featureGeometry.getNumZCells(); - - for(usize z = 0; z < zPoints; z++) - { - const usize zStride = z * xPoints * yPoints; - for(usize y = 0; y < yPoints; y++) - { - const usize yStride = y * xPoints; - for(usize x = 0; x < xPoints; x++) - { - if(shouldCancel) - { - return; - } - - const int32 gNum = featureIds[zStride + yStride + x]; - if(gNum != 0 && !surfaceFeatures[gNum]) - { - if(IsPointASurfaceFeature(Point3D{x, y, z}, xPoints, yPoints, zPoints, markFeature0Neighbors, featureIds)) - { - surfaceFeatures[gNum] = 1; - } - } - } - } - } -} - -void findSurfaceFeatures2D(DataStructure& dataStructure, const DataPath& featureGeometryPathValue, const DataPath& featureIdsArrayPathValue, const DataPath& surfaceFeaturesArrayPathValue, - bool markFeature0Neighbors, const std::atomic_bool& shouldCancel) -{ - const auto& featureGeometry = dataStructure.getDataRefAs(featureGeometryPathValue); - const auto& featureIds = dataStructure.getDataAs(featureIdsArrayPathValue)->getDataStoreRef(); - auto& surfaceFeatures = dataStructure.getDataAs(surfaceFeaturesArrayPathValue)->getDataStoreRef(); - - usize xPoints = 0; - usize yPoints = 0; - - if(featureGeometry.getNumXCells() == 1) - { - xPoints = featureGeometry.getNumYCells(); - yPoints = featureGeometry.getNumZCells(); - } - if(featureGeometry.getNumYCells() == 1) - { - xPoints = featureGeometry.getNumXCells(); - yPoints = featureGeometry.getNumZCells(); - } - if(featureGeometry.getNumZCells() == 1) - { - xPoints = featureGeometry.getNumXCells(); - yPoints = featureGeometry.getNumYCells(); - } - - for(usize y = 0; y < yPoints; y++) - { - const usize yStride = y * xPoints; - - for(usize x = 0; x < xPoints; x++) - { - if(shouldCancel) - { - return; - } - - const int32 gNum = featureIds[yStride + x]; - if(gNum != 0 && surfaceFeatures[gNum] == 0) - { - if(IsPointASurfaceFeature(Point2D{x, y}, xPoints, yPoints, markFeature0Neighbors, featureIds)) - { - surfaceFeatures[gNum] = 1; - } - } - } - } -} -} // namespace - // ----------------------------------------------------------------------------- ComputeSurfaceFeatures::ComputeSurfaceFeatures(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, ComputeSurfaceFeaturesInputValues* inputValues) @@ -196,36 +24,6 @@ ComputeSurfaceFeatures::~ComputeSurfaceFeatures() noexcept = default; // ----------------------------------------------------------------------------- Result<> ComputeSurfaceFeatures::operator()() { - - const auto pMarkFeature0NeighborsValue = m_InputValues->MarkFeature0Neighbors; - const auto pFeatureGeometryPathValue = m_InputValues->InputImageGeometryPath; - const auto pFeatureIdsArrayPathValue = m_InputValues->FeatureIdsPath; - const auto pFeaturesAttributeMatrixPathValue = m_InputValues->FeatureAttributeMatrixPath; - const auto pSurfaceFeaturesArrayPathValue = pFeaturesAttributeMatrixPathValue.createChildPath(m_InputValues->SurfaceFeaturesArrayName); - - // Resize the surface features array to the proper size - const auto& featureIdsArray = m_DataStructure.getDataRefAs(pFeatureIdsArrayPathValue); - - auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, pFeaturesAttributeMatrixPathValue, featureIdsArray, false, m_MessageHandler); - if(validateNumFeatResult.invalid()) - { - return validateNumFeatResult; - } - - // Find surface features - const auto& featureGeometry = m_DataStructure.getDataRefAs(pFeatureGeometryPathValue); - if(const usize geometryDimensionality = featureGeometry.getDimensionality(); geometryDimensionality == 3) - { - findSurfaceFeatures3D(m_DataStructure, pFeatureGeometryPathValue, pFeatureIdsArrayPathValue, pSurfaceFeaturesArrayPathValue, pMarkFeature0NeighborsValue, m_ShouldCancel); - } - else if(geometryDimensionality == 2) - { - findSurfaceFeatures2D(m_DataStructure, pFeatureGeometryPathValue, pFeatureIdsArrayPathValue, pSurfaceFeaturesArrayPathValue, pMarkFeature0NeighborsValue, m_ShouldCancel); - } - else - { - MakeErrorResult(-1000, fmt::format("Image Geometry at path '{}' must be either 3D or 2D", pFeatureGeometryPathValue.toString())); - } - - return {}; + auto* featureIdsArray = m_DataStructure.getDataAs(m_InputValues->FeatureIdsPath); + return DispatchAlgorithm({featureIdsArray}, m_DataStructure, m_MessageHandler, m_ShouldCancel, m_InputValues); } diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp index d7810f2502..172a29a51d 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp @@ -9,6 +9,12 @@ namespace nx::core { struct ComputeSurfaceFeaturesInputValues; +/** + * @class ComputeSurfaceFeaturesDirect + * @brief In-core algorithm for ComputeSurfaceFeatures. Preserves the original 2D/3D branching + * with sequential voxel iteration and face-neighbor surface detection. Selected by + * DispatchAlgorithm when all input arrays are backed by in-memory DataStore. + */ class SIMPLNXCORE_EXPORT ComputeSurfaceFeaturesDirect { public: diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp index a5a6c23bf2..d3be471d12 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp @@ -9,6 +9,13 @@ namespace nx::core { struct ComputeSurfaceFeaturesInputValues; +/** + * @class ComputeSurfaceFeaturesScanline + * @brief Out-of-core algorithm for ComputeSurfaceFeatures. Uses chunk-sequential 3D iteration + * with 2D coordinate remapping for degenerate dimensions, ensuring sequential disk I/O + * on ZarrStore-backed arrays. Selected by DispatchAlgorithm when any input array is + * backed by ZarrStore. + */ class SIMPLNXCORE_EXPORT ComputeSurfaceFeaturesScanline { public: From 0aabf6d9e212a03dca178f70a4e83df3ccc2748e Mon Sep 17 00:00:00 2001 From: Joey Kleingers Date: Thu, 5 Mar 2026 14:57:31 -0500 Subject: [PATCH 4/5] FIX: PR review fixes for Group B face-neighbor filters - Fix missing return before MakeErrorResult in ComputeSurfaceFeaturesDirect for 1D geometry error path (pre-existing bug). Add same error path to ComputeSurfaceFeaturesScanline. - Replace std::minmax_element full-scan validation with deferred max tracking during the main loop in ComputeFeatureNeighborsDirect and Scanline, eliminating a redundant OOC full-scan (1.4x OOC speedup). - Add Doxygen @brief comments to operator()() in all 8 new algorithm classes. --- ...adDataNeighborOrientationCheckScanline.cpp | 7 ++++ ...adDataNeighborOrientationCheckWorklist.cpp | 6 ++++ .../Algorithms/ComputeBoundaryCellsDirect.cpp | 4 +++ .../ComputeBoundaryCellsScanline.cpp | 5 +++ .../ComputeFeatureNeighborsDirect.cpp | 34 +++++++++++------- .../ComputeFeatureNeighborsScanline.cpp | 35 ++++++++++++------- .../ComputeSurfaceAreaToVolumeDirect.cpp | 5 +++ .../ComputeSurfaceAreaToVolumeScanline.cpp | 5 +++ .../ComputeSurfaceFeaturesDirect.cpp | 7 +++- .../ComputeSurfaceFeaturesScanline.cpp | 9 +++++ 10 files changed, 92 insertions(+), 25 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.cpp index f198ff60b0..f7b20b1803 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckScanline.cpp @@ -27,6 +27,13 @@ BadDataNeighborOrientationCheckScanline::BadDataNeighborOrientationCheckScanline BadDataNeighborOrientationCheckScanline::~BadDataNeighborOrientationCheckScanline() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Flips bad voxels to good using chunk-sequential multi-pass scanning. + * OOC path: Phase 1 counts matching good neighbors using chunk-sequential iteration. + * Phase 2 uses repeated chunk-sequential scans (instead of a worklist) to flip + * eligible voxels, with a chunk-skip optimization that avoids loading chunks + * with no eligible voxels. + */ Result<> BadDataNeighborOrientationCheckScanline::operator()() { const float misorientationTolerance = m_InputValues->MisorientationTolerance * numbers::pi_v / 180.0f; diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.cpp index f3b8702b7f..82795f05e4 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/BadDataNeighborOrientationCheckWorklist.cpp @@ -29,6 +29,12 @@ BadDataNeighborOrientationCheckWorklist::BadDataNeighborOrientationCheckWorklist BadDataNeighborOrientationCheckWorklist::~BadDataNeighborOrientationCheckWorklist() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Flips bad voxels to good using a worklist-based propagation algorithm. + * In-core path: Phase 1 counts matching good neighbors per bad voxel. + * Phase 2 iteratively flips eligible voxels using a deque worklist, + * propagating new eligibility to neighbors immediately. + */ Result<> BadDataNeighborOrientationCheckWorklist::operator()() { const float misorientationTolerance = m_InputValues->MisorientationTolerance * numbers::pi_v / 180.0f; diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp index edde51e412..efedbd6348 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsDirect.cpp @@ -22,6 +22,10 @@ ComputeBoundaryCellsDirect::ComputeBoundaryCellsDirect(DataStructure& dataStruct ComputeBoundaryCellsDirect::~ComputeBoundaryCellsDirect() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Counts boundary faces per voxel using direct Z-Y-X iteration. + * In-core path: iterates all voxels sequentially, checking 6 face neighbors. + */ Result<> ComputeBoundaryCellsDirect::operator()() { const auto& imageGeometry = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp index 9a3e5f783e..25a5e90812 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryCellsScanline.cpp @@ -22,6 +22,11 @@ ComputeBoundaryCellsScanline::ComputeBoundaryCellsScanline(DataStructure& dataSt ComputeBoundaryCellsScanline::~ComputeBoundaryCellsScanline() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Counts boundary faces per voxel using chunk-sequential iteration. + * OOC path: iterates chunks in order via loadChunk/getChunkLowerBounds/getChunkUpperBounds, + * then Z-Y-X within each chunk. Same logic as ComputeBoundaryCellsDirect. + */ Result<> ComputeBoundaryCellsScanline::operator()() { const auto& imageGeometry = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp index da14331911..f3fb5d5e0f 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp @@ -9,8 +9,6 @@ #include "simplnx/Utilities/MessageHelper.hpp" #include "simplnx/Utilities/NeighborUtilities.hpp" -#include - using namespace nx::core; // ----------------------------------------------------------------------------- @@ -27,6 +25,11 @@ ComputeFeatureNeighborsDirect::ComputeFeatureNeighborsDirect(DataStructure& data ComputeFeatureNeighborsDirect::~ComputeFeatureNeighborsDirect() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Computes feature neighbor lists using direct sequential iteration. + * In-core path: Phase 1 iterates all voxels to build per-feature neighbor + * counts. Phase 2 deduplicates and computes shared surface areas. + */ Result<> ComputeFeatureNeighborsDirect::operator()() { auto storeBoundaryCells = m_InputValues->StoreBoundaryCells; @@ -58,16 +61,10 @@ Result<> ComputeFeatureNeighborsDirect::operator()() usize totalPoints = featureIds.getNumberOfTuples(); usize totalFeatures = numNeighbors.getNumberOfTuples(); - /* Ensure that we will be able to work with the user selected featureId Array */ - const auto [minFeatureId, maxFeatureId] = std::minmax_element(featureIds.begin(), featureIds.end()); - if(static_cast(*maxFeatureId) >= totalFeatures) - { - std::stringstream out; - out << "Data Array " << featureIdsPath.getTargetName() << " has a maximum value of " << *maxFeatureId << " which is greater than the " - << " number of features from array " << numNeighborsPath.getTargetName() << " which has " << totalFeatures << ". Did you select the " - << " incorrect array for the 'FeatureIds' array?"; - return MakeErrorResult(-24500, out.str()); - } + // Max feature ID validation is deferred to after the main loop to avoid + // a separate full scan through the data store. The loop's + // `feature < neighborlist.size()` guard prevents out-of-bounds access. + int32 observedMaxFeatureId = 0; auto& imageGeom = m_DataStructure.getDataRefAs(imageGeomPath); SizeVec3 uDims = imageGeom.getDimensions(); @@ -123,6 +120,10 @@ Result<> ComputeFeatureNeighborsDirect::operator()() onsurf = 0; feature = featureIds[voxelIndex]; + if(feature > observedMaxFeatureId) + { + observedMaxFeatureId = feature; + } if(feature > 0 && feature < neighborlist.size()) { int64 xIdx = voxelIndex % dims[0]; @@ -168,6 +169,15 @@ Result<> ComputeFeatureNeighborsDirect::operator()() } } + // Validate max feature ID (deferred from before the loop to avoid a separate full scan) + if(static_cast(observedMaxFeatureId) >= totalFeatures) + { + return MakeErrorResult(-24500, + fmt::format("Data Array {} has a maximum value of {} which is greater than the number of features from array {} which has {}. " + "Did you select the incorrect array for the 'FeatureIds' array?", + featureIdsPath.getTargetName(), observedMaxFeatureId, numNeighborsPath.getTargetName(), totalFeatures)); + } + FloatVec3 spacing = imageGeom.getSpacing(); // We do this to create new set of NeighborList objects diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp index c081fed686..44376308ab 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp @@ -9,8 +9,6 @@ #include "simplnx/Utilities/MessageHelper.hpp" #include "simplnx/Utilities/NeighborUtilities.hpp" -#include - using namespace nx::core; // ----------------------------------------------------------------------------- @@ -27,6 +25,12 @@ ComputeFeatureNeighborsScanline::ComputeFeatureNeighborsScanline(DataStructure& ComputeFeatureNeighborsScanline::~ComputeFeatureNeighborsScanline() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Computes feature neighbor lists using chunk-sequential iteration. + * OOC path: Phase 1 iterates chunks in order to build per-feature neighbor + * counts. Phase 2 deduplicates and computes shared surface areas. + * Same logic as ComputeFeatureNeighborsDirect. + */ Result<> ComputeFeatureNeighborsScanline::operator()() { auto storeBoundaryCells = m_InputValues->StoreBoundaryCells; @@ -58,16 +62,10 @@ Result<> ComputeFeatureNeighborsScanline::operator()() usize totalPoints = featureIds.getNumberOfTuples(); usize totalFeatures = numNeighbors.getNumberOfTuples(); - /* Ensure that we will be able to work with the user selected featureId Array */ - const auto [minFeatureId, maxFeatureId] = std::minmax_element(featureIds.begin(), featureIds.end()); - if(static_cast(*maxFeatureId) >= totalFeatures) - { - std::stringstream out; - out << "Data Array " << featureIdsPath.getTargetName() << " has a maximum value of " << *maxFeatureId << " which is greater than the " - << " number of features from array " << numNeighborsPath.getTargetName() << " which has " << totalFeatures << ". Did you select the " - << " incorrect array for the 'FeatureIds' array?"; - return MakeErrorResult(-24500, out.str()); - } + // Max feature ID validation is deferred to after the chunk-sequential loop + // to avoid a separate full scan through OOC data. The loop's + // `feature < neighborlist.size()` guard prevents out-of-bounds access. + int32 observedMaxFeatureId = 0; auto& imageGeom = m_DataStructure.getDataRefAs(imageGeomPath); SizeVec3 uDims = imageGeom.getDimensions(); @@ -140,6 +138,10 @@ Result<> ComputeFeatureNeighborsScanline::operator()() const int64 voxelIndex = kStride + jStride + static_cast(xIdx); onsurf = 0; feature = featureIds[voxelIndex]; + if(feature > observedMaxFeatureId) + { + observedMaxFeatureId = feature; + } if(feature > 0 && feature < neighborlist.size()) { const int64 ix = static_cast(xIdx); @@ -188,6 +190,15 @@ Result<> ComputeFeatureNeighborsScanline::operator()() } } + // Validate max feature ID (deferred from before the loop to avoid a separate OOC scan) + if(static_cast(observedMaxFeatureId) >= totalFeatures) + { + return MakeErrorResult(-24500, + fmt::format("Data Array {} has a maximum value of {} which is greater than the number of features from array {} which has {}. " + "Did you select the incorrect array for the 'FeatureIds' array?", + featureIdsPath.getTargetName(), observedMaxFeatureId, numNeighborsPath.getTargetName(), totalFeatures)); + } + // Phase 2: Build NeighborList objects (operates on small feature-level data, no chunk iteration needed) FloatVec3 spacing = imageGeom.getSpacing(); diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp index a131959ff2..97bb4ddb4a 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp @@ -24,6 +24,11 @@ ComputeSurfaceAreaToVolumeDirect::ComputeSurfaceAreaToVolumeDirect(DataStructure ComputeSurfaceAreaToVolumeDirect::~ComputeSurfaceAreaToVolumeDirect() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Computes surface-area-to-volume ratio using direct Z-Y-X iteration. + * In-core path: accumulates per-feature surface area from face-neighbor + * comparisons, then divides by voxel volume. Optionally computes sphericity. + */ Result<> ComputeSurfaceAreaToVolumeDirect::operator()() { // Input Cell Data diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp index 2a2105ae87..f442854a20 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp @@ -24,6 +24,11 @@ ComputeSurfaceAreaToVolumeScanline::ComputeSurfaceAreaToVolumeScanline(DataStruc ComputeSurfaceAreaToVolumeScanline::~ComputeSurfaceAreaToVolumeScanline() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Computes surface-area-to-volume ratio using chunk-sequential iteration. + * OOC path: iterates chunks in order via loadChunk/getChunkLowerBounds/getChunkUpperBounds. + * Same logic as ComputeSurfaceAreaToVolumeDirect. + */ Result<> ComputeSurfaceAreaToVolumeScanline::operator()() { // Input Cell Data diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.cpp index ba5def33eb..0a0f98ca13 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.cpp @@ -196,6 +196,11 @@ ComputeSurfaceFeaturesDirect::ComputeSurfaceFeaturesDirect(DataStructure& dataSt ComputeSurfaceFeaturesDirect::~ComputeSurfaceFeaturesDirect() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Identifies surface features using direct Z-Y-X iteration. + * In-core path: delegates to findSurfaceFeatures3D or findSurfaceFeatures2D + * depending on geometry dimensionality. + */ Result<> ComputeSurfaceFeaturesDirect::operator()() { const auto pMarkFeature0NeighborsValue = m_InputValues->MarkFeature0Neighbors; @@ -223,7 +228,7 @@ Result<> ComputeSurfaceFeaturesDirect::operator()() } else { - MakeErrorResult(-1000, fmt::format("Image Geometry at path '{}' must be either 3D or 2D", pFeatureGeometryPathValue.toString())); + return MakeErrorResult(-1000, fmt::format("Image Geometry at path '{}' must be either 3D or 2D", pFeatureGeometryPathValue.toString())); } return {}; diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.cpp index 306c955c7a..ace8e906af 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.cpp @@ -110,6 +110,11 @@ ComputeSurfaceFeaturesScanline::ComputeSurfaceFeaturesScanline(DataStructure& da ComputeSurfaceFeaturesScanline::~ComputeSurfaceFeaturesScanline() noexcept = default; // ----------------------------------------------------------------------------- +/** + * @brief Identifies surface features using chunk-sequential iteration. + * OOC path: iterates chunks in order, handling both 3D and 2D geometries + * with coordinate remapping. Same logic as ComputeSurfaceFeaturesDirect. + */ Result<> ComputeSurfaceFeaturesScanline::operator()() { const auto pMarkFeature0NeighborsValue = m_InputValues->MarkFeature0Neighbors; @@ -216,6 +221,10 @@ Result<> ComputeSurfaceFeaturesScanline::operator()() surfaceFeatures[gNum] = 1; } } + else + { + return MakeErrorResult(-1000, fmt::format("Image Geometry at path '{}' must be either 3D or 2D", pFeatureGeometryPathValue.toString())); + } } } } From 6f45328fb824d05b7584533b17dbbea1d374777d Mon Sep 17 00:00:00 2001 From: Joey Kleingers Date: Thu, 5 Mar 2026 22:39:33 -0500 Subject: [PATCH 5/5] STYLE: Run clang-format on Group B PR files --- .../Algorithms/ComputeFeatureNeighborsDirect.cpp | 7 +++---- .../Algorithms/ComputeFeatureNeighborsScanline.cpp | 10 ++++------ .../Algorithms/ComputeSurfaceFeaturesDirect.hpp | 3 +-- .../Algorithms/ComputeSurfaceFeaturesScanline.hpp | 3 +-- 4 files changed, 9 insertions(+), 14 deletions(-) diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp index f3fb5d5e0f..1db4813dc5 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsDirect.cpp @@ -172,10 +172,9 @@ Result<> ComputeFeatureNeighborsDirect::operator()() // Validate max feature ID (deferred from before the loop to avoid a separate full scan) if(static_cast(observedMaxFeatureId) >= totalFeatures) { - return MakeErrorResult(-24500, - fmt::format("Data Array {} has a maximum value of {} which is greater than the number of features from array {} which has {}. " - "Did you select the incorrect array for the 'FeatureIds' array?", - featureIdsPath.getTargetName(), observedMaxFeatureId, numNeighborsPath.getTargetName(), totalFeatures)); + return MakeErrorResult(-24500, fmt::format("Data Array {} has a maximum value of {} which is greater than the number of features from array {} which has {}. " + "Did you select the incorrect array for the 'FeatureIds' array?", + featureIdsPath.getTargetName(), observedMaxFeatureId, numNeighborsPath.getTargetName(), totalFeatures)); } FloatVec3 spacing = imageGeom.getSpacing(); diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp index 44376308ab..bdd454bc58 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureNeighborsScanline.cpp @@ -127,8 +127,7 @@ Result<> ComputeFeatureNeighborsScanline::operator()() const int64 jStride = static_cast(yIdx) * dims[0]; for(usize xIdx = chunkLowerBounds[2]; xIdx <= chunkUpperBounds[2]; xIdx++) { - throttledMessenger.sendThrottledMessage( - [&]() { return fmt::format("Determining Neighbor Lists || Chunk {}/{}", chunkIdx + 1, numChunks); }); + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Determining Neighbor Lists || Chunk {}/{}", chunkIdx + 1, numChunks); }); if(m_ShouldCancel) { @@ -193,10 +192,9 @@ Result<> ComputeFeatureNeighborsScanline::operator()() // Validate max feature ID (deferred from before the loop to avoid a separate OOC scan) if(static_cast(observedMaxFeatureId) >= totalFeatures) { - return MakeErrorResult(-24500, - fmt::format("Data Array {} has a maximum value of {} which is greater than the number of features from array {} which has {}. " - "Did you select the incorrect array for the 'FeatureIds' array?", - featureIdsPath.getTargetName(), observedMaxFeatureId, numNeighborsPath.getTargetName(), totalFeatures)); + return MakeErrorResult(-24500, fmt::format("Data Array {} has a maximum value of {} which is greater than the number of features from array {} which has {}. " + "Did you select the incorrect array for the 'FeatureIds' array?", + featureIdsPath.getTargetName(), observedMaxFeatureId, numNeighborsPath.getTargetName(), totalFeatures)); } // Phase 2: Build NeighborList objects (operates on small feature-level data, no chunk iteration needed) diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp index 172a29a51d..ebacffdb14 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesDirect.hpp @@ -18,8 +18,7 @@ struct ComputeSurfaceFeaturesInputValues; class SIMPLNXCORE_EXPORT ComputeSurfaceFeaturesDirect { public: - ComputeSurfaceFeaturesDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, - const ComputeSurfaceFeaturesInputValues* inputValues); + ComputeSurfaceFeaturesDirect(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, const ComputeSurfaceFeaturesInputValues* inputValues); ~ComputeSurfaceFeaturesDirect() noexcept; ComputeSurfaceFeaturesDirect(const ComputeSurfaceFeaturesDirect&) = delete; diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp index d3be471d12..3e757e7c01 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceFeaturesScanline.hpp @@ -19,8 +19,7 @@ struct ComputeSurfaceFeaturesInputValues; class SIMPLNXCORE_EXPORT ComputeSurfaceFeaturesScanline { public: - ComputeSurfaceFeaturesScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, - const ComputeSurfaceFeaturesInputValues* inputValues); + ComputeSurfaceFeaturesScanline(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, const ComputeSurfaceFeaturesInputValues* inputValues); ~ComputeSurfaceFeaturesScanline() noexcept; ComputeSurfaceFeaturesScanline(const ComputeSurfaceFeaturesScanline&) = delete;