GCC Code Coverage Report

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