BUG: Fix periodic feature-centroid model in Triangle and Feature centroid filters

src/Plugins/SimplnxCore/docs/ComputeFeatureCentroidsFilter.md

@@ -6,7 +6,7 @@ Generic (Morphological)
 
 ## Description
 
-This **Filter** calculates the *centroid* of each **Feature** by determining the average X, Y, and Z position (in physical coordinates) of all the **Cells** belonging to the **Feature**. The per-cell coordinates are accumulated using Kahan compensated summation to limit floating-point round-off on features with large cell counts. An *Is Periodic* option is available: when enabled, a **Feature** that spans the full extent of an axis (touching both opposing faces of the **Image Geometry**) is treated as wrapping around to the opposite face, and its centroid on that axis is shifted by half the physical distance between the first and last cell centers. When *Is Periodic* is disabled, **Features** that intersect the outer surfaces of the sample will still have *centroids* calculated, but they will be *centroids* of the truncated part of the **Feature** that lies inside the sample.
+This **Filter** calculates the *centroid* of each **Feature** by determining the average X, Y, and Z position (in physical coordinates) of all the **Cells** belonging to the **Feature**. The per-cell coordinates are accumulated using Kahan compensated summation to limit floating-point round-off on features with large cell counts. An *Is Periodic* option is available: when enabled, a **Feature** that spans the full extent of an axis (touching both opposing faces of the **Image Geometry**) is treated as wrapping around to the opposite face, and its centroid on that axis is computed as a minimum-image (circular) mean — the cell positions are unwrapped across the boundary before averaging and the result is mapped back into the domain. This depends on where the **Feature**'s cells actually lie, so it is correct for features whose mass wraps asymmetrically and never places the centroid outside the sample; a feature that fills the whole domain along an axis falls back to the ordinary average on that axis. When *Is Periodic* is disabled, **Features** that intersect the outer surfaces of the sample will still have *centroids* calculated, but they will be *centroids* of the truncated part of the **Feature** that lies inside the sample.
 
 A **Feature** with no **Cells** (an unused Feature Id) keeps a centroid of (0, 0, 0).
 

src/Plugins/SimplnxCore/docs/ComputeTriangleGeomCentroidsFilter.md

@@ -14,6 +14,20 @@ using the following algorithm:
    node to be entered once for a given owner*)
 3. For each **Feature**, find the average (x,y,z) coordinate from the set of nodes that bound it
 
+## Is Periodic
+
+When **Is Periodic** is enabled, features whose nodes wrap across opposing boundaries of the mesh are handled with a
+minimum-image (periodic) centroid rather than the plain average. For any axis on which a feature touches both opposing
+faces of the geometry's bounding box, the centroid component on that axis is computed by unwrapping the feature's nodes
+across the boundary (using the largest empty gap in their distribution) before averaging, and mapping the result back
+into the domain. This depends on where the feature's mass actually lies, so it is correct for asymmetrically-wrapped
+features and never places the centroid outside the domain. Axes on which the feature does not wrap keep the ordinary
+average, and a feature that fills the whole domain along an axis falls back to the ordinary average on that axis. When
+**Is Periodic** is disabled (the default), the plain average from step 3 is reported.
+
+> **Note:** Enabling **Is Periodic** assumes the periodic domain equals the geometry's bounding box. An informational
+> message is emitted for each feature that is adjusted, since a manual review may be warranted.
+
 % Auto generated parameter table will be inserted here
 
 ## Example Pipelines

src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeFeatureCentroids.cpp

@@ -4,10 +4,21 @@
 #include "simplnx/DataStructure/DataGroup.hpp"
 #include "simplnx/DataStructure/Geometry/ImageGeom.hpp"
 #include "simplnx/Utilities/DataArrayUtilities.hpp"
-#include "simplnx/Utilities/GeometryHelpers.hpp"
 #include "simplnx/Utilities/ParallelDataAlgorithm.hpp"
 
 #include <algorithm>
+#include <array>
+#include <cmath>
+
+namespace
+{
+// 2*pi, used to map cell-center positions onto the periodic unit circle.
+constexpr double k_TwoPi = 6.283185307179586;
+// A feature whose per-axis unit-vector resultant falls below this length has its mass spread ~uniformly
+// around the domain (a domain-filling feature); its circular mean is not meaningful, so the arithmetic
+// mean is kept instead.
+constexpr double k_DegenerateResultant = 1.0e-6;
+} // namespace
 
 using namespace nx::core;
 
