BUG: C11: float32 accumulation loses precision in 16 filters #1871

Summary

Many filters count or sum values in float32 variables. A float32 has a 24-bit significand. A counter incremented with ++ or += 1.0f stops at 16,777,216 (2^24), and a long sum loses precision. Large Features, large slices, large ranges or coordinates far from the origin then give wrong output without an error. In one case a loop variable stops advancing and the filter hangs.

As of develop 419f9db (2026-10-21) there are 16 open findings in 16 filters: 1 CRASH (hang), 14 WRONG_DATA and 1 MINOR. Each finding is a candidate. It was found and re-checked by reading code, not by running it.

Mechanism

float32 can store every integer only up to 2^24. Then 16777216.0f + 1.0f == 16777216.0f, so a counter stops. A sum of many small values stops growing when each value is smaller than half the spacing between floats at the current total. The forms are:

• Counts in float32. Voxel or sample counters are float32 and are incremented one at a time. A Feature with more than 16,777,216 cells (for example a matrix Feature in a 300^3 Image Geometry) gets a capped count, so volumes, fractions and means are wrong.
• Sums in float32. Centroids, means and threshold sums use float32 accumulators.
• Intermediate math in float32. Histogram increments and bin indices, color-map normalization and triangle volume determinants use float32 for int32, int64 or float64 input, or for coordinates far from the origin (cancellation).
• Loop variable in float32. lineY += lineSpacing stops advancing when the spacing is smaller than half the float spacing at lineY, and the loop never ends.

Example (ComputeBoundaryElementFractions). Both counters are float32. They are incremented with ++ for every cell (lines 70 and 73), so they stop at 2^24 and the fraction is wrong for large Features:

  std::vector<float32> surfVoxCounts(numFeatures, 0);
  std::vector<float32> voxCounts(numFeatures, 0);

Source: src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeBoundaryElementFractions.cpp#L43-L44

Fixed form. PR #1817 (commit 01809c2) changed ComputeSurfaceAreaToVolume to compute and accumulate in float64 and to convert to float32 only when it writes the output:

  const float64 voxelVol = static_cast<float64>(spacing[0]) * spacing[1] * spacing[2];
  const float64 xyFaceArea = static_cast<float64>(spacing[0]) * spacing[1];
  const float64 yzFaceArea = static_cast<float64>(spacing[1]) * spacing[2];
  const float64 zxFaceArea = static_cast<float64>(spacing[2]) * spacing[0];

  // Multiple cells contribute to one feature, so accumulation stays local until
  // the cell scan is complete.
  std::vector<float64> featureSurfaceArea(numFeatures, 0.0);

Source: src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeSurfaceAreaToVolumeDirect.cpp#L53-L60

Rule

• Count in an integer type (usize or uint64). Sum and do intermediate math in float64. Convert to the output type only when you write the output.
• Loops over a real-valued range use an integer step counter k and compute the value as start + k * step in float64.
• Moving accumulators to float64 changes exemplar outputs. Treat this class as V&V-style work: prove each change with an independent oracle. Do NOT regenerate exemplars just to absorb the change. If an exemplar comparison fails after the fix, report the difference and the oracle evidence in the pull request, and let the project lead decide about the exemplar.
• Checks belong in preflight unless they need DataArray values.

How to fix

1. Change each float32 counter to usize or uint64. Compute fractions and means in float64, then cast the result.
2. Change each float32 sum to float64. If the order of summation matters for exact results, keep the existing order.
3. In src/simplnx/Utilities/HistogramUtilities.hpp, compute CalculateIncrement, FillBinRanges and CalculateBin in float64 (lines 21-24 and 84-90 at 419f9db). The per-feature summation at line 395 is also vector<float32>. These helpers are shared, so one change fixes ComputeArrayHistogram, ComputeArrayHistogramByFeature and part of ComputeArrayStatistics. Check every caller.
4. In src/simplnx/Utilities/ColorTableUtilities.hpp line 151, divide in float64, then cast the [0,1] result to float32.
5. In src/simplnx/Utilities/Meshing/TriangleUtilities.cpp (FindTriangleVolume, lines 922-938), compute the determinant in float64 relative to a reference point near the Feature, and accumulate in float64.
6. In CreateAMScanPaths.cpp line 190, iterate an integer k and compute lineY = startLineY + k * lineSpacing in float64. Add a cancel check in the loop.
7. Fix every copy of the code in the same filter: Direct and Scanline algorithms, 2D and 3D paths, and in-core and out-of-core paths (the checklist names the other lines).
8. New error codes (if any): use the filter's existing error-code range. Codes are more negative than -10000, and each filter uses one range of 100 codes. Check uniqueness with git grep -nE -- "-NNNN([^0-9]|$)" (not \b) on develop and on open pull request refs.

