Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
2118 lines (1970 sloc) 77.4 KB
// MSQUARES :: https://github.com/prideout/par
// Converts fp32 grayscale images, or 8-bit color images, into triangles.
//
// For grayscale images, a threshold is specified to determine insideness.
// For color images, an exact color is specified to determine insideness.
// Color images can be r8, rg16, rgb24, or rgba32. For a visual overview of
// the API and all the flags, see:
//
// http://github.prideout.net/marching-squares/
//
// The MIT License
// Copyright (c) 2015 Philip Rideout
#ifndef PAR_MSQUARES_H
#define PAR_MSQUARES_H
#ifdef __cplusplus
extern "C" {
#endif
#include <stdint.h>
// -----------------------------------------------------------------------------
// BEGIN PUBLIC API
// -----------------------------------------------------------------------------
typedef uint8_t par_byte;
typedef struct par_msquares_meshlist_s par_msquares_meshlist;
// Results of a marching squares operation. Triangles are counter-clockwise.
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;
// 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 in a single <path>,
// provided they each terminate with a "Z" and use 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;
// Reverses the "insideness" test.
#define PAR_MSQUARES_INVERT (1 << 0)
// Returns a meshlist with two meshes: one for the inside, one for the outside.
#define PAR_MSQUARES_DUAL (1 << 1)
// Requests that returned meshes have 3-tuple coordinates instead of 2-tuples.
// When using a color-based function, the Z coordinate represents the alpha
// value of the nearest pixel.
#define PAR_MSQUARES_HEIGHTS (1 << 2)
// Applies a step function to the Z coordinates. Requires HEIGHTS and DUAL.
#define PAR_MSQUARES_SNAP (1 << 3)
// Adds extrusion triangles to each mesh other than the lowest mesh. Requires
// the PAR_MSQUARES_HEIGHTS flag to be present.
#define PAR_MSQUARES_CONNECT (1 << 4)
// Enables quick & dirty (not best) simpification of the returned mesh.
#define PAR_MSQUARES_SIMPLIFY (1 << 5)
// Indicates that the "color" argument is ABGR instead of ARGB.
#define PAR_MSQUARES_SWIZZLE (1 << 6)
// Ensures there are no T-junction vertices. (par_msquares_color_multi only)
// Requires the PAR_MSQUARES_SIMPLIFY flag to be disabled.
#define PAR_MSQUARES_CLEAN (1 << 7)
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(par_byte const* data, int width,
int height, int cellsize, uint32_t color, int bpp, int flags);
par_msquares_mesh const* par_msquares_get_mesh(par_msquares_meshlist*, int n);
int par_msquares_get_count(par_msquares_meshlist*);
void par_msquares_free(par_msquares_meshlist*);
void par_msquares_free_boundary(par_msquares_boundary*);
typedef int (*par_msquares_inside_fn)(int, void*);
typedef float (*par_msquares_height_fn)(float, float, void*);
par_msquares_meshlist* par_msquares_function(int width, int height,
int cellsize, int flags, void* context, par_msquares_inside_fn insidefn,
par_msquares_height_fn heightfn);
par_msquares_meshlist* par_msquares_grayscale_multi(float const* data,
int width, int height, int cellsize, float const* thresholds,
int nthresholds, int flags);
par_msquares_meshlist* par_msquares_color_multi(par_byte const* data, int width,
int height, int cellsize, int bpp, int flags);
par_msquares_boundary* par_msquares_extract_boundary(par_msquares_mesh const* );
#ifndef PAR_PI
#define PAR_PI (3.14159265359)
#define PAR_MIN(a, b) (a > b ? b : a)
#define PAR_MAX(a, b) (a > b ? a : b)
#define PAR_CLAMP(v, lo, hi) PAR_MAX(lo, PAR_MIN(hi, v))
#define PAR_SWAP(T, A, B) { T tmp = B; B = A; A = tmp; }
#define PAR_SQR(a) ((a) * (a))
#endif
#ifndef PAR_MALLOC
#define PAR_MALLOC(T, N) ((T*) malloc(N * sizeof(T)))
#define PAR_CALLOC(T, N) ((T*) calloc(N * sizeof(T), 1))
#define PAR_REALLOC(T, BUF, N) ((T*) realloc(BUF, sizeof(T) * (N)))
#define PAR_FREE(BUF) free(BUF)
#endif
#ifdef __cplusplus
}
#endif
// -----------------------------------------------------------------------------
// END PUBLIC API
// -----------------------------------------------------------------------------
#ifdef PAR_MSQUARES_IMPLEMENTATION
#include <stdlib.h>
#include <assert.h>
#include <float.h>
#include <string.h>
typedef struct {
uint16_t* values;
size_t count;
size_t capacity;
} par__uint16list;
typedef struct {
float* points;
int npoints;
uint16_t* triangles;
int ntriangles;
int dim;
uint32_t color;
int nconntriangles;
uint16_t* conntri;
par__uint16list* tjunctions;
} par_msquares__mesh;
struct par_msquares_meshlist_s {
int nmeshes;
par_msquares__mesh** meshes;
};
static int** par_msquares_binary_point_table = 0;
static int** par_msquares_binary_triangle_table = 0;
static int* par_msquares_quaternary_triangle_table[64][4];
static int* par_msquares_quaternary_boundary_table[64][4];
static par_msquares_meshlist* par_msquares__merge(par_msquares_meshlist** lists,
int count, int snap);
static void par_init_tables()
{
char const* BINARY_TABLE =
"0"
"1017"
"1123"
"2023370"
"1756"
"2015560"
"2123756"
"3023035056"
"1345"
"4013034045057"
"2124451"
"3024045057"
"2734467"
"3013034046"
"3124146167"
"2024460";
char const* binary_token = BINARY_TABLE;
par_msquares_binary_point_table = PAR_CALLOC(int*, 16);
par_msquares_binary_triangle_table = PAR_CALLOC(int*, 16);
for (int i = 0; i < 16; i++) {
int ntris = *binary_token - '0';
binary_token++;
par_msquares_binary_triangle_table[i] =
PAR_CALLOC(int, (ntris + 1) * 3);
int* sqrtris = par_msquares_binary_triangle_table[i];
sqrtris[0] = ntris;
int mask = 0;
int* sqrpts = par_msquares_binary_point_table[i] = PAR_CALLOC(int, 7);
sqrpts[0] = 0;
for (int j = 0; j < ntris * 3; j++, binary_token++) {
int midp = *binary_token - '0';
int bit = 1 << midp;
if (!(mask & bit)) {
mask |= bit;
sqrpts[++sqrpts[0]] = midp;
}
sqrtris[j + 1] = midp;
}
}
char const* QUATERNARY_TABLE =
"2024046000"
"3346360301112300"
"3346360301112300"
"3346360301112300"
"3560502523013450"
"2015056212414500"
"4018087785756212313828348450"
"4018087785756212313828348450"
"3560502523013450"
"4018087785756212313828348450"
"2015056212414500"
"4018087785756212313828348450"
"3560502523013450"
"4018087785756212313828348450"
"4018087785756212313828348450"
"2015056212414500"
"3702724745001756"
"2018087212313828348452785756"
"4013034045057112301756"
"4013034045057112301756"
"2023037027347460"
"1701312414616700"
"2018087212313847857568348450"
"2018087212313847857568348450"
"4018087123138028348452785756"
"1701467161262363513450"
"2018087412313883484502785756"
"2018087212313828348452785756"
"4018087123138028348452785756"
"1701467161262363513450"
"2018087212313828348452785756"
"2018087412313883484502785756"
"3702724745001756"
"4013034045057112301756"
"2018087212313828348452785756"
"4013034045057112301756"
"4018087123138028348452785756"
"2018087412313883484502785756"
"1701467161262363513450"
"2018087212313828348452785756"
"2023037027347460"
"2018087212313847857568348450"
"1701312414616700"
"2018087212313847857568348450"
"4018087123138028348452785756"
"2018087212313828348452785756"
"1701467161262363513450"
"2018087412313883484502785756"
"3702724745001756"
"4013034045057112301756"
"4013034045057112301756"
"2018087212313828348452785756"
"4018087123138028348452785756"
"2018087412313883484502785756"
"2018087212313828348452785756"
"1701467161262363513450"
"4018087123138028348452785756"
"2018087212313828348452785756"
"2018087412313883484502785756"
"1701467161262363513450"
"2023037027347460"
"2018087212313847857568348450"
"2018087212313847857568348450"
"1701312414616700";
char const* quaternary_token = QUATERNARY_TABLE;
int* quaternary_values = PAR_CALLOC(int, strlen(QUATERNARY_TABLE));
int* vals = quaternary_values;
for (int i = 0; i < 64; i++) {
int ntris = *quaternary_token++ - '0';
*vals = ntris;
par_msquares_quaternary_triangle_table[i][0] = vals++;
for (int j = 0; j < ntris * 3; j++) {
int pt = *quaternary_token++ - '0';
assert(pt >= 0 && pt < 9);
*vals++ = pt;
}
ntris = *quaternary_token++ - '0';
*vals = ntris;
par_msquares_quaternary_triangle_table[i][1] = vals++;
for (int j = 0; j < ntris * 3; j++) {
int pt = *quaternary_token++ - '0';
assert(pt >= 0 && pt < 9);
*vals++ = pt;
}
ntris = *quaternary_token++ - '0';
*vals = ntris;
par_msquares_quaternary_triangle_table[i][2] = vals++;
for (int j = 0; j < ntris * 3; j++) {
int pt = *quaternary_token++ - '0';
assert(pt >= 0 && pt < 9);
*vals++ = pt;
}
ntris = *quaternary_token++ - '0';
*vals = ntris;
par_msquares_quaternary_triangle_table[i][3] = vals++;
for (int j = 0; j < ntris * 3; j++) {
int pt = *quaternary_token++ - '0';
assert(pt >= 0 && pt < 9);
*vals++ = pt;
}
}
assert(vals == quaternary_values + strlen(QUATERNARY_TABLE));
char const* QUATERNARY_EDGES =
"0000"
"11313100113131001131310013501530"
"115151002188523881258830218852388125883013501530"
"218852388125883011515100218852388125883013501530"
"218852388125883021885238812588301151510015700175"
"2188723881258832788521357131017521357131017513701730"
"11717100218872388127883021887238812788302388702588327885"
"1172713515302188725881027885218872388125883278852388702588327885"
"11727135153021887238812588327885218872588102788515700175"
"213571310175218872388125883278852135713101752388702588327885"
"21887258810278851172713515302188723881258832788513701730"
"21887238812788301171710021887238812788302388702588327885"
"21887238812588327885117271351530218872588102788515700175"
"213571310175213571310175218872388125883278852388702588327885"
"2188725881027885218872388125883278851172713515302388702588327885"
"21887238812588327885218872588102788511727135153013701730"
"2188723881278830218872388127883011717100";
quaternary_token = QUATERNARY_EDGES;
quaternary_values = PAR_CALLOC(int, strlen(QUATERNARY_EDGES));
vals = quaternary_values;
for (int i = 0; i < 64; i++) {
int nedges = *quaternary_token++ - '0';
*vals = nedges;
par_msquares_quaternary_boundary_table[i][0] = vals++;
for (int j = 0; j < nedges * 2; j++) {
int pt = *quaternary_token++ - '0';
assert(pt >= 0 && pt < 9);
*vals++ = pt;
}
nedges = *quaternary_token++ - '0';
*vals = nedges;
par_msquares_quaternary_boundary_table[i][1] = vals++;
for (int j = 0; j < nedges * 2; j++) {
int pt = *quaternary_token++ - '0';
assert(pt >= 0 && pt < 9);
*vals++ = pt;
}
nedges = *quaternary_token++ - '0';
*vals = nedges;
par_msquares_quaternary_boundary_table[i][2] = vals++;
for (int j = 0; j < nedges * 2; j++) {
int pt = *quaternary_token++ - '0';
assert(pt >= 0 && pt < 9);
*vals++ = pt;
}
nedges = *quaternary_token++ - '0';
*vals = nedges;
par_msquares_quaternary_boundary_table[i][3] = vals++;
for (int j = 0; j < nedges * 2; j++) {
int pt = *quaternary_token++ - '0';
assert(pt >= 0 && pt < 9);
*vals++ = pt;
}
}
assert(vals == quaternary_values + strlen(QUATERNARY_EDGES));
}
typedef struct {
float const* data;
float threshold;
float lower_bound;
float upper_bound;
int width;
int height;
} par_gray_context;
static int gray_inside(int location, void* contextptr)
{
par_gray_context* context = (par_gray_context*) contextptr;
return context->data[location] > context->threshold;
}
static int gray_multi_inside(int location, void* contextptr)
{
par_gray_context* context = (par_gray_context*) contextptr;
float val = context->data[location];
float upper = context->upper_bound;
float lower = context->lower_bound;
return val >= lower && val < upper;
}
static float gray_height(float x, float y, void* contextptr)
{
par_gray_context* context = (par_gray_context*) contextptr;
int i = PAR_CLAMP(context->width * x, 0, context->width - 1);
int j = PAR_CLAMP(context->height * y, 0, context->height - 1);
return context->data[i + j * context->width];
}
typedef struct {
par_byte const* data;
par_byte color[4];
int bpp;
int width;
int height;
} par_color_context;
static int color_inside(int location, void* contextptr)
{
par_color_context* context = (par_color_context*) contextptr;
par_byte const* data = context->data + location * context->bpp;
for (int i = 0; i < context->bpp; i++) {
if (data[i] != context->color[i]) {
return 0;
}
}
return 1;
}
static float color_height(float x, float y, void* contextptr)
{
par_color_context* context = (par_color_context*) contextptr;
assert(context->bpp == 4);
int i = PAR_CLAMP(context->width * x, 0, context->width - 1);
int j = PAR_CLAMP(context->height * y, 0, context->height - 1);
int k = i + j * context->width;
return context->data[k * 4 + 3] / 255.0;
}
par_msquares_meshlist* par_msquares_color(par_byte const* data, int width,
int height, int cellsize, uint32_t color, int bpp, int flags)
{
par_color_context context;
context.bpp = bpp;
if (flags & PAR_MSQUARES_SWIZZLE) {
context.color[0] = (color >> 0) & 0xff;
context.color[1] = (color >> 8) & 0xff;
context.color[2] = (color >> 16) & 0xff;
context.color[3] = (color >> 24) & 0xff;
} else {
context.color[0] = (color >> 16) & 0xff;
context.color[1] = (color >> 8) & 0xff;
context.color[2] = (color >> 0) & 0xff;
context.color[3] = (color >> 24) & 0xff;
}
context.data = data;
context.width = width;
context.height = height;
return par_msquares_function(
width, height, cellsize, flags, &context, color_inside, color_height);
}
par_msquares_meshlist* par_msquares_grayscale(float const* data, int width,
int height, int cellsize, float threshold, int flags)
{
par_gray_context context;
context.width = width;
context.height = height;
context.data = data;
context.threshold = threshold;
return par_msquares_function(
width, height, cellsize, flags, &context, gray_inside, gray_height);
}
par_msquares_meshlist* par_msquares_grayscale_multi(float const* data,
int width, int height, int cellsize, float const* thresholds,
int nthresholds, int flags)
{
par_msquares_meshlist* mlists[2];
mlists[0] = PAR_CALLOC(par_msquares_meshlist, 1);
int connect = flags & PAR_MSQUARES_CONNECT;
int snap = flags & PAR_MSQUARES_SNAP;
int heights = flags & PAR_MSQUARES_HEIGHTS;
if (!heights) {
snap = connect = 0;
}
flags &= ~PAR_MSQUARES_INVERT;
flags &= ~PAR_MSQUARES_DUAL;
flags &= ~PAR_MSQUARES_CONNECT;
flags &= ~PAR_MSQUARES_SNAP;
par_gray_context context;
context.width = width;
context.height = height;
context.data = data;
context.lower_bound = -FLT_MAX;
for (int i = 0; i <= nthresholds; i++) {
int mergeconf = i > 0 ? connect : 0;
if (i == nthresholds) {
context.upper_bound = FLT_MAX;
mergeconf |= snap;
} else {
context.upper_bound = thresholds[i];
}
mlists[1] = par_msquares_function(width, height, cellsize, flags,
&context, gray_multi_inside, gray_height);
mlists[0] = par_msquares__merge(mlists, 2, mergeconf);
context.lower_bound = context.upper_bound;
flags |= connect;
}
return mlists[0];
}
par_msquares_mesh const* par_msquares_get_mesh(
par_msquares_meshlist* mlist, int mindex)
{
assert(mlist && mindex < mlist->nmeshes);
return (par_msquares_mesh const*) mlist->meshes[mindex];
}
int par_msquares_get_count(par_msquares_meshlist* mlist)
{
assert(mlist);
return mlist->nmeshes;
}
void par_msquares_free(par_msquares_meshlist* mlist)
{
if (!mlist) {
return;
}
par_msquares__mesh** meshes = mlist->meshes;
for (int i = 0; i < mlist->nmeshes; i++) {
free(meshes[i]->points);
free(meshes[i]->triangles);
free(meshes[i]);
}
free(meshes);
free(mlist);
}
// Combine multiple meshlists by moving mesh pointers, and optionally apply
// a "snap" operation that assigns a single Z value across all verts in each
// mesh. The Z value determined by the mesh's position in the final mesh list.
static par_msquares_meshlist* par_msquares__merge(par_msquares_meshlist** lists,
int count, int snap)
{
par_msquares_meshlist* merged = PAR_CALLOC(par_msquares_meshlist, 1);
merged->nmeshes = 0;
for (int i = 0; i < count; i++) {
merged->nmeshes += lists[i]->nmeshes;
}
merged->meshes = PAR_CALLOC(par_msquares__mesh*, merged->nmeshes);
par_msquares__mesh** pmesh = merged->meshes;
for (int i = 0; i < count; i++) {
par_msquares_meshlist* meshlist = lists[i];
for (int j = 0; j < meshlist->nmeshes; j++) {
*pmesh++ = meshlist->meshes[j];
}
free(meshlist);
}
if (!snap) {
return merged;
}
pmesh = merged->meshes;
float zmin = FLT_MAX;
float zmax = -zmin;
for (int i = 0; i < merged->nmeshes; i++, pmesh++) {
float* pzed = (*pmesh)->points + 2;
for (int j = 0; j < (*pmesh)->npoints; j++, pzed += 3) {
zmin = PAR_MIN(*pzed, zmin);
zmax = PAR_MAX(*pzed, zmax);
}
}
float zextent = zmax - zmin;
pmesh = merged->meshes;
for (int i = 0; i < merged->nmeshes; i++, pmesh++) {
float* pzed = (*pmesh)->points + 2;
float zed = zmin + zextent * i / (merged->nmeshes - 1);
for (int j = 0; j < (*pmesh)->npoints; j++, pzed += 3) {
*pzed = zed;
}
}
if (!(snap & PAR_MSQUARES_CONNECT)) {
return merged;
}
for (int i = 1; i < merged->nmeshes; i++) {
par_msquares__mesh* mesh = merged->meshes[i];
// Find all extrusion points. This is tightly coupled to the
// tessellation code, which generates two "connector" triangles for each
// extruded edge. The first two verts of the second triangle are the
// verts that need to be displaced.
char* markers = PAR_CALLOC(char, mesh->npoints);
int tri = mesh->ntriangles - mesh->nconntriangles;
while (tri < mesh->ntriangles) {
markers[mesh->triangles[tri * 3 + 3]] = 1;
markers[mesh->triangles[tri * 3 + 4]] = 1;
tri += 2;
}
// Displace all extrusion points down to the previous level.
float zed = zmin + zextent * (i - 1) / (merged->nmeshes - 1);
float* pzed = mesh->points + 2;
for (int j = 0; j < mesh->npoints; j++, pzed += 3) {
if (markers[j]) {
*pzed = zed;
}
}
free(markers);
}
return merged;
}
static void par_remove_unreferenced_verts(par_msquares__mesh* mesh)
{
if (mesh->npoints == 0) {
return;
}
char* markers = PAR_CALLOC(char, mesh->npoints);
uint16_t const* ptris = mesh->triangles;
int newnpts = 0;
for (int i = 0; i < mesh->ntriangles * 3; i++, ptris++) {
if (!markers[*ptris]) {
newnpts++;
markers[*ptris] = 1;
}
}
float* newpts = PAR_CALLOC(float, newnpts * mesh->dim);
uint16_t* mapping = PAR_CALLOC(uint16_t, mesh->npoints);
float const* ppts = mesh->points;
float* pnewpts = newpts;
int j = 0;
if (mesh->dim == 3) {
for (int i = 0; i < mesh->npoints; i++, ppts += 3) {
if (markers[i]) {
*pnewpts++ = ppts[0];
*pnewpts++ = ppts[1];
*pnewpts++ = ppts[2];
mapping[i] = j++;
}
}
} else {
for (int i = 0; i < mesh->npoints; i++, ppts += 2) {
if (markers[i]) {
*pnewpts++ = ppts[0];
*pnewpts++ = ppts[1];
mapping[i] = j++;
}
}
}
free(mesh->points);
free(markers);
mesh->points = newpts;
mesh->npoints = newnpts;
for (int i = 0; i < mesh->ntriangles * 3; i++) {
mesh->triangles[i] = mapping[mesh->triangles[i]];
}
free(mapping);
}
par_msquares_meshlist* par_msquares_function(int width, int height,
int cellsize, int flags, void* context, par_msquares_inside_fn insidefn,
par_msquares_height_fn heightfn)
{
assert(width > 0 && width % cellsize == 0);
assert(height > 0 && height % cellsize == 0);
if (flags & PAR_MSQUARES_DUAL) {
int connect = flags & PAR_MSQUARES_CONNECT;
int snap = flags & PAR_MSQUARES_SNAP;
int heights = flags & PAR_MSQUARES_HEIGHTS;
if (!heights) {
snap = connect = 0;
}
flags ^= PAR_MSQUARES_INVERT;
flags &= ~PAR_MSQUARES_DUAL;
flags &= ~PAR_MSQUARES_CONNECT;
par_msquares_meshlist* m[2];
m[0] = par_msquares_function(width, height, cellsize, flags,
context, insidefn, heightfn);
flags ^= PAR_MSQUARES_INVERT;
if (connect) {
flags |= PAR_MSQUARES_CONNECT;
}
m[1] = par_msquares_function(width, height, cellsize, flags,
context, insidefn, heightfn);
return par_msquares__merge(m, 2, snap | connect);
}
int invert = flags & PAR_MSQUARES_INVERT;
// Create the two code tables if we haven't already. These are tables of
// fixed constants, so it's embarassing that we use dynamic memory
// allocation for them. However it's easy and it's one-time-only.
if (!par_msquares_binary_point_table) {
par_init_tables();
}
// Allocate the meshlist and the first mesh.
par_msquares_meshlist* mlist = PAR_CALLOC(par_msquares_meshlist, 1);
mlist->nmeshes = 1;
mlist->meshes = PAR_CALLOC(par_msquares__mesh*, 1);
mlist->meshes[0] = PAR_CALLOC(par_msquares__mesh, 1);
par_msquares__mesh* mesh = mlist->meshes[0];
mesh->dim = (flags & PAR_MSQUARES_HEIGHTS) ? 3 : 2;
int ncols = width / cellsize;
int nrows = height / cellsize;
// Worst case is four triangles and six verts per cell, so allocate that
// much.
int maxtris = ncols * nrows * 4;
int maxpts = ncols * nrows * 6;
int maxedges = ncols * nrows * 2;
// However, if we include extrusion triangles for boundary edges,
// we need space for another 4 triangles and 4 points per cell.
uint16_t* conntris = 0;
int nconntris = 0;
uint16_t* edgemap = 0;
if (flags & PAR_MSQUARES_CONNECT) {
conntris = PAR_CALLOC(uint16_t, maxedges * 6);
maxtris += maxedges * 2;
maxpts += maxedges * 2;
edgemap = PAR_CALLOC(uint16_t, maxpts);
for (int i = 0; i < maxpts; i++) {
edgemap[i] = 0xffff;
}
}
uint16_t* tris = PAR_CALLOC(uint16_t, maxtris * 3);
int ntris = 0;
float* pts = PAR_CALLOC(float, maxpts * mesh->dim);
int npts = 0;
// The "verts" x/y/z arrays are the 4 corners and 4 midpoints around the
// square, in counter-clockwise order. The origin of "triangle space" is at
// the lower-left, although we expect the image data to be in raster order
// (starts at top-left).
float vertsx[8], vertsy[8];
float normalization = 1.0f / PAR_MAX(width, height);
float normalized_cellsize = cellsize * normalization;
int maxrow = (height - 1) * width;
uint16_t* ptris = tris;
uint16_t* pconntris = conntris;
float* ppts = pts;
uint8_t* prevrowmasks = PAR_CALLOC(uint8_t, ncols);
int* prevrowinds = PAR_CALLOC(int, ncols * 3);
// If simplification is enabled, we need to track all 'F' cells and their
// repsective triangle indices.
uint8_t* simplification_codes = 0;
uint16_t* simplification_tris = 0;
uint8_t* simplification_ntris = 0;
if (flags & PAR_MSQUARES_SIMPLIFY) {
simplification_codes = PAR_CALLOC(uint8_t, nrows * ncols);
simplification_tris = PAR_CALLOC(uint16_t, nrows * ncols);
simplification_ntris = PAR_CALLOC(uint8_t, nrows * ncols);
}
// Do the march!
for (int row = 0; row < nrows; row++) {
vertsx[0] = vertsx[6] = vertsx[7] = 0;
vertsx[1] = vertsx[5] = 0.5 * normalized_cellsize;
vertsx[2] = vertsx[3] = vertsx[4] = normalized_cellsize;
vertsy[0] = vertsy[1] = vertsy[2] = normalized_cellsize * (row + 1);
vertsy[4] = vertsy[5] = vertsy[6] = normalized_cellsize * row;
vertsy[3] = vertsy[7] = normalized_cellsize * (row + 0.5);
int northi = row * cellsize * width;
int southi = PAR_MIN(northi + cellsize * width, maxrow);
int northwest = invert ^ insidefn(northi, context);
int southwest = invert ^ insidefn(southi, context);
int previnds[8] = {0};
uint8_t prevmask = 0;
for (int col = 0; col < ncols; col++) {
northi += cellsize;
southi += cellsize;
if (col == ncols - 1) {
northi--;
southi--;
}
int northeast = invert ^ insidefn(northi, context);
int southeast = invert ^ insidefn(southi, context);
int code = southwest | (southeast << 1) | (northwest << 2) |
(northeast << 3);
int const* pointspec = par_msquares_binary_point_table[code];
int ptspeclength = *pointspec++;
int currinds[8] = {0};
uint8_t mask = 0;
uint8_t prevrowmask = prevrowmasks[col];
while (ptspeclength--) {
int midp = *pointspec++;
int bit = 1 << midp;
mask |= bit;
// The following six conditionals perform welding to reduce the
// number of vertices. The first three perform welding with the
// cell to the west; the latter three perform welding with the
// cell to the north.
if (bit == 1 && (prevmask & 4)) {
currinds[midp] = previnds[2];
continue;
}
if (bit == 128 && (prevmask & 8)) {
currinds[midp] = previnds[3];
continue;
}
if (bit == 64 && (prevmask & 16)) {
currinds[midp] = previnds[4];
continue;
}
if (bit == 16 && (prevrowmask & 4)) {
currinds[midp] = prevrowinds[col * 3 + 2];
continue;
}
if (bit == 32 && (prevrowmask & 2)) {
currinds[midp] = prevrowinds[col * 3 + 1];
continue;
}
if (bit == 64 && (prevrowmask & 1)) {
currinds[midp] = prevrowinds[col * 3 + 0];
continue;
}
ppts[0] = vertsx[midp];
ppts[1] = vertsy[midp];
// Adjust the midpoints to a more exact crossing point.
if (midp == 1) {
int begin = southi - cellsize / 2;
int previous = 0;
for (int i = 0; i < cellsize; i++) {
int offset = begin + i / 2 * ((i % 2) ? -1 : 1);
int inside = insidefn(offset, context);
if (i > 0 && inside != previous) {
ppts[0] = normalization *
(col * cellsize + offset - southi + cellsize);
break;
}
previous = inside;
}
} else if (midp == 5) {
int begin = northi - cellsize / 2;
int previous = 0;
for (int i = 0; i < cellsize; i++) {
int offset = begin + i / 2 * ((i % 2) ? -1 : 1);
int inside = insidefn(offset, context);
if (i > 0 && inside != previous) {
ppts[0] = normalization *
(col * cellsize + offset - northi + cellsize);
break;
}
previous = inside;
}
} else if (midp == 3) {
int begin = northi + width * cellsize / 2;
int previous = 0;
for (int i = 0; i < cellsize; i++) {
int offset = begin +
width * (i / 2 * ((i % 2) ? -1 : 1));
int inside = insidefn(offset, context);
if (i > 0 && inside != previous) {
ppts[1] = normalization *
(row * cellsize +
(offset - northi) / (float) width);
break;
}
previous = inside;
}
} else if (midp == 7) {
int begin = northi + width * cellsize / 2 - cellsize;
int previous = 0;
for (int i = 0; i < cellsize; i++) {
int offset = begin +
width * (i / 2 * ((i % 2) ? -1 : 1));
int inside = insidefn(offset, context);
if (i > 0 && inside != previous) {
ppts[1] = normalization *
(row * cellsize +
(offset - northi - cellsize) / (float) width);
break;
}
previous = inside;
}
}
if (mesh->dim == 3) {
ppts[2] = heightfn(ppts[0], ppts[1], context);
}
ppts += mesh->dim;
currinds[midp] = npts++;
}
int const* trianglespec = par_msquares_binary_triangle_table[code];
int trispeclength = *trianglespec++;
if (flags & PAR_MSQUARES_SIMPLIFY) {
simplification_codes[ncols * row + col] = code;
simplification_tris[ncols * row + col] = ntris;
simplification_ntris[ncols * row + col] = trispeclength;
}
// Add triangles.
while (trispeclength--) {
int a = *trianglespec++;
int b = *trianglespec++;
int c = *trianglespec++;
*ptris++ = currinds[c];
*ptris++ = currinds[b];
*ptris++ = currinds[a];
ntris++;
}
// Create two extrusion triangles for each boundary edge.
if (flags & PAR_MSQUARES_CONNECT) {
trianglespec = par_msquares_binary_triangle_table[code];
trispeclength = *trianglespec++;
while (trispeclength--) {
int a = *trianglespec++;
int b = *trianglespec++;
int c = *trianglespec++;
int i = currinds[a];
int j = currinds[b];
int k = currinds[c];
int u = 0, v = 0, w = 0;
if ((a % 2) && (b % 2)) {
u = v = 1;
} else if ((a % 2) && (c % 2)) {
u = w = 1;
} else if ((b % 2) && (c % 2)) {
v = w = 1;
} else {
continue;
}
if (u && edgemap[i] == 0xffff) {
for (int d = 0; d < mesh->dim; d++) {
*ppts++ = pts[i * mesh->dim + d];
}
edgemap[i] = npts++;
}
if (v && edgemap[j] == 0xffff) {
for (int d = 0; d < mesh->dim; d++) {
*ppts++ = pts[j * mesh->dim + d];
}
edgemap[j] = npts++;
}
if (w && edgemap[k] == 0xffff) {
for (int d = 0; d < mesh->dim; d++) {
*ppts++ = pts[k * mesh->dim + d];
}
edgemap[k] = npts++;
}
if ((a % 2) && (b % 2)) {
*pconntris++ = i;
*pconntris++ = j;
*pconntris++ = edgemap[j];
*pconntris++ = edgemap[j];
*pconntris++ = edgemap[i];
*pconntris++ = i;
} else if ((a % 2) && (c % 2)) {
*pconntris++ = edgemap[k];
*pconntris++ = k;
*pconntris++ = i;
*pconntris++ = edgemap[i];
*pconntris++ = edgemap[k];
*pconntris++ = i;
} else if ((b % 2) && (c % 2)) {
*pconntris++ = j;
*pconntris++ = k;
*pconntris++ = edgemap[k];
*pconntris++ = edgemap[k];
*pconntris++ = edgemap[j];
*pconntris++ = j;
}
nconntris += 2;
}
}
// Prepare for the next cell.
prevrowmasks[col] = mask;
prevrowinds[col * 3 + 0] = currinds[0];
prevrowinds[col * 3 + 1] = currinds[1];
prevrowinds[col * 3 + 2] = currinds[2];
prevmask = mask;
northwest = northeast;
southwest = southeast;
for (int i = 0; i < 8; i++) {
previnds[i] = currinds[i];
vertsx[i] += normalized_cellsize;
}
}
}
free(edgemap);
free(prevrowmasks);
free(prevrowinds);
// Perform quick-n-dirty simplification by iterating two rows at a time.
// In no way does this create the simplest possible mesh, but at least it's
// fast and easy.
if (flags & PAR_MSQUARES_SIMPLIFY) {
int in_run = 0, start_run;
// First figure out how many triangles we can eliminate.
int neliminated_triangles = 0;
for (int row = 0; row < nrows - 1; row += 2) {
for (int col = 0; col < ncols; col++) {
int a = simplification_codes[ncols * row + col] == 0xf;
int b = simplification_codes[ncols * row + col + ncols] == 0xf;
if (a && b) {
if (!in_run) {
in_run = 1;
start_run = col;
}
continue;
}
if (in_run) {
in_run = 0;
int run_width = col - start_run;
neliminated_triangles += run_width * 4 - 2;
}
}
if (in_run) {
in_run = 0;
int run_width = ncols - start_run;
neliminated_triangles += run_width * 4 - 2;
}
}
// Build a new index array cell-by-cell. If any given cell is 'F' and
// its neighbor to the south is also 'F', then it's part of a run.
int nnewtris = ntris + nconntris - neliminated_triangles;
uint16_t* newtris = PAR_CALLOC(uint16_t, nnewtris * 3);
uint16_t* pnewtris = newtris;
in_run = 0;
for (int row = 0; row < nrows - 1; row += 2) {
for (int col = 0; col < ncols; col++) {
int cell = ncols * row + col;
int south = cell + ncols;
int a = simplification_codes[cell] == 0xf;
int b = simplification_codes[south] == 0xf;
if (a && b) {
if (!in_run) {
in_run = 1;
start_run = col;
}
continue;
}
if (in_run) {
in_run = 0;
int nw_cell = ncols * row + start_run;
int ne_cell = ncols * row + col - 1;
int sw_cell = nw_cell + ncols;
int se_cell = ne_cell + ncols;
int nw_tri = simplification_tris[nw_cell];
int ne_tri = simplification_tris[ne_cell];
int sw_tri = simplification_tris[sw_cell];
int se_tri = simplification_tris[se_cell];
int nw_corner = nw_tri * 3 + 4;
int ne_corner = ne_tri * 3 + 0;
int sw_corner = sw_tri * 3 + 2;
int se_corner = se_tri * 3 + 1;
*pnewtris++ = tris[se_corner];
*pnewtris++ = tris[sw_corner];
*pnewtris++ = tris[nw_corner];
*pnewtris++ = tris[nw_corner];
*pnewtris++ = tris[ne_corner];
*pnewtris++ = tris[se_corner];
}
int ncelltris = simplification_ntris[cell];
int celltri = simplification_tris[cell];
for (int t = 0; t < ncelltris; t++, celltri++) {
*pnewtris++ = tris[celltri * 3];
*pnewtris++ = tris[celltri * 3 + 1];
*pnewtris++ = tris[celltri * 3 + 2];
}
ncelltris = simplification_ntris[south];
celltri = simplification_tris[south];
for (int t = 0; t < ncelltris; t++, celltri++) {
*pnewtris++ = tris[celltri * 3];
*pnewtris++ = tris[celltri * 3 + 1];
*pnewtris++ = tris[celltri * 3 + 2];
}
}
if (in_run) {
in_run = 0;
int nw_cell = ncols * row + start_run;
int ne_cell = ncols * row + ncols - 1;
int sw_cell = nw_cell + ncols;
int se_cell = ne_cell + ncols;
int nw_tri = simplification_tris[nw_cell];
int ne_tri = simplification_tris[ne_cell];
int sw_tri = simplification_tris[sw_cell];
int se_tri = simplification_tris[se_cell];
int nw_corner = nw_tri * 3 + 4;
int ne_corner = ne_tri * 3 + 0;
int sw_corner = sw_tri * 3 + 2;
int se_corner = se_tri * 3 + 1;
*pnewtris++ = tris[se_corner];
*pnewtris++ = tris[sw_corner];
*pnewtris++ = tris[nw_corner];
*pnewtris++ = tris[nw_corner];
*pnewtris++ = tris[ne_corner];
*pnewtris++ = tris[se_corner];
}
}
ptris = pnewtris;
ntris -= neliminated_triangles;
free(tris);
tris = newtris;
free(simplification_codes);
free(simplification_tris);
free(simplification_ntris);
// Remove unreferenced points.
char* markers = PAR_CALLOC(char, npts);
ptris = tris;
int newnpts = 0;
for (int i = 0; i < ntris * 3; i++, ptris++) {
if (!markers[*ptris]) {
newnpts++;
markers[*ptris] = 1;
}
}
for (int i = 0; i < nconntris * 3; i++) {
if (!markers[conntris[i]]) {
newnpts++;
markers[conntris[i]] = 1;
}
}
float* newpts = PAR_CALLOC(float, newnpts * mesh->dim);
uint16_t* mapping = PAR_CALLOC(uint16_t, npts);
ppts = pts;
float* pnewpts = newpts;
int j = 0;
if (mesh->dim == 3) {
for (int i = 0; i < npts; i++, ppts += 3) {
if (markers[i]) {
*pnewpts++ = ppts[0];
*pnewpts++ = ppts[1];
*pnewpts++ = ppts[2];
mapping[i] = j++;
}
}
} else {
for (int i = 0; i < npts; i++, ppts += 2) {
if (markers[i]) {
*pnewpts++ = ppts[0];
*pnewpts++ = ppts[1];
mapping[i] = j++;
}
}
}
free(pts);
free(markers);
pts = newpts;
npts = newnpts;
for (int i = 0; i < ntris * 3; i++) {
tris[i] = mapping[tris[i]];
}
for (int i = 0; i < nconntris * 3; i++) {
conntris[i] = mapping[conntris[i]];
}
free(mapping);
}
// Append all extrusion triangles to the main triangle array.
// We need them to be last so that they form a contiguous sequence.
pconntris = conntris;
for (int i = 0; i < nconntris; i++) {
*ptris++ = *pconntris++;
*ptris++ = *pconntris++;
*ptris++ = *pconntris++;
ntris++;
}
free(conntris);
// Final cleanup and return.
assert(npts <= maxpts);
assert(ntris <= maxtris);
mesh->npoints = npts;
mesh->points = pts;
mesh->ntriangles = ntris;
mesh->triangles = tris;
mesh->nconntriangles = nconntris;
return mlist;
}
typedef struct {
uint16_t outera;
uint16_t outerb;
uint16_t innera;
uint16_t innerb;
char i;
char j;
par_msquares__mesh* mesh;
int mesh_index;
} par_connector;
static par_connector* par_conn_find(par_connector* conns, int nconns,
char i, char j)
{
for (int c = 0; c < nconns; c++) {
if (conns[c].i == i && conns[c].j == j) {
return conns + c;
}
}
return 0;
}
static int par_msquares_cmp(const void *a, const void *b)
{
uint32_t arg1 = *((uint32_t const*) a);
uint32_t arg2 = *((uint32_t const*) b);
if (arg1 < arg2) return -1;
if (arg1 > arg2) return 1;
return 0;
}
typedef int (*par_msquares_code_fn)(int, int, int, int, void*);
static int par_msquares_multi_code(int sw, int se, int ne, int nw)
{
int code[4];
int ncols = 0;
code[0] = ncols++;
if (se == sw) {
code[1] = code[0];
} else {
code[1] = ncols++;
}
if (ne == se) {
code[2] = code[1];
} else if (ne == sw) {
code[2] = code[0];
} else {
code[2] = ncols++;
}
if (nw == ne) {
code[3] = code[2];
} else if (nw == se) {
code[3] = code[1];
} else if (nw == sw) {
code[3] = code[0];
} else {
code[3] = ncols++;
}
return code[0] | (code[1] << 2) | (code[2] << 4) | (code[3] << 6);
}
static uint32_t par_msquares_argb(par_byte const* pdata, int bpp)
{
uint32_t color = 0;
if (bpp == 4) {
color |= pdata[2];
color |= pdata[1] << 8;
color |= pdata[0] << 16;
color |= pdata[3] << 24;
return color;
}
for (int j = 0; j < bpp; j++) {
color <<= 8;
color |= pdata[j];
}
return color;
}
// Merge connective triangles into the primary triangle list.
static void par_msquares__finalize(par_msquares_meshlist* mlist)
{
if (mlist->nmeshes < 2 || mlist->meshes[1]->nconntriangles == 0) {
return;
}
for (int m = 1; m < mlist->nmeshes; m++) {
par_msquares__mesh* mesh = mlist->meshes[m];
int ntris = mesh->ntriangles + mesh->nconntriangles;
uint16_t* triangles = PAR_CALLOC(uint16_t, ntris * 3);
uint16_t* dst = triangles;
uint16_t const* src = mesh->triangles;
for (int t = 0; t < mesh->ntriangles; t++) {
*dst++ = *src++;
*dst++ = *src++;
*dst++ = *src++;
}
src = mesh->conntri;
for (int t = 0; t < mesh->nconntriangles; t++) {
*dst++ = *src++;
*dst++ = *src++;
*dst++ = *src++;
}
free(mesh->triangles);
free(mesh->conntri);
mesh->triangles = triangles;
mesh->ntriangles = ntris;
mesh->conntri = 0;
mesh->nconntriangles = 0;
}
}
static par__uint16list* par__uint16list_create()
{
par__uint16list* list = PAR_CALLOC(par__uint16list, 1);
list->count = 0;
list->capacity = 32;
list->values = PAR_CALLOC(uint16_t, list->capacity);
return list;
}
static void par__uint16list_add3(par__uint16list* list,
uint16_t a, uint16_t b, uint16_t c)
{
if (list->count + 3 > list->capacity) {
list->capacity *= 2;
list->values = PAR_REALLOC(uint16_t, list->values, list->capacity);
}
list->values[list->count++] = a;
list->values[list->count++] = b;
list->values[list->count++] = c;
}
static void par__uint16list_free(par__uint16list* list)
{
if (list) {
PAR_FREE(list->values);
PAR_FREE(list);
}
}
static void par_msquares__repair_tjunctions(par_msquares_meshlist* mlist)
{
for (int m = 0; m < mlist->nmeshes; m++) {
par_msquares__mesh* mesh = mlist->meshes[m];
par__uint16list* tjunctions = mesh->tjunctions;
int njunctions = (int) tjunctions->count / 3;
if (njunctions == 0) {
continue;
}
int ntriangles = mesh->ntriangles + njunctions;
mesh->triangles = PAR_REALLOC(uint16_t, mesh->triangles,
ntriangles * 3);
uint16_t const* jun = tjunctions->values;
uint16_t* new_triangles = mesh->triangles + mesh->ntriangles * 3;
int ncreated = 0;
for (int j = 0; j < njunctions; j++, jun += 3) {
uint16_t* tri = mesh->triangles;
int t;
for (t = 0; t < mesh->ntriangles; t++, tri += 3) {
int i = -1;
if (tri[0] == jun[0] && tri[1] == jun[1]) {
i = 0;
} else if (tri[1] == jun[0] && tri[2] == jun[1]) {
i = 1;
} else if (tri[2] == jun[0] && tri[0] == jun[1]) {
i = 2;
} else {
continue;
}
new_triangles[0] = tri[(i + 0) % 3];
new_triangles[1] = jun[2];
new_triangles[2] = tri[(i + 2) % 3];
tri[(i + 0) % 3] = jun[2];
new_triangles += 3;
ncreated++;
break;
}
// TODO: Need to investigate the "msquares_multi_diagram.obj" test.
assert(t != mesh->ntriangles &&
"Error with T-Junction repair; please disable the CLEAN flag.");
}
mesh->ntriangles += ncreated;
}
}
par_msquares_meshlist* par_msquares_color_multi(par_byte const* data, int width,
int height, int cellsize, int bpp, int flags)
{
if (!par_msquares_binary_point_table) {
par_init_tables();
}
const int ncols = width / cellsize;
const int nrows = height / cellsize;
const int maxrow = (height - 1) * width;
const int ncells = ncols * nrows;
const int dim = (flags & PAR_MSQUARES_HEIGHTS) ? 3 : 2;
const int west_to_east[9] = { 2, -1, -1, -1, -1, -1, 4, 3, -1 };
const int north_to_south[9] = { -1, -1, -1, -1, 2, 1, 0, -1, -1 };
assert(!(flags & PAR_MSQUARES_HEIGHTS) || bpp == 4);
assert(bpp > 0 && bpp <= 4 && "Bytes per pixel must be 1, 2, 3, or 4.");
assert(!(flags & PAR_MSQUARES_CLEAN) || !(flags & PAR_MSQUARES_SIMPLIFY));
assert(!(flags & PAR_MSQUARES_SNAP) &&
"SNAP is not supported with color_multi");
assert(!(flags & PAR_MSQUARES_INVERT) &&
"INVERT is not supported with color_multi");
assert(!(flags & PAR_MSQUARES_DUAL) &&
"DUAL is not supported with color_multi");
// Find all unique colors and ensure there are no more than 256 colors.
uint32_t colors[256];
int ncolors = 0;
par_byte const* pdata = data;
for (int i = 0; i < width * height; i++, pdata += bpp) {
uint32_t color = par_msquares_argb(pdata, bpp);
if (0 == bsearch(&color, colors, ncolors, 4, par_msquares_cmp)) {
assert(ncolors < 256);
colors[ncolors++] = color;
qsort(colors, ncolors, sizeof(uint32_t), par_msquares_cmp);
}
}
// Convert the color image to grayscale using the mapping table.
par_byte* pixels = PAR_CALLOC(par_byte, width * height);
pdata = data;
for (int i = 0; i < width * height; i++, pdata += bpp) {
uint32_t color = par_msquares_argb(pdata, bpp);
void* result = bsearch(&color, colors, ncolors, 4, par_msquares_cmp);
pixels[i] = (uint32_t*) result - &colors[0];
}
// Allocate 1 mesh for each color.
par_msquares_meshlist* mlist = PAR_CALLOC(par_msquares_meshlist, 1);
mlist->nmeshes = ncolors;
mlist->meshes = PAR_CALLOC(par_msquares__mesh*, ncolors);
par_msquares__mesh* mesh;
int maxtris_per_cell = 6;
int maxpts_per_cell = 9;
if (flags & PAR_MSQUARES_CONNECT) {
maxpts_per_cell += 6;
}
for (int i = 0; i < ncolors; i++) {
mesh = mlist->meshes[i] = PAR_CALLOC(par_msquares__mesh, 1);
mesh->color = colors[i];
mesh->points = PAR_CALLOC(float, ncells * maxpts_per_cell * dim);
mesh->triangles = PAR_CALLOC(uint16_t, ncells * maxtris_per_cell * 3);
mesh->dim = dim;
mesh->tjunctions = par__uint16list_create();
if (flags & PAR_MSQUARES_CONNECT) {
mesh->conntri = PAR_CALLOC(uint16_t, ncells * 8 * 3);
}
}
// The "verts" x/y/z arrays are the 4 corners and 4 midpoints around the
// square, in counter-clockwise order, starting at the lower-left. The
// ninth vert is the center point.
float vertsx[9], vertsy[9];
float normalization = 1.0f / PAR_MAX(width, height);
float normalized_cellsize = cellsize * normalization;
uint8_t cella[256];
uint8_t cellb[256];
uint8_t* currcell = cella;
uint8_t* prevcell = cellb;
uint16_t inds0[256 * 9];
uint16_t inds1[256 * 9];
uint16_t* currinds = inds0;
uint16_t* previnds = inds1;
uint16_t* rowindsa = PAR_CALLOC(uint16_t, ncols * 3 * 256);
uint8_t* rowcellsa = PAR_CALLOC(uint8_t, ncols * 256);
uint16_t* rowindsb = PAR_CALLOC(uint16_t, ncols * 3 * 256);
uint8_t* rowcellsb = PAR_CALLOC(uint8_t, ncols * 256);
uint16_t* prevrowinds = rowindsa;
uint16_t* currrowinds = rowindsb;
uint8_t* prevrowcells = rowcellsa;
uint8_t* currrowcells = rowcellsb;
uint32_t* simplification_words = 0;
if (flags & PAR_MSQUARES_SIMPLIFY) {
simplification_words = PAR_CALLOC(uint32_t, 2 * nrows * ncols);
}
// Do the march!
for (int row = 0; row < nrows; row++) {
vertsx[0] = vertsx[6] = vertsx[7] = 0;
vertsx[1] = vertsx[5] = vertsx[8] = 0.5 * normalized_cellsize;
vertsx[2] = vertsx[3] = vertsx[4] = normalized_cellsize;
vertsy[0] = vertsy[1] = vertsy[2] = normalized_cellsize * (row + 1);
vertsy[4] = vertsy[5] = vertsy[6] = normalized_cellsize * row;
vertsy[3] = vertsy[7] = vertsy[8] = normalized_cellsize * (row + 0.5);
int northi = row * cellsize * width;
int southi = PAR_MIN(northi + cellsize * width, maxrow);
int nwval = pixels[northi];
int swval = pixels[southi];
memset(currrowcells, 0, ncols * 256);
for (int col = 0; col < ncols; col++) {
northi += cellsize;
southi += cellsize;
if (col == ncols - 1) {
northi--;
southi--;
}
// Obtain 8-bit code and grab the four corresponding triangle lists.
int neval = pixels[northi];
int seval = pixels[southi];
int code = par_msquares_multi_code(swval, seval, neval, nwval) >> 2;
int const* trispecs[4] = {
par_msquares_quaternary_triangle_table[code][0],
par_msquares_quaternary_triangle_table[code][1],
par_msquares_quaternary_triangle_table[code][2],
par_msquares_quaternary_triangle_table[code][3]
};
int ntris[4] = {
*trispecs[0]++,
*trispecs[1]++,
*trispecs[2]++,
*trispecs[3]++
};
int const* edgespecs[4] = {
par_msquares_quaternary_boundary_table[code][0],
par_msquares_quaternary_boundary_table[code][1],
par_msquares_quaternary_boundary_table[code][2],
par_msquares_quaternary_boundary_table[code][3]
};
int nedges[4] = {
*edgespecs[0]++,
*edgespecs[1]++,
*edgespecs[2]++,
*edgespecs[3]++
};
int vals[4] = { swval, seval, neval, nwval };
// Gather topology information.
par_connector edges[16];
int ncedges = 0;
for (int c = 0; c < 4; c++) {
int color = vals[c];
par_msquares__mesh* mesh = mlist->meshes[color];
par_connector edge;
for (int e = 0; e < nedges[c]; e++) {
char previndex = edgespecs[c][e * 2];
char currindex = edgespecs[c][e * 2 + 1];
edge.i = previndex;
edge.j = currindex;
edge.mesh_index = color;
edge.mesh = mesh;
edges[ncedges++] = edge;
}
}
assert(ncedges < 16);
// Push triangles and points into the four affected meshes.
for (int m = 0; m < ncolors; m++) {
currcell[m] = 0;
}
uint32_t colors = 0;
uint32_t counts = 0;
uint16_t* conntris_start[4];
for (int c = 0; c < 4; c++) {
int color = vals[c];
colors |= color << (8 * c);
counts |= ntris[c] << (8 * c);
par_msquares__mesh* mesh = mlist->meshes[color];
float height = (mesh->color >> 24) / 255.0;
conntris_start[c] = mesh->conntri + mesh->nconntriangles * 3;
int usedpts[9] = {0};
uint16_t* pcurrinds = currinds + 9 * color;
uint16_t const* pprevinds = previnds + 9 * color;
uint16_t const* pprevrowinds =
prevrowinds + ncols * 3 * color + col * 3;
uint8_t prevrowcell = prevrowcells[color * ncols + col];
float* pdst = mesh->points + mesh->npoints * mesh->dim;
int previndex, prevflag;
for (int t = 0; t < ntris[c] * 3; t++) {
uint16_t index = trispecs[c][t];
if (usedpts[index]) {
continue;
}
usedpts[index] = 1;
if (index < 8) {
currcell[color] |= 1 << index;
}
// Vertical welding.
previndex = north_to_south[index];
prevflag = (previndex > -1) ? (1 << previndex) : 0;
if (row > 0 && (prevrowcell & prevflag)) {
pcurrinds[index] = pprevrowinds[previndex];
continue;
}
// Horizontal welding.
previndex = west_to_east[index];
prevflag = (previndex > -1) ? (1 << previndex) : 0;
if (col > 0 && (prevcell[color] & prevflag)) {
pcurrinds[index] = pprevinds[previndex];
continue;
}
// Insert brand new point.
float* vertex = pdst;
*pdst++ = vertsx[index];
*pdst++ = 1 - vertsy[index];
if (mesh->dim == 3) {
*pdst++ = height;
}
pcurrinds[index] = mesh->npoints++;
// If this is a midpoint, nudge it to the intersection.
if (index == 1) {
int begin = southi - cellsize;
for (int i = 1; i < cellsize + 1; i++) {
int val = pixels[begin + i];
if (val != pixels[begin]) {
vertex[0] = vertsx[0] + normalized_cellsize *
(float) i / cellsize;
break;
}
}
} else if (index == 3) {
int begin = northi;
for (int i = 1; i < cellsize + 1; i++) {
int val = pixels[begin + i * width];
if (val != pixels[begin]) {
vertex[1] = (1 - vertsy[4]) -
normalized_cellsize * (float) i / cellsize;
break;
}
}
}
}
// Look for T junctions and note them for later repairs.
uint8_t prc = prevrowcell;
if (usedpts[4] && !usedpts[5] && usedpts[6] && (prc & 2)) {
// Above cell had a middle vert, current cell straddles it.
par__uint16list_add3(mesh->tjunctions,
pcurrinds[4], pcurrinds[6], pprevrowinds[1]);
} else if ((prc & 1) && !(prc & 2) && (prc & 4) && usedpts[5]) {
// Current cell has a middle vert, above cell straddles it.
par__uint16list_add3(mesh->tjunctions,
pprevrowinds[0], pprevrowinds[2], pcurrinds[5]);
}
uint8_t pcc = col > 0 ? prevcell[color] : 0;
if (usedpts[0] && !usedpts[7] && usedpts[6] && (pcc & 8)) {
// Left cell had a middle vert, current cell straddles it.
par__uint16list_add3(mesh->tjunctions,
pcurrinds[6], pcurrinds[0], pprevinds[3]);
}
if ((pcc & 4) && !(pcc & 8) && (pcc & 16) && usedpts[7]) {
// Current cell has a middle vert, left cell straddles it.
par__uint16list_add3(mesh->tjunctions,
pprevinds[2], pprevinds[4], pcurrinds[7]);
}
// Stamp out the cell's triangle indices for this color.
uint16_t* tdst = mesh->triangles + mesh->ntriangles * 3;
mesh->ntriangles += ntris[c];
for (int t = 0; t < ntris[c] * 3; t++) {
uint16_t index = trispecs[c][t];
*tdst++ = pcurrinds[index];
}
// Add extrusion points and connective triangles if requested.
if (!(flags & PAR_MSQUARES_CONNECT)) {
continue;
}
for (int e = 0; e < nedges[c]; e++) {
int previndex = edgespecs[c][e * 2];
int currindex = edgespecs[c][e * 2 + 1];
par_connector* thisedge = par_conn_find(edges,
ncedges, previndex, currindex);
thisedge->innera = pcurrinds[previndex];
thisedge->innerb = pcurrinds[currindex];
thisedge->outera = mesh->npoints;
thisedge->outerb = mesh->npoints + 1;
par_connector* oppedge = par_conn_find(edges,
ncedges, currindex, previndex);
if (oppedge->mesh_index > color) continue;
*pdst++ = vertsx[previndex];
*pdst++ = 1 - vertsy[previndex];
if (mesh->dim == 3) {
*pdst++ = height;
}
mesh->npoints++;
*pdst++ = vertsx[currindex];
*pdst++ = 1 - vertsy[currindex];
if (mesh->dim == 3) {
*pdst++ = height;
}
mesh->npoints++;
uint16_t i0 = mesh->npoints - 1;
uint16_t i1 = mesh->npoints - 2;
uint16_t i2 = pcurrinds[previndex];
uint16_t i3 = pcurrinds[currindex];
uint16_t* ptr = mesh->conntri +
mesh->nconntriangles * 3;
*ptr++ = i2; *ptr++ = i1; *ptr++ = i0;
*ptr++ = i0; *ptr++ = i3; *ptr++ = i2;
mesh->nconntriangles += 2;
}
}
// Adjust the positions of the extrusion verts.
if (flags & PAR_MSQUARES_CONNECT) {
for (int c = 0; c < 4; c++) {
int color = vals[c];
uint16_t* pconninds = conntris_start[c];
par_msquares__mesh* mesh = mlist->meshes[color];
for (int e = 0; e < nedges[c]; e++) {
int previndex = edgespecs[c][e * 2];
int currindex = edgespecs[c][e * 2 + 1];
uint16_t i1 = pconninds[1];
uint16_t i0 = pconninds[2];
par_connector const* oppedge = par_conn_find(edges,
ncedges, currindex, previndex);
if (oppedge->mesh_index > color) continue;
int d = mesh->dim;
float* dst = mesh->points;
float const* src = oppedge->mesh->points;
dst[i0 * d + 0] = src[oppedge->innera * d + 0];
dst[i0 * d + 1] = src[oppedge->innera * d + 1];
dst[i1 * d + 0] = src[oppedge->innerb * d + 0];
dst[i1 * d + 1] = src[oppedge->innerb * d + 1];
if (d == 3) {
dst[i0 * d + 2] = src[oppedge->innera * d + 2];
dst[i1 * d + 2] = src[oppedge->innerb * d + 2];
}
pconninds += 6;
}
}
}
// Stash the bottom indices for each mesh in this cell to enable
// vertical as-you-go welding.
uint8_t* pcurrrowcells = currrowcells;
uint16_t* pcurrrowinds = currrowinds;
uint16_t const* pcurrinds = currinds;
for (int color = 0; color < ncolors; color++) {
pcurrrowcells[col] = currcell[color];
pcurrrowcells += ncols;
pcurrrowinds[col * 3 + 0] = pcurrinds[0];
pcurrrowinds[col * 3 + 1] = pcurrinds[1];
pcurrrowinds[col * 3 + 2] = pcurrinds[2];
pcurrrowinds += ncols * 3;
pcurrinds += 9;
}
// Stash some information later used by simplification.
if (flags & PAR_MSQUARES_SIMPLIFY) {
int cell = col + row * ncols;
simplification_words[cell * 2] = colors;
simplification_words[cell * 2 + 1] = counts;
}
// Advance the cursor.
nwval = neval;
swval = seval;
for (int i = 0; i < 9; i++) {
vertsx[i] += normalized_cellsize;
}
PAR_SWAP(uint8_t*, prevcell, currcell);
PAR_SWAP(uint16_t*, previnds, currinds);
}
PAR_SWAP(uint8_t*, prevrowcells, currrowcells);
PAR_SWAP(uint16_t*, prevrowinds, currrowinds);
}
free(prevrowinds);
free(prevrowcells);
free(pixels);
if (flags & PAR_MSQUARES_CLEAN) {
par_msquares__repair_tjunctions(mlist);
}
for (int m = 0; m < mlist->nmeshes; m++) {
par_msquares__mesh* mesh = mlist->meshes[m];
par__uint16list_free(mesh->tjunctions);
}
if (!(flags & PAR_MSQUARES_SIMPLIFY)) {
par_msquares__finalize(mlist);
return mlist;
}
uint8_t* simplification_blocks = PAR_CALLOC(uint8_t, nrows * ncols);
uint32_t* simplification_tris = PAR_CALLOC(uint32_t, nrows * ncols);
uint8_t* simplification_ntris = PAR_CALLOC(uint8_t, nrows * ncols);
// Perform quick-n-dirty simplification by iterating two rows at a time.
// In no way does this create the simplest possible mesh, but at least it's
// fast and easy.
for (uint32_t color = 0; color < (uint32_t) ncolors; color++) {
par_msquares__mesh* mesh = mlist->meshes[color];
// Populate the per-mesh info grids.
int ntris = 0;
for (int row = 0; row < nrows; row++) {
for (int col = 0; col < ncols; col++) {
int cell = ncols * row + col;
uint32_t colors = simplification_words[cell * 2];
uint32_t counts = simplification_words[cell * 2 + 1];
int ncelltris = 0;
int ncorners = 0;
if ((colors & 0xff) == color) {
ncelltris = counts & 0xff;
ncorners++;
}
if (((colors >> 8) & 0xff) == color) {
ncelltris += (counts >> 8) & 0xff;
ncorners++;
}
if (((colors >> 16) & 0xff) == color) {
ncelltris += (counts >> 16) & 0xff;
ncorners++;
}
if (((colors >> 24) & 0xff) == color) {
ncelltris += (counts >> 24) & 0xff;
ncorners++;
}
simplification_ntris[cell] = ncelltris;
simplification_tris[cell] = ntris;
simplification_blocks[cell] = ncorners == 4;
ntris += ncelltris;
}
}
// First figure out how many triangles we can eliminate.
int in_run = 0, start_run;
int neliminated_triangles = 0;
for (int row = 0; row < nrows - 1; row += 2) {
for (int col = 0; col < ncols; col++) {
int cell = ncols * row + col;
int a = simplification_blocks[cell];
int b = simplification_blocks[cell + ncols];
if (a && b) {
if (!in_run) {
in_run = 1;
start_run = col;
}
continue;
}
if (in_run) {
in_run = 0;
int run_width = col - start_run;
neliminated_triangles += run_width * 4 - 2;
}
}
if (in_run) {
in_run = 0;
int run_width = ncols - start_run;
neliminated_triangles += run_width * 4 - 2;
}
}
if (neliminated_triangles == 0) {
continue;
}
// Build a new index array cell-by-cell. If any given cell is 'F' and
// its neighbor to the south is also 'F', then it's part of a run.
int nnewtris = mesh->ntriangles - neliminated_triangles;
uint16_t* newtris = PAR_CALLOC(uint16_t, nnewtris * 3);
uint16_t* pnewtris = newtris;
in_run = 0;
uint16_t* tris = mesh->triangles;
for (int row = 0; row < nrows - 1; row += 2) {
for (int col = 0; col < ncols; col++) {
int cell = ncols * row + col;
int south = cell + ncols;
int a = simplification_blocks[cell];
int b = simplification_blocks[south];
if (a && b) {
if (!in_run) {
in_run = 1;
start_run = col;
}
continue;
}
if (in_run) {
in_run = 0;
int nw_cell = ncols * row + start_run;
int ne_cell = ncols * row + col - 1;
int sw_cell = nw_cell + ncols;
int se_cell = ne_cell + ncols;
int nw_tri = simplification_tris[nw_cell];
int ne_tri = simplification_tris[ne_cell];
int sw_tri = simplification_tris[sw_cell];
int se_tri = simplification_tris[se_cell];
int nw_corner = nw_tri * 3 + 5;
int ne_corner = ne_tri * 3 + 2;
int sw_corner = sw_tri * 3 + 0;
int se_corner = se_tri * 3 + 1;
*pnewtris++ = tris[nw_corner];
*pnewtris++ = tris[sw_corner];
*pnewtris++ = tris[se_corner];
*pnewtris++ = tris[se_corner];
*pnewtris++ = tris[ne_corner];
*pnewtris++ = tris[nw_corner];
}
int ncelltris = simplification_ntris[cell];
int celltri = simplification_tris[cell];
for (int t = 0; t < ncelltris; t++, celltri++) {
*pnewtris++ = tris[celltri * 3];
*pnewtris++ = tris[celltri * 3 + 1];
*pnewtris++ = tris[celltri * 3 + 2];
}
ncelltris = simplification_ntris[south];
celltri = simplification_tris[south];
for (int t = 0; t < ncelltris; t++, celltri++) {
*pnewtris++ = tris[celltri * 3];
*pnewtris++ = tris[celltri * 3 + 1];
*pnewtris++ = tris[celltri * 3 + 2];
}
}
if (in_run) {
in_run = 0;
int nw_cell = ncols * row + start_run;
int ne_cell = ncols * row + ncols - 1;
int sw_cell = nw_cell + ncols;
int se_cell = ne_cell + ncols;
int nw_tri = simplification_tris[nw_cell];
int ne_tri = simplification_tris[ne_cell];
int sw_tri = simplification_tris[sw_cell];
int se_tri = simplification_tris[se_cell];
int nw_corner = nw_tri * 3 + 5;
int ne_corner = ne_tri * 3 + 2;
int sw_corner = sw_tri * 3 + 0;
int se_corner = se_tri * 3 + 1;
*pnewtris++ = tris[nw_corner];
*pnewtris++ = tris[sw_corner];
*pnewtris++ = tris[se_corner];
*pnewtris++ = tris[se_corner];
*pnewtris++ = tris[ne_corner];
*pnewtris++ = tris[nw_corner];
}
}
mesh->ntriangles -= neliminated_triangles;
free(mesh->triangles);
mesh->triangles = newtris;
}
free(simplification_blocks);
free(simplification_ntris);
free(simplification_tris);
free(simplification_words);
par_msquares__finalize(mlist);
for (int i = 0; i < mlist->nmeshes; i++) {
par_remove_unreferenced_verts(mlist->meshes[i]);
}
return mlist;
}
void par_msquares_free_boundary(par_msquares_boundary* polygon)
{
free(polygon->points);
free(polygon->chains);
free(polygon->lengths);
free(polygon);
}
typedef struct par__hedge_s {
uint32_t key;
struct par__hvert_s* endvert;
struct par__hedge_s* opposite;
struct par__hedge_s* next;
struct par__hedge_s* prev;
} par__hedge;
typedef struct par__hvert_s {
par__hedge* incoming;
} par__hvert;
typedef struct {
par_msquares_mesh const* mesh;
par__hvert* verts;
par__hedge* edges;
par__hedge** sorted_edges;
} par__hemesh;
static int par__hedge_cmp(const void *arg0, const void *arg1)
{
par__hedge* he0 = *((par__hedge**) arg0);
par__hedge* he1 = *((par__hedge**) arg1);
if (he0->key < he1->key) return -1;
if (he0->key > he1->key) return 1;
return 0;
}
static par__hedge* par__hedge_find(par__hemesh* hemesh, uint32_t key)
{
par__hedge target = {0};
target.key = key;
par__hedge* ptarget = &target;
int nedges = hemesh->mesh->ntriangles * 3;
par__hedge** result = (par__hedge**) bsearch(&ptarget, hemesh->sorted_edges,
nedges, sizeof(par__hedge*), par__hedge_cmp);
return result ? *result : 0;
}
static uint32_t par__hedge_key(par__hvert* a, par__hvert* b, par__hvert* s)
{
uint32_t ai = a - s;
uint32_t bi = b - s;
return (ai << 16) | bi;
}
par_msquares_boundary* par_msquares_extract_boundary(
par_msquares_mesh const* mesh)
{
par_msquares_boundary* result = PAR_CALLOC(par_msquares_boundary, 1);
par__hemesh hemesh = {0};
hemesh.mesh = mesh;
int nedges = mesh->ntriangles * 3;
// Populate all fields of verts and edges, except opposite.
hemesh.edges = PAR_CALLOC(par__hedge, nedges);
par__hvert* hverts = hemesh.verts = PAR_CALLOC(par__hvert, mesh->npoints);
par__hedge* edge = hemesh.edges;
uint16_t const* tri = mesh->triangles;
for (int n = 0; n < mesh->ntriangles; n++, edge += 3, tri += 3) {
edge[0].endvert = hverts + tri[1];
edge[1].endvert = hverts + tri[2];
edge[2].endvert = hverts + tri[0];
hverts[tri[1]].incoming = edge + 0;
hverts[tri[2]].incoming = edge + 1;
hverts[tri[0]].incoming = edge + 2;
edge[0].next = edge + 1;
edge[1].next = edge + 2;
edge[2].next = edge + 0;
edge[0].prev = edge + 2;
edge[1].prev = edge + 0;
edge[2].prev = edge + 1;
edge[0].key = par__hedge_key(edge[2].endvert, edge[0].endvert, hverts);
edge[1].key = par__hedge_key(edge[0].endvert, edge[1].endvert, hverts);
edge[2].key = par__hedge_key(edge[1].endvert, edge[2].endvert, hverts);
}
// Sort the edges according to their key.
hemesh.sorted_edges = PAR_CALLOC(par__hedge*, mesh->ntriangles * 3);
for (int n = 0; n < nedges; n++) {
hemesh.sorted_edges[n] = hemesh.edges + n;
}
qsort(hemesh.sorted_edges, nedges, sizeof(par__hedge*), par__hedge_cmp);
// Populate the "opposite" field in each edge.
for (int n = 0; n < nedges; n++) {
par__hedge* edge = hemesh.edges + n;
par__hedge* prev = edge->prev;
par__hvert* start = edge->endvert;
par__hvert* end = prev->endvert;
uint32_t key = par__hedge_key(start, end, hverts);
edge->opposite = par__hedge_find(&hemesh, key);
}
// Re-use the sorted_edges array, filling it with boundary edges only.
// Also create a mapping table to consolidate all boundary verts.
int nborders = 0;
for (int n = 0; n < nedges; n++) {
par__hedge* edge = hemesh.edges + n;
if (!edge->opposite) {
hemesh.sorted_edges[nborders++] = edge;
}
}
// Allocate for the worst case (all separate triangles).
// We'll adjust the lengths later.
result->nchains = nborders / 3;
result->npoints = nborders + result->nchains;
result->points = PAR_CALLOC(float, 2 * result->npoints);
result->chains = PAR_CALLOC(float*, result->nchains);
result->lengths = PAR_CALLOC(uint16_t, result->nchains);
// Iterate over each polyline.
edge = hemesh.sorted_edges[0];
int pt = 0;
int nwritten = 0;
int nchains = 0;
while (1) {
float* points = result->points;
par__hedge* orig = edge;
uint16_t index = edge->prev->endvert - hverts;
result->chains[nchains] = points + pt;
result->lengths[nchains]++;
points[pt++] = mesh->points[index * mesh->dim];
points[pt++] = mesh->points[index * mesh->dim + 1];
while (1) {
index = edge->endvert - hverts;
edge->key = 0;
nwritten++;
result->lengths[nchains]++;
points[pt++] = mesh->points[index * mesh->dim];
points[pt++] = mesh->points[index * mesh->dim + 1];
par__hedge* next = edge->next;
while (next != edge) {
if (!next->opposite) {
break;
}
next = next->opposite->next;
}
edge = next;
if (edge == orig) {
break;
}
}
nchains++;
if (nwritten >= nborders) {
break;
}
for (int i = 0; i < nborders; i++) {
edge = hemesh.sorted_edges[i];
if (edge->key) {
break;
}
}
}
result->npoints = pt / 2;
result->nchains = nchains;
free(hemesh.verts);
free(hemesh.edges);
free(hemesh.sorted_edges);
return result;
}
#endif // PAR_MSQUARES_IMPLEMENTATION
#endif // PAR_MSQUARES_H