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