commit 4e2d6a0d2e499884a3fb25d97afd8fa85def6282
parent 7d6f14f13c3388d420cd674380d22e4003be4d68
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 19 Nov 2015 16:29:55 +0100
Major update of the geometry distribution
Make configurable the distribution of the geometry parameter. The
following 3 distributions are currently supported: constant, lognormal
and histogram.
Diffstat:
7 files changed, 573 insertions(+), 152 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -66,11 +66,13 @@ set(SCHIFF_FILES_SRC
schiff.c
schiff_args.c
schiff_geometry.c
+ schiff_histogram.c
schiff_mesh.c
schiff_optical_properties.c)
set(SCHIFF_FILES_INC
schiff_args.h
schiff_geometry.h
+ schiff_histogram.h
schiff_mesh.h
schiff_optical_properties.h
schiff_streamline.h)
diff --git a/src/schiff_args.c b/src/schiff_args.c
@@ -33,6 +33,7 @@
#include <rsys/dynamic_array_char.h>
#include <rsys/cstr.h>
+#include <rsys/str.h>
#include <rsys/stretchy_array.h>
#ifdef COMPILER_CL
@@ -70,12 +71,11 @@ printf(
"area. The \"E'\", \"A'\", \"S'\" and \"P'\" values are the standard error of the\n"
"aforementioned estimations.\n\n");
printf(
-" -c R:S:A[:N] the micro organisms are cylindrical meshes discretized in N\n"
-" slices. By default N is %u. The radius `r' and the height `h'\n"
-" of the cylinders are distributed as follow:\n"
-" y = 1/(ln(S)*x*sqrt(2PI)) * exp(-(ln(x)-ln(R))^2 / (2*ln(S)^2))\n"
-" r = y / (3 / (2*A))^1/3\n"
-" h = 2*r/A\n",
+" -c RADIUS,HEIGHT[,N]\n"
+" the micro organisms are cylindrical meshes discretized in N\n"
+" slices. By default N is %u. RADIUS and HEIGHT are the\n"
+" distributions of, respectively, the radius and the height\n"
+" parameters of the cylinder.\n",
SCHIFF_CYLINDER_DEFAULT.nslices);
printf(
" -d DIRS number of sampled directions for each geometry. Default is %u.\n",
@@ -98,18 +98,15 @@ printf(
printf(
" -q do not print the helper message when no FILE is submitted.\n");
printf(
-" -s R:S[:N] the micro organisms are spherical meshes discretized in N\n"
-" slices along 2PI. By default N is %u. Their radius `r' is\n"
-" distributed with respect to a lognormal distribution:\n"
-" r = 1/(ln(S)*x*sqrt(2PI)) * exp(-(ln(x)-ln(R))^2 / (2*ln(S)^2))\n",
+" -s RADIUS[,N] the micro organisms are spherical meshes discretized in N\n"
+" slices along 2PI. By default N is %u. Their radius parameter\n"
+" is distributed with respect to the distribution RADIUS.\n",
SCHIFF_SPHERE_DEFAULT.nslices);
printf(
-" -u M_min:M_max:N2_min:N2_max[:N]\n"
-" the micro organisms are 3D super shapes discretized in N slices\n"
+" -u M,N2[,N] the micro organisms are 3D super shapes discretized in N slices\n"
" along 2PI. By default N is %u. The A, B, N0 and N1 parameters\n"
-" of the super formulas are fixed to 1 while their M and N2\n"
-" variables are uniformally distributed in [M_min, M_max] and\n"
-" [N2_min, N2_max], respectively.\n",
+" of the super formulas are fixed to 1 while the M and N2 control\n"
+" de distribution of the remaining parameters.",
SCHIFF_SUPER_SHAPE_DEFAULT.nslices);
printf(
" -w A[:B]... list of wavelengths in vacuum (expressed in micron) to\n"
@@ -117,47 +114,106 @@ printf(
}
static res_T
+parse_param_distribution(const char* str, struct schiff_param* param)
+{
+ struct str buf;
+ char* tk;
+ res_T res = RES_OK;
+ ASSERT(str && param);
+ param->distribution = SCHIFF_PARAM_NONE;
+
+ str_init(&mem_default_allocator, &buf);
+ res = str_set(&buf, str);
+ if(res != RES_OK) {
+ fprintf(stderr,
+ "Not enough memory. Couldn't parse the parameter distribution `%s'.\n",
+ str);
+ goto error;
+ }
+
+ tk = strtok(str_get(&buf), "=");
+ if(!strcmp(tk, "c")) { /* Constant */
+ param->distribution = SCHIFF_PARAM_CONSTANT;
+ res = cstr_to_double(strtok(NULL, "\0"), ¶m->data.constant);
+ } else if(!strcmp(tk, "h")) { /* Histogram */
+ param->distribution = SCHIFF_PARAM_HISTOGRAM;
+ darray_hentry_init(&mem_default_allocator, ¶m->data.histogram);
+ res = schiff_histogram_load(¶m->data.histogram, strtok(NULL, "\0"));
+ } else if(!strcmp(tk, "l")) { /* lognormal */
+ double list[2];
+ size_t len;
+ param->distribution = SCHIFF_PARAM_LOGNORMAL;
+ res = cstr_to_list_double(strtok(NULL, "\0"), list, &len, 2);
+ if(res == RES_OK) {
+ if(len != 2) {
+ res = RES_BAD_ARG;
+ } else {
+ param->data.lognormal.zeta = list[0];
+ param->data.lognormal.sigma = list[1];
+ }
+ }
+ } else {
+ res = RES_BAD_ARG;
+ }
+
+ if(res != RES_OK) {
+ fprintf(stderr, "Invalid parameter distribution `%s'.\n", str);
+ goto error;
+ }
+
+exit:
+ str_release(&buf);
+ return res;
+error:
+ if(param->distribution == SCHIFF_PARAM_HISTOGRAM)
+ darray_hentry_release(¶m->data.histogram);
+ goto exit;
+}
+
+static res_T
parse_cylinder_distribution(const char* str, struct schiff_args* args)
{
- double list[4];
- size_t len;
- unsigned nslices;
+ struct str buf;
+ char* tk;
+ char* tk_ctx;
res_T res = RES_OK;
ASSERT(args && str);
- /* The distribution was already set to a non cylinder distribution. */
- if(args->geom.type!= SCHIFF_NONE && args->geom.type != SCHIFF_CYLINDER) {
+ str_init(&mem_default_allocator, &buf);
+
+ /* The distribution was already set to a non sphere distribution. */
+ if(args->geom.type != SCHIFF_NONE && args->geom.type != SCHIFF_CYLINDER) {
res = RES_BAD_ARG;
goto error;
-
}
- res = cstr_to_list_double(str, list, &len, 4);
- if(res != RES_OK) goto error; /* Wrong argument formating */
- /* At least the radius the sigma and the aspect_ratio must be provided */
- if(len < 3) {
- res = RES_BAD_ARG;
+ res = str_set(&buf, str);
+ if(res != RES_OK) {
+ fprintf(stderr,
+ "Not enough memory. Couldn't parse the cylinder distribution `%s'.\n",
+ str);
goto error;
}
- if(len == 3) {
- nslices = SCHIFF_CYLINDER_DEFAULT.nslices;
+ tk = strtok_r(str_get(&buf), ",", &tk_ctx);
+ res = parse_param_distribution(tk, &args->geom.data.cylinder.radius);
+ if(res != RES_OK) goto error;
+
+ tk = strtok_r(NULL, ",", &tk_ctx);
+ res = parse_param_distribution(tk, &args->geom.data.cylinder.height);
+ if(res != RES_OK) goto error;
+
+ tk = strtok_r(NULL, "\0", &tk_ctx);
+ if(tk) {
+ args->geom.data.cylinder.nslices = SCHIFF_CYLINDER_DEFAULT.nslices;
} else {
- nslices = (unsigned)list[3];
- if((double)nslices != list[3]) { /* The #slices arg is not an integer */
- res = RES_BAD_ARG;
- goto error;
- }
+ res = cstr_to_uint(tk, &args->geom.data.cylinder.nslices);
+ if(res != RES_OK) goto error;
}
-
- /* Update the distribution */
args->geom.type = SCHIFF_CYLINDER;
- args->geom.data.cylinder.radius = list[0];
- args->geom.data.cylinder.sigma = list[1];
- args->geom.data.cylinder.aspect_ratio = list[2];
- args->geom.data.cylinder.nslices = nslices;
exit:
+ str_release(&buf);
return res;
error:
goto exit;
@@ -166,43 +222,43 @@ error:
static res_T
parse_sphere_distribution(const char* str, struct schiff_args* args)
{
- double list[3];
- size_t len;
- unsigned nslices;
+ struct str buf;
+ char* tk;
+ char* tk_ctx;
res_T res = RES_OK;
ASSERT(args && str);
+ str_init(&mem_default_allocator, &buf);
+
/* The distribution was already set to a non sphere distribution. */
if(args->geom.type != SCHIFF_NONE && args->geom.type != SCHIFF_SPHERE) {
res = RES_BAD_ARG;
goto error;
-
}
- res = cstr_to_list_double(str, list, &len, 3);
- if(res != RES_OK) goto error; /* Wrong argument formating */
- if(len < 2) {/* At least the radius and the sigma argument must be provided */
- res = RES_BAD_ARG;
+ res = str_set(&buf, str);
+ if(res != RES_OK) {
+ fprintf(stderr,
+ "Not enough memory. Couldn't parse the sphere distribution `%s'.\n",
+ str);
goto error;
}
- if(len == 2) {
- nslices = SCHIFF_SPHERE_DEFAULT.nslices;
+ tk = strtok_r(str_get(&buf), ",", &tk_ctx);
+ res = parse_param_distribution(tk, &args->geom.data.sphere.radius);
+ if(res != RES_OK) goto error;
+
+ tk = strtok_r(NULL, "\0", &tk_ctx);
+ if(!tk) {
+ args->geom.data.sphere.nslices = SCHIFF_SPHERE_DEFAULT.nslices;
} else {
- nslices = (unsigned)list[2];
- if((double)nslices != list[2]) { /* The #slices arg is not an integer */
- res = RES_BAD_ARG;
- goto error;
- }
+ res = cstr_to_uint(tk, &args->geom.data.sphere.nslices);
+ if(res != RES_OK) goto error;
}
-
- /* Update the distribution */
args->geom.type = SCHIFF_SPHERE;
- args->geom.data.sphere.radius = list[0];
- args->geom.data.sphere.sigma = list[1];
- args->geom.data.sphere.nslices = nslices;
exit:
+ str_release(&buf);
return res;
error:
goto exit;
@@ -211,49 +267,47 @@ error:
static res_T
parse_super_shape_distribution(const char* str, struct schiff_args* args)
{
- double list[5];
- size_t len;
- unsigned nslices;
+ struct str buf;
+ char* tk;
+ char* tk_ctx;
res_T res = RES_OK;
ASSERT(args && str);
- /* The distribution was already set to a non super_shape distribution. */
+ str_init(&mem_default_allocator, &buf);
+
+ /* The distribution was already set to a non sphere distribution. */
if(args->geom.type != SCHIFF_NONE && args->geom.type != SCHIFF_SUPER_SHAPE) {
res = RES_BAD_ARG;
goto error;
-
}
- res = cstr_to_list_double(str, list, &len, 5);
- if(res != RES_OK) goto error; /* Wrong argument formating */
- if(len < 4) {/* At least the bounds of the M and N2 params must be provided */
- res = RES_BAD_ARG;
- goto error;
- }
- if(list[0] > list[1] || list[2] > list[3]) {
- res = RES_BAD_ARG;
+ res = str_set(&buf, str);
+ if(res != RES_OK) {
+ fprintf(stderr,
+ "Not enough memory. Couldn't parse the super shape distribution `%s'.\n",
+ str);
goto error;
}
- if(len == 4) {
- nslices = SCHIFF_SPHERE_DEFAULT.nslices;
+ tk = strtok_r(str_get(&buf), ",", &tk_ctx);
+ res = parse_param_distribution(tk, &args->geom.data.super_shape.M);
+ if(res != RES_OK) goto error;
+
+ tk = strtok_r(NULL, ",", &tk_ctx);
+ res = parse_param_distribution(tk, &args->geom.data.super_shape.N2);
+ if(res != RES_OK) goto error;
+
+ tk = strtok_r(NULL, "\0", &tk_ctx);
+ if(tk) {
+ args->geom.data.super_shape.nslices = SCHIFF_SUPER_SHAPE_DEFAULT.nslices;
} else {
- nslices = (unsigned)list[4];
- if((double)nslices != list[4]) { /* The #slices arg is not an integer */
- res = RES_BAD_ARG;
- goto error;
- }
+ res = cstr_to_uint(tk, &args->geom.data.super_shape.nslices);
+ if(res != RES_OK) goto error;
}
-
- /* Update the distribution */
args->geom.type = SCHIFF_SUPER_SHAPE;
- args->geom.data.super_shape.lower_M = list[0];
- args->geom.data.super_shape.upper_M = list[1];
- args->geom.data.super_shape.lower_N2 = list[2];
- args->geom.data.super_shape.upper_N2 = list[3];
- args->geom.data.super_shape.nslices = nslices;
exit:
+ str_release(&buf);
return res;
error:
goto exit;
diff --git a/src/schiff_geometry.c b/src/schiff_geometry.c
@@ -34,6 +34,8 @@
#include <star/ssp.h>
#include <star/sschiff.h>
+enum { A, B, M, N0, N1, N2 }; /* Super formula arguments */
+
/* 3D Context of the sphere geometry */
struct sphere_context {
float radius; /* Sphere radius */
@@ -50,12 +52,6 @@ struct super_shape_context {
double formulas[2][6];
};
-/* Distribution context of a sphere geometry */
-struct sphere_distribution_context {
- double log_mean_radius; /* Log of the mean radius */
- double log_sigma; /* Log of the sigma argument of the lognormal distribution*/
-};
-
/* Distribution context of a cylinder geometry */
struct cylinder_distribution_context {
double log_mean_radius; /* Log of the mean radius */
@@ -63,17 +59,10 @@ struct cylinder_distribution_context {
double aspect_ratio; /* aspect ratio of the cylinder distribution */
};
-/* Distribtion context of a super shape geometry */
-struct super_shape_distribution_context {
- double lower_M;
- double upper_M;
- double lower_N2;
- double upper_N2;
-};
-
/* 3D context of a generic geometry */
struct geometry_context {
struct schiff_mesh* mesh; /* Triangular mesh of the geometry */
+ struct schiff_geometry* geometry;
enum schiff_geometry_type type;
union {
struct cylinder_context cylinder;
@@ -84,19 +73,37 @@ struct geometry_context {
/* Distribution context of a generic geometry */
struct geometry_distribution_context {
- struct schiff_mesh* mesh; /* Triangular mesh of the cylindrical mesh */
+ struct schiff_mesh* mesh; /* Triangular mesh of the geometry */
struct schiff_optical_properties* properties; /* Per wavelength properties */
- enum schiff_geometry_type type;
- union {
- struct cylinder_distribution_context cylinder;
- struct sphere_distribution_context sphere;
- struct super_shape_distribution_context super_shape;
- } data;
+ const struct schiff_geometry* geometry; /* Geometry descriptor */
};
/*******************************************************************************
* Helper functions
******************************************************************************/
+static INLINE double
+eval_param(const struct schiff_param* param, struct ssp_rng* rng)
+{
+ double val = 0.0;
+ ASSERT(param && rng);
+
+ switch(param->distribution) {
+ case SCHIFF_PARAM_CONSTANT:
+ val = param->data.constant;
+ break;
+ case SCHIFF_PARAM_LOGNORMAL:
+ val = ssp_ran_lognormal(rng, log(param->data.lognormal.zeta),
+ log(param->data.lognormal.sigma));
+ break;
+ case SCHIFF_PARAM_HISTOGRAM:
+ val = schiff_histogram_sample
+ (¶m->data.histogram, ssp_rng_canonical(rng));
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return val;
+}
+
static void
get_material_property
(void* mtl,
@@ -173,21 +180,30 @@ geometry_sample_cylinder
void* context)
{
struct geometry_distribution_context* distrib = context;
+ const struct schiff_cylinder* cylinder;
struct geometry_context geom;
struct s3d_vertex_data attrib;
size_t nverts, nprims;
- double y;
+ double r;
ASSERT(rng && mtl && shape && context);
- ASSERT(distrib->type == SCHIFF_CYLINDER);
- y = ssp_ran_lognormal
- (rng, distrib->data.cylinder.log_mean_radius, distrib->data.cylinder.log_sigma);
+ cylinder = &distrib->geometry->data.cylinder;
geom.type = SCHIFF_CYLINDER;
geom.mesh = distrib->mesh;
- geom.data.cylinder.radius =
- (float)(y / pow(3.0 / (2.0*distrib->data.cylinder.aspect_ratio), 1.0/3.0));
- geom.data.cylinder.height =
- (float)(2.f * geom.data.cylinder.radius / distrib->data.cylinder.aspect_ratio);
+ switch(distrib->geometry->type) {
+ case SCHIFF_CYLINDER:
+ geom.data.cylinder.radius = (float)eval_param(&cylinder->radius, rng);
+ geom.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 =
+ (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);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
attrib.usage = S3D_POSITION;
attrib.type = S3D_FLOAT3;
@@ -211,16 +227,17 @@ geometry_sample_sphere
void* context)
{
struct geometry_distribution_context* distrib = context;
+ const struct schiff_sphere* sphere;
struct geometry_context geom;
struct s3d_vertex_data attrib;
size_t nverts, nprims;
ASSERT(rng && mtl && shape && context);
- ASSERT(distrib->type == SCHIFF_SPHERE);
+ ASSERT(distrib->geometry->type == SCHIFF_SPHERE);
+ sphere = &distrib->geometry->data.sphere;
geom.mesh = distrib->mesh;
geom.type = SCHIFF_SPHERE;
- geom.data.sphere.radius = (float)ssp_ran_lognormal
- (rng, distrib->data.sphere.log_mean_radius, distrib->data.sphere.log_sigma);
+ geom.data.sphere.radius = (float)eval_param(&sphere->radius, rng);
attrib.usage = S3D_POSITION;
attrib.type = S3D_FLOAT3;
@@ -246,12 +263,12 @@ geometry_sample_super_shape
struct geometry_distribution_context* distrib = context;
struct geometry_context geom;
struct s3d_vertex_data attrib;
- const struct super_shape_distribution_context* sshape;
+ const struct schiff_super_shape* sshape;
size_t nverts, nprims;
ASSERT(rng && mtl && shape && context);
- ASSERT(distrib->type == SCHIFF_SUPER_SHAPE);
+ ASSERT(distrib->geometry->type == SCHIFF_SUPER_SHAPE);
- sshape = &distrib->data.super_shape;
+ sshape = &distrib->geometry->data.super_shape;
geom.mesh = distrib->mesh;
geom.type = SCHIFF_SUPER_SHAPE;
@@ -265,14 +282,10 @@ geometry_sample_super_shape
geom.data.super_shape.formulas[1][N0] = 1.0;
geom.data.super_shape.formulas[1][N1] = 1.0;
- geom.data.super_shape.formulas[0][M] = (float)
- ssp_rng_uniform_double(rng, sshape->lower_M, sshape->upper_M);
- geom.data.super_shape.formulas[1][M] = (float)
- ssp_rng_uniform_double(rng, sshape->lower_M, sshape->upper_M);
- geom.data.super_shape.formulas[0][N2] = (float)
- ssp_rng_uniform_double(rng, sshape->lower_M, sshape->upper_N2);
- geom.data.super_shape.formulas[1][N2] = (float)
- ssp_rng_uniform_double(rng, sshape->lower_M, sshape->upper_N2);
+ geom.data.super_shape.formulas[0][M] = (float)eval_param(&sshape->M, rng);
+ geom.data.super_shape.formulas[1][M] = (float)eval_param(&sshape->M, rng);
+ geom.data.super_shape.formulas[0][N2] = (float)eval_param(&sshape->N2, rng);
+ geom.data.super_shape.formulas[1][N2] = (float)eval_param(&sshape->N2, rng);
attrib.usage = S3D_POSITION;
attrib.type = S3D_FLOAT3;
@@ -319,27 +332,19 @@ schiff_geometry_distribution_init
/* 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);
- ctx->data.cylinder.log_mean_radius = log(geom->data.cylinder.radius);
- ctx->data.cylinder.log_sigma = log(geom->data.cylinder.sigma);
- ctx->data.cylinder.aspect_ratio = geom->data.cylinder.aspect_ratio;
distrib->sample = geometry_sample_cylinder;
break;
case SCHIFF_SPHERE:
res = schiff_mesh_init_sphere
(&mem_default_allocator, ctx->mesh, geom->data.sphere.nslices);
- ctx->data.sphere.log_mean_radius = log(geom->data.sphere.radius);
- ctx->data.sphere.log_sigma = log(geom->data.sphere.sigma);
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);
- ctx->data.super_shape.lower_M = geom->data.super_shape.lower_M;
- ctx->data.super_shape.upper_M = geom->data.super_shape.upper_M;
- ctx->data.super_shape.lower_N2 = geom->data.super_shape.lower_N2;
- ctx->data.super_shape.upper_N2 = geom->data.super_shape.upper_N2;
distrib->sample = geometry_sample_super_shape;
break;
default: FATAL("Unreachable code.\n"); break;
@@ -349,7 +354,7 @@ schiff_geometry_distribution_init
goto error;
}
ctx->properties = properties;
- ctx->type = geom->type;
+ ctx->geometry = geom;
distrib->context = ctx;
exit:
diff --git a/src/schiff_geometry.h b/src/schiff_geometry.h
@@ -30,44 +30,61 @@
#define SCHIFF_GEOMETRY_H
#include <rsys/rsys.h>
+#include "schiff_histogram.h"
+
+enum schiff_param_distribution {
+ SCHIFF_PARAM_CONSTANT,
+ SCHIFF_PARAM_LOGNORMAL,
+ SCHIFF_PARAM_HISTOGRAM,
+ SCHIFF_PARAM_NONE
+};
+
+struct schiff_param {
+ enum schiff_param_distribution distribution;
+ union {
+ double constant;
+ struct { double zeta, sigma; } lognormal;
+ struct darray_hentry histogram;
+ } data;
+};
+#define SCHIFF_PARAM_DEFAULT__ {SCHIFF_PARAM_CONSTANT, {1.0}}
enum schiff_geometry_type {
SCHIFF_CYLINDER,
+ SCHIFF_CYLINDER_AS_SPHERE, /* The cylinder volume is control by a sphere */
SCHIFF_SPHERE,
SCHIFF_SUPER_SHAPE,
SCHIFF_NONE
};
struct schiff_sphere {
- double radius;
- double sigma;
+ struct schiff_param radius;
unsigned nslices;
};
-#define SCHIFF_SPHERE_DEFAULT__ {1.0, 1.0, 64}
+#define SCHIFF_SPHERE_DEFAULT__ {SCHIFF_PARAM_DEFAULT__, 64}
static const struct schiff_sphere SCHIFF_SPHERE_DEFAULT =
SCHIFF_SPHERE_DEFAULT__;
struct schiff_cylinder {
- double radius;
- double sigma;
- double aspect_ratio;
+ struct schiff_param radius;
+ struct schiff_param height; /* In use by SCHIFF_CYLINDER */
+ double aspect_ratio; /* In use by CYLINDER_AS_SPHERE */
unsigned nslices;
};
-#define SCHIFF_CYLINDER_DEFAULT__ {1.0, 1.0, 1.0, 64}
+#define SCHIFF_CYLINDER_DEFAULT__ \
+ {SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, -1.0, 64}
static const struct schiff_cylinder SCHIFF_CYLINDER_DEFAULT =
SCHIFF_CYLINDER_DEFAULT__;
-enum schiff_super_formula { A, B, M, N0, N1, N2 };
struct schiff_super_shape {
- double lower_M;
- double upper_M;
- double lower_N2;
- double upper_N2;
+ struct schiff_param M;
+ struct schiff_param N2;
unsigned nslices;
};
-#define SCHIFF_SUPER_SHAPE_DEFAULT__ {1.0, 9.0, 1.0, 9.0, 64}
+#define SCHIFF_SUPER_SHAPE_DEFAULT__ \
+ {SCHIFF_PARAM_DEFAULT__, SCHIFF_PARAM_DEFAULT__, 64}
static const struct schiff_super_shape SCHIFF_SUPER_SHAPE_DEFAULT =
SCHIFF_SUPER_SHAPE_DEFAULT__;
diff --git a/src/schiff_histogram.c b/src/schiff_histogram.c
@@ -0,0 +1,196 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#include "schiff_histogram.h"
+#include "schiff_streamline.h"
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/stretchy_array.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+parse_schiff_hentry
+ (struct schiff_hentry* entry,
+ const char* filename,
+ const size_t iline,
+ const char* str)
+{
+ char buf[128];
+ double tmp[2];
+ char* tk;
+ int i;
+ res_T res = RES_OK;
+ ASSERT(entry && filename && str);
+
+ if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
+ fprintf(stderr,
+ "%s:%lu: not enough memory: cannot parse the schiff histogram value `%s'.\n",
+ filename, (unsigned long)iline, str);
+ return RES_MEM_ERR;
+ }
+
+ strncpy(buf, str, sizeof(buf));
+ for(i=0, tk=strtok(buf, " \t"); tk && i < 2; ++i, tk=strtok(NULL, " \t")) {
+ res = cstr_to_double(tk, tmp + i);
+ if(res != RES_OK) {
+ fprintf(stderr, "%s:%lu: cannot read the histogram value `%s'.\n",
+ filename, (unsigned long)iline, tk);
+ return res;
+ }
+ }
+ if(i < 2) {
+ fprintf(stderr,
+ "%s:%lu: wrong histogram value `%s'. Expect <VALUE> <PROBABILITY>\n",
+ filename, (unsigned long)iline, str);
+ return RES_BAD_ARG;
+ }
+
+ if(tk) {
+ fprintf(stderr, "%s:%lu: wrong histogram entry `%s'.\n",
+ filename, (unsigned long)iline, str);
+ }
+ entry->value = tmp[0];
+ entry->cdf = tmp[1];
+ return RES_OK;
+}
+
+static INLINE int
+cmp_entry(const void* a, const void* b)
+{
+ double val;
+ const struct schiff_hentry* entry = b;
+ ASSERT(a && b);
+ val = *(const double*)a;
+ return val < entry->cdf ? -1 : val > entry->cdf ? 1 : 0;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+schiff_histogram_load(struct darray_hentry* histo, const char* filename)
+{
+ FILE* fp = NULL;
+ res_T res = RES_OK;
+ ASSERT(histo && filename);
+
+ fp = fopen(filename, "r");
+ if(!fp) {
+ fprintf(stderr, "Cannot open the file of optical properties `%s'.\n",
+ filename);
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ res = schiff_histogram_load_stream(histo, fp, filename);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(fp) fclose(fp);
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+schiff_histogram_load_stream
+ (struct darray_hentry* histo,
+ FILE* stream,
+ const char* stream_name)
+{
+ struct schiff_streamline streamline;
+ size_t iline;
+ char* line;
+ double cdf = 0.0;
+ res_T res = RES_OK;
+ ASSERT(histo && stream_name && stream_name);
+
+ darray_hentry_clear(histo);
+ schiff_streamline_init(&streamline);
+
+ iline = 1;
+ while((res = schiff_streamline_read(&streamline, stream, &line)) != RES_EOF) {
+ struct schiff_hentry entry;
+
+ if(*line == '\0' /*Empty line*/ || *line == '#' /* Comment */)
+ continue;
+
+ res = parse_schiff_hentry(&entry, stream_name, iline, strtok(line, "#"));
+ if(res == RES_OK) { /* *NO* error */
+ cdf += entry.cdf;
+ entry.cdf = cdf;
+ res = darray_hentry_push_back(histo, &entry);
+ if(res != RES_OK) {
+ fprintf(stderr, "%s:%lu: couldn't store the histogram entry.\n",
+ stream_name, (unsigned long)iline);
+ goto error;
+ }
+ }
+ ++iline;
+ }
+ if(res != RES_EOF) goto error;
+ res = RES_OK;
+ if(!eq_eps(cdf, 1, 1.e-6)) {
+ fprintf(stderr, "%s: unormalized histogram.\n", stream_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+
+exit:
+ schiff_streamline_release(&streamline);
+ return res;
+error:
+ darray_hentry_clear(histo);
+ goto exit;
+}
+
+double
+schiff_histogram_sample(const struct darray_hentry* histo, const double u)
+{
+ const struct schiff_hentry* entries;
+ const struct schiff_hentry* find;
+ size_t nentries;
+ ASSERT(histo && u >= 0.0 && u < 1.0);
+
+ entries = darray_hentry_cdata_get(histo);
+ nentries = darray_hentry_size_get(histo);
+ ASSERT(nentries);
+
+ find = search_lower_bound
+ (&u, entries, nentries, sizeof(struct schiff_hentry), cmp_entry);
+ if(!find) { /* Numerical issue */
+ return entries[nentries - 1].value;
+ } else {
+ return find->value;
+ }
+}
+
diff --git a/src/schiff_histogram.h b/src/schiff_histogram.h
@@ -0,0 +1,56 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#ifndef SCHIFF_HISTOGRAM_H
+#define SCHIFF_HISTOGRAM_H
+
+#include <rsys/dynamic_array.h>
+
+struct schiff_hentry { double value, cdf; };
+#define DARRAY_NAME hentry
+#define DARRAY_DATA struct schiff_hentry
+#include <rsys/dynamic_array.h>
+
+extern LOCAL_SYM res_T
+schiff_histogram_load
+ (struct darray_hentry* histo,
+ const char* filename);
+
+extern LOCAL_SYM res_T
+schiff_histogram_load_stream
+ (struct darray_hentry* histo,
+ FILE* stream,
+ const char* stream_name);
+
+extern LOCAL_SYM double
+schiff_histogram_sample
+ (const struct darray_hentry* histo,
+ const double u); /* in [0, 1[ */
+
+#endif /* SCHIFF_HISTOGRAM_H */
+
diff --git a/src/schiff_streambuf.h b/src/schiff_streambuf.h
@@ -0,0 +1,91 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#ifndef SCHIFF_STREAMBUF_H
+#define SHCIFF_STREAMBUF_H
+
+#include <rsys/stretchy_array.h>
+#include <string.h>
+
+struct schiff_streambuf { char* buf; };
+
+res_T
+schiff_streambuf_init(struct schiff_streambuf* sbuf)
+{
+ ASSERT(sbuf);
+ memset(sbuf, 0, sizeof(struct schiff_sbuf));
+ sbuf->buf = sa_add(buf, 128);
+}
+
+res_T
+schiff_streambuf_release(struct schiff_streambuf* sbuf)
+{
+ ASSERT(sbuf);
+ sa_release(sbuf->buf);
+}
+
+res_T
+schiff_streambuf_read_line
+ (struct schiff_streambuf* sbuf,
+ FILE* stream,
+ char** out_line)
+{
+ char* line = NULL;
+ size_t last_char;
+ const size_t chunk = 128;
+ res_T res = RES_OK;
+ ASSERT(sbuf && stream && out_line);
+
+ if(!fgets(sbuf->buf, (int)sa_size(sbuf->buf), stream)) {
+ res = RES_EOF;
+ goto exit;
+ }
+
+ while(!strrchr(buf, '\n')) { /* Ensure that the whole line is read */
+ if(!fgets(sa_add(sbuf->buf, buf_chunk), buf_chunk, stream)) /* EOF */
+ break;
+ }
+
+ /* Remove leading spaces */
+ line = sbuf->buf;
+ while((*line == ' ' || *line == '\t') && *line != '\0') ++line;
+
+ /* Remove newline character(s) */
+ last_char = strlen(line);
+ while(last_char-- && (line[last_char]=='\n' || line[last_char]=='\r'));
+ line[last_char + 1] = '\0';
+exit:
+ *out_line = line;
+ return res;
+error:
+ line = NULL;
+ goto error;
+}
+
+#endif /* SCHIFF_STREAMBUF_H */
+
