ENH: Add slice-by-slice processing option to Isolate Largest Feature. (#1169)

Signed-off-by: Joey Kleingers <[email protected]>

Author: joeykleingers

src/Plugins/SimplnxCore/docs/IdentifySampleFilter.md

@@ -17,6 +17,8 @@ If *Fill Holes* is set to *true*:
 
 *Note:* if there are in fact "holes" in the sample, then this **Filter** will "close" them (if *Fill Holes* is set to true) by calling all the **Cells** "inside" the sample *good*.  If the user wants to reidentify those holes, then reuse the threshold **Filter** with the criteria of *GoodVoxels = 1* and whatever original criteria identified the "holes", as this will limit applying those original criteria to within the sample and not the outer border region.
 
+*Additional Note:* Only completely water-tight, internal holes within the sample are addressed when *Fill Holes* is enabled.  To fill in a contiguous group of good cells that includes holes located along the outer edge of the sample, try enabling *Process Data Slice-By-Slice*.  For each slice of the chosen plane, this will search for the largest contiguous set of *good* **Cells**, set all other *good* **Cells** to be *bad* **Cells**, and (if *Fill Holes* is enabled) fill all water-tight holes PER SLICE instead of the whole 3D volume at once.  This option can be used to allow non water-tight holes to be filled without also accidentally filling the surrounding overscan area.
+
 | Name | Description |
 |------|-------------|
 |![Small IN100 IPF Map](Images/Small_IN100.png) | Good dataset to use this filter |

src/Plugins/SimplnxCore/src/SimplnxCore/Filters/IdentifySampleFilter.cpp

@@ -4,6 +4,7 @@
 #include "simplnx/DataStructure/Geometry/ImageGeom.hpp"
 #include "simplnx/Parameters/ArraySelectionParameter.hpp"
 #include "simplnx/Parameters/BoolParameter.hpp"
+#include "simplnx/Parameters/ChoicesParameter.hpp"
 #include "simplnx/Parameters/DataGroupSelectionParameter.hpp"
 #include "simplnx/Parameters/GeometrySelectionParameter.hpp"
 #include "simplnx/Utilities/FilterUtilities.hpp"
@@ -46,11 +47,8 @@ struct IdentifySampleFunctor
     std::vector<bool> checked(totalPoints, false);
     std::vector<bool> sample(totalPoints, false);
     int64 biggestBlock = 0;
-    usize count = 0;
-    int32 good = 0;
+
     int64 neighbor = 0;
-    int64 column = 0, row = 0, plane = 0;
-    int64 index = 0;
 
     // In this loop over the data we are finding the biggest contiguous set of GoodVoxels and calling that the 'sample'  All GoodVoxels that do not touch the 'sample'
     // are flipped to be called 'bad' voxels or 'not sample'
@@ -70,16 +68,16 @@ struct IdentifySampleFunctor
       if(!checked[i] && goodVoxels.getValue(i))
       {
         currentVList.push_back(i);
-        count = 0;
+        usize count = 0;
         while(count < currentVList.size())
         {
-          index = currentVList[count];
-          column = index % xp;
-          row = (index / xp) % yp;
-          plane = index / (xp * yp);
+          int64 index = currentVList[count];
+          int64 column = index % xp;
+          int64 row = (index / xp) % yp;
+          int64 plane = index / (xp * yp);
           for(int32 j = 0; j < 6; j++)
           {
-            good = 1;
+            int32 good = 1;
             neighbor = index + neighborPoints[j];
             if(j == 0 && plane == 0)
             {
@@ -156,21 +154,21 @@ struct IdentifySampleFunctor
         if(!checked[i] && !goodVoxels.getValue(i))
         {
           currentVList.push_back(i);
-          count = 0;
+          usize count = 0;
           touchesBoundary = false;
           while(count < currentVList.size())
           {
-            index = currentVList[count];
-            column = index % xp;
-            row = (index / xp) % yp;
-            plane = index / (xp * yp);
+            int64 index = currentVList[count];
+            int64 column = index % xp;
+            int64 row = (index / xp) % yp;
+            int64 plane = index / (xp * yp);
             if(column == 0 || column == (xp - 1) || row == 0 || row == (yp - 1) || plane == 0 || plane == (zp - 1))
             {
               touchesBoundary = true;
             }
             for(int32 j = 0; j < 6; j++)
             {
-              good = 1;
+              int32 good = 1;
               neighbor = index + neighborPoints[j];
               if(j == 0 && plane == 0)
               {
@@ -218,6 +216,200 @@ struct IdentifySampleFunctor
     checked.clear();
   }
 };
+
+struct IdentifySampleSliceBySliceFunctor
+{
+  enum class Plane
+  {
+    XY,
+    XZ,
+    YZ
+  };
+
+  template <typename T>
+  void operator()(const ImageGeom* imageGeom, IDataArray* goodVoxelsPtr, bool fillHoles, Plane plane)
+  {
+    auto& goodVoxels = goodVoxelsPtr->template getIDataStoreRefAs<AbstractDataStore<T>>();
+
+    SizeVec3 uDims = imageGeom->getDimensions();
+    const int64 dimX = static_cast<int64>(uDims[0]);
+    const int64 dimY = static_cast<int64>(uDims[1]);
+    const int64 dimZ = static_cast<int64>(uDims[2]);
+
+    int64 planeDim1, planeDim2, fixedDim;
+    int64 stride1, stride2, fixedStride;
+
+    switch(plane)
+    {
+    case Plane::XY:
+      planeDim1 = dimX;
+      planeDim2 = dimY;
+      fixedDim = dimZ;
+      stride1 = 1;
+      stride2 = dimX;
+      fixedStride = dimX * dimY;
+      break;
+
+    case Plane::XZ:
+      planeDim1 = dimX;
+      planeDim2 = dimZ;
+      fixedDim = dimY;
+      stride1 = 1;
+      stride2 = dimX * dimY;
+      fixedStride = dimX;
+      break;
+
+    case Plane::YZ:
+      planeDim1 = dimY;
+      planeDim2 = dimZ;
+      fixedDim = dimX;
+      stride1 = dimX;
+      stride2 = dimX * dimY;
+      fixedStride = 1;
+      break;
+    }
+
+    for(int64 fixedIdx = 0; fixedIdx < fixedDim; ++fixedIdx) // Process each slice
+    {
+      std::vector<bool> checked(planeDim1 * planeDim2, false);
+      std::vector<bool> sample(planeDim1 * planeDim2, false);
+      std::vector<int64> currentVList;
+      int64 biggestBlock = 0;
+
+      // Identify the largest contiguous set of good voxels in the slice
+      for(int64 p2 = 0; p2 < planeDim2; ++p2)
+      {
+        for(int64 p1 = 0; p1 < planeDim1; ++p1)
+        {
+          int64 planeIndex = p2 * planeDim1 + p1;
+          int64 globalIndex = fixedIdx * fixedStride + p2 * stride2 + p1 * stride1;
+
+          if(!checked[planeIndex] && goodVoxels.getValue(globalIndex))
+          {
+            currentVList.push_back(planeIndex);
+            int64 count = 0;
+
+            while(count < currentVList.size())
+            {
+              int64 localIdx = currentVList[count];
+              int64 localP1 = localIdx % planeDim1;
+              int64 localP2 = localIdx / planeDim1;
+
+              for(int j = 0; j < 4; ++j)
+              {
+                int64 dp1[4] = {0, 0, -1, 1};
+                int64 dp2[4] = {-1, 1, 0, 0};
+
+                int64 neighborP1 = localP1 + dp1[j];
+                int64 neighborP2 = localP2 + dp2[j];
+
+                if(neighborP1 >= 0 && neighborP1 < planeDim1 && neighborP2 >= 0 && neighborP2 < planeDim2)
+                {
+                  int64 neighborIdx = neighborP2 * planeDim1 + neighborP1;
+                  int64 globalNeighborIdx = fixedIdx * fixedStride + neighborP2 * stride2 + neighborP1 * stride1;
+
+                  if(!checked[neighborIdx] && goodVoxels.getValue(globalNeighborIdx))
+                  {
+                    currentVList.push_back(neighborIdx);
+                    checked[neighborIdx] = true;
+                  }
+                }
+              }
+              count++;
+            }
+
+            if(static_cast<int64>(currentVList.size()) > biggestBlock)
+            {
+              biggestBlock = currentVList.size();
+              sample.assign(planeDim1 * planeDim2, false);
+              for(int64 idx : currentVList)
+              {
+                sample[idx] = true;
+              }
+            }
+            currentVList.clear();
+          }
+        }
+      }
+
+      for(int64 p2 = 0; p2 < planeDim2; ++p2)
+      {
+        for(int64 p1 = 0; p1 < planeDim1; ++p1)
+        {
+          int64 planeIndex = p2 * planeDim1 + p1;
+          int64 globalIndex = fixedIdx * fixedStride + p2 * stride2 + p1 * stride1;
+
+          if(!sample[planeIndex])
+          {
+            goodVoxels.setValue(globalIndex, false);
+          }
+        }
+      }
+
+      if(fillHoles)
+      {
+        for(int64 p2 = 0; p2 < planeDim2; ++p2)
+        {
+          for(int64 p1 = 0; p1 < planeDim1; ++p1)
+          {
+            int64 planeIndex = p2 * planeDim1 + p1;
+            int64 globalIndex = fixedIdx * fixedStride + p2 * stride2 + p1 * stride1;
+
+            if(!checked[planeIndex] && !goodVoxels.getValue(globalIndex))
+            {
+              currentVList.push_back(planeIndex);
+              int64 count = 0;
+              bool touchesBoundary = false;
+
+              while(count < currentVList.size())
+              {
+                int64 localIdx = currentVList[count];
+                int64 localP1 = localIdx % planeDim1;
+                int64 localP2 = localIdx / planeDim1;
+
+                if(localP1 == 0 || localP1 == planeDim1 - 1 || localP2 == 0 || localP2 == planeDim2 - 1)
+                {
+                  touchesBoundary = true;
+                }
+
+                for(int j = 0; j < 4; ++j)
+                {
+                  int64 dp1[4] = {0, 0, -1, 1};
+                  int64 dp2[4] = {-1, 1, 0, 0};
+
+                  int64 neighborP1 = localP1 + dp1[j];
+                  int64 neighborP2 = localP2 + dp2[j];
+
+                  if(neighborP1 >= 0 && neighborP1 < planeDim1 && neighborP2 >= 0 && neighborP2 < planeDim2)
+                  {
+                    int64 neighborIdx = neighborP2 * planeDim1 + neighborP1;
+                    int64 globalNeighborIdx = fixedIdx * fixedStride + neighborP2 * stride2 + neighborP1 * stride1;
+
+                    if(!checked[neighborIdx] && !goodVoxels.getValue(globalNeighborIdx))
+                    {
+                      currentVList.push_back(neighborIdx);
+                      checked[neighborIdx] = true;
+                    }
+                  }
+                }
+                count++;
+              }
+
+              if(!touchesBoundary)
+              {
+                for(int64 idx : currentVList)
+                {
+                  goodVoxels.setValue(fixedIdx * fixedStride + idx, true);
+                }
+              }
+              currentVList.clear();
+            }
+          }
+        }
+      }
+    }
+  }
+};
 } // namespace
 
 //------------------------------------------------------------------------------
@@ -257,11 +449,26 @@ Parameters IdentifySampleFilter::parameters() const
 
   params.insertSeparator(Parameters::Separator{"Input Parameter(s)"});
   params.insert(std::make_unique<BoolParameter>(k_FillHoles_Key, "Fill Holes in Largest Feature", "Whether to fill holes within sample after it is identified", true));
+  params.insertLinkableParameter(std::make_unique<BoolParameter>(k_SliceBySlice_Key, "Process Data Slice-By-Slice",
+                                                                 "Whether to identify the largest sample (and optionally fill holes) slice-by-slice.  This option is useful if you have a sample that "
+                                                                 "is not water-tight and the holes open up to the overscan section, or if you have holes that sit on a boundary.  The original "
+                                                                 "algorithm will not fill holes that have these characteristics, only holes that are completely enclosed by the sample and "
+                                                                 "water-tight.  If you have holes that are not water-tight or sit on a boundary, choose this option and then pick the plane that will "
+                                                                 "allow the holes to be water-tight on each slice of that plane.",
+                                                                 false));
+  params.insert(
+      std::make_unique<ChoicesParameter>(k_SliceBySlicePlane_Key, "Slice-By-Slice Plane",
+                                         "Set the plane that the data will be processed slice-by-slice.  For example, if you pick the XY plane, the data will be processed in the Z direction.", 0,
+                                         ChoicesParameter::Choices{"XY", "XZ", "YZ"}));
+
   params.insert(std::make_unique<GeometrySelectionParameter>(k_SelectedImageGeometryPath_Key, "Image Geometry", "DataPath to the target ImageGeom", DataPath(),
                                                              GeometrySelectionParameter::AllowedTypes{IGeometry::Type::Image}));
   params.insert(std::make_unique<ArraySelectionParameter>(k_MaskArrayPath_Key, "Mask Array", "DataPath to the mask array defining what is sample and what is not", DataPath(),
                                                           ArraySelectionParameter::AllowedTypes{nx::core::DataType::boolean, nx::core::DataType::uint8},
                                                           ArraySelectionParameter::AllowedComponentShapes{{1}}));
+
+  params.linkParameters(k_SliceBySlice_Key, k_SliceBySlicePlane_Key, true);
+
   return params;
 }
 
@@ -301,13 +508,22 @@ Result<> IdentifySampleFilter::executeImpl(DataStructure& dataStructure, const A
                                            const std::atomic_bool& shouldCancel) const
 {
   const auto fillHoles = args.value<bool>(k_FillHoles_Key);
+  const auto sliceBySlice = args.value<bool>(k_SliceBySlice_Key);
+  const auto sliceBySlicePlane = static_cast<IdentifySampleSliceBySliceFunctor::Plane>(args.value<ChoicesParameter::ValueType>(k_SliceBySlicePlane_Key));
   const auto imageGeomPath = args.value<DataPath>(k_SelectedImageGeometryPath_Key);
   const auto goodVoxelsArrayPath = args.value<DataPath>(k_MaskArrayPath_Key);
 
   auto* inputData = dataStructure.getDataAs<IDataArray>(goodVoxelsArrayPath);
   const auto* imageGeom = dataStructure.getDataAs<ImageGeom>(imageGeomPath);
 
-  ExecuteDataFunction(IdentifySampleFunctor{}, inputData->getDataType(), imageGeom, inputData, fillHoles);
+  if(sliceBySlice)
+  {
+    ExecuteDataFunction(IdentifySampleSliceBySliceFunctor{}, inputData->getDataType(), imageGeom, inputData, fillHoles, sliceBySlicePlane);
+  }
+  else
+  {
+    ExecuteDataFunction(IdentifySampleFunctor{}, inputData->getDataType(), imageGeom, inputData, fillHoles);
+  }
 
   return {};
 }

src/Plugins/SimplnxCore/src/SimplnxCore/Filters/IdentifySampleFilter.hpp

@@ -22,6 +22,8 @@ class SIMPLNXCORE_EXPORT IdentifySampleFilter : public IFilter
 
   // Parameter Keys
   static inline constexpr StringLiteral k_FillHoles_Key = "fill_holes";
+  static inline constexpr StringLiteral k_SliceBySlice_Key = "slice_by_slice";
+  static inline constexpr StringLiteral k_SliceBySlicePlane_Key = "slice_by_slice_plane_index";
   static inline constexpr StringLiteral k_SelectedImageGeometryPath_Key = "input_image_geometry_path";
   static inline constexpr StringLiteral k_MaskArrayPath_Key = "mask_array_path";