BUG: New implementation for ComputeNeighborhoods filter

The previous implementation has issues with the binning system and will only find
that a centroid is within a radius if both centroids fall into the same bin.

New unit test file was created and hand verified by MAJ.

Author: imikejackson

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

@@ -1,45 +1,108 @@
 #include "ComputeNeighborhoods.hpp"
 
+#include "simplnx/Common/Array.hpp"
+#include "simplnx/Common/Range.hpp"
+#include "simplnx/Common/Result.hpp"
+#include "simplnx/Common/Types.hpp"
+#include "simplnx/DataStructure/AbstractDataStore.hpp"
 #include "simplnx/DataStructure/DataArray.hpp"
-#include "simplnx/DataStructure/DataGroup.hpp"
+#include "simplnx/DataStructure/DataStructure.hpp"
 #include "simplnx/DataStructure/Geometry/ImageGeom.hpp"
+#include "simplnx/Filter/IFilter.hpp"
+#include "simplnx/Utilities/MessageHelper.hpp"
 #include "simplnx/Utilities/ParallelDataAlgorithm.hpp"
 
-#include <cmath>
+#include <fmt/format.h>
+
+#include <atomic>
+#include <functional>
+#include <mutex>
+#include <unordered_map>
+#include <vector>
 
 using namespace nx::core;
 namespace
 {
+struct BinKey
+{
+  int64 bx = 0;
+  int64 by = 0;
+  int64 bz = 0;
+
+  BinKey(int64 x, int64 y, int64 z)
+  : bx(x)
+  , by(y)
+  , bz(z)
+  {
+  }
+
+  BinKey(const std::vector<int64>& bins, usize tupleIndex)
+  {
+    bx = bins[3 * tupleIndex + 0];
+    by = bins[3 * tupleIndex + 1];
+    bz = bins[3 * tupleIndex + 2];
+  }
+
+  bool operator==(const BinKey& other) const noexcept
+  {
+    return bx == other.bx && by == other.by && bz == other.bz;
+  }
+};
+
+struct BinKeyHasher
+{
+  usize operator()(const BinKey& key) const noexcept
+  {
+    // simple hash combine
+    const usize h1 = std::hash<int64>{}(key.bx);
+    const usize h2 = std::hash<int64>{}(key.by);
+    const usize h3 = std::hash<int64>{}(key.bz);
+    usize seed = h1;
+    seed ^= h2 + 0x9e3779b97f4a7c15ULL + (seed << 6) + (seed >> 2);
+    seed ^= h3 + 0x9e3779b97f4a7c15ULL + (seed << 6) + (seed >> 2);
+    return seed;
+  }
+};
+
 class ComputeNeighborhoodsImpl
 {
 public:
-  ComputeNeighborhoodsImpl(ComputeNeighborhoods* filter, usize totalFeatures, const std::vector<int64_t>& bins, const std::vector<float>& criticalDistance, const std::atomic_bool& shouldCancel,
-                           ProgressMessageHelper& progressMessageHelper)
+  ComputeNeighborhoodsImpl(ComputeNeighborhoods* filter, const nx::core::AbstractDataStore<float>& centroids, const std::vector<int64>& bins, float32 avgDiam, float32 multiplesOfAverage,
+                           const std::atomic_bool& shouldCancel, ProgressMessageHelper& progressMessageHelper)
   : m_Filter(filter)
-  , m_TotalFeatures(totalFeatures)
+  , m_Centroids(centroids)
   , m_Bins(bins)
-  , m_CriticalDistance(criticalDistance)
+  , m_AvgDiam(avgDiam)
+  , m_MultiplesOfAverage(multiplesOfAverage)
   , m_ShouldCancel(shouldCancel)
   , m_ProgressMessageHelper(progressMessageHelper)
   {
   }
 
   void convert(usize start, usize end) const
   {
-    int64 bin1x, bin2x, bin1y, bin2y, bin1z, bin2z;
-    float32 dBinX, dBinY, dBinZ;
-    float32 criticalDistance1, criticalDistance2;
-
-    auto increment = static_cast<float64>(end - start) / 100.0;
-    float64 incCount = 0.0;
-    // NEVER start at 0.
-    if(start == 0)
+    const auto increment = static_cast<int64>((end - start) / 100.0);
+    int64 incCount = 0.0;
+
+    const usize totalFeatures = m_Centroids.getNumberOfTuples();
+
+    // 1. Build spatial grid: BinKey -> list of features
+    std::unordered_map<BinKey, std::vector<usize>, BinKeyHasher> binToFeatures;
+    binToFeatures.reserve(totalFeatures);
+
+    for(usize i = 1; i < totalFeatures; ++i) // assuming feature 0 is background
     {
-      start = 1;
+      const BinKey key(m_Bins, i);
+      binToFeatures[key].push_back(i);
     }
 
+    // 2. Precompute radius info
+    const float32 radius = m_AvgDiam * m_MultiplesOfAverage / 2.0f;
+    const float32 radiusSq = radius * radius;
+    const int64 k = static_cast<int64>(std::ceil(m_MultiplesOfAverage));
+
     ProgressMessenger progressMessenger = m_ProgressMessageHelper.createProgressMessenger();
-    for(usize featureIdx = start; featureIdx < end; featureIdx++)
+    for(usize i = start; i < end; i++)
     {
       incCount++;
       if(incCount >= increment)
@@ -52,31 +115,56 @@ class ComputeNeighborhoodsImpl
       {
         return;
       }
+      // (a) Get feature's i position
+      const float32 xi = m_Centroids[3 * i + 0];
+      const float32 yi = m_Centroids[3 * i + 1];
+      const float32 zi = m_Centroids[3 * i + 2];
 
-      bin1x = m_Bins[3 * featureIdx];
-      bin1y = m_Bins[3 * featureIdx + 1];
-      bin1z = m_Bins[3 * featureIdx + 2];
-      criticalDistance1 = m_CriticalDistance[featureIdx];
+      // (b) Get its bin
+      const int64 bx0 = m_Bins[3 * i + 0];
+      const int64 by0 = m_Bins[3 * i + 1];
+      const int64 bz0 = m_Bins[3 * i + 2];
 
-      for(usize j = featureIdx + 1; j < m_TotalFeatures; j++)
+      // (c) Scan all bins within +/- k in each dimension
+      for(int64 dbx = -k; dbx <= k; ++dbx)
       {
-        bin2x = m_Bins[3 * j];
-        bin2y = m_Bins[3 * j + 1];
-        bin2z = m_Bins[3 * j + 2];
-        criticalDistance2 = m_CriticalDistance[j];
+        for(int64 dby = -k; dby <= k; ++dby)
+        {
+          for(int64 dbz = -k; dbz <= k; ++dbz)
+          {
+            const BinKey nbKey{bx0 + dbx, by0 + dby, bz0 + dbz};
 
-        dBinX = std::abs(static_cast<float32>(bin2x - bin1x));
-        dBinY = std::abs(static_cast<float32>(bin2y - bin1y));
-        dBinZ = std::abs(static_cast<float32>(bin2z - bin1z));
+            auto it = binToFeatures.find(nbKey);
+            if(it == binToFeatures.end())
+            {
+              continue; // no features in this bin
+            }
 
-        if(dBinX < criticalDistance1 && dBinY < criticalDistance1 && dBinZ < criticalDistance1)
-        {
-          m_Filter->updateNeighborHood(featureIdx, j);
-        }
+            const std::vector<usize>& candidates = it->second;
 
-        if(dBinX < criticalDistance2 && dBinY < criticalDistance2 && dBinZ < criticalDistance2)
-        {
-          m_Filter->updateNeighborHood(j, featureIdx);
+            // (d) Check actual distances to candidates in this bin
+            for(const usize j : candidates)
+            {
+              if(j == i)
+              {
+                continue; // skip self
+              }
+
+              const float32 xj = m_Centroids[3 * j + 0];
+              const float32 yj = m_Centroids[3 * j + 1];
+              const float32 zj = m_Centroids[3 * j + 2];
+
+              const float32 dx = xi - xj;
+              const float32 dy = yi - yj;
+              const float32 dz = zi - zj;
+
+              const float32 distSq = dx * dx + dy * dy + dz * dz;
+              if(distSq <= radiusSq)
+              {
+                m_Filter->updateNeighborHood(i, j);
+              }
+            }
+          }
         }
       }
     }
@@ -90,9 +178,10 @@ class ComputeNeighborhoodsImpl
 
 private:
   ComputeNeighborhoods* m_Filter = nullptr;
-  usize m_TotalFeatures = 0;
+  const nx::core::AbstractDataStore<float>& m_Centroids;
   const std::vector<int64>& m_Bins;
-  const std::vector<float32>& m_CriticalDistance;
+  float32 m_AvgDiam;
+  float32 m_MultiplesOfAverage;
   const std::atomic_bool& m_ShouldCancel;
   ProgressMessageHelper& m_ProgressMessageHelper;
 };
@@ -111,16 +200,10 @@ ComputeNeighborhoods::ComputeNeighborhoods(DataStructure& dataStructure, const I
 // -----------------------------------------------------------------------------
 ComputeNeighborhoods::~ComputeNeighborhoods() noexcept = default;
 
-// -----------------------------------------------------------------------------
-const std::atomic_bool& ComputeNeighborhoods::getCancel()
-{
-  return m_ShouldCancel;
-}
-
 // -----------------------------------------------------------------------------
 void ComputeNeighborhoods::updateNeighborHood(usize sourceIndex, usize destIndex)
 {
-  const std::lock_guard<std::mutex> lock(m_Mutex);
+  const std::scoped_lock lock(m_Mutex);
   (*m_Neighborhoods)[sourceIndex].inc();
   m_LocalNeighborhoodList[sourceIndex].push_back(static_cast<int32_t>(destIndex));
 }
@@ -129,60 +212,59 @@ void ComputeNeighborhoods::updateNeighborHood(usize sourceIndex, usize destIndex
 Result<> ComputeNeighborhoods::operator()()
 {
   // m_ProgressCounter initialized to zero on filter creation
-  std::vector<float32> criticalDistance;
-
   auto multiplesOfAverage = m_InputValues->MultiplesOfAverage;
   const auto& equivalentDiameters = m_DataStructure.getDataAs<Float32Array>(m_InputValues->EquivalentDiametersArrayPath)->getDataStoreRef();
   const auto& centroids = m_DataStructure.getDataAs<Float32Array>(m_InputValues->CentroidsArrayPath)->getDataStoreRef();
 
   m_Neighborhoods = m_DataStructure.getDataAs<Int32Array>(m_InputValues->NeighborhoodsArrayName);
 
-  usize totalFeatures = equivalentDiameters.getNumberOfTuples();
+  const usize totalFeatures = equivalentDiameters.getNumberOfTuples();
 
   ProgressMessageHelper progressMessageHelper = m_MessageHelper.createProgressMessageHelper();
   progressMessageHelper.setMaxProgresss(totalFeatures);
   progressMessageHelper.setProgressMessageTemplate("Finding Feature Neighborhoods: {:.2f}%");
 
   m_LocalNeighborhoodList.resize(totalFeatures);
-  criticalDistance.resize(totalFeatures);
 
-  float32 aveDiam = 0.0f;
+  // (a) This section finds the average equivalent spherical (ESD) diameter of ALL features
+  float32 avgDiameter = 0.0f;
   for(usize i = 1; i < totalFeatures; i++)
   {
     (*m_Neighborhoods)[i] = 0;
-    aveDiam += equivalentDiameters[i];
-    criticalDistance[i] = equivalentDiameters[i] * multiplesOfAverage;
-  }
-  aveDiam /= static_cast<float32>(totalFeatures);
-  for(usize i = 1; i < totalFeatures; i++)
-  {
-    criticalDistance[i] /= aveDiam;
+    avgDiameter += equivalentDiameters[i];
   }
+  avgDiameter /= static_cast<float32>(totalFeatures);
+  m_MessageHandler(IFilter::Message::Type::Info, fmt::format("Feature Average Diameter: '{}'", avgDiameter));
 
+  // (c) We are going to place each feature's centroid into a bin in the normalized 3D space.
+  // The centroid is normalized by the Average Diameter
   std::vector<int64> bins(3 * totalFeatures, 0);
   FloatVec3 origin = m_DataStructure.getDataAs<ImageGeom>(m_InputValues->InputImageGeometry)->getOrigin();
   for(usize i = 1; i < totalFeatures; i++)
   {
-    float32 x = centroids[3 * i];
-    float32 y = centroids[3 * i + 1];
-    float32 z = centroids[3 * i + 2];
-    bins[3 * i] = static_cast<int64>((x - origin[0]) / aveDiam);     // x-Bin
-    bins[3 * i + 1] = static_cast<int64>((y - origin[1]) / aveDiam); // y-Bin
-    bins[3 * i + 2] = static_cast<int64>((z - origin[2]) / aveDiam); // z-Bin
+    const float32 x = centroids[3 * i];
+    const float32 y = centroids[3 * i + 1];
+    const float32 z = centroids[3 * i + 2];
+    bins[3 * i] = static_cast<int64>((x - origin[0]) / avgDiameter);     // x-Bin
+    bins[3 * i + 1] = static_cast<int64>((y - origin[1]) / avgDiameter); // y-Bin
+    bins[3 * i + 2] = static_cast<int64>((z - origin[2]) / avgDiameter); // z-Bin
+  }
+  if(m_ShouldCancel)
+  {
+    return {};
   }
-
   ParallelDataAlgorithm parallelAlgorithm;
   parallelAlgorithm.setRange(Range(0, totalFeatures));
   parallelAlgorithm.setParallelizationEnabled(true);
-  parallelAlgorithm.execute(ComputeNeighborhoodsImpl(this, totalFeatures, bins, criticalDistance, m_ShouldCancel, progressMessageHelper));
+  parallelAlgorithm.execute(ComputeNeighborhoodsImpl(this, centroids, bins, avgDiameter, multiplesOfAverage, m_ShouldCancel, progressMessageHelper));
 
   // Output Variables
   auto& outputNeighborList = m_DataStructure.getDataRefAs<NeighborList<int32>>(m_InputValues->NeighborhoodListArrayName);
   // Set the vector for each list into the NeighborList Object
   for(usize i = 1; i < totalFeatures; i++)
   {
     // Construct a shared vector<int32> through the std::vector<> copy constructor.
-    NeighborList<int32>::SharedVectorType sharedMisOrientationList(new std::vector<int32>(m_LocalNeighborhoodList[i]));
+    const NeighborList<int32>::SharedVectorType sharedMisOrientationList(new std::vector<int32>(m_LocalNeighborhoodList[i]));
     outputNeighborList.setList(static_cast<int32>(i), sharedMisOrientationList);
   }
 

src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/ComputeNeighborhoods.hpp

@@ -41,8 +41,6 @@ class SIMPLNXCORE_EXPORT ComputeNeighborhoods
 
   Result<> operator()();
 
-  const std::atomic_bool& getCancel();
-
   void updateNeighborHood(usize sourceIndex, usize targetIndex);
 
 private:

src/Plugins/SimplnxCore/test/CMakeLists.txt

@@ -273,6 +273,7 @@ if(EXISTS "${DREAM3D_DATA_DIR}" AND SIMPLNX_DOWNLOAD_TEST_FILES)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME ReadZeissTxmFileTest_v2.tar.gz SHA512 2d962655d919480fa7baf0c8720affbcbd1182e8482480f5413b53ead41787d35d0f15487219163ba0dcffad69be90e627730fc780107d80df3d97e14fb15cac)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME segment_features_neighbor_scheme_test.tar.gz SHA512 4700b42eb3a9058dc08520c786d0e6aa20e45bfa06233952462c0c8a428e437406f530125822f1742dfe57405be66e4f0fd5ac2f76db4a13379a7ae189d370cd)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME write_lammps_test.tar.gz SHA512 82bb5360b76e857f3233d37733c602f67fd2ac667e49b24741a70ab649e8046fb7905493df37d142808b740c2771fe7cdccd71c9d70679afafe398529ee5771e)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME compute_feature_neighborhoods.tar.gz SHA512 dda96546de1b924d8145e8d173d00605bbf7d311c8d3719edd9d698bb8a326009fcfdb885d2fe08204b8e8dcbebe2eb699d3255e25fe34cfc36b9f950a8e07da)
 
 endif()