diff --git a/src/Plugins/SimplnxCore/docs/ComputeSurfaceAreaToVolumeFilter.md b/src/Plugins/SimplnxCore/docs/ComputeSurfaceAreaToVolumeFilter.md index efe9c63b12..36ff9b9418 100644 --- a/src/Plugins/SimplnxCore/docs/ComputeSurfaceAreaToVolumeFilter.md +++ b/src/Plugins/SimplnxCore/docs/ComputeSurfaceAreaToVolumeFilter.md @@ -12,32 +12,78 @@ This filter also optionally calculates the [Sphericity](https://en.wikipedia.org ![Equation for Sphericity used in the filter](Images/Sphericity_Equation.png) -This **Filter** determines whether a **Feature** touches an outer *Surface* of the sample volume. A **Feature** is considered touching the *Surface* of the sample if either of the following conditions are met: - -+ Any cell location is xmin, xmax, ymin, ymax, zmin or zmax -+ Any cell has **Feature ID = 0** as a neighbor. +Use [Compute Surface Features](ComputeSurfaceFeaturesFilter.md) when you need a separate flag that identifies features touching the sample surface. ### Required Input Sources - **Cell Feature Ids** -- the per-**Cell** **Feature** label array, typically produced by [Segment Features (Scalar)](ScalarSegmentFeaturesFilter.md) or another segmentation filter. +- **Number of Cells** -- the per-feature cell counts, typically produced by [Compute Feature Sizes](ComputeFeatureSizesFilter.md). + +The label array must have one tuple per cell in the selected geometry. The array can be stored separately from that geometry. Its tuple count, cell ordering, and feature labels must describe the selected image. You are responsible for checking that the cell ordering and feature labels match the selected image. The tuple-count check cannot verify this correspondence. Each Feature Id must index a tuple in the selected feature data; the filter checks these index bounds. + +## Algorithm + +For each cell with a positive Feature Id, the filter examines its six face-sharing neighbors within the image. A face contributes area to the current feature when the neighbor has a different Feature Id, including Feature Id 0. An interface between two positive features contributes to each feature's area. Faces outside the image are omitted. Feature Id 0 is the background; its output tuple remains zero. + +For cell spacing `(dx, dy, dz)`, the area of a face uses the two spacing components perpendicular to its normal: + +- X-normal face: `dy * dz`. +- Y-normal face: `dx * dz`. +- Z-normal face: `dx * dy`. + +![Fig. 1: Generated schematic showing all six voxel faces and the two spacing components used for each face area.](Images/ComputeSurfaceAreaToVolume_VoxelFaceAreas.png) + +The filter calculates each feature's volume as `V = Number of Cells * dx * dy * dz`. It uses the supplied cell count without recounting the labels or checking that the count matches them. The surface-area-to-volume ratio is `S / V`, where `S` is the accumulated face area. Area has units of length squared, and volume has units of length cubed. + +When *Calculate Sphericity* is enabled, the filter also calculates `Psi = pi^(1/3) * (6 * V)^(2/3) / S`. Psi is dimensionless. A mathematical sphere with a complete closed surface and consistent positive volume has sphericity 1. For a voxelized sphere with a complete counted surface and consistent positive volume, the stair-stepped surface has more area than an equal-volume smooth sphere, so its sphericity is below 1. Both output arrays use 32-bit floating-point values; face areas, accumulated areas, and metric calculations use 64-bit floating-point values. + +### Algorithm selection + +The filter automatically uses *Direct* when Cell Feature Ids are stored in memory and *Scanline* when they use out-of-core storage. Both methods use the same face-counting and volume rules. + +### Direct + +Direct reads neighboring labels from the array and accumulates one surface area per feature. Its additional working memory scales as `O(F)`, where `F` is the number of feature tuples, including background. -## Algorithm Details +### Scanline -- First, all the boundary **Cells** are found for each **Feature** (see also [Compute Boundary Cells (Image)](ComputeBoundaryCellsFilter.md)). -- Next, the surface area for each face that is in contact with a different **Feature** is totalled as long as that neighboring *featureId* is > 0. -- This number is divided by the volume of each **Feature**, calculated by taking the number of **Cells** of each **Feature** and multiplying by the volume of a **Cell**. +Scanline reads labels in sequential slices and keeps three adjacent XY slices in memory. It releases these slices before reading feature counts and writing metrics in bounded batches. Each batch has at most `min(65536, max(1, floor(F / 2)))` tuples, where `F` includes background. The count and result buffers together use at most 512 KiB. This is an implementation memory budget, not a limit on the number of features. + +The 64-bit area accumulator still uses eight bytes per feature tuple. Additional working memory scales as `O(X * Y + F)`, where `X` and `Y` are image dimensions in cells. Scanline avoids buffering the full label volume, but the area accumulator can still be large when the feature count approaches the cell count. These memory estimates exclude the input/output stores, their caches, and allocator overhead. + +A bulk read or write error stops execution. If a feature-count read or metric write fails after earlier batches completed, those earlier outputs remain written. Output updates are not transactional. + +### Compatibility with earlier results + +Corrected X/Y face weights can change results when `dx` differs from `dy`. The face-weight correction does not change the measured area when `dx = dy`, including when `dz` differs. Changes to numerical precision and the sphericity exponents can still change results, including isotropic results. + +### Small IN100 comparison + +The generated views below use segmented Small IN100 feature labels in a separate demonstration copy. The original spacing is `(0.25, 0.25, 0.25)` micrometers. The demonstration changes the spacing to `(0.25, 0.50, 0.25)` micrometers so that X-normal and Y-normal faces have different areas. The original dataset is unchanged, and the labels are reused without resampling. Both calculations use identical labels, cell counts derived from those labels, and demonstration spacing. + +The earlier Direct implementation and the current filter each calculate sphericity over the complete three-dimensional volume. The views display those three-dimensional feature values on the same central XY slice, at zero-based Z index 58. Both views use the same dimensionless 0–1 Viridis color scale and physical aspect ratio. The comparison includes corrections to numerical precision and sphericity exponents as well as the X/Y face weights. Features cut by the image boundary can have incomplete counted surfaces, as described below. + +![Small IN100 sphericity before and after the correction](Images/ComputeSurfaceAreaToVolume_SmallIN100_Sphericity.png) + +Fig. 2: Generated before-and-after views of full-volume Small IN100 feature sphericity on a central XY slice. Spacing was altered from 0.25/0.25/0.25 to 0.25/0.50/0.25 micrometers for demonstration; both panels share the 0–1 Viridis color scale and physical aspect ratio. ### WARNING - Aliasing The surface area will be the surface area of the **Cells** in contact with the neighboring **Feature** and will be influenced by the *aliasing* of the structure — the jagged, stair-stepped approximation of a smooth surface that results from representing it on a discrete voxel grid. As a result, the surface area to volume will likely be over-estimated with respect to the *real* **Feature**. -### WARNING - Skewed Results for features touching the surface +### WARNING - Incomplete surfaces and nonfinite results -Because the filter does not include any surface that is touching/connected to a "FeatureId = 0", those features that are in contact with the edge of the virtual volume or in contact with internal features that are labeled as "FeatureId = 0" will have their values skewed. +Image-exterior faces are omitted, so a feature touching the image boundary can have an incomplete measured surface. This can reduce its surface-area-to-volume ratio and increase its sphericity, even above 1. In-image faces against Feature Id 0 are included. + +The filter does not reject zero or negative supplied cell counts or replace nonfinite results. For positive area and zero volume, the ratio is infinite and sphericity is zero. For positive volume and zero measured area, the ratio is zero and sphericity is infinite. Zero area with zero volume produces NaN (not a number). A negative volume can produce a negative ratio and NaN sphericity. + +Zero or negative supplied counts do not give a positive feature volume. Check **Number of Cells** against the feature labels before interpreting either metric. A zero measured area does not support a finite sphericity value, but it can result from omitted image-exterior faces even when the labels and counts are correct. ### Warning - 2D Image Geometry Results -Because even a single slice has *volume* according to DREAM3D-NX, results will still be computed. These results should **NOT** be interpreted as "Boundary Length to Area" values. +For a single XY slice, each cell still has thickness `dz`, and the filter uses three-dimensional cell volumes. Both Z-normal faces lie outside the image and are omitted. The counted side-face area equals the counted in-plane boundary length times `dz`. + +When the supplied counts match the labels, thickness cancels in the ratio. The result then equals counted in-plane boundary length divided by feature area, with units of inverse length. Image-exterior edges remain omitted, so this need not be the feature's full perimeter-to-area ratio. Optional sphericity still uses the three-dimensional formula and depends on the slice thickness; it is not a two-dimensional circularity measure. % Auto generated parameter table will be inserted here diff --git a/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_SmallIN100_Sphericity.png b/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_SmallIN100_Sphericity.png new file mode 100644 index 0000000000..6ac15e3d40 Binary files /dev/null and b/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_SmallIN100_Sphericity.png differ diff --git a/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_SmallIN100_Sphericity.svg b/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_SmallIN100_Sphericity.svg new file mode 100644 index 0000000000..029d214fcc --- /dev/null +++ b/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_SmallIN100_Sphericity.svg @@ -0,0 +1,38 @@ + + + + +Small IN100 demonstration: Sphericity +Generated views • full 3D calculation; central XY slice shows each feature’s 3D metric + +Altered spacing for demonstration: dx = 0.25, dy = 0.50, dz = 0.25 µm. +Original spacing: 0.25, 0.25, 0.25 µm. Original file and labels are unchanged. +Before +Historical Direct • 5aaa879debf9 + + +X extent 47.25 µm; Y extent 100.50 µm +After +Current filter • corrected face weights + + +X extent 47.25 µm; Y extent 100.50 µm +Shared sphericity scale • dimensionless • Viridis + + +0 + +0.25 + +0.5 + +0.75 + +1 + +Background • Feature Id 0 +No scale clipping • Y increases upward • Z index 58 of 0–116 +Whole volume: 2,317 populated positive features; mean |change| 0.028309; max |change| 0.19037 dimensionless. +Both views use identical labels, counts, spacing, physical aspect and full-range color limits. +Comparison includes precision/exponent changes. Out-of-image faces are excluded; edge features may have incomplete surfaces. + \ No newline at end of file diff --git a/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_VoxelFaceAreas.png b/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_VoxelFaceAreas.png new file mode 100644 index 0000000000..b2d3848b4d Binary files /dev/null and b/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_VoxelFaceAreas.png differ diff --git a/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_VoxelFaceAreas.svg b/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_VoxelFaceAreas.svg new file mode 100644 index 0000000000..f318cf6fd3 --- /dev/null +++ b/src/Plugins/SimplnxCore/docs/Images/ComputeSurfaceAreaToVolume_VoxelFaceAreas.svg @@ -0,0 +1,53 @@ + + + + + +Voxel face areas: the face normal determines the two side lengths +Generated schematic • opposite views expose all six faces • lengths are illustrative +View toward the +X, +Y and +Z faces +Opposite view: −X, −Y and −Z faces + ++X + ++Y + ++Z +dx +dy +dz + +−X + +−Y + +−Z +dx +dy +dz + +Face normal +Face area +Demonstration spacing (µm) +Area (µm²) + + +−X and +X +dy × dz +0.50 × 0.25 +0.125 + + +−Y and +Y +dx × dz +0.25 × 0.25 +0.0625 + + +−Z and +Z +dx × dy +0.25 × 0.50 +0.125 +Count a face when its in-bounds neighbor has a different Feature Id. Out-of-image faces are excluded. +Anisotropic Small IN100 demonstration: dx = 0.25, dy = 0.50, dz = 0.25 µm. Original data spacing remains unchanged. + \ No newline at end of file diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp index 0802643fd8..91aedaf4f9 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp @@ -1,10 +1,13 @@ -#include - #include "ComputeSurfaceAreaToVolumeDirect.hpp" +#include +#include +#include + #include "ComputeSurfaceAreaToVolume.hpp" #include "simplnx/Common/Constants.hpp" +#include "simplnx/Common/Types.hpp" #include "simplnx/DataStructure/DataArray.hpp" #include "simplnx/DataStructure/DataGroup.hpp" #include "simplnx/DataStructure/Geometry/ImageGeom.hpp" @@ -39,7 +42,7 @@ Result<> ComputeSurfaceAreaToVolumeDirect::operator()() { return validateNumFeatResult; } - auto numFeatures = static_cast(numCells.getNumberOfTuples()); + const usize numFeatures = numCells.getNumberOfTuples(); SizeVec3 dims = imageGeom.getDimensions(); FloatVec3 spacing = imageGeom.getSpacing(); @@ -47,11 +50,14 @@ Result<> ComputeSurfaceAreaToVolumeDirect::operator()() auto yPoints = static_cast(dims[1]); auto zPoints = static_cast(dims[2]); - float32 voxelVol = spacing[0] * spacing[1] * spacing[2]; + const float64 voxelVol = static_cast(spacing[0]) * spacing[1] * spacing[2]; + const float64 xyFaceArea = static_cast(spacing[0]) * spacing[1]; + const float64 yzFaceArea = static_cast(spacing[1]) * spacing[2]; + const float64 zxFaceArea = static_cast(spacing[2]) * spacing[0]; // Multiple cells contribute to one feature, so accumulation stays local until // the cell scan is complete. - std::vector featureSurfaceArea(static_cast(numFeatures), 0.0f); + std::vector featureSurfaceArea(numFeatures, 0.0); // Flat offsets avoid repeated coordinate-to-index conversion for in-memory // six-neighbor reads. @@ -75,7 +81,7 @@ Result<> ComputeSurfaceAreaToVolumeDirect::operator()() int64 yStride = yIdx * xPoints; for(int64 xIdx = 0; xIdx < xPoints; xIdx++) { - float32 onSurface = 0.0f; + float64 onSurface = 0.0; int32 currentFeatureId = featureIdsStoreRef[zStride + yStride + xIdx]; // Feature 0 is the background tuple. if(currentFeatureId < 1) @@ -116,41 +122,41 @@ Result<> ComputeSurfaceAreaToVolumeDirect::operator()() { if(neighborOffsetIndex == 0 || neighborOffsetIndex == 5) { - onSurface = onSurface + spacing[0] * spacing[1]; + onSurface += xyFaceArea; } if(neighborOffsetIndex == 1 || neighborOffsetIndex == 4) { - onSurface = onSurface + spacing[1] * spacing[2]; + onSurface += zxFaceArea; } if(neighborOffsetIndex == 2 || neighborOffsetIndex == 3) { - onSurface = onSurface + spacing[2] * spacing[0]; + onSurface += yzFaceArea; } } } - int32 featureId = featureIdsStoreRef[zStride + yStride + xIdx]; - featureSurfaceArea[featureId] = featureSurfaceArea[featureId] + onSurface; + featureSurfaceArea[currentFeatureId] += onSurface; } } } - const float32 thirdRootPi = std::pow(nx::core::Constants::k_PiF, 0.333333f); for(usize i = 1; i < numFeatures; i++) { - float32 featureVolume = voxelVol * numCells[i]; - surfaceAreaVolumeRatio[i] = featureSurfaceArea[i] / featureVolume; + const float64 featureVolume = voxelVol * numCells[i]; + surfaceAreaVolumeRatio[i] = static_cast(featureSurfaceArea[i] / featureVolume); } if(m_InputValues->CalculateSphericity) { m_MessageHandler.sendMessage(IFilter::Message::Type::Info, fmt::format("Computing Sphericity")); + const float64 thirdRootPi = std::cbrt(nx::core::Constants::k_PiD); + auto& sphericity = m_DataStructure.getDataAs(m_InputValues->SphericityArrayName)->getDataStoreRef(); // Sphericity is one for a sphere and decreases as the surface becomes less compact. - for(usize i = 1; i < static_cast(numFeatures); i++) + for(usize i = 1; i < numFeatures; i++) { - float32 featureVolume = voxelVol * numCells[i]; - sphericity[i] = (thirdRootPi * std::pow((6.0f * featureVolume), 0.66666f)) / featureSurfaceArea[i]; + const float64 featureVolume = voxelVol * numCells[i]; + sphericity[i] = static_cast((thirdRootPi * std::pow((6.0 * featureVolume), 2.0 / 3.0)) / featureSurfaceArea[i]); } } diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp index e982b8daa2..e0cc80ae7c 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.cpp @@ -1,8 +1,14 @@ #include "ComputeSurfaceAreaToVolumeScanline.hpp" +#include +#include +#include + #include "ComputeSurfaceAreaToVolume.hpp" #include "simplnx/Common/Constants.hpp" +#include "simplnx/Common/Types.hpp" +#include "simplnx/DataStructure/AbstractDataStore.hpp" #include "simplnx/DataStructure/DataArray.hpp" #include "simplnx/DataStructure/DataGroup.hpp" #include "simplnx/DataStructure/Geometry/ImageGeom.hpp" @@ -27,7 +33,7 @@ ComputeSurfaceAreaToVolumeScanline::~ComputeSurfaceAreaToVolumeScanline() noexce Result<> ComputeSurfaceAreaToVolumeScanline::operator()() { auto featureIdsArrayPtr = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath); - auto& featureIdsStore = featureIdsArrayPtr->getDataStoreRef(); + const auto& featureIdsStore = featureIdsArrayPtr->getDataStoreRef(); const auto& numCells = m_DataStructure.getDataAs(m_InputValues->NumCellsArrayPath)->getDataStoreRef(); auto& surfaceAreaVolumeRatio = m_DataStructure.getDataAs(m_InputValues->SurfaceAreaVolumeRatioArrayName)->getDataStoreRef(); @@ -38,7 +44,11 @@ Result<> ComputeSurfaceAreaToVolumeScanline::operator()() { return validateNumFeatResult; } - auto numFeatures = static_cast(numCells.getNumberOfTuples()); + const usize numFeatures = numCells.getNumberOfTuples(); + if(numFeatures == 0) + { + return {}; + } SizeVec3 dims = imageGeom.getDimensions(); FloatVec3 spacing = imageGeom.getSpacing(); @@ -46,163 +56,178 @@ Result<> ComputeSurfaceAreaToVolumeScanline::operator()() auto yPoints = static_cast(dims[1]); auto zPoints = static_cast(dims[2]); - float32 voxelVol = spacing[0] * spacing[1] * spacing[2]; + const float64 voxelVol = static_cast(spacing[0]) * spacing[1] * spacing[2]; // Multiple cells contribute to one feature, so accumulation stays local until // the cell scan is complete. - std::vector featureSurfaceArea(static_cast(numFeatures), 0.0f); + std::vector featureSurfaceArea(numFeatures, 0.0); - const float32 xyFaceArea = spacing[0] * spacing[1]; - const float32 yzFaceArea = spacing[1] * spacing[2]; - const float32 zxFaceArea = spacing[2] * spacing[0]; + const float64 xyFaceArea = static_cast(spacing[0]) * spacing[1]; + const float64 yzFaceArea = static_cast(spacing[1]) * spacing[2]; + const float64 zxFaceArea = static_cast(spacing[2]) * spacing[0]; const usize sliceSize = static_cast(yPoints) * static_cast(xPoints); - // The rolling window turns Z-neighbor reads into sequential slice I/O. - std::vector prevSlice(sliceSize, 0); - std::vector curSlice(sliceSize, 0); - std::vector nextSlice(sliceSize, 0); - - // A bulk-I/O error stops the scan before later slices or output arrays change. - Result<> ioResult = featureIdsStore.copyIntoBuffer(0, nonstd::span(curSlice.data(), sliceSize)); - if(ioResult.invalid()) { - return ioResult; - } - if(zPoints > 1) - { - ioResult = featureIdsStore.copyIntoBuffer(sliceSize, nonstd::span(nextSlice.data(), sliceSize)); + // The rolling window turns Z-neighbor reads into sequential slice I/O. + std::vector prevSlice(sliceSize, 0); + std::vector curSlice(sliceSize, 0); + std::vector nextSlice(sliceSize, 0); + + // A bulk-I/O error stops the scan before later slices or output arrays change. + Result<> ioResult = featureIdsStore.copyIntoBuffer(0, nonstd::span(curSlice.data(), sliceSize)); if(ioResult.invalid()) { return ioResult; } - } - - for(int64 z = 0; z < zPoints; z++) - { - if(m_ShouldCancel) + if(zPoints > 1) { - return {}; + ioResult = featureIdsStore.copyIntoBuffer(sliceSize, nonstd::span(nextSlice.data(), sliceSize)); + if(ioResult.invalid()) + { + return ioResult; + } } - for(int64 y = 0; y < yPoints; y++) + + for(int64 zIndex = 0; zIndex < zPoints; zIndex++) { - for(int64 x = 0; x < xPoints; x++) + if(m_ShouldCancel) { - const usize inSlice = static_cast(y) * static_cast(xPoints) + static_cast(x); - int32 currentFeatureId = curSlice[inSlice]; - if(currentFeatureId < 1) + return {}; + } + for(int64 yIndex = 0; yIndex < yPoints; yIndex++) + { + for(int64 xIndex = 0; xIndex < xPoints; xIndex++) { - continue; - } + const usize inSlice = (static_cast(yIndex) * static_cast(xPoints)) + static_cast(xIndex); + const int32 currentFeatureId = curSlice[inSlice]; + if(currentFeatureId < 1) + { + continue; + } - float32 onSurface = 0.0f; + float64 onSurface = 0.0; - // Z neighbors use adjacent buffers. X and Y neighbors use the current buffer. - if(z > 0) - { - if(prevSlice[inSlice] != currentFeatureId) + // Z neighbors use adjacent buffers. X and Y neighbors use the current buffer. + if(zIndex > 0) { - onSurface += xyFaceArea; + if(prevSlice[inSlice] != currentFeatureId) + { + onSurface += xyFaceArea; + } } - } - if(z < zPoints - 1) - { - if(nextSlice[inSlice] != currentFeatureId) + if(zIndex < zPoints - 1) { - onSurface += xyFaceArea; + if(nextSlice[inSlice] != currentFeatureId) + { + onSurface += xyFaceArea; + } } - } - if(y > 0) - { - if(curSlice[inSlice - static_cast(xPoints)] != currentFeatureId) + if(yIndex > 0) { - onSurface += yzFaceArea; + if(curSlice[inSlice - static_cast(xPoints)] != currentFeatureId) + { + onSurface += zxFaceArea; + } } - } - if(y < yPoints - 1) - { - if(curSlice[inSlice + static_cast(xPoints)] != currentFeatureId) + if(yIndex < yPoints - 1) { - onSurface += yzFaceArea; + if(curSlice[inSlice + static_cast(xPoints)] != currentFeatureId) + { + onSurface += zxFaceArea; + } } - } - if(x > 0) - { - if(curSlice[inSlice - 1] != currentFeatureId) + if(xIndex > 0) { - onSurface += zxFaceArea; + if(curSlice[inSlice - 1] != currentFeatureId) + { + onSurface += yzFaceArea; + } } - } - if(x < xPoints - 1) - { - if(curSlice[inSlice + 1] != currentFeatureId) + if(xIndex < xPoints - 1) { - onSurface += zxFaceArea; + if(curSlice[inSlice + 1] != currentFeatureId) + { + onSurface += yzFaceArea; + } } - } - featureSurfaceArea[currentFeatureId] += onSurface; + featureSurfaceArea[currentFeatureId] += onSurface; + } } - } - // Vector swaps rotate the window without copying a slice. - std::swap(prevSlice, curSlice); - std::swap(curSlice, nextSlice); - if(z + 2 < zPoints) - { - ioResult = featureIdsStore.copyIntoBuffer(static_cast(z + 2) * sliceSize, nonstd::span(nextSlice.data(), sliceSize)); - if(ioResult.invalid()) + // Vector swaps rotate the window without copying a slice. + std::swap(prevSlice, curSlice); + std::swap(curSlice, nextSlice); + if(zIndex + 2 < zPoints) { - return ioResult; + ioResult = featureIdsStore.copyIntoBuffer(static_cast(zIndex + 2) * sliceSize, nonstd::span(nextSlice.data(), sliceSize)); + if(ioResult.invalid()) + { + return ioResult; + } } } } - // Feature-level caches avoid per-value OOC access during metric calculation. - const usize numFeaturesUSize = static_cast(numFeatures); - std::vector localNumCells(numFeaturesUSize); - ioResult = numCells.copyIntoBuffer(0, nonstd::span(localNumCells.data(), numFeaturesUSize)); - if(ioResult.invalid()) - { - return ioResult; - } - - std::vector localSurfaceAreaVolumeRatio(numFeaturesUSize, 0.0f); - - const float32 thirdRootPi = std::pow(nx::core::Constants::k_PiF, 0.333333f); - for(usize i = 1; i < numFeaturesUSize; i++) - { - float32 featureVolume = voxelVol * localNumCells[i]; - localSurfaceAreaVolumeRatio[i] = featureSurfaceArea[i] / featureVolume; - } - ioResult = surfaceAreaVolumeRatio.copyFromBuffer(0, nonstd::span(localSurfaceAreaVolumeRatio.data(), numFeaturesUSize)); - if(ioResult.invalid()) - { - return ioResult; - } + // The slice buffers are released before these batches are allocated. The cap + // gives paired count/result buffers a 512 KiB budget, not a feature-count limit. + // For F >= 2, B <= F/2 also limits their combined payload to 4*F bytes. + constexpr usize maxBatchTuples = 65536; + const usize batchCapacity = std::min(maxBatchTuples, std::max(usize{1}, numFeatures / 2)); + std::vector localNumCells(batchCapacity); + std::vector localMetrics(batchCapacity); + AbstractDataStore* sphericityPtr = nullptr; if(m_InputValues->CalculateSphericity) { m_MessageHandler.sendMessage(IFilter::Message::Type::Info, fmt::format("Computing Sphericity")); + sphericityPtr = &m_DataStructure.getDataAs(m_InputValues->SphericityArrayName)->getDataStoreRef(); + } + const float64 thirdRootPi = m_InputValues->CalculateSphericity ? std::cbrt(nx::core::Constants::k_PiD) : 0.0; + for(usize offset = 0; offset < numFeatures;) + { + const usize batchSize = std::min(batchCapacity, numFeatures - offset); + auto ioResult = numCells.copyIntoBuffer(offset, nonstd::span(localNumCells.data(), batchSize)); + if(ioResult.invalid()) + { + return ioResult; + } - auto& sphericity = m_DataStructure.getDataAs(m_InputValues->SphericityArrayName)->getDataStoreRef(); - // Sphericity is one for a sphere and decreases as the surface becomes less compact. - std::vector localSphericity(numFeaturesUSize, 0.0f); - for(usize i = 1; i < numFeaturesUSize; i++) + // Background stays zero. Every other tuple uses the supplied cell count. + const usize firstFeature = offset == 0 ? 1 : 0; + localMetrics[0] = 0.0f; + for(usize batchIndex = firstFeature; batchIndex < batchSize; batchIndex++) { - float32 featureVolume = voxelVol * localNumCells[i]; - localSphericity[i] = (thirdRootPi * std::pow((6.0f * featureVolume), 0.66666f)) / featureSurfaceArea[i]; + const float64 featureVolume = voxelVol * localNumCells[batchIndex]; + localMetrics[batchIndex] = static_cast(featureSurfaceArea[offset + batchIndex] / featureVolume); } - ioResult = sphericity.copyFromBuffer(0, nonstd::span(localSphericity.data(), numFeaturesUSize)); + ioResult = surfaceAreaVolumeRatio.copyFromBuffer(offset, nonstd::span(localMetrics.data(), batchSize)); if(ioResult.invalid()) { return ioResult; } + + if(sphericityPtr != nullptr) + { + // Reuse the ratio buffer after its write; no second output buffer is needed. + for(usize batchIndex = firstFeature; batchIndex < batchSize; batchIndex++) + { + const float64 featureVolume = voxelVol * localNumCells[batchIndex]; + localMetrics[batchIndex] = static_cast((thirdRootPi * std::pow((6.0 * featureVolume), 2.0 / 3.0)) / featureSurfaceArea[offset + batchIndex]); + } + ioResult = sphericityPtr->copyFromBuffer(offset, nonstd::span(localMetrics.data(), batchSize)); + if(ioResult.invalid()) + { + return ioResult; + } + } + offset += batchSize; } return {}; diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp index 94bacfc293..003def673e 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp @@ -13,9 +13,13 @@ struct ComputeSurfaceAreaToVolumeInputValues; * @class ComputeSurfaceAreaToVolumeScanline * @brief Computes surface metrics with a three-slice Feature Id window. * - * The window converts disk-backed neighbor reads to sequential slice I/O. Memory - * is three cell-scale slices plus feature-scale caches. Feature caches can be - * large when the feature count approaches the cell count. + * The window converts disk-backed neighbor reads to sequential slice I/O. The + * three label slices are released before count and metric batches are allocated. + * Let F be the number of feature tuples and S the number of cells per XY slice. + * For F > 0, these buffers require 8*F + max(12*S, 8*B) bytes. + * The batch capacity is B = min(65536, max(1, floor(F/2))). + * The 8*F-byte area accumulator can still be large. This estimate excludes + * stores, caches, allocator overhead, and the separate Feature Id validation. * * @see ComputeSurfaceAreaToVolumeDirect for the in-core variant. * @see ComputeSurfaceAreaToVolume for the dispatcher. @@ -46,10 +50,11 @@ class SIMPLNXCORE_EXPORT ComputeSurfaceAreaToVolumeScanline /** * @brief Calculates surface metrics with bulk I/O. - * @return Error from Feature Id validation, or success after cancellation. + * @return Error from Feature Id validation or bulk I/O; success on completion or cancellation. * - * Cancellation during the cell scan leaves metric outputs unchanged. Current - * bulk-I/O Result values are not inspected. + * At each slice start, cancellation stops processing before metric writes. + * A bulk-I/O error stops processing. A later count read or metric write error + * can leave earlier batches written. Output updates are not transactional. */ Result<> operator()(); diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/ComputeSurfaceAreaToVolumeFilter.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/ComputeSurfaceAreaToVolumeFilter.cpp index 6da89db6e4..e3c4d65dc1 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/ComputeSurfaceAreaToVolumeFilter.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/ComputeSurfaceAreaToVolumeFilter.cpp @@ -2,18 +2,27 @@ #include "SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolume.hpp" +#include "simplnx/Common/Result.hpp" #include "simplnx/DataStructure/AttributeMatrix.hpp" #include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/DataObject.hpp" #include "simplnx/DataStructure/DataPath.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" #include "simplnx/Filter/Actions/CreateArrayAction.hpp" +#include "simplnx/Filter/Output.hpp" #include "simplnx/Parameters/ArraySelectionParameter.hpp" #include "simplnx/Parameters/BoolParameter.hpp" #include "simplnx/Parameters/DataObjectNameParameter.hpp" +#include "simplnx/Utilities/DataArrayUtilities.hpp" #include "simplnx/Utilities/SIMPLConversion.hpp" #include "simplnx/Parameters/GeometrySelectionParameter.hpp" +#include + +#include + using namespace nx::core; namespace nx::core @@ -95,6 +104,24 @@ IFilter::PreflightResult ComputeSurfaceAreaToVolumeFilter::preflightImpl(const D auto pFeatureIdsArrayPathValue = filterArgs.value(k_CellFeatureIdsArrayPath_Key); auto pNumCellsArrayPathValue = filterArgs.value(k_NumCellsArrayPath_Key); auto pCalculateSphericityValue = filterArgs.value(k_CalculateSphericity_Key); + auto pImageGeometryPathValue = filterArgs.value(k_SelectedImageGeometryPath_Key); + + const auto& imageGeom = dataStructure.getDataRefAs(pImageGeometryPathValue); + // ArraySelectionParameter also accepts StringArray and NeighborList for numeric selections. + // Keep the local type checks until the parameter validates supported array kinds. + // TODO(#1833): Revisit these guards after the shared fix preserves generic array support. + const auto* featureIdsPtr = dataStructure.getDataAs(pFeatureIdsArrayPathValue); + if(featureIdsPtr == nullptr) + { + const auto& selectedObject = dataStructure.getDataRefAs(pFeatureIdsArrayPathValue); + return MakePreflightErrorResult(-12803, fmt::format("The selected Feature Ids object '{}' has type '{}', but a scalar int32 DataArray is required. Select a scalar int32 DataArray.", + pFeatureIdsArrayPathValue.toString(), selectedObject.getTypeName())); + } + auto tupleCountResult = ValidateCellTupleCount(imageGeom, pImageGeometryPathValue, featureIdsPtr->getNumberOfTuples(), pFeatureIdsArrayPathValue, -12804); + if(tupleCountResult.invalid()) + { + return {.outputActions = ConvertInvalidResult(std::move(tupleCountResult))}; + } PreflightResult preflightResult; diff --git a/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTest.cpp b/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTest.cpp index 8e54ad89f7..9710f541b8 100644 --- a/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTest.cpp +++ b/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTest.cpp @@ -1,18 +1,46 @@ +#include "SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolume.hpp" +#include "SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeScanline.hpp" #include "SimplnxCore/Filters/ComputeSurfaceAreaToVolumeFilter.hpp" #include "SimplnxCore/SimplnxCore_test_dirs.hpp" +#include "ComputeSurfaceAreaToVolumeTestUtils.hpp" + +#include "simplnx/Common/Aliases.hpp" +#include "simplnx/Common/Array.hpp" +#include "simplnx/Common/Result.hpp" +#include "simplnx/Common/Types.hpp" #include "simplnx/DataStructure/AttributeMatrix.hpp" #include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/DataPath.hpp" +#include "simplnx/DataStructure/DataStore.hpp" +#include "simplnx/DataStructure/DataStructure.hpp" #include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/DataStructure/IDataArray.hpp" +#include "simplnx/DataStructure/NeighborList.hpp" +#include "simplnx/DataStructure/StringArray.hpp" +#include "simplnx/Filter/Arguments.hpp" +#include "simplnx/Filter/IFilter.hpp" #include "simplnx/Parameters/ArrayCreationParameter.hpp" #include "simplnx/Parameters/BoolParameter.hpp" +#include "simplnx/UnitTest/AlgorithmTestScope.hpp" #include "simplnx/UnitTest/UnitTestCommon.hpp" #include "simplnx/Utilities/AlgorithmDispatch.hpp" #include "simplnx/Utilities/DataStoreUtilities.hpp" #include - +#include +#include + +#include +#include +#include +#include +#include #include +#include +#include +#include +#include namespace fs = std::filesystem; @@ -32,6 +60,71 @@ constexpr usize k_TotalVoxels = k_Dim * k_Dim * k_Dim; constexpr int32 k_NumOctantFeatures = 8; const std::string k_GeomName = "ImageGeom"; +// Resident test stores observe Scanline transfers without adding a production fault hook. +template +class SurfaceMetricTransferStore : public DataStore +{ +public: + SurfaceMetricTransferStore(usize numTuples, T initialValue) + : DataStore(numTuples, initialValue) + { + } + + /** + * @brief Selects the read that returns an injected error before copying data. + * @param readNumber One-based call number; zero disables the injected error. + */ + void failReadOn(usize readNumber) + { + m_FailRead = readNumber; + } + + /** + * @brief Selects the write that returns an injected error before copying data. + * @param writeNumber One-based call number; zero disables the injected error. + */ + void failWriteOn(usize writeNumber) + { + m_FailWrite = writeNumber; + } + + const std::vector>& getReads() const + { + return m_Reads; + } + + const std::vector>& getWrites() const + { + return m_Writes; + } + + Result<> copyIntoBuffer(usize offset, nonstd::span buffer) const override + { + m_Reads.push_back({offset, buffer.size()}); + if(m_Reads.size() == m_FailRead) + { + return MakeErrorResult(-94801, "Injected surface metric input read failure"); + } + return DataStore::copyIntoBuffer(offset, buffer); + } + + Result<> copyFromBuffer(usize offset, nonstd::span buffer) override + { + m_Writes.push_back({offset, buffer.size()}); + if(m_Writes.size() == m_FailWrite) + { + return MakeErrorResult(-94802, "Injected surface metric output write failure"); + } + return DataStore::copyFromBuffer(offset, buffer); + } + +private: + mutable std::vector> m_Reads; + std::vector> m_Writes; + usize m_FailRead = 0; + usize m_FailWrite = 0; +}; + void BuildOctantWithNumCells(DataStructure& ds) { const ShapeType cellTupleShape = {k_Dim, k_Dim, k_Dim}; @@ -86,6 +179,380 @@ void BuildOctantWithNumCells(DataStructure& ds) } } // namespace +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeFilter: sparse feature table batch boundaries", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + const auto scenario = GENERATE(from_range(UnitTest::SelectAlgorithmTestScenariosForInMemoryStores())); + const auto numFeatures = GENERATE(usize{1}, usize{2}, usize{3}, usize{7}, usize{131071}, usize{131072}, usize{131073}); + const bool calculateSphericity = GENERATE(false, true); + CAPTURE(scenario, numFeatures, calculateSphericity); + UnitTest::AlgorithmTestScope scope(scenario); + surface_area_to_volume_test::CheckSparseFeatureSurfaceMetrics(scope, numFeatures, calculateSphericity); +} + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeScanline: bounded metric transfers", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + // This targets Scanline directly: the test contract is its bulk-transfer bound, not dispatch. + const auto featureAndBatchCounts = GENERATE(std::pair{0, 1}, std::pair{1, 1}, std::pair{2, 1}, std::pair{3, 1}, std::pair{7, 3}, + std::pair{131071, 65535}, std::pair{131072, 65536}, std::pair{131073, 65536}); + const auto [numFeatures, batchCapacity] = featureAndBatchCounts; + const bool calculateSphericity = GENERATE(false, true); + CAPTURE(numFeatures, batchCapacity, calculateSphericity); + DataStructure dataStructure; + const auto occupiedFeatures = surface_area_to_volume_test::BuildSparseFeatureSurfaceData(dataStructure, numFeatures); + const DataPath featurePath({"ImageGeom", "CellFeatureData"}); + const DataPath countsPath = featurePath.createChildPath("NumCells"); + REQUIRE_NOTHROW(dataStructure.getDataRefAs(featurePath)); + const auto featureId = dataStructure.getDataRefAs(featurePath).getId(); + auto countsStorePtr = std::make_shared>(numFeatures, 0); + for(const usize feature : occupiedFeatures) + { + countsStorePtr->setValue(feature, 1); + } + REQUIRE_NOTHROW(dataStructure.getDataRefAs(countsPath)); + const auto replacement = dataStructure.getDataRefAs(countsPath).setDataStore(countsStorePtr); + SIMPLNX_RESULT_REQUIRE_VALID(replacement); + auto ratioStorePtr = std::make_shared>(numFeatures, -47.0f); + REQUIRE(Float32Array::Create(dataStructure, "SurfaceAreaVolumeRatio", ratioStorePtr, featureId) != nullptr); + auto sphericityStorePtr = std::make_shared>(numFeatures, -47.0f); + if(calculateSphericity) + { + REQUIRE(Float32Array::Create(dataStructure, "Sphericity", sphericityStorePtr, featureId) != nullptr); + } + const ComputeSurfaceAreaToVolumeInputValues inputValues = {.FeatureIdsArrayPath = DataPath({"ImageGeom", "CellData", "FeatureIds"}), + .NumCellsArrayPath = countsPath, + .SurfaceAreaVolumeRatioArrayName = featurePath.createChildPath("SurfaceAreaVolumeRatio"), + .CalculateSphericity = calculateSphericity, + .SphericityArrayName = featurePath.createChildPath("Sphericity"), + .InputImageGeometry = DataPath({"ImageGeom"})}; + const std::atomic_bool shouldCancel = false; + const IFilter::MessageHandler messageHandler; + const auto result = ComputeSurfaceAreaToVolumeScanline(dataStructure, messageHandler, shouldCancel, &inputValues)(); + SIMPLNX_RESULT_REQUIRE_VALID(result); + + std::vector> expectedTransfers; + for(usize offset = 0; offset < numFeatures; offset += batchCapacity) + { + expectedTransfers.push_back({offset, std::min(batchCapacity, numFeatures - offset)}); + } + CHECK(countsStorePtr->getReads() == expectedTransfers); + CHECK(ratioStorePtr->getWrites() == expectedTransfers); + if(calculateSphericity) + { + CHECK(sphericityStorePtr->getWrites() == expectedTransfers); + } + else + { + CHECK(sphericityStorePtr->getWrites().empty()); + } + if(numFeatures > 0) + { + CHECK(ratioStorePtr->getValue(0) == 0.0f); + if(calculateSphericity) + { + CHECK(sphericityStorePtr->getValue(0) == 0.0f); + } + } + UnitTest::CheckArraysInheritTupleDims(dataStructure); +} + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeScanline: later metric I/O failures preserve prior batches", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + // Exercise each failed operation on the second batch, after the first batch was written. + const int failure = GENERATE(0, 1, 2); // NumCells read, ratio write, sphericity write. + CAPTURE(failure); + constexpr usize numFeatures = 7; + constexpr float32 unwritten = -47.0f; + DataStructure dataStructure; + const auto occupiedFeatures = surface_area_to_volume_test::BuildSparseFeatureSurfaceData(dataStructure, numFeatures); + const DataPath featurePath({"ImageGeom", "CellFeatureData"}); + const DataPath countsPath = featurePath.createChildPath("NumCells"); + REQUIRE_NOTHROW(dataStructure.getDataRefAs(featurePath)); + const auto featureId = dataStructure.getDataRefAs(featurePath).getId(); + auto countsStorePtr = std::make_shared>(numFeatures, 0); + countsStorePtr->failReadOn(failure == 0 ? 2 : 0); + for(const usize feature : occupiedFeatures) + { + countsStorePtr->setValue(feature, 1); + } + REQUIRE_NOTHROW(dataStructure.getDataRefAs(countsPath)); + const auto replacement = dataStructure.getDataRefAs(countsPath).setDataStore(countsStorePtr); + SIMPLNX_RESULT_REQUIRE_VALID(replacement); + auto ratioStorePtr = std::make_shared>(numFeatures, unwritten); + ratioStorePtr->failWriteOn(failure == 1 ? 2 : 0); + auto sphericityStorePtr = std::make_shared>(numFeatures, unwritten); + sphericityStorePtr->failWriteOn(failure == 2 ? 2 : 0); + REQUIRE(Float32Array::Create(dataStructure, "SurfaceAreaVolumeRatio", ratioStorePtr, featureId) != nullptr); + REQUIRE(Float32Array::Create(dataStructure, "Sphericity", sphericityStorePtr, featureId) != nullptr); + const ComputeSurfaceAreaToVolumeInputValues inputValues = {.FeatureIdsArrayPath = DataPath({"ImageGeom", "CellData", "FeatureIds"}), + .NumCellsArrayPath = countsPath, + .SurfaceAreaVolumeRatioArrayName = featurePath.createChildPath("SurfaceAreaVolumeRatio"), + .CalculateSphericity = true, + .SphericityArrayName = featurePath.createChildPath("Sphericity"), + .InputImageGeometry = DataPath({"ImageGeom"})}; + const std::atomic_bool shouldCancel = false; + const IFilter::MessageHandler messageHandler; + const auto result = ComputeSurfaceAreaToVolumeScanline(dataStructure, messageHandler, shouldCancel, &inputValues)(); + SIMPLNX_RESULT_REQUIRE_INVALID(result); + REQUIRE(result.errors().size() == 1); + CHECK(result.errors().front().code == (failure == 0 ? -94801 : -94802)); + CHECK(countsStorePtr->getReads().size() == 2); + CHECK(ratioStorePtr->getWrites().size() == (failure == 0 ? 1 : 2)); + CHECK(sphericityStorePtr->getWrites().size() == (failure == 2 ? 2 : 1)); + + const double piValue = std::acos(-1.0); + const double sphereRadius = std::cbrt(3.0 * 30.0 / (4.0 * piValue)); + const double expectedSphericity = 4.0 * piValue * sphereRadius * sphereRadius / 62.0; + CHECK(ratioStorePtr->getValue(0) == 0.0f); + CHECK(sphericityStorePtr->getValue(0) == 0.0f); + for(usize feature = 1; feature < numFeatures; ++feature) + { + CAPTURE(feature); + const bool occupied = std::ranges::binary_search(occupiedFeatures, feature); + if(feature < (failure == 2 ? 6 : 3)) + { + if(occupied) + { + CHECK(ratioStorePtr->getValue(feature) == Approx(62.0 / 30.0).epsilon(1.0e-6)); + } + else + { + CHECK(std::isnan(ratioStorePtr->getValue(feature))); + } + } + else + { + CHECK(ratioStorePtr->getValue(feature) == unwritten); + } + if(feature < 3) + { + CHECK(sphericityStorePtr->getValue(feature) == Approx(expectedSphericity).epsilon(1.0e-6)); + } + else + { + CHECK(sphericityStorePtr->getValue(feature) == unwritten); + } + } + UnitTest::CheckArraysInheritTupleDims(dataStructure); +} + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeFilter: geometry and FeatureIds tuple cardinality", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + const auto labelDimensions = GENERATE(SizeVec3{2, 2, 2}, SizeVec3{6, 5, 4}, SizeVec3{6, 5, 2}); + const auto calculateSphericity = GENERATE(false, true); + const ShapeType labelShape = {labelDimensions[2], labelDimensions[1], labelDimensions[0]}; + const usize labelCount = labelDimensions[0] * labelDimensions[1] * labelDimensions[2]; + CAPTURE(labelCount, calculateSphericity); + + DataStructure dataStructure; + const DataPath geomPath({"SelectedGeometry"}); + const DataPath labelGeomPath({"LabelGeometry"}); + const DataPath idsPath = labelGeomPath.createChildPath("CellData").createChildPath("FeatureIds"); + const DataPath countsPath = geomPath.createChildPath("CellFeatureData").createChildPath("NumCells"); + const DataPath ratioPath = countsPath.replaceName("SurfaceAreaVolumeRatio"); + const DataPath sphericityPath = countsPath.replaceName("Sphericity"); + + auto* imageGeomPtr = ImageGeom::Create(dataStructure, geomPath.getTargetName()); + REQUIRE(imageGeomPtr != nullptr); + imageGeomPtr->setDimensions({4, 3, 5}); + auto* cellAttributeMatrixPtr = AttributeMatrix::Create(dataStructure, "CellData", {5, 3, 4}, imageGeomPtr->getId()); + REQUIRE(cellAttributeMatrixPtr != nullptr); + imageGeomPtr->setCellData(*cellAttributeMatrixPtr); + auto* featureAttributeMatrixPtr = AttributeMatrix::Create(dataStructure, "CellFeatureData", {2}, imageGeomPtr->getId()); + REQUIRE(featureAttributeMatrixPtr != nullptr); + + auto* labelGeomPtr = ImageGeom::Create(dataStructure, labelGeomPath.getTargetName()); + REQUIRE(labelGeomPtr != nullptr); + labelGeomPtr->setDimensions(labelDimensions); + auto* labelAttributeMatrixPtr = AttributeMatrix::Create(dataStructure, "CellData", labelShape, labelGeomPtr->getId()); + REQUIRE(labelAttributeMatrixPtr != nullptr); + labelGeomPtr->setCellData(*labelAttributeMatrixPtr); + + auto idsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, idsPath, labelShape, {1}); + auto countsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, countsPath, {2}, {1}); + REQUIRE(idsStorePtr != nullptr); + REQUIRE(countsStorePtr != nullptr); + REQUIRE(Int32Array::Create(dataStructure, "FeatureIds", idsStorePtr, labelAttributeMatrixPtr->getId()) != nullptr); + REQUIRE(Int32Array::Create(dataStructure, "NumCells", countsStorePtr, featureAttributeMatrixPtr->getId()) != nullptr); + + const ComputeSurfaceAreaToVolumeFilter filter; + Arguments args = filter.getDefaultArguments(); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SelectedImageGeometryPath_Key, std::make_any(geomPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CellFeatureIdsArrayPath_Key, std::make_any(idsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_NumCellsArrayPath_Key, std::make_any(countsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CalculateSphericity_Key, std::make_any(calculateSphericity)); + + auto preflight = filter.preflight(dataStructure, args); + if(labelCount == 60) + { + SIMPLNX_RESULT_REQUIRE_VALID(preflight.outputActions); + CHECK(preflight.outputActions.value().actions.size() == (calculateSphericity ? 2 : 1)); + } + else + { + // An invalid output-actions result exposes no actions that could allocate outputs for the unsafe input. + SIMPLNX_RESULT_REQUIRE_INVALID(preflight.outputActions); + REQUIRE(preflight.outputActions.errors().size() == 1); + const auto& error = preflight.outputActions.errors().front(); + CHECK(error.code == -12804); + CHECK(error.message.find(idsPath.toString()) != std::string::npos); + CHECK(error.message.find(geomPath.toString()) != std::string::npos); + CHECK(error.message.find(fmt::format("{} tuples", labelCount)) != std::string::npos); + CHECK(error.message.find("60 cells") != std::string::npos); + } + CHECK(dataStructure.getDataAs(ratioPath) == nullptr); + CHECK(dataStructure.getDataAs(sphericityPath) == nullptr); + UnitTest::CheckArraysInheritTupleDims(dataStructure); +} + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeFilter: rejects non-DataArray selections", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + const bool replaceFeatureIds = GENERATE(true, false); + const bool useStringArray = GENERATE(true, false); + CAPTURE(replaceFeatureIds, useStringArray); + + DataStructure dataStructure; + const DataPath geomPath({"ImageGeom"}); + const DataPath cellPath = geomPath.createChildPath("CellData"); + const DataPath featurePath = geomPath.createChildPath("CellFeatureData"); + const DataPath idsPath = cellPath.createChildPath("FeatureIds"); + const DataPath countsPath = featurePath.createChildPath("NumCells"); + const DataPath ratioPath = featurePath.createChildPath("SurfaceAreaVolumeRatio"); + const DataPath sphericityPath = featurePath.createChildPath("Sphericity"); + const ShapeType cellShape = {2, 2, 2}; + const ShapeType featureShape = {2}; + const std::vector labels(8, 1); + const std::vector counts = {0, 8}; + + auto* imageGeomPtr = ImageGeom::Create(dataStructure, geomPath.getTargetName()); + REQUIRE(imageGeomPtr != nullptr); + imageGeomPtr->setDimensions({2, 2, 2}); + auto* cellAttributeMatrixPtr = AttributeMatrix::Create(dataStructure, cellPath.getTargetName(), cellShape, imageGeomPtr->getId()); + REQUIRE(cellAttributeMatrixPtr != nullptr); + imageGeomPtr->setCellData(*cellAttributeMatrixPtr); + auto* featureAttributeMatrixPtr = AttributeMatrix::Create(dataStructure, featurePath.getTargetName(), featureShape, imageGeomPtr->getId()); + REQUIRE(featureAttributeMatrixPtr != nullptr); + auto idsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, idsPath, cellShape, {1}); + auto countsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, countsPath, featureShape, {1}); + REQUIRE(idsStorePtr != nullptr); + REQUIRE(countsStorePtr != nullptr); + REQUIRE(Int32Array::Create(dataStructure, idsPath.getTargetName(), idsStorePtr, cellAttributeMatrixPtr->getId()) != nullptr); + REQUIRE(Int32Array::Create(dataStructure, countsPath.getTargetName(), countsStorePtr, featureAttributeMatrixPtr->getId()) != nullptr); + const auto idsWrite = idsStorePtr->copyFromBuffer(0, nonstd::span(labels.data(), labels.size())); + SIMPLNX_RESULT_REQUIRE_VALID(idsWrite); + const auto countsWrite = countsStorePtr->copyFromBuffer(0, nonstd::span(counts.data(), counts.size())); + SIMPLNX_RESULT_REQUIRE_VALID(countsWrite); + + const auto* parentPtr = replaceFeatureIds ? cellAttributeMatrixPtr : featureAttributeMatrixPtr; + const DataPath parentPath = replaceFeatureIds ? cellPath : featurePath; + const auto& values = replaceFeatureIds ? labels : counts; + const std::string malformedName = useStringArray ? "StringSelection" : "NeighborListSelection"; + const DataPath malformedPath = parentPath.createChildPath(malformedName); + REQUIRE(values.size() == parentPtr->getNumberOfTuples()); + if(useStringArray) + { + std::vector strings; + strings.reserve(values.size()); + for(const int32 value : values) + { + strings.push_back(std::to_string(value)); + } + auto* stringsPtr = StringArray::CreateWithValues(dataStructure, malformedName, parentPtr->getShape(), std::move(strings), parentPtr->getId()); + REQUIRE(stringsPtr != nullptr); + } + else + { + auto* listsPtr = NeighborList::Create(dataStructure, malformedName, parentPtr->getShape(), parentPtr->getId()); + REQUIRE(listsPtr != nullptr); + for(usize index = 0; index < values.size(); ++index) + { + listsPtr->setList(static_cast(index), std::vector{values[index]}); + } + CHECK(listsPtr->getComponentShape() == ShapeType{1}); + } + + const ComputeSurfaceAreaToVolumeFilter filter; + Arguments args = filter.getDefaultArguments(); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SelectedImageGeometryPath_Key, std::make_any(geomPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CellFeatureIdsArrayPath_Key, std::make_any(replaceFeatureIds ? malformedPath : idsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_NumCellsArrayPath_Key, std::make_any(replaceFeatureIds ? countsPath : malformedPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CalculateSphericity_Key, std::make_any(true)); + + IFilter::PreflightResult preflight; + REQUIRE_NOTHROW(preflight = filter.preflight(dataStructure, args)); + SIMPLNX_RESULT_REQUIRE_INVALID(preflight.outputActions); + REQUIRE(preflight.outputActions.errors().size() == 1); + const auto& error = preflight.outputActions.errors().front(); + CHECK(error.code == (replaceFeatureIds ? -12803 : -12801)); + CHECK(error.message.find(malformedPath.toString()) != std::string::npos); + if(replaceFeatureIds) + { + CHECK(error.message.find(useStringArray ? "StringArray" : "NeighborList") != std::string::npos); + CHECK(error.message.find("scalar int32 DataArray") != std::string::npos); + CHECK(error.message.find("Select a scalar int32 DataArray") != std::string::npos); + } + CHECK(dataStructure.getData(ratioPath) == nullptr); + CHECK(dataStructure.getData(sphericityPath) == nullptr); + UnitTest::CheckArraysInheritTupleDims(dataStructure); +} + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeFilter: analytical interior pair", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + const auto scenario = GENERATE(from_range(UnitTest::SelectAlgorithmTestScenariosForInMemoryStores())); + const auto pairAxis = GENERATE(usize{0}, usize{1}, usize{2}); + const auto spacing = GENERATE(FloatVec3{2.0f, 3.0f, 5.0f}, FloatVec3{1.0f, 1.0f, 1.0f}); + const auto calculateSphericity = GENERATE(false, true); + CAPTURE(scenario); + UnitTest::AlgorithmTestScope scope(scenario); + surface_area_to_volume_test::CheckInteriorPairSurfaceMetrics(scope, pairAxis, spacing, calculateSphericity); +} + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeFilter: analytical single-slice pair", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + const auto scenario = GENERATE(from_range(UnitTest::SelectAlgorithmTestScenariosForInMemoryStores())); + const auto pairAxis = GENERATE(usize{0}, usize{1}); + const auto calculateSphericity = GENERATE(false, true); + CAPTURE(scenario); + UnitTest::AlgorithmTestScope scope(scenario); + surface_area_to_volume_test::CheckInteriorPairSurfaceMetrics(scope, pairAxis, FloatVec3{2.0f, 3.0f, 5.0f}, calculateSphericity, true); +} + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeFilter: decimal connected comb", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + const auto scenario = GENERATE(from_range(UnitTest::SelectAlgorithmTestScenariosForInMemoryStores())); + const auto spacing = GENERATE(FloatVec3{0.13f, 0.27f, 1.1f}, FloatVec3{0.13f, 0.13f, 0.13f}); + const auto calculateSphericity = GENERATE(false, true); + CAPTURE(scenario); + UnitTest::AlgorithmTestScope scope(scenario); + surface_area_to_volume_test::CheckConnectedCombSurfaceMetrics(scope, spacing, calculateSphericity); +} + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeFilter: scale-invariant sphericity", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + const auto scenario = GENERATE(from_range(UnitTest::SelectAlgorithmTestScenariosForInMemoryStores())); + CAPTURE(scenario); + UnitTest::AlgorithmTestScope scope(scenario); + surface_area_to_volume_test::CheckScaleInvariantSphericity(scope); +} + +TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolumeFilter: nonpositive supplied pair counts", "[SimplnxCore][ComputeSurfaceAreaToVolumeFilter]") +{ + UnitTest::LoadPlugins(); + const auto scenario = GENERATE(from_range(UnitTest::SelectAlgorithmTestScenariosForInMemoryStores())); + const auto suppliedPairCount = GENERATE(int32{0}, int32{-2}); + CAPTURE(scenario, suppliedPairCount); + UnitTest::AlgorithmTestScope scope(scenario); + surface_area_to_volume_test::CheckInteriorPairSurfaceMetrics(scope, 0, FloatVec3{2.0f, 3.0f, 5.0f}, true, false, suppliedPairCount); +} + TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolume", "[SimplnxCore][ComputeSurfaceAreaToVolume]") { UnitTest::LoadPlugins(); @@ -135,6 +602,7 @@ TEST_CASE("SimplnxCore::ComputeSurfaceAreaToVolume", "[SimplnxCore][ComputeSurfa calculatedPath = DataPath({k_DataContainer, k_CellFeatureData, k_SphericityArrayNameNX}); REQUIRE_NOTHROW(dataStructure.getDataRefAs(exemplarPath)); REQUIRE_NOTHROW(dataStructure.getDataRefAs(calculatedPath)); + // The archived exemplar predates the sphericity exponent correction. Its absolute tolerance can hide small output differences; analytical cases establish the corrected exponent behavior. CompareDataArrays(dataStructure.getDataRefAs(exemplarPath), dataStructure.getDataRefAs(calculatedPath)); } diff --git a/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTestUtils.hpp b/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTestUtils.hpp new file mode 100644 index 0000000000..f21f90ccf2 --- /dev/null +++ b/src/Plugins/SimplnxCore/test/ComputeSurfaceAreaToVolumeTestUtils.hpp @@ -0,0 +1,533 @@ +#pragma once + +#include "SimplnxCore/Filters/ComputeSurfaceAreaToVolumeFilter.hpp" + +#include "simplnx/Common/Aliases.hpp" +#include "simplnx/Common/Array.hpp" +#include "simplnx/Common/Result.hpp" +#include "simplnx/Common/Types.hpp" +#include "simplnx/DataStructure/AttributeMatrix.hpp" +#include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/DataPath.hpp" +#include "simplnx/DataStructure/DataStructure.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/DataStructure/IDataArray.hpp" +#include "simplnx/Filter/Arguments.hpp" +#include "simplnx/UnitTest/UnitTestCommon.hpp" +#include "simplnx/Utilities/DataStoreUtilities.hpp" + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace surface_area_to_volume_test +{ +/** + * @brief Builds isolated cells that exercise large feature tables in a small image. + * @param dataStructure Receives the geometry, labels, and supplied cell counts. + * @param numFeatures Number of feature tuples, including background. Zero creates empty inputs. + * @return Sorted feature identifiers for the isolated one-cell boxes. + */ +inline std::vector BuildSparseFeatureSurfaceData(nx::core::DataStructure& dataStructure, nx::core::usize numFeatures) +{ + using namespace nx::core; + const DataPath geomPath({"ImageGeom"}); + const DataPath cellPath = geomPath.createChildPath("CellData"); + const DataPath featurePath = geomPath.createChildPath("CellFeatureData"); + const DataPath idsPath = cellPath.createChildPath("FeatureIds"); + const DataPath countsPath = featurePath.createChildPath("NumCells"); + const ShapeType cellShape = numFeatures == 0 ? ShapeType{0, 0, 0} : ShapeType{3, 3, 19}; + const ShapeType featureShape = {numFeatures}; + auto* geomPtr = ImageGeom::Create(dataStructure, geomPath.getTargetName()); + REQUIRE(geomPtr != nullptr); + geomPtr->setDimensions({cellShape[2], cellShape[1], cellShape[0]}); + geomPtr->setSpacing({2.0f, 3.0f, 5.0f}); + auto* cellAMPtr = AttributeMatrix::Create(dataStructure, cellPath.getTargetName(), cellShape, geomPtr->getId()); + REQUIRE(cellAMPtr != nullptr); + geomPtr->setCellData(*cellAMPtr); + auto* featureAMPtr = AttributeMatrix::Create(dataStructure, featurePath.getTargetName(), featureShape, geomPtr->getId()); + REQUIRE(featureAMPtr != nullptr); + + auto idsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, idsPath, cellShape, {1}); + auto countsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, countsPath, featureShape, {1}); + REQUIRE(idsStorePtr != nullptr); + REQUIRE(countsStorePtr != nullptr); + REQUIRE(Int32Array::Create(dataStructure, idsPath.getTargetName(), idsStorePtr, cellAMPtr->getId()) != nullptr); + REQUIRE(Int32Array::Create(dataStructure, countsPath.getTargetName(), countsStorePtr, featureAMPtr->getId()) != nullptr); + if(numFeatures == 0) + { + return {}; + } + + // Probe both sides of the half-table and fixed-cap boundaries, plus the last tuple. + std::vector occupiedFeatures = {1, numFeatures / 2, (numFeatures / 2) + 1, numFeatures - 1, 65535, 65536, 65537}; + if(numFeatures > 3) + { + occupiedFeatures.push_back((numFeatures / 2) - 1); + } + std::erase_if(occupiedFeatures, [numFeatures](usize feature) { return feature == 0 || feature >= numFeatures; }); + std::ranges::sort(occupiedFeatures); + occupiedFeatures.erase(std::ranges::unique(occupiedFeatures).begin(), occupiedFeatures.end()); + REQUIRE(occupiedFeatures.size() <= 8); + std::vector labels(usize{19} * 3 * 3, 0); + std::vector counts(numFeatures, 0); + for(usize index = 0; index < occupiedFeatures.size(); ++index) + { + // Background separates every one-cell box from its neighbors and the image exterior. + labels[(((1 * 3) + 1) * 19) + 1 + (2 * index)] = static_cast(occupiedFeatures[index]); + counts[occupiedFeatures[index]] = 1; + } + counts[0] = static_cast(labels.size() - occupiedFeatures.size()); + const auto idsWrite = idsStorePtr->copyFromBuffer(0, nonstd::span(labels.data(), labels.size())); + SIMPLNX_RESULT_REQUIRE_VALID(idsWrite); + const auto countsWrite = countsStorePtr->copyFromBuffer(0, nonstd::span(counts.data(), counts.size())); + SIMPLNX_RESULT_REQUIRE_VALID(countsWrite); + return occupiedFeatures; +} + +/** + * @brief Checks all sparse-table metrics against box and equal-volume sphere geometry. + * @tparam ScopeT Verifies the selected algorithm and actual array storage. + * @param scope Active algorithm and storage scenario. + * @param numFeatures Number of feature tuples, including background. Must be positive. + * @param calculateSphericity Enables the optional sphere-area ratio check. + */ +template +void CheckSparseFeatureSurfaceMetrics(ScopeT& scope, nx::core::usize numFeatures, bool calculateSphericity) +{ + using namespace nx::core; + REQUIRE(numFeatures > 0); + DataStructure dataStructure; + const auto occupiedFeatures = BuildSparseFeatureSurfaceData(dataStructure, numFeatures); + const DataPath geomPath({"ImageGeom"}); + const DataPath idsPath({"ImageGeom", "CellData", "FeatureIds"}); + const DataPath featurePath({"ImageGeom", "CellFeatureData"}); + const DataPath countsPath = featurePath.createChildPath("NumCells"); + const DataPath ratioPath = featurePath.createChildPath("SurfaceAreaVolumeRatio"); + const DataPath sphericityPath = featurePath.createChildPath("Sphericity"); + REQUIRE_NOTHROW(dataStructure.getDataRefAs(idsPath)); + REQUIRE_NOTHROW(dataStructure.getDataRefAs(countsPath)); + scope.requireExpectedStore(dataStructure.getDataRefAs(idsPath)); + scope.requireExpectedStore(dataStructure.getDataRefAs(countsPath)); + + const ComputeSurfaceAreaToVolumeFilter filter; + Arguments args = filter.getDefaultArguments(); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SelectedImageGeometryPath_Key, std::make_any(geomPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CellFeatureIdsArrayPath_Key, std::make_any(idsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_NumCellsArrayPath_Key, std::make_any(countsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CalculateSphericity_Key, std::make_any(calculateSphericity)); + const auto result = scope.executeFilter(filter, dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(result.result); + + REQUIRE_NOTHROW(dataStructure.getDataRefAs(ratioPath)); + const auto& ratios = dataStructure.getDataRefAs(ratioPath); + scope.requireExpectedStore(ratios); + CHECK(ratios.getTupleShape() == ShapeType{numFeatures}); + CHECK(ratios.getComponentShape() == ShapeType{1}); + std::vector values(numFeatures); + const auto ratioRead = ratios.getDataStoreRef().copyIntoBuffer(0, nonstd::span(values.data(), values.size())); + SIMPLNX_RESULT_REQUIRE_VALID(ratioRead); + CHECK(values[0] == 0.0f); + // An isolated 2-by-3-by-5 box has area 62 and volume 30. + for(usize feature = 1; feature < numFeatures; ++feature) + { + CAPTURE(feature); + if(std::ranges::binary_search(occupiedFeatures, feature)) + { + REQUIRE(values[feature] == Approx(62.0 / 30.0).epsilon(1.0e-6)); + } + else + { + REQUIRE(std::isnan(values[feature])); + } + } + if(calculateSphericity) + { + REQUIRE_NOTHROW(dataStructure.getDataRefAs(sphericityPath)); + const auto& sphericity = dataStructure.getDataRefAs(sphericityPath); + scope.requireExpectedStore(sphericity); + CHECK(sphericity.getTupleShape() == ShapeType{numFeatures}); + CHECK(sphericity.getComponentShape() == ShapeType{1}); + const auto sphericityRead = sphericity.getDataStoreRef().copyIntoBuffer(0, nonstd::span(values.data(), values.size())); + SIMPLNX_RESULT_REQUIRE_VALID(sphericityRead); + CHECK(values[0] == 0.0f); + // Independently calculate the equal-volume sphere from its radius. + const double piValue = std::acos(-1.0); + const double radius = std::cbrt(3.0 * 30.0 / (4.0 * piValue)); + const double expected = 4.0 * piValue * radius * radius / 62.0; + for(usize feature = 1; feature < numFeatures; ++feature) + { + CAPTURE(feature); + if(std::ranges::binary_search(occupiedFeatures, feature)) + { + REQUIRE(values[feature] == Approx(expected).epsilon(1.0e-6)); + } + else + { + REQUIRE(std::isnan(values[feature])); + } + } + } + else + { + CHECK(dataStructure.getData(sphericityPath) == nullptr); + } + UnitTest::CheckArraysInheritTupleDims(dataStructure); +} + +// A wrong X/Y face weight changes the physical area of this interior box. +template +std::array CheckInteriorPairSurfaceMetrics(ScopeT& scope, nx::core::usize pairAxis, const nx::core::FloatVec3& spacing, bool calculateSphericity, bool singleSlice = false, + nx::core::int32 suppliedPairCount = 2) +{ + using namespace nx::core; + const SizeVec3 dimensions = singleSlice ? SizeVec3{5, 5, 1} : SizeVec3{5, 5, 5}; + DataStructure dataStructure; + const DataPath geomPath({"ImageGeom"}); + const DataPath cellPath = geomPath.createChildPath("CellData"); + const DataPath featurePath = geomPath.createChildPath("CellFeatureData"); + const DataPath idsPath = cellPath.createChildPath("FeatureIds"); + const DataPath countsPath = featurePath.createChildPath("NumCells"); + const DataPath ratioPath = featurePath.createChildPath("SurfaceAreaVolumeRatio"); + const DataPath sphericityPath = featurePath.createChildPath("Sphericity"); + const ShapeType cellShape = {dimensions[2], dimensions[1], dimensions[0]}; + const ShapeType featureShape = {3}; + const ShapeType componentShape = {1}; + const FloatVec3 origin = {-1.0f, 4.0f, 7.0f}; + REQUIRE(pairAxis < (singleSlice ? 2 : 3)); + constexpr usize pairFeature = 2; + constexpr int32 enclosingFeature = 1; + const usize cellCount = dimensions[0] * dimensions[1] * dimensions[2]; + + auto* imageGeomPtr = ImageGeom::Create(dataStructure, "ImageGeom"); + REQUIRE(imageGeomPtr != nullptr); + imageGeomPtr->setDimensions(dimensions); + imageGeomPtr->setSpacing(spacing); + imageGeomPtr->setOrigin(origin); + auto* cellAttributeMatrixPtr = AttributeMatrix::Create(dataStructure, "CellData", cellShape, imageGeomPtr->getId()); + REQUIRE(cellAttributeMatrixPtr != nullptr); + imageGeomPtr->setCellData(*cellAttributeMatrixPtr); + auto* featureAttributeMatrixPtr = AttributeMatrix::Create(dataStructure, "CellFeatureData", featureShape, imageGeomPtr->getId()); + REQUIRE(featureAttributeMatrixPtr != nullptr); + + auto idsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, idsPath, cellShape, componentShape); + auto countsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, countsPath, featureShape, componentShape); + REQUIRE(idsStorePtr != nullptr); + REQUIRE(countsStorePtr != nullptr); + REQUIRE(Int32Array::Create(dataStructure, "FeatureIds", idsStorePtr, cellAttributeMatrixPtr->getId()) != nullptr); + REQUIRE(Int32Array::Create(dataStructure, "NumCells", countsStorePtr, featureAttributeMatrixPtr->getId()) != nullptr); + std::vector labels(cellCount, enclosingFeature); + const std::array first = {2, 2, singleSlice ? usize{0} : usize{2}}; + auto second = first; + second.at(pairAxis)++; + labels[(((first[2] * dimensions[1]) + first[1]) * dimensions[0]) + first[0]] = static_cast(pairFeature); + labels[(((second[2] * dimensions[1]) + second[1]) * dimensions[0]) + second[0]] = static_cast(pairFeature); + const std::array counts = {0, singleSlice ? 23 : 123, suppliedPairCount}; + auto idsWrite = idsStorePtr->copyFromBuffer(0, nonstd::span(labels.data(), labels.size())); + SIMPLNX_RESULT_REQUIRE_VALID(idsWrite); + auto countsWrite = countsStorePtr->copyFromBuffer(0, nonstd::span(counts.data(), counts.size())); + SIMPLNX_RESULT_REQUIRE_VALID(countsWrite); + + REQUIRE_NOTHROW(dataStructure.getDataRefAs(idsPath)); + REQUIRE_NOTHROW(dataStructure.getDataRefAs(countsPath)); + const auto& ids = dataStructure.getDataRefAs(idsPath); + const auto& numCells = dataStructure.getDataRefAs(countsPath); + scope.requireExpectedStore(ids); + scope.requireExpectedStore(numCells); + const auto idsFormat = ids.getDataFormat(); + const auto countsFormat = numCells.getDataFormat(); + + const ComputeSurfaceAreaToVolumeFilter filter; + Arguments args; + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SelectedImageGeometryPath_Key, std::make_any(geomPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CellFeatureIdsArrayPath_Key, std::make_any(idsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_NumCellsArrayPath_Key, std::make_any(countsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CalculateSphericity_Key, std::make_any(calculateSphericity)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SurfaceAreaVolumeRatioArrayName_Key, std::make_any("SurfaceAreaVolumeRatio")); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SphericityArrayName_Key, std::make_any("Sphericity")); + auto preflight = filter.preflight(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(preflight.outputActions); + auto result = scope.executeFilter(filter, dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(result.result); + + // Independent rectangular-box oracle: the pair doubles one physical side. + std::array sides = {spacing[0], spacing[1], spacing[2]}; + sides.at(pairAxis) *= 2.0; + // Exterior Z faces are excluded for a single slice, as in the production boundary policy. + const double area = 2.0 * ((singleSlice ? 0.0 : sides[0] * sides[1]) + (sides[0] * sides[2]) + (sides[1] * sides[2])); + const double voxelVolume = static_cast(spacing[0]) * spacing[1] * spacing[2]; + std::array volumes = {0.0, 0.0, 0.0}; + volumes[pairFeature] = suppliedPairCount * voxelVolume; + volumes[enclosingFeature] = (singleSlice ? 23.0 : 123.0) * voxelVolume; + CAPTURE(pairAxis, spacing[0], spacing[1], spacing[2], calculateSphericity, singleSlice, suppliedPairCount, area); + REQUIRE(area > 0.0); + + REQUIRE_NOTHROW(dataStructure.getDataRefAs(ratioPath)); + const auto& ratios = dataStructure.getDataRefAs(ratioPath); + scope.requireExpectedStore(ratios); + CHECK(ratios.getTupleShape() == featureShape); + CHECK(ratios.getComponentShape() == componentShape); + std::array ratioValues = {}; + auto ratiosRead = ratios.getDataStoreRef().copyIntoBuffer(0, nonstd::span(ratioValues.data(), ratioValues.size())); + SIMPLNX_RESULT_REQUIRE_VALID(ratiosRead); + CHECK(ratioValues[0] == 0.0f); + for(usize feature = 1; feature < 3; feature++) + { + CAPTURE(feature); + if(volumes.at(feature) == 0.0) + { + CHECK(std::isinf(ratioValues.at(feature))); + CHECK(ratioValues.at(feature) > 0.0f); + } + else + { + CHECK(ratioValues.at(feature) == Approx(area / volumes.at(feature)).epsilon(1.0e-6)); + } + } + std::array sphericityValues = {}; + if(calculateSphericity) + { + REQUIRE_NOTHROW(dataStructure.getDataRefAs(sphericityPath)); + const auto& sphericity = dataStructure.getDataRefAs(sphericityPath); + scope.requireExpectedStore(sphericity); + CHECK(sphericity.getTupleShape() == featureShape); + CHECK(sphericity.getComponentShape() == componentShape); + auto sphericityRead = sphericity.getDataStoreRef().copyIntoBuffer(0, nonstd::span(sphericityValues.data(), sphericityValues.size())); + SIMPLNX_RESULT_REQUIRE_VALID(sphericityRead); + CHECK(sphericityValues[0] == 0.0f); + for(usize feature = 1; feature < 3; feature++) + { + CAPTURE(feature); + if(volumes.at(feature) < 0.0) + { + CHECK(std::isnan(sphericityValues.at(feature))); + } + else if(volumes.at(feature) == 0.0) + { + CHECK(sphericityValues.at(feature) == 0.0f); + } + else + { + // Derive the equal-volume sphere's area from its radius independently of the implementation formula. + const double piValue = std::acos(-1.0); + const double sphereRadius = std::cbrt(3.0 * volumes.at(feature) / (4.0 * piValue)); + const double expected = 4.0 * piValue * sphereRadius * sphereRadius / area; + CHECK(sphericityValues.at(feature) == Approx(expected).epsilon(1.0e-6)); + } + } + } + else + { + CHECK(dataStructure.getDataAs(sphericityPath) == nullptr); + } + + std::vector actualLabels(cellCount); + std::array actualCounts = {}; + auto idsRead = ids.getDataStoreRef().copyIntoBuffer(0, nonstd::span(actualLabels.data(), actualLabels.size())); + SIMPLNX_RESULT_REQUIRE_VALID(idsRead); + auto countsRead = numCells.getDataStoreRef().copyIntoBuffer(0, nonstd::span(actualCounts.data(), actualCounts.size())); + SIMPLNX_RESULT_REQUIRE_VALID(countsRead); + CHECK(actualLabels == labels); + CHECK(actualCounts == counts); + CHECK(ids.getTupleShape() == cellShape); + CHECK(ids.getComponentShape() == componentShape); + CHECK(numCells.getTupleShape() == featureShape); + CHECK(numCells.getComponentShape() == componentShape); + CHECK(ids.getDataFormat() == idsFormat); + CHECK(numCells.getDataFormat() == countsFormat); + CHECK(&ids.getDataStoreRef() == idsStorePtr.get()); + CHECK(&numCells.getDataStoreRef() == countsStorePtr.get()); + scope.requireExpectedStore(ids); + scope.requireExpectedStore(numCells); + REQUIRE_NOTHROW(dataStructure.getDataRefAs(geomPath)); + const auto& finalGeom = dataStructure.getDataRefAs(geomPath); + CHECK(finalGeom.getDimensions() == dimensions); + CHECK(finalGeom.getSpacing() == spacing); + CHECK(finalGeom.getOrigin() == origin); + CHECK(cellAttributeMatrixPtr->getShape() == cellShape); + CHECK(featureAttributeMatrixPtr->getShape() == featureShape); + UnitTest::CheckArraysInheritTupleDims(dataStructure); + return sphericityValues; +} + +template +void CheckScaleInvariantSphericity(ScopeT& scope) +{ + using namespace nx::core; + const auto reference = CheckInteriorPairSurfaceMetrics(scope, 0, FloatVec3{2.0f, 3.0f, 5.0f}, true); + for(const int exponent : {-20, 20}) + { + CAPTURE(exponent); + const FloatVec3 spacing = {std::ldexp(2.0f, exponent), std::ldexp(3.0f, exponent), std::ldexp(5.0f, exponent)}; + // The pair helper independently checks each scale against the equal-volume sphere oracle. + const auto scaled = CheckInteriorPairSurfaceMetrics(scope, 0, spacing, true); + for(usize feature = 1; feature < reference.size(); feature++) + { + CAPTURE(feature); + CHECK(scaled.at(feature) == Approx(reference.at(feature)).epsilon(1.0e-6)); + } + } +} +// A connected plate with long, separated teeth exposes repeated decimal-area accumulation. +template +void CheckConnectedCombSurfaceMetrics(ScopeT& scope, const nx::core::FloatVec3& spacing, bool calculateSphericity) +{ + using namespace nx::core; + DataStructure dataStructure; + const SizeVec3 dimensions = {137, 113, 37}; + const ShapeType cellShape = {dimensions[2], dimensions[1], dimensions[0]}; + const ShapeType featureShape = {2}; + const ShapeType componentShape = {1}; + const DataPath geomPath({"ImageGeom"}); + const DataPath cellPath = geomPath.createChildPath("CellData"); + const DataPath featurePath = geomPath.createChildPath("CellFeatureData"); + const DataPath idsPath = cellPath.createChildPath("FeatureIds"); + const DataPath countsPath = featurePath.createChildPath("NumCells"); + const DataPath ratioPath = featurePath.createChildPath("SurfaceAreaVolumeRatio"); + const DataPath sphericityPath = featurePath.createChildPath("Sphericity"); + const FloatVec3 origin = {-1.0f, 4.0f, 7.0f}; + const usize cellCount = dimensions[0] * dimensions[1] * dimensions[2]; + + // The x=1 backing plate is 111 by 35. Each odd (y,z) location extends a tooth + // another 134 cells along X. All teeth connect through the plate and the entire + // comb has a one-cell background border, so no exterior-domain faces enter the oracle. + constexpr int32 plateY = 111; + constexpr int32 plateZ = 35; + constexpr int32 toothLength = 134; + constexpr int32 teethY = 56; + constexpr int32 teethZ = 18; + constexpr int32 combCellCount = (plateY * plateZ) + (toothLength * teethY * teethZ); + constexpr std::array normalFaceCounts = {2 * plateY * plateZ, (2 * plateZ) + (2 * toothLength * teethY * teethZ), (2 * plateY) + (2 * toothLength * teethY * teethZ)}; + STATIC_REQUIRE(combCellCount == 138957); + STATIC_REQUIRE(normalFaceCounts[0] == 7770); + STATIC_REQUIRE(normalFaceCounts[1] == 270214); + STATIC_REQUIRE(normalFaceCounts[2] == 270366); + const std::array counts = {static_cast(cellCount) - combCellCount, combCellCount}; + std::vector labels(cellCount, 0); + for(usize zIndex = 1; zIndex <= plateZ; zIndex++) + { + for(usize yIndex = 1; yIndex <= plateY; yIndex++) + { + for(usize xIndex = 1; xIndex <= toothLength + 1; xIndex++) + { + if(xIndex == 1 || (yIndex % 2 == 1 && zIndex % 2 == 1)) + { + labels[(((zIndex * dimensions[1]) + yIndex) * dimensions[0]) + xIndex] = 1; + } + } + } + } + REQUIRE(std::count(labels.begin(), labels.end(), int32{1}) == combCellCount); + + auto* imageGeomPtr = ImageGeom::Create(dataStructure, "ImageGeom"); + REQUIRE(imageGeomPtr != nullptr); + imageGeomPtr->setDimensions(dimensions); + imageGeomPtr->setSpacing(spacing); + imageGeomPtr->setOrigin(origin); + auto* cellAttributeMatrixPtr = AttributeMatrix::Create(dataStructure, "CellData", cellShape, imageGeomPtr->getId()); + REQUIRE(cellAttributeMatrixPtr != nullptr); + imageGeomPtr->setCellData(*cellAttributeMatrixPtr); + auto* featureAttributeMatrixPtr = AttributeMatrix::Create(dataStructure, "CellFeatureData", featureShape, imageGeomPtr->getId()); + REQUIRE(featureAttributeMatrixPtr != nullptr); + auto idsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, idsPath, cellShape, componentShape); + auto countsStorePtr = DataStoreUtilities::CreateDataStore(dataStructure, countsPath, featureShape, componentShape); + REQUIRE(idsStorePtr != nullptr); + REQUIRE(countsStorePtr != nullptr); + REQUIRE(Int32Array::Create(dataStructure, "FeatureIds", idsStorePtr, cellAttributeMatrixPtr->getId()) != nullptr); + REQUIRE(Int32Array::Create(dataStructure, "NumCells", countsStorePtr, featureAttributeMatrixPtr->getId()) != nullptr); + auto idsWrite = idsStorePtr->copyFromBuffer(0, nonstd::span(labels.data(), labels.size())); + SIMPLNX_RESULT_REQUIRE_VALID(idsWrite); + auto countsWrite = countsStorePtr->copyFromBuffer(0, nonstd::span(counts.data(), counts.size())); + SIMPLNX_RESULT_REQUIRE_VALID(countsWrite); + + REQUIRE_NOTHROW(dataStructure.getDataRefAs(idsPath)); + REQUIRE_NOTHROW(dataStructure.getDataRefAs(countsPath)); + const auto& ids = dataStructure.getDataRefAs(idsPath); + const auto& numCells = dataStructure.getDataRefAs(countsPath); + scope.requireExpectedStore(ids); + scope.requireExpectedStore(numCells); + const auto idsFormat = ids.getDataFormat(); + const auto countsFormat = numCells.getDataFormat(); + + const ComputeSurfaceAreaToVolumeFilter filter; + Arguments args; + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SelectedImageGeometryPath_Key, std::make_any(geomPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CellFeatureIdsArrayPath_Key, std::make_any(idsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_NumCellsArrayPath_Key, std::make_any(countsPath)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_CalculateSphericity_Key, std::make_any(calculateSphericity)); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SurfaceAreaVolumeRatioArrayName_Key, std::make_any("SurfaceAreaVolumeRatio")); + args.insertOrAssign(ComputeSurfaceAreaToVolumeFilter::k_SphericityArrayName_Key, std::make_any("Sphericity")); + auto preflight = filter.preflight(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(preflight.outputActions); + auto result = scope.executeFilter(filter, dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(result.result); + + // Count entire plate and tooth faces analytically; do not repeat the production neighbor scan. + const std::array physicalSpacing = {spacing[0], spacing[1], spacing[2]}; + const double area = (normalFaceCounts[0] * physicalSpacing[1] * physicalSpacing[2]) + (normalFaceCounts[1] * physicalSpacing[0] * physicalSpacing[2]) + + (normalFaceCounts[2] * physicalSpacing[0] * physicalSpacing[1]); + const double volume = combCellCount * physicalSpacing[0] * physicalSpacing[1] * physicalSpacing[2]; + CAPTURE(spacing[0], spacing[1], spacing[2], calculateSphericity, area, volume); + REQUIRE_NOTHROW(dataStructure.getDataRefAs(ratioPath)); + const auto& ratios = dataStructure.getDataRefAs(ratioPath); + scope.requireExpectedStore(ratios); + CHECK(ratios.getTupleShape() == featureShape); + CHECK(ratios.getComponentShape() == componentShape); + std::array ratioValues = {}; + auto ratiosRead = ratios.getDataStoreRef().copyIntoBuffer(0, nonstd::span(ratioValues.data(), ratioValues.size())); + SIMPLNX_RESULT_REQUIRE_VALID(ratiosRead); + CHECK(ratioValues[0] == 0.0f); + CHECK(ratioValues[1] == Approx(area / volume).epsilon(1.0e-6)); + if(calculateSphericity) + { + REQUIRE_NOTHROW(dataStructure.getDataRefAs(sphericityPath)); + const auto& sphericity = dataStructure.getDataRefAs(sphericityPath); + scope.requireExpectedStore(sphericity); + CHECK(sphericity.getTupleShape() == featureShape); + CHECK(sphericity.getComponentShape() == componentShape); + std::array values = {}; + auto sphericityRead = sphericity.getDataStoreRef().copyIntoBuffer(0, nonstd::span(values.data(), values.size())); + SIMPLNX_RESULT_REQUIRE_VALID(sphericityRead); + CHECK(values[0] == 0.0f); + const double expected = std::cbrt(std::acos(-1.0)) * std::pow(6.0 * volume, 2.0 / 3.0) / area; + CHECK(values[1] == Approx(expected).epsilon(1.0e-6)); + } + else + { + CHECK(dataStructure.getDataAs(sphericityPath) == nullptr); + } + + std::vector actualLabels(cellCount); + std::array actualCounts = {}; + auto idsRead = ids.getDataStoreRef().copyIntoBuffer(0, nonstd::span(actualLabels.data(), actualLabels.size())); + SIMPLNX_RESULT_REQUIRE_VALID(idsRead); + auto countsRead = numCells.getDataStoreRef().copyIntoBuffer(0, nonstd::span(actualCounts.data(), actualCounts.size())); + SIMPLNX_RESULT_REQUIRE_VALID(countsRead); + CHECK(actualLabels == labels); + CHECK(actualCounts == counts); + CHECK(ids.getTupleShape() == cellShape); + CHECK(ids.getComponentShape() == componentShape); + CHECK(numCells.getTupleShape() == featureShape); + CHECK(numCells.getComponentShape() == componentShape); + CHECK(ids.getDataFormat() == idsFormat); + CHECK(numCells.getDataFormat() == countsFormat); + CHECK(&ids.getDataStoreRef() == idsStorePtr.get()); + CHECK(&numCells.getDataStoreRef() == countsStorePtr.get()); + scope.requireExpectedStore(ids); + scope.requireExpectedStore(numCells); + REQUIRE_NOTHROW(dataStructure.getDataRefAs(geomPath)); + const auto& finalGeom = dataStructure.getDataRefAs(geomPath); + CHECK(finalGeom.getDimensions() == dimensions); + CHECK(finalGeom.getSpacing() == spacing); + CHECK(finalGeom.getOrigin() == origin); + CHECK(cellAttributeMatrixPtr->getShape() == cellShape); + CHECK(featureAttributeMatrixPtr->getShape() == featureShape); + UnitTest::CheckArraysInheritTupleDims(dataStructure); +} +} // namespace surface_area_to_volume_test diff --git a/src/simplnx/Utilities/DataArrayUtilities.cpp b/src/simplnx/Utilities/DataArrayUtilities.cpp index e78a403b85..30f1445b24 100644 --- a/src/simplnx/Utilities/DataArrayUtilities.cpp +++ b/src/simplnx/Utilities/DataArrayUtilities.cpp @@ -1,8 +1,10 @@ #include "DataArrayUtilities.hpp" +#include "simplnx/Common/Result.hpp" #include "simplnx/Common/Types.hpp" #include "simplnx/Common/TypesUtility.hpp" #include "simplnx/Core/Application.hpp" +#include "simplnx/DataStructure/Geometry/IGridGeometry.hpp" #include "simplnx/Filter/Actions/CreateArrayAction.hpp" #include "simplnx/Filter/Actions/CreateNeighborListAction.hpp" #include "simplnx/Filter/Actions/CreateStringArrayAction.hpp" @@ -13,6 +15,7 @@ #include "simplnx/Utilities/StringInterpretationUtilities.hpp" #include "simplnx/Utilities/TemplateHelpers.hpp" +#include #include using namespace nx::core; @@ -141,6 +144,17 @@ struct CreateDefaultValueNeighborListFunctor namespace nx::core { +Result<> ValidateCellTupleCount(const IGridGeometry& geometry, const DataPath& geometryPath, usize tupleCount, const DataPath& dataPath, int32 errorCode) +{ + const usize cellCount = geometry.getNumberOfCells(); + if(tupleCount != cellCount) + { + return MakeErrorResult(errorCode, fmt::format("The selected data '{}' has {} tuples, but the selected grid geometry '{}' has {} cells. Select data with one tuple per geometry cell.", + dataPath.toString(), tupleCount, geometryPath.toString(), cellCount)); + } + return {}; +} + bool IsChildOfAttributeMatrix(const DataStructure& dataStructure, const DataPath& arrayPath, const AttributeMatrix& attributeMatrix) { if(arrayPath.getLength() < 2) diff --git a/src/simplnx/Utilities/DataArrayUtilities.hpp b/src/simplnx/Utilities/DataArrayUtilities.hpp index 4ed4a436c3..487dad7759 100644 --- a/src/simplnx/Utilities/DataArrayUtilities.hpp +++ b/src/simplnx/Utilities/DataArrayUtilities.hpp @@ -6,9 +6,11 @@ #include "simplnx/Common/Bit.hpp" #include "simplnx/Common/Result.hpp" #include "simplnx/Common/ScratchBuffer.hpp" +#include "simplnx/Common/Types.hpp" #include "simplnx/DataStructure/AbstractListStore.hpp" #include "simplnx/DataStructure/AttributeMatrix.hpp" #include "simplnx/DataStructure/DataArray.hpp" +#include "simplnx/DataStructure/DataPath.hpp" #include "simplnx/DataStructure/NeighborList.hpp" #include "simplnx/DataStructure/StringArray.hpp" #include "simplnx/Filter/IFilter.hpp" @@ -48,6 +50,21 @@ */ namespace nx::core { +class IGridGeometry; + +/** + * @brief Checks that data has one tuple per geometry cell. + * @param geometry Supplies the cell count. + * @param geometryPath Identifies the geometry in an error message. + * @param tupleCount Number of tuples in the selected data. + * @param dataPath Identifies the selected data in an error message. + * @param errorCode Error code to return if the counts differ. + * @return Empty valid result if the counts match; otherwise, an error with both paths and counts. + * + * This check compares counts only. It does not read array values or require matching tuple shapes or a common parent. + */ +SIMPLNX_EXPORT Result<> ValidateCellTupleCount(const IGridGeometry& geometry, const DataPath& geometryPath, usize tupleCount, const DataPath& dataPath, int32 errorCode); + /** * @brief Tests whether a path's parent identifies the specified Attribute Matrix. * @param dataStructure Resolves path identifiers. diff --git a/test/DataArrayTest.cpp b/test/DataArrayTest.cpp index 908ead1b8f..1ac964fdda 100644 --- a/test/DataArrayTest.cpp +++ b/test/DataArrayTest.cpp @@ -1,16 +1,21 @@ #include "simplnx/DataStructure/DataArray.hpp" #include "simplnx/Common/Array.hpp" #include "simplnx/Common/ScopeGuard.hpp" +#include "simplnx/Common/Types.hpp" #include "simplnx/Common/TypesUtility.hpp" #include "simplnx/Core/Application.hpp" #include "simplnx/Core/Preferences.hpp" #include "simplnx/DataStructure/AttributeMatrix.hpp" #include "simplnx/DataStructure/DataGroup.hpp" +#include "simplnx/DataStructure/DataPath.hpp" #include "simplnx/DataStructure/DataStore.hpp" #include "simplnx/DataStructure/DataStructure.hpp" #include "simplnx/DataStructure/EmptyDataStore.hpp" #include "simplnx/DataStructure/EmptyListStore.hpp" #include "simplnx/DataStructure/EmptyStringStore.hpp" +#include "simplnx/DataStructure/Geometry/IGridGeometry.hpp" +#include "simplnx/DataStructure/Geometry/ImageGeom.hpp" +#include "simplnx/DataStructure/Geometry/RectGridGeom.hpp" #include "simplnx/DataStructure/IO/Generic/DataIOCollection.hpp" #include "simplnx/DataStructure/IO/Generic/IDataIOManager.hpp" #include "simplnx/DataStructure/IO/Generic/IDataStoreFormatResolver.hpp" @@ -27,6 +32,7 @@ #include "DataStructObserver.hpp" #include +#include #include #include @@ -35,11 +41,59 @@ #include #include #include +#include #include #include using namespace nx::core; +TEST_CASE("DataArrayUtilities::ValidateCellTupleCount", "[simplnx][DataArrayUtilities]") +{ + const bool useRectGrid = GENERATE(false, true); + const bool emptyGeometry = GENERATE(false, true); + const usize tupleCount = GENERATE(usize{0}, usize{1}, usize{59}, usize{60}, usize{61}); + const int32 errorCode = GENERATE(-1234, -5678); + const usize cellCount = emptyGeometry ? 0 : 60; + CAPTURE(useRectGrid, emptyGeometry, tupleCount, errorCode); + + DataStructure dataStructure; + auto* parentPtr = DataGroup::Create(dataStructure, "Geometry Group"); + REQUIRE(parentPtr != nullptr); + const DataPath geometryPath({"Geometry Group", "Selected Geometry"}); + const DataPath dataPath({"Other Geometry", "Cell Data", "Values"}); + IGridGeometry* gridGeomPtr = nullptr; + if(useRectGrid) + { + gridGeomPtr = RectGridGeom::Create(dataStructure, geometryPath.getTargetName(), parentPtr->getId()); + } + else + { + gridGeomPtr = ImageGeom::Create(dataStructure, geometryPath.getTargetName(), parentPtr->getId()); + } + REQUIRE(gridGeomPtr != nullptr); + gridGeomPtr->setDimensions({emptyGeometry ? usize{0} : usize{4}, 3, 5}); + + auto result = ValidateCellTupleCount(*gridGeomPtr, geometryPath, tupleCount, dataPath, errorCode); + if(tupleCount == cellCount) + { + REQUIRE(result.valid()); + } + else + { + REQUIRE(result.invalid()); + REQUIRE(result.errors().size() == 1); + const auto& error = result.errors().front(); + CHECK(error.code == errorCode); + CHECK(error.message.find(geometryPath.toString()) != std::string::npos); + CHECK(error.message.find(dataPath.toString()) != std::string::npos); + CHECK(error.message.find(fmt::format("{} tuples", tupleCount)) != std::string::npos); + CHECK(error.message.find(fmt::format("{} cells", cellCount)) != std::string::npos); + CHECK(error.message.find("Select data with one tuple per geometry cell.") != std::string::npos); + } + CHECK(result.warnings().empty()); + UnitTest::CheckArraysInheritTupleDims(dataStructure); +} + TEST_CASE("DataArrayUtilities::IsChildOfAttributeMatrix", "[simplnx][DataArrayUtilities]") { DataStructure dataStructure;