WIP #14. Initial MarchingSquares algorithm. Working on connected component...

WIP #14. Initial MarchingSquares algorithm. Working on connected component labeling detect multiple groups.

radixalgorithm/marchingsquares.hh

 #ifndef RADIX_RADIXALGORITHM_MARCHINGSQUARES_HH_
 #define RADIX_RADIXALGORITHM_MARCHINGSQUARES_HH_
 
+#include <limits>
+#include <memory>
+#include <numeric>
 #include <vector>
 
+#include "radixbug/bug.hh"
+
 namespace radix
 {
-template <typename data_type, typename sign_type = int>
-void marching_squares(const std::vector<data_type>& points,
-                      std::vector<data_type>& out, size_t rows, size_t cols,
-                      data_type isovalue, sign_type negative = 0,
-                      sign_type positive = 1);
+template <typename data_type>
+class MarchingSquares
+{
+ protected:
+  std::vector<data_type> mData;
+  std::vector<int> mBit;
+  std::vector<int> mComponent;
+  size_t mRows;
+  size_t mColumns;
+
+  enum class StepDirection
+  {
+    Up,
+    Down,
+    Right,
+    Left,
+    None
+  };
+  StepDirection next_step;
+  StepDirection prev_step;
+  /**
+   * @brief starting_point find a starting position which matches or exceeds
+   * isovalue
+   * @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;
+  /**
+   * @brief step Performs a step in the march
+   * @param r row aka y index
+   * @param c column aka x index
+   */
+  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;
+  }
+
+ public:
+  /**
+   * @brief MarchingSquares
+   * @param data const std::vector<data_type> 1D representation of 2D data
+   * @param rows number of rows, aka y size
+   * @param columns number of columns, aka x size
+   */
+  MarchingSquares(const std::vector<data_type>& data, size_t rows,
+                  size_t columns)
+      : mData(data)
+      , mBit(data.size(), 0)
+      , mComponent(data.size())
+      , mRows(rows)
+      , mColumns(columns)
+  {
+    radix_line("MarchingSquares data size(" << data.size() << ")");
+    radix_check(data.size() == (rows * columns));
+  }
+
+  /**
+   * @brief march March the data on the isovalue threshold. Can be called
+   * multiple times for multiple grouping at threshold.
+   * @param isovalue the threshold value of the data
+   * @param wash_bit the value that replace all data matchings isovalue
+   * @param wash_threshold the threshold to which washing shouldn't apply.
+   * Primary use in multi-contour requests.
+   * @return std::vector<std::pair<int,int>> listing of all x,y indices, aka
+   * a closed polygon
+   */
+  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());
+};  // class
 }  // namespace
 
 /** Include implementation file */

radixalgorithm/marchingsquares.i.hh

 
+#include <array>
+#include <unordered_set>
 #include <vector>
 
+#include "radixalgorithm/marchingsquares.hh"
 #include "radixbug/bug.hh"
 
 namespace radix
 {
-template <typename data_type, typename sign_type = int>
-void marching_squares(const std::vector<data_type>& points,
-                      std::vector<data_type>& out, size_t rows, size_t cols,
-                      data_type isovalue, sign_type negative,
-                      sign_type positive)
+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)
 {
-  radix_check(points.size() == (rows * cols));
+  std::pair<size_t, size_t> start;
+  std::vector<std::pair<int, int>> out;
 
-  out.clear();
-  out.resize(points.size());
   //
-  // Create a 2d view of the points
-  auto view = [&points, cols](size_t row, size_t column) {
-    return points[cols * row + column];
-  };
-  for (size_t p_i = 0; p_i < points.size(); ++p_i)
+  // 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)
   {
-    // select data as 0 below value or 1 greater than
-    if (points[p_i] < isovalue)
+    radix_tagged_line("Did not find a starting point.");
+    // return empty listing
+    return out;
+  }
+
+  radix_tagged_line("Starting point[" << start.first << "," << start.second
+                                      << "]");
+  //
+  // Initialize bit field
+  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)
     {
-      out[p_i] = negative;
+      mBit[p_i] = 1;
     }
