Skip to content

VV/BUG: Compute Feature Sizes Full V&V and Thread Summation Patch #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

Open
nyoungbq wants to merge 3 commits into
base: develop
Choose a base branch
from
Open

VV/BUG: Compute Feature Sizes Full V&V and Thread Summation Patch #1638

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
35c46a9
Remove Legacy Test
nyoungbq Jun 5, 2026
7d8e90d
provenance file
nyoungbq Jun 5, 2026
aa689d3
- patched thread kahan summation bug
nyoungbq Jun 11, 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
6 changes: 3 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 @@ -344,13 +344,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
74 changes: 1 addition & 73 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,18 +8,11 @@

#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
{
// Geometry Level
Expand Down Expand Up @@ -71,7 +64,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 Down Expand Up @@ -647,71 +640,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
106 changes: 106 additions & 0 deletions src/Plugins/SimplnxCore/vv/ComputeFeatureSizesFilter.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,106 @@
# V&V Report: ComputeFeatureSizes

| | |
|---|---|
| Plugin | `SimplnxCore` |
| SIMPLNX UUID | `c666ee17-ca58-4969-80d0-819986c72485` |
| DREAM3D 6.5.171 equivalent | `FindSizes` — UUID `656f144c-a120-5c3b-bee5-06deab438588` |
| Verified commit | *pending* |
| Status | COMPLETE |
| Sign-off | Nathan Young, June 10th, 2026 |

## At a glance

| Aspect | State |
|---|---|
| Algorithm relationship | **Port** of both `FindSizes::execute()` (ImageGeom) and `FindSizes::findSizesUnstructured()` (RectGridGeom) |
| Oracle | **Class 1 (Analytical)** — all test data inlined, hand-derived |
| Code paths | **14 of 19** exercised; 5 gaps (cancel ×2, INT32_MAX overflow ×2, other-geom no-op ×1) |
| Tests | **9 TEST_CASEs** — all pass |
| External archive | None |
| Deviations | **1 active** (`ComputeFeatureSizes-D1`): Kahan vs naive summation in RectGridGeom — see deviations file |
| Open bugs | **1 open** (`Bug-1`): 2D area formula uses all 3 spacings (= volume) instead of the 2 non-flat spacings |

## Summary

`ComputeFeatureSizesFilter` produces three arrays per feature: `NumElements` (voxel count, `int32`), `Volumes`/`Areas` (float32), and `EquivalentDiameters` (ESD or ECD, float32). For ImageGeom it uses voxel count × voxel size; for RectGridGeom it sums per-cell element sizes per feature. Both paths are parallelized via `ParallelDataAlgorithm` + TBB `tbb::combinable`. All tests use Class 1 oracles (hand-constructed fixtures with first-principles expected values). Source-inspection confirms both geometry paths are ports of legacy `FindSizes`. One precision deviation (`D1`) and one open bug (`Bug-1`) are documented below.

## Algorithm Relationship

**Port** of `FindSizes` (both ImageGeom and RectGridGeom paths).

- **ImageGeom** (`ProcessImageGeom`): voxel count × voxel volume/area → ESD/ECD. Identical to `FindSizes::execute()`. Formulas: `ESD = 2·∛(V / (4π/3))`; `ECD = 2·√(A / π)`.
- **RectGridGeom** (`ProcessRectGridGeom`): per-cell `ΔxΔyΔz` sizes, summed per feature. Port of `FindSizes::findSizesUnstructured()`. SIMPLNX departs in precision: float32 element sizes are promoted to float64 and Kahan summation is applied at two levels (per-thread and post-reduction). See deviation `D1`.

Port-time additions not in DREAM3D 6.5.171:
- **Parallel execution** — `ParallelDataAlgorithm` + `tbb::combinable`; legacy was serial.
- **Execute-time FeatureId bounds check** — `ValidateFeatureIdsToFeatureAttributeMatrixIndexing`; legacy accessed out-of-bounds silently.
- **Preflight dimension guard** — rejects ImageGeom with two or more dimensions equal to 1. This is because a single empty dimension converts the calculation from Volume to Area and "1D" Images would not make sense to get an Area formula.

## Oracle

**Class 1 (Analytical)** — all expected values are derived from first principles and asserted in `test/ComputeFeatureSizesTest.cpp`.

| Fixture | Geometry | Derived quantity |
|---|---|---|
| `Create2DImageDataStructure()` | 5×5×1 ImageGeom, spacing 20.2×0.1×1.0 | `voxelArea = 20.2 × 0.1 × 1.0 = 2.02`; areas = count × 2.02; ECDs from circle formula |
| `Create3DImageDataStructure()` | 5×5×5 ImageGeom, spacing 1.2×0.9×2.1 | `voxelVol = 1.2 × 0.9 × 2.1 = 2.268`; volumes = count × 2.268; ESDs from sphere formula |
| `CreateRectGridDataStructure()` | 4×4×4 non-uniform RectGrid | Per-cell `ΔxΔyΔz` summed per feature; step-by-step trace in provenance sidecar |

Negative tests use degenerate inputs (out-of-bounds FeatureId; degenerate dims) and assert an error result is returned.

Second-engineer review pending: verify hand-derivations for all three fixtures and the tolerance choice (`std::numeric_limits<float32>::epsilon()`).

## Code path coverage

14 of 19 paths exercised. Source: `Algorithms/ComputeFeatureSizes.cpp`.

| # | Path | Exercised by |
|---|---|---|
| 1 | Preflight: `emptyDimCount > 1` → error | `Invalid: Preflight Failure` (4 sub-checks) |
| 2 | Preflight: valid dims → pass | All 6 positive tests |
| 3 | Execute validate: FeatureId > numFeatures → error | `Invalid: Execution Failure` |
| 4 | Execute validate: passes → geometry dispatch | All 6 positive tests |
| 5 | Dispatch: ImageGeom → `ProcessImageGeom` | `Image 2D *`, `Image Stack 3D *` |
| 6 | Dispatch: RectGridGeom → `ProcessRectGridGeom` | `Rectilinear Grid *` |
| 7 | Dispatch: other geometry → no-op | *Not tested. `GeometrySelectionParameter` prevents this at runtime; gap acceptable.* |
| 8 | Image: any dim == 1 → area + ECD | `Image 2D *` |
| 9 | Image: all dims > 1 → volume + ESD | `Image Stack 3D *` |
| 10 | Image: voxelCount > `k_MaxVoxelCount` → error | *Not tested. Requires >2³¹ voxels in one feature; impractical. Gap acceptable.* |
| 11 | Image: `SaveElementSizes = false` | `Image 2D`, `Image Stack 3D` |
| 12 | Image: `SaveElementSizes = true` | `Image 2D with Element Sizes`, `Image Stack 3D with Element Size` |
| 13 | Image: `shouldCancel` in loop → early return | *Not tested. Cancel disregard would cause hangs in any test; low-value gap.* |
| 14 | RectGrid: `findElementSizes` → per-cell volumes available | `Rectilinear Grid *` |
| 15 | RectGrid: Kahan summation in `RectGridSummationImpl` + post-reduction | `Rectilinear Grid *` (non-uniform spacing exercises summation) |
| 16 | RectGrid: voxelCount > `k_MaxVoxelCount` → error | *Not tested. Same rationale as path 10.* |
| 17 | RectGrid: `SaveElementSizes = false` → `deleteElementSizes()` | `Rectilinear Grid` |
| 18 | RectGrid: `SaveElementSizes = true` → element sizes retained | `Rectilinear Grid with Element Size` |
| 19 | RectGrid: `shouldCancel` in loop → early return | *Not tested. Same rationale as path 13.* |

## Test inventory

| Test case | Notes |
|---|---|
| `Valid: Image 2D` | Class 1; 5×5×1, spacing 20.2×0.1×1.0; SaveElementSizes=false |
| `Valid: Image 2D with Element Sizes` | Same fixture; SaveElementSizes=true |
| `Valid: Image Stack 3D` | Class 1; 5×5×5, spacing 1.2×0.9×2.1; SaveElementSizes=false |
| `Valid: Image Stack 3D with Element Size` | Same fixture; SaveElementSizes=true |
| `Valid: Rectilinear Grid` | Class 1; 4×4×4 non-uniform; SaveElementSizes=false; Kahan path exercised |
| `Valid: Rectilinear Grid with Element Size` | Same fixture; SaveElementSizes=true |
| `Invalid: Execution Failure` | FeatureId 10 in a 4-feature AM; asserts execute result invalid |
| `Invalid: Preflight Failure` | 4 degenerate dim configs; asserts preflight result invalid |
| `SIMPL Backwards Compatibility` | `DYNAMIC_SECTION` over SIMPL 6.4 + 6.5; validates UUID + arg-key + param-value decoding |

All 9 TEST_CASEs pass.

## Exemplar archive

None — all fixtures constructed in C++ at test time. Provenance in `vv/provenance/ComputeFeatureSizesFilter.md`.

## Deviations from DREAM3D 6.5.171

**ImageGeom path:** No deviations. Formulas and logic identical to `FindSizes::execute()`. The float64 intermediate in the ESD/ECD computation differs from likely float32 in legacy by ≤1 ULP for typical EBSD spacings; non-material.

**RectGridGeom path — `ComputeFeatureSizes-D1`:** Per-feature volumes differ from `FindSizes::findSizesUnstructured()` output due to two precision improvements in SIMPLNX: (1) element sizes promoted from float32 to float64 before accumulation; (2) Kahan compensated summation applied per-thread and during TBB post-reduction, reducing accumulated error from O(N·ε_float32) to O(ε_float64). SIMPLNX is more accurate. Users migrating from DREAM3D 6.5.171 should expect small shifts in per-feature volumes for RectGridGeom; largest for features with many cells on grids with high cell-volume variation.

Full entry in `vv/deviations/ComputeFeatureSizesFilter.md`.
47 changes: 47 additions & 0 deletions src/Plugins/SimplnxCore/vv/deviations/ComputeFeatureSizesFilter.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
# Deviations from DREAM3D 6.5.171: ComputeFeatureSizes

This file lists every documented behavioral difference between this SIMPLNX filter and its DREAM3D 6.5.171 equivalent.

Entries are referenced by stable ID (`ComputeFeatureSizes-D<N>`) from the V&V report and from public migration guidance. The ID is stable across renames; the Filter UUID field is the permanent cross-reference anchor.

---

## ImageGeom path — no deviations

Source-inspection comparison of `ProcessImageGeom` against DREAM3D 6.5.171 `FindSizes::execute()` confirmed identical algorithmic structure:

- Per-feature voxel count: identical parallel accumulation logic (serial in legacy).
- 3D volume formula: `volume = voxelCount × voxelVolume`; ESD: `2·∛(volume / (4π/3))` — both use the standard sphere formula.
- 2D area formula (any dim == 1): `area = voxelCount × voxelArea`; ECD: `2·√(area / π)` — both use the standard circle formula.

**Precision non-deviation (not flagged):** SIMPLNX promotes `voxelVolume` to `float64` and evaluates `cbrt` / `sqrt` in `float64` before casting to `float32` for storage. The legacy filter likely used `float32` throughout. The difference is within ~1 ULP of `float32` for typical EBSD spacings and is non-material for downstream morphological statistics.

---

## ComputeFeatureSizes-D1

| Field | Value |
|---|---|
| **Deviation ID** | `ComputeFeatureSizes-D1` |
| **Filter UUID** | `c666ee17-ca58-4969-80d0-819986c72485` |
| **Status** | Active — precision improvement; SIMPLNX output is more accurate |

**Symptom:** Per-feature volumes for features in `RectGridGeom` geometries differ from DREAM3D 6.5.171 `FindSizes::findSizesUnstructured()` output. The discrepancy grows with feature size and with the variation in cell volumes across the grid, scaling approximately as `O(N · ε_float32)` for the legacy result vs `O(ε_float64)` for the SIMPLNX result, where N is the number of cells in the feature.

**Root cause:** Precision — two compounding improvements in SIMPLNX's `ProcessRectGridGeom` over the legacy `findSizesUnstructured`:

1. **float32 → float64 promotion.** Element sizes are stored as `float32` (the output of `findElementSizes`). The legacy `findSizesUnstructured` accumulated these directly in `float32`, so each per-cell addition carried a relative error of ~`ε_float32 ≈ 1.2×10⁻⁷`. SIMPLNX promotes each element size to `float64` before accumulation (`static_cast<float64>(m_ElemSizes.getValue(voxelIdx))`), reducing the per-operation rounding error to `ε_float64 ≈ 2.2×10⁻¹⁶`.

2. **Kahan summation.** SIMPLNX applies Kahan compensated summation at two levels:
- *Per-thread, per-cell* in `RectGridSummationImpl::convert()` (lines 274–283): standard Kahan with `y = elemSize - c`, `t = sum + y`, `c′ = (t - sum) - y`. For a feature spanning N cells, naive float64 accumulation still has O(N · ε_float64) error; Kahan reduces this to O(ε_float64) by carrying a compensation term `c` that recovers the low-order bits lost in each `sum + y` operation.
- *Post-reduction, per thread-local vector* in `threadLocalVolumes.combine_each()` (lines 341–358): the same Kahan scheme is applied when combining the TBB thread-local partial sums into the final per-feature volume.

For a non-uniform rectilinear grid where cell volumes span several orders of magnitude (e.g., a grid with fine near-surface resolution and coarse interior), naive summation of N cells can lose ~`log₂(V_max / V_min)` bits of precision in the smaller cell contributions. Kahan summation recovers those bits regardless of N or the dynamic range of cell volumes.

**Affected users:** Any workflow that runs `ComputeFeatureSizes` on a `RectGridGeom` and then compares per-feature volume or ESD values against DREAM3D 6.5.171 output. The deviation is largest for large features (high N) on grids with high volume variation. On uniform RectGrids (all cell volumes equal), Kahan has no practical effect and the only deviation is the float32→float64 promotion.

**Recommendation:** Trust SIMPLNX. The legacy `findSizesUnstructured` accumulated unnecessary floating-point error that grew with feature size and grid non-uniformity. The SIMPLNX result is strictly more accurate. Users migrating pipelines should expect small positive or negative shifts in per-feature volumes; ESD values are affected proportionally through the `cbrt` formula.

---

*If a future comparison run against DREAM3D 6.5.171 output reveals additional deviations, add them here as `ComputeFeatureSizes-D2` etc.*
Loading
Loading