Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Marching Squares in C++
#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;
}
}
}
#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
@danelliottster

This comment has been minimized.

Copy link

@danelliottster danelliottster commented Sep 4, 2019

What does DouglasPeucker::simplify(poly, epsilon); do?

@mieko

This comment has been minimized.

Copy link
Owner Author

@mieko mieko commented Sep 4, 2019

Hi, Dan.

It removes points that are nearly-colinear on each edge of a polygon. The epsilon is the “tolerance”. Basically noise reduction on extraneous points.

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.

@danelliottster

This comment has been minimized.

Copy link

@danelliottster danelliottster commented Sep 4, 2019

@mieko

This comment has been minimized.

Copy link
Owner Author

@mieko mieko commented Sep 10, 2019

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

@danelliottster

This comment has been minimized.

Copy link

@danelliottster danelliottster commented Sep 10, 2019

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment