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