How to test

Write a unit test that fails on develop first. Then fix the code.

• Counter saturation: a test with more than 2^24 elements in one Feature or slice. A uniform input (for example a 257 x 257 x 257 Image Geometry with Feature Id 1 everywhere, 16,974,593 cells) has an exact oracle: the count is the number of cells, and the mean of a constant is the constant. If the memory cost is too high for the normal test suite, test the extracted helper with a synthetic stream of 2^24 + N values, or mark the test as a long test. Write the choice in the pull request.
• Sum and math precision: use inputs with a closed-form answer (for example vertices offset by 1e5 to 1e6 from the origin around a unit cube, whose volume is 1 and whose centroid is known). Compare with a tolerance that is derived from float64, not from the old output.
• Histogram: an int32 DataArray {0, 100000000} with 10 bins. The oracle is the bin edges from integer arithmetic, and the maximum value is in the last bin.
• Hang: an Edge Geometry with coordinates near 1e5 and Hatch Spacing < 0.004. Assert that the filter finishes and that the line count equals floor((maxY - startY) / spacing) + 1.
• If the filter has Direct and Scanline algorithms, test both through the existing algorithm-selection test scope.
• Value checks need an independent expected result (oracle), not the output of the code under test.

Conventions

• One pull request for the whole class. One commit per filter.
• Commit prefix BUG:. Add a Signed-off-by line. Do not add AI attribution.
• Every fix starts with a unit test that fails on develop. If you cannot reproduce a finding, write that in the pull request and on this issue with the reason. Do not force a fix.
• Value checks need an independent expected result (oracle).
• Run ctest serially. Never use -j; the tests depend on each other.
• Format with clang-format 19.1.1 and the repository .clang-format.
• Do not add per-filter documentation sentences about input requirements. The generated parameter table documents them.
• If a fix touches a filter that has a V&V report (src/Plugins/<Plugin>/vv/<Filter>.md), update every file of that filter's V&V package (report, deviations, provenance, test tables). V&V documents reference DREAM3D 6.5.171 only.
• Rebase on develop. Never merge develop into the branch.

Checklist