@@ -21,7 +32,8 @@ class ComputeFeatureCentroidsImpl1
 {
 public:
   ComputeFeatureCentroidsImpl1(Float64AbstractDataStore& sum, Float64AbstractDataStore& compensation, UInt64AbstractDataStore& count, std::array<size_t, 3> dims, const nx::core::ImageGeom& imageGeom,
-                               const Int32AbstractDataStore& featureIds, UInt64AbstractDataStore& rangeXStoreRef, UInt64AbstractDataStore& rangeYStoreRef, UInt64AbstractDataStore& rangeZStoreRef)
+                               const Int32AbstractDataStore& featureIds, UInt64AbstractDataStore& rangeXStoreRef, UInt64AbstractDataStore& rangeYStoreRef, UInt64AbstractDataStore& rangeZStoreRef,
+                               bool isPeriodic, Float64AbstractDataStore& sumCos, Float64AbstractDataStore& sumSin, std::array<double, 3> origin, std::array<double, 3> domainLength)
   : m_Sum(sum)
   , m_Compensation(compensation)
   , m_Count(count)
@@ -31,6 +43,11 @@ class ComputeFeatureCentroidsImpl1
   , m_RangeXStoreRef(rangeXStoreRef)
   , m_RangeYStoreRef(rangeYStoreRef)
   , m_RangeZStoreRef(rangeZStoreRef)
+  , m_IsPeriodic(isPeriodic)
+  , m_SumCos(sumCos)
+  , m_SumSin(sumSin)
+  , m_Origin(origin)
+  , m_DomainLength(domainLength)
   {
   }
   ~ComputeFeatureCentroidsImpl1() = default;
@@ -85,6 +102,21 @@ class ComputeFeatureCentroidsImpl1
           m_Compensation[featureId_idx] = (temp - m_Sum[featureId_idx]) - componentValue;
           m_Sum[featureId_idx] = temp;
           m_Count[featureId_idx].inc();
+
+          // For periodic runs, also accumulate the unit vector of each cell center's angular position
+          // around the domain on each axis. The per-feature (cos, sin) sums produce a minimum-image
+          // centroid at finalize for features that wrap the boundary.
+          if(m_IsPeriodic)
+          {
+            const nx::core::Point3Dd center = {voxel_center[0], voxel_center[1], voxel_center[2]};
+            for(size_t axis = 0; axis < 3; axis++)
+            {
+              const double phase = k_TwoPi * (center[axis] - m_Origin[axis]) / m_DomainLength[axis];
+              const size_t axisIdx = featureId * 3ULL + axis;
+              m_SumCos[axisIdx] = m_SumCos.getValue(axisIdx) + std::cos(phase);
+              m_SumSin[axisIdx] = m_SumSin.getValue(axisIdx) + std::sin(phase);
+            }
+          }
         }
       }
     }
@@ -105,6 +137,11 @@ class ComputeFeatureCentroidsImpl1
   UInt64AbstractDataStore& m_RangeXStoreRef;
   UInt64AbstractDataStore& m_RangeYStoreRef;
   UInt64AbstractDataStore& m_RangeZStoreRef;
+  bool m_IsPeriodic = false;
+  Float64AbstractDataStore& m_SumCos;
+  Float64AbstractDataStore& m_SumSin;
+  std::array<double, 3> m_Origin = {0.0, 0.0, 0.0};
+  std::array<double, 3> m_DomainLength = {0.0, 0.0, 0.0};
 };
 
 } // namespace
@@ -169,6 +206,21 @@ Result<> ComputeFeatureCentroids::operator()()
   compensation.fill(0.0);
   count.fill(0.0);
 
