Skip to content

VV/BUG: Compute Feature Sizes Full V&V + precision and 2D area fixes #1638

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
imikejackson merged 9 commits into from
Jul 29, 2026
Merged

VV/BUG: Compute Feature Sizes Full V&V + precision and 2D area fixes #1638

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
68bec0a
Remove Legacy Test
nyoungbq Jun 5, 2026
7c7c9cd
provenance file
nyoungbq Jun 5, 2026
9f5c246
- patched thread kahan summation bug
nyoungbq Jun 11, 2026
f3a3711
DOC: ComputeFeatureSizes V&V — add D2 (ESD evaluation) and A/B verifi…
imikejackson Jun 27, 2026
6e70ccd
BUG: Fix ComputeFeatureSizes 2D area formula; review fixes
imikejackson Jul 8, 2026
0df5f2c
TEST: Use relative tolerance for ComputeFeatureSizes oracle assertions
imikejackson Jul 16, 2026
d9ff3d3
DOC: Correct ComputeFeatureSizes V&V claims against verified sources
imikejackson Jul 16, 2026
ec9086d
Revert Code change to the PR #1590 conventions for 2D Image Geometries.
imikejackson Jul 17, 2026
5ca3f12
TEST: Add missing [!shouldfail] tag to 2D flat-spacing characterizati…
imikejackson Jul 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 19 additions & 3 deletions src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureSizes.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,12 @@ Result<> ProcessImageGeom(ImageGeom& imageGeom, Float32AbstractDataStore& volume
FeatureVoxelCountsT threadLocalVoxelCounts([numFeatures] { return std::vector<uint64>(numFeatures, 0); });
ParallelDataAlgorithm dataAlg;
dataAlg.setRange(0, dims[2]);
// The worker reads featureIds concurrently across slices; per the project thread-safety policy a
// DataStore is only safe for concurrent access when it is resident in memory, so gate parallelization
// on that (an out-of-core store falls back to serial execution).
IParallelAlgorithm::AlgorithmStores algStores;
algStores.push_back(&featureIds);
dataAlg.requireStoresInMemory(algStores);
dataAlg.execute(ImageSummationImpl(threadLocalVoxelCounts, dims, featureIds, shouldCancel));

if(shouldCancel)
Expand Down Expand Up @@ -142,6 +148,11 @@ Result<> ProcessImageGeom(ImageGeom& imageGeom, Float32AbstractDataStore& volume
* For these two cases the following code would BREAK, so do not enable.
**/

// OPEN DESIGN QUESTION: this includes the flat dimension's spacing, matching the PR #1590
// "slab" convention used by ImageGeom::findElementSizes, but it diverges from legacy DREAM3D
// 6.5.171 (FindSizes::findSizesImage uses only the two non-flat resolutions) whenever the flat
// dimension's spacing != 1. See the "[2DFlatSpacing]" test case, which characterizes the
// divergence, and the ComputeFeatureSizesFilter V&V deviations entry.
// Calculate the area of a single voxel
const float64 voxelArea = static_cast<float64>(spacing[0]) * static_cast<float64>(spacing[1]) * static_cast<float64>(spacing[2]);

Expand Down Expand Up @@ -326,6 +337,11 @@ Result<> ProcessRectGridGeom(RectGridGeom& rectGridGeom, Float32AbstractDataStor
FeatureVolumesT threadLocalVolumes([numFeatures] { return std::vector<float64>(numFeatures, 0); });
ParallelDataAlgorithm dataAlg;
dataAlg.setRange(0, dims[2]);
// The worker reads featureIds concurrently across slices; gate parallelization on the store being
// resident in memory (see the ProcessImageGeom note; project thread-safety policy).
IParallelAlgorithm::AlgorithmStores algStores;
algStores.push_back(&featureIds);
dataAlg.requireStoresInMemory(algStores);
dataAlg.execute(RectGridSummationImpl(threadLocalVoxelCounts, threadLocalVolumes, dims, featureIds, elemSizes, shouldCancel));

if(shouldCancel)
Expand All @@ -344,13 +360,13 @@ Result<> ProcessRectGridGeom(RectGridGeom& rectGridGeom, Float32AbstractDataStor
// Use Kahan summation to determine overall volume

// Attempt to recover low order into the value. The first instance is 0
const float64 value = featureVolumes[featureIdx] - featureCompensators[featureIdx];
const float64 value = localVolumes[featureIdx] - featureCompensators[featureIdx];

// low order may be lost
const float64 volSum = localVolumes[featureIdx] + value;
const float64 volSum = featureVolumes[featureIdx] + value;

// recover and cache low order
featureCompensators[featureIdx] = (volSum - localVolumes[featureIdx]) - value;
featureCompensators[featureIdx] = (volSum - featureVolumes[featureIdx]) - value;

// store volumes
featureVolumes[featureIdx] = volSum;
Expand Down
180 changes: 87 additions & 93 deletions src/Plugins/SimplnxCore/test/ComputeFeatureSizesTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,20 +8,18 @@

#include <catch2/catch.hpp>
#include <filesystem>
#include <fstream>

using namespace nx::core;

namespace fs = std::filesystem;

namespace LegacyTest
{
const std::string k_Volumes("Volumes");
const std::string k_EquivalentDiameters("EquivalentDiameters");
} // namespace LegacyTest

namespace Test
{
// Relative tolerance for comparing float32 results against the hand-derived oracle values. A few
// float32 ULPs (~1.2e-7 each) of slack so the pins survive platform and TBB reduction-order
// differences without demanding bit-exact equality; anything beyond this is a real deviation.
constexpr float64 k_RelativeTolerance = 1.0e-6;

// Geometry Level
const std::string k_ImageGeomName = "Image";
const DataPath k_ImageGeomPath = DataPath({k_ImageGeomName});
Expand Down Expand Up @@ -71,7 +69,7 @@ DataStructure Create2DImageDataStructure()
1, 1, 1, 1, 1,
1, 1, 1, 3, 3,
3, 3, 1, 1, 3,
3, 3, 3, 3, 3,
3, 3, 3, 3, 3
};
// clang-format on

Expand All @@ -94,14 +92,14 @@ void Validate2DImageDataStructure(const DataStructure& dataStructure)
REQUIRE(numElements.getValue(3) == 13);
// areas: 0.0 22.22 2.02 26.26
const auto& areas = dataStructure.getDataRefAs<Float32Array>(k_VolumesPath);
REQUIRE(std::abs(areas.getValue(1) - 22.220001f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(areas.getValue(2) - 2.0200002f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(areas.getValue(3) - 26.260002f) < std::numeric_limits<float32>::epsilon());
REQUIRE(areas.getValue(1) == Approx(22.220001f).epsilon(k_RelativeTolerance));
REQUIRE(areas.getValue(2) == Approx(2.0200002f).epsilon(k_RelativeTolerance));
REQUIRE(areas.getValue(3) == Approx(26.260002f).epsilon(k_RelativeTolerance));
// eqDiameters: 0.0 5.318964 1.603728 5.78232
const auto& equivalentDiameters = dataStructure.getDataRefAs<Float32Array>(k_EquivalentDiametersPath);
REQUIRE(std::abs(equivalentDiameters.getValue(1) - 5.3189644f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(equivalentDiameters.getValue(2) - 1.60372818f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(equivalentDiameters.getValue(3) - 5.7823243f) < std::numeric_limits<float32>::epsilon());
REQUIRE(equivalentDiameters.getValue(1) == Approx(5.3189644f).epsilon(k_RelativeTolerance));
REQUIRE(equivalentDiameters.getValue(2) == Approx(1.60372818f).epsilon(k_RelativeTolerance));
REQUIRE(equivalentDiameters.getValue(3) == Approx(5.7823243f).epsilon(k_RelativeTolerance));
}

DataStructure Create3DImageDataStructure()
Expand Down Expand Up @@ -178,14 +176,14 @@ void Validate3DImageDataStructure(const DataStructure& dataStructure)
REQUIRE(numElements.getValue(3) == 23);
// volumes: 0.0 165.564 65.772 52.164
const auto& volumes = dataStructure.getDataRefAs<Float32Array>(k_VolumesPath);
REQUIRE(std::abs(volumes.getValue(1) - 165.564f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(volumes.getValue(2) - 65.771995f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(volumes.getValue(3) - 52.163997f) < std::numeric_limits<float32>::epsilon());
REQUIRE(volumes.getValue(1) == Approx(165.564f).epsilon(k_RelativeTolerance));
REQUIRE(volumes.getValue(2) == Approx(65.771995f).epsilon(k_RelativeTolerance));
REQUIRE(volumes.getValue(3) == Approx(52.163997f).epsilon(k_RelativeTolerance));
// eqDiameters: 0.0 6.81275 5.00819 4.63579
const auto& equivalentDiameters = dataStructure.getDataRefAs<Float32Array>(k_EquivalentDiametersPath);
REQUIRE(std::abs(equivalentDiameters.getValue(1) - 6.8127493f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(equivalentDiameters.getValue(2) - 5.0081901f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(equivalentDiameters.getValue(3) - 4.6357936f) < std::numeric_limits<float32>::epsilon());
REQUIRE(equivalentDiameters.getValue(1) == Approx(6.8127493f).epsilon(k_RelativeTolerance));
REQUIRE(equivalentDiameters.getValue(2) == Approx(5.0081901f).epsilon(k_RelativeTolerance));
REQUIRE(equivalentDiameters.getValue(3) == Approx(4.6357936f).epsilon(k_RelativeTolerance));
}

DataStructure CreateRectGridDataStructure()
Expand Down Expand Up @@ -232,8 +230,8 @@ DataStructure CreateRectGridDataStructure()
// clang-format off
// Expected Outputs:
// numElements: 0 39 15 10
// volumes: 0.0 2358.834 352.462 15.104
// eqDiameters: 0.0 16.516 8.764 3.0994
// volumes: 0.0 2362.434 352.462 15.104
// eqDiameters: 0.0 16.5242 8.76404 3.06689
const std::array<uint8, 64> featureIdsArray = {
1, 2, 2, 2,
1, 1, 1, 1,
Expand Down Expand Up @@ -275,14 +273,14 @@ void ValidateRectGridDataStructure(const DataStructure& dataStructure)
REQUIRE(numElements.getValue(3) == 10);
// volumes: 0.0 2362.434 352.462 15.104
const auto& volumes = dataStructure.getDataRefAs<Float32Array>(k_VolumesPath);
REQUIRE(std::abs(volumes.getValue(1) - 2362.43384f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(volumes.getValue(2) - 352.461884f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(volumes.getValue(3) - 15.1039925f) < std::numeric_limits<float32>::epsilon());
REQUIRE(volumes.getValue(1) == Approx(2362.43384f).epsilon(k_RelativeTolerance));
REQUIRE(volumes.getValue(2) == Approx(352.461884f).epsilon(k_RelativeTolerance));
REQUIRE(volumes.getValue(3) == Approx(15.1039925f).epsilon(k_RelativeTolerance));
// eqDiameters: 0.0 16.5242 8.76404 3.06689
const auto& equivalentDiameters = dataStructure.getDataRefAs<Float32Array>(k_EquivalentDiametersPath);
REQUIRE(std::abs(equivalentDiameters.getValue(1) - 16.5241966f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(equivalentDiameters.getValue(2) - 8.7640428f) < std::numeric_limits<float32>::epsilon());
REQUIRE(std::abs(equivalentDiameters.getValue(3) - 3.0668866f) < std::numeric_limits<float32>::epsilon());
REQUIRE(equivalentDiameters.getValue(1) == Approx(16.5241966f).epsilon(k_RelativeTolerance));
REQUIRE(equivalentDiameters.getValue(2) == Approx(8.7640428f).epsilon(k_RelativeTolerance));
REQUIRE(equivalentDiameters.getValue(3) == Approx(3.0668866f).epsilon(k_RelativeTolerance));
}
} // namespace Test

Expand Down Expand Up @@ -326,6 +324,67 @@ TEST_CASE("SimplnxCore::ComputeFeatureSizes: Valid: Image 2D", "[SimplnxCore][Co
UnitTest::CheckArraysInheritTupleDims(dataStructure);
}

// Characterization test for an OPEN BUG / pending design decision on the 2D area formula.
//
// The assertions below encode the legacy DREAM3D 6.5.171 semantics (FindSizes::findSizesImage):
// a 2D area is the product of the two NON-flat spacings only, so the single feature's area is
// 4 voxels * (2.0 * 3.0) = 24.0. The current implementation follows the PR #1590 "slab" convention
// (matching ImageGeom::findElementSizes) and multiplies ALL THREE spacings, yielding 120.0 whenever
// the flat dimension's spacing != 1. The existing "Valid: Image 2D" fixture uses a flat-dimension
// spacing of exactly 1.0, which is why it cannot see the divergence; here the flat spacing is
// deliberately 5.0, and all three flat orientations are GENERATEd so no axis-specific special case
// can satisfy it.
//
// While the #1590 convention is in place, the assertions below fail; the [!shouldfail] tag makes
// Catch2 report that as an expected failure (so CI stays green) and flips to a hard failure the
// moment the divergence is resolved, forcing this test to be updated alongside the design decision.
TEST_CASE("SimplnxCore::ComputeFeatureSizes: 2D area excludes the flat-dimension spacing", "[SimplnxCore][ComputeFeatureSizes][2DFlatSpacing][!shouldfail]")
{
auto [label, dims, spacing] = GENERATE(std::make_tuple("flat Z", SizeVec3{2, 2, 1}, FloatVec3{2.0f, 3.0f, 5.0f}), std::make_tuple("flat X", SizeVec3{1, 2, 2}, FloatVec3{5.0f, 2.0f, 3.0f}),
std::make_tuple("flat Y", SizeVec3{2, 1, 2}, FloatVec3{2.0f, 5.0f, 3.0f}));

DYNAMIC_SECTION(label)
{
// Four cells, all feature 1 (feature 0 unused). Non-flat spacings are always 2.0 and 3.0, so the
// per-voxel area is 6.0, and the single feature's area is 4 * 6.0 = 24.0 regardless of which axis is
// flat. The flat dimension's spacing is 5.0 and must NOT appear in the product.
DataStructure dataStructure;
auto* imageGeom = ImageGeom::Create(dataStructure, Test::k_ImageGeomName);
imageGeom->setSpacing(spacing);
imageGeom->setOrigin(FloatVec3{0.0f, 0.0f, 0.0f});
imageGeom->setDimensions(dims);

const ShapeType tupleShape{dims[2], dims[1], dims[0]};
auto* cellData = AttributeMatrix::Create(dataStructure, Test::k_CellAMName, tupleShape, imageGeom->getId());
imageGeom->setCellData(*cellData);
auto* featureIds = Int32Array::CreateWithStore<Int32DataStore>(dataStructure, Test::k_FeatureIdsName, cellData->getShape(), ShapeType{1}, cellData->getId());
featureIds->fill(1);
AttributeMatrix::Create(dataStructure, Test::k_FeatureAMName, ShapeType{2}, imageGeom->getId());

ComputeFeatureSizesFilter filter;
Arguments args;
args.insert(ComputeFeatureSizesFilter::k_GeometryPath_Key, std::make_any<DataPath>(Test::k_ImageGeomPath));
args.insert(ComputeFeatureSizesFilter::k_SaveElementSizes_Key, std::make_any<bool>(false));
args.insert(ComputeFeatureSizesFilter::k_CellFeatureIdsArrayPath_Key, std::make_any<DataPath>(Test::k_FeatureIdsPath));
args.insert(ComputeFeatureSizesFilter::k_CellFeatureAttributeMatrixPath_Key, std::make_any<DataPath>(Test::k_FeatureAMPath));
args.insert(ComputeFeatureSizesFilter::k_VolumesName_Key, std::make_any<std::string>(Test::k_VolumesName));
args.insert(ComputeFeatureSizesFilter::k_EquivalentDiametersName_Key, std::make_any<std::string>(Test::k_EquivalentDiametersName));
args.insert(ComputeFeatureSizesFilter::k_NumElementsName_Key, std::make_any<std::string>(Test::k_NumElementsName));

auto preflightResult = filter.preflight(dataStructure, args);
SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions);
auto executeResult = filter.execute(dataStructure, args);
SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result);

const auto& areas = dataStructure.getDataRefAs<Float32Array>(Test::k_VolumesPath);
// Legacy 6.5.171 expectation: 4 voxels * (2.0 * 3.0) = 24.0. The current #1590 slab convention
// produces 4 * (2.0 * 3.0 * 5.0) = 120.0, so this assertion fails (expected, see [!shouldfail]).
REQUIRE(areas.getValue(1) == Approx(24.0f));

UnitTest::CheckArraysInheritTupleDims(dataStructure);
}
}

TEST_CASE("SimplnxCore::ComputeFeatureSizes: Valid: Image 2D with Element Sizes", "[SimplnxCore][ComputeFeatureSizes]")
{
DataStructure dataStructure = Test::Create2DImageDataStructure();
Expand Down Expand Up @@ -647,71 +706,6 @@ TEST_CASE("SimplnxCore::ComputeFeatureSizes: Invalid: Preflight Failure", "[Simp
#endif
}

TEST_CASE("SimplnxCore::ComputeFeatureSizes: Legacy: Small IN100 Test", "[SimplnxCore][ComputeFeatureSizes]")
{
UnitTest::LoadPlugins();

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
auto baseDataFilePath = fs::path(fmt::format("{}/6_6_stats_test_v2.dream3d", unit_test::k_TestFilesDir));
DataStructure dataStructure = UnitTest::LoadDataStructure(baseDataFilePath);
DataPath smallIn100Group({Constants::k_DataContainer});
DataPath cellDataPath = smallIn100Group.createChildPath(Constants::k_CellData);
DataPath cellPhasesPath = cellDataPath.createChildPath(Constants::k_Phases);
DataPath featureIdsPath = cellDataPath.createChildPath(Constants::k_FeatureIds);
DataPath featureGroup = smallIn100Group.createChildPath(Constants::k_CellFeatureData);
std::string volumesName = "computed_volumes";
std::string numElementsName = "computed_NumElements";
std::string EquivalentDiametersName = "computed_EquivalentDiameters";

std::vector<std::string> featureNames = {LegacyTest::k_Volumes, LegacyTest::k_EquivalentDiameters, Constants::k_NumElements};

{
ComputeFeatureSizesFilter filter;
Arguments args;

args.insert(ComputeFeatureSizesFilter::k_GeometryPath_Key, std::make_any<DataPath>(smallIn100Group));
args.insert(ComputeFeatureSizesFilter::k_SaveElementSizes_Key, std::make_any<bool>(false));
args.insert(ComputeFeatureSizesFilter::k_CellFeatureIdsArrayPath_Key, std::make_any<DataPath>(featureIdsPath));
args.insert(ComputeFeatureSizesFilter::k_CellFeatureAttributeMatrixPath_Key, std::make_any<DataPath>(featureGroup));
args.insert(ComputeFeatureSizesFilter::k_VolumesName_Key, std::make_any<std::string>(volumesName));
args.insert(ComputeFeatureSizesFilter::k_EquivalentDiametersName_Key, std::make_any<std::string>(EquivalentDiametersName));
args.insert(ComputeFeatureSizesFilter::k_NumElementsName_Key, std::make_any<std::string>(numElementsName));

// Preflight the filter and check result
auto preflightResult = filter.preflight(dataStructure, args);
SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions);

// Execute the filter and check the result
auto executeResult = filter.execute(dataStructure, args);
SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result);
}

// Compare Outputs
{
DataPath exemplaryDataPath = featureGroup.createChildPath(LegacyTest::k_Volumes);
UnitTest::CompareArrays<float32>(dataStructure, exemplaryDataPath, featureGroup.createChildPath(volumesName));
}

{
DataPath exemplaryDataPath = featureGroup.createChildPath(LegacyTest::k_EquivalentDiameters);
UnitTest::CompareArrays<float32>(dataStructure, exemplaryDataPath, featureGroup.createChildPath(EquivalentDiametersName));
}

{
DataPath exemplaryDataPath = featureGroup.createChildPath(Constants::k_NumElements);
UnitTest::CompareArrays<int32>(dataStructure, exemplaryDataPath, featureGroup.createChildPath(numElementsName));
}

// Write the DataStructure out to the file system
#ifdef SIMPLNX_WRITE_TEST_OUTPUT
WriteTestDataStructure(dataStructure, fs::path(fmt::format("{}/calculate_feature_sizes/legacy_test.dream3d", unit_test::k_BinaryTestOutputDir)));
#endif

UnitTest::CheckArraysInheritTupleDims(dataStructure);
}

TEST_CASE("SimplnxCore::ComputeFeatureSizesFilter: SIMPL Backwards Compatibility", "[SimplnxCore][ComputeFeatureSizesFilter][BackwardsCompatibility]")
{
auto app = Application::GetOrCreateInstance();
Expand Down
Loading
Loading