schiff

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

commit d3148f07ab4939da6ba21a9bdf43e16b37099dae
parent 936a8b08f826d6d5f6313d6fe7de96596aa12de3
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 11 Mar 2016 16:29:15 +0100

Finalize the YAML parsing of geometry distributions

Use the updated Star-Schiff API. Remove outdated documentations of the
command help. Handle multiple geometry distributions.

Diffstat:
Msrc/schiff.c | 38++++++++++++++++++--------------------
Msrc/schiff_args.c | 349++++++++++++++++++++++++++++++++++++++++++++++++++-----------------------------
Msrc/schiff_args.h | 31+++++++++++++++++++++----------
Msrc/schiff_geometry.c | 303+++++++++++++++++++++++++++++++++++++++++++++----------------------------------
Msrc/schiff_geometry.h | 4++++
Msrc/schiff_mesh.h | 22+++++++++++-----------
6 files changed, 448 insertions(+), 299 deletions(-)

diff --git a/src/schiff.c b/src/schiff.c @@ -49,7 +49,6 @@ dump_geometries struct s3d_device* s3d = NULL; struct s3d_scene* scn = NULL; struct s3d_shape* shape = NULL; - struct sschiff_material mtl_dummy; unsigned i; res_T res = RES_OK; ASSERT(args && distrib && rng && stream); @@ -78,7 +77,7 @@ dump_geometries FOR_EACH(i, 0, args->ngeoms_dump) { unsigned ivert, nverts; unsigned itri, ntris; - res = distrib->sample(rng, &mtl_dummy, shape, distrib->context); + res = distrib->sample(rng, shape, distrib->context); if(res != RES_OK) { fprintf(stderr, "Couldn't sample the micro organism geometry.\n"); goto error; @@ -123,12 +122,11 @@ run_integration { struct sschiff_device* sschiff = NULL; struct sschiff_estimator* estimator = NULL; - struct sschiff_estimator_state* states = NULL; size_t iwlen, nwlens; res_T res = RES_OK; ASSERT(args && sa_size(args->wavelengths) && rng && stream); - res = sschiff_device_create(NULL, NULL, args->nthreads, 0, NULL, &sschiff); + res = sschiff_device_create(NULL, NULL, args->nthreads, 1, NULL, &sschiff); if(res != RES_OK) { fprintf(stderr, "Couldn't create the Star Schiff device.\n"); goto error; @@ -137,37 +135,36 @@ run_integration /* Invoke the Schiff integration */ nwlens = sa_size(args->wavelengths); - res = sschiff_integrate(sschiff, rng, distrib, args->wavelengths, nwlens, 1, - args->ngeoms, args->ndirs, &estimator); + res = sschiff_integrate(sschiff, rng, distrib, args->wavelengths, nwlens, + sschiff_uniform_scattering_angles, args->nangles, args->nrealisations, + args->ndirs, &estimator); if(res != RES_OK) { fprintf(stderr, "Schiff integration error.\n"); goto error; } - /* Retrieve estimation results */ - states = sa_add(states, nwlens); - SSCHIFF(estimator_get_states(estimator, states)); - - /* Print the estimation results */ + /* Print the estimated cross sections */ FOR_EACH(iwlen, 0, nwlens) { - const struct sschiff_estimator_value* val; + const struct sschiff_state* val; + struct sschiff_cross_section cross_section; + + SSCHIFF(estimator_get_cross_section(estimator, iwlen, &cross_section)); - fprintf(stream, "%g ", states[iwlen].wavelength); + fprintf(stream, "%g ", args->wavelengths[iwlen]); - val = states[iwlen].values + SSCHIFF_EXTINCTION_CROSS_SECTION; + val = &cross_section.extinction; fprintf(stream, "%g %g ", val->E, val->SE); - val = states[iwlen].values + SSCHIFF_ABSORPTION_CROSS_SECTION; + val = &cross_section.absorption; fprintf(stream, "%g %g ", val->E, val->SE); - val = states[iwlen].values + SSCHIFF_SCATTERING_CROSS_SECTION; + val = &cross_section.scattering; fprintf(stream, "%g %g ", val->E, val->SE); - val = states[iwlen].values + SSCHIFF_AVERAGE_PROJECTED_AREA; + val = &cross_section.average_projected_area; fprintf(stream, "%g %g\n", val->E, val->SE); } exit: if(estimator) SSCHIFF(estimator_ref_put(estimator)); if(sschiff) SSCHIFF(device_ref_put(sschiff)); - sa_release(states); return res; error: goto exit; @@ -198,8 +195,9 @@ run(const struct schiff_args* args) goto error; } - res = schiff_geometry_distribution_init - (&distrib, &args->geom, args->properties); + res = schiff_geometry_distribution_init(&distrib, args->geoms, + sa_size(args->geoms), args->caracteristic_length, args->ran_geoms, + args->properties); if(res != RES_OK) { fprintf(stderr, "Couldn't initialise the Star Schiff geometry distribution.\n"); diff --git a/src/schiff_args.c b/src/schiff_args.c @@ -54,129 +54,94 @@ print_help(const char* binary) { printf( "Usage: %s [OPTIONS] [FILE]\n" -"Estimate the radiative properties of soft particles with an \"Approximation\n" -"Method for Short Wavelength or High Energy Scattering\" (L. Schiff, 1956).\n" +"Estimate the radiative properties of soft particles with an \"Approximation\n" +"Method for Short Wavelength or High Energy Scattering\" (L. Schiff, 1956).\n" "Assume a small contrast on the relative refractive index of the particles.\n\n", binary); printf( -"FILE lists the per wavelength optical properties of the soft particles. Each\n" -"line must be formatted as \"W Nr Kr Ne\" where \"W\" is the wavelength in vacuum\n" -"expressed in micron, \"Nr\" and \"Kr\" are the real and imaginary parts,\n" -"respectively, of the relative refractive index, and \"Ne\" the refractive index\n" -"of the medium. With no FILE, read optical properties from standard input.\n\n"); +"FILE lists the per wavelength optical properties of the soft particles. Each\n" +"line must be formatted as \"W Nr Kr Ne\" where \"W\" is the wavelength in\n" +"vacuum expressed in micron, \"Nr\" and \"Kr\" are the real and imaginary parts,\n" +"respectively, of the relative refractive index, and \"Ne\" the refractive index\n" +"of the medium. With no FILE, read optical properties from standard input.\n\n"); printf( -"The estimated results are written to OUTPUT or to standard output whether the\n" -"\"-o\" option is defined or not, respectively. The line format of the outputed\n" -"results is \"W E E' A A' S S' P P'\" with \"W\" the wavelength in vacuum (expressed\n" -"in micron), \"E\", \"A\" and \"S\" the estimation of the extinction, absorption and\n" -"scattering cross section, respectively, and \"P\" the estimated average projected\n" -"area. The \"E'\", \"A'\", \"S'\" and \"P'\" values are the standard error of the\n" -"aforementioned estimations.\n\n"); - printf( -"The shapes of the soft particles are controlled by a geometric distribution\n" -"whose parameters are distributed according to an user defined distribution. The\n" -"available distributions of the parameters and the syntax to control them are\n" -"described below:\n\n"); - printf( -" c=VAL parameter is the constant VAL.\n"); - printf( -" h=HISTOGRAM The first valid line of the HISTOGRAM file should be formatted\n" -" as \"LOWER UPPER NINTERVALS\" where \"LOWER\" and \"UPPER\" are,\n" -" respectively, the lower and the upper bounds of the histogram,\n" -" and \"NINTERVALS\" the number of intervals. The remaining valid\n" -" lines should list the probability of the interval upper bound\n" -" until all intervals are defined.\n"); - printf( -" l=ZETA:SIGMA parameter is distributed with respect to the lognormal\n" -" distribution:\n" -" P(x) dx = 1/(log(SIGMA)*x*sqrt(2*PI))\n" -" * exp(-(ln(x)-log(ZETA))^2 / (2*log(SIGMA)^2)) dx\n\n"); +"The estimated results are written to OUTPUT or to standard output whether the\n" +"\"-o\" option is defined or not, respectively. The line format of the outputed\n" +"results is \"W E E' A A' S S' P P'\" with \"W\" the wavelength in vacuum\n" +"(expressed in micron), \"E\", \"A\" and \"S\" the estimation of the extinction,\n" +"absorption and scattering cross section, respectively, and \"P\" the estimated\n" +"average projected area. The \"E'\", \"A'\", \"S'\" and \"P'\" values are the\n" +"standard error of the aforementioned estimations.\n\n"); printf( "Options:\n\n"); printf( -" -c RDISTRIB,HDISTRIB[,N]\n" -" the soft particles are cylindrical meshes discretized in N\n" -" slices. By default N is %u. RDISTRIB and HDISTRIB are the\n" -" distributions of, respectively, the radius and the height\n" -" parameters of the cylinder.\n", - SCHIFF_CYLINDER_DEFAULT.nslices); - printf( -" -C ASPECT_RATIO,RDISTRIB[,N]\n" -" the soft particles are cylindrical meshes discretized in N\n" -" slices. By default N is %u. The volume of the cylinders is\n" -" defined according to the volume of a sphere whose radius \"R\" is\n" -" distributed with respect to the RDISTRIB distribution. The\n" -" ratio between the cylinder radius \"r\" and its height \"h\" is\n" -" fixed by the ASPECT_RATIO floating point parameter.\n", - SCHIFF_CYLINDER_DEFAULT.nslices); - printf( -" r = R * ((2*ASPECT_RATIO)/3)^1/3\n" -" h = r * 2 / ASPECT_RATIO\n"); - +" -a ANGLES number of phase function scattering angles. Default is %u.\n", + SCHIFF_ARGS_NULL.nangles); printf( -" -d DIRS number of sampled directions for each geometry. Default is %u.\n", +" -d DIRS number of sampled directions for each geometry.\n" +" Default is %u.\n", SCHIFF_ARGS_NULL.ndirs); printf( -" -g GOEMS number of sampled geometries. This is actually the number of\n" -" realisations. Default is %u.\n", - SCHIFF_ARGS_NULL.ngeoms); +" -g GOEMS number of sampled geometries. This is actually the number of\n" +" realisations. Default is %u.\n", + SCHIFF_ARGS_NULL.nrealisations); printf( -" -G NUM sampled `NUM' geometries with respect to the defined\n" -" distribution, dump their data and exit. The data are written\n" -" to OUTPUT or the standard output whether the `-o' option is\n" -" defined or not, respectively. The outputted data followed the\n" +" -G NUM sampled `NUM' geometries with respect to the defined\n" +" distribution, dump their data and exit. The data are written to\n" +" OUTPUT or the standard output whether the `-o' option is\n" +" defined or not, respectively. The outputted data followed the\n" " Alias Wavefront obj file format.\n"); printf( " -h display this help and exit.\n"); printf( -" -n NTHREADS hint on the number of threads to use during the integration.\n" +" -i DISTRIB YAML file that defines the geometry distributions of the soft\n" +" particles.\n"); + printf( +" -l LENGTH caracteristic length of the soft particles.\n"); + printf( +" -n NTHREADS hint on the number of threads to use during the integration.\n" " By default use as many threads as CPU cores.\n"); printf( -" -o OUTPUT write results to OUTPUT. If not defined, write results to\n" +" -o OUTPUT write results to OUTPUT. If not defined, write results to\n" " standard output.\n"); printf( " -q do not print the helper message when no FILE is submitted.\n"); printf( -" -s RDISTRIB[,N]\n" -" the soft particles are spherical meshes discretized in N slices\n" -" along 2PI. By default N is %u. Their radius parameter is\n" -" distributed with respect to the RDISTRIB distribution.\n", - SCHIFF_SPHERE_DEFAULT.nslices); - printf( -" -u FORMULA0#FORMULA1[#N]\n" -" the soft particles are 3D super shapes discretized in N slices\n" -" along 2PI. By default N is %u. The FORMULA<0|1> arguments\n" -" control the parameter distribution of the superformulas. Each\n" -" FORMULA is formated as \"A,B,M,N0,N1,N2\" where A, B, M, N0, N1\n" -" and N2 are the parameter distributions of the corresponding\n" -" superformula parmater.\n", - SCHIFF_SUPER_SHAPE_DEFAULT.nslices); - printf( -" -w A[:B]... list of wavelengths in vacuum (expressed in micron) to\n" +" -w A[:B]... list of wavelengths in vacuum (expressed in micron) to\n" " integrate.\n\n"); - printf( -"Examples:\n\n"); - printf( -" Spheric soft particles whose radius is distributed according to a lognormal\n" -" distribution and their optical properties are defined in the \"optical_props\"\n" -" file. The estimation is evaluated for the wavelength 0.6 micron:\n" -" %s -s l=1:1.18 -w 0.6 optical_props\n\n", - binary); - printf( -" Soft particles are cylinders discretized in 128 slices. Their radius is a\n" -" constant and their height follows the distribution described in the \"histo\"\n" -" file. Optical properties are read from standard input. Estimation is\n" -" performed for the wavelengths 0.5, 0.6 and 0.65 micron:\n" -" %s -c c=4.1e-3,h=hist,128 -w 0.5:0.6:0.65\n\n", - binary); - printf( -" Cylindrical soft particles whith a fixed aspect ratio of 1.33. Their volume\n" -" is defined with respect to a sphere whose radius is distributed according to\n" -" a lognormal distribution. Optical properties are read from standard input.\n" -" The estimation is evaluated for the wavelength 0.6 micron and is stored in\n" -" the \"output\" file:\n" -" %s -C 1.33,l=1.2:1.834 -w 0.6 -o output\n", - binary); +} + +static INLINE void +param_release(struct schiff_param* param) +{ + ASSERT(param); + + if(param->distribution == SCHIFF_PARAM_HISTOGRAM + && param->data.histogram.entries) { + sa_release(param->data.histogram.entries); + } + param->distribution = SCHIFF_PARAM_NONE; +} + +static void +geometry_release(struct schiff_geometry* geom) +{ + ASSERT(geom); + switch(geom->type) { + 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); + break; + case SCHIFF_SPHERE: + param_release(&geom->data.sphere.radius); + break; + case SCHIFF_NONE: /* Do nothing */ break; + default: FATAL("Unreachable code\n"); break; + } + geom->type = SCHIFF_NONE; } static res_T @@ -252,6 +217,26 @@ parse_yaml_double(const char* filename, const yaml_node_t* node, double* value) } static res_T +parse_yaml_uint(const char* filename, const yaml_node_t* node, unsigned* value) +{ + res_T res = RES_OK; + ASSERT(filename && node && value); + + if(node->type != YAML_SCALAR_NODE) { + log_err(filename, node, "expecting an unsigned integer.\n"); + return RES_BAD_ARG; + } + res = cstr_to_uint((char*)node->data.scalar.value, value); + if(res != RES_OK) { + log_err(filename, node, "invalid unsigned integer `%s'.\n", + node->data.scalar.value); + return RES_BAD_ARG; + } + + return RES_OK; +} + +static res_T parse_yaml_param_lognormal (const char* filename, yaml_document_t* doc, @@ -527,19 +512,21 @@ parse_yaml_cylinder (const char* filename, yaml_document_t* doc, const yaml_node_t* node, - struct schiff_geometry* geom) + struct schiff_geometry* geom, + double* geom_proba) { enum { PROBA = BIT(0), RADIUS = BIT(1), HEIGHT = BIT(2), - ASPECT_RATIO = BIT(3) + ASPECT_RATIO = BIT(3), + SLICES = BIT(4) }; int mask = 0; /* Register the parsed attributes */ size_t nattrs; size_t i; res_T res = RES_OK; - ASSERT(filename && doc && node && geom); + ASSERT(filename && doc && node && geom && geom_proba); if(node->type != YAML_MAPPING_NODE) { log_err(filename, node, "expecting a mapping of cylinder attributes.\n"); @@ -559,13 +546,12 @@ parse_yaml_cylinder /* Distribution probability */ if(!strcmp((char*)key->data.scalar.value, "proba")) { - double proba; /* TODO use this proba */ if(mask & PROBA) { log_err(filename, key, "the cylinder proba is already defined.\n"); return RES_BAD_ARG; } mask |= PROBA; - res = parse_yaml_double(filename, val, &proba); + res = parse_yaml_double(filename, val, geom_proba); /* Cylinder radius */ } else if(!strcmp((char*)key->data.scalar.value, "radius")) { @@ -594,6 +580,14 @@ parse_yaml_cylinder mask |= ASPECT_RATIO; res = parse_yaml_double (filename, val, &geom->data.cylinder.aspect_ratio); + } else if(!strcmp((char*)key->data.scalar.value, "slices")) { + if(mask & SLICES) { + log_err(filename, key, + "the cylinder number of slices is already defined.\n"); + return RES_BAD_ARG; + } + mask |= SLICES; + res = parse_yaml_uint(filename, val, &geom->data.cylinder.nslices); } else { log_err(filename, key, "unkown cylinder attribute `%s'.\n", key->data.scalar.value); @@ -606,9 +600,6 @@ parse_yaml_cylinder if(!(mask & RADIUS)) { log_err(filename, node, "missing the radius attribute.\n"); return RES_BAD_ARG; - } else if(!(mask & PROBA)) { - log_err(filename, node, "missing the proba 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; @@ -619,6 +610,13 @@ parse_yaml_cylinder return RES_BAD_ARG; } + 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) { geom->type = SCHIFF_CYLINDER; } else { @@ -633,17 +631,19 @@ parse_yaml_sphere (const char* filename, yaml_document_t* doc, const yaml_node_t* node, - struct schiff_geometry* geom) + struct schiff_geometry* geom, + double* geom_proba) { enum { PROBA = BIT(0), - RADIUS = BIT(1) + RADIUS = BIT(1), + SLICES = BIT(2) }; int mask = 0; /* Register the parsed attributes */ size_t nattrs; size_t i; res_T res = RES_OK; - ASSERT(filename && doc && node && geom); + ASSERT(filename && doc && node && geom && geom_proba); if(node->type != YAML_MAPPING_NODE) { log_err(filename, node, "expecting a mapping of sphere attributes.\n"); @@ -662,13 +662,12 @@ parse_yaml_sphere ASSERT(key->type == YAML_SCALAR_NODE); if(!strcmp((char*)key->data.scalar.value, "proba")) { - double proba; /* TODO use this proba */ if(mask & PROBA) { log_err(filename, node, "the sphere proba is already defined.\n"); return RES_BAD_ARG; } mask |= PROBA; - res = parse_yaml_double(filename, val, &proba); + res = parse_yaml_double(filename, val, geom_proba); } else if(!strcmp((char*)key->data.scalar.value, "radius")) { if(mask & RADIUS) { log_err(filename, key, "the sphere radius is already defined.\n"); @@ -677,6 +676,14 @@ parse_yaml_sphere mask |= RADIUS; res = parse_yaml_param_distribution (filename, doc, val, &geom->data.sphere.radius); + } else if(!strcmp((char*)key->data.scalar.value, "slices")) { + if(mask & SLICES) { + log_err(filename, key, + "the sphere number of slices is already defined.\n"); + return RES_BAD_ARG; + } + mask |= SLICES; + res = parse_yaml_uint(filename, val, &geom->data.sphere.nslices); } else { log_err(filename, key, "unkown sphere attribute `%s'.\n", key->data.scalar.value); @@ -689,26 +696,35 @@ parse_yaml_sphere if(!(mask & RADIUS)) { log_err(filename, node, "missing the radius attribute.\n"); return RES_BAD_ARG; - } else if(!(mask & PROBA)) { - log_err(filename, node, "missing the proba attribute.\n"); - return RES_BAD_ARG; } - + if(!(mask & PROBA)) { /* Default proba */ + *geom_proba = 1.0; + } + if(!(mask & SLICES)) { /* Default number of slices */ + geom->data.sphere.nslices = SCHIFF_SPHERE_DEFAULT.nslices; + } + geom->type = SCHIFF_SPHERE; return RES_OK; } static res_T -parse_yaml(const char* filename) +parse_yaml + (const char* filename, + struct schiff_geometry** out_geoms, + struct ssp_ran_discrete** out_ran) { - struct schiff_geometry geom; yaml_parser_t parser; yaml_document_t doc; yaml_node_t* root; size_t ndistribs; size_t idistrib; + struct schiff_geometry* geoms = NULL; + double* probas = NULL; + double accum_proba = 0; + struct ssp_ran_discrete* ran = NULL; FILE* file = NULL; res_T res = RES_OK; - ASSERT(filename); + ASSERT(filename && out_geoms && out_ran); if(!yaml_parser_initialize(&parser)) { fprintf(stderr, "Couldn't intialise the YAML parser.\n"); @@ -744,6 +760,32 @@ parse_yaml(const char* filename) ndistribs = (size_t) (root->data.sequence.items.top - root->data.sequence.items.start); + + if(!sa_add(geoms, ndistribs)) { + log_err(filename, root, + "couldn't allocate up to %lu geometry distributions.\n", + (unsigned long)ndistribs); + res = RES_MEM_ERR; + goto error; + } + FOR_EACH(idistrib, 0, ndistribs) + geoms[idistrib] = SCHIFF_GEOMETRY_NULL; + + if(!sa_add(probas, ndistribs)) { + log_err(filename, root, + "couldn't allocate the list of geometry distribution probabilities.\n"); + res = RES_MEM_ERR; + goto error; + } + + res = ssp_ran_discrete_create(&mem_default_allocator, &ran); + if(res != RES_OK) { + log_err(filename, root, + "couldn't allocate the random variate of geometry distributions.\n"); + goto error; + } + + accum_proba = 0; FOR_EACH(idistrib, 0, ndistribs) { yaml_node_t* distrib, *key, *val; @@ -762,9 +804,11 @@ parse_yaml(const char* filename) ASSERT(key->type == YAML_SCALAR_NODE); if(!strcmp((char*)key->data.scalar.value, "cylinder")) { - res = parse_yaml_cylinder(filename, &doc, val, &geom); + res = parse_yaml_cylinder + (filename, &doc, val, &geoms[idistrib], &probas[idistrib]); } else if(!strcmp((char*)key->data.scalar.value, "sphere")) { - res = parse_yaml_sphere(filename, &doc, val, &geom); + res = parse_yaml_sphere + (filename, &doc, val, &geoms[idistrib], &probas[idistrib]); } else { log_err(filename, key, "unknown distribution `%s'.\n", key->data.scalar.value); @@ -772,14 +816,40 @@ parse_yaml(const char* filename) goto error; } if(res != RES_OK) goto error; + accum_proba += probas[idistrib]; + } + + FOR_EACH(idistrib, 0, ndistribs) { + probas[idistrib] /= accum_proba; + } + + res = ssp_ran_discrete_setup(ran, probas, ndistribs); + if(res != RES_OK) { + log_err(filename, root, + "couldn't setup the discrete geometry distributions.\n"); + goto error; } exit: yaml_parser_delete(&parser); yaml_document_delete(&doc); if(file) fclose(file); + if(probas) sa_release(probas); + *out_geoms = geoms; + *out_ran = ran; return res; error: + if(ran) { + SSP(ran_discrete_ref_put(ran)); + ran = NULL; + } + if(geoms) { + FOR_EACH(idistrib, 0, ndistribs) { + geometry_release(&geoms[idistrib]); + } + sa_release(geoms); + geoms = NULL; + } goto exit; } @@ -798,17 +868,23 @@ schiff_args_init ASSERT(argc && argv && args); *args = SCHIFF_ARGS_NULL; - while((opt = getopt(argc, argv, "d:g:G:hi:n:o:qw:")) != -1) { + while((opt = getopt(argc, argv, "a:d:g:G:hi:l:n:o:qw:")) != -1) { switch(opt) { + case 'a': res = cstr_to_uint(optarg, &args->nangles); break; case 'd': res = cstr_to_uint(optarg, &args->ndirs); break; - case 'g': res = cstr_to_uint(optarg, &args->ngeoms); break; + case 'g': res = cstr_to_uint(optarg, &args->nrealisations); break; case 'G': res = cstr_to_uint(optarg, &args->ngeoms_dump); break; case 'h': print_help(argv[0]); schiff_args_release(args); return RES_OK; case 'i': - res = parse_yaml(optarg); + res = parse_yaml(optarg, &args->geoms, &args->ran_geoms); + break; + case 'l': + res = cstr_to_double(optarg, &args->caracteristic_length); + if(res == RES_OK && args->caracteristic_length <= 0.0) + res = RES_BAD_ARG; break; case 'n': res = cstr_to_uint(optarg, &args->nthreads); @@ -829,9 +905,9 @@ schiff_args_init } } - if(args->geom.type == SCHIFF_NONE) { + if(args->geoms == NULL) { fprintf(stderr, - "%s: missing geometry distribution.\nTry '%s -h' for more information.\n", + "%s: missing geometry distribution.\nTry '%s -h' for more informations.\n", argv[0], argv[0]); res = RES_BAD_ARG; goto error; @@ -841,7 +917,15 @@ schiff_args_init if(!args->wavelengths) { fprintf(stderr, - "Usage: %s [OPTIONS] [FILE]\nTry '%s -h' for more information.\n", + "%s: missing wavelengths.\nTry '%s -h' for more informations.\n", + argv[0], argv[0]); + res = RES_BAD_ARG; + goto error; + } + + if(args->caracteristic_length < 0.0) { + fprintf(stderr, +"%s: missing the caracteristic length.\nTry '%s -h' for more informations\n", argv[0], argv[0]); res = RES_BAD_ARG; goto error; @@ -883,9 +967,16 @@ error: void schiff_args_release(struct schiff_args* args) { + size_t i, count; ASSERT(args); sa_release(args->properties); sa_release(args->wavelengths); + + count = sa_size(args->geoms); + FOR_EACH(i, 0, count) geometry_release(&args->geoms[i]); + sa_release(args->geoms); + if(args->ran_geoms) SSP(ran_discrete_ref_put(args->ran_geoms)); + args->geoms = NULL; *args = SCHIFF_ARGS_NULL; } diff --git a/src/schiff_args.h b/src/schiff_args.h @@ -34,24 +34,35 @@ #include <star/sschiff.h> struct schiff_args { - const char* output_filename; - struct schiff_optical_properties* properties; - double* wavelengths; - struct schiff_geometry geom; - unsigned ngeoms_dump; - unsigned ngeoms; - unsigned ndirs; - unsigned nthreads; + const char* output_filename; /* File in which output data are stored */ + struct schiff_optical_properties* properties; /* List of properties */ + double* wavelengths; /* List of wavelengths to integrate */ + + double caracteristic_length; + + /* List of geometry and its associated random variates */ + struct schiff_geometry* geoms; + struct ssp_ran_discrete* ran_geoms; + + unsigned ngeoms_dump; /* # geometries to dump */ + unsigned nrealisations; /* # realisation */ + unsigned ndirs; /* Number of directions to sample per realisation */ + unsigned nangles; /* Number of scattering angles */ + + unsigned nthreads; /* Hint on the number of thread to use */ }; static const struct schiff_args SCHIFF_ARGS_NULL = { NULL, /* Output filename */ NULL, /* List of optical properties */ NULL, /* List of wavelength to integrate */ - SCHIFF_GEOMETRY_NULL__, /* List of Schiff geometries */ + -1.0, /* Caracteristic length */ + NULL, /* List of Schiff geometries */ + NULL, /* Schiff geometry random variates */ 0, /* # Dumped geometries */ 1000, /* # Sampled geometries */ - 100, /* # Sampled directions per gemetry */ + 100, /* # Sampled directions per geometry */ + 1000, /* # Scattering angles */ SSCHIFF_NTHREADS_DEFAULT }; diff --git a/src/schiff_geometry.c b/src/schiff_geometry.c @@ -61,9 +61,9 @@ struct cylinder_distribution_context { }; /* 3D context of a generic geometry */ -struct geometry_context { - struct schiff_mesh* mesh; /* Triangular mesh of the geometry */ - struct schiff_geometry* geometry; +struct mesh_context { + const struct schiff_mesh* mesh; /* Triangular mesh of the geometry */ + const struct schiff_geometry* geometry; enum schiff_geometry_type type; union { struct cylinder_context cylinder; @@ -74,9 +74,10 @@ struct geometry_context { /* Distribution context of a generic geometry */ struct geometry_distribution_context { - struct schiff_mesh* mesh; /* Triangular mesh of the geometry */ struct schiff_optical_properties* properties; /* Per wavelength properties */ - const struct schiff_geometry* geometry; /* Geometry descriptor */ + struct schiff_mesh* meshes; /* List of triangular meshes */ + struct schiff_geometry const ** geometries; /* List of geometries */ + struct ssp_ran_discrete* ran_geometries; /* Geometries random variates */ }; /******************************************************************************* @@ -164,47 +165,47 @@ get_material_property static void geometry_get_indices(const unsigned itri, unsigned ids[3], void* ctx) { - struct geometry_context* geom = ctx; - ASSERT(geom); - schiff_mesh_get_indices(geom->mesh, itri, ids); + struct mesh_context* mesh_ctx = ctx; + ASSERT(ctx); + schiff_mesh_get_indices(mesh_ctx->mesh, itri, ids); } static void cylinder_get_position(const unsigned ivert, float vertex[3], void* ctx) { - struct geometry_context* geom = ctx; - ASSERT(geom && geom->type == SCHIFF_CYLINDER); - schiff_mesh_get_cartesian_position(geom->mesh, ivert, vertex); - vertex[0] *= geom->data.cylinder.radius; - vertex[1] *= geom->data.cylinder.radius; - vertex[2] *= geom->data.cylinder.height; + struct mesh_context* mesh_ctx = ctx; + ASSERT(mesh_ctx && mesh_ctx->type == SCHIFF_CYLINDER); + schiff_mesh_get_cartesian_position(mesh_ctx->mesh, ivert, vertex); + vertex[0] *= mesh_ctx->data.cylinder.radius; + vertex[1] *= mesh_ctx->data.cylinder.radius; + vertex[2] *= mesh_ctx->data.cylinder.height; } static void sphere_get_position(const unsigned ivert, float vertex[3], void* ctx) { - struct geometry_context* geom = ctx; - ASSERT(geom && geom->type == SCHIFF_SPHERE); - schiff_mesh_get_cartesian_position(geom->mesh, ivert, vertex); - vertex[0] *= geom->data.sphere.radius; - vertex[1] *= geom->data.sphere.radius; - vertex[2] *= geom->data.sphere.radius; + struct mesh_context* mesh_ctx = ctx; + ASSERT(mesh_ctx && mesh_ctx->type == SCHIFF_SPHERE); + schiff_mesh_get_cartesian_position(mesh_ctx->mesh, ivert, vertex); + vertex[0] *= mesh_ctx->data.sphere.radius; + vertex[1] *= mesh_ctx->data.sphere.radius; + vertex[2] *= mesh_ctx->data.sphere.radius; } static void super_shape_get_position(const unsigned ivert, float vertex[3], void* ctx) { - struct geometry_context* geom = ctx; + struct mesh_context* mesh_ctx = ctx; struct sin_cos angles[2]; double uv[2]; int iform; - ASSERT(geom && geom->type == SCHIFF_SUPER_SHAPE); + ASSERT(mesh_ctx && mesh_ctx->type == SCHIFF_SUPER_SHAPE); - schiff_mesh_get_polar_position(geom->mesh, ivert, angles); + schiff_mesh_get_polar_position(mesh_ctx->mesh, ivert, angles); FOR_EACH(iform, 0, 2) { double m, k, g; - double* form = geom->data.super_shape.formulas[iform]; + double* form = mesh_ctx->data.super_shape.formulas[iform]; m = cos(form[M] * angles[iform].angle / 4.0); m = fabs(m) / form[A]; @@ -222,32 +223,31 @@ super_shape_get_position(const unsigned ivert, float vertex[3], void* ctx) static res_T geometry_sample_cylinder (struct ssp_rng* rng, - struct sschiff_material* mtl, - struct s3d_shape* shape, - void* context) + const struct schiff_geometry* geom, + const struct schiff_mesh* mesh, + struct s3d_shape* shape) { - struct geometry_distribution_context* distrib = context; const struct schiff_cylinder* cylinder; - struct geometry_context geom; + struct mesh_context mesh_ctx; struct s3d_vertex_data attrib; size_t nverts, nprims; double r; - ASSERT(rng && mtl && shape && context); + ASSERT(rng && geom && mesh && shape); - cylinder = &distrib->geometry->data.cylinder; - geom.type = SCHIFF_CYLINDER; - geom.mesh = distrib->mesh; - switch(distrib->geometry->type) { + cylinder = &geom->data.cylinder; + mesh_ctx.type = SCHIFF_CYLINDER; + mesh_ctx.mesh = mesh; + switch(geom->type) { case SCHIFF_CYLINDER: - geom.data.cylinder.radius = (float)eval_param(&cylinder->radius, rng); - geom.data.cylinder.height = (float)eval_param(&cylinder->height, rng); + 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); - geom.data.cylinder.radius = + mesh_ctx.data.cylinder.radius = (float)(r / pow(3.0 / (2.0*cylinder->aspect_ratio), 1.0/3.0)); - geom.data.cylinder.height = - (float)(2.f * geom.data.cylinder.radius / cylinder->aspect_ratio); + mesh_ctx.data.cylinder.height = + (float)(2.f * mesh_ctx.data.cylinder.radius / cylinder->aspect_ratio); break; default: FATAL("Unreachable code.\n"); break; } @@ -256,73 +256,65 @@ geometry_sample_cylinder attrib.type = S3D_FLOAT3; attrib.get = cylinder_get_position; - mtl->get_property = get_material_property; - mtl->material = distrib->properties; - - nverts = darray_float_size_get(&distrib->mesh->vertices.cartesian) / 3/*#coords*/; - nprims = darray_uint_size_get(&distrib->mesh->indices) / 3/*#indices per prim*/; + nverts = darray_float_size_get(&mesh->vertices.cartesian) / 3/*#coords*/; + nprims = darray_uint_size_get(&mesh->indices) / 3/*#indices per prim*/; return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, - geometry_get_indices, (unsigned)nverts, &attrib, 1, &geom); + geometry_get_indices, (unsigned)nverts, &attrib, 1, &mesh_ctx); } static res_T geometry_sample_sphere (struct ssp_rng* rng, - struct sschiff_material* mtl, - struct s3d_shape* shape, - void* context) + const struct schiff_geometry* geom, + const struct schiff_mesh* mesh, + struct s3d_shape* shape) { - struct geometry_distribution_context* distrib = context; const struct schiff_sphere* sphere; - struct geometry_context geom; + struct mesh_context mesh_ctx; struct s3d_vertex_data attrib; size_t nverts, nprims; - ASSERT(rng && mtl && shape && context); - ASSERT(distrib->geometry->type == SCHIFF_SPHERE); + ASSERT(rng && geom && mesh && shape); + ASSERT(geom->type == SCHIFF_SPHERE); - sphere = &distrib->geometry->data.sphere; - geom.mesh = distrib->mesh; - geom.type = SCHIFF_SPHERE; - geom.data.sphere.radius = (float)eval_param(&sphere->radius, rng); + sphere = &geom->data.sphere; + mesh_ctx.mesh = mesh; + mesh_ctx.type = SCHIFF_SPHERE; + mesh_ctx.data.sphere.radius = (float)eval_param(&sphere->radius, rng); attrib.usage = S3D_POSITION; attrib.type = S3D_FLOAT3; attrib.get = sphere_get_position; - mtl->get_property = get_material_property; - mtl->material = distrib->properties; - - nverts = darray_float_size_get(&distrib->mesh->vertices.cartesian) / 3/*#coords*/; - nprims = darray_uint_size_get(&distrib->mesh->indices) / 3/*#indices*/; + nverts = darray_float_size_get(&mesh->vertices.cartesian) / 3/*#coords*/; + nprims = darray_uint_size_get(&mesh->indices) / 3/*#indices*/; return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, - geometry_get_indices, (unsigned)nverts, &attrib, 1, &geom); + geometry_get_indices, (unsigned)nverts, &attrib, 1, &mesh_ctx); } static res_T geometry_sample_super_shape (struct ssp_rng* rng, - struct sschiff_material* mtl, - struct s3d_shape* shape, - void* context) + const struct schiff_geometry* geom, + const struct schiff_mesh* mesh, + struct s3d_shape* shape) { - struct geometry_distribution_context* distrib = context; - struct geometry_context geom; - struct s3d_vertex_data attrib; const struct schiff_super_shape* sshape; + struct mesh_context mesh_ctx; + struct s3d_vertex_data attrib; size_t nverts, nprims; int iform, iattr; - ASSERT(rng && mtl && shape && context); - ASSERT(distrib->geometry->type == SCHIFF_SUPER_SHAPE); + ASSERT(rng && geom && mesh && shape); + ASSERT(geom->type == SCHIFF_SUPER_SHAPE); - sshape = &distrib->geometry->data.super_shape; - geom.mesh = distrib->mesh; - geom.type = SCHIFF_SUPER_SHAPE; + sshape = &geom->data.super_shape; + mesh_ctx.mesh = mesh; + mesh_ctx.type = SCHIFF_SUPER_SHAPE; FOR_EACH(iform, 0, 2) { FOR_EACH(iattr, 0, 6) { - geom.data.super_shape.formulas[iform][iattr] = + mesh_ctx.data.super_shape.formulas[iform][iattr] = (float)eval_param(&sshape->formulas[iform][iattr], rng); }} @@ -330,95 +322,148 @@ geometry_sample_super_shape attrib.type = S3D_FLOAT3; attrib.get = super_shape_get_position; - mtl->get_property = get_material_property; - mtl->material = distrib->properties; - - nverts = darray_uint_size_get(&distrib->mesh->vertices.polar) / 2/*#theta/phi*/; - nprims = darray_uint_size_get(&distrib->mesh->indices) / 3/*#indices*/; + nverts = darray_uint_size_get(&mesh->vertices.polar) / 2/*#theta/phi*/; + nprims = darray_uint_size_get(&mesh->indices) / 3/*#indices*/; return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, - geometry_get_indices, (unsigned)nverts, &attrib, 1, &geom); + geometry_get_indices, (unsigned)nverts, &attrib, 1, &mesh_ctx); +} + +static res_T +geometry_sample + (struct ssp_rng* rng, + struct s3d_shape* shape, + void* ctx) +{ + struct geometry_distribution_context* distrib = ctx; + const struct schiff_geometry* geom; + const struct schiff_mesh* mesh; + size_t isamp; + res_T res = RES_OK; + ASSERT(rng && shape && ctx); + + isamp = ssp_ran_discrete(rng, distrib->ran_geometries); + geom = distrib->geometries[isamp]; + mesh = distrib->meshes + isamp; + + switch(geom->type) { + case SCHIFF_CYLINDER: + case SCHIFF_CYLINDER_AS_SPHERE: + res = geometry_sample_cylinder(rng, geom, mesh, shape); + break; + case SCHIFF_SPHERE: + res = geometry_sample_sphere(rng, geom, mesh, shape); + break; + case SCHIFF_SUPER_SHAPE: + res = geometry_sample_super_shape(rng, geom, mesh, shape); + break; + default: FATAL("Unreachable code.\n"); break; + } + if(res != RES_OK) goto error; + +exit: + return res; +error: + goto exit; } /******************************************************************************* * Local functions ******************************************************************************/ +void +schiff_geometry_distribution_release + (struct sschiff_geometry_distribution* distrib) +{ + struct geometry_distribution_context* ctx; + size_t i, count; + ASSERT(distrib); + ctx = distrib->context; + + count = sa_size(ctx->meshes); + FOR_EACH(i, 0, count) schiff_mesh_release(&ctx->meshes[i]); + sa_release(ctx->meshes); + sa_release(ctx->geometries); + SSP(ran_discrete_ref_put(ctx->ran_geometries)); + mem_rm(ctx); +} + res_T schiff_geometry_distribution_init (struct sschiff_geometry_distribution* distrib, - const struct schiff_geometry* geom, + const struct schiff_geometry* geoms, + const size_t ngeoms, + const double caracteristic_length, + struct ssp_ran_discrete* ran_geoms, struct schiff_optical_properties* properties) { struct geometry_distribution_context* ctx = NULL; + size_t i; res_T res = RES_OK; - ASSERT(distrib && geom); + ASSERT(distrib && geoms && ngeoms && ran_geoms && properties); + + *distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION; ctx = mem_calloc(1, sizeof(struct geometry_distribution_context)); if(!ctx) { fprintf(stderr, -"Not enough memory: couldn't allocate the geometry distribution context\n"); + "Couldn't allocate the geometry distribution context\n"); res = RES_MEM_ERR; goto error; } - ctx->mesh = mem_calloc(1, sizeof(struct schiff_mesh)); - if(!ctx->mesh) { + + if(!sa_add(ctx->meshes, ngeoms)) { fprintf(stderr, -"Not enough memory: couldn't allocate the mesh of the geometry.\n"); + "Couldn't allocate the meshes of the geometry distributions.\n"); res = RES_MEM_ERR; goto error; } - /* Generate the mesh template and setup its distribution context */ - switch(geom->type) { - case SCHIFF_CYLINDER: - case SCHIFF_CYLINDER_AS_SPHERE: - res = schiff_mesh_init_cylinder - (&mem_default_allocator, ctx->mesh, geom->data.cylinder.nslices); - distrib->sample = geometry_sample_cylinder; - break; - case SCHIFF_SPHERE: - res = schiff_mesh_init_sphere - (&mem_default_allocator, ctx->mesh, geom->data.sphere.nslices); - distrib->sample = geometry_sample_sphere; - break; - case SCHIFF_SUPER_SHAPE: - res = schiff_mesh_init_sphere_polar - (&mem_default_allocator, ctx->mesh, geom->data.super_shape.nslices); - distrib->sample = geometry_sample_super_shape; - break; - default: FATAL("Unreachable code.\n"); break; - } - if(res != RES_OK) { - fprintf(stderr, "Couldn't initialise the mesh template.\n"); + if(!sa_add(ctx->geometries, ngeoms)) { + fprintf(stderr, +"Couldn't allocate the array that reference the geometry distributions.\n"); + res = RES_MEM_ERR; goto error; } + + FOR_EACH(i, 0, ngeoms) { + + /* Generate the mesh template and setup its distribution context */ + switch(geoms[i].type) { + case SCHIFF_CYLINDER: + case SCHIFF_CYLINDER_AS_SPHERE: + res = schiff_mesh_init_cylinder(&mem_default_allocator, &ctx->meshes[i], + geoms[i].data.cylinder.nslices); + break; + case SCHIFF_SPHERE: + res = schiff_mesh_init_sphere(&mem_default_allocator, &ctx->meshes[i], + geoms[i].data.sphere.nslices); + break; + case SCHIFF_SUPER_SHAPE: + res = schiff_mesh_init_sphere_polar(&mem_default_allocator, + &ctx->meshes[i], geoms[i].data.super_shape.nslices); + break; + default: FATAL("Unreachable code.\n"); break; + } + if(res != RES_OK) { + fprintf(stderr, "Couldn't initialise the mesh template.\n"); + goto error; + } + ctx->geometries[i] = geoms + i; + } ctx->properties = properties; - ctx->geometry = geom; + SSP(ran_discrete_ref_get(ran_geoms)); + ctx->ran_geometries = ran_geoms; + + distrib->material.get_property = get_material_property; + distrib->material.material = properties; distrib->context = ctx; + distrib->sample = geometry_sample; + distrib->caracteristic_length = caracteristic_length; exit: return res; error: - if(ctx) { - if(ctx->mesh) { - schiff_mesh_release(ctx->mesh); - mem_rm(ctx->mesh); - } - mem_rm(ctx); - ctx = NULL; - } + schiff_geometry_distribution_release(distrib); goto exit; } -void -schiff_geometry_distribution_release - (struct sschiff_geometry_distribution* distrib) -{ - struct geometry_distribution_context* ctx; - ASSERT(distrib); - ctx = distrib->context; - schiff_mesh_release(ctx->mesh); - mem_rm(ctx->mesh); - mem_rm(ctx); -} - diff --git a/src/schiff_geometry.h b/src/schiff_geometry.h @@ -30,6 +30,7 @@ #define SCHIFF_GEOMETRY_H #include <rsys/rsys.h> +#include <star/ssp.h> enum schiff_param_distribution { SCHIFF_PARAM_CONSTANT, @@ -121,6 +122,9 @@ extern LOCAL_SYM res_T schiff_geometry_distribution_init (struct sschiff_geometry_distribution* distrib, const struct schiff_geometry* geometry, + const size_t ngeoms, + const double caracteristic_length, + struct ssp_ran_discrete* ran_geoms, struct schiff_optical_properties* properties); extern LOCAL_SYM void diff --git a/src/schiff_mesh.h b/src/schiff_mesh.h @@ -80,7 +80,7 @@ schiff_mesh_release static INLINE void schiff_mesh_get_indices - (struct schiff_mesh* mesh, + (const struct schiff_mesh* mesh, const unsigned itri, unsigned ids[3]) { @@ -88,14 +88,14 @@ schiff_mesh_get_indices ASSERT(mesh); ASSERT(darray_uint_size_get(&mesh->indices) % 3 == 0); ASSERT(itri < darray_uint_size_get(&mesh->indices) / 3); - ids[0] = darray_uint_data_get(&mesh->indices)[i + 0]; - ids[1] = darray_uint_data_get(&mesh->indices)[i + 1]; - ids[2] = darray_uint_data_get(&mesh->indices)[i + 2]; + ids[0] = darray_uint_cdata_get(&mesh->indices)[i + 0]; + ids[1] = darray_uint_cdata_get(&mesh->indices)[i + 1]; + ids[2] = darray_uint_cdata_get(&mesh->indices)[i + 2]; } static INLINE void schiff_mesh_get_cartesian_position - (struct schiff_mesh* mesh, + (const struct schiff_mesh* mesh, const unsigned ivert, float vertex[3]) { @@ -103,14 +103,14 @@ schiff_mesh_get_cartesian_position ASSERT(mesh && vertex && mesh->coordinates == SCHIFF_CARTESIAN); ASSERT(darray_float_size_get(&mesh->vertices.cartesian) % 3 == 0); ASSERT(ivert < darray_float_size_get(&mesh->vertices.cartesian) / 3); - vertex[0] = darray_float_data_get(&mesh->vertices.cartesian)[i + 0]; - vertex[1] = darray_float_data_get(&mesh->vertices.cartesian)[i + 1]; - vertex[2] = darray_float_data_get(&mesh->vertices.cartesian)[i + 2]; + vertex[0] = darray_float_cdata_get(&mesh->vertices.cartesian)[i + 0]; + vertex[1] = darray_float_cdata_get(&mesh->vertices.cartesian)[i + 1]; + vertex[2] = darray_float_cdata_get(&mesh->vertices.cartesian)[i + 2]; } static INLINE void schiff_mesh_get_polar_position - (struct schiff_mesh* mesh, + (const struct schiff_mesh* mesh, const unsigned ivert, struct sin_cos angles[2]) { @@ -121,8 +121,8 @@ schiff_mesh_get_polar_position ASSERT(ivert < darray_uint_size_get(&mesh->vertices.polar) / 2); iangles = darray_uint_cdata_get(&mesh->vertices.polar) + i; - angles[0] = darray_sincos_data_get(&mesh->thetas)[iangles[0]]; - angles[1] = darray_sincos_data_get(&mesh->phis)[iangles[1]]; + angles[0] = darray_sincos_cdata_get(&mesh->thetas)[iangles[0]]; + angles[1] = darray_sincos_cdata_get(&mesh->phis)[iangles[1]]; } #endif /* SBOX_SCHIFF_MESH_H */