+  // Per-feature/per-axis unit-vector sums for the periodic (circular-mean) centroid. Only populated when
+  // Is Periodic is enabled; harmless (all zero) otherwise.
+  auto sumCosPtr = DataStoreUtilities::CreateDataStore<float64>(tupleShape, componentShape, IDataAction::Mode::Execute);
+  auto sumSinPtr = DataStoreUtilities::CreateDataStore<float64>(tupleShape, componentShape, IDataAction::Mode::Execute);
+  Float64AbstractDataStore& sumCos = *sumCosPtr.get();
+  Float64AbstractDataStore& sumSin = *sumSinPtr.get();
+  sumCos.fill(0.0);
+  sumSin.fill(0.0);
+
+  const auto geomOrigin = imageGeom.getOrigin();
+  const auto geomSpacing = imageGeom.getSpacing();
+  const std::array<double, 3> origin = {static_cast<double>(geomOrigin[0]), static_cast<double>(geomOrigin[1]), static_cast<double>(geomOrigin[2])};
+  // Periodic domain length on each axis = number of cells * spacing (the full physical extent).
+  const std::array<double, 3> domainLength = {static_cast<double>(xPoints) * geomSpacing[0], static_cast<double>(yPoints) * geomSpacing[1], static_cast<double>(zPoints) * geomSpacing[2]};
+
   // Create data stores to check if feature IDs are periodic
   componentShape[0] = 2;
   auto rangeXStorePtr = DataStoreUtilities::CreateDataStore<uint64>(tupleShape, componentShape, IDataAction::Mode::Execute);
@@ -187,7 +239,8 @@ Result<> ComputeFeatureCentroids::operator()()
   // by the total number of cores/threads and do a ParallelTask Algorithm instead
   // we might see some speedup.
   dataAlg.setParallelizationEnabled(false);
-  dataAlg.execute(ComputeFeatureCentroidsImpl1(sum, compensation, count, {xPoints, yPoints, zPoints}, imageGeom, featureIdsStoreRef, rangeXStoreRef, rangeYStoreRef, rangeZStoreRef));
+  dataAlg.execute(ComputeFeatureCentroidsImpl1(sum, compensation, count, {xPoints, yPoints, zPoints}, imageGeom, featureIdsStoreRef, rangeXStoreRef, rangeYStoreRef, rangeZStoreRef,
+                                               m_InputValues->IsPeriodic, sumCos, sumSin, origin, domainLength));
 
   // Here we are only looping over the number of features so let this just go in serial mode.
   // The count store carries the same voxel count in all three components of a feature; a feature with
@@ -216,9 +269,48 @@ Result<> ComputeFeatureCentroids::operator()()
   if(m_InputValues->IsPeriodic)
   {
     m_MessageHandler({IFilter::Message::Type::Info, "Checking for periodic data."});
-    if(GeometryHelpers::Topology::AdjustCentroidsForPeriodicFaces(imageGeom, rangeXStoreRef, rangeYStoreRef, rangeZStoreRef, centroids))
+
+    // For each axis on which a feature spans the full extent (a cell at index 0 and at the last index), the
+    // naive arithmetic centroid lands in the empty middle of the wrapped feature. Replace that component with
+    // the circular (minimum-image) mean derived from the per-feature unit-vector sums. A feature whose mass
+    // is spread ~uniformly around the domain (near-zero resultant) keeps its arithmetic mean.
+    const std::array<UInt64AbstractDataStore*, 3> rangeStores = {&rangeXStoreRef, &rangeYStoreRef, &rangeZStoreRef};
+    const std::array<size_t, 3> dims = {xPoints, yPoints, zPoints};
+    bool anyAdjusted = false;
+    for(size_t featureId = 0; featureId < totalFeatures; featureId++)
+    {
+      for(size_t axis = 0; axis < 3; axis++)
+      {
+        const size_t axisIdx = featureId * 3 + axis;
+        if(count[axisIdx] == 0)
+        {
+          continue;
+        }
+        const UInt64AbstractDataStore& rangeStore = *rangeStores[axis];
+        const bool spansExtent = (rangeStore.getValue(featureId * 2 + 0) == 0 && rangeStore.getValue(featureId * 2 + 1) == dims[axis] - 1);
+        if(!spansExtent)
+        {
+          continue;
+        }
+        const double sumCosValue = sumCos.getValue(axisIdx);
+        const double sumSinValue = sumSin.getValue(axisIdx);
+        const double resultant = std::sqrt(sumCosValue * sumCosValue + sumSinValue * sumSinValue) / static_cast<double>(count[axisIdx]);
+        if(resultant < k_DegenerateResultant)
+        {
+          continue;
+        }
+        double phase = std::atan2(sumSinValue, sumCosValue);
+        if(phase < 0.0)
+        {
+          phase += k_TwoPi;
+        }
+        centroids[axisIdx] = static_cast<float>(origin[axis] + (phase / k_TwoPi) * domainLength[axis]);
+        anyAdjusted = true;
+      }
+    }
+    if(anyAdjusted)
     {
-      m_MessageHandler({IFilter::Message::Type::Info, "ComputeFeatureCentroids found Non-Contiguous Features. Centroids may require additional checks."});
+      m_MessageHandler({IFilter::Message::Type::Info, "ComputeFeatureCentroids adjusted centroids of features that wrap the periodic boundary."});
     }
   }
 

