ssol_shape.c (47231B)
1 /* Copyright (C) 2018-2026 |Meso|Star> (contact@meso-star.com) 2 * Copyright (C) 2016, 2018 CNRS 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #define _POSIX_C_SOURCE 200112L /* copysign support */ 18 19 #include "ssol.h" 20 #include "ssol_c.h" 21 #include "ssol_device_c.h" 22 #include "ssol_shape_c.h" 23 24 #include <rsys/double2.h> 25 #include <rsys/double3.h> 26 #include <rsys/double33.h> 27 #include <rsys/dynamic_array_double.h> 28 #include <rsys/dynamic_array_size_t.h> 29 #include <rsys/float3.h> 30 #include <rsys/mem_allocator.h> 31 #include <rsys/ref_count.h> 32 #include <rsys/rsys.h> 33 #include <rsys/math.h> 34 35 #include <star/scpr.h> 36 #include <star/s3dut.h> 37 38 #include <limits.h> /* UINT_MAX constant */ 39 #include <math.h> /* copysign function */ 40 41 struct mesh_context { 42 const double* coords; 43 const size_t* ids; 44 }; 45 46 struct quadric_mesh_context { 47 const double* coords; 48 const size_t* ids; 49 const union private_data* data; 50 enum ssol_quadric_type quadric_type; 51 const double* transform; /* 3x4 column major matrix */ 52 double lower[2]; 53 double upper[2]; 54 }; 55 56 struct get_ctx { 57 size_t nbvert; 58 double two_pi_over_nbvert; 59 double radius; 60 }; 61 62 /******************************************************************************* 63 * Helper functions 64 ******************************************************************************/ 65 static FINLINE int 66 check_plane(const struct ssol_quadric_plane* plane) 67 { 68 return plane != NULL; 69 } 70 71 static FINLINE int 72 check_parabol(const struct ssol_quadric_parabol* parabol) 73 { 74 return parabol && parabol->focal > 0; 75 } 76 77 static FINLINE int 78 check_hyperbol(const struct ssol_quadric_hyperbol* hyperbol) 79 { 80 return hyperbol && hyperbol->img_focal > 0 && hyperbol->real_focal > 0; 81 } 82 83 static FINLINE int 84 check_parabolic_cylinder 85 (const struct ssol_quadric_parabolic_cylinder* parabolic_cylinder) 86 { 87 return parabolic_cylinder && parabolic_cylinder->focal > 0; 88 } 89 90 static FINLINE int 91 check_hemisphere(const struct ssol_quadric_hemisphere* hemisphere) 92 { 93 return hemisphere && hemisphere->radius > 0; 94 } 95 96 static INLINE int 97 check_quadric(const struct ssol_quadric* quadric) 98 { 99 if(!quadric) return RES_BAD_ARG; 100 101 switch (quadric->type) { 102 case SSOL_QUADRIC_PLANE: 103 return check_plane(&quadric->data.plane); 104 case SSOL_QUADRIC_PARABOL: 105 return check_parabol(&quadric->data.parabol); 106 case SSOL_QUADRIC_HYPERBOL: 107 return check_hyperbol(&quadric->data.hyperbol); 108 case SSOL_QUADRIC_PARABOLIC_CYLINDER: 109 return check_parabolic_cylinder(&quadric->data.parabolic_cylinder); 110 case SSOL_QUADRIC_HEMISPHERE: 111 return check_hemisphere(&quadric->data.hemisphere); 112 default: return 0; 113 } 114 } 115 116 static INLINE int 117 check_carving(const struct ssol_carving* polygon) 118 { 119 /* We don't check that the polygon defines a not empty area in such case, the 120 * quadric is valid but can have zero surface */ 121 return polygon && polygon->get && polygon->nb_vertices > 0; 122 } 123 124 static INLINE int 125 check_punched_surface(const struct ssol_punched_surface* punched_surface) 126 { 127 size_t i; 128 129 if(!punched_surface 130 || (punched_surface->nb_carvings == 0 131 && punched_surface->quadric->type != SSOL_QUADRIC_HEMISPHERE) 132 || (punched_surface->nb_carvings && !punched_surface->carvings) 133 || !punched_surface->quadric 134 || !check_quadric(punched_surface->quadric)) 135 return 0; 136 137 FOR_EACH(i, 0, punched_surface->nb_carvings) { 138 if(!check_carving(&punched_surface->carvings[i])) 139 return 0; 140 } 141 /* We don't check that carvings define a non empty area in such case, the 142 * quadric is valid but has zero surface */ 143 return 1; 144 } 145 146 static INLINE int 147 check_shape(const struct ssol_shape* shape) 148 { 149 return shape && shape->dev && (unsigned)shape->type < SHAPE_TYPES_COUNT__; 150 } 151 152 static INLINE enum scpr_operation 153 ssol_to_scpr_clip_op(const enum ssol_clipping_op clip_op) 154 { 155 enum scpr_operation op; 156 switch(clip_op) { 157 case SSOL_AND: op = SCPR_AND; break; 158 case SSOL_SUB: op = SCPR_SUB; break; 159 default: FATAL("Unreachable code.\n"); break; 160 } 161 return op; 162 } 163 164 static void 165 mesh_get_ids(const size_t itri, size_t ids[3], void* ctx) 166 { 167 const size_t i = itri*3/*#ids per triangle*/; 168 const struct mesh_context* msh = ctx; 169 ASSERT(ids && ctx); 170 ids[0] = msh->ids[i+0]; 171 ids[1] = msh->ids[i+1]; 172 ids[2] = msh->ids[i+2]; 173 } 174 175 static void 176 mesh_get_pos(const size_t ivert, double pos[2], void* ctx) 177 { 178 const size_t i = ivert*2/*#coords per vertex*/; 179 const struct mesh_context* msh = ctx; 180 ASSERT(pos && ctx); 181 pos[0] = msh->coords[i+0]; 182 pos[1] = msh->coords[i+1]; 183 } 184 185 static void 186 quadric_mesh_get_uv(const unsigned ivert, float uv[2], void* ctx) 187 { 188 const size_t i = ivert*2/*#coords per vertex*/; 189 const struct quadric_mesh_context* msh = ctx; 190 double tmp[2]; 191 ASSERT(uv && ctx); 192 tmp[0] = (msh->coords[i+0] - msh->lower[0]) / (msh->upper[0] - msh->lower[0]); 193 tmp[1] = (msh->coords[i+1] - msh->lower[1]) / (msh->upper[1] - msh->lower[1]); 194 195 uv[0] = (float)tmp[0]; 196 uv[1] = (float)tmp[1]; 197 } 198 199 static void 200 quadric_mesh_get_ids(const unsigned itri, unsigned ids[3], void* ctx) 201 { 202 const size_t i = itri*3/*#ids per triangle*/; 203 const struct quadric_mesh_context* msh = ctx; 204 ASSERT(ids && ctx); 205 ids[0] = (unsigned)msh->ids[i+0]; 206 ids[1] = (unsigned)msh->ids[i+1]; 207 ids[2] = (unsigned)msh->ids[i+2]; 208 } 209 210 static void 211 quadric_mesh_plane_get_pos(const unsigned ivert, float pos[3], void* ctx) 212 { 213 const size_t i = ivert*2/*#coords per vertex*/; 214 const struct quadric_mesh_context* msh = ctx; 215 double p[3]; /* Temporary quadric space position */ 216 ASSERT(pos && ctx); 217 p[0] = (float)msh->coords[i+0]; 218 p[1] = (float)msh->coords[i+1]; 219 p[2] = 0.f; 220 221 /* Transform the position in object space */ 222 d33_muld3(p, msh->transform, p); 223 d3_add(p, p, msh->transform+9); 224 225 f3_set_d3(pos, p); 226 } 227 228 static FINLINE double 229 hyperbol_z 230 (const double p[2], 231 const struct priv_hyperbol_data* hyperbol) 232 { 233 const double z0 = hyperbol->g_square + hyperbol->abs_b; 234 const double r2 = p[0] * p[0] + p[1] * p[1]; 235 return hyperbol->abs_b * sqrt(1 + r2 * hyperbol->one_over_a_square) 236 + hyperbol->g_square - z0; 237 } 238 239 static FINLINE double 240 parabol_z 241 (const double p[2], 242 const struct priv_parabol_data* parabol) 243 { 244 const double r2 = p[0] * p[0] + p[1] * p[1]; 245 return r2 * parabol->one_over_4focal; 246 } 247 248 static FINLINE double 249 parabolic_cylinder_z 250 (const double p[2], 251 const struct priv_pcylinder_data* pcyl) 252 { 253 return (p[1] * p[1]) * pcyl->one_over_4focal; 254 } 255 256 static FINLINE double 257 hemisphere_z 258 (const double p[2], 259 const struct priv_hemisphere_data* hemisphere) 260 { 261 const double r2 = p[0] * p[0] + p[1] * p[1]; 262 const double z2 = hemisphere->sqr_radius - r2; 263 /* manage numerical unaccuracy */ 264 ASSERT(z2 >= -hemisphere->sqr_radius * FLT_EPSILON); 265 return (z2 > 0) ? -sqrt(z2) + hemisphere->radius : hemisphere->radius; 266 } 267 268 static void 269 quadric_mesh_parabol_get_pos(const unsigned ivert, float pos[3], void* ctx) 270 { 271 const size_t i = ivert*2/*#coords per vertex*/; 272 const struct quadric_mesh_context* msh = ctx; 273 double p[3]; /* Temporary quadric space position */ 274 ASSERT(pos && ctx); 275 p[0] = msh->coords[i+0]; 276 p[1] = msh->coords[i+1]; 277 p[2] = parabol_z(p, &msh->data->parabol); 278 279 /* Transform the position in object space */ 280 d33_muld3(p, msh->transform, p); 281 d3_add(p, p, msh->transform+9); 282 283 f3_set_d3(pos, p); 284 } 285 286 static void 287 quadric_mesh_hyperbol_get_pos(const unsigned ivert, float pos[3], void* ctx) 288 { 289 const size_t i = ivert * 2/*#coords per vertex*/; 290 const struct quadric_mesh_context* msh = ctx; 291 double p[3]; /* Temporary quadric space position */ 292 ASSERT(pos && ctx); 293 p[0] = msh->coords[i+0]; 294 p[1] = msh->coords[i+1]; 295 p[2] = hyperbol_z(p, &msh->data->hyperbol); 296 297 /* Transform the position in object space */ 298 d33_muld3(p, msh->transform, p); 299 d3_add(p, p, msh->transform+9); 300 301 f3_set_d3(pos, p); 302 } 303 304 static void 305 quadric_mesh_parabolic_cylinder_get_pos 306 (const unsigned ivert, float pos[3], void* ctx) 307 { 308 const size_t i = ivert*2/*#coords per vertex*/; 309 const struct quadric_mesh_context* msh = ctx; 310 double p[3]; /* Temporary quadric space position */ 311 ASSERT(pos && ctx); 312 p[0] = msh->coords[i+0]; 313 p[1] = msh->coords[i+1]; 314 p[2] = parabolic_cylinder_z(p, &msh->data->pcylinder); 315 316 /* Transform the position in object space */ 317 d33_muld3(p, msh->transform, p); 318 d3_add(p, p, msh->transform+9); 319 320 f3_set_d3(pos, p); 321 } 322 323 static void 324 quadric_mesh_hemisphere_get_pos 325 (const unsigned ivert, float pos[3], void* ctx) 326 { 327 const size_t i = ivert * 2/*#coords per vertex*/; 328 const struct quadric_mesh_context* msh = ctx; 329 double p[3]; /* Temporary quadric space position */ 330 ASSERT(pos && ctx); 331 p[0] = msh->coords[i + 0]; 332 p[1] = msh->coords[i + 1]; 333 p[2] = hemisphere_z(p, &msh->data->hemisphere); 334 335 /* Transform the position in object space */ 336 d33_muld3(p, msh->transform, p); 337 d3_add(p, p, msh->transform + 9); 338 339 f3_set_d3(pos, p); 340 } 341 342 static FINLINE int 343 aabb_is_degenerated(const double lower[2], const double upper[2]) 344 { 345 ASSERT(lower && upper); 346 return lower[0] >= upper[0] || lower[1] >= upper[1]; 347 } 348 349 static void 350 carvings_compute_aabb 351 (const struct ssol_carving* carvings, 352 const size_t ncarvings, 353 double lower[2], 354 double upper[2]) 355 { 356 size_t icarving; 357 ASSERT(carvings && ncarvings && lower && upper); 358 359 d2_splat(lower, DBL_MAX); 360 d2_splat(upper,-DBL_MAX); 361 362 FOR_EACH(icarving, 0, ncarvings) { 363 size_t ivert; 364 FOR_EACH(ivert, 0, carvings[icarving].nb_vertices) { 365 double pos[2]; 366 /* Discard the polygons to subtract */ 367 if(carvings[icarving].operation == SSOL_SUB) continue; 368 369 carvings[icarving].get(ivert, pos, carvings[icarving].context); 370 d2_min(lower, lower, pos); 371 d2_max(upper, upper, pos); 372 } 373 } 374 } 375 376 static double 377 carvings_compute_radius 378 (const struct ssol_carving* carvings, const size_t ncarvings) 379 { 380 size_t icarving; 381 double r2 = -DBL_MAX; 382 ASSERT(carvings); 383 384 if(!ncarvings) return DBL_MAX; 385 386 FOR_EACH(icarving, 0, ncarvings) { 387 size_t ivert; 388 FOR_EACH(ivert, 0, carvings[icarving].nb_vertices) { 389 double pos[2]; 390 /* Discard the polygons to subtract */ 391 if (carvings[icarving].operation == SSOL_SUB) continue; 392 393 carvings[icarving].get(ivert, pos, carvings[icarving].context); 394 r2 = MMAX(r2, d2_dot(pos, pos)); 395 } 396 } 397 return r2 >= 0 ? sqrt(r2) : DBL_MAX; 398 } 399 400 static res_T 401 build_triangulated_disk 402 (struct darray_double* coords, 403 struct darray_size_t* ids, 404 const double radius, 405 const size_t nsteps) 406 { 407 struct s3dut_mesh_data data; 408 struct s3dut_mesh* mesh = NULL; 409 double *c_ptr = NULL; 410 size_t* i_ptr = NULL; 411 size_t i; 412 res_T res = RES_OK; 413 ASSERT(coords && ids && nsteps && radius > 0); 414 ASSERT(nsteps < UINT_MAX); 415 416 s3dut_create_hemisphere 417 (coords->allocator, radius, (unsigned)nsteps, (unsigned)nsteps, &mesh); 418 if (res != RES_OK) { 419 fprintf(stderr, "Could not create the hemisphere 3D data.\n"); 420 goto error; 421 } 422 423 S3DUT(mesh_get_data(mesh, &data)); 424 if (!data.nprimitives || !data.nvertices) { 425 res = RES_BAD_ARG; 426 goto error; 427 } 428 429 darray_double_clear(coords); 430 darray_size_t_clear(ids); 431 432 /* Reserve the memory space for the plane vertices */ 433 res = darray_double_resize(coords, data.nvertices * 2/*#coords per vertex*/); 434 if (res != RES_OK) goto error; 435 436 /* Reserve the memory space for the plane indices */ 437 res = darray_size_t_resize(ids, data.nprimitives * 3/*#ids per triangle*/); 438 if (res != RES_OK) goto error; 439 440 c_ptr = darray_double_data_get(coords); 441 FOR_EACH(i, 0, data.nvertices) { 442 d2_set(c_ptr + i * 2, data.positions + i * 3); /* don't get z */ 443 } 444 i_ptr = darray_size_t_data_get(ids); 445 FOR_EACH(i, 0, data.nprimitives * 3) i_ptr[i] = data.indices[i]; 446 447 exit: 448 if(mesh) S3DUT(mesh_ref_put(mesh)); 449 return res; 450 error: 451 darray_double_clear(coords); 452 darray_size_t_clear(ids); 453 goto exit; 454 } 455 456 static res_T 457 build_triangulated_plane 458 (struct darray_double* coords, 459 struct darray_size_t* ids, 460 const double lower[2], 461 const double upper[2], 462 const size_t nsteps) 463 { 464 size_t nsteps2[2]; 465 size_t nverts[2]; 466 size_t ix, iy; 467 double size[2]; 468 double size_min; 469 double delta; 470 res_T res = RES_OK; 471 ASSERT(coords && ids && lower && upper && nsteps); 472 ASSERT(!aabb_is_degenerated(lower, upper)); 473 474 darray_double_clear(coords); 475 darray_size_t_clear(ids); 476 477 d2_sub(size, upper, lower); 478 size_min = MMIN(size[0], size[1]); 479 480 if(eq_eps(size_min, 0, 1.e-6)) { 481 res = RES_BAD_ARG; 482 goto error; 483 } 484 485 delta = size_min / (double)nsteps; 486 nsteps2[0] = (size_t)ceil(size[0] / delta); 487 nsteps2[1] = (size_t)ceil(size[1] / delta); 488 nverts[0] = nsteps2[0] + 1; 489 nverts[1] = nsteps2[1] + 1; 490 491 /* Reserve the memory space for the plane vertices */ 492 res = darray_double_reserve(coords, 493 nverts[0]*nverts[1]*2/*#coords per vertex*/); 494 if(res != RES_OK) goto error; 495 496 /* Reserve the memory space for the plane indices */ 497 res = darray_size_t_reserve(ids, 498 nsteps2[0] * nsteps2[1] * 2/*#triangle per step*/*3/*#ids per triangle*/); 499 if(res != RES_OK) goto error; 500 501 /* Setup the plane vertices */ 502 FOR_EACH(ix, 0, nverts[0]) { 503 double x = lower[0] + (double)ix*delta; 504 x = MMIN(x, upper[0]); 505 FOR_EACH(iy, 0, nverts[1]) { 506 double y = lower[1] + (double)iy*delta; 507 y = MMIN(y, upper[1]); 508 darray_double_push_back(coords, &x); 509 darray_double_push_back(coords, &y); 510 } 511 } 512 513 /* Setup the plane indices */ 514 FOR_EACH(ix, 0, nsteps2[0]) { 515 const size_t offset0 = ix*nverts[1]; 516 const size_t offset1 = (ix+1)*nverts[1]; 517 518 FOR_EACH(iy, 0, nsteps2[1]) { 519 const size_t id0 = offset0 + iy; 520 const size_t id1 = offset1 + iy; 521 const size_t id2 = offset0 + iy + 1; 522 const size_t id3 = offset1 + iy + 1; 523 524 darray_size_t_push_back(ids, &id0); 525 darray_size_t_push_back(ids, &id3); 526 darray_size_t_push_back(ids, &id1); 527 528 darray_size_t_push_back(ids, &id0); 529 darray_size_t_push_back(ids, &id2); 530 darray_size_t_push_back(ids, &id3); 531 } 532 } 533 534 exit: 535 return res; 536 error: 537 darray_double_clear(coords); 538 darray_size_t_clear(ids); 539 goto exit; 540 } 541 542 static void 543 get_carving_count 544 (const size_t ic, 545 size_t *count, 546 void* ctx) 547 { 548 const struct ssol_carving* carving = (const struct ssol_carving*)ctx; 549 ASSERT(ic == 0 && count); (void)ic; 550 *count = carving->nb_vertices; 551 } 552 553 static void 554 get_carving_pos 555 (const size_t ic, 556 const size_t iv, 557 double pos[2], 558 void* ctx) 559 { 560 struct ssol_carving* carving = (struct ssol_carving*)ctx; 561 ASSERT(ic == 0); (void)ic; 562 carving->get(iv, pos, carving->context); 563 } 564 565 static res_T 566 clip_triangulated_sheet 567 (struct darray_double* coords, 568 struct darray_size_t* ids, 569 struct scpr_mesh* mesh, 570 const struct ssol_carving* carvings, 571 const size_t ncarvings, 572 struct scpr_device* scpr) 573 { 574 struct mesh_context msh; 575 size_t nverts; 576 size_t ntris; 577 size_t icarving; 578 size_t i; 579 struct scpr_polygon* polygon = NULL; 580 res_T res = RES_OK; 581 ASSERT(coords && ids && carvings && ncarvings); 582 ASSERT(darray_double_size_get(coords) % 2 == 0); 583 ASSERT(darray_size_t_size_get(ids) % 3 == 0); 584 585 nverts = darray_double_size_get(coords)/2; 586 ntris = darray_size_t_size_get(ids)/3; 587 if(!nverts || !ntris) goto exit; 588 589 /* Setup the Star-CliPpeR mesh */ 590 msh.coords = darray_double_cdata_get(coords); 591 msh.ids = darray_size_t_cdata_get(ids); 592 res = scpr_mesh_setup_indexed_vertices 593 (mesh, ntris, mesh_get_ids, nverts, mesh_get_pos, &msh); 594 if(res != RES_OK) goto error; 595 596 /* Apply each carving operation to the Star-CliPpeR mesh */ 597 FOR_EACH(icarving, 0, ncarvings) { 598 enum scpr_operation op = ssol_to_scpr_clip_op(carvings[icarving].operation); 599 600 SCPR(polygon_create(scpr, &polygon)); 601 res = scpr_polygon_setup_indexed_vertices(polygon, 1, 602 get_carving_count, 603 get_carving_pos, 604 (void*)&carvings[icarving]); 605 if(res != RES_OK) goto error; 606 607 res = scpr_mesh_clip(mesh, op, polygon); 608 if(res != RES_OK) goto error; 609 SCPR(polygon_ref_put(polygon)); 610 polygon = NULL; 611 } 612 613 /* Reserve the output index/vertex buffer memory space */ 614 SCPR(mesh_get_vertices_count(mesh, &nverts)); 615 SCPR(mesh_get_triangles_count(mesh, &ntris)); 616 darray_double_clear(coords); 617 darray_size_t_clear(ids); 618 res = darray_double_reserve(coords, nverts*2/*#coords per vertex*/); 619 if(res != RES_OK) goto error; 620 res = darray_size_t_reserve(ids, ntris*3/*#ids per triangle*/); 621 if(res != RES_OK) goto error; 622 623 /* Save the coordinates of the clipped mesh */ 624 FOR_EACH(i, 0, nverts) { 625 double pos[2]; 626 SCPR(mesh_get_position(mesh, i, pos)); 627 darray_double_push_back(coords, pos+0); 628 darray_double_push_back(coords, pos+1); 629 } 630 631 /* Save the indices of the clipped mesh */ 632 FOR_EACH(i, 0, ntris) { 633 size_t tri[3]; 634 SCPR(mesh_get_indices(mesh, i, tri)); 635 darray_size_t_push_back(ids, tri+0); 636 darray_size_t_push_back(ids, tri+1); 637 darray_size_t_push_back(ids, tri+2); 638 } 639 640 exit: 641 return res; 642 error: 643 if(polygon) { 644 SCPR(polygon_ref_put(polygon)); 645 } 646 goto exit; 647 } 648 649 static double 650 mesh_compute_area 651 (const unsigned ntris, 652 void (*get_indices)(const unsigned itri, unsigned ids[3], void* data), 653 const unsigned nverts, 654 void (*get_position)(const unsigned ivert, float position[3], void* data), 655 void* ctx) 656 { 657 unsigned itri; 658 double area = 0; 659 (void)nverts; 660 661 FOR_EACH(itri, 0, ntris) { 662 float v0[3], v1[3], v2[3]; 663 double E0[3], E1[3], N[3]; 664 double V0[3], V1[3], V2[3]; 665 unsigned IDS[3]; 666 667 get_indices(itri, IDS, ctx); 668 ASSERT(IDS[0] < nverts); 669 ASSERT(IDS[1] < nverts); 670 ASSERT(IDS[2] < nverts); 671 672 get_position(IDS[0], v0, ctx); 673 get_position(IDS[1], v1, ctx); 674 get_position(IDS[2], v2, ctx); 675 d3_set_f3(V0, v0); 676 d3_set_f3(V1, v1); 677 d3_set_f3(V2, v2); 678 d3_sub(E0, V1, V0); 679 d3_sub(E1, V2, V0); 680 681 area += d3_len(d3_cross(N, E0, E1)); 682 } 683 return area * 0.5; 684 } 685 686 /* Setup the Star-3D shape of the quadric to ray-trace, i.e. the clipped 2D 687 * profile of the quadric whose vertices are displaced with respect to the 688 * quadric equation */ 689 static res_T 690 quadric_setup_s3d_shape_rt 691 (const struct ssol_shape* shape, 692 const struct darray_double* coords, 693 const struct darray_size_t* ids, 694 const double lower[2], 695 const double upper[2], 696 struct s3d_shape* s3dshape, 697 double* rt_area) 698 { 699 struct quadric_mesh_context ctx; 700 struct s3d_vertex_data vdata[2]; 701 unsigned nverts; 702 unsigned ntris; 703 res_T res; 704 ASSERT(shape && coords && ids && lower && upper && s3dshape && rt_area); 705 ASSERT(darray_double_size_get(coords)%2 == 0); 706 ASSERT(darray_size_t_size_get(ids)%3 == 0); 707 ASSERT(darray_double_size_get(coords)/2 <= UINT_MAX); 708 ASSERT(darray_size_t_size_get(ids)/3 <= UINT_MAX); 709 ASSERT(!aabb_is_degenerated(lower, upper)); 710 711 nverts = (unsigned)darray_double_size_get(coords) / 2/*#coords per vertex*/; 712 ntris = (unsigned)darray_size_t_size_get(ids) / 3/*#ids per triangle*/; 713 ctx.coords = darray_double_cdata_get(coords); 714 ctx.ids = darray_size_t_cdata_get(ids); 715 ctx.transform = shape->transform; 716 d2_set(ctx.lower, lower); 717 d2_set(ctx.upper, upper); 718 719 vdata[0].usage = S3D_POSITION; 720 vdata[0].type = S3D_FLOAT3; 721 vdata[0].get = NULL; 722 723 vdata[1].usage = SSOL_TO_S3D_TEXCOORD; 724 vdata[1].type = S3D_FLOAT2; 725 vdata[1].get = quadric_mesh_get_uv; 726 727 ctx.data = &shape->private_data; 728 ctx.quadric_type = shape->quadric_type; 729 switch (shape->quadric_type) { 730 case SSOL_QUADRIC_PARABOL: 731 vdata[0].get = quadric_mesh_parabol_get_pos; 732 break; 733 case SSOL_QUADRIC_HYPERBOL: 734 vdata[0].get = quadric_mesh_hyperbol_get_pos; 735 break; 736 case SSOL_QUADRIC_PARABOLIC_CYLINDER: 737 vdata[0].get = quadric_mesh_parabolic_cylinder_get_pos; 738 break; 739 case SSOL_QUADRIC_PLANE: 740 vdata[0].get = quadric_mesh_plane_get_pos; 741 break; 742 case SSOL_QUADRIC_HEMISPHERE: 743 vdata[0].get = quadric_mesh_hemisphere_get_pos; 744 break; 745 default: FATAL("Unreachable code.\n"); break; 746 } 747 748 res = s3d_mesh_setup_indexed_vertices 749 (s3dshape, ntris, quadric_mesh_get_ids, nverts, vdata, 2, &ctx); 750 if(res != RES_OK) return res; 751 752 *rt_area = mesh_compute_area 753 (ntris, quadric_mesh_get_ids, nverts, vdata[0].get, &ctx); 754 return RES_OK; 755 } 756 757 /* Setup the Star-3D shape of the quadric to sample, i.e. the clipped 2D 758 * profile of the quadric */ 759 static res_T 760 quadric_setup_s3d_shape_samp 761 (const struct ssol_quadric* quadric, 762 const struct darray_double* coords, 763 const struct darray_size_t* ids, 764 const double lower[2], 765 const double upper[2], 766 struct s3d_shape* shape, 767 double *samp_area) 768 { 769 struct quadric_mesh_context ctx; 770 struct s3d_vertex_data vdata[2]; 771 unsigned nverts; 772 unsigned ntris; 773 res_T res; 774 ASSERT(coords && ids && shape && ids && lower && samp_area); 775 ASSERT(darray_double_size_get(coords)%2 == 0); 776 ASSERT(darray_size_t_size_get(ids)%3 == 0); 777 ASSERT(darray_double_size_get(coords)/2 <= UINT_MAX); 778 ASSERT(darray_size_t_size_get(ids)/3 <= UINT_MAX); 779 ASSERT(!aabb_is_degenerated(lower, upper)); 780 781 nverts = (unsigned)darray_double_size_get(coords) / 2/*#coords per vertex*/; 782 ntris = (unsigned)darray_size_t_size_get(ids) / 3/*#ids per triangle*/; 783 ctx.coords = darray_double_cdata_get(coords); 784 ctx.ids = darray_size_t_cdata_get(ids); 785 ctx.transform = quadric->transform; 786 d2_set(ctx.lower, lower); 787 d2_set(ctx.upper, upper); 788 789 vdata[0].usage = S3D_POSITION; 790 vdata[0].type = S3D_FLOAT3; 791 vdata[0].get = quadric_mesh_plane_get_pos; 792 793 vdata[1].usage = SSOL_TO_S3D_TEXCOORD; 794 vdata[1].type = S3D_FLOAT2; 795 vdata[1].get = quadric_mesh_get_uv; 796 797 res = s3d_mesh_setup_indexed_vertices 798 (shape, ntris, quadric_mesh_get_ids, nverts, vdata, 2, &ctx); 799 if(res != RES_OK) return res; 800 *samp_area = mesh_compute_area 801 (ntris, quadric_mesh_get_ids, nverts, quadric_mesh_plane_get_pos, &ctx); 802 return RES_OK; 803 } 804 805 static res_T 806 shape_create 807 (struct ssol_device* dev, 808 struct ssol_shape** out_shape, 809 enum shape_type type) 810 { 811 struct ssol_shape* shape = NULL; 812 res_T res = RES_OK; 813 814 if(!dev || !out_shape || type >= SHAPE_TYPES_COUNT__) { 815 res = RES_BAD_ARG; 816 goto error; 817 } 818 819 shape = MEM_CALLOC(dev->allocator, 1, sizeof(struct ssol_shape)); 820 if(!shape) { 821 res = RES_MEM_ERR; 822 goto error; 823 } 824 SSOL(device_ref_get(dev)); 825 shape->dev = dev; 826 shape->type = type; 827 ref_init(&shape->ref); 828 829 /* Create the s3d_shape to ray-trace */ 830 res = s3d_shape_create_mesh(dev->s3d, &shape->shape_rt); 831 if(res != RES_OK) goto error; 832 res = s3d_mesh_set_hit_filter_function 833 (shape->shape_rt, hit_filter_function, NULL); 834 if(res != RES_OK) goto error; 835 836 /* Create the s3d_shape to sample */ 837 res = s3d_shape_create_mesh(dev->s3d, &shape->shape_samp); 838 if(res != RES_OK) goto error; 839 840 exit: 841 if(out_shape) *out_shape = shape; 842 return res; 843 error: 844 if(shape) { 845 SSOL(shape_ref_put(shape)); 846 shape = NULL; 847 } 848 goto exit; 849 } 850 851 /* Solve a 2nd degree equation. "hint" is used to select among the 2 solutions 852 * (if applies) the selected solution is then the closest to hint ans is 853 * returned in dist[0]. 854 * If there is a second solution, it is returned in dist[1]. 855 * Returns the number of roots. */ 856 static int 857 solve_second 858 (const double a, 859 const double b, 860 const double c, 861 const double hint, 862 double dist[2]) 863 { 864 ASSERT(dist && hint >= 0); 865 if(a == 0) { 866 if(b != 0) { 867 dist[0] = -c / b; /* Degenerated case: 1st degree only */ 868 return 1; 869 } 870 return 0; /* 0 distance determined */ 871 } else { /* Standard case: 2nd degree */ 872 const double delta = b*b - 4*a*c; 873 874 if(delta == 0) { 875 dist[0] = -b / (2*a); 876 return 1; 877 } else { 878 const double sqrt_delta = sqrt(delta); 879 /* Precise formula */ 880 const double t1 = (-b - copysign(sqrt_delta, b)) / (2*a); 881 const double t2 = c / (a*t1); 882 /* dist[0] is the closest value to hint */ 883 dist[0] = fabs(t1 - hint) < fabs(t2 - hint) ? t1 : t2; 884 dist[1] = fabs(t1 - hint) < fabs(t2 - hint) ? t2 : t1; 885 return 2; 886 } 887 } 888 } 889 890 static FINLINE void 891 quadric_plane_gradient_local(double grad[3]) 892 { 893 ASSERT(grad); 894 grad[0] = 0; 895 grad[1] = 0; 896 grad[2] = 1; 897 } 898 899 static FINLINE void 900 quadric_parabol_gradient_local 901 (const struct priv_parabol_data* quad, 902 const double pt[3], 903 double grad[3]) 904 { 905 ASSERT(quad && pt && grad); 906 grad[0] = -pt[0]; 907 grad[1] = -pt[1]; 908 grad[2] = 2 * quad->focal; 909 } 910 911 static FINLINE void 912 quadric_hyperbol_gradient_local 913 (const struct priv_hyperbol_data* quad, 914 const double pt[3], 915 double grad[3]) 916 { 917 ASSERT(quad && pt && grad); 918 { 919 const double z0 = quad->g_square + quad->abs_b; 920 grad[0] = pt[0]; 921 grad[1] = pt[1]; 922 grad[2] = -(pt[2] + z0 - quad->g_square) * quad->a_square_over_b_square; 923 } 924 } 925 926 static FINLINE void 927 quadric_parabolic_cylinder_gradient_local 928 (const struct priv_pcylinder_data* quad, 929 const double pt[3], 930 double grad[3]) 931 { 932 ASSERT(quad && pt && grad); 933 grad[0] = 0; 934 grad[1] = -pt[1]; 935 grad[2] = 2 * quad->focal; 936 } 937 938 static FINLINE void 939 quadric_hemisphere_gradient_local 940 (const struct priv_hemisphere_data* quad, 941 const double pt[3], 942 double grad[3]) 943 { 944 ASSERT(pt && grad); 945 grad[0] = -pt[0]; 946 grad[1] = -pt[1]; 947 grad[2] = quad->radius - pt[2]; 948 } 949 950 static FINLINE int 951 quadric_plane_intersect_local 952 (const double org[3], 953 const double dir[3], 954 const double hint, 955 double hit_pt[3], 956 double grad[3], 957 double* dist) 958 { 959 /* Define 0 z^2 + z + 0 = 0 */ 960 const double a = 0; 961 const double b = dir[2]; 962 const double c = org[2]; 963 double dst[2]; 964 const int n = solve_second(a, b, c, hint, dst); 965 966 if(!n) return 0; 967 d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); 968 quadric_plane_gradient_local(grad); 969 *dist = *dst; 970 return 1; 971 } 972 973 static FINLINE int 974 quadric_parabol_intersect_local 975 (const struct priv_parabol_data* quad, 976 const double org[3], 977 const double dir[3], 978 const double hint, 979 double hit_pt[3], 980 double grad[3], 981 double* dist) /* in/out: */ 982 { 983 /* Define x^2 + y^2 - 4*focal*z = 0 */ 984 double dst[2]; 985 const double a = dir[0] * dir[0] + dir[1] * dir[1]; 986 const double b = 987 2 * org[0] * dir[0] + 2 * org[1] * dir[1] - 4 * quad->focal * dir[2]; 988 const double c = org[0] * org[0] + org[1] * org[1] - 4 * quad->focal * org[2]; 989 const int n = solve_second(a, b, c, hint, dst); 990 991 if(!n) return 0; 992 d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); 993 quadric_parabol_gradient_local(quad, hit_pt, grad); 994 *dist = *dst; 995 return 1; 996 } 997 998 static FINLINE int 999 quadric_hyperbol_intersect_local 1000 (const struct priv_hyperbol_data* quad, 1001 const double org[3], 1002 const double dir[3], 1003 const double hint, 1004 double hit_pt[3], 1005 double grad[3], 1006 double* dist) 1007 { 1008 double dst[2]; 1009 const double b2 = quad->abs_b * quad->abs_b; 1010 const double b2_a2 = b2 * quad->one_over_a_square; 1011 const double z0 = quad->g_square + quad->abs_b; 1012 const double a = 1013 b2_a2 * (dir[0] * dir[0] + dir[1] * dir[1]) - dir[2] * dir[2]; 1014 const double b = 1015 2 * (b2_a2 * (org[0] * dir[0] + org[1] * dir[1]) 1016 - (org[2] + z0 - quad->g_square) * dir[2]); 1017 const double c = b2_a2 * (org[0] * org[0] + org[1] * org[1]) + b2 1018 - (org[2] + z0 - quad->g_square) * (org[2] + z0 - quad->g_square); 1019 const int n = solve_second(a, b, c, hint, dst); 1020 1021 if(!n) return 0; 1022 d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); 1023 quadric_hyperbol_gradient_local(quad, hit_pt, grad); 1024 *dist = *dst; 1025 return 1; 1026 } 1027 1028 static FINLINE int 1029 quadric_hemisphere_intersect_local 1030 (const struct priv_hemisphere_data* quad, 1031 const double org[3], 1032 const double dir[3], 1033 const double hint, 1034 double hit_pt[3], 1035 double grad[3], 1036 double* dist) 1037 { 1038 double dst[2]; 1039 double z0 = -quad->radius; 1040 const double a = dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]; 1041 const double b = 2 * (org[0] * dir[0] + org[1] * dir[1] + org[2] * dir[2] + z0 * dir[2]); 1042 const double c = 1043 org[0] * org[0] + org[1] * org[1] + org[2] * org[2] - quad->sqr_radius 1044 + 2 * z0 * org[2] + z0 * z0; 1045 const int n = solve_second(a, b, c, hint, dst); 1046 1047 if(!n) return 0; 1048 d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); 1049 quadric_hemisphere_gradient_local(quad, hit_pt, grad); 1050 *dist = *dst; 1051 return 1; 1052 } 1053 1054 static FINLINE int 1055 quadric_parabolic_cylinder_intersect_local 1056 (const struct priv_pcylinder_data* quad, 1057 const double org[3], 1058 const double dir[3], 1059 const double hint, 1060 double hit_pt[3], 1061 double grad[3], 1062 double* dist) 1063 { 1064 double dst[2]; 1065 const double a = dir[1] * dir[1]; 1066 const double b = 2 * org[1] * dir[1] - 4 * quad->focal * dir[2]; 1067 const double c = org[1] * org[1] - 4 * quad->focal * org[2]; 1068 const int n = solve_second(a, b, c, hint, dst); 1069 1070 if(!n) return 0; 1071 d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst)); 1072 quadric_parabolic_cylinder_gradient_local(quad, hit_pt, grad); 1073 *dist = *dst; 1074 return 1; 1075 } 1076 1077 static FINLINE void 1078 punched_shape_set_z_local(const struct ssol_shape* shape, double pt[3]) 1079 { 1080 ASSERT(shape && pt); 1081 ASSERT(shape->type == SHAPE_PUNCHED); 1082 switch (shape->quadric_type) { 1083 case SSOL_QUADRIC_PLANE: 1084 pt[2] = 0; 1085 break; 1086 case SSOL_QUADRIC_PARABOLIC_CYLINDER: 1087 pt[2] = parabolic_cylinder_z(pt, &shape->private_data.pcylinder); 1088 break; 1089 case SSOL_QUADRIC_PARABOL: 1090 pt[2] = parabol_z(pt, &shape->private_data.parabol); 1091 break; 1092 case SSOL_QUADRIC_HYPERBOL: 1093 pt[2] = hyperbol_z(pt, &shape->private_data.hyperbol); 1094 break; 1095 case SSOL_QUADRIC_HEMISPHERE: 1096 pt[2] = hemisphere_z(pt, &shape->private_data.hemisphere); 1097 break; 1098 default: FATAL("Unreachable code\n"); break; 1099 } 1100 } 1101 1102 static FINLINE void 1103 punched_shape_set_normal_local 1104 (const struct ssol_shape* shape, 1105 const double pt[3], 1106 double normal[3]) 1107 { 1108 ASSERT(shape && pt); 1109 ASSERT(shape->type == SHAPE_PUNCHED); 1110 switch (shape->quadric_type) { 1111 case SSOL_QUADRIC_PLANE: 1112 quadric_plane_gradient_local(normal); 1113 break; 1114 case SSOL_QUADRIC_PARABOLIC_CYLINDER: 1115 quadric_parabolic_cylinder_gradient_local 1116 (&shape->private_data.pcylinder, pt, normal); 1117 break; 1118 case SSOL_QUADRIC_PARABOL: 1119 quadric_parabol_gradient_local 1120 (&shape->private_data.parabol, pt, normal); 1121 break; 1122 case SSOL_QUADRIC_HYPERBOL: 1123 quadric_hyperbol_gradient_local 1124 (&shape->private_data.hyperbol, pt, normal); 1125 break; 1126 case SSOL_QUADRIC_HEMISPHERE: 1127 quadric_hemisphere_gradient_local 1128 (&shape->private_data.hemisphere, pt, normal); 1129 break; 1130 default: FATAL("Unreachable code\n"); break; 1131 } 1132 } 1133 1134 int 1135 punched_shape_intersect_local 1136 (const struct ssol_shape* shape, 1137 const double org[3], 1138 const double dir[3], 1139 const double hint, 1140 double pt[3], 1141 double N[3], 1142 double* dist) 1143 { 1144 int hit; 1145 ASSERT(shape && org && dir && hint >= 0 && pt && N && dist); 1146 ASSERT(shape->type == SHAPE_PUNCHED); 1147 ASSERT(dir[0] || dir[1] || dir[2]); 1148 1149 /* Hits on quadrics must be recomputed more accurately */ 1150 switch (shape->quadric_type) { 1151 case SSOL_QUADRIC_PLANE: 1152 hit = quadric_plane_intersect_local(org, dir, hint, pt, N, dist); 1153 break; 1154 case SSOL_QUADRIC_PARABOLIC_CYLINDER: 1155 hit = quadric_parabolic_cylinder_intersect_local 1156 (&shape->private_data.pcylinder, org, dir, hint, pt, N, dist); 1157 break; 1158 case SSOL_QUADRIC_PARABOL: 1159 hit = quadric_parabol_intersect_local 1160 (&shape->private_data.parabol, org, dir, hint, pt, N, dist); 1161 break; 1162 case SSOL_QUADRIC_HYPERBOL: 1163 hit = quadric_hyperbol_intersect_local 1164 (&shape->private_data.hyperbol, org, dir, hint, pt, N, dist); 1165 break; 1166 case SSOL_QUADRIC_HEMISPHERE: 1167 hit = quadric_hemisphere_intersect_local 1168 (&shape->private_data.hemisphere, org, dir, hint, pt, N, dist); 1169 break; 1170 default: FATAL("Unreachable code\n"); break; 1171 } 1172 return hit; 1173 } 1174 1175 static void 1176 shape_release(ref_T* ref) 1177 { 1178 struct ssol_device* dev; 1179 struct ssol_shape* shape = CONTAINER_OF(ref, struct ssol_shape, ref); 1180 ASSERT(ref); 1181 dev = shape->dev; 1182 ASSERT(dev && dev->allocator); 1183 if(shape->shape_rt) S3D(shape_ref_put(shape->shape_rt)); 1184 if(shape->shape_samp) S3D(shape_ref_put(shape->shape_samp)); 1185 MEM_RM(dev->allocator, shape); 1186 SSOL(device_ref_put(dev)); 1187 } 1188 1189 /* Return the parabol discretisation parameter */ 1190 static FINLINE void 1191 priv_parabol_data_setup 1192 (struct priv_parabol_data* data, 1193 const struct ssol_quadric_parabol* parabol) 1194 { 1195 ASSERT(data && parabol); 1196 data->focal = parabol->focal; 1197 data->one_over_4focal = 1 / (4.0 * parabol->focal); 1198 } 1199 1200 static FINLINE void 1201 priv_hyperbol_data_setup 1202 (struct priv_hyperbol_data* data, 1203 const struct ssol_quadric_hyperbol* hyperbol) 1204 { 1205 double g, f, a2; 1206 ASSERT(data && hyperbol); 1207 1208 /* Re-dimensionalize */ 1209 g = hyperbol->real_focal + hyperbol->img_focal; 1210 f = hyperbol->real_focal / g; 1211 a2 = g * g * (f - f * f); 1212 1213 data->g_square = g * 0.5; 1214 data->abs_b = g * fabs(f - 0.5); 1215 data->a_square_over_b_square = a2 / (data->abs_b * data->abs_b); 1216 data->one_over_a_square = 1 / a2; 1217 } 1218 1219 static FINLINE void 1220 priv_parabolic_cylinder_data_setup 1221 (struct priv_pcylinder_data* data, 1222 const struct ssol_quadric_parabolic_cylinder* parabolic_cylinder) 1223 { 1224 ASSERT(data && parabolic_cylinder); 1225 data->focal = parabolic_cylinder->focal; 1226 data->one_over_4focal = 1 / (4.0 * parabolic_cylinder->focal); 1227 } 1228 1229 static FINLINE void 1230 priv_hemisphere_data_setup 1231 (struct priv_hemisphere_data* data, 1232 const struct ssol_quadric_hemisphere* hemisphere) 1233 { 1234 ASSERT(data && hemisphere); 1235 data->radius = hemisphere->radius; 1236 data->sqr_radius = hemisphere->radius * hemisphere->radius; 1237 } 1238 1239 static INLINE void 1240 priv_quadric_data_setup 1241 (union private_data* priv_data, 1242 const struct ssol_quadric* quadric) 1243 { 1244 ASSERT(priv_data && quadric); 1245 switch(quadric->type) { 1246 case SSOL_QUADRIC_PLANE: /* Do nothing */ break; 1247 case SSOL_QUADRIC_PARABOL: 1248 priv_parabol_data_setup 1249 (&priv_data->parabol, &quadric->data.parabol); 1250 break; 1251 case SSOL_QUADRIC_HYPERBOL: 1252 priv_hyperbol_data_setup 1253 (&priv_data->hyperbol, &quadric->data.hyperbol); 1254 break; 1255 case SSOL_QUADRIC_PARABOLIC_CYLINDER: 1256 priv_parabolic_cylinder_data_setup 1257 (&priv_data->pcylinder, &quadric->data.parabolic_cylinder); 1258 break; 1259 case SSOL_QUADRIC_HEMISPHERE: 1260 priv_hemisphere_data_setup 1261 (&priv_data->hemisphere, &quadric->data.hemisphere); 1262 break; 1263 default: FATAL("Unreachable code\n"); break; 1264 } 1265 } 1266 1267 static INLINE size_t 1268 priv_quadric_data_compute_slices_count 1269 (const enum ssol_quadric_type type, 1270 const union private_data* priv_data, 1271 const double lower[2], 1272 const double upper[2]) 1273 { 1274 size_t nslices; 1275 double max_z; 1276 ASSERT(priv_data && lower && upper); 1277 1278 switch(type) { 1279 case SSOL_QUADRIC_PLANE: nslices = 1; break; 1280 case SSOL_QUADRIC_PARABOL: 1281 max_z = MMAX 1282 (parabol_z(lower, &priv_data->parabol), 1283 parabol_z(upper, &priv_data->parabol)); 1284 nslices = MMIN(50, (size_t)(3 + sqrt(max_z) * 6)); 1285 break; 1286 case SSOL_QUADRIC_HYPERBOL: 1287 max_z = MMAX 1288 (hyperbol_z(lower, &priv_data->hyperbol), 1289 hyperbol_z(upper, &priv_data->hyperbol)); 1290 nslices = MMIN(50, (size_t)(3 + sqrt(max_z) * 6)); 1291 break; 1292 case SSOL_QUADRIC_PARABOLIC_CYLINDER: 1293 max_z = MMAX 1294 (parabolic_cylinder_z(lower, &priv_data->pcylinder), 1295 parabolic_cylinder_z(upper, &priv_data->pcylinder)); 1296 nslices = MMIN(50, (size_t)(3 + sqrt(max_z) * 6)); 1297 break; 1298 default: FATAL("Unreachable code\n"); break; 1299 } 1300 return nslices; 1301 } 1302 1303 static INLINE size_t 1304 hemisphere_compute_slices_count 1305 (const struct priv_hemisphere_data* hemisphere, const double radius) 1306 { 1307 size_t nslices; 1308 ASSERT(hemisphere && radius > 0 && hemisphere->radius >= radius); 1309 /* default ranging from 5 to 16 */ 1310 nslices = (size_t)(5.5 + 11 * radius / hemisphere->radius); 1311 return nslices; 1312 } 1313 1314 /******************************************************************************* 1315 * Local functions 1316 ******************************************************************************/ 1317 void 1318 punched_shape_project_point 1319 (struct ssol_shape* shape, 1320 const double transform[12], /* Shape to world space transformation */ 1321 const double pos[3], /* World space position near of the quadric */ 1322 double pos_quadric[3], /* World space position onto the quadric */ 1323 double N_quadric[3]) /* World space normal onto the quadric */ 1324 { 1325 double R[9]; /* Quadric to world rotation matrix */ 1326 double R_invtrans[9]; /* Inverse transpose of R */ 1327 double T[3]; /* Quadric to world translation vector */ 1328 double T_inv[3]; /* Inverse of T */ 1329 double pos_local[3]; 1330 double N_local[3]; 1331 ASSERT(shape && transform && pos && pos_quadric && N_quadric); 1332 ASSERT(shape->type == SHAPE_PUNCHED); 1333 1334 /* Compute world<->quadric space transformations */ 1335 d33_muld33(R, transform, shape->transform); 1336 d33_muld3(T, transform, shape->transform+9); 1337 d3_add(T, T, transform + 9); 1338 d33_invtrans(R_invtrans, R); 1339 d3_minus(T_inv, T); 1340 1341 /* Transform pos in quadric space */ 1342 d3_add(pos_local, pos, T_inv); 1343 d3_muld33(pos_local, pos_local, R_invtrans); 1344 1345 /* Project pos_local onto the quadric and compute its associated normal */ 1346 punched_shape_set_z_local(shape, pos_local); 1347 punched_shape_set_normal_local(shape, pos_local, N_local); 1348 1349 /* Transform the local position in world space */ 1350 d33_muld3(pos_quadric, R, pos_local); 1351 d3_add(pos_quadric, pos_quadric, T); 1352 1353 /* Transform the quadric normal in world space */ 1354 d33_muld3(N_quadric, R_invtrans, N_local); 1355 d3_normalize(N_quadric, N_quadric); 1356 } 1357 1358 double 1359 shape_trace_ray 1360 (struct ssol_shape* shape, 1361 const double transform[12], /* Shape to world space transformation */ 1362 const double org[3], /* World space position near of the ray origin */ 1363 const double dir[3], /* World space ray direction */ 1364 const double hint_dst, /* Hint on the hit distance */ 1365 double N_shape[3], /* World space normal onto the shape */ 1366 intersect_local_fn local) /* the intersection function for this shape */ 1367 { 1368 double R[9]; /* Quadric to world rotation matrix */ 1369 double R_invtrans[9]; /* Inverse transpose of R */ 1370 double T[3]; /* Quadric to world translation vector */ 1371 double T_inv[3]; /* Inverse of T */ 1372 double dir_local[3]; 1373 double org_local[3]; 1374 double hit_local[3]; 1375 double N_local[3]; 1376 double dst; /* Hit distance */ 1377 int valid; 1378 ASSERT(shape && transform && org && N_shape); 1379 1380 /* Compute world<->quadric space transformations */ 1381 d33_muld33(R, transform, shape->transform); 1382 d33_muld3(T, transform, shape->transform+9); 1383 d3_add(T, T, transform + 9); 1384 d33_invtrans(R_invtrans, R); 1385 d3_minus(T_inv, T); 1386 1387 /* Transform pos in quadric space */ 1388 d3_add(org_local, org, T_inv); 1389 d3_muld33(org_local, org_local, R_invtrans); 1390 1391 /* Transform dir in quadric space */ 1392 d3_muld33(dir_local, dir, R_invtrans); 1393 1394 /* Project pos_local onto the shape and compute its associated normal */ 1395 valid = local 1396 (shape, org_local, dir_local, hint_dst, hit_local, N_local, &dst); 1397 if(!valid) return INF; 1398 1399 /* Transform the shape normal in world space */ 1400 d33_muld3(N_shape, R_invtrans, N_local); 1401 d3_normalize(N_shape, N_shape); 1402 return dst; 1403 } 1404 1405 res_T 1406 shape_fetched_raw_vertex_attrib 1407 (const struct ssol_shape* shape, 1408 const unsigned ivert, 1409 const enum ssol_attrib_usage usage, 1410 double value[3]) 1411 { 1412 struct s3d_attrib s3d_attr; 1413 enum s3d_attrib_usage s3d_usage; 1414 res_T res = RES_OK; 1415 1416 ASSERT(shape && value); 1417 s3d_usage = ssol_to_s3d_attrib_usage(usage); 1418 1419 res = s3d_mesh_get_vertex_attrib 1420 (shape->shape_rt, ivert, s3d_usage, &s3d_attr); 1421 if(res != RES_OK) return res; 1422 1423 d3_splat(value, 1); 1424 1425 switch(s3d_attr.type) { 1426 case S3D_FLOAT3: value[2] = (double)s3d_attr.value[2]; FALLTHROUGH; 1427 case S3D_FLOAT2: value[1] = (double)s3d_attr.value[1]; FALLTHROUGH; 1428 case S3D_FLOAT: value[0] = (double)s3d_attr.value[0]; break; 1429 default: FATAL("Unexpected vertex attrib type\n"); break; 1430 } 1431 return RES_OK; 1432 } 1433 1434 /******************************************************************************* 1435 * Exported ssol_shape functions 1436 ******************************************************************************/ 1437 res_T 1438 ssol_shape_create_mesh 1439 (struct ssol_device* dev, 1440 struct ssol_shape** out_shape) 1441 { 1442 return shape_create(dev, out_shape, SHAPE_MESH); 1443 } 1444 1445 res_T 1446 ssol_shape_create_punched_surface 1447 (struct ssol_device* dev, 1448 struct ssol_shape** out_shape) 1449 { 1450 return shape_create(dev, out_shape, SHAPE_PUNCHED); 1451 } 1452 1453 res_T 1454 ssol_shape_ref_get(struct ssol_shape* shape) 1455 { 1456 if(!shape) return RES_BAD_ARG; 1457 ref_get(&shape->ref); 1458 return RES_OK; 1459 } 1460 1461 res_T 1462 ssol_shape_ref_put(struct ssol_shape* shape) 1463 { 1464 if(!shape) return RES_BAD_ARG; 1465 ref_put(&shape->ref, shape_release); 1466 return RES_OK; 1467 } 1468 1469 res_T 1470 ssol_shape_get_vertices_count 1471 (const struct ssol_shape* shape, unsigned* nverts) 1472 { 1473 if(!shape || !nverts) return RES_BAD_ARG; 1474 return s3d_mesh_get_vertices_count(shape->shape_rt, nverts); 1475 } 1476 1477 res_T 1478 ssol_shape_get_vertex_attrib 1479 (const struct ssol_shape* shape, 1480 const unsigned ivert, 1481 const enum ssol_attrib_usage usage, 1482 double value[]) 1483 { 1484 res_T res = RES_OK; 1485 if(!shape || (unsigned)usage >= SSOL_ATTRIBS_COUNT__ || !value) 1486 return RES_BAD_ARG; 1487 1488 res = shape_fetched_raw_vertex_attrib(shape, ivert, usage, value); 1489 if(res != RES_OK) return res; 1490 1491 /* Transform the fetch attrib */ 1492 if(shape->type == SHAPE_PUNCHED) { 1493 if(usage == SSOL_POSITION) { 1494 d33_muld3(value, shape->transform, value); 1495 d3_add(value, shape->transform + 9, value); 1496 } else if(usage == SSOL_NORMAL) { 1497 double R_invtrans[9]; 1498 d33_invtrans(R_invtrans, shape->transform); 1499 d33_muld3(value, R_invtrans, value); 1500 } 1501 } 1502 return RES_OK; 1503 } 1504 1505 res_T 1506 ssol_shape_get_triangles_count(const struct ssol_shape* shape, unsigned* ntris) 1507 { 1508 if(!shape || !ntris) return RES_BAD_ARG; 1509 return s3d_mesh_get_triangles_count(shape->shape_rt, ntris); 1510 } 1511 1512 res_T 1513 ssol_shape_get_triangle_indices 1514 (const struct ssol_shape* shape, const unsigned itri, unsigned ids[3]) 1515 { 1516 if(!shape || !ids) return RES_BAD_ARG; 1517 return s3d_mesh_get_triangle_indices(shape->shape_rt, itri, ids); 1518 } 1519 1520 res_T 1521 ssol_punched_surface_setup 1522 (struct ssol_shape* shape, 1523 const struct ssol_punched_surface* psurf) 1524 { 1525 double lower[2], upper[2]; /* Carvings AABB */ 1526 double radius = -1; 1527 struct darray_double coords; 1528 struct darray_size_t ids; 1529 size_t nslices; 1530 res_T res = RES_OK; 1531 1532 darray_double_init(shape->dev->allocator, &coords); 1533 darray_size_t_init(shape->dev->allocator, &ids); 1534 1535 if(!check_shape(shape) 1536 || !check_punched_surface(psurf) 1537 || shape->type != SHAPE_PUNCHED) { 1538 res = RES_BAD_ARG; 1539 goto error; 1540 } 1541 1542 /* Save quadric for further object instancing */ 1543 d33_set(shape->transform, psurf->quadric->transform); 1544 d3_set(shape->transform+9, psurf->quadric->transform+9); 1545 shape->quadric_type = psurf->quadric->type; 1546 1547 if(psurf->quadric->type == SSOL_QUADRIC_HEMISPHERE) { 1548 radius = carvings_compute_radius(psurf->carvings, psurf->nb_carvings); 1549 radius = MMIN(radius, psurf->quadric->data.hemisphere.radius); 1550 lower[0] = lower[1] = -radius; 1551 upper[0] = upper[1] = +radius; 1552 } else { 1553 carvings_compute_aabb(psurf->carvings, psurf->nb_carvings, lower, upper); 1554 if(aabb_is_degenerated(lower, upper)) { 1555 log_error(shape->dev, 1556 "%s: infinite or null punched surface.\n", 1557 FUNC_NAME); 1558 res = RES_BAD_ARG; 1559 goto error; 1560 } 1561 } 1562 1563 /* Setup internal data */ 1564 priv_quadric_data_setup(&shape->private_data, psurf->quadric); 1565 1566 /* Define the #slices of the discreet quadric */ 1567 if(psurf->quadric->slices_count_hint != SIZE_MAX) { 1568 nslices = psurf->quadric->slices_count_hint; 1569 } else if(psurf->quadric->type == SSOL_QUADRIC_HEMISPHERE) { 1570 nslices = hemisphere_compute_slices_count 1571 (&shape->private_data.hemisphere, radius); 1572 } else { 1573 nslices = priv_quadric_data_compute_slices_count 1574 (shape->quadric_type, &shape->private_data, lower, upper); 1575 } 1576 1577 /* Build the quadric mesh */ 1578 if(psurf->quadric->type == SSOL_QUADRIC_HEMISPHERE) { 1579 res = build_triangulated_disk(&coords, &ids, radius, nslices); 1580 } else { 1581 res = build_triangulated_plane(&coords, &ids, lower, upper, nslices); 1582 } 1583 if(res != RES_OK) goto error; 1584 1585 /* Clip the quadric mesh if necessary */ 1586 if(psurf->nb_carvings) { 1587 res = clip_triangulated_sheet 1588 (&coords, &ids, shape->dev->scpr_mesh, psurf->carvings, psurf->nb_carvings, 1589 shape->dev->scpr_device); 1590 if(res != RES_OK) goto error; 1591 } 1592 1593 /* Setup the Star-3D shape to ray-trace */ 1594 res = quadric_setup_s3d_shape_rt(shape, &coords, &ids, lower, upper, 1595 shape->shape_rt, &shape->shape_rt_area); 1596 if(res != RES_OK) goto error; 1597 1598 /* Setup the Star-3D shape to sample */ 1599 res = quadric_setup_s3d_shape_samp(psurf->quadric, &coords, &ids, lower, 1600 upper, shape->shape_samp, &shape->shape_samp_area); 1601 if(res != RES_OK) goto error; 1602 1603 exit: 1604 darray_double_release(&coords); 1605 darray_size_t_release(&ids); 1606 return res; 1607 error: 1608 goto exit; 1609 } 1610 1611 res_T 1612 ssol_mesh_setup 1613 (struct ssol_shape* shape, 1614 const unsigned ntris, 1615 void(*get_indices)(const unsigned itri, unsigned ids[3], void* ctx), 1616 const unsigned nverts, 1617 const struct ssol_vertex_data attribs [], 1618 const unsigned nattribs, 1619 void* data) 1620 { 1621 struct s3d_vertex_data attrs[SSOL_ATTRIBS_COUNT__]; 1622 void (*get_position)(const unsigned ivert, float position[3], void* data) = NULL; 1623 res_T res = RES_OK; 1624 unsigned i; 1625 1626 if(!check_shape(shape) 1627 || shape->type != SHAPE_MESH 1628 || !get_indices 1629 || !ntris 1630 || !nverts 1631 || !attribs 1632 || !nattribs) { 1633 res = RES_BAD_ARG; 1634 goto error; 1635 } 1636 1637 if(nattribs > SSOL_ATTRIBS_COUNT__) { 1638 res = RES_MEM_ERR; 1639 goto error; 1640 } 1641 1642 FOR_EACH(i, 0, nattribs) { 1643 attrs[i].get = attribs[i].get; 1644 switch (attribs[i].usage) { 1645 case SSOL_POSITION: 1646 attrs[i].usage = SSOL_TO_S3D_POSITION; 1647 attrs[i].type = S3D_FLOAT3; 1648 ASSERT(!get_position); 1649 get_position = attrs[i].get; 1650 break; 1651 case SSOL_NORMAL: 1652 attrs[i].usage = SSOL_TO_S3D_NORMAL; 1653 attrs[i].type = S3D_FLOAT3; 1654 break; 1655 case SSOL_TEXCOORD: 1656 attrs[i].usage = SSOL_TO_S3D_TEXCOORD; 1657 attrs[i].type = S3D_FLOAT2; 1658 break; 1659 default: FATAL("Unreachable code.\n"); break; 1660 } 1661 } 1662 ASSERT(get_position); 1663 1664 res = s3d_mesh_setup_indexed_vertices 1665 (shape->shape_rt, ntris, get_indices, nverts, attrs, nattribs, data); 1666 if(res != RES_OK) goto error; 1667 shape->shape_rt_area = 1668 mesh_compute_area(ntris, get_indices, nverts, get_position, data); 1669 1670 /* The Star-3D shape to sample is the same that the one to ray-trace */ 1671 res = s3d_mesh_copy(shape->shape_rt, shape->shape_samp); 1672 if(res != RES_OK) goto error; 1673 shape->shape_samp_area = shape->shape_rt_area; 1674 1675 exit: 1676 return res; 1677 error: 1678 goto exit; 1679 }