-    else
+  }
+
+  // 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++)
+    {
+      union_coords(row_i, col_i, row_i + 1, col_i);
+      union_coords(row_i, col_i, row_i, col_i + 1);
+    }
+  radix_block(
+      std::unordered_set<int> groups(mComponent.begin(), mComponent.end()));
+  radix_tagged_line("Number of groups(" << groups.size() << ")");
+
+  //
+  // Walk the perimeter
+  size_t row = start.first, column = start.second;
+  do
+  {
+    step(row, column);
+    // If our current point is within our image
+    // add it to the list of points
+    if (column >= 0 && column < mColumns && row >= 0 && row < mRows)
+      out.push_back({row, column});
+
+    switch (next_step)
     {
-      out[p_i] = positive;
+      case StepDirection::Up:
+        --row;
+        break;
+      case StepDirection::Left:
+        --column;
+        break;
+      case StepDirection::Down:
+        ++row;
+        break;
+      case StepDirection::Right:
+        ++column;
+        break;
+      default:
+        break;
     }
+  } while (row != start.first || column != start.second);
+
+  //
+  // Now we need to wash out selected group within contour, leaving those at or
+  // above threhold
+
+  return out;
+}  // march
+
+template <typename data_type>
+bool MarchingSquares<data_type>::starting_point(
+    std::pair<size_t, size_t>& pos) const
+{
+  for (size_t i = 0; i < mData.size(); ++i)
+  {
+    if (mBit[i] == 1)
+    {
+      // save the row
+      pos.first = i / mColumns;
+      // save the column
+      pos.second = i % mColumns;
+      return true;
+    }
+  }
+  return false;
+}
+template <typename data_type>
+bool MarchingSquares<data_type>::accepts(size_t r, size_t c) const
+{
+  // Make sure we don't pick a point out of bounds
+  if (c < 0 || r < 0 || c >= mColumns || r >= mRows) return false;
+
+  // Check the data value
+  if (mData[mColumns * r + c] > 0) return true;
+
+  return false;
+}
+
+template <typename data_type>
+void MarchingSquares<data_type>::step(size_t r, size_t c)
+{
+  //
+  // The meat of the marching squares algorithm
+  // See https://en.wikipedia.org/wiki/Marching_squares
+  // for specifics on the algorithm
+  // upper relations
+  bool u_left  = accepts(r - 1, c - 1);
+  bool u_right = accepts(r - 1, c);
+  // lower relations
+  bool l_left  = accepts(r, c - 1);
+  bool l_right = accepts(r, c);
+
+  prev_step = next_step;
+  int state = 0;
+
+  if (u_left) state |= 1;
+  if (u_right) state |= 2;
+  if (l_left) state |= 4;
+  if (l_right) state |= 8;
+  // State now contains a number between 0 and 15
+  // representing our state.
+  // In binary, it looks like 0000-1111 (in binary)
+
+  // An example. Let's say the top two pixels are filled,
+  // and the bottom two are empty.
+  // Stepping through the if statements above with a state
+  // of 0b0000 initially produces:
+  // upper Left == true ==>  0b0001
+  // upper Right == true ==> 0b0011
+  // The others are false, so 0b0011 is our state
+  // (That's a result 3)
+
+  // Looking at the chart above, we see that state
+  // corresponds to a move right, so in our switch statement
+  // below, we add a case for 3, and assign Right as the
+  // direction of the next step. We repeat this process
+  // for all 16 states.
+
+  // So we can use a switch statement to determine our
+  // next direction based on
+  switch (state)
+  {
+    case 1:
+      next_step = StepDirection::Up;
+      break;
+    case 2:
+      next_step = StepDirection::Right;
+      break;
+    case 3:
+      next_step = StepDirection::Right;
+      break;
+    case 4:
+      next_step = StepDirection::Left;
+      break;
+    case 5:
+      next_step = StepDirection::Up;
+      break;
+    case 6:
+      if (prev_step == StepDirection::Up)
+      {
+        next_step = StepDirection::Left;
+      }
+      else
+      {
+        next_step = StepDirection::Right;
+      }
+      break;
+    case 7:
+      next_step = StepDirection::Right;
+      break;
+    case 8:
+      next_step = StepDirection::Down;
+      break;
+    case 9:
+      if (prev_step == StepDirection::Right)
+      {
+        next_step = StepDirection::Up;
+      }
+      else
+      {
+        next_step = StepDirection::Down;
+      }
+      break;
+    case 10:
+      next_step = StepDirection::Down;
+      break;
+    case 11:
+      next_step = StepDirection::Down;
+      break;
+    case 12:
+      next_step = StepDirection::Left;
+      break;
+    case 13:
+      next_step = StepDirection::Up;
+      break;
+    case 14:
+      next_step = StepDirection::Left;
+      break;
+    default:
+      next_step = StepDirection::None;
+      break;
   }
 }
 }  // namespace

radixalgorithm/tests/CMakeLists.txt

 INCLUDE(GoogleTest)
 
 ADD_GOOGLE_TEST(tstOrdering.cc NP 1)
+ADD_GOOGLE_TEST(tstMarchingSquares.cc NP 1)

radixalgorithm/tests/MarchingSquares.Simple.txt

+0,0,0,0,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,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,2,2,0,0,2,2,2,0,0,0,0,0,0,0,0,0,0,0,
+0,0,2,0,0,0,0,0,2,2,2,0,0,0,0,0,0,0,0,0,
+0,0,2,0,0,0,0,0,0,0,2,2,2,0,0,0,0,0,0,0,
+0,0,2,0,0,0,0,0,0,0,0,2,2,0,0,0,0,0,0,0,
+0,0,2,2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
\ No newline at end of file