src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeTriangleGeomCentroids.cpp

@@ -58,6 +58,10 @@ Result<> ComputeTriangleGeomCentroids::operator()()
 
   for(MeshIndexType i = 0; i < numTriangles; i++)
   {
+    if(m_ShouldCancel)
+    {
+      return {};
+    }
     const int32 faceLabel0 = faceLabels[2 * i + 0];
     const int32 faceLabel1 = faceLabels[2 * i + 1];
     if(faceLabel0 > 0)
@@ -78,7 +82,11 @@ Result<> ComputeTriangleGeomCentroids::operator()()
 
   for(MeshIndexType i = 0; i < numFeatures; i++)
   {
-    std::set<MeshIndexType> vertexSet = vertexSets[i];
+    if(m_ShouldCancel)
+    {
+      return {};
+    }
+    const std::set<MeshIndexType>& vertexSet = vertexSets[i];
     auto periodicFaces = GeometryHelpers::Topology::FindElementPeriodicFaces(boundingBox, vertexCoords, vertexSet);
 
     for(const auto& vert : vertexSets[i])
@@ -95,10 +103,11 @@ Result<> ComputeTriangleGeomCentroids::operator()()
 
       if(m_InputValues->IsPeriodic)
       {
-        if(GeometryHelpers::Topology::AdjustCentroidsForPeriodicFaces(boundingBox, periodicFaces, centroids, i))
+        if(GeometryHelpers::Topology::AdjustCentroidsForPeriodicFaces(boundingBox, periodicFaces, vertexCoords, vertexSets[i], centroids, i))
         {
-          IFilter::Message warningMsg{IFilter::Message::Type::Info, fmt::format("Feature ID {} may be periodic. Manual review may be necessary.", i)};
-          m_MessageHandler.m_Callback(warningMsg);
+          // Use the MessageHandler's call operator (which guards against an empty callback) rather than
+          // invoking m_Callback directly, which throws std::bad_function_call when no handler is installed.
+          m_MessageHandler(IFilter::Message{IFilter::Message::Type::Info, fmt::format("Feature ID {} may be periodic. Manual review may be necessary.", i)});
         }
       }
     }

src/Plugins/SimplnxCore/test/ComputeFeatureCentroidsTest.cpp

@@ -160,37 +160,49 @@ TEST_CASE("SimplnxCore::ComputeFeatureCentroidsFilter: Class 1/4 - Periodic Boun
   using namespace CentroidToy;
   UnitTest::LoadPlugins();
 
-  SECTION("Fixture D - periodic wrap, unit spacing")
+  SECTION("Fixture D - periodic wrap, unit spacing (asymmetric)")
   {
     // 4x1x1; FeatureIds [1,2,1,1]; feature 1 spans x=0 and x=3 (full extent), feature 2 (x=1) does not.
-    // feature 1 cells x=0,2,3 -> centers 0.5,2.5,3.5 -> mean 6.5/3 = 2.16667
+    // feature 1 cells x=0,2,3 -> centers 0.5,2.5,3.5 -> non-periodic mean 6.5/3 = 2.16667
     auto sNP = Build(4, 1, 1, {1.0f, 1.0f, 1.0f}, {0.0f, 0.0f, 0.0f}, 3, {1, 2, 1, 1});
     auto nonPeriodic = Run(sNP, false);
     RequireCentroid(nonPeriodic, 1, 6.5f / 3.0f, 0.5f, 0.5f);
     RequireCentroid(nonPeriodic, 2, 1.5f, 0.5f, 0.5f); // single cell x=1
 
     auto sP = Build(4, 1, 1, {1.0f, 1.0f, 1.0f}, {0.0f, 0.0f, 0.0f}, 3, {1, 2, 1, 1});
     auto periodic = Run(sP, true);
-    // Class 1 (given the DREAM3D periodic-shift model): spanning feature 1 gets +(dim-1)/2 = +1.5.
-    RequireCentroid(periodic, 1, 6.5f / 3.0f + 1.5f, 0.5f, 0.5f);
+    // Class 1 minimum-image oracle: domain L=4 (seam at x=0). The cell at center 0.5 wraps to 4.5, so the
+    // contiguous cluster is {2.5, 3.5, 4.5}; mean = 10.5/3 = 3.5, mapped back into [0,4) -> 3.5.
+    // (The old constant-offset model gave 6.5/3 + 1.5 = 3.6667, which is wrong.)
+    RequireCentroid(periodic, 1, 3.5f, 0.5f, 0.5f);
     // Class 4 invariant: the non-spanning feature 2 must be unchanged by the periodic pass.
     RequireCentroid(periodic, 2, 1.5f, 0.5f, 0.5f);
   }
 
-  SECTION("Fixture E - periodic wrap, NON-unit spacing (regression pin for the spacing fix)")
+  SECTION("Fixture E - periodic wrap, NON-unit spacing (asymmetric, interior result)")
   {
-    // 4x1x1, spacing (2,1,1), origin (10,0,0); FeatureIds [1,2,2,1]
-    // x-centers: idx0 -> 11.0, idx3 -> 17.0; non-periodic centroid.x = 14.0
-    auto sNP = Build(4, 1, 1, {2.0f, 1.0f, 1.0f}, {10.0f, 0.0f, 0.0f}, 3, {1, 2, 2, 1});
+    // 4x1x1, spacing (2,1,1), origin (10,0,0); FeatureIds [1,1,2,1]
+    // feature 1 cells idx 0,1,3 -> x-centers 11.0, 13.0, 17.0; non-periodic centroid.x = 41/3 = 13.6667
+    auto sNP = Build(4, 1, 1, {2.0f, 1.0f, 1.0f}, {10.0f, 0.0f, 0.0f}, 3, {1, 1, 2, 1});
     auto nonPeriodic = Run(sNP, false);
-    RequireCentroid(nonPeriodic, 1, 14.0f, 0.5f, 0.5f); // unambiguous analytical value
+    RequireCentroid(nonPeriodic, 1, 41.0f / 3.0f, 0.5f, 0.5f);
 
-    auto sP = Build(4, 1, 1, {2.0f, 1.0f, 1.0f}, {10.0f, 0.0f, 0.0f}, 3, {1, 2, 2, 1});
+    auto sP = Build(4, 1, 1, {2.0f, 1.0f, 1.0f}, {10.0f, 0.0f, 0.0f}, 3, {1, 1, 2, 1});
     auto periodic = Run(sP, true);
-    // The periodic offset must scale with spacing: (dim-1)*spacing_x/2 = 3*2/2 = 3.0 -> 14.0 + 3.0 = 17.0.
-    // Before the fix this returned 15.5 (offset 1.5 in cell units, ignoring spacing). Regression pin for
-    // the AdjustCentroidsForPeriodicFaces spacing fix. See vv/ComputeFeatureCentroidsFilter.md Phase 6.
-    RequireCentroid(periodic, 1, 17.0f, 0.5f, 0.5f);
+    // Class 1 minimum-image oracle: domain L=8 (seam at x=10 == x=18). Positions relative to origin are
+    // {1, 3, 7}; the largest empty gap is 3->7, so 1 and 3 wrap to 9 and 11. Cluster {7, 9, 11}; mean = 9,
+    // mapped back into [10,18) -> x = 11.0. (The old constant-offset model gave 13.6667 + 3.0 = 16.6667.)
+    RequireCentroid(periodic, 1, 11.0f, 0.5f, 0.5f);
+  }
+
+  SECTION("Fixture F - domain-filling feature (arithmetic-mean fallback)")
+  {
+    // 4x1x1, unit spacing; a single feature fills every cell (centers 0.5,1.5,2.5,3.5). It spans the full
+    // extent, so the periodic branch fires, but the unit vectors cancel (zero resultant): the circular mean
+    // is undefined, so the arithmetic mean 2.0 is kept. Class 4 fallback path.
+    auto sP = Build(4, 1, 1, {1.0f, 1.0f, 1.0f}, {0.0f, 0.0f, 0.0f}, 2, {1, 1, 1, 1});
+    auto periodic = Run(sP, true);
+    RequireCentroid(periodic, 1, 2.0f, 0.5f, 0.5f);
   }
 }