Initial working version of connect components labeling within MarchingSquares.

radixalgorithm/marchingsquares.hh

@@ -4,6 +4,7 @@
 #include <limits>
 #include <memory>
 #include <numeric>
+#include <ostream>
 #include <vector>
 
 #include "radixbug/bug.hh"
@@ -13,6 +14,10 @@ namespace radix
 template <typename data_type>
 class MarchingSquares
 {
+ private:
+  std::vector<int> dx;
+  std::vector<int> dy;
+
  protected:
   std::vector<data_type> mData;
   std::vector<int> mBit;
@@ -36,7 +41,8 @@ class MarchingSquares
    * @param pos std::pair<size_t, size_t> <row, column> aka <y index, xindex>
    * @return bool true if a position was found
    */
-  bool starting_point(std::pair<size_t, size_t>& pos) const;
+  bool starting_point(data_type isovalue, data_type wash_threshold,
+                      std::pair<size_t, size_t>& pos) const;
   /**
    * @brief step Performs a step in the march
    * @param r row aka y index
@@ -44,21 +50,8 @@ class MarchingSquares
    */
   void step(size_t r, size_t c);
   bool accepts(size_t r, size_t c) const;
-
-  void union_coords(size_t r1, size_t c1, size_t r2, size_t c2)
-  {
-    size_t first  = mColumns * c1 + r1;
-    size_t second = mColumns * c2 + r2;
-    if (r2 < mRows && c2 < mColumns && mBit[first] && mBit[second])
-      do_union(int(first), int(second));
-  }
-  void do_union(int a, int b)
-  {
-    // get the root component of a and b, and set the one's parent to the other
-    while (mComponent[a] != a) a = mComponent[a];
-    while (mComponent[b] != b) b = mComponent[b];
-    mComponent[b] = a;
-  }
+  void find_components();
+  void dfs(int x, int y, size_t current_label);
 
  public:
   /**
@@ -71,12 +64,13 @@ class MarchingSquares
                   size_t columns)
       : mData(data)
       , mBit(data.size(), 0)
-      , mComponent(data.size())
       , mRows(rows)
       , mColumns(columns)
   {
-    radix_line("MarchingSquares data size(" << data.size() << ")");
+    //    radix_line("MarchingSquares data size(" << data.size() << ")");
     radix_check(data.size() == (rows * columns));
+    dx = {+1, 0, -1, 0};
+    dy = {0, +1, 0, -1};
   }
 
   /**
@@ -92,8 +86,10 @@ class MarchingSquares
   std::vector<std::pair<int, int>> march(
       data_type isovalue, data_type wash_bit = {0},
       data_type wash_threshold = std::numeric_limits<data_type>::max());
+
+  void dump_component_map(std::ostream& os) const;
 };  // class
-}  // namespace
+}  // namespace radix
 
 /** Include implementation file */
 #include "radixalgorithm/marchingsquares.i.hh"

radixalgorithm/marchingsquares.i.hh

@@ -8,9 +8,43 @@
 
 namespace radix
 {
+template <typename data_type>
+void MarchingSquares<data_type>::dfs(int x, int y, size_t current_label)
+{
+  if (x < 0 || x == mColumns) return;  // out of bounds
+  if (y < 0 || y == mRows) return;     // out of bounds
+  size_t c_i = mColumns * y + x;
+
+  if (mComponent[c_i] || !mBit[c_i])
+    return;  // already labeled or not marked with 1 in m
+
+  // mark the current cell
+  mComponent[c_i] = current_label;
+
+  // recursively mark the neighbors
+  for (int direction = 0; direction < 4; ++direction)
+    dfs(x + dx[direction], y + dy[direction], current_label);
+}
+
+template <typename data_type>
+void MarchingSquares<data_type>::find_components()
+{
+  mComponent.resize(mBit.size(), 0);
+  size_t component = 0;
+  for (size_t c_i = 0; c_i < mComponent.size(); ++c_i)
+  {
+    // save the row
+    size_t y = c_i / mColumns;
+    // save the column
+    size_t x = c_i % mColumns;
+    // check if it is already labeled
+    if (!mComponent[c_i] && mBit[c_i]) dfs(x, y, ++component);
+  }
+}
+
 template <typename data_type>
 std::vector<std::pair<int, int>> MarchingSquares<data_type>::march(
-    data_type isovalue, data_type wash_bit, data_type wash_thresold)
+    data_type isovalue, data_type wash_bit, data_type wash_threshold)
 {
   std::pair<size_t, size_t> start;
   std::vector<std::pair<int, int>> out;
@@ -18,7 +52,7 @@ std::vector<std::pair<int, int>> MarchingSquares<data_type>::march(
   //
   // Get a starting point first, because if there isn't a place to start
   // then we don't care about going any further
-  if (starting_point(start) == false)
+  if (starting_point(isovalue, wash_threshold, start) == false)
   {
     radix_tagged_line("Did not find a starting point.");
     // return empty listing
@@ -32,22 +66,23 @@ std::vector<std::pair<int, int>> MarchingSquares<data_type>::march(
   for (size_t p_i = 0; p_i < mData.size(); ++p_i)
   {
     // select data as 1 if greater than isovalue
-    if (mData[p_i] >= isovalue)
+    if (mData[p_i] >= isovalue && mData[p_i] < wash_threshold)
     {
       mBit[p_i] = 1;
     }
-  }
-
-  // We required connect commponent labeling (CCL)
-  // to deal with multiple groups
-  // initialize all cells as unique component 0 -> data.size()
-  std::iota(std::begin(mComponent), std::end(mComponent), 0);
-  for (size_t row_i = 0; row_i < mRows; row_i++)
-    for (size_t col_i = 0; col_i < mColumns; col_i++)
+    else if (mData[p_i] >= wash_threshold)
     {
-      union_coords(row_i, col_i, row_i + 1, col_i);
-      union_coords(row_i, col_i, row_i, col_i + 1);
+      mBit[p_i] = 2;
     }
+  }
+  //
+  // connected component labeling is required for multiple groups
+  if (mComponent.empty()) find_components();
+  size_t group_i = mColumns * start.first + start.second;
+  radix_line("Starting point is a member of component(" << mComponent[group_i]
+                                                        << ")");
+
+  std::cerr << std::endl;
   radix_block(
       std::unordered_set<int> groups(mComponent.begin(), mComponent.end()));
   radix_tagged_line("Number of groups(" << groups.size() << ")");
@@ -85,17 +120,25 @@ std::vector<std::pair<int, int>> MarchingSquares<data_type>::march(
   //
   // Now we need to wash out selected group within contour, leaving those at or
   // above threhold
-
+  // Get the group the starting point is a member of
+  for (size_t i = 0; i < mData.size(); ++i)
+  {
+    if (mComponent[i] == mComponent[group_i])
+    {
+      mData[i] = wash_bit;
+    }
+  }
   return out;
 }  // march
 
 template <typename data_type>
 bool MarchingSquares<data_type>::starting_point(
+    data_type isovalue, data_type wash_threshold,
     std::pair<size_t, size_t>& pos) const
 {
   for (size_t i = 0; i < mData.size(); ++i)
   {
-    if (mBit[i] == 1)
+    if (mData[i] >= isovalue && mData[i] < wash_threshold)
     {
       // save the row
       pos.first = i / mColumns;
@@ -223,4 +266,15 @@ void MarchingSquares<data_type>::step(size_t r, size_t c)
       break;
   }
 }
-}  // namespace
+
+template <typename data_type>
+void MarchingSquares<data_type>::dump_component_map(std::ostream& os) const
+{
+  for (size_t i = 0; i < mComponent.size(); ++i)
+  {
+    os << std::setw(4) << mComponent[i] << ", ";
+    if (((i + 1) % mColumns) == 0) os << std::endl;
+  }
+  os << std::endl;
+}
+}  // namespace radix

radixalgorithm/tests/MarchingSquares.4EdgeGroups.txt

+0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,5,5,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,5,5,0,0,0,
+2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,4,
+2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,4,4,
+2,2,0,0,0,0,3,3,0,0,0,0,0,0,0,0,0,4,4,0,
+2,2,0,0,0,0,3,3,0,0,0,0,0,0,0,0,0,4,4,0,
+0,0,0,0,0,0,3,3,0,0,0,0,0,0,0,0,0,0,4,4,
+0,0,0,3,3,3,3,3,0,0,0,0,0,0,0,0,0,0,0,0

radixalgorithm/tests/tstMarchingSquares.cc

@@ -5,7 +5,7 @@
 #include "radixbug/bug.hh"
 
 using namespace radix;
-TEST(MarchingSquares, DISABLED_Simple)
+TEST(MarchingSquares, Simple)
 {
   //
   // define 10x20 grid of values
@@ -55,7 +55,7 @@ TEST(MarchingSquares, DISABLED_Simple)
   }
 }
 
-TEST(MarchingSquares, DISABLED_TwoGroups)
+TEST(MarchingSquares, TwoGroups)
 {
   //
   // define 10x20 grid of values
@@ -93,20 +93,73 @@ TEST(MarchingSquares, DISABLED_TwoGroups)
     }
   }
 
-  {
+  {  // get the second grouping
     std::vector<std::pair<int, int>> second_contour = ms.march(2.);
-    std::cout << "[ row , col ]" << std::endl;
+    //    std::cout << "[ row , col ]" << std::endl;
+    //    std::cout << "{" << std::endl;
+    //    for (size_t i = 0; i < second_contour.size(); ++i)
+    //    {
+    //      const auto& pixel = second_contour[i];
+    //      std::cout << "{" << pixel.first << "," << pixel.second << "}";
+    //      if (i != second_contour.size() - 1)
+    //      {
+    //        std::cout << ",";
+    //      }
+    //      std::cout << std::endl;
+    //    }
+    //    std::cout << "}" << std::endl;
+    std::vector<std::pair<int, int>> second_blessed{
+        {7, 16}, {8, 16}, {9, 16}, {9, 17}, {9, 18}, {8, 18}, {7, 18}, {7, 17}};
+    ASSERT_EQ(second_blessed.size(), second_contour.size());
+    for (size_t i = 0; i < second_blessed.size(); ++i)
+    {
+      EXPECT_EQ(second_blessed[i], second_contour[i]);
+    }
+  }
+  //
+  // Expect that there is no
+  EXPECT_EQ(0, ms.march(2.).size());
+}
+
+TEST(MarchingSquares, 4EdgeGroups)
+{
+  //
+  // define 10x20 grid of values
+  size_t rows{10};
+  size_t columns{20};
+  std::vector<double> grid{
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 5, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 5, 0, 0, 0, 2, 2, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 4, 4, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 4, 4, 4, 2, 2, 0, 0, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,
+      4, 0, 2, 2, 0, 0, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 0, 0,
+      0, 0, 0, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 0, 0, 0, 3,
+      3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  MarchingSquares<double> ms(grid, rows, columns);
+
+  {
+    std::vector<std::pair<int, int>> contour = ms.march(2., 0., 3.);
+    ms.dump_component_map(std::cout);
     std::cout << "{" << std::endl;
-    for (size_t i = 0; i < second_contour.size(); ++i)
+    for (size_t i = 0; i < contour.size(); ++i)
     {
-      const auto& pixel = second_contour[i];
+      const auto& pixel = contour[i];
       std::cout << "{" << pixel.first << "," << pixel.second << "}";
-      if (i != second_contour.size() - 1)
+      if (i != contour.size() - 1)
       {
         std::cout << ",";
       }
       std::cout << std::endl;
     }
     std::cout << "}" << std::endl;
+    std::vector<std::pair<int, int>> blessed;
+
+    //    ASSERT_EQ(blessed.size(), contour.size());
+    for (size_t i = 0; i < blessed.size(); ++i)
+    {
+      EXPECT_EQ(blessed[i], contour[i]);
+    }
   }
 }