git.sesse.net Git - pistorm/blob - raylib_pi4_test/external/par_shapes.h

   1 // SHAPES :: https://github.com/prideout/par
   2 // Simple C library for creation and manipulation of triangle meshes.
   3 //
   4 // The API is divided into three sections:
   5 //
   6 //   - Generators.  Create parametric surfaces, platonic solids, etc.
   7 //   - Queries.     Ask a mesh for its axis-aligned bounding box, etc.
   8 //   - Transforms.  Rotate a mesh, merge it with another, add normals, etc.
   9 //
  10 // In addition to the comment block above each function declaration, the API
  11 // has informal documentation here:
  12 //
  13 //     https://prideout.net/shapes
  14 //
  15 // For our purposes, a "mesh" is a list of points and a list of triangles; the
  16 // former is a flattened list of three-tuples (32-bit floats) and the latter is
  17 // also a flattened list of three-tuples (16-bit uints).  Triangles are always
  18 // oriented such that their front face winds counter-clockwise.
  19 //
  20 // Optionally, meshes can contain 3D normals (one per vertex), and 2D texture
  21 // coordinates (one per vertex).  That's it!  If you need something fancier,
  22 // look elsewhere.
  23 //
  24 // The MIT License
  25 // Copyright (c) 2015 Philip Rideout
  26
  27 #ifndef PAR_SHAPES_H
  28 #define PAR_SHAPES_H
  29
  30 #ifdef __cplusplus
  31 extern "C" {
  32 #endif
  33
  34 #include <stdint.h>
  35
  36 // Ray: commented to avoid conflict with raylib bool
  37 /*
  38 #if !defined(_MSC_VER)
  39 # include <stdbool.h>
  40 #else // MSVC
  41 # if _MSC_VER >= 1800
  42 #  include <stdbool.h>
  43 # else // stdbool.h missing prior to MSVC++ 12.0 (VS2013)
  44 //#  define bool int
  45 //#  define true 1
  46 //#  define false 0
  47 # endif
  48 #endif
  49 */
  50
  51 #ifndef PAR_SHAPES_T
  52 #define PAR_SHAPES_T uint16_t
  53 #endif
  54
  55 typedef struct par_shapes_mesh_s {
  56     float* points;           // Flat list of 3-tuples (X Y Z X Y Z...)
  57     int npoints;             // Number of points
  58     PAR_SHAPES_T* triangles; // Flat list of 3-tuples (I J K I J K...)
  59     int ntriangles;          // Number of triangles
  60     float* normals;          // Optional list of 3-tuples (X Y Z X Y Z...)
  61     float* tcoords;          // Optional list of 2-tuples (U V U V U V...)
  62 } par_shapes_mesh;
  63
  64 void par_shapes_free_mesh(par_shapes_mesh*);
  65
  66 // Generators ------------------------------------------------------------------
  67
  68 // Instance a cylinder that sits on the Z=0 plane using the given tessellation
  69 // levels across the UV domain.  Think of "slices" like a number of pizza
  70 // slices, and "stacks" like a number of stacked rings.  Height and radius are
  71 // both 1.0, but they can easily be changed with par_shapes_scale.
  72 par_shapes_mesh* par_shapes_create_cylinder(int slices, int stacks);
  73
  74 // Create a donut that sits on the Z=0 plane with the specified inner radius.
  75 // The outer radius can be controlled with par_shapes_scale.
  76 par_shapes_mesh* par_shapes_create_torus(int slices, int stacks, float radius);
  77
  78 // Create a sphere with texture coordinates and small triangles near the poles.
  79 par_shapes_mesh* par_shapes_create_parametric_sphere(int slices, int stacks);
  80
  81 // Approximate a sphere with a subdivided icosahedron, which produces a nice
  82 // distribution of triangles, but no texture coordinates.  Each subdivision
  83 // level scales the number of triangles by four, so use a very low number.
  84 par_shapes_mesh* par_shapes_create_subdivided_sphere(int nsubdivisions);
  85
  86 // More parametric surfaces.
  87 par_shapes_mesh* par_shapes_create_klein_bottle(int slices, int stacks);
  88 par_shapes_mesh* par_shapes_create_trefoil_knot(int slices, int stacks,
  89     float radius);
  90 par_shapes_mesh* par_shapes_create_hemisphere(int slices, int stacks);
  91 par_shapes_mesh* par_shapes_create_plane(int slices, int stacks);
  92
  93 // Create a parametric surface from a callback function that consumes a 2D
  94 // point in [0,1] and produces a 3D point.
  95 typedef void (*par_shapes_fn)(float const*, float*, void*);
  96 par_shapes_mesh* par_shapes_create_parametric(par_shapes_fn, int slices,
  97     int stacks, void* userdata);
  98
  99 // Generate points for a 20-sided polyhedron that fits in the unit sphere.
 100 // Texture coordinates and normals are not generated.
 101 par_shapes_mesh* par_shapes_create_icosahedron();
 102
 103 // Generate points for a 12-sided polyhedron that fits in the unit sphere.
 104 // Again, texture coordinates and normals are not generated.
 105 par_shapes_mesh* par_shapes_create_dodecahedron();
 106
 107 // More platonic solids.
 108 par_shapes_mesh* par_shapes_create_octahedron();
 109 par_shapes_mesh* par_shapes_create_tetrahedron();
 110 par_shapes_mesh* par_shapes_create_cube();
 111
 112 // Generate an orientable disk shape in 3-space.  Does not include normals or
 113 // texture coordinates.
 114 par_shapes_mesh* par_shapes_create_disk(float radius, int slices,
 115     float const* center, float const* normal);
 116
 117 // Create an empty shape.  Useful for building scenes with merge_and_free.
 118 par_shapes_mesh* par_shapes_create_empty();
 119
 120 // Generate a rock shape that sits on the Y=0 plane, and sinks into it a bit.
 121 // This includes smooth normals but no texture coordinates.  Each subdivision
 122 // level scales the number of triangles by four, so use a very low number.
 123 par_shapes_mesh* par_shapes_create_rock(int seed, int nsubdivisions);
 124
 125 // Create trees or vegetation by executing a recursive turtle graphics program.
 126 // The program is a list of command-argument pairs.  See the unit test for
 127 // an example.  Texture coordinates and normals are not generated.
 128 par_shapes_mesh* par_shapes_create_lsystem(char const* program, int slices,
 129     int maxdepth);
 130
 131 // Queries ---------------------------------------------------------------------
 132
 133 // Dump out a text file conforming to the venerable OBJ format.
 134 void par_shapes_export(par_shapes_mesh const*, char const* objfile);
 135
 136 // Take a pointer to 6 floats and set them to min xyz, max xyz.
 137 void par_shapes_compute_aabb(par_shapes_mesh const* mesh, float* aabb);
 138
 139 // Make a deep copy of a mesh.  To make a brand new copy, pass null to "target".
 140 // To avoid memory churn, pass an existing mesh to "target".
 141 par_shapes_mesh* par_shapes_clone(par_shapes_mesh const* mesh,
 142     par_shapes_mesh* target);
 143
 144 // Transformations -------------------------------------------------------------
 145
 146 void par_shapes_merge(par_shapes_mesh* dst, par_shapes_mesh const* src);
 147 void par_shapes_translate(par_shapes_mesh*, float x, float y, float z);
 148 void par_shapes_rotate(par_shapes_mesh*, float radians, float const* axis);
 149 void par_shapes_scale(par_shapes_mesh*, float x, float y, float z);
 150 void par_shapes_merge_and_free(par_shapes_mesh* dst, par_shapes_mesh* src);
 151
 152 // Reverse the winding of a run of faces.  Useful when drawing the inside of
 153 // a Cornell Box.  Pass 0 for nfaces to reverse every face in the mesh.
 154 void par_shapes_invert(par_shapes_mesh*, int startface, int nfaces);
 155
 156 // Remove all triangles whose area is less than minarea.
 157 void par_shapes_remove_degenerate(par_shapes_mesh*, float minarea);
 158
 159 // Dereference the entire index buffer and replace the point list.
 160 // This creates an inefficient structure, but is useful for drawing facets.
 161 // If create_indices is true, a trivial "0 1 2 3..." index buffer is generated.
 162 void par_shapes_unweld(par_shapes_mesh* mesh, bool create_indices);
 163
 164 // Merge colocated verts, build a new index buffer, and return the
 165 // optimized mesh.  Epsilon is the maximum distance to consider when
 166 // welding vertices. The mapping argument can be null, or a pointer to
 167 // npoints integers, which gets filled with the mapping from old vertex
 168 // indices to new indices.
 169 par_shapes_mesh* par_shapes_weld(par_shapes_mesh const*, float epsilon,
 170     PAR_SHAPES_T* mapping);
 171
 172 // Compute smooth normals by averaging adjacent facet normals.
 173 void par_shapes_compute_normals(par_shapes_mesh* m);
 174
 175 #ifndef PAR_PI
 176 #define PAR_PI (3.14159265359)
 177 #define PAR_MIN(a, b) (a > b ? b : a)
 178 #define PAR_MAX(a, b) (a > b ? a : b)
 179 #define PAR_CLAMP(v, lo, hi) PAR_MAX(lo, PAR_MIN(hi, v))
 180 #define PAR_SWAP(T, A, B) { T tmp = B; B = A; A = tmp; }
 181 #define PAR_SQR(a) ((a) * (a))
 182 #endif
 183
 184 #ifndef PAR_MALLOC
 185 #define PAR_MALLOC(T, N) ((T*) malloc(N * sizeof(T)))
 186 #define PAR_CALLOC(T, N) ((T*) calloc(N * sizeof(T), 1))
 187 #define PAR_REALLOC(T, BUF, N) ((T*) realloc(BUF, sizeof(T) * (N)))
 188 #define PAR_FREE(BUF) free(BUF)
 189 #endif
 190
 191 #ifdef __cplusplus
 192 }
 193 #endif
 194
 195 // -----------------------------------------------------------------------------
 196 // END PUBLIC API
 197 // -----------------------------------------------------------------------------
 198
 199 #ifdef PAR_SHAPES_IMPLEMENTATION
 200 #include <stdlib.h>
 201 #include <stdio.h>
 202 #include <assert.h>
 203 #include <float.h>
 204 #include <string.h>
 205 #include <math.h>
 206 #include <errno.h>
 207
 208 static void par_shapes__sphere(float const* uv, float* xyz, void*);
 209 static void par_shapes__hemisphere(float const* uv, float* xyz, void*);
 210 static void par_shapes__plane(float const* uv, float* xyz, void*);
 211 static void par_shapes__klein(float const* uv, float* xyz, void*);
 212 static void par_shapes__cylinder(float const* uv, float* xyz, void*);
 213 static void par_shapes__torus(float const* uv, float* xyz, void*);
 214 static void par_shapes__trefoil(float const* uv, float* xyz, void*);
 215
 216 struct osn_context;
 217 static int par__simplex_noise(int64_t seed, struct osn_context** ctx);
 218 static void par__simplex_noise_free(struct osn_context* ctx);
 219 static double par__simplex_noise2(struct osn_context* ctx, double x, double y);
 220
 221 static void par_shapes__copy3(float* result, float const* a)
 222 {
 223     result[0] = a[0];
 224     result[1] = a[1];
 225     result[2] = a[2];
 226 }
 227
 228 static float par_shapes__dot3(float const* a, float const* b)
 229 {
 230     return b[0] * a[0] + b[1] * a[1] + b[2] * a[2];
 231 }
 232
 233 static void par_shapes__transform3(float* p, float const* x, float const* y,
 234     float const* z)
 235 {
 236     float px = par_shapes__dot3(p, x);
 237     float py = par_shapes__dot3(p, y);
 238     float pz = par_shapes__dot3(p, z);
 239     p[0] = px;
 240     p[1] = py;
 241     p[2] = pz;
 242 }
 243
 244 static void par_shapes__cross3(float* result, float const* a, float const* b)
 245 {
 246     float x = (a[1] * b[2]) - (a[2] * b[1]);
 247     float y = (a[2] * b[0]) - (a[0] * b[2]);
 248     float z = (a[0] * b[1]) - (a[1] * b[0]);
 249     result[0] = x;
 250     result[1] = y;
 251     result[2] = z;
 252 }
 253
 254 static void par_shapes__mix3(float* d, float const* a, float const* b, float t)
 255 {
 256     float x = b[0] * t + a[0] * (1 - t);
 257     float y = b[1] * t + a[1] * (1 - t);
 258     float z = b[2] * t + a[2] * (1 - t);
 259     d[0] = x;
 260     d[1] = y;
 261     d[2] = z;
 262 }
 263
 264 static void par_shapes__scale3(float* result, float a)
 265 {
 266     result[0] *= a;
 267     result[1] *= a;
 268     result[2] *= a;
 269 }
 270
 271 static void par_shapes__normalize3(float* v)
 272 {
 273     float lsqr = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
 274     if (lsqr > 0) {
 275         par_shapes__scale3(v, 1.0f / lsqr);
 276     }
 277 }
 278
 279 static void par_shapes__subtract3(float* result, float const* a)
 280 {
 281     result[0] -= a[0];
 282     result[1] -= a[1];
 283     result[2] -= a[2];
 284 }
 285
 286 static void par_shapes__add3(float* result, float const* a)
 287 {
 288     result[0] += a[0];
 289     result[1] += a[1];
 290     result[2] += a[2];
 291 }
 292
 293 static float par_shapes__sqrdist3(float const* a, float const* b)
 294 {
 295     float dx = a[0] - b[0];
 296     float dy = a[1] - b[1];
 297     float dz = a[2] - b[2];
 298     return dx * dx + dy * dy + dz * dz;
 299 }
 300
 301 static void par_shapes__compute_welded_normals(par_shapes_mesh* m)
 302 {
 303     m->normals = PAR_MALLOC(float, m->npoints * 3);
 304     PAR_SHAPES_T* weldmap = PAR_MALLOC(PAR_SHAPES_T, m->npoints);
 305     par_shapes_mesh* welded = par_shapes_weld(m, 0.01, weldmap);
 306     par_shapes_compute_normals(welded);
 307     float* pdst = m->normals;
 308     for (int i = 0; i < m->npoints; i++, pdst += 3) {
 309         int d = weldmap[i];
 310         float const* pnormal = welded->normals + d * 3;
 311         pdst[0] = pnormal[0];
 312         pdst[1] = pnormal[1];
 313         pdst[2] = pnormal[2];
 314     }
 315     PAR_FREE(weldmap);
 316     par_shapes_free_mesh(welded);
 317 }
 318
 319 par_shapes_mesh* par_shapes_create_cylinder(int slices, int stacks)
 320 {
 321     if (slices < 3 || stacks < 1) {
 322         return 0;
 323     }
 324     return par_shapes_create_parametric(par_shapes__cylinder, slices,
 325         stacks, 0);
 326 }
 327
 328 par_shapes_mesh* par_shapes_create_parametric_sphere(int slices, int stacks)
 329 {
 330     if (slices < 3 || stacks < 3) {
 331         return 0;
 332     }
 333     par_shapes_mesh* m = par_shapes_create_parametric(par_shapes__sphere,
 334         slices, stacks, 0);
 335     par_shapes_remove_degenerate(m, 0.0001);
 336     return m;
 337 }
 338
 339 par_shapes_mesh* par_shapes_create_hemisphere(int slices, int stacks)
 340 {
 341     if (slices < 3 || stacks < 3) {
 342         return 0;
 343     }
 344     par_shapes_mesh* m = par_shapes_create_parametric(par_shapes__hemisphere,
 345         slices, stacks, 0);
 346     par_shapes_remove_degenerate(m, 0.0001);
 347     return m;
 348 }
 349
 350 par_shapes_mesh* par_shapes_create_torus(int slices, int stacks, float radius)
 351 {
 352     if (slices < 3 || stacks < 3) {
 353         return 0;
 354     }
 355     assert(radius <= 1.0 && "Use smaller radius to avoid self-intersection.");
 356     assert(radius >= 0.1 && "Use larger radius to avoid self-intersection.");
 357     void* userdata = (void*) &radius;
 358     return par_shapes_create_parametric(par_shapes__torus, slices,
 359         stacks, userdata);
 360 }
 361
 362 par_shapes_mesh* par_shapes_create_klein_bottle(int slices, int stacks)
 363 {
 364     if (slices < 3 || stacks < 3) {
 365         return 0;
 366     }
 367     par_shapes_mesh* mesh = par_shapes_create_parametric(
 368         par_shapes__klein, slices, stacks, 0);
 369     int face = 0;
 370     for (int stack = 0; stack < stacks; stack++) {
 371         for (int slice = 0; slice < slices; slice++, face += 2) {
 372             if (stack < 27 * stacks / 32) {
 373                 par_shapes_invert(mesh, face, 2);
 374             }
 375         }
 376     }
 377     par_shapes__compute_welded_normals(mesh);
 378     return mesh;
 379 }
 380
 381 par_shapes_mesh* par_shapes_create_trefoil_knot(int slices, int stacks,
 382     float radius)
 383 {
 384     if (slices < 3 || stacks < 3) {
 385         return 0;
 386     }
 387     assert(radius <= 3.0 && "Use smaller radius to avoid self-intersection.");
 388     assert(radius >= 0.5 && "Use larger radius to avoid self-intersection.");
 389     void* userdata = (void*) &radius;
 390     return par_shapes_create_parametric(par_shapes__trefoil, slices,
 391         stacks, userdata);
 392 }
 393
 394 par_shapes_mesh* par_shapes_create_plane(int slices, int stacks)
 395 {
 396     if (slices < 1 || stacks < 1) {
 397         return 0;
 398     }
 399     return par_shapes_create_parametric(par_shapes__plane, slices,
 400         stacks, 0);
 401 }
 402
 403 par_shapes_mesh* par_shapes_create_parametric(par_shapes_fn fn,
 404     int slices, int stacks, void* userdata)
 405 {
 406     par_shapes_mesh* mesh = PAR_CALLOC(par_shapes_mesh, 1);
 407
 408     // Generate verts.
 409     mesh->npoints = (slices + 1) * (stacks + 1);
 410     mesh->points = PAR_CALLOC(float, 3 * mesh->npoints);
 411     float uv[2];
 412     float xyz[3];
 413     float* points = mesh->points;
 414     for (int stack = 0; stack < stacks + 1; stack++) {
 415         uv[0] = (float) stack / stacks;
 416         for (int slice = 0; slice < slices + 1; slice++) {
 417             uv[1] = (float) slice / slices;
 418             fn(uv, xyz, userdata);
 419             *points++ = xyz[0];
 420             *points++ = xyz[1];
 421             *points++ = xyz[2];
 422         }
 423     }
 424
 425     // Generate texture coordinates.
 426     mesh->tcoords = PAR_CALLOC(float, 2 * mesh->npoints);
 427     float* uvs = mesh->tcoords;
 428     for (int stack = 0; stack < stacks + 1; stack++) {
 429         uv[0] = (float) stack / stacks;
 430         for (int slice = 0; slice < slices + 1; slice++) {
 431             uv[1] = (float) slice / slices;
 432             *uvs++ = uv[0];
 433             *uvs++ = uv[1];
 434         }
 435     }
 436
 437     // Generate faces.
 438     mesh->ntriangles = 2 * slices * stacks;
 439     mesh->triangles = PAR_CALLOC(PAR_SHAPES_T, 3 * mesh->ntriangles);
 440     int v = 0;
 441     PAR_SHAPES_T* face = mesh->triangles;
 442     for (int stack = 0; stack < stacks; stack++) {
 443         for (int slice = 0; slice < slices; slice++) {
 444             int next = slice + 1;
 445             *face++ = v + slice + slices + 1;
 446             *face++ = v + next;
 447             *face++ = v + slice;
 448             *face++ = v + slice + slices + 1;
 449             *face++ = v + next + slices + 1;
 450             *face++ = v + next;
 451         }
 452         v += slices + 1;
 453     }
 454
 455     par_shapes__compute_welded_normals(mesh);
 456     return mesh;
 457 }
 458
 459 void par_shapes_free_mesh(par_shapes_mesh* mesh)
 460 {
 461     PAR_FREE(mesh->points);
 462     PAR_FREE(mesh->triangles);
 463     PAR_FREE(mesh->normals);
 464     PAR_FREE(mesh->tcoords);
 465     PAR_FREE(mesh);
 466 }
 467
 468 void par_shapes_export(par_shapes_mesh const* mesh, char const* filename)
 469 {
 470     FILE* objfile = fopen(filename, "wt");
 471     float const* points = mesh->points;
 472     float const* tcoords = mesh->tcoords;
 473     float const* norms = mesh->normals;
 474     PAR_SHAPES_T const* indices = mesh->triangles;
 475     if (tcoords && norms) {
 476         for (int nvert = 0; nvert < mesh->npoints; nvert++) {
 477             fprintf(objfile, "v %f %f %f\n", points[0], points[1], points[2]);
 478             fprintf(objfile, "vt %f %f\n", tcoords[0], tcoords[1]);
 479             fprintf(objfile, "vn %f %f %f\n", norms[0], norms[1], norms[2]);
 480             points += 3;
 481             norms += 3;
 482             tcoords += 2;
 483         }
 484         for (int nface = 0; nface < mesh->ntriangles; nface++) {
 485             int a = 1 + *indices++;
 486             int b = 1 + *indices++;
 487             int c = 1 + *indices++;
 488             fprintf(objfile, "f %d/%d/%d %d/%d/%d %d/%d/%d\n",
 489                 a, a, a, b, b, b, c, c, c);
 490         }
 491     } else if (norms) {
 492         for (int nvert = 0; nvert < mesh->npoints; nvert++) {
 493             fprintf(objfile, "v %f %f %f\n", points[0], points[1], points[2]);
 494             fprintf(objfile, "vn %f %f %f\n", norms[0], norms[1], norms[2]);
 495             points += 3;
 496             norms += 3;
 497         }
 498         for (int nface = 0; nface < mesh->ntriangles; nface++) {
 499             int a = 1 + *indices++;
 500             int b = 1 + *indices++;
 501             int c = 1 + *indices++;
 502             fprintf(objfile, "f %d//%d %d//%d %d//%d\n", a, a, b, b, c, c);
 503         }
 504     } else if (tcoords) {
 505         for (int nvert = 0; nvert < mesh->npoints; nvert++) {
 506             fprintf(objfile, "v %f %f %f\n", points[0], points[1], points[2]);
 507             fprintf(objfile, "vt %f %f\n", tcoords[0], tcoords[1]);
 508             points += 3;
 509             tcoords += 2;
 510         }
 511         for (int nface = 0; nface < mesh->ntriangles; nface++) {
 512             int a = 1 + *indices++;
 513             int b = 1 + *indices++;
 514             int c = 1 + *indices++;
 515             fprintf(objfile, "f %d/%d %d/%d %d/%d\n", a, a, b, b, c, c);
 516         }
 517     } else {
 518         for (int nvert = 0; nvert < mesh->npoints; nvert++) {
 519             fprintf(objfile, "v %f %f %f\n", points[0], points[1], points[2]);
 520             points += 3;
 521         }
 522         for (int nface = 0; nface < mesh->ntriangles; nface++) {
 523             int a = 1 + *indices++;
 524             int b = 1 + *indices++;
 525             int c = 1 + *indices++;
 526             fprintf(objfile, "f %d %d %d\n", a, b, c);
 527         }
 528     }
 529     fclose(objfile);
 530 }
 531
 532 static void par_shapes__sphere(float const* uv, float* xyz, void* userdata)
 533 {
 534     float phi = uv[0] * PAR_PI;
 535     float theta = uv[1] * 2 * PAR_PI;
 536     xyz[0] = cosf(theta) * sinf(phi);
 537     xyz[1] = sinf(theta) * sinf(phi);
 538     xyz[2] = cosf(phi);
 539 }
 540
 541 static void par_shapes__hemisphere(float const* uv, float* xyz, void* userdata)
 542 {
 543     float phi = uv[0] * PAR_PI;
 544     float theta = uv[1] * PAR_PI;
 545     xyz[0] = cosf(theta) * sinf(phi);
 546     xyz[1] = sinf(theta) * sinf(phi);
 547     xyz[2] = cosf(phi);
 548 }
 549
 550 static void par_shapes__plane(float const* uv, float* xyz, void* userdata)
 551 {
 552     xyz[0] = uv[0];
 553     xyz[1] = uv[1];
 554     xyz[2] = 0;
 555 }
 556
 557 static void par_shapes__klein(float const* uv, float* xyz, void* userdata)
 558 {
 559     float u = uv[0] * PAR_PI;
 560     float v = uv[1] * 2 * PAR_PI;
 561     u = u * 2;
 562     if (u < PAR_PI) {
 563         xyz[0] = 3 * cosf(u) * (1 + sinf(u)) + (2 * (1 - cosf(u) / 2)) *
 564             cosf(u) * cosf(v);
 565         xyz[2] = -8 * sinf(u) - 2 * (1 - cosf(u) / 2) * sinf(u) * cosf(v);
 566     } else {
 567         xyz[0] = 3 * cosf(u) * (1 + sinf(u)) + (2 * (1 - cosf(u) / 2)) *
 568             cosf(v + PAR_PI);
 569         xyz[2] = -8 * sinf(u);
 570     }
 571     xyz[1] = -2 * (1 - cosf(u) / 2) * sinf(v);
 572 }
 573
 574 static void par_shapes__cylinder(float const* uv, float* xyz, void* userdata)
 575 {
 576     float theta = uv[1] * 2 * PAR_PI;
 577     xyz[0] = sinf(theta);
 578     xyz[1] = cosf(theta);
 579     xyz[2] = uv[0];
 580 }
 581
 582 static void par_shapes__torus(float const* uv, float* xyz, void* userdata)
 583 {
 584     float major = 1;
 585     float minor = *((float*) userdata);
 586     float theta = uv[0] * 2 * PAR_PI;
 587     float phi = uv[1] * 2 * PAR_PI;
 588     float beta = major + minor * cosf(phi);
 589     xyz[0] = cosf(theta) * beta;
 590     xyz[1] = sinf(theta) * beta;
 591     xyz[2] = sinf(phi) * minor;
 592 }
 593
 594 static void par_shapes__trefoil(float const* uv, float* xyz, void* userdata)
 595 {
 596     float minor = *((float*) userdata);
 597     const float a = 0.5f;
 598     const float b = 0.3f;
 599     const float c = 0.5f;
 600     const float d = minor * 0.1f;
 601     const float u = (1 - uv[0]) * 4 * PAR_PI;
 602     const float v = uv[1] * 2 * PAR_PI;
 603     const float r = a + b * cos(1.5f * u);
 604     const float x = r * cos(u);
 605     const float y = r * sin(u);
 606     const float z = c * sin(1.5f * u);
 607     float q[3];
 608     q[0] =
 609         -1.5f * b * sin(1.5f * u) * cos(u) - (a + b * cos(1.5f * u)) * sin(u);
 610     q[1] =
 611         -1.5f * b * sin(1.5f * u) * sin(u) + (a + b * cos(1.5f * u)) * cos(u);
 612     q[2] = 1.5f * c * cos(1.5f * u);
 613     par_shapes__normalize3(q);
 614     float qvn[3] = {q[1], -q[0], 0};
 615     par_shapes__normalize3(qvn);
 616     float ww[3];
 617     par_shapes__cross3(ww, q, qvn);
 618     xyz[0] = x + d * (qvn[0] * cos(v) + ww[0] * sin(v));
 619     xyz[1] = y + d * (qvn[1] * cos(v) + ww[1] * sin(v));
 620     xyz[2] = z + d * ww[2] * sin(v);
 621 }
 622
 623 void par_shapes_merge(par_shapes_mesh* dst, par_shapes_mesh const* src)
 624 {
 625     PAR_SHAPES_T offset = dst->npoints;
 626     int npoints = dst->npoints + src->npoints;
 627     int vecsize = sizeof(float) * 3;
 628     dst->points = PAR_REALLOC(float, dst->points, 3 * npoints);
 629     memcpy(dst->points + 3 * dst->npoints, src->points, vecsize * src->npoints);
 630     dst->npoints = npoints;
 631     if (src->normals || dst->normals) {
 632         dst->normals = PAR_REALLOC(float, dst->normals, 3 * npoints);
 633         if (src->normals) {
 634             memcpy(dst->normals + 3 * offset, src->normals,
 635                 vecsize * src->npoints);
 636         }
 637     }
 638     if (src->tcoords || dst->tcoords) {
 639         int uvsize = sizeof(float) * 2;
 640         dst->tcoords = PAR_REALLOC(float, dst->tcoords, 2 * npoints);
 641         if (src->tcoords) {
 642             memcpy(dst->tcoords + 2 * offset, src->tcoords,
 643                 uvsize * src->npoints);
 644         }
 645     }
 646     int ntriangles = dst->ntriangles + src->ntriangles;
 647     dst->triangles = PAR_REALLOC(PAR_SHAPES_T, dst->triangles, 3 * ntriangles);
 648     PAR_SHAPES_T* ptriangles = dst->triangles + 3 * dst->ntriangles;
 649     PAR_SHAPES_T const* striangles = src->triangles;
 650     for (int i = 0; i < src->ntriangles; i++) {
 651         *ptriangles++ = offset + *striangles++;
 652         *ptriangles++ = offset + *striangles++;
 653         *ptriangles++ = offset + *striangles++;
 654     }
 655     dst->ntriangles = ntriangles;
 656 }
 657
 658 par_shapes_mesh* par_shapes_create_disk(float radius, int slices,
 659     float const* center, float const* normal)
 660 {
 661     par_shapes_mesh* mesh = PAR_CALLOC(par_shapes_mesh, 1);
 662     mesh->npoints = slices + 1;
 663     mesh->points = PAR_MALLOC(float, 3 * mesh->npoints);
 664     float* points = mesh->points;
 665     *points++ = 0;
 666     *points++ = 0;
 667     *points++ = 0;
 668     for (int i = 0; i < slices; i++) {
 669         float theta = i * PAR_PI * 2 / slices;
 670         *points++ = radius * cos(theta);
 671         *points++ = radius * sin(theta);
 672         *points++ = 0;
 673     }
 674     float nnormal[3] = {normal[0], normal[1], normal[2]};
 675     par_shapes__normalize3(nnormal);
 676     mesh->normals = PAR_MALLOC(float, 3 * mesh->npoints);
 677     float* norms = mesh->normals;
 678     for (int i = 0; i < mesh->npoints; i++) {
 679         *norms++ = nnormal[0];
 680         *norms++ = nnormal[1];
 681         *norms++ = nnormal[2];
 682     }
 683     mesh->ntriangles = slices;
 684     mesh->triangles = PAR_MALLOC(PAR_SHAPES_T, 3 * mesh->ntriangles);
 685     PAR_SHAPES_T* triangles = mesh->triangles;
 686     for (int i = 0; i < slices; i++) {
 687         *triangles++ = 0;
 688         *triangles++ = 1 + i;
 689         *triangles++ = 1 + (i + 1) % slices;
 690     }
 691     float k[3] = {0, 0, -1};
 692     float axis[3];
 693     par_shapes__cross3(axis, nnormal, k);
 694     par_shapes__normalize3(axis);
 695     par_shapes_rotate(mesh, acos(nnormal[2]), axis);
 696     par_shapes_translate(mesh, center[0], center[1], center[2]);
 697     return mesh;
 698 }
 699
 700 par_shapes_mesh* par_shapes_create_empty()
 701 {
 702     return PAR_CALLOC(par_shapes_mesh, 1);
 703 }
 704
 705 void par_shapes_translate(par_shapes_mesh* m, float x, float y, float z)
 706 {
 707     float* points = m->points;
 708     for (int i = 0; i < m->npoints; i++) {
 709         *points++ += x;
 710         *points++ += y;
 711         *points++ += z;
 712     }
 713 }
 714
 715 void par_shapes_rotate(par_shapes_mesh* mesh, float radians, float const* axis)
 716 {
 717     float s = sinf(radians);
 718     float c = cosf(radians);
 719     float x = axis[0];
 720     float y = axis[1];
 721     float z = axis[2];
 722     float xy = x * y;
 723     float yz = y * z;
 724     float zx = z * x;
 725     float oneMinusC = 1.0f - c;
 726     float col0[3] = {
 727         (((x * x) * oneMinusC) + c),
 728         ((xy * oneMinusC) + (z * s)), ((zx * oneMinusC) - (y * s))
 729     };
 730     float col1[3] = {
 731         ((xy * oneMinusC) - (z * s)),
 732         (((y * y) * oneMinusC) + c), ((yz * oneMinusC) + (x * s))
 733     };
 734     float col2[3] = {
 735         ((zx * oneMinusC) + (y * s)),
 736         ((yz * oneMinusC) - (x * s)), (((z * z) * oneMinusC) + c)
 737     };
 738     float* p = mesh->points;
 739     for (int i = 0; i < mesh->npoints; i++, p += 3) {
 740         float x = col0[0] * p[0] + col1[0] * p[1] + col2[0] * p[2];
 741         float y = col0[1] * p[0] + col1[1] * p[1] + col2[1] * p[2];
 742         float z = col0[2] * p[0] + col1[2] * p[1] + col2[2] * p[2];
 743         p[0] = x;
 744         p[1] = y;
 745         p[2] = z;
 746     }
 747     p = mesh->normals;
 748     if (p) {
 749         for (int i = 0; i < mesh->npoints; i++, p += 3) {
 750             float x = col0[0] * p[0] + col1[0] * p[1] + col2[0] * p[2];
 751             float y = col0[1] * p[0] + col1[1] * p[1] + col2[1] * p[2];
 752             float z = col0[2] * p[0] + col1[2] * p[1] + col2[2] * p[2];
 753             p[0] = x;
 754             p[1] = y;
 755             p[2] = z;
 756         }
 757     }
 758 }
 759
 760 void par_shapes_scale(par_shapes_mesh* m, float x, float y, float z)
 761 {
 762     float* points = m->points;
 763     for (int i = 0; i < m->npoints; i++) {
 764         *points++ *= x;
 765         *points++ *= y;
 766         *points++ *= z;
 767     }
 768 }
 769
 770 void par_shapes_merge_and_free(par_shapes_mesh* dst, par_shapes_mesh* src)
 771 {
 772     par_shapes_merge(dst, src);
 773     par_shapes_free_mesh(src);
 774 }
 775
 776 void par_shapes_compute_aabb(par_shapes_mesh const* m, float* aabb)
 777 {
 778     float* points = m->points;
 779     aabb[0] = aabb[3] = points[0];
 780     aabb[1] = aabb[4] = points[1];
 781     aabb[2] = aabb[5] = points[2];
 782     points += 3;
 783     for (int i = 1; i < m->npoints; i++, points += 3) {
 784         aabb[0] = PAR_MIN(points[0], aabb[0]);
 785         aabb[1] = PAR_MIN(points[1], aabb[1]);
 786         aabb[2] = PAR_MIN(points[2], aabb[2]);
 787         aabb[3] = PAR_MAX(points[0], aabb[3]);
 788         aabb[4] = PAR_MAX(points[1], aabb[4]);
 789         aabb[5] = PAR_MAX(points[2], aabb[5]);
 790     }
 791 }
 792
 793 void par_shapes_invert(par_shapes_mesh* m, int face, int nfaces)
 794 {
 795     nfaces = nfaces ? nfaces : m->ntriangles;
 796     PAR_SHAPES_T* tri = m->triangles + face * 3;
 797     for (int i = 0; i < nfaces; i++) {
 798         PAR_SWAP(PAR_SHAPES_T, tri[0], tri[2]);
 799         tri += 3;
 800     }
 801 }
 802
 803 par_shapes_mesh* par_shapes_create_icosahedron()
 804 {
 805     static float verts[] = {
 806         0.000,  0.000,  1.000,
 807         0.894,  0.000,  0.447,
 808         0.276,  0.851,  0.447,
 809         -0.724,  0.526,  0.447,
 810         -0.724, -0.526,  0.447,
 811         0.276, -0.851,  0.447,
 812         0.724,  0.526, -0.447,
 813         -0.276,  0.851, -0.447,
 814         -0.894,  0.000, -0.447,
 815         -0.276, -0.851, -0.447,
 816         0.724, -0.526, -0.447,
 817         0.000,  0.000, -1.000
 818     };
 819     static PAR_SHAPES_T faces[] = {
 820         0,1,2,
 821         0,2,3,
 822         0,3,4,
 823         0,4,5,
 824         0,5,1,
 825         7,6,11,
 826         8,7,11,
 827         9,8,11,
 828         10,9,11,
 829         6,10,11,
 830         6,2,1,
 831         7,3,2,
 832         8,4,3,
 833         9,5,4,
 834         10,1,5,
 835         6,7,2,
 836         7,8,3,
 837         8,9,4,
 838         9,10,5,
 839         10,6,1
 840     };
 841     par_shapes_mesh* mesh = PAR_CALLOC(par_shapes_mesh, 1);
 842     mesh->npoints = sizeof(verts) / sizeof(verts[0]) / 3;
 843     mesh->points = PAR_MALLOC(float, sizeof(verts) / 4);
 844     memcpy(mesh->points, verts, sizeof(verts));
 845     mesh->ntriangles = sizeof(faces) / sizeof(faces[0]) / 3;
 846     mesh->triangles = PAR_MALLOC(PAR_SHAPES_T, sizeof(faces) / 2);
 847     memcpy(mesh->triangles, faces, sizeof(faces));
 848     return mesh;
 849 }
 850
 851 par_shapes_mesh* par_shapes_create_dodecahedron()
 852 {
 853     static float verts[20 * 3] = {
 854         0.607, 0.000, 0.795,
 855         0.188, 0.577, 0.795,
 856         -0.491, 0.357, 0.795,
 857         -0.491, -0.357, 0.795,
 858         0.188, -0.577, 0.795,
 859         0.982, 0.000, 0.188,
 860         0.304, 0.934, 0.188,
 861         -0.795, 0.577, 0.188,
 862         -0.795, -0.577, 0.188,
 863         0.304, -0.934, 0.188,
 864         0.795, 0.577, -0.188,
 865         -0.304, 0.934, -0.188,
 866         -0.982, 0.000, -0.188,
 867         -0.304, -0.934, -0.188,
 868         0.795, -0.577, -0.188,
 869         0.491, 0.357, -0.795,
 870         -0.188, 0.577, -0.795,
 871         -0.607, 0.000, -0.795,
 872         -0.188, -0.577, -0.795,
 873         0.491, -0.357, -0.795,
 874     };
 875     static PAR_SHAPES_T pentagons[12 * 5] = {
 876         0,1,2,3,4,
 877         5,10,6,1,0,
 878         6,11,7,2,1,
 879         7,12,8,3,2,
 880         8,13,9,4,3,
 881         9,14,5,0,4,
 882         15,16,11,6,10,
 883         16,17,12,7,11,
 884         17,18,13,8,12,
 885         18,19,14,9,13,
 886         19,15,10,5,14,
 887         19,18,17,16,15
 888     };
 889     int npentagons = sizeof(pentagons) / sizeof(pentagons[0]) / 5;
 890     par_shapes_mesh* mesh = PAR_CALLOC(par_shapes_mesh, 1);
 891     int ncorners = sizeof(verts) / sizeof(verts[0]) / 3;
 892     mesh->npoints = ncorners;
 893     mesh->points = PAR_MALLOC(float, mesh->npoints * 3);
 894     memcpy(mesh->points, verts, sizeof(verts));
 895     PAR_SHAPES_T const* pentagon = pentagons;
 896     mesh->ntriangles = npentagons * 3;
 897     mesh->triangles = PAR_MALLOC(PAR_SHAPES_T, mesh->ntriangles * 3);
 898     PAR_SHAPES_T* tris = mesh->triangles;
 899     for (int p = 0; p < npentagons; p++, pentagon += 5) {
 900         *tris++ = pentagon[0];
 901         *tris++ = pentagon[1];
 902         *tris++ = pentagon[2];
 903         *tris++ = pentagon[0];
 904         *tris++ = pentagon[2];
 905         *tris++ = pentagon[3];
 906         *tris++ = pentagon[0];
 907         *tris++ = pentagon[3];
 908         *tris++ = pentagon[4];
 909     }
 910     return mesh;
 911 }
 912
 913 par_shapes_mesh* par_shapes_create_octahedron()
 914 {
 915     static float verts[6 * 3] = {
 916         0.000, 0.000, 1.000,
 917         1.000, 0.000, 0.000,
 918         0.000, 1.000, 0.000,
 919         -1.000, 0.000, 0.000,
 920         0.000, -1.000, 0.000,
 921         0.000, 0.000, -1.000
 922     };
 923     static PAR_SHAPES_T triangles[8 * 3] = {
 924         0,1,2,
 925         0,2,3,
 926         0,3,4,
 927         0,4,1,
 928         2,1,5,
 929         3,2,5,
 930         4,3,5,
 931         1,4,5,
 932     };
 933     int ntris = sizeof(triangles) / sizeof(triangles[0]) / 3;
 934     par_shapes_mesh* mesh = PAR_CALLOC(par_shapes_mesh, 1);
 935     int ncorners = sizeof(verts) / sizeof(verts[0]) / 3;
 936     mesh->npoints = ncorners;
 937     mesh->points = PAR_MALLOC(float, mesh->npoints * 3);
 938     memcpy(mesh->points, verts, sizeof(verts));
 939     PAR_SHAPES_T const* triangle = triangles;
 940     mesh->ntriangles = ntris;
 941     mesh->triangles = PAR_MALLOC(PAR_SHAPES_T, mesh->ntriangles * 3);
 942     PAR_SHAPES_T* tris = mesh->triangles;
 943     for (int p = 0; p < ntris; p++) {
 944         *tris++ = *triangle++;
 945         *tris++ = *triangle++;
 946         *tris++ = *triangle++;
 947     }
 948     return mesh;
 949 }
 950
 951 par_shapes_mesh* par_shapes_create_tetrahedron()
 952 {
 953     static float verts[4 * 3] = {
 954         0.000, 1.333, 0,
 955         0.943, 0, 0,
 956         -0.471, 0, 0.816,
 957         -0.471, 0, -0.816,
 958     };
 959     static PAR_SHAPES_T triangles[4 * 3] = {
 960         2,1,0,
 961         3,2,0,
 962         1,3,0,
 963         1,2,3,
 964     };
 965     int ntris = sizeof(triangles) / sizeof(triangles[0]) / 3;
 966     par_shapes_mesh* mesh = PAR_CALLOC(par_shapes_mesh, 1);
 967     int ncorners = sizeof(verts) / sizeof(verts[0]) / 3;
 968     mesh->npoints = ncorners;
 969     mesh->points = PAR_MALLOC(float, mesh->npoints * 3);
 970     memcpy(mesh->points, verts, sizeof(verts));
 971     PAR_SHAPES_T const* triangle = triangles;
 972     mesh->ntriangles = ntris;
 973     mesh->triangles = PAR_MALLOC(PAR_SHAPES_T, mesh->ntriangles * 3);
 974     PAR_SHAPES_T* tris = mesh->triangles;
 975     for (int p = 0; p < ntris; p++) {
 976         *tris++ = *triangle++;
 977         *tris++ = *triangle++;
 978         *tris++ = *triangle++;
 979     }
 980     return mesh;
 981 }
 982
 983 par_shapes_mesh* par_shapes_create_cube()
 984 {
 985     static float verts[8 * 3] = {
 986         0, 0, 0, // 0
 987         0, 1, 0, // 1
 988         1, 1, 0, // 2
 989         1, 0, 0, // 3
 990         0, 0, 1, // 4
 991         0, 1, 1, // 5
 992         1, 1, 1, // 6
 993         1, 0, 1, // 7
 994     };
 995     static PAR_SHAPES_T quads[6 * 4] = {
 996         7,6,5,4, // front
 997         0,1,2,3, // back
 998         6,7,3,2, // right
 999         5,6,2,1, // top
1000         4,5,1,0, // left
1001         7,4,0,3, // bottom
1002     };
1003     int nquads = sizeof(quads) / sizeof(quads[0]) / 4;
1004     par_shapes_mesh* mesh = PAR_CALLOC(par_shapes_mesh, 1);
1005     int ncorners = sizeof(verts) / sizeof(verts[0]) / 3;
1006     mesh->npoints = ncorners;
1007     mesh->points = PAR_MALLOC(float, mesh->npoints * 3);
1008     memcpy(mesh->points, verts, sizeof(verts));
1009     PAR_SHAPES_T const* quad = quads;
1010     mesh->ntriangles = nquads * 2;
1011     mesh->triangles = PAR_MALLOC(PAR_SHAPES_T, mesh->ntriangles * 3);
1012     PAR_SHAPES_T* tris = mesh->triangles;
1013     for (int p = 0; p < nquads; p++, quad += 4) {
1014         *tris++ = quad[0];
1015         *tris++ = quad[1];
1016         *tris++ = quad[2];
1017         *tris++ = quad[2];
1018         *tris++ = quad[3];
1019         *tris++ = quad[0];
1020     }
1021     return mesh;
1022 }
1023
1024 typedef struct {
1025     char* cmd;
1026     char* arg;
1027 } par_shapes__command;
1028
1029 typedef struct {
1030     char const* name;
1031     int weight;
1032     int ncommands;
1033     par_shapes__command* commands;
1034 } par_shapes__rule;
1035
1036 typedef struct {
1037     int pc;
1038     float position[3];
1039     float scale[3];
1040     par_shapes_mesh* orientation;
1041     par_shapes__rule* rule;
1042 } par_shapes__stackframe;
1043
1044 static par_shapes__rule* par_shapes__pick_rule(const char* name,
1045     par_shapes__rule* rules, int nrules)
1046 {
1047     par_shapes__rule* rule = 0;
1048     int total = 0;
1049     for (int i = 0; i < nrules; i++) {
1050         rule = rules + i;
1051         if (!strcmp(rule->name, name)) {
1052             total += rule->weight;
1053         }
1054     }
1055     float r = (float) rand() / RAND_MAX;
1056     float t = 0;
1057     for (int i = 0; i < nrules; i++) {
1058         rule = rules + i;
1059         if (!strcmp(rule->name, name)) {
1060             t += (float) rule->weight / total;
1061             if (t >= r) {
1062                 return rule;
1063             }
1064         }
1065     }
1066     return rule;
1067 }
1068
1069 static par_shapes_mesh* par_shapes__create_turtle()
1070 {
1071     const float xaxis[] = {1, 0, 0};
1072     const float yaxis[] = {0, 1, 0};
1073     const float zaxis[] = {0, 0, 1};
1074     par_shapes_mesh* turtle = PAR_CALLOC(par_shapes_mesh, 1);
1075     turtle->npoints = 3;
1076     turtle->points = PAR_CALLOC(float, turtle->npoints * 3);
1077     par_shapes__copy3(turtle->points + 0, xaxis);
1078     par_shapes__copy3(turtle->points + 3, yaxis);
1079     par_shapes__copy3(turtle->points + 6, zaxis);
1080     return turtle;
1081 }
1082
1083 static par_shapes_mesh* par_shapes__apply_turtle(par_shapes_mesh* mesh,
1084     par_shapes_mesh* turtle, float const* pos, float const* scale)
1085 {
1086     par_shapes_mesh* m = par_shapes_clone(mesh, 0);
1087     for (int p = 0; p < m->npoints; p++) {
1088         float* pt = m->points + p * 3;
1089         pt[0] *= scale[0];
1090         pt[1] *= scale[1];
1091         pt[2] *= scale[2];
1092         par_shapes__transform3(pt,
1093             turtle->points + 0, turtle->points + 3, turtle->points + 6);
1094         pt[0] += pos[0];
1095         pt[1] += pos[1];
1096         pt[2] += pos[2];
1097     }
1098     return m;
1099 }
1100
1101 static void par_shapes__connect(par_shapes_mesh* scene,
1102     par_shapes_mesh* cylinder, int slices)
1103 {
1104     int stacks = 1;
1105     int npoints = (slices + 1) * (stacks + 1);
1106     assert(scene->npoints >= npoints && "Cannot connect to empty scene.");
1107
1108     // Create the new point list.
1109     npoints = scene->npoints + (slices + 1);
1110     float* points = PAR_MALLOC(float, npoints * 3);
1111     memcpy(points, scene->points, sizeof(float) * scene->npoints * 3);
1112     float* newpts = points + scene->npoints * 3;
1113     memcpy(newpts, cylinder->points + (slices + 1) * 3,
1114         sizeof(float) * (slices + 1) * 3);
1115     PAR_FREE(scene->points);
1116     scene->points = points;
1117
1118     // Create the new triangle list.
1119     int ntriangles = scene->ntriangles + 2 * slices * stacks;
1120     PAR_SHAPES_T* triangles = PAR_MALLOC(PAR_SHAPES_T, ntriangles * 3);
1121     memcpy(triangles, scene->triangles, 2 * scene->ntriangles * 3);
1122     int v = scene->npoints - (slices + 1);
1123     PAR_SHAPES_T* face = triangles + scene->ntriangles * 3;
1124     for (int stack = 0; stack < stacks; stack++) {
1125         for (int slice = 0; slice < slices; slice++) {
1126             int next = slice + 1;
1127             *face++ = v + slice + slices + 1;
1128             *face++ = v + next;
1129             *face++ = v + slice;
1130             *face++ = v + slice + slices + 1;
1131             *face++ = v + next + slices + 1;
1132             *face++ = v + next;
1133         }
1134         v += slices + 1;
1135     }
1136     PAR_FREE(scene->triangles);
1137     scene->triangles = triangles;
1138
1139     scene->npoints = npoints;
1140     scene->ntriangles = ntriangles;
1141 }
1142
1143 par_shapes_mesh* par_shapes_create_lsystem(char const* text, int slices,
1144     int maxdepth)
1145 {
1146     char* program;
1147     program = PAR_MALLOC(char, strlen(text) + 1);
1148
1149     // The first pass counts the number of rules and commands.
1150     strcpy(program, text);
1151     char *cmd = strtok(program, " ");
1152     int nrules = 1;
1153     int ncommands = 0;
1154     while (cmd) {
1155         char *arg = strtok(0, " ");
1156         if (!arg) {
1157             //puts("lsystem error: unexpected end of program.");
1158             break;
1159         }
1160         if (!strcmp(cmd, "rule")) {
1161             nrules++;
1162         } else {
1163             ncommands++;
1164         }
1165         cmd = strtok(0, " ");
1166     }
1167
1168     // Allocate space.
1169     par_shapes__rule* rules = PAR_MALLOC(par_shapes__rule, nrules);
1170     par_shapes__command* commands = PAR_MALLOC(par_shapes__command, ncommands);
1171
1172     // Initialize the entry rule.
1173     par_shapes__rule* current_rule = &rules[0];
1174     par_shapes__command* current_command = &commands[0];
1175     current_rule->name = "entry";
1176     current_rule->weight = 1;
1177     current_rule->ncommands = 0;
1178     current_rule->commands = current_command;
1179
1180     // The second pass fills in the structures.
1181     strcpy(program, text);
1182     cmd = strtok(program, " ");
1183     while (cmd) {
1184         char *arg = strtok(0, " ");
1185         if (!strcmp(cmd, "rule")) {
1186             current_rule++;
1187
1188             // Split the argument into a rule name and weight.
1189             char* dot = strchr(arg, '.');
1190             if (dot) {
1191                 current_rule->weight = atoi(dot + 1);
1192                 *dot = 0;
1193             } else {
1194                 current_rule->weight = 1;
1195             }
1196
1197             current_rule->name = arg;
1198             current_rule->ncommands = 0;
1199             current_rule->commands = current_command;
1200         } else {
1201             current_rule->ncommands++;
1202             current_command->cmd = cmd;
1203             current_command->arg = arg;
1204             current_command++;
1205         }
1206         cmd = strtok(0, " ");
1207     }
1208
1209     // For testing purposes, dump out the parsed program.
1210     #ifdef TEST_PARSE
1211     /*
1212     for (int i = 0; i < nrules; i++) {
1213         par_shapes__rule rule = rules[i];
1214         printf("rule %s.%d\n", rule.name, rule.weight);
1215         for (int c = 0; c < rule.ncommands; c++) {
1216             par_shapes__command cmd = rule.commands[c];
1217             printf("\t%s %s\n", cmd.cmd, cmd.arg);
1218         }
1219     }
1220     */
1221     #endif
1222
1223     // Instantiate the aggregated shape and the template shapes.
1224     par_shapes_mesh* scene = PAR_CALLOC(par_shapes_mesh, 1);
1225     par_shapes_mesh* tube = par_shapes_create_cylinder(slices, 1);
1226     par_shapes_mesh* turtle = par_shapes__create_turtle();
1227
1228     // We're not attempting to support texture coordinates and normals
1229     // with L-systems, so remove them from the template shape.
1230     PAR_FREE(tube->normals);
1231     PAR_FREE(tube->tcoords);
1232     tube->normals = 0;
1233     tube->tcoords = 0;
1234
1235     const float xaxis[] = {1, 0, 0};
1236     const float yaxis[] = {0, 1, 0};
1237     const float zaxis[] = {0, 0, 1};
1238     const float units[] = {1, 1, 1};
1239
1240     // Execute the L-system program until the stack size is 0.
1241     par_shapes__stackframe* stack =
1242         PAR_CALLOC(par_shapes__stackframe, maxdepth);
1243     int stackptr = 0;
1244     stack[0].orientation = turtle;
1245     stack[0].rule = &rules[0];
1246     par_shapes__copy3(stack[0].scale, units);
1247     while (stackptr >= 0) {
1248         par_shapes__stackframe* frame = &stack[stackptr];
1249         par_shapes__rule* rule = frame->rule;
1250         par_shapes_mesh* turtle = frame->orientation;
1251         float* position = frame->position;
1252         float* scale = frame->scale;
1253         if (frame->pc >= rule->ncommands) {
1254             par_shapes_free_mesh(turtle);
1255             stackptr--;
1256             continue;
1257         }
1258
1259         par_shapes__command* cmd = rule->commands + (frame->pc++);
1260         #ifdef DUMP_TRACE
1261         //printf("%5s %5s %5s:%d  %03d\n", cmd->cmd, cmd->arg, rule->name, frame->pc - 1, stackptr);
1262         #endif
1263
1264         float value;
1265         if (!strcmp(cmd->cmd, "shape")) {
1266             par_shapes_mesh* m = par_shapes__apply_turtle(tube, turtle,
1267                 position, scale);
1268             if (!strcmp(cmd->arg, "connect")) {
1269                 par_shapes__connect(scene, m, slices);
1270             } else {
1271                 par_shapes_merge(scene, m);
1272             }
1273             par_shapes_free_mesh(m);
1274         } else if (!strcmp(cmd->cmd, "call") && stackptr < maxdepth - 1) {
1275             rule = par_shapes__pick_rule(cmd->arg, rules, nrules);
1276             frame = &stack[++stackptr];
1277             frame->rule = rule;
1278             frame->orientation = par_shapes_clone(turtle, 0);
1279             frame->pc = 0;
1280             par_shapes__copy3(frame->scale, scale);
1281             par_shapes__copy3(frame->position, position);
1282             continue;
1283         } else {
1284             value = atof(cmd->arg);
1285             if (!strcmp(cmd->cmd, "rx")) {
1286                 par_shapes_rotate(turtle, value * PAR_PI / 180.0, xaxis);
1287             } else if (!strcmp(cmd->cmd, "ry")) {
1288                 par_shapes_rotate(turtle, value * PAR_PI / 180.0, yaxis);
1289             } else if (!strcmp(cmd->cmd, "rz")) {
1290                 par_shapes_rotate(turtle, value * PAR_PI / 180.0, zaxis);
1291             } else if (!strcmp(cmd->cmd, "tx")) {
1292                 float vec[3] = {value, 0, 0};
1293                 float t[3] = {
1294                     par_shapes__dot3(turtle->points + 0, vec),
1295                     par_shapes__dot3(turtle->points + 3, vec),
1296                     par_shapes__dot3(turtle->points + 6, vec)
1297                 };
1298                 par_shapes__add3(position, t);
1299             } else if (!strcmp(cmd->cmd, "ty")) {
1300                 float vec[3] = {0, value, 0};
1301                 float t[3] = {
1302                     par_shapes__dot3(turtle->points + 0, vec),
1303                     par_shapes__dot3(turtle->points + 3, vec),
1304                     par_shapes__dot3(turtle->points + 6, vec)
1305                 };
1306                 par_shapes__add3(position, t);
1307             } else if (!strcmp(cmd->cmd, "tz")) {
1308                 float vec[3] = {0, 0, value};
1309                 float t[3] = {
1310                     par_shapes__dot3(turtle->points + 0, vec),
1311                     par_shapes__dot3(turtle->points + 3, vec),
1312                     par_shapes__dot3(turtle->points + 6, vec)
1313                 };
1314                 par_shapes__add3(position, t);
1315             } else if (!strcmp(cmd->cmd, "sx")) {
1316                 scale[0] *= value;
1317             } else if (!strcmp(cmd->cmd, "sy")) {
1318                 scale[1] *= value;
1319             } else if (!strcmp(cmd->cmd, "sz")) {
1320                 scale[2] *= value;
1321             } else if (!strcmp(cmd->cmd, "sa")) {
1322                 scale[0] *= value;
1323                 scale[1] *= value;
1324                 scale[2] *= value;
1325             }
1326         }
1327     }
1328     PAR_FREE(stack);
1329     PAR_FREE(program);
1330     PAR_FREE(rules);
1331     PAR_FREE(commands);
1332     return scene;
1333 }
1334
1335 void par_shapes_unweld(par_shapes_mesh* mesh, bool create_indices)
1336 {
1337     int npoints = mesh->ntriangles * 3;
1338     float* points = PAR_MALLOC(float, 3 * npoints);
1339     float* dst = points;
1340     PAR_SHAPES_T const* index = mesh->triangles;
1341     for (int i = 0; i < npoints; i++) {
1342         float const* src = mesh->points + 3 * (*index++);
1343         *dst++ = src[0];
1344         *dst++ = src[1];
1345         *dst++ = src[2];
1346     }
1347     PAR_FREE(mesh->points);
1348     mesh->points = points;
1349     mesh->npoints = npoints;
1350     if (create_indices) {
1351         PAR_SHAPES_T* tris = PAR_MALLOC(PAR_SHAPES_T, 3 * mesh->ntriangles);
1352         PAR_SHAPES_T* index = tris;
1353         for (int i = 0; i < mesh->ntriangles * 3; i++) {
1354             *index++ = i;
1355         }
1356         PAR_FREE(mesh->triangles);
1357         mesh->triangles = tris;
1358     }
1359 }
1360
1361 void par_shapes_compute_normals(par_shapes_mesh* m)
1362 {
1363     PAR_FREE(m->normals);
1364     m->normals = PAR_CALLOC(float, m->npoints * 3);
1365     PAR_SHAPES_T const* triangle = m->triangles;
1366     float next[3], prev[3], cp[3];
1367     for (int f = 0; f < m->ntriangles; f++, triangle += 3) {
1368         float const* pa = m->points + 3 * triangle[0];
1369         float const* pb = m->points + 3 * triangle[1];
1370         float const* pc = m->points + 3 * triangle[2];
1371         par_shapes__copy3(next, pb);
1372         par_shapes__subtract3(next, pa);
1373         par_shapes__copy3(prev, pc);
1374         par_shapes__subtract3(prev, pa);
1375         par_shapes__cross3(cp, next, prev);
1376         par_shapes__add3(m->normals + 3 * triangle[0], cp);
1377         par_shapes__copy3(next, pc);
1378         par_shapes__subtract3(next, pb);
1379         par_shapes__copy3(prev, pa);
1380         par_shapes__subtract3(prev, pb);
1381         par_shapes__cross3(cp, next, prev);
1382         par_shapes__add3(m->normals + 3 * triangle[1], cp);
1383         par_shapes__copy3(next, pa);
1384         par_shapes__subtract3(next, pc);
1385         par_shapes__copy3(prev, pb);
1386         par_shapes__subtract3(prev, pc);
1387         par_shapes__cross3(cp, next, prev);
1388         par_shapes__add3(m->normals + 3 * triangle[2], cp);
1389     }
1390     float* normal = m->normals;
1391     for (int p = 0; p < m->npoints; p++, normal += 3) {
1392         par_shapes__normalize3(normal);
1393     }
1394 }
1395
1396 static void par_shapes__subdivide(par_shapes_mesh* mesh)
1397 {
1398     assert(mesh->npoints == mesh->ntriangles * 3 && "Must be unwelded.");
1399     int ntriangles = mesh->ntriangles * 4;
1400     int npoints = ntriangles * 3;
1401     float* points = PAR_CALLOC(float, npoints * 3);
1402     float* dpoint = points;
1403     float const* spoint = mesh->points;
1404     for (int t = 0; t < mesh->ntriangles; t++, spoint += 9, dpoint += 3) {
1405         float const* a = spoint;
1406         float const* b = spoint + 3;
1407         float const* c = spoint + 6;
1408         float const* p0 = dpoint;
1409         float const* p1 = dpoint + 3;
1410         float const* p2 = dpoint + 6;
1411         par_shapes__mix3(dpoint, a, b, 0.5);
1412         par_shapes__mix3(dpoint += 3, b, c, 0.5);
1413         par_shapes__mix3(dpoint += 3, a, c, 0.5);
1414         par_shapes__add3(dpoint += 3, a);
1415         par_shapes__add3(dpoint += 3, p0);
1416         par_shapes__add3(dpoint += 3, p2);
1417         par_shapes__add3(dpoint += 3, p0);
1418         par_shapes__add3(dpoint += 3, b);
1419         par_shapes__add3(dpoint += 3, p1);
1420         par_shapes__add3(dpoint += 3, p2);
1421         par_shapes__add3(dpoint += 3, p1);
1422         par_shapes__add3(dpoint += 3, c);
1423     }
1424     PAR_FREE(mesh->points);
1425     mesh->points = points;
1426     mesh->npoints = npoints;
1427     mesh->ntriangles = ntriangles;
1428 }
1429
1430 par_shapes_mesh* par_shapes_create_subdivided_sphere(int nsubd)
1431 {
1432     par_shapes_mesh* mesh = par_shapes_create_icosahedron();
1433     par_shapes_unweld(mesh, false);
1434     PAR_FREE(mesh->triangles);
1435     mesh->triangles = 0;
1436     while (nsubd--) {
1437         par_shapes__subdivide(mesh);
1438     }
1439     for (int i = 0; i < mesh->npoints; i++) {
1440         par_shapes__normalize3(mesh->points + i * 3);
1441     }
1442     mesh->triangles = PAR_MALLOC(PAR_SHAPES_T, 3 * mesh->ntriangles);
1443     for (int i = 0; i < mesh->ntriangles * 3; i++) {
1444         mesh->triangles[i] = i;
1445     }
1446     par_shapes_mesh* tmp = mesh;
1447     mesh = par_shapes_weld(mesh, 0.01, 0);
1448     par_shapes_free_mesh(tmp);
1449     par_shapes_compute_normals(mesh);
1450     return mesh;
1451 }
1452
1453 par_shapes_mesh* par_shapes_create_rock(int seed, int subd)
1454 {
1455     par_shapes_mesh* mesh = par_shapes_create_subdivided_sphere(subd);
1456     struct osn_context* ctx;
1457     par__simplex_noise(seed, &ctx);
1458     for (int p = 0; p < mesh->npoints; p++) {
1459         float* pt = mesh->points + p * 3;
1460         float a = 0.25, f = 1.0;
1461         double n = a * par__simplex_noise2(ctx, f * pt[0], f * pt[2]);
1462         a *= 0.5; f *= 2;
1463         n += a * par__simplex_noise2(ctx, f * pt[0], f * pt[2]);
1464         pt[0] *= 1 + 2 * n;
1465         pt[1] *= 1 + n;
1466         pt[2] *= 1 + 2 * n;
1467         if (pt[1] < 0) {
1468             pt[1] = -pow(-pt[1], 0.5) / 2;
1469         }
1470     }
1471     par__simplex_noise_free(ctx);
1472     par_shapes_compute_normals(mesh);
1473     return mesh;
1474 }
1475
1476 par_shapes_mesh* par_shapes_clone(par_shapes_mesh const* mesh,
1477     par_shapes_mesh* clone)
1478 {
1479     if (!clone) {
1480         clone = PAR_CALLOC(par_shapes_mesh, 1);
1481     }
1482     clone->npoints = mesh->npoints;
1483     clone->points = PAR_REALLOC(float, clone->points, 3 * clone->npoints);
1484     memcpy(clone->points, mesh->points, sizeof(float) * 3 * clone->npoints);
1485     clone->ntriangles = mesh->ntriangles;
1486     clone->triangles = PAR_REALLOC(PAR_SHAPES_T, clone->triangles, 3 *
1487         clone->ntriangles);
1488     memcpy(clone->triangles, mesh->triangles,
1489         sizeof(PAR_SHAPES_T) * 3 * clone->ntriangles);
1490     if (mesh->normals) {
1491         clone->normals = PAR_REALLOC(float, clone->normals, 3 * clone->npoints);
1492         memcpy(clone->normals, mesh->normals,
1493             sizeof(float) * 3 * clone->npoints);
1494     }
1495     if (mesh->tcoords) {
1496         clone->tcoords = PAR_REALLOC(float, clone->tcoords, 2 * clone->npoints);
1497         memcpy(clone->tcoords, mesh->tcoords,
1498             sizeof(float) * 2 * clone->npoints);
1499     }
1500     return clone;
1501 }
1502
1503 static struct {
1504     float const* points;
1505     int gridsize;
1506 } par_shapes__sort_context;
1507
1508 static int par_shapes__cmp1(const void *arg0, const void *arg1)
1509 {
1510     const int g = par_shapes__sort_context.gridsize;
1511
1512     // Convert arg0 into a flattened grid index.
1513     PAR_SHAPES_T d0 = *(const PAR_SHAPES_T*) arg0;
1514     float const* p0 = par_shapes__sort_context.points + d0 * 3;
1515     int i0 = (int) p0[0];
1516     int j0 = (int) p0[1];
1517     int k0 = (int) p0[2];
1518     int index0 = i0 + g * j0 + g * g * k0;
1519
1520     // Convert arg1 into a flattened grid index.
1521     PAR_SHAPES_T d1 = *(const PAR_SHAPES_T*) arg1;
1522     float const* p1 = par_shapes__sort_context.points + d1 * 3;
1523     int i1 = (int) p1[0];
1524     int j1 = (int) p1[1];
1525     int k1 = (int) p1[2];
1526     int index1 = i1 + g * j1 + g * g * k1;
1527
1528     // Return the ordering.
1529     if (index0 < index1) return -1;
1530     if (index0 > index1) return 1;
1531     return 0;
1532 }
1533
1534 static void par_shapes__sort_points(par_shapes_mesh* mesh, int gridsize,
1535     PAR_SHAPES_T* sortmap)
1536 {
1537     // Run qsort over a list of consecutive integers that get deferenced
1538     // within the comparator function; this creates a reorder mapping.
1539     for (int i = 0; i < mesh->npoints; i++) {
1540         sortmap[i] = i;
1541     }
1542     par_shapes__sort_context.gridsize = gridsize;
1543     par_shapes__sort_context.points = mesh->points;
1544     qsort(sortmap, mesh->npoints, sizeof(PAR_SHAPES_T), par_shapes__cmp1);
1545
1546     // Apply the reorder mapping to the XYZ coordinate data.
1547     float* newpts = PAR_MALLOC(float, mesh->npoints * 3);
1548     PAR_SHAPES_T* invmap = PAR_MALLOC(PAR_SHAPES_T, mesh->npoints);
1549     float* dstpt = newpts;
1550     for (int i = 0; i < mesh->npoints; i++) {
1551         invmap[sortmap[i]] = i;
1552         float const* srcpt = mesh->points + 3 * sortmap[i];
1553         *dstpt++ = *srcpt++;
1554         *dstpt++ = *srcpt++;
1555         *dstpt++ = *srcpt++;
1556     }
1557     PAR_FREE(mesh->points);
1558     mesh->points = newpts;
1559
1560     // Apply the inverse reorder mapping to the triangle indices.
1561     PAR_SHAPES_T* newinds = PAR_MALLOC(PAR_SHAPES_T, mesh->ntriangles * 3);
1562     PAR_SHAPES_T* dstind = newinds;
1563     PAR_SHAPES_T const* srcind = mesh->triangles;
1564     for (int i = 0; i < mesh->ntriangles * 3; i++) {
1565         *dstind++ = invmap[*srcind++];
1566     }
1567     PAR_FREE(mesh->triangles);
1568     mesh->triangles = newinds;
1569
1570     // Cleanup.
1571     memcpy(sortmap, invmap, sizeof(PAR_SHAPES_T) * mesh->npoints);
1572     PAR_FREE(invmap);
1573 }
1574
1575 static void par_shapes__weld_points(par_shapes_mesh* mesh, int gridsize,
1576     float epsilon, PAR_SHAPES_T* weldmap)
1577 {
1578     // Each bin contains a "pointer" (really an index) to its first point.
1579     // We add 1 because 0 is reserved to mean that the bin is empty.
1580     // Since the points are spatially sorted, there's no need to store
1581     // a point count in each bin.
1582     PAR_SHAPES_T* bins = PAR_CALLOC(PAR_SHAPES_T,
1583         gridsize * gridsize * gridsize);
1584     int prev_binindex = -1;
1585     for (int p = 0; p < mesh->npoints; p++) {
1586         float const* pt = mesh->points + p * 3;
1587         int i = (int) pt[0];
1588         int j = (int) pt[1];
1589         int k = (int) pt[2];
1590         int this_binindex = i + gridsize * j + gridsize * gridsize * k;
1591         if (this_binindex != prev_binindex) {
1592             bins[this_binindex] = 1 + p;
1593         }
1594         prev_binindex = this_binindex;
1595     }
1596
1597     // Examine all bins that intersect the epsilon-sized cube centered at each
1598     // point, and check for colocated points within those bins.
1599     float const* pt = mesh->points;
1600     int nremoved = 0;
1601     for (int p = 0; p < mesh->npoints; p++, pt += 3) {
1602
1603         // Skip if this point has already been welded.
1604         if (weldmap[p] != p) {
1605             continue;
1606         }
1607
1608         // Build a list of bins that intersect the epsilon-sized cube.
1609         int nearby[8];
1610         int nbins = 0;
1611         int minp[3], maxp[3];
1612         for (int c = 0; c < 3; c++) {
1613             minp[c] = (int) (pt[c] - epsilon);
1614             maxp[c] = (int) (pt[c] + epsilon);
1615         }
1616         for (int i = minp[0]; i <= maxp[0]; i++) {
1617             for (int j = minp[1]; j <= maxp[1]; j++) {
1618                 for (int k = minp[2]; k <= maxp[2]; k++) {
1619                     int binindex = i + gridsize * j + gridsize * gridsize * k;
1620                     PAR_SHAPES_T binvalue = *(bins + binindex);
1621                     if (binvalue > 0) {
1622                         if (nbins == 8) {
1623                             //printf("Epsilon value is too large.\n");
1624                             break;
1625                         }
1626                         nearby[nbins++] = binindex;
1627                     }
1628                 }
1629             }
1630         }
1631
1632         // Check for colocated points in each nearby bin.
1633         for (int b = 0; b < nbins; b++) {
1634             int binindex = nearby[b];
1635             PAR_SHAPES_T binvalue = *(bins + binindex);
1636             PAR_SHAPES_T nindex = binvalue - 1;
1637             while (true) {
1638
1639                 // If this isn't "self" and it's colocated, then weld it!
1640                 if (nindex != p && weldmap[nindex] == nindex) {
1641                     float const* thatpt = mesh->points + nindex * 3;
1642                     float dist2 = par_shapes__sqrdist3(thatpt, pt);
1643                     if (dist2 < epsilon) {
1644                         weldmap[nindex] = p;
1645                         nremoved++;
1646                     }
1647                 }
1648
1649                 // Advance to the next point if possible.
1650                 if (++nindex >= mesh->npoints) {
1651                     break;
1652                 }
1653
1654                 // If the next point is outside the bin, then we're done.
1655                 float const* nextpt = mesh->points + nindex * 3;
1656                 int i = (int) nextpt[0];
1657                 int j = (int) nextpt[1];
1658                 int k = (int) nextpt[2];
1659                 int nextbinindex = i + gridsize * j + gridsize * gridsize * k;
1660                 if (nextbinindex != binindex) {
1661                     break;
1662                 }
1663             }
1664         }
1665     }
1666     PAR_FREE(bins);
1667
1668     // Apply the weldmap to the vertices.
1669     int npoints = mesh->npoints - nremoved;
1670     float* newpts = PAR_MALLOC(float, 3 * npoints);
1671     float* dst = newpts;
1672     PAR_SHAPES_T* condensed_map = PAR_MALLOC(PAR_SHAPES_T, mesh->npoints);
1673     PAR_SHAPES_T* cmap = condensed_map;
1674     float const* src = mesh->points;
1675     int ci = 0;
1676     for (int p = 0; p < mesh->npoints; p++, src += 3) {
1677         if (weldmap[p] == p) {
1678             *dst++ = src[0];
1679             *dst++ = src[1];
1680             *dst++ = src[2];
1681             *cmap++ = ci++;
1682         } else {
1683             *cmap++ = condensed_map[weldmap[p]];
1684         }
1685     }
1686     assert(ci == npoints);
1687     PAR_FREE(mesh->points);
1688     memcpy(weldmap, condensed_map, mesh->npoints * sizeof(PAR_SHAPES_T));
1689     PAR_FREE(condensed_map);
1690     mesh->points = newpts;
1691     mesh->npoints = npoints;
1692
1693     // Apply the weldmap to the triangle indices and skip the degenerates.
1694     PAR_SHAPES_T const* tsrc = mesh->triangles;
1695     PAR_SHAPES_T* tdst = mesh->triangles;
1696     int ntriangles = 0;
1697     for (int i = 0; i < mesh->ntriangles; i++, tsrc += 3) {
1698         PAR_SHAPES_T a = weldmap[tsrc[0]];
1699         PAR_SHAPES_T b = weldmap[tsrc[1]];
1700         PAR_SHAPES_T c = weldmap[tsrc[2]];
1701         if (a != b && a != c && b != c) {
1702             *tdst++ = a;
1703             *tdst++ = b;
1704             *tdst++ = c;
1705             ntriangles++;
1706         }
1707     }
1708     mesh->ntriangles = ntriangles;
1709 }
1710
1711 par_shapes_mesh* par_shapes_weld(par_shapes_mesh const* mesh, float epsilon,
1712     PAR_SHAPES_T* weldmap)
1713 {
1714     par_shapes_mesh* clone = par_shapes_clone(mesh, 0);
1715     float aabb[6];
1716     int gridsize = 20;
1717     float maxcell = gridsize - 1;
1718     par_shapes_compute_aabb(clone, aabb);
1719     float scale[3] = {
1720         aabb[3] == aabb[0] ? 1.0f : maxcell / (aabb[3] - aabb[0]),
1721         aabb[4] == aabb[1] ? 1.0f : maxcell / (aabb[4] - aabb[1]),
1722         aabb[5] == aabb[2] ? 1.0f : maxcell / (aabb[5] - aabb[2]),
1723     };
1724     par_shapes_translate(clone, -aabb[0], -aabb[1], -aabb[2]);
1725     par_shapes_scale(clone, scale[0], scale[1], scale[2]);
1726     PAR_SHAPES_T* sortmap = PAR_MALLOC(PAR_SHAPES_T, mesh->npoints);
1727     par_shapes__sort_points(clone, gridsize, sortmap);
1728     bool owner = false;
1729     if (!weldmap) {
1730         owner = true;
1731         weldmap = PAR_MALLOC(PAR_SHAPES_T, mesh->npoints);
1732     }
1733     for (int i = 0; i < mesh->npoints; i++) {
1734         weldmap[i] = i;
1735     }
1736     par_shapes__weld_points(clone, gridsize, epsilon, weldmap);
1737     if (owner) {
1738         PAR_FREE(weldmap);
1739     } else {
1740         PAR_SHAPES_T* newmap = PAR_MALLOC(PAR_SHAPES_T, mesh->npoints);
1741         for (int i = 0; i < mesh->npoints; i++) {
1742             newmap[i] = weldmap[sortmap[i]];
1743         }
1744         memcpy(weldmap, newmap, sizeof(PAR_SHAPES_T) * mesh->npoints);
1745         PAR_FREE(newmap);
1746     }
1747     PAR_FREE(sortmap);
1748     par_shapes_scale(clone, 1.0 / scale[0], 1.0 / scale[1], 1.0 / scale[2]);
1749     par_shapes_translate(clone, aabb[0], aabb[1], aabb[2]);
1750     return clone;
1751 }
1752
1753 // -----------------------------------------------------------------------------
1754 // BEGIN OPEN SIMPLEX NOISE
1755 // -----------------------------------------------------------------------------
1756
1757 #define STRETCH_CONSTANT_2D (-0.211324865405187)  // (1 / sqrt(2 + 1) - 1 ) / 2;
1758 #define SQUISH_CONSTANT_2D (0.366025403784439)  // (sqrt(2 + 1) -1) / 2;
1759 #define STRETCH_CONSTANT_3D (-1.0 / 6.0)  // (1 / sqrt(3 + 1) - 1) / 3;
1760 #define SQUISH_CONSTANT_3D (1.0 / 3.0)  // (sqrt(3+1)-1)/3;
1761 #define STRETCH_CONSTANT_4D (-0.138196601125011)  // (1 / sqrt(4 + 1) - 1) / 4;
1762 #define SQUISH_CONSTANT_4D (0.309016994374947)  // (sqrt(4 + 1) - 1) / 4;
1763
1764 #define NORM_CONSTANT_2D (47.0)
1765 #define NORM_CONSTANT_3D (103.0)
1766 #define NORM_CONSTANT_4D (30.0)
1767
1768 #define DEFAULT_SEED (0LL)
1769
1770 struct osn_context {
1771     int16_t* perm;
1772     int16_t* permGradIndex3D;
1773 };
1774
1775 #define ARRAYSIZE(x) (sizeof((x)) / sizeof((x)[0]))
1776
1777 /*
1778  * Gradients for 2D. They approximate the directions to the
1779  * vertices of an octagon from the center.
1780  */
1781 static const int8_t gradients2D[] = {
1782     5, 2, 2, 5, -5, 2, -2, 5, 5, -2, 2, -5, -5, -2, -2, -5,
1783 };
1784
1785 /*
1786  * Gradients for 3D. They approximate the directions to the
1787  * vertices of a rhombicuboctahedron from the center, skewed so
1788  * that the triangular and square facets can be inscribed inside
1789  * circles of the same radius.
1790  */
1791 static const signed char gradients3D[] = {
1792     -11, 4, 4, -4, 11, 4, -4, 4, 11, 11, 4, 4, 4, 11, 4, 4, 4, 11, -11, -4, 4,
1793     -4, -11, 4, -4, -4, 11, 11, -4, 4, 4, -11, 4, 4, -4, 11, -11, 4, -4, -4, 11,
1794     -4, -4, 4, -11, 11, 4, -4, 4, 11, -4, 4, 4, -11, -11, -4, -4, -4, -11, -4,
1795     -4, -4, -11, 11, -4, -4, 4, -11, -4, 4, -4, -11,
1796 };
1797
1798 /*
1799  * Gradients for 4D. They approximate the directions to the
1800  * vertices of a disprismatotesseractihexadecachoron from the center,
1801  * skewed so that the tetrahedral and cubic facets can be inscribed inside
1802  * spheres of the same radius.
1803  */
1804 static const signed char gradients4D[] = {
1805     3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3, -3, 1, 1, 1, -1, 3, 1, 1,
1806     -1, 1, 3, 1, -1, 1, 1, 3, 3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3, 1, 1, -1, 1,
1807     3, -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3, 3, 1, -1, 1, 1,
1808     3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3, -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3,
1809     1, -1, 1, -1, 3, 3, -1, -1, 1, 1, -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3, -3,
1810     -1, -1, 1, -1, -3, -1, 1, -1, -1, -3, 1, -1, -1, -1, 3, 3, 1, 1, -1, 1, 3,
1811     1, -1, 1, 1, 3, -1, 1, 1, 1, -3, -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3, -1,
1812     -1, 1, 1, -3, 3, -1, 1, -1, 1, -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3, -3,
1813     -1, 1, -1, -1, -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3, 3, 1, -1, -1, 1, 3,
1814     -1, -1, 1, 1, -3, -1, 1, 1, -1, -3, -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3,
1815     -1, -1, 1, -1, -3, 3, -1, -1, -1, 1, -3, -1, -1, 1, -1, -3, -1, 1, -1, -1,
1816     -3, -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3,
1817 };
1818
1819 static double extrapolate2(
1820     struct osn_context* ctx, int xsb, int ysb, double dx, double dy)
1821 {
1822     int16_t* perm = ctx->perm;
1823     int index = perm[(perm[xsb & 0xFF] + ysb) & 0xFF] & 0x0E;
1824     return gradients2D[index] * dx + gradients2D[index + 1] * dy;
1825 }
1826
1827 static inline int fastFloor(double x)
1828 {
1829     int xi = (int) x;
1830     return x < xi ? xi - 1 : xi;
1831 }
1832
1833 static int allocate_perm(struct osn_context* ctx, int nperm, int ngrad)
1834 {
1835     PAR_FREE(ctx->perm);
1836     PAR_FREE(ctx->permGradIndex3D);
1837     ctx->perm = PAR_MALLOC(int16_t, nperm);
1838     if (!ctx->perm) {
1839         return -ENOMEM;
1840     }
1841     ctx->permGradIndex3D = PAR_MALLOC(int16_t, ngrad);
1842     if (!ctx->permGradIndex3D) {
1843         PAR_FREE(ctx->perm);
1844         return -ENOMEM;
1845     }
1846     return 0;
1847 }
1848
1849 static int par__simplex_noise(int64_t seed, struct osn_context** ctx)
1850 {
1851     int rc;
1852     int16_t source[256];
1853     int i;
1854     int16_t* perm;
1855     int16_t* permGradIndex3D;
1856     *ctx = PAR_MALLOC(struct osn_context, 1);
1857     if (!(*ctx)) {
1858         return -ENOMEM;
1859     }
1860     (*ctx)->perm = NULL;
1861     (*ctx)->permGradIndex3D = NULL;
1862     rc = allocate_perm(*ctx, 256, 256);
1863     if (rc) {
1864         PAR_FREE(*ctx);
1865         return rc;
1866     }
1867     perm = (*ctx)->perm;
1868     permGradIndex3D = (*ctx)->permGradIndex3D;
1869     for (i = 0; i < 256; i++) {
1870         source[i] = (int16_t) i;
1871     }
1872     seed = seed * 6364136223846793005LL + 1442695040888963407LL;
1873     seed = seed * 6364136223846793005LL + 1442695040888963407LL;
1874     seed = seed * 6364136223846793005LL + 1442695040888963407LL;
1875     for (i = 255; i >= 0; i--) {
1876         seed = seed * 6364136223846793005LL + 1442695040888963407LL;
1877         int r = (int) ((seed + 31) % (i + 1));
1878         if (r < 0)
1879             r += (i + 1);
1880         perm[i] = source[r];
1881         permGradIndex3D[i] =
1882             (short) ((perm[i] % (ARRAYSIZE(gradients3D) / 3)) * 3);
1883         source[r] = source[i];
1884     }
1885     return 0;
1886 }
1887
1888 static void par__simplex_noise_free(struct osn_context* ctx)
1889 {
1890     if (!ctx)
1891         return;
1892     if (ctx->perm) {
1893         PAR_FREE(ctx->perm);
1894         ctx->perm = NULL;
1895     }
1896     if (ctx->permGradIndex3D) {
1897         PAR_FREE(ctx->permGradIndex3D);
1898         ctx->permGradIndex3D = NULL;
1899     }
1900     PAR_FREE(ctx);
1901 }
1902
1903 static double par__simplex_noise2(struct osn_context* ctx, double x, double y)
1904 {
1905     // Place input coordinates onto grid.
1906     double stretchOffset = (x + y) * STRETCH_CONSTANT_2D;
1907     double xs = x + stretchOffset;
1908     double ys = y + stretchOffset;
1909
1910     // Floor to get grid coordinates of rhombus (stretched square) super-cell
1911     // origin.
1912     int xsb = fastFloor(xs);
1913     int ysb = fastFloor(ys);
1914
1915     // Skew out to get actual coordinates of rhombus origin. We'll need these
1916     // later.
1917     double squishOffset = (xsb + ysb) * SQUISH_CONSTANT_2D;
1918     double xb = xsb + squishOffset;
1919     double yb = ysb + squishOffset;
1920
1921     // Compute grid coordinates relative to rhombus origin.
1922     double xins = xs - xsb;
1923     double yins = ys - ysb;
1924
1925     // Sum those together to get a value that determines which region we're in.
1926     double inSum = xins + yins;
1927
1928     // Positions relative to origin point.
1929     double dx0 = x - xb;
1930     double dy0 = y - yb;
1931
1932     // We'll be defining these inside the next block and using them afterwards.
1933     double dx_ext, dy_ext;
1934     int xsv_ext, ysv_ext;
1935
1936     double value = 0;
1937
1938     // Contribution (1,0)
1939     double dx1 = dx0 - 1 - SQUISH_CONSTANT_2D;
1940     double dy1 = dy0 - 0 - SQUISH_CONSTANT_2D;
1941     double attn1 = 2 - dx1 * dx1 - dy1 * dy1;
1942     if (attn1 > 0) {
1943         attn1 *= attn1;
1944         value += attn1 * attn1 * extrapolate2(ctx, xsb + 1, ysb + 0, dx1, dy1);
1945     }
1946
1947     // Contribution (0,1)
1948     double dx2 = dx0 - 0 - SQUISH_CONSTANT_2D;
1949     double dy2 = dy0 - 1 - SQUISH_CONSTANT_2D;
1950     double attn2 = 2 - dx2 * dx2 - dy2 * dy2;
1951     if (attn2 > 0) {
1952         attn2 *= attn2;
1953         value += attn2 * attn2 * extrapolate2(ctx, xsb + 0, ysb + 1, dx2, dy2);
1954     }
1955
1956     if (inSum <= 1) {  // We're inside the triangle (2-Simplex) at (0,0)
1957         double zins = 1 - inSum;
1958         if (zins > xins || zins > yins) {
1959             if (xins > yins) {
1960                 xsv_ext = xsb + 1;
1961                 ysv_ext = ysb - 1;
1962                 dx_ext = dx0 - 1;
1963                 dy_ext = dy0 + 1;
1964             } else {
1965                 xsv_ext = xsb - 1;
1966                 ysv_ext = ysb + 1;
1967                 dx_ext = dx0 + 1;
1968                 dy_ext = dy0 - 1;
1969             }
1970         } else {  //(1,0) and (0,1) are the closest two vertices.
1971             xsv_ext = xsb + 1;
1972             ysv_ext = ysb + 1;
1973             dx_ext = dx0 - 1 - 2 * SQUISH_CONSTANT_2D;
1974             dy_ext = dy0 - 1 - 2 * SQUISH_CONSTANT_2D;
1975         }
1976     } else {  // We're inside the triangle (2-Simplex) at (1,1)
1977         double zins = 2 - inSum;
1978         if (zins < xins || zins < yins) {
1979             if (xins > yins) {
1980                 xsv_ext = xsb + 2;
1981                 ysv_ext = ysb + 0;
1982                 dx_ext = dx0 - 2 - 2 * SQUISH_CONSTANT_2D;
1983                 dy_ext = dy0 + 0 - 2 * SQUISH_CONSTANT_2D;
1984             } else {
1985                 xsv_ext = xsb + 0;
1986                 ysv_ext = ysb + 2;
1987                 dx_ext = dx0 + 0 - 2 * SQUISH_CONSTANT_2D;
1988                 dy_ext = dy0 - 2 - 2 * SQUISH_CONSTANT_2D;
1989             }
1990         } else {  //(1,0) and (0,1) are the closest two vertices.
1991             dx_ext = dx0;
1992             dy_ext = dy0;
1993             xsv_ext = xsb;
1994             ysv_ext = ysb;
1995         }
1996         xsb += 1;
1997         ysb += 1;
1998         dx0 = dx0 - 1 - 2 * SQUISH_CONSTANT_2D;
1999         dy0 = dy0 - 1 - 2 * SQUISH_CONSTANT_2D;
2000     }
2001
2002     // Contribution (0,0) or (1,1)
2003     double attn0 = 2 - dx0 * dx0 - dy0 * dy0;
2004     if (attn0 > 0) {
2005         attn0 *= attn0;
2006         value += attn0 * attn0 * extrapolate2(ctx, xsb, ysb, dx0, dy0);
2007     }
2008
2009     // Extra Vertex
2010     double attn_ext = 2 - dx_ext * dx_ext - dy_ext * dy_ext;
2011     if (attn_ext > 0) {
2012         attn_ext *= attn_ext;
2013         value += attn_ext * attn_ext *
2014             extrapolate2(ctx, xsv_ext, ysv_ext, dx_ext, dy_ext);
2015     }
2016
2017     return value / NORM_CONSTANT_2D;
2018 }
2019
2020 void par_shapes_remove_degenerate(par_shapes_mesh* mesh, float mintriarea)
2021 {
2022     int ntriangles = 0;
2023     PAR_SHAPES_T* triangles = PAR_MALLOC(PAR_SHAPES_T, mesh->ntriangles * 3);
2024     PAR_SHAPES_T* dst = triangles;
2025     PAR_SHAPES_T const* src = mesh->triangles;
2026     float next[3], prev[3], cp[3];
2027     float mincplen2 = (mintriarea * 2) * (mintriarea * 2);
2028     for (int f = 0; f < mesh->ntriangles; f++, src += 3) {
2029         float const* pa = mesh->points + 3 * src[0];
2030         float const* pb = mesh->points + 3 * src[1];
2031         float const* pc = mesh->points + 3 * src[2];
2032         par_shapes__copy3(next, pb);
2033         par_shapes__subtract3(next, pa);
2034         par_shapes__copy3(prev, pc);
2035         par_shapes__subtract3(prev, pa);
2036         par_shapes__cross3(cp, next, prev);
2037         float cplen2 = par_shapes__dot3(cp, cp);
2038         if (cplen2 >= mincplen2) {
2039             *dst++ = src[0];
2040             *dst++ = src[1];
2041             *dst++ = src[2];
2042             ntriangles++;
2043         }
2044     }
2045     mesh->ntriangles = ntriangles;
2046     PAR_FREE(mesh->triangles);
2047     mesh->triangles = triangles;
2048 }
2049
2050 #endif // PAR_SHAPES_IMPLEMENTATION
2051 #endif // PAR_SHAPES_H