WIP. starting on marchingsquares algorithm.

radixalgorithm/CMakeLists.txt

@@ -5,6 +5,8 @@ ordering.cc
 )
 SET(HEADERS
 ordering.hh
+marchingsquares.hh
+marchingsquares.i.hh
 )
 
 TRIBITS_ADD_LIBRARY(radixalgorithmlib

radixalgorithm/marchingsquares.hh

+#ifndef RADIX_RADIXALGORITHM_MARCHINGSQUARES_HH_
+#define RADIX_RADIXALGORITHM_MARCHINGSQUARES_HH_
+
+#include <vector>
+
+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);
+}  // namespace
+
+/** Include implementation file */
+#include "radixalgorithm/marchingsquares.i.hh"
+#endif /**  RADIX_RADIXALGORITHM_MARCHINGSQUARES_HH_ */

radixalgorithm/marchingsquares.i.hh

+
+#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,
+                      sign_type positive)
+{
+  radix_check(points.size() == (rows * cols));
+
+  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)
+  {
+    // select data as 0 below value or 1 greater than
+    if (points[p_i] < isovalue)
+    {
+      out[p_i] = negative;
+    }
+    else
+    {
+      out[p_i] = positive;
+    }
+  }
+}
+}  // namespace