schiff.c (13169B)
1 /* Copyright (C) 2015-2016 CNRS 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "schiff_args.h" 17 #include "schiff_geometry.h" 18 #include "schiff_optical_properties.h" 19 20 #include <rsys/stretchy_array.h> 21 22 #include <star/s3d.h> 23 #include <star/sschiff.h> 24 #include <star/ssp.h> 25 26 /******************************************************************************* 27 * Helper functions 28 ******************************************************************************/ 29 static res_T 30 dump_geometries 31 (const struct schiff_args* args, 32 struct s3d_device* s3d, 33 struct sschiff_geometry_distribution* distrib, 34 struct ssp_rng* rng, 35 FILE* stream) 36 { 37 struct s3d_scene* scn = NULL; 38 struct s3d_shape* shape = NULL; 39 unsigned i; 40 res_T res = RES_OK; 41 ASSERT(args && distrib && rng && stream); 42 43 res = s3d_scene_create(s3d, &scn); 44 if(res != RES_OK) { 45 fprintf(stderr, "Couldn't create the Star-3D scene.\n"); 46 goto error; 47 } 48 res = s3d_shape_create_mesh(s3d, &shape); 49 if(res != RES_OK) { 50 fprintf(stderr, "Couldn't create the Star-3D shape.\n"); 51 goto error; 52 } 53 res = s3d_scene_attach_shape(scn, shape); 54 if(res != RES_OK) { 55 fprintf(stderr, "Couldn't attach the Star-3D shape to the Star-3D scene.\n"); 56 goto error; 57 } 58 59 FOR_EACH(i, 0, args->ngeoms_dump) { 60 unsigned ivert, nverts; 61 unsigned itri, ntris; 62 double volume_scaling; 63 double dist_scaling; 64 res = distrib->sample(rng, shape, &volume_scaling, distrib->context); 65 if(res != RES_OK) { 66 fprintf(stderr, "Couldn't sample the micro organism geometry.\n"); 67 goto error; 68 } 69 70 dist_scaling = pow(volume_scaling, 1.0/3.0); 71 72 fprintf(stream, "g schiff_geometry_%u\n", i); 73 74 /* Dump vertex position */ 75 fprintf(stream, "# List of vertex positions\n"); 76 S3D(mesh_get_vertices_count(shape, &nverts)); 77 FOR_EACH(ivert, 0, nverts) { 78 struct s3d_attrib attr; 79 S3D(mesh_get_vertex_attrib(shape, ivert, S3D_POSITION, &attr)); 80 f3_mulf(attr.value, attr.value, (float)dist_scaling); 81 fprintf(stream, "v %f %f %f\n", SPLIT3(attr.value)); 82 } 83 84 /* Dump triangle indices */ 85 fprintf(stream, "# Vertex indices of the triangular faces\n"); 86 S3D(mesh_get_triangles_count(shape, &ntris)); 87 FOR_EACH(itri, 0, ntris) { 88 unsigned ids[3]; 89 S3D(mesh_get_triangle_indices(shape, itri, ids)); 90 fprintf(stream, "f %u %u %u\n", ids[0]+1, ids[1]+1, ids[2]+1); 91 } 92 } 93 94 exit: 95 if(scn) S3D(scene_ref_put(scn)); 96 if(shape) S3D(shape_ref_put(shape)); 97 return res; 98 error: 99 goto exit; 100 } 101 102 static void 103 write_cross_sections 104 (FILE* stream, 105 struct sschiff_estimator* estimator, 106 const double* wlens, /* List of wavelengths in microns */ 107 const size_t nwlens) /* #wavelengths */ 108 { 109 size_t iwlen; 110 ASSERT(stream && estimator && wlens && nwlens); 111 112 FOR_EACH(iwlen, 0, nwlens) { 113 const struct sschiff_state* val; 114 struct sschiff_cross_section cross_section; 115 116 SSCHIFF(estimator_get_cross_section(estimator, iwlen, &cross_section)); 117 fprintf(stream, "%g ", wlens[iwlen]); 118 119 val = &cross_section.extinction; 120 fprintf(stream, "%g %g ", val->E, val->SE); 121 val = &cross_section.absorption; 122 fprintf(stream, "%g %g ", val->E, val->SE); 123 val = &cross_section.scattering; 124 fprintf(stream, "%g %g ", val->E, val->SE); 125 val = &cross_section.average_projected_area; 126 fprintf(stream, "%g %g\n", val->E, val->SE); 127 } 128 fprintf(stream, "\n"); 129 } 130 131 static void 132 write_phase_function_descriptors 133 (FILE* stream, 134 struct sschiff_estimator* estimator, 135 const double* wlens, /* List of wavelengths in microns */ 136 const size_t nwlens, /* #wavelenghts */ 137 const size_t nangles_inv) /* Number of inversed cumulative entries */ 138 { 139 size_t iwlen; 140 const double* angles; 141 size_t nangles; 142 res_T res = RES_OK; 143 ASSERT(stream && estimator && wlens && nwlens && nangles_inv > 1); 144 145 SSCHIFF(estimator_get_scattering_angles(estimator, &angles, &nangles)); 146 147 FOR_EACH(iwlen, 0, nwlens) { 148 size_t iangle; 149 struct sschiff_state cdf, pf; 150 double n; 151 152 res = sschiff_estimator_get_limit_scattering_angle_index 153 (estimator, iwlen, &iangle); 154 155 if(res == RES_BAD_OP) { 156 /* No valid limit angle => no phase function. Print -1 for all the data */ 157 fprintf(stderr, 158 "The phase function couldn't be computed.\n" 159 "wavelength = %g microns\n", wlens[iwlen]); 160 fprintf(stream, "-1 -1 -1 -1 -1 -1 -1 -1 -1\n"); 161 162 } else { 163 /* Retrieve the phase function descriptor data */ 164 ASSERT(res == RES_OK); 165 SSCHIFF(estimator_get_wide_scattering_angle_model_parameter 166 (estimator, iwlen, &n)); 167 SSCHIFF(estimator_get_differential_cross_section 168 (estimator, iwlen, iangle, &pf)); 169 SSCHIFF(estimator_get_differential_cross_section_cumulative 170 (estimator, iwlen, iangle, &cdf)); 171 172 /* Write the phase function descriptor */ 173 fprintf(stream, "%g %g %g %g %g %g %g %lu %lu\n", 174 wlens[iwlen], angles[iangle], pf.E, pf.SE, cdf.E, cdf.SE, n, 175 (unsigned long)nangles, 176 (unsigned long)nangles_inv); 177 } 178 } 179 fprintf(stream, "\n"); /* Empty line */ 180 } 181 182 static void 183 write_phase_functions 184 (FILE* stream, 185 const struct sschiff_estimator* estimator, 186 const double* wlens, 187 const size_t nwlens) 188 { 189 size_t iwlen; 190 const double* angles; 191 const struct sschiff_state* phase_funcs; 192 size_t iangle, nangles; 193 ASSERT(stream && estimator && wlens && nwlens); 194 (void)wlens; 195 196 SSCHIFF(estimator_get_scattering_angles(estimator, &angles, &nangles)); 197 FOR_EACH(iwlen, 0, nwlens) { 198 const res_T res = sschiff_estimator_get_phase_function 199 (estimator, iwlen, &phase_funcs); 200 if(res == RES_BAD_OP) { 201 /* The phase function couldn't be computed. Write default data. */ 202 FOR_EACH(iangle, 0, nangles) { 203 fprintf(stream, "-1 -1 -1\n"); 204 } 205 } else { 206 ASSERT(res == RES_OK); 207 FOR_EACH(iangle, 0, nangles) { 208 fprintf(stream, "%g %g %g\n", 209 angles[iangle], 210 phase_funcs[iangle].E, 211 phase_funcs[iangle].SE); 212 } 213 } 214 fprintf(stream, "\n"); 215 } 216 } 217 218 static void 219 write_phase_functions_cum 220 (FILE* stream, 221 struct sschiff_estimator* estimator, 222 const double* wlens, 223 const size_t nwlens) 224 { 225 size_t iwlen; 226 const double* angles; 227 const struct sschiff_state* phase_funcs_cum; 228 size_t iangle, nangles; 229 ASSERT(stream && estimator && wlens && nwlens); 230 (void)wlens; 231 232 SSCHIFF(estimator_get_scattering_angles(estimator, &angles, &nangles)); 233 FOR_EACH(iwlen, 0, nwlens) { 234 const res_T res = sschiff_estimator_get_phase_function_cumulative 235 (estimator, iwlen, &phase_funcs_cum); 236 if(res == RES_BAD_OP) { 237 /* The cumulative phase func couldn't be computed. Write default data */ 238 FOR_EACH(iangle, 0, nangles) { 239 fprintf(stream, "-1 -1 -1\n"); 240 } 241 } else { 242 ASSERT(res == RES_OK); 243 FOR_EACH(iangle, 0, nangles) { 244 fprintf(stream, "%g %g %g\n", 245 angles[iangle], 246 phase_funcs_cum[iangle].E, 247 phase_funcs_cum[iangle].SE); 248 } 249 } 250 fprintf(stream, "\n"); 251 } 252 } 253 254 static void 255 write_phase_functions_invcum 256 (FILE* stream, 257 struct sschiff_estimator* estimator, 258 const double* wlens, 259 const size_t nwlens, 260 const size_t nangles_inv) 261 { 262 double* invcum = NULL; 263 size_t iangle; 264 size_t iwlen; 265 double step; 266 ASSERT(stream && estimator && wlens && nwlens && nangles_inv > 1); 267 268 step = 1.0 / (double)(nangles_inv - 1); 269 270 if(!sa_add(invcum, nangles_inv)) { 271 fprintf(stderr, 272 "Couldn't allocate the list of inverse cumulative phase function angles.\n"); 273 274 /* Print the default -1 value */ 275 FOR_EACH(iwlen, 0, nwlens) { 276 FOR_EACH(iangle, 0, nangles_inv) { 277 fprintf(stream, "-1 -1\n"); 278 } 279 fprintf(stream, "\n"); 280 } 281 } else { 282 FOR_EACH(iwlen, 0, nwlens) { 283 const res_T res = sschiff_estimator_inverse_cumulative_phase_function 284 (estimator, iwlen, invcum, nangles_inv); 285 if(res != RES_OK) { 286 /* The inverse cumulative cannot be computed. Write -1 */ 287 fprintf(stderr, 288 "Couldn't inverse the cumulative phase function.\n" 289 "wavelength = %g microns\n", wlens[iwlen]); 290 FOR_EACH(iangle, 0, nangles_inv) { 291 fprintf(stream, "-1 -1\n"); 292 } 293 } else { 294 /* Write the inverse cumulative values */ 295 FOR_EACH(iangle, 0, nangles_inv-1) { 296 fprintf(stream, "%g %g\n", (double)iangle*step, invcum[iangle]); 297 } 298 /* Handle precision issues */ 299 fprintf(stream, "1 %g\n", invcum[nangles_inv-1]); 300 } 301 fprintf(stream, "\n"); /* Empty line */ 302 } 303 } 304 305 if(invcum) 306 sa_release(invcum); 307 } 308 309 static res_T 310 run_integration 311 (const struct schiff_args* args, 312 struct s3d_device* s3d, 313 struct sschiff_geometry_distribution* distrib, 314 struct ssp_rng* rng, 315 FILE* stream) 316 { 317 struct sschiff_device* sschiff = NULL; 318 struct sschiff_estimator* estimator = NULL; 319 size_t nwlens; 320 res_T res = RES_OK; 321 ASSERT(args && sa_size(args->wavelengths) && rng && stream); 322 323 res = sschiff_device_create(NULL, NULL, args->nthreads, 1, s3d, &sschiff); 324 if(res != RES_OK) { 325 fprintf(stderr, "Couldn't create the Star Schiff device.\n"); 326 goto error; 327 } 328 329 /* Invoke the Schiff integration */ 330 nwlens = sa_size(args->wavelengths); 331 res = sschiff_integrate(sschiff, rng, distrib, args->wavelengths, nwlens, 332 sschiff_uniform_scattering_angles, args->nangles, args->nrealisations, 333 args->ndirs, &estimator); 334 if(res != RES_OK) { 335 fprintf(stderr, "Schiff integration error.\n"); 336 goto error; 337 } 338 339 /* Write results */ 340 write_cross_sections(stream, estimator, args->wavelengths, nwlens); 341 write_phase_function_descriptors 342 (stream, estimator, args->wavelengths, nwlens, args->nangles_inv); 343 write_phase_functions(stream, estimator, args->wavelengths, nwlens); 344 write_phase_functions_cum(stream, estimator, args->wavelengths, nwlens); 345 write_phase_functions_invcum 346 (stream, estimator, args->wavelengths, nwlens, args->nangles_inv); 347 348 exit: 349 if(estimator) SSCHIFF(estimator_ref_put(estimator)); 350 if(sschiff) SSCHIFF(device_ref_put(sschiff)); 351 return res; 352 error: 353 goto exit; 354 } 355 356 static res_T 357 run(const struct schiff_args* args) 358 { 359 FILE* fp = stdout; 360 struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION; 361 struct ssp_rng* rng = NULL; 362 struct s3d_device* s3d = NULL; 363 res_T res = RES_OK; 364 ASSERT(args); 365 366 if(args->output_filename) { 367 fp = fopen(args->output_filename, "w"); 368 if(!fp) { 369 fprintf(stderr, 370 "Couldn't open the output file `%s'.\n", args->output_filename); 371 res = RES_IO_ERR; 372 goto error; 373 } 374 } 375 376 res = ssp_rng_create(&mem_default_allocator, &ssp_rng_threefry, &rng); 377 if(res != RES_OK) { 378 fprintf(stderr, "Couldn't create the Random Number Generator.\n"); 379 goto error; 380 } 381 382 res = s3d_device_create(NULL, &mem_default_allocator, 0, &s3d); 383 if(res != RES_OK) { 384 fprintf(stderr, "Couldn't create the Star-3D device.\n"); 385 goto error; 386 } 387 388 res = schiff_geometry_distribution_init(&distrib, s3d, args->geoms, 389 sa_size(args->geoms), args->characteristic_length, args->ran_geoms, 390 args->properties); 391 if(res != RES_OK) { 392 fprintf(stderr, 393 "Couldn't initialise the Star Schiff geometry distribution.\n"); 394 goto error; 395 } 396 397 if(args->ngeoms_dump) { 398 res = dump_geometries(args, s3d, &distrib, rng, fp); 399 } else { 400 res = run_integration(args, s3d, &distrib, rng, fp); 401 } 402 403 /* Release before the check of the integration error code since if an error 404 * occurred the function will return without releasing the distribution. */ 405 schiff_geometry_distribution_release(&distrib); 406 407 if(res != RES_OK) 408 goto error; 409 410 exit: 411 if(fp && fp != stdout) fclose(fp); 412 if(rng) SSP(rng_ref_put(rng)); 413 if(s3d) S3D(device_ref_put(s3d)); 414 return res; 415 error: 416 goto exit; 417 } 418 419 /******************************************************************************* 420 * Entry point 421 ******************************************************************************/ 422 int 423 main(int argc, char** argv) 424 { 425 struct schiff_args args = SCHIFF_ARGS_NULL; 426 int err = 0; 427 res_T res; 428 429 res = schiff_args_init(&args, argc, argv); 430 if(res != RES_OK) goto error; 431 if(!args.ngeoms_dump && !args.wavelengths) goto exit; 432 res = run(&args); 433 if(res != RES_OK) goto error; 434 435 exit: 436 schiff_args_release(&args); 437 if(mem_allocated_size() != 0) { 438 fprintf(stderr, "Memory leaks: %lu Bytes\n", 439 (unsigned long) mem_allocated_size()); 440 } 441 return err; 442 error: 443 err = 1; 444 goto exit; 445 } 446