schiff_geometry.c (35225B)
1 /* Copyright (C) 2015, 2016, 2026 Centre National de la Recherche Scientifique 2 * Copyright (C) 2026 Clermont Auvergne INP 3 * Copyright (C) 2026 Institut Mines Télécom Albi-Carmaux 4 * Copyright (C) 2017, 2019-2021, 2026 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2026 Université de Lorraine 6 * Copyright (C) 2026 Université de Toulouse 7 * 8 * This program is free software: you can redistribute it and/or modify 9 * it under the terms of the GNU General Public License as published by 10 * the Free Software Foundation, either version 3 of the License, or 11 * (at your option) any later version. 12 * 13 * This program is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 * GNU General Public License for more details. 17 * 18 * You should have received a copy of the GNU General Public License 19 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 20 21 #include "schiff_geometry.h" 22 #include "schiff_mesh.h" 23 #include "schiff_optical_properties.h" 24 25 #include <rsys/algorithm.h> 26 #include <rsys/stretchy_array.h> 27 28 #include <star/s3d.h> 29 #include <star/ssp.h> 30 #include <star/sschiff.h> 31 32 /* 3D Context of an ellipsoid */ 33 struct ellipsoid_context { 34 /* Precomputed values used to define the radius of the ellipsoid for a 35 * {theta, phi} position : 36 * r = (a*b*c) 37 * / sqrt(b^2*c^2*cos(theta)^2*cos(phi)^2 38 * + a^2*c^2*sin(theta)^2*cos(phi)^2 39 * + a^2*b^2*sin(phi)^2) */ 40 double abc; /* a*b*c */ 41 double b2c2; /* b^2*c^2 */ 42 double a2c2; /* a^2*c^2 */ 43 double a2b2; /* a^2*b^2 */ 44 }; 45 46 /* 3D Context of the sphere geometry */ 47 struct sphere_context { 48 float radius; /* Sphere radius */ 49 }; 50 51 /* 3D Context of the cylinder geometry */ 52 struct cylinder_context { 53 float radius; 54 float height; 55 }; 56 57 /* 3D Context of the supershape geometry */ 58 struct supershape_context { 59 double formulas[2][6]; 60 }; 61 62 /* Distribution context of a cylinder geometry */ 63 struct cylinder_distribution_context { 64 double log_mean_radius; /* Log of the mean radius */ 65 double log_sigma; /* Log of the sigma argument of the lognormal distribution*/ 66 double aspect_ratio; /* aspect ratio of the cylinder distribution */ 67 }; 68 69 /* 3D context of a generic geometry */ 70 struct mesh_context { 71 const struct schiff_mesh* mesh; /* Triangular mesh of the geometry */ 72 const struct schiff_geometry* geometry; 73 enum schiff_geometry_type type; 74 union { 75 struct ellipsoid_context ellipsoid; 76 struct cylinder_context cylinder; 77 struct sphere_context sphere; 78 struct supershape_context supershape; 79 } data; 80 }; 81 82 /* Specific mesh context for the helical pipe */ 83 struct helical_pipe_mesh_context { 84 const unsigned* indices; 85 const float* vertices; 86 }; 87 88 89 struct shape { 90 const struct schiff_geometry* geometry; /* Pointer onto the geometry */ 91 struct s3d_shape* shape; /* May be NULL => No volume distribution */ 92 struct schiff_mesh mesh; /* Triangular mesh of the shape */ 93 double volume; /* Volume of the shape */ 94 }; 95 96 /* Distribution context of a generic geometry */ 97 struct geometry_distribution_context { 98 struct schiff_optical_properties* properties; /* Per wavelength properties */ 99 struct shape* shapes; /* List of shapes */ 100 struct ssp_ranst_discrete* ran_geometries; /* Geometries random variates */ 101 }; 102 103 /******************************************************************************* 104 * Helper functions 105 ******************************************************************************/ 106 static int 107 cmp_double(const void* op0, const void* op1) 108 { 109 const double* a = (const double*)op0; 110 const double* b = (const double*)op1; 111 return *a < *b ? -1 : (*a > *b ? 1 : 0); 112 } 113 114 static double 115 histogram_sample 116 (const double* entries, 117 const size_t nentries, 118 const double lower, 119 const double upper, 120 const double u) 121 { 122 const double* find; 123 double step; 124 double from, to; 125 double v; 126 ASSERT(entries && nentries && lower < upper && u >= 0 && u < 1.0); 127 128 /* Search for the first entry that is gequal than the canonical variable u */ 129 find = search_lower_bound(&u, entries, nentries, sizeof(double), cmp_double); 130 131 if(!find) /* Handle numerical issue */ 132 find = entries + (nentries - 1); 133 134 step = (upper - lower) / (double)nentries; /* Size of an histogram interval */ 135 from = lower + (double)(find - entries) * step; 136 to = from + step; 137 138 /* Linear interpolation */ 139 v = CLAMP((u - from) / step, 0.0, 1.0); 140 if(eq_eps(v, 0, 1.e-6)) { 141 return from; 142 } else if(eq_eps(v, 1, 1.e-6)) { 143 return to; 144 } else { 145 return v*to + (1 - v)*from; 146 } 147 } 148 149 static FINLINE double 150 ran_gaussian_truncated 151 (struct ssp_rng* rng, 152 const double mu, 153 const double sigma, 154 const double range[2]) 155 { 156 double val; 157 ASSERT(rng && mu > 0 && mu >= range[0] && mu <= range[1]); 158 do { 159 val = ssp_ran_gaussian(rng, mu, sigma); 160 } while(val < range[0] || val > range[1] ); 161 return val; 162 } 163 164 static INLINE double 165 eval_param(const struct schiff_param* param, struct ssp_rng* rng) 166 { 167 double val = 0.0; 168 ASSERT(param); 169 170 switch(param->distribution) { 171 case SCHIFF_PARAM_CONSTANT: 172 val = param->data.constant; 173 break; 174 case SCHIFF_PARAM_LOGNORMAL: 175 val = ssp_ran_lognormal(rng, log(param->data.lognormal.mu), 176 log(param->data.lognormal.sigma)); 177 break; 178 case SCHIFF_PARAM_GAUSSIAN: 179 val = ran_gaussian_truncated(rng, param->data.gaussian.mu, 180 param->data.gaussian.sigma, param->data.gaussian.range); 181 break; 182 case SCHIFF_PARAM_HISTOGRAM: 183 val = histogram_sample 184 (param->data.histogram.entries, 185 sa_size(param->data.histogram.entries), 186 param->data.histogram.lower, 187 param->data.histogram.upper, 188 ssp_rng_canonical(rng)); 189 break; 190 default: FATAL("Unreachable code.\n"); break; 191 } 192 return val; 193 } 194 195 static void 196 get_material_property 197 (void* mtl, 198 const double wavelength, 199 struct sschiff_material_properties* props) 200 { 201 struct schiff_optical_properties* properties = mtl; 202 schiff_optical_properties_fetch(properties, wavelength, props); 203 } 204 205 static INLINE void 206 geometry_get_indices(const unsigned itri, unsigned ids[3], void* ctx) 207 { 208 struct mesh_context* mesh_ctx = ctx; 209 ASSERT(ctx); 210 schiff_mesh_get_indices(mesh_ctx->mesh, itri, ids); 211 } 212 213 static INLINE void 214 geometry_get_position(const unsigned ivert, float vertex[3], void* ctx) 215 { 216 struct mesh_context* mesh_ctx = ctx; 217 ASSERT(ctx); 218 schiff_mesh_get_cartesian_position(mesh_ctx->mesh, ivert, vertex); 219 } 220 221 static void 222 ellipsoid_get_position(const unsigned ivert, float vertex[3], void* ctx) 223 { 224 struct mesh_context* mesh_ctx = ctx; 225 struct sin_cos angles[2]; 226 double sqr_cos_theta; 227 double sqr_cos_phi; 228 double sqr_sin_theta; 229 double sqr_sin_phi; 230 double radius; 231 double tmp; 232 ASSERT(mesh_ctx && mesh_ctx->type == SCHIFF_ELLIPSOID); 233 234 schiff_mesh_get_polar_position(mesh_ctx->mesh, ivert, angles); 235 sqr_cos_theta = angles[0].cosine * angles[0].cosine; 236 sqr_sin_theta = angles[0].sinus * angles[0].sinus; 237 sqr_cos_phi = angles[1].cosine * angles[1].cosine; 238 sqr_sin_phi = angles[1].sinus * angles[1].sinus; 239 240 tmp = mesh_ctx->data.ellipsoid.b2c2 * sqr_cos_theta * sqr_cos_phi 241 + mesh_ctx->data.ellipsoid.a2c2 * sqr_sin_theta * sqr_cos_phi 242 + mesh_ctx->data.ellipsoid.a2b2 * sqr_sin_phi; 243 radius = mesh_ctx->data.ellipsoid.abc / sqrt(tmp); 244 245 vertex[0] = (float)(angles[0].cosine * angles[1].cosine * radius); 246 vertex[1] = (float)(angles[0].sinus * angles[1].cosine * radius); 247 vertex[2] = (float)(angles[1].sinus * radius); 248 } 249 250 static void 251 cylinder_get_position(const unsigned ivert, float vertex[3], void* ctx) 252 { 253 struct mesh_context* mesh_ctx = ctx; 254 ASSERT(mesh_ctx && mesh_ctx->type == SCHIFF_CYLINDER); 255 schiff_mesh_get_cartesian_position(mesh_ctx->mesh, ivert, vertex); 256 vertex[0] *= mesh_ctx->data.cylinder.radius; 257 vertex[1] *= mesh_ctx->data.cylinder.radius; 258 vertex[2] *= mesh_ctx->data.cylinder.height; 259 } 260 261 static void 262 helical_pipe_get_indices(const unsigned itri, unsigned ids[3], void* ctx) 263 { 264 size_t i = itri * 3/*#indices per triangle*/; 265 const struct helical_pipe_mesh_context* mesh_ctx = ctx; 266 ASSERT(mesh_ctx); 267 ids[0] = mesh_ctx->indices[i + 0]; 268 ids[1] = mesh_ctx->indices[i + 1]; 269 ids[2] = mesh_ctx->indices[i + 2]; 270 } 271 272 static void 273 helical_pipe_get_position(const unsigned ivert, float vertex[3], void* ctx) 274 { 275 const size_t i = ivert * 3/*#coords*/; 276 const struct helical_pipe_mesh_context* mesh_ctx = ctx; 277 ASSERT(mesh_ctx); 278 vertex[0] = mesh_ctx->vertices[i + 0]; 279 vertex[1] = mesh_ctx->vertices[i + 1]; 280 vertex[2] = mesh_ctx->vertices[i + 2]; 281 } 282 283 static void 284 sphere_get_position(const unsigned ivert, float vertex[3], void* ctx) 285 { 286 struct mesh_context* mesh_ctx = ctx; 287 ASSERT(mesh_ctx && mesh_ctx->type == SCHIFF_SPHERE); 288 schiff_mesh_get_cartesian_position(mesh_ctx->mesh, ivert, vertex); 289 vertex[0] *= mesh_ctx->data.sphere.radius; 290 vertex[1] *= mesh_ctx->data.sphere.radius; 291 vertex[2] *= mesh_ctx->data.sphere.radius; 292 } 293 294 static void 295 supershape_get_position(const unsigned ivert, float vertex[3], void* ctx) 296 { 297 struct mesh_context* mesh_ctx = ctx; 298 struct sin_cos angles[2]; 299 double uv[2]; 300 int iform; 301 ASSERT(mesh_ctx && mesh_ctx->type == SCHIFF_SUPERSHAPE); 302 303 schiff_mesh_get_polar_position(mesh_ctx->mesh, ivert, angles); 304 305 FOR_EACH(iform, 0, 2) { 306 double m, k, g; 307 double* form = mesh_ctx->data.supershape.formulas[iform]; 308 309 m = cos(form[M] * angles[iform].angle / 4.0); 310 m = fabs(m) / form[A]; 311 k = sin(form[M] * angles[iform].angle / 4.0); 312 k = fabs(k) / form[B]; 313 g = pow(m, form[N1]) + pow(k, form[N2]); 314 uv[iform] = pow(g, (-1.0/form[N0])); 315 } 316 317 vertex[0] = (float)(uv[0] * angles[0].cosine * uv[1] * angles[1].cosine); 318 vertex[1] = (float)(uv[0] * angles[0].sinus * uv[1] * angles[1].cosine); 319 vertex[2] = (float)(uv[1] * angles[1].sinus); 320 } 321 322 static res_T 323 compute_s3d_shape_volume 324 (struct s3d_device* s3d, struct s3d_shape* shape, double* out_volume) 325 { 326 struct s3d_scene* scn = NULL; 327 struct s3d_scene_view* view = NULL; 328 float volume = 0.f; 329 res_T res = RES_OK; 330 ASSERT(s3d && shape && out_volume); 331 332 if(RES_OK != (res = s3d_scene_create(s3d, &scn))) goto error; 333 if(RES_OK != (res = s3d_scene_attach_shape(scn, shape))) goto error; 334 if(RES_OK != (res = s3d_scene_view_create(scn, 0, &view))) goto error; 335 if(RES_OK != (res = s3d_scene_view_compute_volume(view, &volume))) goto error; 336 337 exit: 338 if(scn) S3D(scene_ref_put(scn)); 339 if(view) S3D(scene_view_ref_put(view)); 340 /* The volume may be negative if the faces are not correctly oriented */ 341 *out_volume = absf(volume); 342 return res; 343 error: 344 goto exit; 345 } 346 347 static void 348 shape_release(struct shape* shape) 349 { 350 ASSERT(shape); 351 schiff_mesh_release(&shape->mesh); 352 if(shape->shape) S3D(shape_ref_put(shape->shape)); 353 memset(shape, 0, sizeof(*shape)); 354 } 355 356 static res_T 357 shape_cylinder_init 358 (struct shape* shape, 359 struct s3d_device* s3d, 360 const struct schiff_geometry* geometry) 361 { 362 struct mesh_context ctx; 363 struct s3d_vertex_data attrib; 364 const struct schiff_cylinder* cylinder; 365 double radius, height; 366 size_t nverts, nprims; 367 res_T res = RES_OK; 368 ASSERT(shape && geometry); 369 ASSERT(geometry->type == SCHIFF_CYLINDER 370 || geometry->type == SCHIFF_CYLINDER_AS_SPHERE); 371 372 memset(shape, 0, sizeof(*shape)); 373 cylinder = &geometry->data.cylinder; 374 375 376 res = schiff_mesh_init_cylinder 377 (&mem_default_allocator, &shape->mesh, cylinder->nslices); 378 if(res != RES_OK) goto error; 379 380 shape->geometry = geometry; 381 382 if(geometry->type == SCHIFF_CYLINDER) /* That's all folks */ 383 goto exit; 384 385 /* Create the Star-3D cylinder shape */ 386 res = s3d_shape_create_mesh(s3d, &shape->shape); 387 if(res != RES_OK) goto error; 388 389 /* Cylinder parameters must be constant */ 390 ASSERT(cylinder->radius.distribution == SCHIFF_PARAM_CONSTANT); 391 ASSERT(cylinder->height.distribution == SCHIFF_PARAM_CONSTANT); 392 radius = cylinder->radius.data.constant; 393 height = cylinder->height.data.constant; 394 395 /* Setup the context of the cylinder */ 396 ctx.type = SCHIFF_CYLINDER; 397 ctx.mesh = &shape->mesh; 398 ctx.data.cylinder.radius = (float)radius; 399 ctx.data.cylinder.height = (float)height; 400 401 /* Define the position vertex attribute */ 402 attrib.usage = S3D_POSITION; 403 attrib.type = S3D_FLOAT3; 404 attrib.get = cylinder_get_position; 405 406 nverts = darray_float_size_get(&shape->mesh.vertices.cartesian) / 3/*#coords*/; 407 nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices per prim*/; 408 409 res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, 410 geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); 411 if(res != RES_OK) goto error; 412 413 shape->volume = PI*height*radius*radius; /* Analytic volume */ 414 415 schiff_mesh_release(&shape->mesh); 416 417 exit: 418 return res; 419 error: 420 shape_release(shape); 421 goto exit; 422 } 423 424 static res_T 425 shape_cylinder_generate_s3d_shape 426 (const struct shape* shape, 427 struct ssp_rng* rng, 428 struct s3d_shape* s3d_shape) 429 { 430 struct s3d_vertex_data attrib; 431 struct mesh_context ctx; 432 size_t nverts, nprims; 433 ASSERT(shape && shape->geometry->type == SCHIFF_CYLINDER); 434 435 ctx.type = SCHIFF_CYLINDER; 436 ctx.mesh = &shape->mesh; 437 ctx.data.cylinder.radius = (float)eval_param 438 (&shape->geometry->data.cylinder.radius, rng); 439 ctx.data.cylinder.height = (float)eval_param 440 (&shape->geometry->data.cylinder.height, rng); 441 442 attrib.usage = S3D_POSITION; 443 attrib.type = S3D_FLOAT3; 444 attrib.get = cylinder_get_position; 445 446 nverts = darray_float_size_get(&shape->mesh.vertices.cartesian) / 3/*#coords*/; 447 nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices per prim*/; 448 449 return s3d_mesh_setup_indexed_vertices(s3d_shape, (unsigned)nprims, 450 geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); 451 } 452 453 static res_T 454 shape_cylinder_as_sphere_generate_s3d_shape 455 (const struct shape* shape, 456 struct ssp_rng* rng, 457 struct s3d_shape* s3d_shape) 458 { 459 ASSERT(shape && shape->geometry->type == SCHIFF_CYLINDER_AS_SPHERE); 460 (void)rng; 461 return s3d_mesh_copy(shape->shape, s3d_shape); 462 } 463 464 static res_T 465 shape_cylinder_as_sphere_sample_volume_scaling 466 (const struct shape* shape, 467 struct ssp_rng* rng, 468 double* volume_scaling) 469 { 470 double radius, sphere_volume; 471 ASSERT(shape && volume_scaling); 472 ASSERT(shape->geometry->type == SCHIFF_CYLINDER_AS_SPHERE); 473 474 radius = eval_param(&shape->geometry->data.cylinder.radius_sphere, rng); 475 sphere_volume = 4.0/3.0 * PI * radius*radius*radius; 476 *volume_scaling = sphere_volume / shape->volume; 477 return RES_OK; 478 } 479 480 static res_T 481 shape_ellipsoid_init 482 (struct shape* shape, 483 struct s3d_device* s3d, 484 const struct schiff_geometry* geometry) 485 { 486 const struct schiff_ellipsoid* ellipsoid; 487 struct mesh_context ctx; 488 struct s3d_vertex_data attrib; 489 double a, b, c, a2, b2, c2; 490 size_t nverts, nprims; 491 res_T res = RES_OK; 492 ASSERT(shape && geometry); 493 ASSERT(geometry->type == SCHIFF_ELLIPSOID 494 || geometry->type == SCHIFF_ELLIPSOID_AS_SPHERE); 495 496 ellipsoid = &geometry->data.ellipsoid; 497 memset(shape, 0, sizeof(*shape)); 498 499 res = schiff_mesh_init_sphere_polar 500 (&mem_default_allocator, &shape->mesh, ellipsoid->nslices); 501 if(res != RES_OK) goto error; 502 503 shape->geometry = geometry; 504 505 if(geometry->type == SCHIFF_ELLIPSOID) /* That's all folks */ 506 goto exit; 507 508 res = s3d_shape_create_mesh(s3d, &shape->shape); 509 if(res != RES_OK) goto error; 510 511 /* Expecting constant parameters */ 512 ASSERT(ellipsoid->a.distribution == SCHIFF_PARAM_CONSTANT); 513 ASSERT(ellipsoid->c.distribution == SCHIFF_PARAM_CONSTANT); 514 a = ellipsoid->a.data.constant; 515 c = ellipsoid->c.data.constant; 516 b = a; 517 a2 = a*a; 518 b2 = b*b; 519 c2 = c*c; 520 521 ctx.mesh = &shape->mesh; 522 ctx.type = SCHIFF_ELLIPSOID; 523 ctx.data.ellipsoid.abc = a*b*c; 524 ctx.data.ellipsoid.b2c2 = b2*c2; 525 ctx.data.ellipsoid.a2c2 = a2*c2; 526 ctx.data.ellipsoid.a2b2 = a2*b2; 527 528 attrib.usage = S3D_POSITION; 529 attrib.type = S3D_FLOAT3; 530 attrib.get = ellipsoid_get_position; 531 532 nverts = darray_uint_size_get(&shape->mesh.vertices.polar) / 2/*#theta/phi*/; 533 nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; 534 535 res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, 536 geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); 537 if(res != RES_OK) goto error; 538 539 shape->volume = 4.0/3.0 * PI * ctx.data.ellipsoid.abc; 540 541 schiff_mesh_release(&shape->mesh); 542 543 exit: 544 return res; 545 error: 546 shape_release(shape); 547 goto exit; 548 } 549 550 static res_T 551 shape_ellipsoid_generate_s3d_shape 552 (const struct shape* shape, 553 struct ssp_rng* rng, 554 struct s3d_shape* s3d_shape) 555 { 556 struct s3d_vertex_data attrib; 557 struct mesh_context ctx; 558 size_t nverts, nprims; 559 double a, b, c, a2, b2, c2; 560 ASSERT(shape && shape->geometry->type == SCHIFF_ELLIPSOID); 561 562 ctx.type = SCHIFF_ELLIPSOID; 563 ctx.mesh = &shape->mesh; 564 a = eval_param(&shape->geometry->data.ellipsoid.a, rng); 565 c = eval_param(&shape->geometry->data.ellipsoid.c, rng); 566 567 b = a; 568 a2 = a*a; 569 b2 = b*b; 570 c2 = c*c; 571 572 ctx.mesh = &shape->mesh; 573 ctx.type = SCHIFF_ELLIPSOID; 574 ctx.data.ellipsoid.abc = a*b*c; 575 ctx.data.ellipsoid.b2c2 = b2*c2; 576 ctx.data.ellipsoid.a2c2 = a2*c2; 577 ctx.data.ellipsoid.a2b2 = a2*b2; 578 579 attrib.usage = S3D_POSITION; 580 attrib.type = S3D_FLOAT3; 581 attrib.get = ellipsoid_get_position; 582 583 nverts = darray_uint_size_get(&shape->mesh.vertices.polar) / 2/*#theta/phi*/; 584 nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; 585 586 return s3d_mesh_setup_indexed_vertices(s3d_shape, (unsigned)nprims, 587 geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); 588 } 589 590 static res_T 591 shape_ellipsoid_as_sphere_generate_s3d_shape 592 (const struct shape* shape, 593 struct ssp_rng* rng, 594 struct s3d_shape* s3d_shape) 595 { 596 ASSERT(shape && shape->geometry->type == SCHIFF_ELLIPSOID_AS_SPHERE); 597 (void)rng; 598 return s3d_mesh_copy(shape->shape, s3d_shape); 599 } 600 601 static res_T 602 shape_ellipsoid_as_sphere_sample_volume_scaling 603 (const struct shape* shape, 604 struct ssp_rng* rng, 605 double* volume_scaling) 606 { 607 double radius, sphere_volume; 608 ASSERT(shape && volume_scaling); 609 ASSERT(shape->geometry->type == SCHIFF_ELLIPSOID_AS_SPHERE); 610 611 radius = eval_param(&shape->geometry->data.ellipsoid.radius_sphere, rng); 612 sphere_volume = 4.0/3.0 * PI * radius*radius*radius; 613 *volume_scaling = sphere_volume / shape->volume; 614 return RES_OK; 615 } 616 617 static res_T 618 shape_helical_pipe_init 619 (struct shape* shape, 620 struct s3d_device* s3d, 621 const struct schiff_geometry* geometry) 622 { 623 const struct schiff_helical_pipe* helical_pipe; 624 float* vertices = NULL; 625 struct helical_pipe_mesh_context ctx; 626 struct s3d_vertex_data attrib; 627 size_t nprims, nverts; 628 double pitch, height, hradius, cradius; 629 res_T res = RES_OK; 630 631 ASSERT(geometry); 632 ASSERT(geometry->type == SCHIFF_HELICAL_PIPE 633 || geometry->type == SCHIFF_HELICAL_PIPE_AS_SPHERE); 634 635 helical_pipe = &geometry->data.helical_pipe; 636 memset(shape, 0, sizeof(*shape)); 637 638 res = schiff_mesh_init_helical_pipe(&mem_default_allocator, &shape->mesh, 639 helical_pipe->nslices_helicoid, helical_pipe->nslices_circle); 640 if(res != RES_OK) goto error; 641 642 shape->geometry = geometry; 643 644 if(geometry->type == SCHIFF_HELICAL_PIPE) /* That's all folks */ 645 goto exit; 646 647 res = s3d_shape_create_mesh(s3d, &shape->shape); 648 if(res != RES_OK) goto error; 649 650 /* Expecting constant parameters */ 651 ASSERT(helical_pipe->radius_helicoid.distribution == SCHIFF_PARAM_CONSTANT); 652 ASSERT(helical_pipe->radius_circle.distribution == SCHIFF_PARAM_CONSTANT); 653 ASSERT(helical_pipe->pitch.distribution == SCHIFF_PARAM_CONSTANT); 654 ASSERT(helical_pipe->height.distribution == SCHIFF_PARAM_CONSTANT); 655 hradius = helical_pipe->radius_helicoid.data.constant; 656 cradius = helical_pipe->radius_circle.data.constant; 657 pitch = helical_pipe->pitch.data.constant; 658 height = helical_pipe->height.data.constant; 659 660 /* Setup the mesh vertices */ 661 res = schiff_mesh_helical_pipe_create_vertices(&shape->mesh, pitch, height, 662 hradius, cradius, helical_pipe->nslices_helicoid, 663 helical_pipe->nslices_circle, &nverts, &vertices); 664 if(res != RES_OK) return res; 665 666 ctx.vertices = vertices; 667 ctx.indices = darray_uint_cdata_get(&shape->mesh.indices); 668 669 attrib.usage = S3D_POSITION; 670 attrib.type = S3D_FLOAT3; 671 attrib.get = helical_pipe_get_position; 672 673 nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices per prim*/; 674 675 res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, 676 helical_pipe_get_indices, (unsigned)nverts, &attrib, 1, &ctx); 677 if(res != RES_OK) goto error; 678 679 res = compute_s3d_shape_volume(s3d, shape->shape, &shape->volume); 680 if(res != RES_OK) goto error; 681 682 schiff_mesh_release(&shape->mesh); 683 684 exit: 685 if(vertices) { 686 schiff_mesh_helical_pipe_destroy_vertices(&shape->mesh, vertices); 687 } 688 return res; 689 error: 690 shape_release(shape); 691 goto exit; 692 } 693 694 static res_T 695 shape_helical_pipe_generate_s3d_shape 696 (struct shape* shape, 697 struct ssp_rng* rng, 698 struct s3d_shape* s3d_shape) 699 { 700 const struct schiff_helical_pipe* helical_pipe; 701 float* vertices = NULL; 702 struct helical_pipe_mesh_context ctx; 703 struct s3d_vertex_data attrib; 704 size_t nprims, nverts; 705 double pitch, height, hradius, cradius; 706 res_T res = RES_OK; 707 ASSERT(shape && shape->geometry->type == SCHIFF_HELICAL_PIPE); 708 709 helical_pipe = &shape->geometry->data.helical_pipe; 710 711 /* Sample the helical pipe attributes */ 712 hradius = eval_param(&helical_pipe->radius_helicoid, rng); 713 cradius = eval_param(&helical_pipe->radius_circle, rng); 714 pitch = eval_param(&helical_pipe->pitch, rng); 715 height = eval_param(&helical_pipe->height, rng); 716 717 /* Setup the mesh vertices */ 718 res = schiff_mesh_helical_pipe_create_vertices(&shape->mesh, pitch, height, 719 hradius, cradius, helical_pipe->nslices_helicoid, 720 helical_pipe->nslices_circle, &nverts, &vertices); 721 if(res != RES_OK) goto error; 722 723 ctx.vertices = vertices; 724 ctx.indices = darray_uint_cdata_get(&shape->mesh.indices); 725 726 attrib.usage = S3D_POSITION; 727 attrib.type = S3D_FLOAT3; 728 attrib.get = helical_pipe_get_position; 729 730 nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices per prim*/; 731 732 res = s3d_mesh_setup_indexed_vertices(s3d_shape, (unsigned)nprims, 733 helical_pipe_get_indices, (unsigned)nverts, &attrib, 1, &ctx); 734 if(res != RES_OK) goto error; 735 736 exit: 737 if(vertices) { 738 schiff_mesh_helical_pipe_destroy_vertices(&shape->mesh, vertices); 739 } 740 return res; 741 error: 742 goto exit; 743 } 744 745 static res_T 746 shape_helical_pipe_as_sphere_generate_s3d_shape 747 (struct shape* shape, 748 struct ssp_rng* rng, 749 struct s3d_shape* s3d_shape) 750 { 751 ASSERT(shape && shape->geometry->type == SCHIFF_HELICAL_PIPE_AS_SPHERE); 752 (void)rng; 753 return s3d_mesh_copy(shape->shape, s3d_shape); 754 } 755 756 static res_T 757 shape_helical_pipe_as_sphere_sample_volume_scaling 758 (const struct shape* shape, 759 struct ssp_rng* rng, 760 double* volume_scaling) 761 { 762 double radius, sphere_volume; 763 ASSERT(shape && volume_scaling); 764 ASSERT(shape->geometry->type == SCHIFF_HELICAL_PIPE_AS_SPHERE); 765 766 radius = eval_param(&shape->geometry->data.helical_pipe.radius_sphere, rng); 767 sphere_volume = 4.0/3.0 * PI * radius*radius*radius; 768 *volume_scaling = sphere_volume / shape->volume; 769 return RES_OK; 770 } 771 772 static res_T 773 shape_sphere_init 774 (struct shape* shape, 775 struct s3d_device* s3d, 776 const struct schiff_geometry* geometry) 777 { 778 struct mesh_context ctx; 779 struct s3d_vertex_data attrib; 780 const struct schiff_sphere* sphere; 781 size_t nverts, nprims; 782 res_T res = RES_OK; 783 ASSERT(shape && geometry && geometry->type == SCHIFF_SPHERE); 784 785 sphere = &geometry->data.sphere; 786 memset(shape, 0, sizeof(*shape)); 787 788 /* Generate the sphere mesh */ 789 res = schiff_mesh_init_sphere 790 (&mem_default_allocator, &shape->mesh, sphere->nslices); 791 if(res != RES_OK) goto error; 792 793 shape->geometry = geometry; 794 795 /* Create the Star-3D sphere shape */ 796 res = s3d_shape_create_mesh(s3d, &shape->shape); 797 if(res != RES_OK) goto error; 798 799 /* Setup the context of the sphere*/ 800 ctx.type = SCHIFF_SPHERE; 801 ctx.mesh = &shape->mesh; 802 ctx.data.sphere.radius = 1.f; 803 804 /* Define the position vertex attribute */ 805 attrib.usage = S3D_POSITION; 806 attrib.type = S3D_FLOAT3; 807 attrib.get = sphere_get_position; 808 809 nverts = darray_float_size_get(&shape->mesh.vertices.cartesian) / 3/*#coords*/; 810 nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; 811 812 res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, 813 geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); 814 if(res != RES_OK) goto error; 815 816 shape->volume = 4.0/3.0*PI; /* Analytic volume. The radius is implicitly 1 */ 817 818 exit: 819 schiff_mesh_release(&shape->mesh); 820 return res; 821 error: 822 shape_release(shape); 823 goto exit; 824 } 825 826 static res_T 827 shape_sphere_generate_s3d_shape 828 (const struct shape* shape, 829 struct ssp_rng* rng, 830 struct s3d_shape* s3d_shape) 831 { 832 ASSERT(shape && shape->geometry->type == SCHIFF_SPHERE); 833 (void)rng; 834 return s3d_mesh_copy(shape->shape, s3d_shape); 835 } 836 837 static res_T 838 shape_sphere_sample_volume_scaling 839 (const struct shape* shape, 840 struct ssp_rng* rng, 841 double* volume_scaling) 842 { 843 double radius, sphere_volume; 844 ASSERT(shape && volume_scaling); 845 ASSERT(shape->geometry->type == SCHIFF_SPHERE); 846 847 radius = eval_param(&shape->geometry->data.sphere.radius, rng); 848 sphere_volume = 4.0/3.0 * PI * radius*radius*radius; 849 *volume_scaling = sphere_volume / shape->volume; 850 return RES_OK; 851 } 852 853 static res_T 854 shape_supershape_init 855 (struct shape* shape, 856 struct s3d_device* s3d, 857 const struct schiff_geometry* geometry) 858 { 859 const struct schiff_supershape* sshape; 860 struct mesh_context ctx; 861 struct s3d_vertex_data attrib; 862 size_t nverts, nprims; 863 int iform, iattr; 864 res_T res = RES_OK; 865 ASSERT(geometry); 866 ASSERT(geometry->type == SCHIFF_SUPERSHAPE 867 || geometry->type == SCHIFF_SUPERSHAPE_AS_SPHERE); 868 869 sshape = &geometry->data.supershape; 870 memset(shape, 0, sizeof(*shape)); 871 872 res = schiff_mesh_init_sphere_polar(&mem_default_allocator, 873 &shape->mesh, geometry->data.supershape.nslices); 874 if(res != RES_OK) goto error; 875 876 shape->geometry = geometry; 877 878 if(geometry->type == SCHIFF_SUPERSHAPE) /* That's all folks */ 879 goto exit; 880 881 res = s3d_shape_create_mesh(s3d, &shape->shape); 882 if(res != RES_OK) goto error; 883 884 ctx.type = SCHIFF_SUPERSHAPE; 885 ctx.mesh = &shape->mesh; 886 FOR_EACH(iform, 0, 2) { 887 FOR_EACH(iattr, 0, 6) { 888 /* Expecting constant parameters */ 889 ASSERT(sshape->formulas[iform][iattr].distribution == SCHIFF_PARAM_CONSTANT); 890 ctx.data.supershape.formulas[iform][iattr] = 891 sshape->formulas[iform][iattr].data.constant; 892 }} 893 894 attrib.usage = S3D_POSITION; 895 attrib.type = S3D_FLOAT3; 896 attrib.get = supershape_get_position; 897 898 nverts = darray_uint_size_get(&shape->mesh.vertices.polar) / 2/*#theta/phi*/; 899 nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; 900 901 res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, 902 geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); 903 if(res != RES_OK) goto error; 904 905 res = compute_s3d_shape_volume(s3d, shape->shape, &shape->volume); 906 if(res != RES_OK) goto error; 907 908 schiff_mesh_release(&shape->mesh); 909 910 exit: 911 return res; 912 error: 913 shape_release(shape); 914 goto exit; 915 } 916 917 static res_T 918 shape_supershape_generate_s3d_shape 919 (const struct shape* shape, 920 struct ssp_rng* rng, 921 struct s3d_shape* s3d_shape) 922 { 923 const struct schiff_supershape* sshape; 924 struct mesh_context ctx; 925 struct s3d_vertex_data attrib; 926 size_t nverts, nprims; 927 int iform, iattr; 928 ASSERT(shape && shape->geometry->type == SCHIFF_SUPERSHAPE); 929 930 sshape = &shape->geometry->data.supershape; 931 ctx.mesh = &shape->mesh; 932 ctx.type = SCHIFF_SUPERSHAPE; 933 934 FOR_EACH(iform, 0, 2) { 935 FOR_EACH(iattr, 0, 6) { 936 ctx.data.supershape.formulas[iform][iattr] = 937 (float)eval_param(&sshape->formulas[iform][iattr], rng); 938 }} 939 940 attrib.usage = S3D_POSITION; 941 attrib.type = S3D_FLOAT3; 942 attrib.get = supershape_get_position; 943 944 nverts = darray_uint_size_get(&shape->mesh.vertices.polar) / 2/*#theta/phi*/; 945 nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; 946 947 return s3d_mesh_setup_indexed_vertices(s3d_shape, (unsigned)nprims, 948 geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); 949 } 950 951 static res_T 952 shape_supershape_as_sphere_generate_s3d_shape 953 (const struct shape* shape, 954 struct ssp_rng* rng, 955 struct s3d_shape* s3d_shape) 956 { 957 ASSERT(shape && shape->geometry->type == SCHIFF_SUPERSHAPE_AS_SPHERE); 958 (void)rng; 959 return s3d_mesh_copy(shape->shape, s3d_shape); 960 } 961 962 static res_T 963 shape_supershape_as_sphere_sample_volume_scaling 964 (const struct shape* shape, 965 struct ssp_rng* rng, 966 double* volume_scaling) 967 { 968 double radius, sphere_volume; 969 ASSERT(shape && volume_scaling); 970 ASSERT(shape->geometry->type == SCHIFF_SUPERSHAPE_AS_SPHERE); 971 972 radius = eval_param(&shape->geometry->data.supershape.radius_sphere, rng); 973 sphere_volume = 4.0/3.0 * PI * radius*radius*radius; 974 *volume_scaling = sphere_volume / shape->volume; 975 return RES_OK; 976 } 977 978 static res_T 979 geometry_sample_shape 980 (struct ssp_rng* rng, 981 void** shape_data, 982 struct s3d_shape* s3d_shape, 983 void* ctx) 984 { 985 struct geometry_distribution_context* distrib = ctx; 986 struct shape* shape; 987 size_t isamp; 988 res_T res = RES_OK; 989 ASSERT(rng && shape_data && ctx); 990 991 isamp = ssp_ranst_discrete_get(rng, distrib->ran_geometries); 992 shape = distrib->shapes + isamp; 993 994 switch(shape->geometry->type) { 995 case SCHIFF_ELLIPSOID: 996 res = shape_ellipsoid_generate_s3d_shape(shape, rng, s3d_shape); 997 break; 998 case SCHIFF_ELLIPSOID_AS_SPHERE: 999 res = shape_ellipsoid_as_sphere_generate_s3d_shape(shape, rng, s3d_shape); 1000 break; 1001 case SCHIFF_CYLINDER: 1002 res = shape_cylinder_generate_s3d_shape(shape, rng, s3d_shape); 1003 break; 1004 case SCHIFF_CYLINDER_AS_SPHERE: 1005 res = shape_cylinder_as_sphere_generate_s3d_shape(shape, rng, s3d_shape); 1006 break; 1007 case SCHIFF_HELICAL_PIPE: 1008 res = shape_helical_pipe_generate_s3d_shape(shape, rng, s3d_shape); 1009 break; 1010 case SCHIFF_HELICAL_PIPE_AS_SPHERE: 1011 res = shape_helical_pipe_as_sphere_generate_s3d_shape(shape, rng, s3d_shape); 1012 break; 1013 case SCHIFF_SPHERE: 1014 res = shape_sphere_generate_s3d_shape(shape, rng, s3d_shape); 1015 break; 1016 case SCHIFF_SUPERSHAPE: 1017 res = shape_supershape_generate_s3d_shape(shape, rng, s3d_shape); 1018 break; 1019 case SCHIFF_SUPERSHAPE_AS_SPHERE: 1020 res = shape_supershape_as_sphere_generate_s3d_shape(shape, rng, s3d_shape); 1021 break; 1022 default: FATAL("Unreachable code.\n"); break; 1023 } 1024 if(res != RES_OK) goto error; 1025 1026 exit: 1027 *shape_data = shape; 1028 return res; 1029 error: 1030 goto exit; 1031 } 1032 1033 static res_T 1034 geometry_sample_volume_scaling 1035 (struct ssp_rng* rng, 1036 void* shape_data, 1037 double* volume_scaling, 1038 void* ctx) 1039 { 1040 struct shape* shape = shape_data; 1041 res_T res = RES_OK; 1042 ASSERT(rng && shape_data && volume_scaling && ctx); 1043 (void)ctx; 1044 1045 switch(shape->geometry->type) { 1046 case SCHIFF_ELLIPSOID: 1047 case SCHIFF_CYLINDER: 1048 case SCHIFF_HELICAL_PIPE: 1049 case SCHIFF_SUPERSHAPE: 1050 *volume_scaling = 1; 1051 break; 1052 case SCHIFF_ELLIPSOID_AS_SPHERE: 1053 res = shape_ellipsoid_as_sphere_sample_volume_scaling 1054 (shape, rng, volume_scaling); 1055 break; 1056 case SCHIFF_CYLINDER_AS_SPHERE: 1057 res = shape_cylinder_as_sphere_sample_volume_scaling 1058 (shape, rng, volume_scaling); 1059 break; 1060 case SCHIFF_HELICAL_PIPE_AS_SPHERE: 1061 res = shape_helical_pipe_as_sphere_sample_volume_scaling 1062 (shape, rng, volume_scaling); 1063 break; 1064 case SCHIFF_SPHERE: 1065 res = shape_sphere_sample_volume_scaling 1066 (shape, rng, volume_scaling); 1067 break; 1068 case SCHIFF_SUPERSHAPE_AS_SPHERE: 1069 res = shape_supershape_as_sphere_sample_volume_scaling 1070 (shape, rng, volume_scaling); 1071 break; 1072 default: FATAL("Unreachable code.\n"); break; 1073 } 1074 if(res != RES_OK) goto error; 1075 ASSERT(*volume_scaling > 0); 1076 1077 exit: 1078 return res; 1079 error: 1080 goto exit; 1081 } 1082 1083 /******************************************************************************* 1084 * Local functions 1085 ******************************************************************************/ 1086 void 1087 schiff_geometry_distribution_release 1088 (struct sschiff_geometry_distribution* distrib) 1089 { 1090 struct geometry_distribution_context* ctx; 1091 size_t i, count; 1092 ASSERT(distrib); 1093 if(!distrib->context) return; 1094 1095 ctx = distrib->context; 1096 if(ctx->shapes) { 1097 count = sa_size(ctx->shapes); 1098 FOR_EACH(i, 0, count) 1099 shape_release(&ctx->shapes[i]); 1100 sa_release(ctx->shapes); 1101 } 1102 if(ctx->ran_geometries) 1103 SSP(ranst_discrete_ref_put(ctx->ran_geometries)); 1104 mem_rm(ctx); 1105 } 1106 1107 res_T 1108 schiff_geometry_distribution_init 1109 (struct sschiff_geometry_distribution* distrib, 1110 struct s3d_device* s3d, 1111 const struct schiff_geometry* geoms, 1112 const size_t ngeoms, 1113 const double characteristic_length, 1114 struct ssp_ranst_discrete* ran_geoms, 1115 struct schiff_optical_properties* properties) 1116 { 1117 struct geometry_distribution_context* ctx = NULL; 1118 size_t i; 1119 res_T res = RES_OK; 1120 ASSERT(distrib && geoms && ngeoms && ran_geoms); 1121 /* Note that properties may be NULL if the "-G" option is defined */ 1122 1123 *distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION; 1124 1125 ctx = mem_calloc(1, sizeof(struct geometry_distribution_context)); 1126 if(!ctx) { 1127 fprintf(stderr, 1128 "Couldn't allocate the geometry distribution context\n"); 1129 res = RES_MEM_ERR; 1130 goto error; 1131 } 1132 /* Directly setup the distribution context to avoid memory leaks on error */ 1133 distrib->context = ctx; 1134 1135 if(!sa_add(ctx->shapes, ngeoms)) { 1136 fprintf(stderr, 1137 "Couldn't allocate the shapes of the geometry distributions.\n"); 1138 res = RES_MEM_ERR; 1139 goto error; 1140 } 1141 memset(ctx->shapes, 0, sizeof(ctx->shapes[0])*ngeoms); 1142 1143 FOR_EACH(i, 0, ngeoms) { 1144 /* Generate the mesh template and setup its distribution context */ 1145 switch(geoms[i].type) { 1146 case SCHIFF_ELLIPSOID: 1147 case SCHIFF_ELLIPSOID_AS_SPHERE: 1148 res = shape_ellipsoid_init(&ctx->shapes[i], s3d, &geoms[i]); 1149 break; 1150 case SCHIFF_CYLINDER: 1151 case SCHIFF_CYLINDER_AS_SPHERE: 1152 res = shape_cylinder_init(&ctx->shapes[i], s3d, &geoms[i]); 1153 break; 1154 case SCHIFF_HELICAL_PIPE: 1155 case SCHIFF_HELICAL_PIPE_AS_SPHERE: 1156 res = shape_helical_pipe_init(&ctx->shapes[i], s3d, &geoms[i]); 1157 break; 1158 case SCHIFF_SPHERE: 1159 res = shape_sphere_init(&ctx->shapes[i], s3d, &geoms[i]); 1160 break; 1161 case SCHIFF_SUPERSHAPE: 1162 case SCHIFF_SUPERSHAPE_AS_SPHERE: 1163 res = shape_supershape_init(&ctx->shapes[i], s3d, &geoms[i]); 1164 break; 1165 default: FATAL("Unreachable code.\n"); break; 1166 } 1167 if(res != RES_OK) { 1168 fprintf(stderr, 1169 "Couldn't initialise the shape of the geometry distribution.\n"); 1170 goto error; 1171 } 1172 } 1173 ctx->properties = properties; 1174 SSP(ranst_discrete_ref_get(ran_geoms)); 1175 ctx->ran_geometries = ran_geoms; 1176 1177 distrib->material.get_property = get_material_property; 1178 distrib->material.material = properties; 1179 distrib->sample = geometry_sample_shape; 1180 distrib->sample_volume_scaling = geometry_sample_volume_scaling; 1181 distrib->characteristic_length = characteristic_length; 1182 1183 exit: 1184 return res; 1185 error: 1186 schiff_geometry_distribution_release(distrib); 1187 goto exit; 1188 } 1189