Marching Squares

In this post I’ll walk through some features of the par_msquares.h library. The API has three entry points:

par_msquares_meshlist* par_msquares_grayscale(float const* data, int width,
    int height, int cellsize, float threshold, int flags);

par_msquares_meshlist* par_msquares_color(uint8_t const* data, int width,
    int height, int cellsize, uint32_t color, int bpp, int flags);

par_msquares_meshlist* par_msquares_color_multi(uint8_t const* data, int width,
    int height, int cellsize, int bpp, int flags);

All of these functions consume a packed array of image data and produce one or more triangle meshes. The passed-in cell size determines the density of the generated meshes. The results are high quality in the following sense:

During the march, meshes are constructed such that neighboring triangles share vertices as much as possible.
Generated edges are nudged at pixel increments to better match the original contour in the image.

The msquares_grayscale function consumes floating-point data (one 32-bit float per pixel), while the msquares_color function consumes color data (one to four bytes per pixel).

For the grayscale function, a threshold is passed in to determine insideness. For the color-based functions, color determines insideness. If you’ve got an image with less than 256 colors, you can generate a separate mesh for each color using msquares_color_multi.

The library also proffers a lower-level function that takes a callback instead of image data. Look at the source to see how to use it.

All imperative functions return an opaque mesh list pointer. Clients can peek at read-only mesh structures by passing the list into a query function:

par_msquares_mesh const* par_msquares_get_mesh(par_msquares_meshlist*, int n);

int par_msquares_get_count(par_msquares_meshlist*);

typedef struct {
    float* points;        // pointer to XY (or XYZ) vertex coordinates
    int npoints;          // number of vertex coordinates
    uint16_t* triangles;  // pointer to 3-tuples of vertex indices
    int ntriangles;       // number of 3-tuples
    int dim;              // number of floats per point (either 2 or 3)
    uint32_t color;       // used only with par_msquares_color_multi
} par_msquares_mesh;

See the unit tests for an example of how to export the above data structure to an OBJ mesh.

The polygonal contour of the mesh can be obtained using the following utility function.

par_msquares_boundary* par_msquares_extract_boundary(par_msquares_mesh const* );

// Polyline boundary extracted from a mesh, composed of one or more chains.
// Counterclockwise chains are solid, clockwise chains are holes.  So, when
// serializing to SVG, all chains can be aggregated into a single <path>,
// provided each one terminates with a "Z" and uses the default fill rule.
typedef struct {
    float* points;        // list of XY vertex coordinates
    int npoints;          // number of vertex coordinates
    float** chains;       // list of pointers to the start of each chain
    uint16_t* lengths;    // list of chain lengths
    int nchains;          // number of chains
} par_msquares_boundary;

See the unit tests for an example of how to export the above data structure to a SVG file.

Examples with a Grayscale Image

Let’s say your source data is a floating-point image that looks like this. It is shown with a perspective tilt for reasons that will be apparent later.

The following SVG image was created by calling par_msquares_grayscale with the above input data. It was invoked four times using various threshold values. After each invocation, the contours were obtained using par_msquares_extract_boundary and appended to the SVG file.

Next we’ll have a visual walkthrough of the flags argument, which takes a bitfield that can be configured in many different ways.

The simplest thing to do is to use the default behavior. Pass 0 for the flags argument, pick a threshold, and pick a cell size.

int flags = 0;
mlist = par_msquares_grayscale(graydata, width, height, cellsize, thresh,
    flags);

The above image appears squashed because we’re applying some perspective — stay tuned.

The INVERT flag does what it sounds like.

mlist = par_msquares_grayscale(graydata, width, height, cellsize, thresh,
    PAR_MSQUARES_INVERT);

The SIMPLIFY flag performs some really basic simplification with no loss to the quality of the boundary. It won’t produce the simplest possible mesh, but it’s fast and simple.

mlist = par_msquares_grayscale(graydata, width, height, cellsize, thresh,
    PAR_MSQUARES_SIMPLIFY);

The DUAL flag causes the returned meshlist to have two entries instead of one. The two meshes are disjoint; the boundary verts of each mesh are perfectly colocated.

mlist = par_msquares_grayscale(graydata, width, height, cellsize, thresh,
    PAR_MSQUARES_DUAL);

The HEIGHTS flag generates a Z value at each vert by sampling from the nearest pixel in the source image.

mlist = par_msquares_grayscale(graydata, width, height, cellsize, thresh,
    PAR_MSQUARES_HEIGHTS);

When the HEIGHTS flag is combined with the SNAP flag, a step function is applied to the Z values such that the number of steps is equal to the number of meshes in the generated mesh list.

mlist = par_msquares_grayscale(graydata, width, height, cellsize, thresh,
    PAR_MSQUARES_DUAL | PAR_MSQUARES_HEIGHTS | PAR_MSQUARES_SNAP);

The CONNECT flag adds triangles that connect the disjoint verts along the boundary.

mlist = par_msquares_grayscale(graydata, width, height, cellsize, thresh,
    PAR_MSQUARES_DUAL | PAR_MSQUARES_HEIGHTS | PAR_MSQUARES_SNAP |
    PAR_MSQUARES_CONNECT);

Color Source Images

Grayscale images use a threshold value to determine insideness, but color images use a simple color selector. For example, suppose we have an ARGB image that looks like this:

The transparent pixels on the outside of the island have an ARGB value of 0x00214562, so we can invoke marching squares like this:

mlist = par_msquares_color(pixels, width, height, cellsize, 0x214562, 4, 0);

Let’s try combining INVERT and HEIGHTS. The alpha values in the source image are wired into the resulting Z values.

mlist = par_msquares_color(pixels, width, height, cellsize, 0x214562, 4,
    PAR_MSQUARES_INVERT | PAR_MSQUARES_HEIGHTS);

Next, let’s try combining a ton of flags.

mlist = par_msquares_color(pixels, width, height, cellsize, 0x214562, 4,
    PAR_MSQUARES_DUAL | PAR_MSQUARES_HEIGHTS | PAR_MSQUARES_SNAP |
    PAR_MSQUARES_CONNECT | PAR_MSQUARES_SIMPLIFY | PAR_MSQUARES_INVERT);

One more example, this time using the color_multi entry point.

mlist = par_msquares_color_multi(pixels, width, height, cellsize, 4,
    PAR_MSQUARES_HEIGHTS | PAR_MSQUARES_CONNECT | PAR_MSQUARES_SIMPLIFY);

Footnote

C is my favorite language nowadays, especially the C99 and C11 variants. I try not to be stubborn about it – it’s obviously not the best choice in many situations. But, I definitely enjoy working with it more than any other language.

Inspired by Sean Barrett’s sweet stb libraries, I’ve started to create my own little collection of single-header libraries, and that’s where the msquares library lives: