Marching Squares in C++

marching_squares.cc

#include "marching-squares.hh"

#include "point.hh"
#include "polygon.hh"
#include "douglas-peucker.hh"
#include <ps/canvas/canvas.hh>

#include <cstring>
#include <vector>

namespace Ps {

  /* Implementation based on: "Better Know An Algorithm 1: Marching Squares" */
  /* "http://devblog.phillipspiess.com/2010/02/23/"
   *  "better-know-an-algorithm-1-marching-squares/" 
   */
  
  namespace {

    typedef char Direction;

    /* We build these as char[4]'s, but can then compare them as ints */
    union Key {
      char c[4];
      int  i;

      Key();
      Key(const char* s);
    };

    inline Key::Key()
      : i(0)
    {
      ;
    }

    inline Key::Key(const char* s)
      : i(0)
    {
      assert(std::strlen(s) >= sizeof(c));
      std::memcpy(c, s, sizeof(c));
    }


    struct Transition {
      Transition() = default;
      
      /* implicit */
      Transition(const char* desc);

      Key       key = Key();
      Direction dir = 'X';
    };

    Transition::Transition(const char* desc)
      : key(desc), dir('X')
    {
      assert(desc != nullptr && strlen(desc) == 6);
      dir = desc[6];
    }

    static const Transition transitions[] = {
      "0000-R", "1100-R", "0011-L", "1110-R",
      "0100-R", "1000-U", "1011-U", "1001-U",
      "0101-D", "1010-U", "0111-L", "0110-L",
      "0001-D", "0010-L", "1101-D", "1111-R"
    };
  }



  struct MarchingSquares::Impl
  {
    Impl(MarchingSquares* inst);
    
  public:
    MarchingSquares*  inst = nullptr;
    size_t            w    = 0;
    size_t            h    = 0;

    /* Full knowledge of vector<bool> specialization */
    std::vector<bool> data;
  };

  MarchingSquares::Impl::Impl(MarchingSquares* inst)
    : inst(inst)
  {
    ;
  }


  MarchingSquares::MarchingSquares(Canvas* canvas)
    : impl(new Impl(this))
  {
    impl->w = canvas->size().w();
    impl->h = canvas->size().h();
    impl->data.reserve(impl->w * impl->h);

    for (size_t y = 0; y != impl->h; ++y) {
      for (size_t x = 0; x != impl->w; ++x) {
        Canvas::Pixel p = canvas->pixel_at(Point(x, y));
        impl->data[y * impl->w + x] = ((p & 0xff000000) >> 24) > 0;
      }
    }
  }

  MarchingSquares::~MarchingSquares()
  {
    delete impl;
  }

  Polygon MarchingSquares::extract_simple(Float epsilon)
  {
    Polygon poly;

    size_t x = 0;
    size_t y = 0;

    /* Find first edge */
    for (int i = 0; i != impl->data.size(); ++i) {
      if (impl->data[i]) {
        y = i / impl->w;
        x = i % impl->w;
        break;
      }
    }

    /* If we got to the bottom-right, we're done: no result */
    if (x == impl->w && y == impl->h) {
      return Polygon();
    }

    /* save terminal condition */
    size_t fx = x, fy = y;

    Transition window;
    char       last_dir = 'X';

    do {
      /* Lets us treat image boundary as 0. */
      #define VIEW(ox, oy) ((ox < 0 || oy < 0 || ox == impl->w || oy == impl->h) \
                             ? '0' : impl->data[(oy) * impl->w + (ox)])
      window.key.c[0] = VIEW(x-1, y-1);
      window.key.c[1] = VIEW(  x, y-1);
      window.key.c[2] = VIEW(x-1,   y);
      window.key.c[3] = VIEW(  x,   y);
      #undef VIEW

      /* Based on the window we're looking at, find the next step */
      Direction dir = 'X';
      for (auto& st: transitions) {
        if (window.key.i == st.key.i) {
          dir = st.dir;
          break;
        }
      }

      /* These two handle special cases that can cause an infinite loop */
      if (window.key.i == transitions[11].key.i) {
        assert(dir == 'L');
        dir = (last_dir == 'U' ? 'L' : 'R');
      } else if (window.key.i == transitions[7].key.i) {
        assert(dir == 'U');
        dir = (last_dir == 'R' ? 'U' : 'D');
      }

      assert(dir != 'X');
      if (dir != last_dir) {
        poly.push_back(Point(x, y));
      }

      switch (dir) {
        case 'U': --y; break;
        case 'D': ++y; break;
        case 'L': --x; break;
        case 'R': ++x; break;
        default:
          assert(false);
      }

      last_dir = dir;
    } while (! (y == fy && x == fx));

    assert(epsilon >= 0.0f);
    if (epsilon > 0.0f) {
      return DouglasPeucker::simplify(poly, epsilon);
    } else {
      return poly;
    }

  }
}

marching_squares.hh

#ifndef PS_GEOM_MARCHING_SQUARES_HH
#define PS_GEOM_MARCHING_SQUARES_HH

namespace Ps {
  class Canvas;
  class Polygon;
  class Mesh2d;

  class MarchingSquares
  {
  public:
    MarchingSquares(Canvas* surface);
    virtual ~MarchingSquares();

    /* Extracts a polygon, then simplifies to simplification_epsilon */
    Polygon extract_simple(Float simplification_epsilon = PS_PIXEL_EPSILON);

  private:
    struct Impl;
    Impl* impl;
  };
}

#endif

Oh, I could see how that could help clean up the contours... Thanks for the reply. If you ever do post that code, I'll gladly take a look at it. Thinking about using your code for a project. - dan

…

On Wed, Sep 4, 2019 at 8:14 AM Mike Owens ***@***.***> wrote: Hi, Daniel. It removes points that are nearly-colinear on each edge of a polygon. The epsilon is the “tolerance”. I obviously had a C++ implementation at one point, but I see I have a old Ruby gist with the ugly “textbook” implementation: https://gist.github.com/mieko/7c1f5a3b51423c4322e2 I’ll see if I can dig up the C++ version later. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <https://gist.github.com/3c17c41b8be9f2598c12?email_source=notifications&email_token=AAED27M2EIO6LQZA2LMKW4DQH6YCZA5CNFSM4ITRC5V2YY3PNVWWK3TUL52HS4DFVNDWS43UINXW23LFNZ2KUY3PNVWWK3TUL5UWJTQAFYDO2#gistcomment-3016429>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAED27MKV4ZJFWQ3WMSW7DTQH6YCZANCNFSM4ITRC5VQ> .

I've resurrected the C++ implementation out of a really-old archive, available here:

https://gist.github.com/mieko/0275f2f4a3b18388ed5131b3364179fb

Not drop-in ready (uses references geometry types undefined), but the missing parts are trivial.

-M

Thank you, sir. I'll check it out.