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