WIP: addressing boundary condition in connected components isovalue elimination.

radixalgorithm/marchingsquares.i.hh

 
-#include <queue>
 #include <set>
 #include <vector>
 
@@ -336,7 +335,6 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
 
   //
   // 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)
   {
@@ -357,12 +355,8 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
 
   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();
-    }
+  auto poly_comparator = [](const Polygon& a, const Polygon& b) {
+    return (a.size() < b.size());
   };
 
   for (size_t isoi = 0; isoi < isovalues_size; ++isoi)
@@ -378,7 +372,7 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
 
     // reset previous for each contour search
     prev = {0, 0};
-    std::priority_queue<Polygon, std::vector<Polygon>, PolyComparator> pqueue;
+    std::vector<Polygon> pvector;
     while (starting_point(label, start, prev))
     {
       prev = start;
@@ -426,11 +420,28 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
       //
       // Finish connecting the contour by making the last=first
       polygon.push_back(start);
-      pqueue.push(polygon);
+      if (pvector.size() == max_contour_polygons)
+      {
+        // compare against the last element
+        // if it is larger then replace it.
+        if (pvector[pvector.size() - 1].size() < polygon.size())
+        {
+          pvector[pvector.size() - 1] = polygon;
+          std::sort(pvector.begin(), pvector.end(), poly_comparator);
+        }
+        // we ignore as this polygon is not large enough to keep
+      }
+      else
+      {
+        pvector.push_back(polygon);
+        std::sort(pvector.begin(), pvector.end(), poly_comparator);
+      }
 
       radix_tagged_line("clearing label ("
                         << mBit[mColumns * start.first + start.second]
-                        << ") threshold (" << ordered_isovalues[isoi] << ")");
+                        << ") threshold (" << ordered_isovalues[isoi]
+                        << ") starting at [" << start.first << ","
+                        << start.second << "]");
 
       clear_isovalue_label(start.first, start.second,
                            mBit[mColumns * start.first + start.second],
@@ -438,30 +449,28 @@ MarchingSquares<data_type>::march(const std::vector<data_type>& isovalues,
     }
     //
     // save the top N polygons
-    size_t count = std::min(max_contour_polygons, pqueue.size());
+    size_t count = pvector.size();
     contours[contouri].resize(count);
     for (size_t i = 0; i < count; ++i)
     {
-      auto& top      = pqueue.top();
+      auto& polygon  = pvector[i];
       auto& cpolygon = contours[contouri][i];
-      if (top.size() > max_polygon_points)
+      if (polygon.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;
+        size_t polygon_size = polygon.size() / 2;
         for (size_t pi = 0; pi < polygon_size; ++pi)
         {
-          cpolygon.push_back(top[pi * 2]);
+          cpolygon.push_back(polygon[pi * 2]);
         }
       }
       else
       {
-        cpolygon.resize(top.size());
-        std::copy(top.begin(), top.end(), cpolygon.begin());
+        cpolygon.resize(polygon.size());
+        std::copy(polygon.begin(), polygon.end(), cpolygon.begin());
       }
-      // remove the top element
-      pqueue.pop();
     }
   }
   return contours;
@@ -472,8 +481,9 @@ 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
+  // we can allow for row=mRows and col=mColumns
+  if (row < 0 || row > mRows) return;     // out of bounds
+  if (col < 0 || col > mColumns) return;  // out of bounds
 
   std::set<size_t> list;
   list.insert(mColumns * row + col);
@@ -488,25 +498,6 @@ void MarchingSquares<data_type>::clear_isovalue_label(
     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)
     {
@@ -521,10 +512,36 @@ void MarchingSquares<data_type>::clear_isovalue_label(
         // don't add it again
         if (list.find(nc_i) == list.end())
         {
+          radix_tagged_line("\t[" << nr << "," << nc << "]");
           list.insert(nc_i);
         }
       }
     }
+    // re-initilize bit field
+    if (label > 1)  // not the last isovalue
+    {
+      data_type value = mData[c_i];
+      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
+    radix_tagged_line("Comparing mBit[" << c_i << "](" << mBit[c_i] << ") to "
+                                        << label);
+    if (mBit[c_i] == label)
+    {
+      radix_tagged_line("\tzeroing [" << row << "," << col << "]");
+      mBit[c_i] = 0;
+    }
     list.erase(it);
   }
 }

radixalgorithm/tests/tstMarchingSquares.cc

@@ -276,7 +276,7 @@ TEST(MarchingSquares, NestedMultiIsovalues)
   {
     std::string line;
     std::getline(stream, line);
-    std::vector<std::string> fields = radix::split_string(",", line);
+    std::vector<std::string> fields = radix::split_string(",", line, true);
     for (size_t i = 0; i < fields.size(); ++i)
     {
       double value = std::atof(fields[i].c_str());
@@ -287,6 +287,7 @@ TEST(MarchingSquares, NestedMultiIsovalues)
   }
   size_t rows    = 256;
   size_t columns = rows;
+  radix_tagged_line("Grid size(" << grid.size() << ")");
   MarchingSquares<double> ms(grid, rows, columns);
 
   std::vector<std::vector<std::vector<std::pair<float, float>>>> contours;