Merge branch 'ccl' into 'master'

Ccl

See merge request !56

radixalgorithm/marchingsquares.hh

@@ -44,6 +44,7 @@ class MarchingSquares
    */
   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
@@ -88,7 +89,33 @@ class MarchingSquares
       data_type isovalue, data_type wash_bit = {0},
       data_type wash_threshold = std::numeric_limits<data_type>::max());
 
+  /**
+   * @brief march March on the data using multiple isovalue thresholds.
+   * @param isotvalues List of threshold values of the data.
+   * @return list of polygons for each isovalue [isovalue i][polygon i][point i]
+   */
+  std::vector<std::vector<std::vector<std::pair<float, float>>>> march(
+      const std::vector<data_type>& isovalues, size_t max_contour_polygons = 10,
+      size_t max_polygon_points = 5000);
+
   void dump_component_map(std::ostream& os) const;
+
+ private:
+  /**
+   * @brief starting_point assumes labeling as occurred and finds a starting
+   * point which matches this label.
+   * @param pos std::pair<size_t, size_t> <row, column> aka <y index, x index>
+   * @param prev previous starting position. assumes clear_connected_components
+   * was called on prev
+   * @return bool true if a position was found
+   */
+  bool starting_point(short label, std::pair<size_t, size_t>& pos,
+                      std::pair<size_t, size_t> prev = {size_t(0),
+                                                        size_t(0)}) const;
+  void step(size_t r, size_t c, short label);
+  bool accepts(size_t r, size_t c, short label) const;
+  void clear_isovalue_label(int row, int col, short label,
+                            const std::vector<data_type>& ordered_isovalues);
 };  // class
 }  // namespace radix
 

radixalgorithm/marchingsquares.i.hh

 
+#include <queue>
 #include <set>
 #include <vector>
 
 #include "radixalgorithm/marchingsquares.hh"
+#include "radixalgorithm/ordering.hh"
 #include "radixbug/bug.hh"
 
 namespace radix
@@ -153,6 +155,28 @@ bool MarchingSquares<data_type>::starting_point(
   }
   return false;
 }
