1 // MSQUARES :: https://github.com/prideout/par
2 // Converts fp32 grayscale images, or 8-bit color images, into triangles.
3 //
4 // For grayscale images, a threshold is specified to determine insideness.
5 // For color images, an exact color is specified to determine insideness.
6 // Color images can be r8, rg16, rgb24, or rgba32. For a visual overview of
7 // the API and all the flags, see:
8 //
9 // http://github.prideout.net/marching-cubes/
10 //
11 // The MIT License
12 // Copyright (c) 2015 Philip Rideout
13
14 #include <stdint.h>
15
16 // -----------------------------------------------------------------------------
17 // BEGIN PUBLIC API
18 // -----------------------------------------------------------------------------
19
20 typedef uint8_t par_byte;
21
22 typedef struct par_msquares_meshlist_s par_msquares_meshlist;
23
24 // Encapsulates the results of a marching squares operation.
25 typedef struct {
26 float* points; // pointer to XY (or XYZ) vertex coordinates
27 int npoints; // number of vertex coordinates
28 uint16_t* triangles; // pointer to 3-tuples of vertex indices
29 int ntriangles; // number of 3-tuples
30 int dim; // number of floats per point (either 2 or 3)
31 int nconntriangles; // internal use only
32 } par_msquares_mesh;
33
34 // Reverses the "insideness" test.
35 #define PAR_MSQUARES_INVERT (1 << 0)
36
37 // Returns a meshlist with two meshes: one for the inside, one for the outside.
38 #define PAR_MSQUARES_DUAL (1 << 1)
39
40 // Returned meshes have 3-tuple coordinates instead of 2-tuples. When using
41 // from_color, the Z coordinate represents the alpha value of the color. With
42 // from_grayscale, the Z coordinate represents the value of the nearest pixel in
43 // the source image.
44 #define PAR_MSQUARES_HEIGHTS (1 << 2)
45
46 // Applies a step function to the Z coordinates. Requires HEIGHTS and DUAL.
47 #define PAR_MSQUARES_SNAP (1 << 3)
48
49 // Adds extrusion triangles to each mesh other than the lowest mesh. Requires
50 // the PAR_MSQUARES_HEIGHTS flag to be present.
51 #define PAR_MSQUARES_CONNECT (1 << 4)
52
53 // Enables quick & dirty (not best) simpification of the returned mesh.
54 #define PAR_MSQUARES_SIMPLIFY (1 << 5)
55
56 par_msquares_meshlist* par_msquares_from_grayscale(float const* data, int width,
57 int height, int cellsize, float threshold, int flags);
58
59 par_msquares_meshlist* par_msquares_from_color(par_byte const* data, int width,
60 int height, int cellsize, uint32_t color, int bpp, int flags);
61
62 typedef int (*par_msquares_inside_fn)(int, void*);
63 typedef float (*par_msquares_height_fn)(float, float, void*);
64
65 par_msquares_meshlist* par_msquares_from_function(int width, int height,
66 int cellsize, int flags, void* context, par_msquares_inside_fn insidefn,
67 par_msquares_height_fn heightfn);
68
69 par_msquares_mesh const* par_msquares_get_mesh(par_msquares_meshlist*, int n);
70
71 int par_msquares_get_count(par_msquares_meshlist*);
72
73 void par_msquares_free(par_msquares_meshlist*);
74
75 // -----------------------------------------------------------------------------
76 // END PUBLIC API
77 // -----------------------------------------------------------------------------
78
79 #ifdef PAR_MSQUARES_IMPLEMENTATION
80
81 #include <stdlib.h>
82 #include <assert.h>
83 #include <float.h>
84
85 #define MIN(a, b) (a > b ? b : a)
86 #define MAX(a, b) (a > b ? a : b)
87 #define CLAMP(v, lo, hi) MAX(lo, MIN(hi, v))
88
89 struct par_msquares_meshlist_s {
90 int nmeshes;
91 par_msquares_mesh** meshes;
92 };
93
94 static int** point_table = 0;
95 static int** triangle_table = 0;
96
par_init_tables()97 static void par_init_tables()
98 {
99 point_table = (int**)malloc(16 * sizeof(int*));
100 triangle_table = (int**)malloc(16 * sizeof(int*));
101
102 char const* CODE_TABLE =
103 "0 0\n"
104 "1 1 0 1 7\n"
105 "2 1 1 2 3\n"
106 "3 2 0 2 3 3 7 0\n"
107 "4 1 7 5 6\n"
108 "5 2 0 1 5 5 6 0\n"
109 "6 2 1 2 3 7 5 6\n"
110 "7 3 0 2 3 0 3 5 0 5 6\n"
111 "8 1 3 4 5\n"
112 "9 4 0 1 3 0 3 4 0 4 5 0 5 7\n"
113 "a 2 1 2 4 4 5 1\n"
114 "b 3 0 2 4 0 4 5 0 5 7\n"
115 "c 2 7 3 4 4 6 7\n"
116 "d 3 0 1 3 0 3 4 0 4 6\n"
117 "e 3 1 2 4 1 4 6 1 6 7\n"
118 "f 2 0 2 4 4 6 0\n";
119
120 char const* table_token = CODE_TABLE;
121 for (int i = 0; i < 16; i++) {
122 char code = *table_token;
123 assert(i == code - (code >= 'a' ? ('a' - 0xa) : '0'));
124 table_token += 2;
125 int ntris = *table_token - '0';
126 table_token += 2;
127 int* sqrtris = triangle_table[i] =
128 (int*)malloc((ntris + 1) * 3 * sizeof(int));
129 sqrtris[0] = ntris;
130 int mask = 0;
131 int* sqrpts = point_table[i] = (int*)malloc(7 * sizeof(int));
132 sqrpts[0] = 0;
133 for (int j = 0; j < ntris * 3; j++, table_token += 2) {
134 int midp = *table_token - '0';
135 int bit = 1 << midp;
136 if (!(mask & bit)) {
137 mask |= bit;
138 sqrpts[++sqrpts[0]] = midp;
139 }
140 sqrtris[j + 1] = midp;
141 }
142 }
143 }
144
145 typedef struct {
146 float const* data;
147 float threshold;
148 int width;
149 int height;
150 } gray_context;
151
gray_inside(int location,void * contextptr)152 static int gray_inside(int location, void* contextptr)
153 {
154 gray_context* context = (gray_context*) contextptr;
155 return context->data[location] > context->threshold;
156 }
157
gray_height(float x,float y,void * contextptr)158 static float gray_height(float x, float y, void* contextptr)
159 {
160 gray_context* context = (gray_context*) contextptr;
161 int i = CLAMP(context->width * x, 0, context->width - 1);
162 int j = CLAMP(context->height * y, 0, context->height - 1);
163 return context->data[i + j * context->width];
164 }
165
166 typedef struct {
167 par_byte const* data;
168 par_byte color[4];
169 int bpp;
170 int width;
171 int height;
172 } color_context;
173
color_inside(int location,void * contextptr)174 static int color_inside(int location, void* contextptr)
175 {
176 color_context* context = (color_context*) contextptr;
177 par_byte const* data = context->data + location * context->bpp;
178 for (int i = 0; i < context->bpp; i++) {
179 if (data[i] != context->color[i]) {
180 return 0;
181 }
182 }
183 return 1;
184 }
185
color_height(float x,float y,void * contextptr)186 static float color_height(float x, float y, void* contextptr)
187 {
188 color_context* context = (color_context*) contextptr;
189 assert(context->bpp == 4);
190 int i = CLAMP(context->width * x, 0, context->width - 1);
191 int j = CLAMP(context->height * y, 0, context->height - 1);
192 int k = i + j * context->width;
193 return context->data[k * 4 + 3] / 255.0;
194 }
195
par_msquares_from_color(par_byte const * data,int width,int height,int cellsize,uint32_t color,int bpp,int flags)196 par_msquares_meshlist* par_msquares_from_color(par_byte const* data, int width,
197 int height, int cellsize, uint32_t color, int bpp, int flags)
198 {
199 color_context context;
200 context.bpp = bpp;
201 context.color[0] = (color >> 16) & 0xff;
202 context.color[1] = (color >> 8) & 0xff;
203 context.color[2] = (color & 0xff);
204 context.color[3] = (color >> 24) & 0xff;
205 context.data = data;
206 context.width = width;
207 context.height = height;
208 return par_msquares_from_function(
209 width, height, cellsize, flags, &context, color_inside, color_height);
210 }
211
par_msquares_from_grayscale(float const * data,int width,int height,int cellsize,float threshold,int flags)212 par_msquares_meshlist* par_msquares_from_grayscale(float const* data, int width,
213 int height, int cellsize, float threshold, int flags)
214 {
215 gray_context context;
216 context.width = width;
217 context.height = height;
218 context.data = data;
219 context.threshold = threshold;
220 return par_msquares_from_function(
221 width, height, cellsize, flags, &context, gray_inside, gray_height);
222 }
223
par_msquares_get_mesh(par_msquares_meshlist * mlist,int mindex)224 par_msquares_mesh const* par_msquares_get_mesh(
225 par_msquares_meshlist* mlist, int mindex)
226 {
227 assert(mlist && mindex < mlist->nmeshes);
228 return mlist->meshes[mindex];
229 }
230
par_msquares_get_count(par_msquares_meshlist * mlist)231 int par_msquares_get_count(par_msquares_meshlist* mlist)
232 {
233 assert(mlist);
234 return mlist->nmeshes;
235 }
236
par_msquares_free(par_msquares_meshlist * mlist)237 void par_msquares_free(par_msquares_meshlist* mlist)
238 {
239 par_msquares_mesh** meshes = mlist->meshes;
240 for (int i = 0; i < mlist->nmeshes; i++) {
241 free(meshes[i]->points);
242 free(meshes[i]->triangles);
243 free(meshes[i]);
244 }
245 free(meshes);
246 free(mlist);
247 }
248
249 // Combine multiple meshlists by moving mesh pointers, and optionally apply
250 // a "snap" operation that assigns a single Z value across all verts in each
251 // mesh. The Z value determined by the mesh's position in the final mesh list.
par_msquares_merge(par_msquares_meshlist ** lists,int count,int snap)252 static par_msquares_meshlist* par_msquares_merge(par_msquares_meshlist** lists,
253 int count, int snap)
254 {
255 par_msquares_meshlist* merged = (par_msquares_meshlist*)malloc(sizeof(par_msquares_meshlist));
256 merged->nmeshes = 0;
257 for (int i = 0; i < count; i++) {
258 merged->nmeshes += lists[i]->nmeshes;
259 }
260 merged->meshes = (par_msquares_mesh**)malloc(sizeof(par_msquares_mesh*) * merged->nmeshes);
261 par_msquares_mesh** pmesh = merged->meshes;
262 for (int i = 0; i < count; i++) {
263 par_msquares_meshlist* meshlist = lists[i];
264 for (int j = 0; j < meshlist->nmeshes; j++) {
265 *pmesh++ = meshlist->meshes[j];
266 }
267 free(meshlist);
268 }
269 if (!snap) {
270 return merged;
271 }
272 pmesh = merged->meshes;
273 float zmin = FLT_MAX;
274 float zmax = -zmin;
275 for (int i = 0; i < merged->nmeshes; i++, pmesh++) {
276 float* pzed = (*pmesh)->points + 2;
277 for (int j = 0; j < (*pmesh)->npoints; j++, pzed += 3) {
278 zmin = MIN(*pzed, zmin);
279 zmax = MAX(*pzed, zmax);
280 }
281 }
282 float zextent = zmax - zmin;
283 pmesh = merged->meshes;
284 for (int i = 0; i < merged->nmeshes; i++, pmesh++) {
285 float* pzed = (*pmesh)->points + 2;
286 float zed = zmin + zextent * i / (merged->nmeshes - 1);
287 for (int j = 0; j < (*pmesh)->npoints; j++, pzed += 3) {
288 *pzed = zed;
289 }
290 }
291 if (!(snap & PAR_MSQUARES_CONNECT)) {
292 return merged;
293 }
294 for (int i = 1; i < merged->nmeshes; i++) {
295 par_msquares_mesh* mesh = merged->meshes[i];
296
297 // Find all extrusion points. This is tightly coupled to the
298 // tessellation code, which generates two "connector" triangles for each
299 // extruded edge. The first two verts of the second triangle are the
300 // verts that need to be displaced.
301 char* markers = (char*)calloc(mesh->npoints, 1);
302 int tri = mesh->ntriangles - mesh->nconntriangles;
303 while (tri < mesh->ntriangles) {
304 markers[mesh->triangles[tri * 3 + 3]] = 1;
305 markers[mesh->triangles[tri * 3 + 4]] = 1;
306 tri += 2;
307 }
308
309 // Displace all extrusion points down to the previous level.
310 float zed = zmin + zextent * (i - 1) / (merged->nmeshes - 1);
311 float* pzed = mesh->points + 2;
312 for (int j = 0; j < mesh->npoints; j++, pzed += 3) {
313 if (markers[j]) {
314 *pzed = zed;
315 }
316 }
317 free(markers);
318 }
319 return merged;
320 }
321
par_msquares_from_function(int width,int height,int cellsize,int flags,void * context,par_msquares_inside_fn insidefn,par_msquares_height_fn heightfn)322 par_msquares_meshlist* par_msquares_from_function(int width, int height,
323 int cellsize, int flags, void* context, par_msquares_inside_fn insidefn,
324 par_msquares_height_fn heightfn)
325 {
326 assert(width > 0 && width % cellsize == 0);
327 assert(height > 0 && height % cellsize == 0);
328
329 if (flags & PAR_MSQUARES_DUAL) {
330 int connect = flags & PAR_MSQUARES_CONNECT;
331 int snap = flags & PAR_MSQUARES_SNAP;
332 int heights = flags & PAR_MSQUARES_HEIGHTS;
333 if (!heights) {
334 snap = connect = 0;
335 }
336 flags ^= PAR_MSQUARES_INVERT;
337 flags &= ~PAR_MSQUARES_DUAL;
338 flags &= ~PAR_MSQUARES_CONNECT;
339 par_msquares_meshlist* m[2];
340 m[0] = par_msquares_from_function(width, height, cellsize, flags,
341 context, insidefn, heightfn);
342 flags ^= PAR_MSQUARES_INVERT;
343 if (connect) {
344 flags |= PAR_MSQUARES_CONNECT;
345 }
346 m[1] = par_msquares_from_function(width, height, cellsize, flags,
347 context, insidefn, heightfn);
348 return par_msquares_merge(m, 2, snap | connect);
349 }
350
351 int invert = flags & PAR_MSQUARES_INVERT;
352
353 // Create the two code tables if we haven't already. These are tables of
354 // fixed constants, so it's embarassing that we use dynamic memory
355 // allocation for them. However it's easy and it's one-time-only.
356
357 if (!point_table) {
358 par_init_tables();
359 }
360
361 // Allocate the meshlist and the first mesh.
362
363 par_msquares_meshlist* mlist = (par_msquares_meshlist*)malloc(sizeof(par_msquares_meshlist));
364 mlist->nmeshes = 1;
365 mlist->meshes = (par_msquares_mesh**)malloc(sizeof(par_msquares_mesh*));
366 mlist->meshes[0] = (par_msquares_mesh*)malloc(sizeof(par_msquares_mesh));
367 par_msquares_mesh* mesh = mlist->meshes[0];
368 mesh->dim = (flags & PAR_MSQUARES_HEIGHTS) ? 3 : 2;
369 int ncols = width / cellsize;
370 int nrows = height / cellsize;
371
372 // Worst case is four triangles and six verts per cell, so allocate that
373 // much.
374
375 int maxtris = ncols * nrows * 4;
376 int maxpts = ncols * nrows * 6;
377 int maxedges = ncols * nrows * 2;
378
379 // However, if we include extrusion triangles for boundary edges,
380 // we need space for another 4 triangles and 4 points per cell.
381
382 uint16_t* conntris = 0;
383 int nconntris = 0;
384 uint16_t* edgemap = 0;
385 if (flags & PAR_MSQUARES_CONNECT) {
386 conntris = (uint16_t*)malloc(maxedges * 6 * sizeof(uint16_t));
387 maxtris += maxedges * 2;
388 maxpts += maxedges * 2;
389 edgemap = (uint16_t*)malloc(maxpts * sizeof(uint16_t));
390 for (int i = 0; i < maxpts; i++) {
391 edgemap[i] = 0xffff;
392 }
393 }
394
395 uint16_t* tris = (uint16_t*)malloc(maxtris * 3 * sizeof(uint16_t));
396 int ntris = 0;
397 float* pts = (float*)malloc(maxpts * mesh->dim * sizeof(float));
398 int npts = 0;
399
400 // The "verts" x/y/z arrays are the 4 corners and 4 midpoints around the
401 // square, in counter-clockwise order. The origin of "triangle space" is at
402 // the lower-left, although we expect the image data to be in raster order
403 // (starts at top-left).
404
405 float normalization = 1.0f / MAX(width, height);
406 float normalized_cellsize = cellsize * normalization;
407 int maxrow = (height - 1) * width;
408 uint16_t* ptris = tris;
409 uint16_t* pconntris = conntris;
410 float* ppts = pts;
411 float vertsx[8], vertsy[8];
412 uint8_t* prevrowmasks = (uint8_t*)calloc(ncols, 1);
413 int* prevrowinds = (int*)calloc(sizeof(int) * ncols * 3, 1);
414
415 // If simplification is enabled, we need to track all 'F' cells and their
416 // repsective triangle indices.
417 uint8_t* simplification_codes = 0;
418 uint16_t* simplification_tris = 0;
419 uint8_t* simplification_ntris = 0;
420 if (flags & PAR_MSQUARES_SIMPLIFY) {
421 simplification_codes = (uint8_t*)malloc(nrows * ncols);
422 simplification_tris = (uint16_t*)malloc(nrows * ncols * sizeof(uint16_t));
423 simplification_ntris = (uint8_t*)malloc(nrows * ncols);
424 }
425
426 // Do the march!
427 for (int row = 0; row < nrows; row++) {
428 vertsx[0] = vertsx[6] = vertsx[7] = 0;
429 vertsx[1] = vertsx[5] = 0.5 * normalized_cellsize;
430 vertsx[2] = vertsx[3] = vertsx[4] = normalized_cellsize;
431 vertsy[0] = vertsy[1] = vertsy[2] = normalized_cellsize * (row + 1);
432 vertsy[4] = vertsy[5] = vertsy[6] = normalized_cellsize * row;
433 vertsy[3] = vertsy[7] = normalized_cellsize * (row + 0.5);
434
435 int northi = row * cellsize * width;
436 int southi = MIN(northi + cellsize * width, maxrow);
437 int northwest = invert ^ insidefn(northi, context);
438 int southwest = invert ^ insidefn(southi, context);
439 int previnds[8] = {0};
440 uint8_t prevmask = 0;
441
442 for (int col = 0; col < ncols; col++) {
443 northi += cellsize;
444 southi += cellsize;
445 if (col == ncols - 1) {
446 northi--;
447 southi--;
448 }
449
450 int northeast = invert ^ insidefn(northi, context);
451 int southeast = invert ^ insidefn(southi, context);
452 int code = southwest | (southeast << 1) | (northwest << 2) |
453 (northeast << 3);
454
455 int const* pointspec = point_table[code];
456 int ptspeclength = *pointspec++;
457 int currinds[8] = {0};
458 uint8_t mask = 0;
459 uint8_t prevrowmask = prevrowmasks[col];
460 while (ptspeclength--) {
461 int midp = *pointspec++;
462 int bit = 1 << midp;
463 mask |= bit;
464
465 // The following six conditionals perform welding to reduce the
466 // number of vertices. The first three perform welding with the
467 // cell to the west; the latter three perform welding with the
468 // cell to the north.
469
470 if (bit == 1 && (prevmask & 4)) {
471 currinds[midp] = previnds[2];
472 continue;
473 }
474 if (bit == 128 && (prevmask & 8)) {
475 currinds[midp] = previnds[3];
476 continue;
477 }
478 if (bit == 64 && (prevmask & 16)) {
479 currinds[midp] = previnds[4];
480 continue;
481 }
482 if (bit == 16 && (prevrowmask & 4)) {
483 currinds[midp] = prevrowinds[col * 3 + 2];
484 continue;
485 }
486 if (bit == 32 && (prevrowmask & 2)) {
487 currinds[midp] = prevrowinds[col * 3 + 1];
488 continue;
489 }
490 if (bit == 64 && (prevrowmask & 1)) {
491 currinds[midp] = prevrowinds[col * 3 + 0];
492 continue;
493 }
494
495 ppts[0] = vertsx[midp];
496 ppts[1] = vertsy[midp];
497
498 // Adjust the midpoints to a more exact crossing point.
499 if (midp == 1) {
500 int begin = southi - cellsize / 2;
501 int previous = 0;
502 for (int i = 0; i < cellsize; i++) {
503 int offset = begin + i / 2 * ((i % 2) ? -1 : 1);
504 int inside = insidefn(offset, context);
505 if (i > 0 && inside != previous) {
506 ppts[0] = normalization *
507 (col * cellsize + offset - southi + cellsize);
508 break;
509 }
510 previous = inside;
511 }
512 } else if (midp == 5) {
513 int begin = northi - cellsize / 2;
514 int previous = 0;
515 for (int i = 0; i < cellsize; i++) {
516 int offset = begin + i / 2 * ((i % 2) ? -1 : 1);
517 int inside = insidefn(offset, context);
518 if (i > 0 && inside != previous) {
519 ppts[0] = normalization *
520 (col * cellsize + offset - northi + cellsize);
521 break;
522 }
523 previous = inside;
524 }
525 } else if (midp == 3) {
526 int begin = northi + width * cellsize / 2;
527 int previous = 0;
528 for (int i = 0; i < cellsize; i++) {
529 int offset = begin +
530 width * (i / 2 * ((i % 2) ? -1 : 1));
531 int inside = insidefn(offset, context);
532 if (i > 0 && inside != previous) {
533 ppts[1] = normalization *
534 (row * cellsize +
535 (offset - northi) / (float) width);
536 break;
537 }
538 previous = inside;
539 }
540 } else if (midp == 7) {
541 int begin = northi + width * cellsize / 2 - cellsize;
542 int previous = 0;
543 for (int i = 0; i < cellsize; i++) {
544 int offset = begin +
545 width * (i / 2 * ((i % 2) ? -1 : 1));
546 int inside = insidefn(offset, context);
547 if (i > 0 && inside != previous) {
548 ppts[1] = normalization *
549 (row * cellsize +
550 (offset - northi - cellsize) / (float) width);
551 break;
552 }
553 previous = inside;
554 }
555 }
556
557 if (mesh->dim == 3) {
558 ppts[2] = heightfn(ppts[0], ppts[1], context);
559 }
560
561 ppts += mesh->dim;
562 currinds[midp] = npts++;
563 }
564
565 int const* trianglespec = triangle_table[code];
566 int trispeclength = *trianglespec++;
567
568 if (flags & PAR_MSQUARES_SIMPLIFY) {
569 simplification_codes[ncols * row + col] = code;
570 simplification_tris[ncols * row + col] = ntris;
571 simplification_ntris[ncols * row + col] = trispeclength;
572 }
573
574 // Add triangles.
575 while (trispeclength--) {
576 int a = *trianglespec++;
577 int b = *trianglespec++;
578 int c = *trianglespec++;
579 *ptris++ = currinds[c];
580 *ptris++ = currinds[b];
581 *ptris++ = currinds[a];
582 ntris++;
583 }
584
585 // Create two extrusion triangles for each boundary edge.
586 if (flags & PAR_MSQUARES_CONNECT) {
587 trianglespec = triangle_table[code];
588 trispeclength = *trianglespec++;
589 while (trispeclength--) {
590 int a = *trianglespec++;
591 int b = *trianglespec++;
592 int c = *trianglespec++;
593 int i = currinds[a];
594 int j = currinds[b];
595 int k = currinds[c];
596 int u = 0, v = 0, w = 0;
597 if ((a % 2) && (b % 2)) {
598 u = v = 1;
599 } else if ((a % 2) && (c % 2)) {
600 u = w = 1;
601 } else if ((b % 2) && (c % 2)) {
602 v = w = 1;
603 } else {
604 continue;
605 }
606 if (u && edgemap[i] == 0xffff) {
607 for (int d = 0; d < mesh->dim; d++) {
608 *ppts++ = pts[i * mesh->dim + d];
609 }
610 edgemap[i] = npts++;
611 }
612 if (v && edgemap[j] == 0xffff) {
613 for (int d = 0; d < mesh->dim; d++) {
614 *ppts++ = pts[j * mesh->dim + d];
615 }
616 edgemap[j] = npts++;
617 }
618 if (w && edgemap[k] == 0xffff) {
619 for (int d = 0; d < mesh->dim; d++) {
620 *ppts++ = pts[k * mesh->dim + d];
621 }
622 edgemap[k] = npts++;
623 }
624 if ((a % 2) && (b % 2)) {
625 *pconntris++ = i;
626 *pconntris++ = j;
627 *pconntris++ = edgemap[j];
628 *pconntris++ = edgemap[j];
629 *pconntris++ = edgemap[i];
630 *pconntris++ = i;
631 } else if ((a % 2) && (c % 2)) {
632 *pconntris++ = edgemap[k];
633 *pconntris++ = k;
634 *pconntris++ = i;
635 *pconntris++ = edgemap[i];
636 *pconntris++ = edgemap[k];
637 *pconntris++ = i;
638 } else if ((b % 2) && (c % 2)) {
639 *pconntris++ = j;
640 *pconntris++ = k;
641 *pconntris++ = edgemap[k];
642 *pconntris++ = edgemap[k];
643 *pconntris++ = edgemap[j];
644 *pconntris++ = j;
645 }
646 nconntris += 2;
647 }
648 }
649
650 // Prepare for the next cell.
651 prevrowmasks[col] = mask;
652 prevrowinds[col * 3 + 0] = currinds[0];
653 prevrowinds[col * 3 + 1] = currinds[1];
654 prevrowinds[col * 3 + 2] = currinds[2];
655 prevmask = mask;
656 northwest = northeast;
657 southwest = southeast;
658 for (int i = 0; i < 8; i++) {
659 previnds[i] = currinds[i];
660 vertsx[i] += normalized_cellsize;
661 }
662 }
663 }
664 free(edgemap);
665 free(prevrowmasks);
666 free(prevrowinds);
667
668 // Perform quick-n-dirty simplification by iterating two rows at a time.
669 // In no way does this create the simplest possible mesh, but at least it's
670 // fast and easy.
671 if (flags & PAR_MSQUARES_SIMPLIFY) {
672 int in_run = 0, start_run;
673
674 // First figure out how many triangles we can eliminate.
675 int neliminated_triangles = 0;
676 for (int row = 0; row < nrows - 1; row += 2) {
677 for (int col = 0; col < ncols; col++) {
678 int a = simplification_codes[ncols * row + col] == 0xf;
679 int b = simplification_codes[ncols * row + col + ncols] == 0xf;
680 if (a && b) {
681 if (!in_run) {
682 in_run = 1;
683 start_run = col;
684 }
685 continue;
686 }
687 if (in_run) {
688 in_run = 0;
689 int run_width = col - start_run;
690 neliminated_triangles += run_width * 4 - 2;
691 }
692 }
693 if (in_run) {
694 in_run = 0;
695 int run_width = ncols - start_run;
696 neliminated_triangles += run_width * 4 - 2;
697 }
698 }
699
700 // Build a new index array cell-by-cell. If any given cell is 'F' and
701 // its neighbor to the south is also 'F', then it's part of a run.
702 int nnewtris = ntris + nconntris - neliminated_triangles;
703 uint16_t* newtris = (uint16_t*)malloc(nnewtris * 3 * sizeof(uint16_t));
704 uint16_t* pnewtris = newtris;
705 in_run = 0;
706 for (int row = 0; row < nrows - 1; row += 2) {
707 for (int col = 0; col < ncols; col++) {
708 int cell = ncols * row + col;
709 int south = cell + ncols;
710 int a = simplification_codes[cell] == 0xf;
711 int b = simplification_codes[south] == 0xf;
712 if (a && b) {
713 if (!in_run) {
714 in_run = 1;
715 start_run = col;
716 }
717 continue;
718 }
719 if (in_run) {
720 in_run = 0;
721 int nw_cell = ncols * row + start_run;
722 int ne_cell = ncols * row + col - 1;
723 int sw_cell = nw_cell + ncols;
724 int se_cell = ne_cell + ncols;
725 int nw_tri = simplification_tris[nw_cell];
726 int ne_tri = simplification_tris[ne_cell];
727 int sw_tri = simplification_tris[sw_cell];
728 int se_tri = simplification_tris[se_cell];
729 int nw_corner = nw_tri * 3 + 4;
730 int ne_corner = ne_tri * 3 + 0;
731 int sw_corner = sw_tri * 3 + 2;
732 int se_corner = se_tri * 3 + 1;
733 *pnewtris++ = tris[se_corner];
734 *pnewtris++ = tris[sw_corner];
735 *pnewtris++ = tris[nw_corner];
736 *pnewtris++ = tris[nw_corner];
737 *pnewtris++ = tris[ne_corner];
738 *pnewtris++ = tris[se_corner];
739 }
740 int ncelltris = simplification_ntris[cell];
741 int celltri = simplification_tris[cell];
742 for (int t = 0; t < ncelltris; t++, celltri++) {
743 *pnewtris++ = tris[celltri * 3];
744 *pnewtris++ = tris[celltri * 3 + 1];
745 *pnewtris++ = tris[celltri * 3 + 2];
746 }
747 ncelltris = simplification_ntris[south];
748 celltri = simplification_tris[south];
749 for (int t = 0; t < ncelltris; t++, celltri++) {
750 *pnewtris++ = tris[celltri * 3];
751 *pnewtris++ = tris[celltri * 3 + 1];
752 *pnewtris++ = tris[celltri * 3 + 2];
753 }
754 }
755 if (in_run) {
756 in_run = 0;
757 int nw_cell = ncols * row + start_run;
758 int ne_cell = ncols * row + ncols - 1;
759 int sw_cell = nw_cell + ncols;
760 int se_cell = ne_cell + ncols;
761 int nw_tri = simplification_tris[nw_cell];
762 int ne_tri = simplification_tris[ne_cell];
763 int sw_tri = simplification_tris[sw_cell];
764 int se_tri = simplification_tris[se_cell];
765 int nw_corner = nw_tri * 3 + 4;
766 int ne_corner = ne_tri * 3 + 0;
767 int sw_corner = sw_tri * 3 + 2;
768 int se_corner = se_tri * 3 + 1;
769 *pnewtris++ = tris[se_corner];
770 *pnewtris++ = tris[sw_corner];
771 *pnewtris++ = tris[nw_corner];
772 *pnewtris++ = tris[nw_corner];
773 *pnewtris++ = tris[ne_corner];
774 *pnewtris++ = tris[se_corner];
775 }
776 }
777 ptris = pnewtris;
778 ntris -= neliminated_triangles;
779 free(tris);
780 tris = newtris;
781 free(simplification_codes);
782 free(simplification_tris);
783 free(simplification_ntris);
784 }
785
786 // Append all extrusion triangles to the main triangle array.
787 // We need them to be last so that they form a contiguous sequence.
788 pconntris = conntris;
789 for (int i = 0; i < nconntris; i++) {
790 *ptris++ = *pconntris++;
791 *ptris++ = *pconntris++;
792 *ptris++ = *pconntris++;
793 ntris++;
794 }
795 free(conntris);
796
797 // Final cleanup and return.
798 assert(npts <= maxpts);
799 assert(ntris <= maxtris);
800 mesh->npoints = npts;
801 mesh->points = pts;
802 mesh->ntriangles = ntris;
803 mesh->triangles = tris;
804 mesh->nconntriangles = nconntris;
805 return mlist;
806 }
807
808 #undef MIN
809 #undef MAX
810 #undef CLAMP
811 #endif
812