schiff

Estimate the radiative properties of soft particless
git clone git://git.meso-star.com/schiff.git
Log | Files | Refs | README | LICENSE

commit 684c105d47317bce6b24d432af0c5a8c881fb1e3
parent 6d287ab6e6b945496cfe75bc563118824d1ab106
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue,  5 Apr 2016 16:31:08 +0200

Handle the Star-Schiff volume scaling attribute

The equivalent sphere distributions rely on the volume scaling attribute
of the geometry distribution. Change the YAML file format to declare an
equivalent sphere distribution of the cylinder, and the ellipsoid.  Add
the equivalent sphere distribution for the helical pipe and the
supershape.

Diffstat:
Msrc/schiff.c | 32++++++++++++++++++++------------
Msrc/schiff_args.c | 186++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
Msrc/schiff_geometry.c | 747+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------
Msrc/schiff_geometry.h | 50+++++++++++++++++++++++++++++++++-----------------
Msrc/schiff_mesh.c | 1+
5 files changed, 741 insertions(+), 275 deletions(-)

diff --git a/src/schiff.c b/src/schiff.c @@ -29,22 +29,17 @@ static res_T dump_geometries (const struct schiff_args* args, + struct s3d_device* s3d, struct sschiff_geometry_distribution* distrib, struct ssp_rng* rng, FILE* stream) { - struct s3d_device* s3d = NULL; struct s3d_scene* scn = NULL; struct s3d_shape* shape = NULL; unsigned i; res_T res = RES_OK; ASSERT(args && distrib && rng && stream); - res = s3d_device_create(NULL, NULL, 0, &s3d); - if(res != RES_OK) { - fprintf(stderr, "Couldn't create the Star-3D device.\n"); - goto error; - } res = s3d_scene_create(s3d, &scn); if(res != RES_OK) { fprintf(stderr, "Couldn't create the Star-3D scene.\n"); @@ -64,12 +59,16 @@ dump_geometries FOR_EACH(i, 0, args->ngeoms_dump) { unsigned ivert, nverts; unsigned itri, ntris; - res = distrib->sample(rng, shape, distrib->context); + double volume_scaling; + double dist_scaling; + res = distrib->sample(rng, shape, &volume_scaling, distrib->context); if(res != RES_OK) { fprintf(stderr, "Couldn't sample the micro organism geometry.\n"); goto error; } + dist_scaling = pow(volume_scaling, 1.0/3.0); + fprintf(stream, "g schiff_geometry_%u\n", i); /* Dump vertex position */ @@ -78,6 +77,7 @@ dump_geometries FOR_EACH(ivert, 0, nverts) { struct s3d_attrib attr; S3D(mesh_get_vertex_attrib(shape, ivert, S3D_POSITION, &attr)); + f3_mulf(attr.value, attr.value, (float)dist_scaling); fprintf(stream, "v %f %f %f\n", SPLIT3(attr.value)); } @@ -92,7 +92,6 @@ dump_geometries } exit: - if(s3d) S3D(device_ref_put(s3d)); if(scn) S3D(scene_ref_put(scn)); if(shape) S3D(shape_ref_put(shape)); return res; @@ -310,6 +309,7 @@ write_phase_functions_invcum static res_T run_integration (const struct schiff_args* args, + struct s3d_device* s3d, struct sschiff_geometry_distribution* distrib, struct ssp_rng* rng, FILE* stream) @@ -320,7 +320,7 @@ run_integration res_T res = RES_OK; ASSERT(args && sa_size(args->wavelengths) && rng && stream); - res = sschiff_device_create(NULL, NULL, args->nthreads, 1, NULL, &sschiff); + res = sschiff_device_create(NULL, NULL, args->nthreads, 1, s3d, &sschiff); if(res != RES_OK) { fprintf(stderr, "Couldn't create the Star Schiff device.\n"); goto error; @@ -359,6 +359,7 @@ run(const struct schiff_args* args) FILE* fp = stdout; struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION; struct ssp_rng* rng = NULL; + struct s3d_device* s3d = NULL; res_T res = RES_OK; ASSERT(args); @@ -378,7 +379,13 @@ run(const struct schiff_args* args) goto error; } - res = schiff_geometry_distribution_init(&distrib, args->geoms, + res = s3d_device_create(NULL, &mem_default_allocator, 0, &s3d); + if(res != RES_OK) { + fprintf(stderr, "Couldn't create the Star-3D device.\n"); + goto error; + } + + res = schiff_geometry_distribution_init(&distrib, s3d, args->geoms, sa_size(args->geoms), args->characteristic_length, args->ran_geoms, args->properties); if(res != RES_OK) { @@ -388,9 +395,9 @@ run(const struct schiff_args* args) } if(args->ngeoms_dump) { - res = dump_geometries(args, &distrib, rng, fp); + res = dump_geometries(args, s3d, &distrib, rng, fp); } else { - res = run_integration(args, &distrib, rng, fp); + res = run_integration(args, s3d, &distrib, rng, fp); } /* Release before the check of the integration error code since if an error @@ -403,6 +410,7 @@ run(const struct schiff_args* args) exit: if(fp && fp != stdout) fclose(fp); if(rng) SSP(rng_ref_put(rng)); + if(s3d) S3D(device_ref_put(s3d)); return res; error: goto exit; diff --git a/src/schiff_args.c b/src/schiff_args.c @@ -116,32 +116,35 @@ geometry_release(struct schiff_geometry* geom) ASSERT(geom); switch(geom->type) { case SCHIFF_ELLIPSOID: + case SCHIFF_ELLIPSOID_AS_SPHERE: param_release(&geom->data.ellipsoid.a); param_release(&geom->data.ellipsoid.c); - case SCHIFF_ELLIPSOID_AS_SPHERE: - param_release(&geom->data.ellipsoid.radius); + param_release(&geom->data.ellipsoid.radius_sphere); break; + case SCHIFF_CYLINDER_AS_SPHERE: case SCHIFF_CYLINDER: param_release(&geom->data.cylinder.radius); param_release(&geom->data.cylinder.height); - break; - case SCHIFF_CYLINDER_AS_SPHERE: - param_release(&geom->data.cylinder.radius); + param_release(&geom->data.cylinder.radius_sphere); break; case SCHIFF_HELICAL_PIPE: + case SCHIFF_HELICAL_PIPE_AS_SPHERE: param_release(&geom->data.helical_pipe.pitch); param_release(&geom->data.helical_pipe.height); - param_release(&geom->data.helical_pipe.helicoid_radius); - param_release(&geom->data.helical_pipe.circle_radius); + param_release(&geom->data.helical_pipe.radius_helicoid); + param_release(&geom->data.helical_pipe.radius_circle); + param_release(&geom->data.helical_pipe.radius_sphere); break; case SCHIFF_SPHERE: param_release(&geom->data.sphere.radius); break; case SCHIFF_SUPERSHAPE: + case SCHIFF_SUPERSHAPE_AS_SPHERE: FOR_EACH(i, 0, 6) { param_release(&geom->data.supershape.formulas[0][i]); param_release(&geom->data.supershape.formulas[1][i]); } + param_release(&geom->data.supershape.radius_sphere); break; case SCHIFF_NONE: /* Do nothing */ break; default: FATAL("Unreachable code\n"); break; @@ -639,9 +642,8 @@ parse_yaml_ellipsoid PROBA = BIT(0), LENGTH_a = BIT(1), LENGTH_c = BIT(2), - RADIUS = BIT(3), - ASPECT_RATIO = BIT(4), - SLICES = BIT(5) + RADIUS_SPHERE = BIT(3), + SLICES = BIT(4) }; int mask = 0; /* Register the parsed histogram attributes */ size_t nattrs; @@ -685,15 +687,10 @@ parse_yaml_ellipsoid (filename, val, 4, 32768, &geom->data.ellipsoid.nslices); /* equivalent sphere radius */ - } else if(!strcmp((char*)key->data.scalar.value, "radius")) { - SETUP_MASK(RADIUS, "equivalent sphere radius of the ellipsoid"); - res = parse_yaml_param_distribution - (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.ellipsoid.radius); - - } else if(!strcmp((char*)key->data.scalar.value, "aspect_ratio")) { - SETUP_MASK(ASPECT_RATIO, "a/b aspect ratio of the ellipsoid"); - res = parse_yaml_double - (filename, val, DBL_MIN, DBL_MAX, &geom->data.ellipsoid.a_over_c); + } else if(!strcmp((char*)key->data.scalar.value, "radius_sphere")) { + SETUP_MASK(RADIUS_SPHERE, "equivalent sphere radius of the ellipsoid"); + res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, + &geom->data.ellipsoid.radius_sphere); /* Length of the ellipsoid "a" semi-principal axis */ } else if(!strcmp((char*)key->data.scalar.value, "a")) { @@ -719,44 +716,35 @@ parse_yaml_ellipsoid } /* Ensure the completeness of the ellipsoid distribution */ - if((mask & (LENGTH_a | LENGTH_c))) { - if(mask & (RADIUS | ASPECT_RATIO)) { - log_err(filename, node, - "the ellipsoid semi-principal axis and its equivalent sphere cannot " - "be defined simultaneously.\n"); - return RES_BAD_ARG; - } else if(!(mask & LENGTH_a)) { - log_err(filename, node, - "missing the length of the semi principal axis \"a\".\n"); - return RES_BAD_ARG; - } else if(!(mask & LENGTH_c)) { - log_err(filename, node, - "missing the length of the semi principal axis \"c\".\n"); - return RES_BAD_ARG; - } - geom->type = SCHIFF_ELLIPSOID; - } else if (mask & (RADIUS | ASPECT_RATIO)) { - if(mask & (LENGTH_a | LENGTH_c)) { - log_err(filename, node, - "the ellipsoid semi-principal axis and its equivalent sphere cannot " - "be defined simultaneously.\n"); - return RES_BAD_ARG; - } else if(!(mask & RADIUS)) { - log_err(filename, node, "missing the equivalent sphere radius.\n"); - return RES_BAD_ARG; - } else if(!(mask & ASPECT_RATIO)) { - log_err(filename, node, "missing the a/c ratio of the ellipsoid.\n"); - return RES_BAD_ARG; - } - geom->type = SCHIFF_ELLIPSOID_AS_SPHERE; + if(!(mask & LENGTH_a)) { + log_err(filename, node, + "missing the length of the semi principal axis \"a\".\n"); + return RES_BAD_ARG; + } else if(!(mask & LENGTH_c)) { + log_err(filename, node, + "missing the length of the semi principal axis \"c\".\n"); + return RES_BAD_ARG; } - + /* Setup the default values if required */ if(!(mask & PROBA)) { /* Default proba */ *geom_proba = 1.0; } if(!(mask & SLICES)) { /* Default number of slices */ geom->data.ellipsoid.nslices = SCHIFF_ELLIPSOID_DEFAULT.nslices; } + /* Define the geometry type */ + if(!(mask & RADIUS_SPHERE)) { + geom->type = SCHIFF_ELLIPSOID; + } else { + if(geom->data.ellipsoid.a.distribution != SCHIFF_PARAM_CONSTANT + || geom->data.ellipsoid.c.distribution != SCHIFF_PARAM_CONSTANT) { + log_err(filename, node, + "the radius_sphere parameter cannot be defined with non constant " + "ellipsoid attributes.\n"); + return RES_BAD_ARG; + } + geom->type = SCHIFF_ELLIPSOID_AS_SPHERE; + } return RES_OK; } @@ -772,7 +760,7 @@ parse_yaml_cylinder PROBA = BIT(0), RADIUS = BIT(1), HEIGHT = BIT(2), - ASPECT_RATIO = BIT(3), + RADIUS_SPHERE = BIT(3), SLICES = BIT(4) }; int mask = 0; /* Register the parsed attributes */ @@ -822,11 +810,11 @@ parse_yaml_cylinder res = parse_yaml_param_distribution (filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.cylinder.height); - /* Cylinder aspect ratio */ - } else if(!strcmp((char*)key->data.scalar.value, "aspect_ratio")) { - SETUP_MASK(ASPECT_RATIO, "cylinder height/radius aspect_ratio"); - res = parse_yaml_double - (filename, val, DBL_MIN, DBL_MAX, &geom->data.cylinder.aspect_ratio); + /* Equivalent sphere radius */ + } else if(!strcmp((char*)key->data.scalar.value, "radius_sphere")) { + SETUP_MASK(RADIUS_SPHERE, "equivalent sphere radius of the cylinder"); + res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, + &geom->data.cylinder.radius_sphere); /* # slices used to discretized the cylinder */ } else if(!strcmp((char*)key->data.scalar.value, "slices")) { @@ -849,27 +837,28 @@ parse_yaml_cylinder if(!(mask & RADIUS)) { log_err(filename, node, "missing the radius attribute.\n"); return RES_BAD_ARG; - } else if(!(mask & (HEIGHT|ASPECT_RATIO))) { - log_err(filename, node, "missing the height or aspect_ratio attribute.\n"); - return RES_BAD_ARG; - } else if((mask & HEIGHT) && (mask & ASPECT_RATIO)) { - log_err(filename, node, - "the height and the aspect_ratio attributes cannot be defined " - "simultanously.\n"); + } else if(!(mask & HEIGHT)) { + log_err(filename, node, "missing the height attribute.\n"); return RES_BAD_ARG; } - + /* Setup the default values if required */ if(!(mask & PROBA)) { /* Default proba */ *geom_proba = 1.0; } if(!(mask & SLICES)) { /* Default number of slices */ geom->data.cylinder.nslices = SCHIFF_CYLINDER_DEFAULT.nslices; } - - if(mask & HEIGHT) { + /* Define the geometry type */ + if(!(mask & RADIUS_SPHERE)) { geom->type = SCHIFF_CYLINDER; } else { - ASSERT(mask & ASPECT_RATIO); + if(geom->data.cylinder.radius.distribution != SCHIFF_PARAM_CONSTANT + || geom->data.cylinder.height.distribution != SCHIFF_PARAM_CONSTANT) { + log_err(filename, node, + "the radius_sphere parameter cannot be defined with non constant " + "cylinder attributes.\n"); + return RES_BAD_ARG; + } geom->type = SCHIFF_CYLINDER_AS_SPHERE; } return RES_OK; @@ -890,7 +879,8 @@ parse_yaml_helical_pipe RADIUS_HELICOID = BIT(3), RADIUS_CIRCLE = BIT(4), SLICES_HELICOID = BIT(5), - SLICES_CIRCLE = BIT(6) + SLICES_CIRCLE = BIT(6), + RADIUS_SPHERE = BIT(7) }; int mask = 0; /* Register the parsed attributes */ size_t nattrs; @@ -943,13 +933,13 @@ parse_yaml_helical_pipe } else if(!strcmp((char*)key->data.scalar.value, "radius_helicoid")) { SETUP_MASK(RADIUS_HELICOID, "helicoid radius"); res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, - &geom->data.helical_pipe.helicoid_radius); + &geom->data.helical_pipe.radius_helicoid); /* Radius of the meridian circle */ } else if(!strcmp((char*)key->data.scalar.value, "radius_circle")) { SETUP_MASK(RADIUS_CIRCLE, "circle radius of the helical pipe"); res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, - &geom->data.helical_pipe.circle_radius); + &geom->data.helical_pipe.radius_circle); /* # slices of the helicoid */ } else if(!strcmp((char*)key->data.scalar.value, "slices_helicoid")) { @@ -963,6 +953,12 @@ parse_yaml_helical_pipe res = parse_yaml_uint (filename, val, 4, 32768, &geom->data.helical_pipe.nslices_circle); + /* Equivalent sphere radius */ + } else if(!strcmp((char*)key->data.scalar.value, "radius_sphere")) { + SETUP_MASK(RADIUS_SPHERE, "equivalent sphere radius of the helical pipe"); + res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, + &geom->data.helical_pipe.radius_sphere); + /* Error */ } else { log_err(filename, key, "unknown helical pipe attribute `%s'.\n", @@ -987,8 +983,7 @@ parse_yaml_helical_pipe log_err(filename, node, "missing the radius of the meridian circle.\n"); return RES_BAD_ARG; } - - /* Setup default value if required */ + /* Setup the default values if required */ if(!(mask & PROBA)) { *geom_proba = 1.0; } @@ -1000,7 +995,22 @@ parse_yaml_helical_pipe geom->data.helical_pipe.nslices_circle = SCHIFF_HELICAL_PIPE_DEFAULT.nslices_circle; } - geom->type = SCHIFF_HELICAL_PIPE; + /* Define the geometry type */ + if(!(mask & RADIUS_SPHERE)) { + geom->type = SCHIFF_HELICAL_PIPE; + } else { + const struct schiff_helical_pipe* hpipe = &geom->data.helical_pipe; + if(hpipe->pitch.distribution != SCHIFF_PARAM_CONSTANT + || hpipe->height.distribution != SCHIFF_PARAM_CONSTANT + || hpipe->radius_helicoid.distribution != SCHIFF_PARAM_CONSTANT + || hpipe->radius_circle.distribution != SCHIFF_PARAM_CONSTANT) { + log_err(filename, node, + "the radius_sphere parameter cannot be defined with non constant " + "helical pipe attributes.\n"); + return RES_BAD_ARG; + } + geom->type = SCHIFF_HELICAL_PIPE_AS_SPHERE; + } return RES_OK; } @@ -1101,8 +1111,9 @@ parse_yaml_supershape enum { FORMULA0 = BIT(0), FORMULA1 = BIT(1), - PROBA = BIT(2), - SLICES = BIT(3) + RADIUS_SPHERE = BIT(2), + PROBA = BIT(3), + SLICES = BIT(4) }; int mask = 0; /* Register the parsed attributes */ size_t nattrs; @@ -1151,6 +1162,12 @@ parse_yaml_supershape res = parse_yaml_superformula (filename, doc, val, geom->data.supershape.formulas[1]); + /* Equivalent sphere radius */ + } else if(!strcmp((char*)key->data.scalar.value, "radius_sphere")) { + SETUP_MASK(RADIUS_SPHERE, "equivalent sphere radius of the supershape"); + res = parse_yaml_param_distribution(filename, doc, val, DBL_MIN, DBL_MAX, + &geom->data.supershape.radius_sphere); + /* # slices used to discretized the super shape */ } else if(!strcmp((char*)key->data.scalar.value, "slices")) { SETUP_MASK(SLICES, "supershape number of slices"); @@ -1176,14 +1193,29 @@ parse_yaml_supershape log_err(filename, node, "missing the formula1 attribute.\n"); return RES_BAD_ARG; } - + /* Setup the default values if required */ if(!(mask & PROBA)) { /* Default proba */ *geom_proba = 1.0; } if(!(mask & SLICES)) { /* Default number of slices */ geom->data.supershape.nslices = SCHIFF_SUPERSHAPE_DEFAULT.nslices; } - geom->type = SCHIFF_SUPERSHAPE; + /* Define the geometry type */ + if(!(mask & RADIUS_SPHERE)) { + geom->type = SCHIFF_SUPERSHAPE; + } else { + const struct schiff_supershape* sshape = &geom->data.supershape; + FOR_EACH(i, 0, 6) { + if(sshape->formulas[0][i].distribution != SCHIFF_PARAM_CONSTANT + || sshape->formulas[1][i].distribution != SCHIFF_PARAM_CONSTANT) { + log_err(filename, node, + "the radius_sphere parameter cannot be defined with non constant " + "supershape attributes.\n"); + return RES_BAD_ARG; + } + } + geom->type = SCHIFF_SUPERSHAPE_AS_SPHERE; + } return RES_OK; } diff --git a/src/schiff_geometry.c b/src/schiff_geometry.c @@ -51,8 +51,8 @@ struct cylinder_context { /* 3D Context of an helical pipe */ struct helical_pipe_context { - double helicoid_radius; - double circle_radius; + double radius_helicoid; + double radius_circle; double pitch; double height; }; @@ -83,11 +83,17 @@ struct mesh_context { } data; }; +struct shape { + const struct schiff_geometry* geometry; /* Pointer onto the geometry */ + struct s3d_shape* shape; /* May be NULL => No volume distribution */ + struct schiff_mesh mesh; /* Triangular mesh of the shape */ + double volume; /* Volume of the shape */ +}; + /* Distribution context of a generic geometry */ struct geometry_distribution_context { struct schiff_optical_properties* properties; /* Per wavelength properties */ - struct schiff_mesh* meshes; /* List of triangular meshes */ - struct schiff_geometry const ** geometries; /* List of geometries */ + struct shape* shapes; /* List of shapes */ struct ssp_ran_discrete* ran_geometries; /* Geometries random variates */ }; @@ -156,7 +162,7 @@ static INLINE double eval_param(const struct schiff_param* param, struct ssp_rng* rng) { double val = 0.0; - ASSERT(param && rng); + ASSERT(param); switch(param->distribution) { case SCHIFF_PARAM_CONSTANT: @@ -289,199 +295,579 @@ supershape_get_position(const unsigned ivert, float vertex[3], void* ctx) } static res_T -geometry_sample_ellipsoid - (struct ssp_rng* rng, - const struct schiff_geometry* geom, - const struct schiff_mesh* mesh, - struct s3d_shape* shape) +compute_s3d_shape_volume + (struct s3d_device* s3d, struct s3d_shape* shape, double* out_volume) +{ + struct s3d_scene* scn = NULL; + float volume = 0.f; + int mask; + res_T res = RES_OK; + ASSERT(s3d && shape && out_volume); + + if(RES_OK != (res = s3d_scene_create(s3d, &scn))) goto error; + if(RES_OK != (res = s3d_scene_attach_shape(scn, shape))) goto error; + if(RES_OK != (res = s3d_scene_begin_session(scn, S3D_TRACE))) goto error; + if(RES_OK != (res = s3d_scene_compute_volume(scn, &volume))) goto error; + S3D(scene_end_session(scn)); + +exit: + if(scn) S3D(scene_ref_put(scn)); + /* The volume may be negative if the faces are not correctly oriented */ + *out_volume = absf(volume); + return res; +error: + S3D(scene_get_session_mask(scn, &mask)); + if(mask) S3D(scene_end_session(scn)); + goto exit; +} + +static void +shape_release(struct shape* shape) +{ + ASSERT(shape); + schiff_mesh_release(&shape->mesh); + if(shape->shape) S3D(shape_ref_put(shape->shape)); +} + +static res_T +shape_cylinder_init + (struct shape* shape, + struct s3d_device* s3d, + const struct schiff_geometry* geometry) +{ + struct mesh_context ctx; + struct s3d_vertex_data attrib; + const struct schiff_cylinder* cylinder; + double radius, height; + size_t nverts, nprims; + res_T res = RES_OK; + ASSERT(shape && geometry); + ASSERT(geometry->type == SCHIFF_CYLINDER + || geometry->type == SCHIFF_CYLINDER_AS_SPHERE); + + memset(shape, 0, sizeof(*shape)); + cylinder = &geometry->data.cylinder; + + + res = schiff_mesh_init_cylinder + (&mem_default_allocator, &shape->mesh, cylinder->nslices); + if(res != RES_OK) goto error; + + shape->geometry = geometry; + + if(geometry->type == SCHIFF_CYLINDER) /* That's all folks */ + goto exit; + + /* Create the Star-3D cylinder shape */ + res = s3d_shape_create_mesh(s3d, &shape->shape); + if(res != RES_OK) goto error; + + /* Cylinder parameters must be constant */ + ASSERT(cylinder->radius.distribution == SCHIFF_PARAM_CONSTANT); + ASSERT(cylinder->height.distribution == SCHIFF_PARAM_CONSTANT); + radius = cylinder->radius.data.constant; + height = cylinder->height.data.constant; + + /* Setup the context of the cylinder */ + ctx.type = SCHIFF_CYLINDER; + ctx.mesh = &shape->mesh; + ctx.data.cylinder.radius = (float)radius; + ctx.data.cylinder.height = (float)height; + + /* Define the position vertex attribute */ + attrib.usage = S3D_POSITION; + attrib.type = S3D_FLOAT3; + attrib.get = cylinder_get_position; + + nverts = darray_float_size_get(&shape->mesh.vertices.cartesian) / 3/*#coords*/; + nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices per prim*/; + + res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, + geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); + if(res != RES_OK) goto error; + + shape->volume = PI*height*radius*radius; /* Analytic volume */ + + schiff_mesh_release(&shape->mesh); + +exit: + return res; +error: + shape_release(shape); + goto exit; +} + +static res_T +shape_cylinder_generate_s3d_shape + (const struct shape* shape, + struct ssp_rng* rng, + double* volume_scaling, + struct s3d_shape* s3d_shape) +{ + struct s3d_vertex_data attrib; + struct mesh_context ctx; + size_t nverts, nprims; + ASSERT(shape && volume_scaling); + ASSERT(shape->geometry->type == SCHIFF_CYLINDER); + + ctx.type = SCHIFF_CYLINDER; + ctx.mesh = &shape->mesh; + ctx.data.cylinder.radius = (float)eval_param + (&shape->geometry->data.cylinder.radius, rng); + ctx.data.cylinder.height = (float)eval_param + (&shape->geometry->data.cylinder.height, rng); + + attrib.usage = S3D_POSITION; + attrib.type = S3D_FLOAT3; + attrib.get = cylinder_get_position; + + nverts = darray_float_size_get(&shape->mesh.vertices.cartesian) / 3/*#coords*/; + nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices per prim*/; + + *volume_scaling = 1.0; + return s3d_mesh_setup_indexed_vertices(s3d_shape, (unsigned)nprims, + geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); +} + +static res_T +shape_cylinder_as_sphere_generate_s3d_shape + (const struct shape* shape, + struct ssp_rng* rng, + double* volume_scaling, + struct s3d_shape* s3d_shape) +{ + double radius, sphere_volume; + ASSERT(shape && volume_scaling); + ASSERT(shape->geometry->type == SCHIFF_CYLINDER_AS_SPHERE); + + radius = eval_param(&shape->geometry->data.cylinder.radius_sphere, rng); + sphere_volume = 4.0/3.0 * PI * radius*radius*radius; + *volume_scaling = sphere_volume / shape->volume; + return s3d_mesh_copy(shape->shape, s3d_shape); +} + +static res_T +shape_ellipsoid_init + (struct shape* shape, + struct s3d_device* s3d, + const struct schiff_geometry* geometry) { const struct schiff_ellipsoid* ellipsoid; - struct mesh_context mesh_ctx; + struct mesh_context ctx; struct s3d_vertex_data attrib; - double a, b, c, a2, b2, c2, r_eq; + double a, b, c, a2, b2, c2; size_t nverts, nprims; - ASSERT(rng && geom && mesh && shape); + res_T res = RES_OK; + ASSERT(shape && geometry); + ASSERT(geometry->type == SCHIFF_ELLIPSOID + || geometry->type == SCHIFF_ELLIPSOID_AS_SPHERE); - ellipsoid = &geom->data.ellipsoid; - switch(geom->type) { - case SCHIFF_ELLIPSOID: - a = eval_param(&ellipsoid->a, rng); - c = eval_param(&ellipsoid->c, rng); - break; - case SCHIFF_ELLIPSOID_AS_SPHERE: - r_eq = eval_param(&ellipsoid->radius, rng); - /* a_over_c^(1/3) could be precomputed */ - a = r_eq * pow(ellipsoid->a_over_c, 1.0/3.0); - c = a / ellipsoid->a_over_c; - break; - default: FATAL("Unreachable code.\n"); break; - } + ellipsoid = &geometry->data.ellipsoid; + memset(shape, 0, sizeof(*shape)); + + res = schiff_mesh_init_sphere_polar + (&mem_default_allocator, &shape->mesh, ellipsoid->nslices); + if(res != RES_OK) goto error; + + shape->geometry = geometry; - /* Schiff exposes ellipsoid whose a and b semi-principal axis are equals */ + if(geometry->type == SCHIFF_ELLIPSOID) /* That's all folks */ + goto exit; + + res = s3d_shape_create_mesh(s3d, &shape->shape); + if(res != RES_OK) goto error; + + /* Expecting constant parameters */ + ASSERT(ellipsoid->a.distribution == SCHIFF_PARAM_CONSTANT); + ASSERT(ellipsoid->c.distribution == SCHIFF_PARAM_CONSTANT); + a = ellipsoid->a.data.constant; + c = ellipsoid->c.data.constant; b = a; a2 = a*a; b2 = b*b; c2 = c*c; - mesh_ctx.mesh = mesh; - mesh_ctx.type = SCHIFF_ELLIPSOID; - mesh_ctx.data.ellipsoid.abc = a*b*c; - mesh_ctx.data.ellipsoid.b2c2 = b2*c2; - mesh_ctx.data.ellipsoid.a2c2 = a2*c2; - mesh_ctx.data.ellipsoid.a2b2 = a2*b2; + ctx.mesh = &shape->mesh; + ctx.type = SCHIFF_ELLIPSOID; + ctx.data.ellipsoid.abc = a*b*c; + ctx.data.ellipsoid.b2c2 = b2*c2; + ctx.data.ellipsoid.a2c2 = a2*c2; + ctx.data.ellipsoid.a2b2 = a2*b2; attrib.usage = S3D_POSITION; attrib.type = S3D_FLOAT3; attrib.get = ellipsoid_get_position; - nverts = darray_uint_size_get(&mesh->vertices.polar) / 2/*#theta/phi*/; - nprims = darray_uint_size_get(&mesh->indices) / 3/*#indices*/; + nverts = darray_uint_size_get(&shape->mesh.vertices.polar) / 2/*#theta/phi*/; + nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; + + res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, + geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); + if(res != RES_OK) goto error; + + shape->volume = 4.0/3.0 * PI * ctx.data.ellipsoid.abc; + + schiff_mesh_release(&shape->mesh); - return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, - geometry_get_indices, (unsigned)nverts, &attrib, 1, &mesh_ctx); +exit: + return res; +error: + shape_release(shape); + goto exit; } static res_T -geometry_sample_cylinder - (struct ssp_rng* rng, - const struct schiff_geometry* geom, - const struct schiff_mesh* mesh, - struct s3d_shape* shape) +shape_ellipsoid_generate_s3d_shape + (const struct shape* shape, + struct ssp_rng* rng, + double* volume_scaling, + struct s3d_shape* s3d_shape) { - const struct schiff_cylinder* cylinder; - struct mesh_context mesh_ctx; struct s3d_vertex_data attrib; + struct mesh_context ctx; size_t nverts, nprims; - double r; - ASSERT(rng && geom && mesh && shape); + double a, b, c, a2, b2, c2; + ASSERT(shape && volume_scaling); + ASSERT(shape->geometry->type == SCHIFF_ELLIPSOID); - cylinder = &geom->data.cylinder; - mesh_ctx.type = SCHIFF_CYLINDER; - mesh_ctx.mesh = mesh; - switch(geom->type) { - case SCHIFF_CYLINDER: - mesh_ctx.data.cylinder.radius = (float)eval_param(&cylinder->radius, rng); - mesh_ctx.data.cylinder.height = (float)eval_param(&cylinder->height, rng); - break; - case SCHIFF_CYLINDER_AS_SPHERE: - r = eval_param(&cylinder->radius, rng); - mesh_ctx.data.cylinder.radius = /* The denominator can be pre-computed */ - (float)(r / pow(3.0 / (2.0*cylinder->aspect_ratio), 1.0/3.0)); - mesh_ctx.data.cylinder.height = - (float)(2.f * mesh_ctx.data.cylinder.radius / cylinder->aspect_ratio); - break; - default: FATAL("Unreachable code.\n"); break; - } + ctx.type = SCHIFF_ELLIPSOID; + ctx.mesh = &shape->mesh; + a = eval_param(&shape->geometry->data.ellipsoid.a, rng); + c = eval_param(&shape->geometry->data.ellipsoid.c, rng); + + b = a; + a2 = a*a; + b2 = b*b; + c2 = c*c; + + ctx.mesh = &shape->mesh; + ctx.type = SCHIFF_ELLIPSOID; + ctx.data.ellipsoid.abc = a*b*c; + ctx.data.ellipsoid.b2c2 = b2*c2; + ctx.data.ellipsoid.a2c2 = a2*c2; + ctx.data.ellipsoid.a2b2 = a2*b2; attrib.usage = S3D_POSITION; attrib.type = S3D_FLOAT3; - attrib.get = cylinder_get_position; + attrib.get = ellipsoid_get_position; - nverts = darray_float_size_get(&mesh->vertices.cartesian) / 3/*#coords*/; - nprims = darray_uint_size_get(&mesh->indices) / 3/*#indices per prim*/; + nverts = darray_uint_size_get(&shape->mesh.vertices.polar) / 2/*#theta/phi*/; + nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; - return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, - geometry_get_indices, (unsigned)nverts, &attrib, 1, &mesh_ctx); + *volume_scaling = 1.0; + return s3d_mesh_setup_indexed_vertices(s3d_shape, (unsigned)nprims, + geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); } static res_T -geometry_sample_helical_pipe - (struct ssp_rng* rng, - const struct schiff_geometry* geom, - struct schiff_mesh* mesh, - struct s3d_shape* shape) +shape_ellipsoid_as_sphere_generate_s3d_shape + (const struct shape* shape, + struct ssp_rng* rng, + double* volume_scaling, + struct s3d_shape* s3d_shape) +{ + double radius, sphere_volume; + ASSERT(shape && volume_scaling); + ASSERT(shape->geometry->type == SCHIFF_ELLIPSOID_AS_SPHERE); + + radius = eval_param(&shape->geometry->data.ellipsoid.radius_sphere, rng); + sphere_volume = 4.0/3.0 * PI * radius*radius*radius; + *volume_scaling = sphere_volume / shape->volume; + return s3d_mesh_copy(shape->shape, s3d_shape); +} + +static res_T +shape_helical_pipe_init + (struct shape* shape, + struct s3d_device* s3d, + const struct schiff_geometry* geometry) { const struct schiff_helical_pipe* helical_pipe; - struct mesh_context mesh_ctx; + struct mesh_context ctx; struct s3d_vertex_data attrib; - size_t nprims; - size_t nverts; - double pitch; - double height; - double hradius; - double cradius; + size_t nprims, nverts; + double pitch, height, hradius, cradius; + res_T res = RES_OK; + + ASSERT(geometry); + ASSERT(geometry->type == SCHIFF_HELICAL_PIPE + || geometry->type == SCHIFF_HELICAL_PIPE_AS_SPHERE); + + helical_pipe = &geometry->data.helical_pipe; + memset(shape, 0, sizeof(*shape)); + + res = schiff_mesh_init_helical_pipe(&mem_default_allocator, &shape->mesh, + helical_pipe->nslices_helicoid, helical_pipe->nslices_circle); + if(res != RES_OK) goto error; + + shape->geometry = geometry; + + if(geometry->type == SCHIFF_HELICAL_PIPE) /* That's all folks */ + goto exit; + + res = s3d_shape_create_mesh(s3d, &shape->shape); + if(res != RES_OK) goto error; + + /* Expecting constant parameters */ + ASSERT(helical_pipe->radius_helicoid.distribution == SCHIFF_PARAM_CONSTANT); + ASSERT(helical_pipe->radius_circle.distribution == SCHIFF_PARAM_CONSTANT); + ASSERT(helical_pipe->pitch.distribution == SCHIFF_PARAM_CONSTANT); + ASSERT(helical_pipe->height.distribution == SCHIFF_PARAM_CONSTANT); + hradius = helical_pipe->radius_helicoid.data.constant; + cradius = helical_pipe->radius_circle.data.constant; + pitch = helical_pipe->pitch.data.constant; + height = helical_pipe->height.data.constant; + + /* Setup the mesh vertices */ + res = schiff_mesh_setup_helical_pipe(&shape->mesh, pitch, height, hradius, + cradius, helical_pipe->nslices_helicoid, helical_pipe->nslices_circle); + if(res != RES_OK) return res; + + ctx.type = SCHIFF_HELICAL_PIPE; + ctx.mesh = &shape->mesh; + + attrib.usage = S3D_POSITION; + attrib.type = S3D_FLOAT3; + attrib.get = geometry_get_position; + + nverts = darray_float_size_get(&shape->mesh.vertices.cartesian) / 3/*#coords*/; + nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices per prim*/; + + res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, + geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); + if(res != RES_OK) goto error; + + res = compute_s3d_shape_volume(s3d, shape->shape, &shape->volume); + if(res != RES_OK) goto error; + + schiff_mesh_release(&shape->mesh); + +exit: + return res; +error: + shape_release(shape); + goto exit; +} + +static res_T +shape_helical_pipe_generate_s3d_shape + (struct shape* shape, + struct ssp_rng* rng, + double* volume_scaling, + struct s3d_shape* s3d_shape) +{ + const struct schiff_helical_pipe* helical_pipe; + struct mesh_context ctx; + struct s3d_vertex_data attrib; + size_t nprims, nverts; + double pitch, height, hradius, cradius; res_T res = RES_OK; - ASSERT(geom && geom->type == SCHIFF_HELICAL_PIPE); - helical_pipe = &geom->data.helical_pipe; + ASSERT(shape && volume_scaling); + ASSERT(shape->geometry->type == SCHIFF_HELICAL_PIPE); + + helical_pipe = &shape->geometry->data.helical_pipe; /* Sample the helical pipe attributes */ - hradius = eval_param(&helical_pipe->helicoid_radius, rng); - cradius = eval_param(&helical_pipe->circle_radius, rng); + hradius = eval_param(&helical_pipe->radius_helicoid, rng); + cradius = eval_param(&helical_pipe->radius_circle, rng); pitch = eval_param(&helical_pipe->pitch, rng); height = eval_param(&helical_pipe->height, rng); /* Setup the mesh vertices */ - res = schiff_mesh_setup_helical_pipe(mesh, pitch, height, hradius, cradius, - helical_pipe->nslices_helicoid, helical_pipe->nslices_circle); + res = schiff_mesh_setup_helical_pipe(&shape->mesh, pitch, height, hradius, + cradius, helical_pipe->nslices_helicoid, helical_pipe->nslices_circle); if(res != RES_OK) return res; - mesh_ctx.mesh = mesh; - mesh_ctx.type = SCHIFF_HELICAL_PIPE; + ctx.mesh = &shape->mesh; + ctx.type = SCHIFF_HELICAL_PIPE; attrib.usage = S3D_POSITION; attrib.type = S3D_FLOAT3; attrib.get = geometry_get_position; - nverts = darray_float_size_get(&mesh->vertices.cartesian) / 3/*#coords*/; - nprims = darray_uint_size_get(&mesh->indices) / 3/*#indices per prim*/; + nverts = darray_float_size_get(&shape->mesh.vertices.cartesian) / 3/*#coords*/; + nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices per prim*/; - return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, - geometry_get_indices, (unsigned)nverts, &attrib, 1, &mesh_ctx); + *volume_scaling = 1.0; + return s3d_mesh_setup_indexed_vertices(s3d_shape, (unsigned)nprims, + geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); } static res_T -geometry_sample_sphere - (struct ssp_rng* rng, - const struct schiff_geometry* geom, - const struct schiff_mesh* mesh, - struct s3d_shape* shape) +shape_helical_pipe_as_sphere_generate_s3d_shape + (const struct shape* shape, + struct ssp_rng* rng, + double* volume_scaling, + struct s3d_shape* s3d_shape) { - const struct schiff_sphere* sphere; - struct mesh_context mesh_ctx; + double radius, sphere_volume; + ASSERT(shape && volume_scaling); + ASSERT(shape->geometry->type == SCHIFF_HELICAL_PIPE_AS_SPHERE); + + radius = eval_param(&shape->geometry->data.helical_pipe.radius_sphere, rng); + sphere_volume = 4.0/3.0 * PI * radius*radius*radius; + *volume_scaling = sphere_volume / shape->volume; + return s3d_mesh_copy(shape->shape, s3d_shape); +} + +static res_T +shape_sphere_init + (struct shape* shape, + struct s3d_device* s3d, + const struct schiff_geometry* geometry) +{ + struct mesh_context ctx; struct s3d_vertex_data attrib; + const struct schiff_sphere* sphere; size_t nverts, nprims; - ASSERT(geom && mesh && shape); - ASSERT(geom->type == SCHIFF_SPHERE); + res_T res = RES_OK; + ASSERT(shape && geometry && geometry->type == SCHIFF_SPHERE); + + sphere = &geometry->data.sphere; + memset(shape, 0, sizeof(*shape)); + + /* Generate the sphere mesh */ + res = schiff_mesh_init_sphere + (&mem_default_allocator, &shape->mesh, sphere->nslices); + if(res != RES_OK) goto error; + + shape->geometry = geometry; - sphere = &geom->data.sphere; - mesh_ctx.mesh = mesh; - mesh_ctx.type = SCHIFF_SPHERE; - mesh_ctx.data.sphere.radius = (float)eval_param(&sphere->radius, rng); + /* Create the Star-3D sphere shape */ + res = s3d_shape_create_mesh(s3d, &shape->shape); + if(res != RES_OK) goto error; + + /* Setup the context of the sphere*/ + ctx.type = SCHIFF_SPHERE; + ctx.mesh = &shape->mesh; + ctx.data.sphere.radius = 1.f; + /* Define the position vertex attribute */ attrib.usage = S3D_POSITION; attrib.type = S3D_FLOAT3; attrib.get = sphere_get_position; - nverts = darray_float_size_get(&mesh->vertices.cartesian) / 3/*#coords*/; - nprims = darray_uint_size_get(&mesh->indices) / 3/*#indices*/; + nverts = darray_float_size_get(&shape->mesh.vertices.cartesian) / 3/*#coords*/; + nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; + + res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, + geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); + if(res != RES_OK) goto error; + + shape->volume = 4.0/3.0*PI; /* Analytic volume. The radius is implicitly 1 */ + +exit: + schiff_mesh_release(&shape->mesh); + return res; +error: + shape_release(shape); + goto exit; +} + +static res_T +shape_sphere_generate_s3d_shape + (const struct shape* shape, + struct ssp_rng* rng, + double* volume_scaling, + struct s3d_shape* s3d_shape) +{ + double radius, sphere_volume; + ASSERT(shape && volume_scaling); + ASSERT(shape->geometry->type == SCHIFF_SPHERE); + + radius = eval_param(&shape->geometry->data.sphere.radius, rng); + sphere_volume = 4.0/3.0 * PI * radius*radius*radius; + *volume_scaling = sphere_volume / shape->volume; + return s3d_mesh_copy(shape->shape, s3d_shape); +} + +static res_T +shape_supershape_init + (struct shape* shape, + struct s3d_device* s3d, + const struct schiff_geometry* geometry) +{ + const struct schiff_supershape* sshape; + struct mesh_context ctx; + struct s3d_vertex_data attrib; + size_t nverts, nprims; + int iform, iattr; + res_T res = RES_OK; + ASSERT(geometry); + ASSERT(geometry->type == SCHIFF_SUPERSHAPE + || geometry->type == SCHIFF_SUPERSHAPE_AS_SPHERE); + + sshape = &geometry->data.supershape; + memset(shape, 0, sizeof(*shape)); + + res = schiff_mesh_init_sphere_polar(&mem_default_allocator, + &shape->mesh, geometry->data.supershape.nslices); + if(res != RES_OK) goto error; + + shape->geometry = geometry; + + if(geometry->type == SCHIFF_SUPERSHAPE) /* That's all folks */ + goto exit; + + res = s3d_shape_create_mesh(s3d, &shape->shape); + if(res != RES_OK) goto error; + + ctx.type = SCHIFF_SUPERSHAPE; + ctx.mesh = &shape->mesh; + FOR_EACH(iform, 0, 2) { + FOR_EACH(iattr, 0, 6) { + /* Expecting constant parameters */ + ASSERT(sshape->formulas[iform][iattr].distribution == SCHIFF_PARAM_CONSTANT); + ctx.data.supershape.formulas[iform][iattr] = + sshape->formulas[iform][iattr].data.constant; + }} + + attrib.usage = S3D_POSITION; + attrib.type = S3D_FLOAT3; + attrib.get = supershape_get_position; + + nverts = darray_uint_size_get(&shape->mesh.vertices.polar) / 2/*#theta/phi*/; + nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; + + res = s3d_mesh_setup_indexed_vertices(shape->shape, (unsigned)nprims, + geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); + if(res != RES_OK) goto error; + + res = compute_s3d_shape_volume(s3d, shape->shape, &shape->volume); + if(res != RES_OK) goto error; + + schiff_mesh_release(&shape->mesh); - return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, - geometry_get_indices, (unsigned)nverts, &attrib, 1, &mesh_ctx); +exit: + return res; +error: + shape_release(shape); + goto exit; } static res_T -geometry_sample_supershape - (struct ssp_rng* rng, - const struct schiff_geometry* geom, - const struct schiff_mesh* mesh, - struct s3d_shape* shape) +shape_supershape_generate_s3d_shape + (const struct shape* shape, + struct ssp_rng* rng, + double* volume_scaling, + struct s3d_shape* s3d_shape) { const struct schiff_supershape* sshape; - struct mesh_context mesh_ctx; + struct mesh_context ctx; struct s3d_vertex_data attrib; size_t nverts, nprims; int iform, iattr; - ASSERT(geom && mesh && shape); - ASSERT(geom->type == SCHIFF_SUPERSHAPE); + ASSERT(shape && volume_scaling); + ASSERT(shape->geometry->type == SCHIFF_SUPERSHAPE); - sshape = &geom->data.supershape; - mesh_ctx.mesh = mesh; - mesh_ctx.type = SCHIFF_SUPERSHAPE; + sshape = &shape->geometry->data.supershape; + ctx.mesh = &shape->mesh; + ctx.type = SCHIFF_SUPERSHAPE; FOR_EACH(iform, 0, 2) { FOR_EACH(iattr, 0, 6) { - mesh_ctx.data.supershape.formulas[iform][iattr] = + ctx.data.supershape.formulas[iform][iattr] = (float)eval_param(&sshape->formulas[iform][iattr], rng); }} @@ -489,47 +875,83 @@ geometry_sample_supershape attrib.type = S3D_FLOAT3; attrib.get = supershape_get_position; - nverts = darray_uint_size_get(&mesh->vertices.polar) / 2/*#theta/phi*/; - nprims = darray_uint_size_get(&mesh->indices) / 3/*#indices*/; + nverts = darray_uint_size_get(&shape->mesh.vertices.polar) / 2/*#theta/phi*/; + nprims = darray_uint_size_get(&shape->mesh.indices) / 3/*#indices*/; - return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, - geometry_get_indices, (unsigned)nverts, &attrib, 1, &mesh_ctx); + *volume_scaling = 1.0; + return s3d_mesh_setup_indexed_vertices(s3d_shape, (unsigned)nprims, + geometry_get_indices, (unsigned)nverts, &attrib, 1, &ctx); +} + +static res_T +shape_supershape_as_sphere_generate_s3d_shape + (const struct shape* shape, + struct ssp_rng* rng, + double* volume_scaling, + struct s3d_shape* s3d_shape) +{ + double radius, sphere_volume; + ASSERT(shape && volume_scaling); + ASSERT(shape->geometry->type == SCHIFF_SUPERSHAPE_AS_SPHERE); + + radius = eval_param(&shape->geometry->data.supershape.radius_sphere, rng); + sphere_volume = 4.0/3.0 * PI * radius*radius*radius; + *volume_scaling = sphere_volume / shape->volume; + return s3d_mesh_copy(shape->shape, s3d_shape); } static res_T geometry_sample (struct ssp_rng* rng, - struct s3d_shape* shape, + struct s3d_shape* s3d_shape, + double* volume_scaling, void* ctx) { struct geometry_distribution_context* distrib = ctx; - const struct schiff_geometry* geom; - struct schiff_mesh* mesh; + struct shape* shape; size_t isamp; res_T res = RES_OK; - ASSERT(rng && shape && ctx); + ASSERT(rng && ctx); isamp = ssp_ran_discrete(rng, distrib->ran_geometries); - geom = distrib->geometries[isamp]; - mesh = distrib->meshes + isamp; + shape = distrib->shapes + isamp; - switch(geom->type) { + switch(shape->geometry->type) { case SCHIFF_ELLIPSOID: + res = shape_ellipsoid_generate_s3d_shape + (shape, rng, volume_scaling, s3d_shape); + break; case SCHIFF_ELLIPSOID_AS_SPHERE: - res = geometry_sample_ellipsoid(rng, geom, mesh, shape); + res = shape_ellipsoid_as_sphere_generate_s3d_shape + (shape, rng, volume_scaling, s3d_shape); break; case SCHIFF_CYLINDER: + res = shape_cylinder_generate_s3d_shape + (shape, rng, volume_scaling, s3d_shape); + break; case SCHIFF_CYLINDER_AS_SPHERE: - res = geometry_sample_cylinder(rng, geom, mesh, shape); + res = shape_cylinder_as_sphere_generate_s3d_shape + (shape, rng, volume_scaling, s3d_shape); break; case SCHIFF_HELICAL_PIPE: - res = geometry_sample_helical_pipe(rng, geom, mesh, shape); + res = shape_helical_pipe_generate_s3d_shape + (shape, rng, volume_scaling, s3d_shape); + break; + case SCHIFF_HELICAL_PIPE_AS_SPHERE: + res = shape_helical_pipe_as_sphere_generate_s3d_shape + (shape, rng, volume_scaling, s3d_shape); break; case SCHIFF_SPHERE: - res = geometry_sample_sphere(rng, geom, mesh, shape); + res = shape_sphere_generate_s3d_shape + (shape, rng, volume_scaling, s3d_shape); break; case SCHIFF_SUPERSHAPE: - res = geometry_sample_supershape(rng, geom, mesh, shape); + res = shape_supershape_generate_s3d_shape + (shape, rng, volume_scaling, s3d_shape); + break; + case SCHIFF_SUPERSHAPE_AS_SPHERE: + res = shape_supershape_as_sphere_generate_s3d_shape + (shape, rng, volume_scaling, s3d_shape); break; default: FATAL("Unreachable code.\n"); break; } @@ -554,13 +976,12 @@ schiff_geometry_distribution_release if(!distrib->context) return; ctx = distrib->context; - if(ctx->meshes) { - count = sa_size(ctx->meshes); - FOR_EACH(i, 0, count) schiff_mesh_release(&ctx->meshes[i]); - sa_release(ctx->meshes); + if(ctx->shapes) { + count = sa_size(ctx->shapes); + FOR_EACH(i, 0, count) + shape_release(&ctx->shapes[i]); + sa_release(ctx->shapes); } - if(ctx->geometries) - sa_release(ctx->geometries); if(ctx->ran_geometries) SSP(ran_discrete_ref_put(ctx->ran_geometries)); mem_rm(ctx); @@ -569,6 +990,7 @@ schiff_geometry_distribution_release res_T schiff_geometry_distribution_init (struct sschiff_geometry_distribution* distrib, + struct s3d_device* s3d, const struct schiff_geometry* geoms, const size_t ngeoms, const double characteristic_length, @@ -593,56 +1015,43 @@ schiff_geometry_distribution_init /* Directly setyp the distribution context to avoid memory leaks on error */ distrib->context = ctx; - if(!sa_add(ctx->meshes, ngeoms)) { + if(!sa_add(ctx->shapes, ngeoms)) { fprintf(stderr, - "Couldn't allocate the meshes of the geometry distributions.\n"); - res = RES_MEM_ERR; - goto error; - } - memset(ctx->meshes, 0, sizeof(ctx->meshes[0])*ngeoms); - - if(!sa_add(ctx->geometries, ngeoms)) { - fprintf(stderr, -"Couldn't allocate the array that reference the geometry distributions.\n"); + "Couldn't allocate the shapes of the geometry distributions.\n"); res = RES_MEM_ERR; goto error; } + memset(ctx->shapes, 0, sizeof(ctx->shapes[0])*ngeoms); FOR_EACH(i, 0, ngeoms) { - /* Generate the mesh template and setup its distribution context */ switch(geoms[i].type) { case SCHIFF_ELLIPSOID: case SCHIFF_ELLIPSOID_AS_SPHERE: - res = schiff_mesh_init_sphere_polar(&mem_default_allocator, - &ctx->meshes[i], geoms[i].data.ellipsoid.nslices); + res = shape_ellipsoid_init(&ctx->shapes[i], s3d, &geoms[i]); break; case SCHIFF_CYLINDER: case SCHIFF_CYLINDER_AS_SPHERE: - res = schiff_mesh_init_cylinder(&mem_default_allocator, &ctx->meshes[i], - geoms[i].data.cylinder.nslices); + res = shape_cylinder_init(&ctx->shapes[i], s3d, &geoms[i]); break; case SCHIFF_HELICAL_PIPE: - res = schiff_mesh_init_helical_pipe - (&mem_default_allocator, &ctx->meshes[i], - geoms[i].data.helical_pipe.nslices_helicoid, - geoms[i].data.helical_pipe.nslices_circle); + case SCHIFF_HELICAL_PIPE_AS_SPHERE: + res = shape_helical_pipe_init(&ctx->shapes[i], s3d, &geoms[i]); break; case SCHIFF_SPHERE: - res = schiff_mesh_init_sphere(&mem_default_allocator, &ctx->meshes[i], - geoms[i].data.sphere.nslices); + res = shape_sphere_init(&ctx->shapes[i], s3d, &geoms[i]); break; case SCHIFF_SUPERSHAPE: - res = schiff_mesh_init_sphere_polar(&mem_default_allocator, - &ctx->meshes[i], geoms[i].data.supershape.nslices); + case SCHIFF_SUPERSHAPE_AS_SPHERE: + res = shape_supershape_init(&ctx->shapes[i], s3d, &geoms[i]); break; default: FATAL("Unreachable code.\n"); break; } if(res != RES_OK) { - fprintf(stderr, "Couldn't initialise the mesh template.\n"); + fprintf(stderr, + "Couldn't initialise the shape of the geometry distribution.\n"); goto error; } - ctx->geometries[i] = geoms + i; } ctx->properties = properties; SSP(ran_discrete_ref_get(ran_geoms)); diff --git a/src/schiff_geometry.h b/src/schiff_geometry.h @@ -39,32 +39,34 @@ struct schiff_param { #define SCHIFF_PARAM_DEFAULT__ {SCHIFF_PARAM_CONSTANT, {1.0}} enum schiff_geometry_type { - SCHIFF_ELLIPSOID, - SCHIFF_ELLIPSOID_AS_SPHERE,/* The ellipsoid volume is controled by a sphere */ SCHIFF_CYLINDER, - SCHIFF_CYLINDER_AS_SPHERE,/* The cylinder volume is controled by a sphere */ + SCHIFF_ELLIPSOID, SCHIFF_HELICAL_PIPE, SCHIFF_SPHERE, SCHIFF_SUPERSHAPE, + + /* Volume is controlled by a sphere */ + SCHIFF_CYLINDER_AS_SPHERE, + SCHIFF_ELLIPSOID_AS_SPHERE, + SCHIFF_HELICAL_PIPE_AS_SPHERE, + SCHIFF_SUPERSHAPE_AS_SPHERE, + SCHIFF_NONE }; /* (x/a)^2 + (y/a)^2 + (z/c)^2 = 1 */ struct schiff_ellipsoid { - /* In use by SCHIFF_ELLIPSOID */ struct schiff_param a; struct schiff_param c; - /* In use by SCHIFF_ELLIPSOID_AS_SPHERE - * a = radius * a_over_c^1/3; c = a / a_over_c */ - struct schiff_param radius; /* Radius of an equivalent sphere */ - double a_over_c; /* Fixed aspect ratio */ + /* In use by SCHIFF_ELLIPSOID_AS_SPHERE */ + struct schiff_param radius_sphere; unsigned nslices; }; #define SCHIFF_ELLIPSOID_DEFAULT__ \ - {SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, 0, 64} + {SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, 64} static const struct schiff_ellipsoid SCHIFF_ELLIPSOID_DEFAULT = SCHIFF_ELLIPSOID_DEFAULT__; @@ -76,18 +78,22 @@ struct schiff_sphere { struct schiff_helical_pipe { struct schiff_param pitch; /* Elevation distance of a full revolution */ struct schiff_param height; /* Total heigh of the helical pipe */ - struct schiff_param helicoid_radius; /* Radius of the helicoid */ - struct schiff_param circle_radius; /* Radius of the meridian circle */ + struct schiff_param radius_helicoid; /* Radius of the helicoid */ + struct schiff_param radius_circle; /* Radius of the meridian circle */ + + /* In use by SCHIFF_HELICAL_PIPE_AS_SPHERE */ + struct schiff_param radius_sphere; unsigned nslices_helicoid; /* # Discrete steps of the helicoid */ unsigned nslices_circle; /* # Discrete steps along 2PI */ }; -#define SCHIFF_HELICAL_PIPE_DEFAULT__ \ +#define SCHIFF_HELICAL_PIPE_DEFAULT__ \ {SCHIFF_PARAM_DEFAULT__, \ SCHIFF_PARAM_DEFAULT__, \ SCHIFF_PARAM_DEFAULT__, \ SCHIFF_PARAM_DEFAULT__, \ + SCHIFF_PARAM_DEFAULT__, \ 128, 64} static const struct schiff_helical_pipe SCHIFF_HELICAL_PIPE_DEFAULT = SCHIFF_HELICAL_PIPE_DEFAULT__; @@ -98,13 +104,16 @@ static const struct schiff_sphere SCHIFF_SPHERE_DEFAULT = struct schiff_cylinder { struct schiff_param radius; - struct schiff_param height; /* In use by SCHIFF_CYLINDER */ - double aspect_ratio; /* In use by CYLINDER_AS_SPHERE */ + struct schiff_param height; + + /* In use by SCHIFF_CYLINDER_AS_SPHERE */ + struct schiff_param radius_sphere; + unsigned nslices; }; #define SCHIFF_CYLINDER_DEFAULT__ \ - {SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, -1.0, 64} + {SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, 64} static const struct schiff_cylinder SCHIFF_CYLINDER_DEFAULT = SCHIFF_CYLINDER_DEFAULT__; @@ -112,6 +121,8 @@ enum { A, B, M, N0, N1, N2 }; /* Super formula arguments */ struct schiff_supershape { struct schiff_param formulas[2][6]; + /* In use by SCHIFF_SUPERSHAPE_AS_SPHERE */ + struct schiff_param radius_sphere; unsigned nslices; }; #define SCHIFF_SUPERSHAPE_DEFAULT__ \ @@ -126,12 +137,15 @@ struct schiff_supershape { SCHIFF_PARAM_DEFAULT__, \ SCHIFF_PARAM_DEFAULT__, \ SCHIFF_PARAM_DEFAULT__, \ - SCHIFF_PARAM_DEFAULT__}}, 64} + SCHIFF_PARAM_DEFAULT__}}, \ + SCHIFF_PARAM_DEFAULT__, \ + 64} static const struct schiff_supershape SCHIFF_SUPERSHAPE_DEFAULT = SCHIFF_SUPERSHAPE_DEFAULT__; struct schiff_geometry { + /* Shape of the geometry */ enum schiff_geometry_type type; union { struct schiff_ellipsoid ellipsoid; @@ -142,17 +156,19 @@ struct schiff_geometry { } data; }; -#define SCHIFF_GEOMETRY_NULL__ { SCHIFF_NONE, { SCHIFF_ELLIPSOID_DEFAULT__ } } +#define SCHIFF_GEOMETRY_NULL__ {SCHIFF_NONE, { SCHIFF_ELLIPSOID_DEFAULT__ }} static const struct schiff_geometry SCHIFF_GEOMETRY_NULL = SCHIFF_GEOMETRY_NULL__; /* Forward declarations */ +struct s3d_device; struct schiff_optical_properties; struct sschiff_geometry_distribution; extern LOCAL_SYM res_T schiff_geometry_distribution_init (struct sschiff_geometry_distribution* distrib, + struct s3d_device* s3d, const struct schiff_geometry* geometry, const size_t ngeoms, const double characteristic_length, diff --git a/src/schiff_mesh.c b/src/schiff_mesh.c @@ -515,5 +515,6 @@ schiff_mesh_release(struct schiff_mesh* mesh) } darray_uint_release(&mesh->indices); mesh->coordinates = SCHIFF_NO_COORDINATE; + mesh->is_init = 0; }