+
+template <typename data_type>
+bool MarchingSquares<data_type>::starting_point(
+    short label, std::pair<size_t, size_t>& pos,
+    std::pair<size_t, size_t> prev) const
+{
+  size_t data_size = mBit.size();
+  size_t starti    = mColumns * prev.first + prev.second;
+  for (size_t i = starti; i < data_size; ++i)
+  {
+    short value = mBit[i];
+    if (value == label)
+    {
+      // 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
 {
@@ -167,17 +191,24 @@ bool MarchingSquares<data_type>::accepts(size_t r, size_t c) const
 
 template <typename data_type>
 void MarchingSquares<data_type>::step(size_t r, size_t c)
+{
+  // step using the default label of 1
+  step(r, c, 1);
+}
+
+template <typename data_type>
+void MarchingSquares<data_type>::step(size_t r, size_t c, short label)
 {
   //
   // 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);
+  bool u_left  = accepts(r - 1, c - 1, label);
+  bool u_right = accepts(r - 1, c, label);
   // lower relations
-  bool l_left  = accepts(r, c - 1);
-  bool l_right = accepts(r, c);
+  bool l_left  = accepts(r, c - 1, label);
+  bool l_right = accepts(r, c, label);
 
   prev_step = next_step;
   int state = 0;
@@ -271,5 +302,228 @@ void MarchingSquares<data_type>::step(size_t r, size_t c)
                          : (next_step == StepDirection::Left) ? "Left"
                                                               : "Right")));
 }
+template <typename data_type>
+bool MarchingSquares<data_type>::accepts(size_t r, size_t c, short label) 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 (mBit[mColumns * r + c] >= label) return true;
+
+  return false;
+}
+template <typename data_type>
+std::vector<std::vector<std::vector<std::pair<float, float>>>>
+MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
+                                  size_t max_contour_polygons,
+                                  size_t max_polygon_points)
+{
+  std::pair<size_t, size_t> start, prev;
+
+  size_t isovalues_size = isovalues.size();
+  // [isovalue i][polygon i][point i]
+  std::vector<std::vector<std::vector<std::pair<float, float>>>> contours(
+      isovalues_size);
+  //
+  // Determine order to ensure largest to smallest ordering
+  // This ensures that we can deal with nested scenarios
+  auto comparator          = [](data_type a, data_type b) { return (a > b); };
+  std::vector<size_t> perm = sort_permutation(isovalues, comparator);
+  // apply ordering so a > b > c > d ...
+  std::vector<data_type> ordered_isovalues = isovalues;
+  apply_permutation(ordered_isovalues, perm);
+
+  //
+  // Initialize bit field
+  std::vector<size_t> label_counts(isovalues_size + 1, 0);
+  size_t data_size = mData.size();
+  for (size_t p_i = 0; p_i < data_size; ++p_i)
+  {
+    // select data as 1 if greater than isovalue
+    data_type value = mData[p_i];
+    for (size_t isoi = 0; isoi < isovalues_size; ++isoi)
+    {
+      if (value >= ordered_isovalues[isoi])
+      {
+        //
+        // save the category (aka there are only number of isovalue categories)
+        // specifically 1-N, 0 is background
+        mBit[p_i] = isovalues_size - isoi;
+        break;
+      }
+    }
+  }
+
+  typedef std::pair<float, float> Point;
+  typedef std::vector<Point> Polygon;
+  struct PolyComparator
+  {
+    bool operator()(const Polygon& a, const Polygon& b)
+    {
+      return a.size() < b.size();
+    }
+  };
+
+  for (size_t isoi = 0; isoi < isovalues_size; ++isoi)
+  {
+    //
+    // since we are using permutation as an interface, to ensure a > b
+    // the label will be reversed to the size
+    short label     = short(isovalues_size - isoi);
+    size_t contouri = size_t(label - 1);
+    //
+    // Get a starting point first, because if there isn't a place to start
+    // then we don't care about going any further for this isovalue
+
+    // reset previous for each contour search
+    prev = {0, 0};
+    std::priority_queue<Polygon, std::vector<Polygon>, PolyComparator> pqueue;
+    while (starting_point(label, start, prev))
+    {
+      prev = start;
+      // we have a new polygon
+      std::vector<std::pair<float, float>> polygon;
+
+      radix_tagged_line("Starting point("
+                        << ordered_isovalues[isoi] << ") label (" << label
+                        << ")[" << start.first << "," << start.second << "]");
+
+      //
+      // Walk the perimeter
+      size_t row = start.first, column = start.second;
+      do
+      {
+        step(row, column, label);
+        // If our current point is within our image
+        // add it to the list of points
+        // We have to allow for row and column being equal to the number
+        // of rows and columns to allow for traveling the boundaries
+        if (column >= 0 && column <= mColumns && row >= 0 && row <= mRows)
+        {
+          polygon.push_back({row, column});
+        }
 
+        switch (next_step)
+        {
+          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);
+
+      //
+      // Finish connecting the contour by making the last=first
+      polygon.push_back(start);
+      pqueue.push(polygon);
+
+      radix_tagged_line("clearing label ("
+                        << mBit[mColumns * start.first + start.second]
+                        << ") threshold (" << ordered_isovalues[isoi] << ")");
+
+      clear_isovalue_label(start.first, start.second,
+                           mBit[mColumns * start.first + start.second],
+                           ordered_isovalues);
+    }
+    //
+    // save the top N polygons
+    size_t count = std::min(max_contour_polygons, pqueue.size());
+    contours[contouri].resize(count);
+    for (size_t i = 0; i < count; ++i)
+    {
+      auto& top      = pqueue.top();
+      auto& cpolygon = contours[contouri][i];
+      if (top.size() > max_polygon_points)
+      {
+        // std::cout << "Taking every other point given size(" << top.size()
+        //          << ") exceeds " << max_polygon_points << std::endl;
+        // take every other point
+        size_t polygon_size = top.size() / 2;
+        for (size_t pi = 0; pi < polygon_size; ++pi)
+        {
+          cpolygon.push_back(top[pi * 2]);
+        }
+      }
+      else
+      {
+        cpolygon.resize(top.size());
+        std::copy(top.begin(), top.end(), cpolygon.begin());
+      }
+    }
+  }
+  return contours;
+}  // march
+
+template <typename data_type>
+void MarchingSquares<data_type>::clear_isovalue_label(
+    int row, int col, short label,
+    const std::vector<data_type>& ordered_isovalues)
+{
+  if (row < 0 || row == mColumns) return;  // out of bounds
+  if (col < 0 || col == mRows) return;     // out of bounds
+
+  std::set<size_t> list;
+  list.insert(mColumns * row + col);
+  size_t isovalues_size = ordered_isovalues.size();
+  size_t isoi           = isovalues_size - size_t(label) + 1;
+  while (!list.empty())
+  {
+    const auto& it = list.begin();
+    size_t c_i     = *it;
+
+    // update the row
+    row = c_i / mColumns;
+    // upate the column
+    col = c_i % mColumns;
+    // re-initilize bit field
+    data_type value = mData[c_i];
+    if (label > 1)  // not the last isovalue
+    {
+      // radix_tagged_line("Comparing " << value
+      //                               << " >= " << ordered_isovalues[isoi]);
+      if (value >= ordered_isovalues[isoi])
+      {
+        //
+        // save the category (aka there are only number of isovalue categories)
+        // specifically 1-N, 0 is background
+        mBit[c_i] = label - 1;
+        // radix_tagged_line("Found new isovalue (" << ordered_isovalues[isoi]
+        //                                         << ") for [" << row << ","
+        //                                         << col << "]");
+      }
+    }
+    // we didn't find a lower level contour, so zero it out
+    if (mBit[c_i] == label) mBit[c_i] = 0;
+    // search neighbors
+    for (int direction = 0; direction < 4; ++direction)
+    {
+      int nc = col + dx[direction];
+      int nr = row + dy[direction];
+      if (nc < 0 || nc >= mColumns) continue;  // out of bounds
+      if (nr < 0 || nr >= mRows) continue;     // out of bounds
+      size_t nc_i = mColumns * nr + nc;
+      if (mBit[nc_i] == label)
+      {
+        // if we already have this cell in the list to look at
+        // don't add it again
+        if (list.find(nc_i) == list.end())
+        {
+          list.insert(nc_i);
+        }
+      }
+    }
+    list.erase(it);
+  }
+}
 }  // namespace radix