1 // SHAPES :: https://github.com/prideout/par
2 // Simple C library for creation and manipulation of triangle meshes.
4 // The API is divided into three sections:
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.
10 // In addition to the comment block above each function declaration, the API
11 // has informal documentation here:
13 // https://prideout.net/shapes
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.
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,
25 // Copyright (c) 2015 Philip Rideout
36 // Ray: commented to avoid conflict with raylib bool
38 #if !defined(_MSC_VER)
43 # else // stdbool.h missing prior to MSVC++ 12.0 (VS2013)
52 #define PAR_SHAPES_T uint16_t
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...)
64 void par_shapes_free_mesh(par_shapes_mesh*);
66 // Generators ------------------------------------------------------------------
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);
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);
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);
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);
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,
90 par_shapes_mesh* par_shapes_create_hemisphere(int slices, int stacks);
91 par_shapes_mesh* par_shapes_create_plane(int slices, int stacks);
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);
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();
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();
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();
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);
117 // Create an empty shape. Useful for building scenes with merge_and_free.
118 par_shapes_mesh* par_shapes_create_empty();
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);
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,
131 // Queries ---------------------------------------------------------------------
133 // Dump out a text file conforming to the venerable OBJ format.
134 void par_shapes_export(par_shapes_mesh const*, char const* objfile);
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);
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);
144 // Transformations -------------------------------------------------------------
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);
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);
156 // Remove all triangles whose area is less than minarea.
157 void par_shapes_remove_degenerate(par_shapes_mesh*, float minarea);
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);
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);
172 // Compute smooth normals by averaging adjacent facet normals.
173 void par_shapes_compute_normals(par_shapes_mesh* m);
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))
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)
195 // -----------------------------------------------------------------------------
197 // -----------------------------------------------------------------------------
199 #ifdef PAR_SHAPES_IMPLEMENTATION
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*);
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);
221 static void par_shapes__copy3(float* result, float const* a)
228 static float par_shapes__dot3(float const* a, float const* b)
230 return b[0] * a[0] + b[1] * a[1] + b[2] * a[2];
233 static void par_shapes__transform3(float* p, float const* x, float const* y,
236 float px = par_shapes__dot3(p, x);
237 float py = par_shapes__dot3(p, y);
238 float pz = par_shapes__dot3(p, z);
244 static void par_shapes__cross3(float* result, float const* a, float const* b)
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]);
254 static void par_shapes__mix3(float* d, float const* a, float const* b, float t)
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);
264 static void par_shapes__scale3(float* result, float a)
271 static void par_shapes__normalize3(float* v)
273 float lsqr = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
275 par_shapes__scale3(v, 1.0f / lsqr);
279 static void par_shapes__subtract3(float* result, float const* a)
286 static void par_shapes__add3(float* result, float const* a)
293 static float par_shapes__sqrdist3(float const* a, float const* b)
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;
301 static void par_shapes__compute_welded_normals(par_shapes_mesh* m)
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) {
310 float const* pnormal = welded->normals + d * 3;
311 pdst[0] = pnormal[0];
312 pdst[1] = pnormal[1];
313 pdst[2] = pnormal[2];
316 par_shapes_free_mesh(welded);
319 par_shapes_mesh* par_shapes_create_cylinder(int slices, int stacks)
321 if (slices < 3 || stacks < 1) {
324 return par_shapes_create_parametric(par_shapes__cylinder, slices,
328 par_shapes_mesh* par_shapes_create_parametric_sphere(int slices, int stacks)
330 if (slices < 3 || stacks < 3) {
333 par_shapes_mesh* m = par_shapes_create_parametric(par_shapes__sphere,
335 par_shapes_remove_degenerate(m, 0.0001);
339 par_shapes_mesh* par_shapes_create_hemisphere(int slices, int stacks)
341 if (slices < 3 || stacks < 3) {
344 par_shapes_mesh* m = par_shapes_create_parametric(par_shapes__hemisphere,
346 par_shapes_remove_degenerate(m, 0.0001);
350 par_shapes_mesh* par_shapes_create_torus(int slices, int stacks, float radius)
352 if (slices < 3 || stacks < 3) {
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,
362 par_shapes_mesh* par_shapes_create_klein_bottle(int slices, int stacks)
364 if (slices < 3 || stacks < 3) {
367 par_shapes_mesh* mesh = par_shapes_create_parametric(
368 par_shapes__klein, slices, stacks, 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);
377 par_shapes__compute_welded_normals(mesh);
381 par_shapes_mesh* par_shapes_create_trefoil_knot(int slices, int stacks,
384 if (slices < 3 || stacks < 3) {
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,
394 par_shapes_mesh* par_shapes_create_plane(int slices, int stacks)
396 if (slices < 1 || stacks < 1) {
399 return par_shapes_create_parametric(par_shapes__plane, slices,
403 par_shapes_mesh* par_shapes_create_parametric(par_shapes_fn fn,
404 int slices, int stacks, void* userdata)
406 par_shapes_mesh* mesh = PAR_CALLOC(par_shapes_mesh, 1);
409 mesh->npoints = (slices + 1) * (stacks + 1);
410 mesh->points = PAR_CALLOC(float, 3 * mesh->npoints);
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);
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;
438 mesh->ntriangles = 2 * slices * stacks;
439 mesh->triangles = PAR_CALLOC(PAR_SHAPES_T, 3 * mesh->ntriangles);
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;
448 *face++ = v + slice + slices + 1;
449 *face++ = v + next + slices + 1;
455 par_shapes__compute_welded_normals(mesh);
459 void par_shapes_free_mesh(par_shapes_mesh* mesh)
461 PAR_FREE(mesh->points);
462 PAR_FREE(mesh->triangles);
463 PAR_FREE(mesh->normals);
464 PAR_FREE(mesh->tcoords);
468 void par_shapes_export(par_shapes_mesh const* mesh, char const* filename)
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]);
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);
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]);
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);
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]);
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);
518 for (int nvert = 0; nvert < mesh->npoints; nvert++) {
519 fprintf(objfile, "v %f %f %f\n", points[0], points[1], points[2]);
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);
532 static void par_shapes__sphere(float const* uv, float* xyz, void* userdata)
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);
541 static void par_shapes__hemisphere(float const* uv, float* xyz, void* userdata)
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);
550 static void par_shapes__plane(float const* uv, float* xyz, void* userdata)
557 static void par_shapes__klein(float const* uv, float* xyz, void* userdata)
559 float u = uv[0] * PAR_PI;
560 float v = uv[1] * 2 * PAR_PI;
563 xyz[0] = 3 * cosf(u) * (1 + sinf(u)) + (2 * (1 - cosf(u) / 2)) *
565 xyz[2] = -8 * sinf(u) - 2 * (1 - cosf(u) / 2) * sinf(u) * cosf(v);
567 xyz[0] = 3 * cosf(u) * (1 + sinf(u)) + (2 * (1 - cosf(u) / 2)) *
569 xyz[2] = -8 * sinf(u);
571 xyz[1] = -2 * (1 - cosf(u) / 2) * sinf(v);
574 static void par_shapes__cylinder(float const* uv, float* xyz, void* userdata)
576 float theta = uv[1] * 2 * PAR_PI;
577 xyz[0] = sinf(theta);
578 xyz[1] = cosf(theta);
582 static void par_shapes__torus(float const* uv, float* xyz, void* userdata)
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;
594 static void par_shapes__trefoil(float const* uv, float* xyz, void* userdata)
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);
609 -1.5f * b * sin(1.5f * u) * cos(u) - (a + b * cos(1.5f * u)) * sin(u);
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);
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);
623 void par_shapes_merge(par_shapes_mesh* dst, par_shapes_mesh const* src)
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);
634 memcpy(dst->normals + 3 * offset, src->normals,
635 vecsize * src->npoints);
638 if (src->tcoords || dst->tcoords) {
639 int uvsize = sizeof(float) * 2;
640 dst->tcoords = PAR_REALLOC(float, dst->tcoords, 2 * npoints);
642 memcpy(dst->tcoords + 2 * offset, src->tcoords,
643 uvsize * src->npoints);
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++;
655 dst->ntriangles = ntriangles;
658 par_shapes_mesh* par_shapes_create_disk(float radius, int slices,
659 float const* center, float const* normal)
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;
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);
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];
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++) {
688 *triangles++ = 1 + i;
689 *triangles++ = 1 + (i + 1) % slices;
691 float k[3] = {0, 0, -1};
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]);
700 par_shapes_mesh* par_shapes_create_empty()
702 return PAR_CALLOC(par_shapes_mesh, 1);
705 void par_shapes_translate(par_shapes_mesh* m, float x, float y, float z)
707 float* points = m->points;
708 for (int i = 0; i < m->npoints; i++) {
715 void par_shapes_rotate(par_shapes_mesh* mesh, float radians, float const* axis)
717 float s = sinf(radians);
718 float c = cosf(radians);
725 float oneMinusC = 1.0f - c;
727 (((x * x) * oneMinusC) + c),
728 ((xy * oneMinusC) + (z * s)), ((zx * oneMinusC) - (y * s))
731 ((xy * oneMinusC) - (z * s)),
732 (((y * y) * oneMinusC) + c), ((yz * oneMinusC) + (x * s))
735 ((zx * oneMinusC) + (y * s)),
736 ((yz * oneMinusC) - (x * s)), (((z * z) * oneMinusC) + c)
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];
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];
760 void par_shapes_scale(par_shapes_mesh* m, float x, float y, float z)
762 float* points = m->points;
763 for (int i = 0; i < m->npoints; i++) {
770 void par_shapes_merge_and_free(par_shapes_mesh* dst, par_shapes_mesh* src)
772 par_shapes_merge(dst, src);
773 par_shapes_free_mesh(src);
776 void par_shapes_compute_aabb(par_shapes_mesh const* m, float* aabb)
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];
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]);
793 void par_shapes_invert(par_shapes_mesh* m, int face, int nfaces)
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]);
803 par_shapes_mesh* par_shapes_create_icosahedron()
805 static float verts[] = {
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,
819 static PAR_SHAPES_T faces[] = {
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));
851 par_shapes_mesh* par_shapes_create_dodecahedron()
853 static float verts[20 * 3] = {
856 -0.491, 0.357, 0.795,
857 -0.491, -0.357, 0.795,
858 0.188, -0.577, 0.795,
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,
875 static PAR_SHAPES_T pentagons[12 * 5] = {
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];
913 par_shapes_mesh* par_shapes_create_octahedron()
915 static float verts[6 * 3] = {
919 -1.000, 0.000, 0.000,
920 0.000, -1.000, 0.000,
923 static PAR_SHAPES_T triangles[8 * 3] = {
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++;
951 par_shapes_mesh* par_shapes_create_tetrahedron()
953 static float verts[4 * 3] = {
959 static PAR_SHAPES_T triangles[4 * 3] = {
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++;
983 par_shapes_mesh* par_shapes_create_cube()
985 static float verts[8 * 3] = {
995 static PAR_SHAPES_T quads[6 * 4] = {
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) {
1027 } par_shapes__command;
1033 par_shapes__command* commands;
1040 par_shapes_mesh* orientation;
1041 par_shapes__rule* rule;
1042 } par_shapes__stackframe;
1044 static par_shapes__rule* par_shapes__pick_rule(const char* name,
1045 par_shapes__rule* rules, int nrules)
1047 par_shapes__rule* rule = 0;
1049 for (int i = 0; i < nrules; i++) {
1051 if (!strcmp(rule->name, name)) {
1052 total += rule->weight;
1055 float r = (float) rand() / RAND_MAX;
1057 for (int i = 0; i < nrules; i++) {
1059 if (!strcmp(rule->name, name)) {
1060 t += (float) rule->weight / total;
1069 static par_shapes_mesh* par_shapes__create_turtle()
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);
1083 static par_shapes_mesh* par_shapes__apply_turtle(par_shapes_mesh* mesh,
1084 par_shapes_mesh* turtle, float const* pos, float const* scale)
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;
1092 par_shapes__transform3(pt,
1093 turtle->points + 0, turtle->points + 3, turtle->points + 6);
1101 static void par_shapes__connect(par_shapes_mesh* scene,
1102 par_shapes_mesh* cylinder, int slices)
1105 int npoints = (slices + 1) * (stacks + 1);
1106 assert(scene->npoints >= npoints && "Cannot connect to empty scene.");
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;
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;
1129 *face++ = v + slice;
1130 *face++ = v + slice + slices + 1;
1131 *face++ = v + next + slices + 1;
1136 PAR_FREE(scene->triangles);
1137 scene->triangles = triangles;
1139 scene->npoints = npoints;
1140 scene->ntriangles = ntriangles;
1143 par_shapes_mesh* par_shapes_create_lsystem(char const* text, int slices,
1147 program = PAR_MALLOC(char, strlen(text) + 1);
1149 // The first pass counts the number of rules and commands.
1150 strcpy(program, text);
1151 char *cmd = strtok(program, " ");
1155 char *arg = strtok(0, " ");
1157 //puts("lsystem error: unexpected end of program.");
1160 if (!strcmp(cmd, "rule")) {
1165 cmd = strtok(0, " ");
1169 par_shapes__rule* rules = PAR_MALLOC(par_shapes__rule, nrules);
1170 par_shapes__command* commands = PAR_MALLOC(par_shapes__command, ncommands);
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;
1180 // The second pass fills in the structures.
1181 strcpy(program, text);
1182 cmd = strtok(program, " ");
1184 char *arg = strtok(0, " ");
1185 if (!strcmp(cmd, "rule")) {
1188 // Split the argument into a rule name and weight.
1189 char* dot = strchr(arg, '.');
1191 current_rule->weight = atoi(dot + 1);
1194 current_rule->weight = 1;
1197 current_rule->name = arg;
1198 current_rule->ncommands = 0;
1199 current_rule->commands = current_command;
1201 current_rule->ncommands++;
1202 current_command->cmd = cmd;
1203 current_command->arg = arg;
1206 cmd = strtok(0, " ");
1209 // For testing purposes, dump out the parsed program.
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);
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();
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);
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};
1240 // Execute the L-system program until the stack size is 0.
1241 par_shapes__stackframe* stack =
1242 PAR_CALLOC(par_shapes__stackframe, maxdepth);
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);
1259 par_shapes__command* cmd = rule->commands + (frame->pc++);
1261 //printf("%5s %5s %5s:%d %03d\n", cmd->cmd, cmd->arg, rule->name, frame->pc - 1, stackptr);
1265 if (!strcmp(cmd->cmd, "shape")) {
1266 par_shapes_mesh* m = par_shapes__apply_turtle(tube, turtle,
1268 if (!strcmp(cmd->arg, "connect")) {
1269 par_shapes__connect(scene, m, slices);
1271 par_shapes_merge(scene, m);
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];
1278 frame->orientation = par_shapes_clone(turtle, 0);
1280 par_shapes__copy3(frame->scale, scale);
1281 par_shapes__copy3(frame->position, position);
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};
1294 par_shapes__dot3(turtle->points + 0, vec),
1295 par_shapes__dot3(turtle->points + 3, vec),
1296 par_shapes__dot3(turtle->points + 6, vec)
1298 par_shapes__add3(position, t);
1299 } else if (!strcmp(cmd->cmd, "ty")) {
1300 float vec[3] = {0, value, 0};
1302 par_shapes__dot3(turtle->points + 0, vec),
1303 par_shapes__dot3(turtle->points + 3, vec),
1304 par_shapes__dot3(turtle->points + 6, vec)
1306 par_shapes__add3(position, t);
1307 } else if (!strcmp(cmd->cmd, "tz")) {
1308 float vec[3] = {0, 0, value};
1310 par_shapes__dot3(turtle->points + 0, vec),
1311 par_shapes__dot3(turtle->points + 3, vec),
1312 par_shapes__dot3(turtle->points + 6, vec)
1314 par_shapes__add3(position, t);
1315 } else if (!strcmp(cmd->cmd, "sx")) {
1317 } else if (!strcmp(cmd->cmd, "sy")) {
1319 } else if (!strcmp(cmd->cmd, "sz")) {
1321 } else if (!strcmp(cmd->cmd, "sa")) {
1335 void par_shapes_unweld(par_shapes_mesh* mesh, bool create_indices)
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++);
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++) {
1356 PAR_FREE(mesh->triangles);
1357 mesh->triangles = tris;
1361 void par_shapes_compute_normals(par_shapes_mesh* m)
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);
1390 float* normal = m->normals;
1391 for (int p = 0; p < m->npoints; p++, normal += 3) {
1392 par_shapes__normalize3(normal);
1396 static void par_shapes__subdivide(par_shapes_mesh* mesh)
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);
1424 PAR_FREE(mesh->points);
1425 mesh->points = points;
1426 mesh->npoints = npoints;
1427 mesh->ntriangles = ntriangles;
1430 par_shapes_mesh* par_shapes_create_subdivided_sphere(int nsubd)
1432 par_shapes_mesh* mesh = par_shapes_create_icosahedron();
1433 par_shapes_unweld(mesh, false);
1434 PAR_FREE(mesh->triangles);
1435 mesh->triangles = 0;
1437 par_shapes__subdivide(mesh);
1439 for (int i = 0; i < mesh->npoints; i++) {
1440 par_shapes__normalize3(mesh->points + i * 3);
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;
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);
1453 par_shapes_mesh* par_shapes_create_rock(int seed, int subd)
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]);
1463 n += a * par__simplex_noise2(ctx, f * pt[0], f * pt[2]);
1468 pt[1] = -pow(-pt[1], 0.5) / 2;
1471 par__simplex_noise_free(ctx);
1472 par_shapes_compute_normals(mesh);
1476 par_shapes_mesh* par_shapes_clone(par_shapes_mesh const* mesh,
1477 par_shapes_mesh* clone)
1480 clone = PAR_CALLOC(par_shapes_mesh, 1);
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 *
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);
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);
1504 float const* points;
1506 } par_shapes__sort_context;
1508 static int par_shapes__cmp1(const void *arg0, const void *arg1)
1510 const int g = par_shapes__sort_context.gridsize;
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;
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;
1528 // Return the ordering.
1529 if (index0 < index1) return -1;
1530 if (index0 > index1) return 1;
1534 static void par_shapes__sort_points(par_shapes_mesh* mesh, int gridsize,
1535 PAR_SHAPES_T* sortmap)
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++) {
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);
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++;
1557 PAR_FREE(mesh->points);
1558 mesh->points = newpts;
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++];
1567 PAR_FREE(mesh->triangles);
1568 mesh->triangles = newinds;
1571 memcpy(sortmap, invmap, sizeof(PAR_SHAPES_T) * mesh->npoints);
1575 static void par_shapes__weld_points(par_shapes_mesh* mesh, int gridsize,
1576 float epsilon, PAR_SHAPES_T* weldmap)
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;
1594 prev_binindex = this_binindex;
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;
1601 for (int p = 0; p < mesh->npoints; p++, pt += 3) {
1603 // Skip if this point has already been welded.
1604 if (weldmap[p] != p) {
1608 // Build a list of bins that intersect the epsilon-sized cube.
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);
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);
1623 //printf("Epsilon value is too large.\n");
1626 nearby[nbins++] = binindex;
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;
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;
1649 // Advance to the next point if possible.
1650 if (++nindex >= mesh->npoints) {
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) {
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;
1676 for (int p = 0; p < mesh->npoints; p++, src += 3) {
1677 if (weldmap[p] == p) {
1683 *cmap++ = condensed_map[weldmap[p]];
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;
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;
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) {
1708 mesh->ntriangles = ntriangles;
1711 par_shapes_mesh* par_shapes_weld(par_shapes_mesh const* mesh, float epsilon,
1712 PAR_SHAPES_T* weldmap)
1714 par_shapes_mesh* clone = par_shapes_clone(mesh, 0);
1717 float maxcell = gridsize - 1;
1718 par_shapes_compute_aabb(clone, aabb);
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]),
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);
1731 weldmap = PAR_MALLOC(PAR_SHAPES_T, mesh->npoints);
1733 for (int i = 0; i < mesh->npoints; i++) {
1736 par_shapes__weld_points(clone, gridsize, epsilon, weldmap);
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]];
1744 memcpy(weldmap, newmap, sizeof(PAR_SHAPES_T) * mesh->npoints);
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]);
1753 // -----------------------------------------------------------------------------
1754 // BEGIN OPEN SIMPLEX NOISE
1755 // -----------------------------------------------------------------------------
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;
1764 #define NORM_CONSTANT_2D (47.0)
1765 #define NORM_CONSTANT_3D (103.0)
1766 #define NORM_CONSTANT_4D (30.0)
1768 #define DEFAULT_SEED (0LL)
1770 struct osn_context {
1772 int16_t* permGradIndex3D;
1775 #define ARRAYSIZE(x) (sizeof((x)) / sizeof((x)[0]))
1778 * Gradients for 2D. They approximate the directions to the
1779 * vertices of an octagon from the center.
1781 static const int8_t gradients2D[] = {
1782 5, 2, 2, 5, -5, 2, -2, 5, 5, -2, 2, -5, -5, -2, -2, -5,
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.
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,
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.
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,
1819 static double extrapolate2(
1820 struct osn_context* ctx, int xsb, int ysb, double dx, double dy)
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;
1827 static inline int fastFloor(double x)
1830 return x < xi ? xi - 1 : xi;
1833 static int allocate_perm(struct osn_context* ctx, int nperm, int ngrad)
1835 PAR_FREE(ctx->perm);
1836 PAR_FREE(ctx->permGradIndex3D);
1837 ctx->perm = PAR_MALLOC(int16_t, nperm);
1841 ctx->permGradIndex3D = PAR_MALLOC(int16_t, ngrad);
1842 if (!ctx->permGradIndex3D) {
1843 PAR_FREE(ctx->perm);
1849 static int par__simplex_noise(int64_t seed, struct osn_context** ctx)
1852 int16_t source[256];
1855 int16_t* permGradIndex3D;
1856 *ctx = PAR_MALLOC(struct osn_context, 1);
1860 (*ctx)->perm = NULL;
1861 (*ctx)->permGradIndex3D = NULL;
1862 rc = allocate_perm(*ctx, 256, 256);
1867 perm = (*ctx)->perm;
1868 permGradIndex3D = (*ctx)->permGradIndex3D;
1869 for (i = 0; i < 256; i++) {
1870 source[i] = (int16_t) i;
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));
1880 perm[i] = source[r];
1881 permGradIndex3D[i] =
1882 (short) ((perm[i] % (ARRAYSIZE(gradients3D) / 3)) * 3);
1883 source[r] = source[i];
1888 static void par__simplex_noise_free(struct osn_context* ctx)
1893 PAR_FREE(ctx->perm);
1896 if (ctx->permGradIndex3D) {
1897 PAR_FREE(ctx->permGradIndex3D);
1898 ctx->permGradIndex3D = NULL;
1903 static double par__simplex_noise2(struct osn_context* ctx, double x, double y)
1905 // Place input coordinates onto grid.
1906 double stretchOffset = (x + y) * STRETCH_CONSTANT_2D;
1907 double xs = x + stretchOffset;
1908 double ys = y + stretchOffset;
1910 // Floor to get grid coordinates of rhombus (stretched square) super-cell
1912 int xsb = fastFloor(xs);
1913 int ysb = fastFloor(ys);
1915 // Skew out to get actual coordinates of rhombus origin. We'll need these
1917 double squishOffset = (xsb + ysb) * SQUISH_CONSTANT_2D;
1918 double xb = xsb + squishOffset;
1919 double yb = ysb + squishOffset;
1921 // Compute grid coordinates relative to rhombus origin.
1922 double xins = xs - xsb;
1923 double yins = ys - ysb;
1925 // Sum those together to get a value that determines which region we're in.
1926 double inSum = xins + yins;
1928 // Positions relative to origin point.
1929 double dx0 = x - xb;
1930 double dy0 = y - yb;
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;
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;
1944 value += attn1 * attn1 * extrapolate2(ctx, xsb + 1, ysb + 0, dx1, dy1);
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;
1953 value += attn2 * attn2 * extrapolate2(ctx, xsb + 0, ysb + 1, dx2, dy2);
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) {
1970 } else { //(1,0) and (0,1) are the closest two vertices.
1973 dx_ext = dx0 - 1 - 2 * SQUISH_CONSTANT_2D;
1974 dy_ext = dy0 - 1 - 2 * SQUISH_CONSTANT_2D;
1976 } else { // We're inside the triangle (2-Simplex) at (1,1)
1977 double zins = 2 - inSum;
1978 if (zins < xins || zins < yins) {
1982 dx_ext = dx0 - 2 - 2 * SQUISH_CONSTANT_2D;
1983 dy_ext = dy0 + 0 - 2 * SQUISH_CONSTANT_2D;
1987 dx_ext = dx0 + 0 - 2 * SQUISH_CONSTANT_2D;
1988 dy_ext = dy0 - 2 - 2 * SQUISH_CONSTANT_2D;
1990 } else { //(1,0) and (0,1) are the closest two vertices.
1998 dx0 = dx0 - 1 - 2 * SQUISH_CONSTANT_2D;
1999 dy0 = dy0 - 1 - 2 * SQUISH_CONSTANT_2D;
2002 // Contribution (0,0) or (1,1)
2003 double attn0 = 2 - dx0 * dx0 - dy0 * dy0;
2006 value += attn0 * attn0 * extrapolate2(ctx, xsb, ysb, dx0, dy0);
2010 double attn_ext = 2 - dx_ext * dx_ext - dy_ext * dy_ext;
2012 attn_ext *= attn_ext;
2013 value += attn_ext * attn_ext *
2014 extrapolate2(ctx, xsv_ext, ysv_ext, dx_ext, dy_ext);
2017 return value / NORM_CONSTANT_2D;
2020 void par_shapes_remove_degenerate(par_shapes_mesh* mesh, float mintriarea)
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) {
2045 mesh->ntriangles = ntriangles;
2046 PAR_FREE(mesh->triangles);
2047 mesh->triangles = triangles;
2050 #endif // PAR_SHAPES_IMPLEMENTATION
2051 #endif // PAR_SHAPES_H