radixalgorithm/tests/MarchingSquares.TwoGroups.txt

+0,0,0,0,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,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,2,2,0,0,2,2,2,0,0,0,0,0,0,0,0,0,0,0,
+0,0,2,0,0,0,0,0,2,2,2,0,0,0,0,0,0,0,0,0,
+0,0,2,0,0,0,0,0,0,0,2,2,2,0,0,0,0,0,0,0,
+0,0,2,0,0,0,0,0,0,0,0,2,2,0,0,0,0,0,0,0,
+0,0,2,2,2,2,2,2,2,2,2,2,2,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,0,0,0,0
\ No newline at end of file

radixalgorithm/tests/tstMarchingSquares.cc

+#include <iostream>
+#include "gtest/gtest.h"
+
+#include "radixalgorithm/marchingsquares.hh"
+#include "radixbug/bug.hh"
+
+using namespace radix;
+TEST(MarchingSquares, DISABLED_Simple)
+{
+  //
+  // 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, 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, 2,
+      2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 2, 2, 2,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2, 2, 2, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 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.);
+  //  std::cout << "[ row , col ]" << std::endl;
+  //  std::cout << "{" << std::endl;
+  //  for (size_t i = 0; i < contour.size(); ++i)
+  //  {
+  //    const auto& pixel = contour[i];
+  //    std::cout << "{" << pixel.first << "," << pixel.second << "}";
+  //    if (i != contour.size() - 1)
+  //    {
+  //      std::cout << ",";
+  //    }
+  //    std::cout << std::endl;
+  //  }
+  //  std::cout << "}" << std::endl;
+
+  std::vector<std::pair<int, int>> blessed{
+      {1, 2},  {2, 2},  {3, 2},  {4, 2},  {5, 2},  {6, 2},  {7, 2},  {8, 2},
+      {8, 3},  {8, 4},  {8, 5},  {8, 6},  {8, 7},  {8, 8},  {8, 9},  {8, 10},
+      {8, 11}, {8, 12}, {8, 13}, {7, 13}, {6, 13}, {5, 13}, {5, 12}, {5, 11},
+      {4, 11}, {4, 10}, {4, 9},  {3, 9},  {3, 8},  {3, 7},  {2, 7},  {2, 6},
+      {2, 5},  {3, 5},  {3, 6},  {4, 6},  {4, 7},  {4, 8},  {5, 8},  {5, 9},
+      {5, 10}, {6, 10}, {6, 11}, {7, 11}, {7, 10}, {7, 9},  {7, 8},  {7, 7},
+      {7, 6},  {7, 5},  {7, 4},  {7, 3},  {6, 3},  {5, 3},  {4, 3},  {4, 4},
+      {3, 4},  {2, 4},  {1, 4},  {1, 3}};
+
+  ASSERT_EQ(blessed.size(), contour.size());
+  for (size_t i = 0; i < blessed.size(); ++i)
+  {
+    EXPECT_EQ(blessed[i], contour[i]);
+  }
+}
+
+TEST(MarchingSquares, DISABLED_TwoGroups)
+{
+  //
+  // 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, 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, 2,
+      2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 2, 2, 2,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2, 2, 2, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 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, 0, 0, 0, 0};
+  MarchingSquares<double> ms(grid, rows, columns);
+
+  {
+    std::vector<std::pair<int, int>> first_contour = ms.march(2.);
+
+    std::vector<std::pair<int, int>> first_blessed{
+        {1, 2},  {2, 2},  {3, 2},  {4, 2},  {5, 2},  {6, 2},  {7, 2},  {8, 2},
+        {8, 3},  {8, 4},  {8, 5},  {8, 6},  {8, 7},  {8, 8},  {8, 9},  {8, 10},
+        {8, 11}, {8, 12}, {8, 13}, {7, 13}, {6, 13}, {5, 13}, {5, 12}, {5, 11},
+        {4, 11}, {4, 10}, {4, 9},  {3, 9},  {3, 8},  {3, 7},  {2, 7},  {2, 6},
+        {2, 5},  {3, 5},  {3, 6},  {4, 6},  {4, 7},  {4, 8},  {5, 8},  {5, 9},
+        {5, 10}, {6, 10}, {6, 11}, {7, 11}, {7, 10}, {7, 9},  {7, 8},  {7, 7},
+        {7, 6},  {7, 5},  {7, 4},  {7, 3},  {6, 3},  {5, 3},  {4, 3},  {4, 4},
+        {3, 4},  {2, 4},  {1, 4},  {1, 3}};
+
+    ASSERT_EQ(first_blessed.size(), first_contour.size());
+    for (size_t i = 0; i < first_blessed.size(); ++i)
+    {
+      EXPECT_EQ(first_blessed[i], first_contour[i]);
+    }
+  }
+
+  {
+    std::vector<std::pair<int, int>> second_contour = ms.march(2.);
+    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;
+  }
+}