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