• CRASH · CreateAMScanPathsFilter · CreateAMScanPaths.cpp:190 — float lineY += lineSpacing: with coordinates near 1e5 and Hatch Spacing < 0.004, lineY stops advancing and the loop never ends (no cancel check). Otherwise the spacing is quantized. CreateAMScanPathsFilter#6
• WRONG_DATA · AlignGeometriesFilter · AlignGeometries.cpp:205 — With Alignment Type = Centroid, vertex centroids are summed in float32 (also lines 220-225, 239-244, 258-263), so a mesh with millions of vertices far from the origin gets a slightly wrong translation. AlignGeometriesFilter#5
• WRONG_DATA · AlignSectionsFeatureCentroidFilter · AlignSectionsFeatureCentroid.cpp:105 — xCentroid/yCentroid are float32 sums over the slice (also the out-of-core path near line 322), so on a 2000 x 2000 slice with a nearly full mask the shift can change by one cell. AlignSectionsFeatureCentroidFilter#6
• WRONG_DATA · AlignSectionsMisorientationFilter · AlignSectionsMisorientation.cpp:137 — disorientation and count are float32 and incremented with ++ (also lines 260-261 and 522-523), so they stop at 2^24 and bias the mismatch fraction on XY slices of more than about 268M cells. AlignSectionsMisorientationFilter#2
• WRONG_DATA · AlignSectionsMutualInformationFilter · AlignSectionsMutualInformation.cpp:371 — The histograms (lines 309-311) and count (lines 371 and 513) are float32 and incremented with ++, so they stop at 2^24 and skew the mutual information on very large slices. AlignSectionsMutualInformationFilter#4
• WRONG_DATA · ComputeArrayHistogramByFeatureFilter · HistogramUtilities.hpp:23 — CalculateIncrement returns float32 and FillBinRanges/CalculateBin use it for int32, int64 and float64 too, so for values over [0, 2^30) the bin edges and the assigned bins disagree near the edges. ComputeArrayHistogramByFeatureFilter#11
• WRONG_DATA · ComputeArrayHistogramFilter · HistogramUtilities.hpp:23 — Same float32 increment and bin index (lines 21-24, 88-89): for an int32 DataArray {0, 100000000} with 10 bins, the maximum lands in overflow and the stored upper edges are wrong. ComputeArrayHistogramFilter#3
• WRONG_DATA · ComputeArrayStatisticsFilter · ComputeArrayStatistics.cpp:2036 — Per-group sums are float32 (line 1907; also HistogramUtilities.hpp:395), so with Compute Statistics Per Feature/Ensemble on and a group of more than 2^24 tuples the Mean is wrong (0.62 instead of 1.0). ComputeArrayStatisticsFilter#4
• WRONG_DATA · ComputeBoundaryElementFractionsFilter · ComputeBoundaryElementFractions.cpp:43 — Cell counts are vector<float32> incremented with ++, so for a Feature of more than 16,777,216 cells the count stops and the fraction is overestimated. ComputeBoundaryElementFractionsFilter#3
• WRONG_DATA · ComputeBoundingBoxStatsFilter · ComputeBoundingBoxStatsDirect.cpp:588 — float32 input is summed into a float32 accumulator, so a box of more than about 1e7 cells (for example all 0.1) drifts. Scanline lines 593 and 651 are the same. ComputeBoundingBoxStatsFilter#6
• WRONG_DATA · ComputeLargestCrossSectionsFilter · ComputeLargestCrossSectionsDirect.cpp:199 — featureCounts is vector<float32> incremented with ++, so it stops at 2^24 per plane (for example 4200 x 4200 XY planes). Scanline line 46 is the same. ComputeLargestCrossSectionsFilter#2
• WRONG_DATA · ComputeShapesFilter · ComputeShapes.cpp:312 — featureVoxelCounts[gnum] += 1.0f on a float32 vector stops at 2^24, which caps Volumes and corrupts Omega3s for large Features. Lines 470 and 515 (second path) are the same. ComputeShapesFilter#3
• WRONG_DATA · ComputeTriangleGeomVolumesFilter · TriangleUtilities.cpp:922 — FindTriangleVolume builds the matrix relative to the origin in float32 and volumes are accumulated in float32, so a small closed Feature at coordinates around 1e5-1e6 loses its volume to cancellation. ComputeTriangleGeomVolumesFilter#6
• WRONG_DATA · CreateColorMapFilter · ColorTableUtilities.hpp:151 — The non-integer path divides in float32, so a float64 range larger than FLT_MAX (for example [0, 1e30, 1e40]) becomes inf, then 0 or NaN, and every value gets the first color. CreateColorMapFilter#4
• WRONG_DATA · RobustAutomaticThresholdFilter · RobustAutomaticThreshold.cpp:47 — numerator and denominator are float and summed per tuple (lines 69-70), so for 20M tuples the denominator stops at 2^24 and the threshold at line 75 is wrong. RobustAutomaticThresholdFilter#1
• MINOR · TriangleCentroidFilter · TriangleCentroid.cpp:40 — Three float32 coordinates are added and divided by 3.0F before the value goes into the float64 output, so the result loses about 1 ulp of float32 precision (for example X = 16777216, 1, 1). TriangleCentroidFilter#1

Already fixed or in progress

• PR BUG: Fix defects in seven SimplnxCore filters #1763 (merged, commit 8bed04c): ComputeKMeansFilter#7, SliceTriangleGeometryFilter#3.
• PR BUG: Correct surface-area and sphericity calculations #1817 (merged, commit 01809c2): ComputeSurfaceAreaToVolumeFilter#4.

Out of scope

• Findings in filters that have a V&V report are excluded from this issue: 7 findings, in ComputeAvgCAxesFilter, ComputeAvgOrientationsFilter, ComputeFeatureReferenceCAxisMisorientationsFilter, ComputeFeatureReferenceMisorientationsFilter, ComputeTriangleGeomCentroidsFilter and DBSCANFilter. A separate V&V branch handles them.
• Changes to exemplar archives only to make a test pass. See Rule.