commit 90e3f632604097d50432ae0ed1d03ea212cdf9ba
parent 19c0d487c13b6bc3e504d0fed5776e30aee82d79
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 6 Nov 2015 13:53:03 +0100
Add the supershape distribution
Add the -u option that defines the lower/upper bounds of the super shape
parameters.
Diffstat:
4 files changed, 268 insertions(+), 29 deletions(-)
diff --git a/src/schiff_args.c b/src/schiff_args.c
@@ -31,6 +31,7 @@
#include "schiff_args.h"
#include "schiff_optical_properties.h"
+#include <rsys/dynamic_array_char.h>
#include <rsys/cstr.h>
#include <rsys/stretchy_array.h>
@@ -52,11 +53,11 @@ print_help(const char* binary)
"Method for Short Wavelength or High Energy Scattering\" (L. Schiff, 1956).\n\n",
binary);
printf(
-"FILE lists the per wavelength optical properties of the micro organisms.\n"
-"Each line must be formatted as \"W Nr Kr Ne\" where \"W\" is the wavelength\n"
-"in micron, \"Nr\" and \"Kr\" are the real and imaginary parts, respectively, of\n"
-"the relative refractive index, and \"Ne\" the refractive index of the medium.\n"
-"With no FILE, read optical properties from standard input.\n\n");
+"FILE lists the per wavelength optical properties of the micro organisms. 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(
" -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"
@@ -90,7 +91,8 @@ print_help(const char* binary)
" r = 1/(ln(S)*x*sqrt(2PI)) * exp(-(ln(x)-ln(R))^2 / (2*ln(S)^2))\n",
SCHIFF_SPHERE_DEFAULT.nslices);
printf(
-" -w A[:B]... list of wavelengths (in micron) to integrate.\n");
+" -w A[:B]... list of wavelengths in vacuum (expressed in micron) to\n"
+" integrate.\n");
}
static res_T
@@ -186,6 +188,125 @@ error:
}
static res_T
+parse_super_shape(const char* str, double formulas[2][6])
+{
+ struct darray_char tmp;
+ char* tk;
+ char* tk_ctx;
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(str && formulas);
+
+ darray_char_init(&mem_default_allocator, &tmp);
+
+ /* Locallly copy the input string */
+ res = darray_char_resize(&tmp, strlen(str) + 1 /*NULL char*/);
+ if(res != RES_OK) goto error;
+ strcpy(darray_char_data_get(&tmp), str);
+
+ tk = strtok_r(darray_char_data_get(&tmp), ",", &tk_ctx);
+ if(!tk) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = cstr_to_list_double(tk, formulas[0], &len, 6);
+ if(res != RES_OK) goto error;
+ if(len != 6) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ tk = strtok_r(NULL, ",", &tk_ctx);
+ if(!tk) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = cstr_to_list_double(tk, formulas[1], &len, 6);
+ if(res != RES_OK) goto error;
+ if(len != 6) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ darray_char_release(&tmp);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_super_shape_distribution(const char* str, struct schiff_args* args)
+{
+ struct darray_char tmp;
+ struct schiff_super_shape* sshape;
+ char* tk;
+ char* tk_ctx;
+ size_t str_len;
+ int i;
+ res_T res = RES_OK;
+
+ darray_char_init(&mem_default_allocator, &tmp);
+
+ if(args->geom.type != SCHIFF_NONE && args->geom.type != SCHIFF_SUPER_SHAPE) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Locallly copy the input string */
+ str_len = strlen(str);
+ res = darray_char_resize(&tmp, str_len + 1 /*NULL char*/);
+ if(res != RES_OK) goto error;
+ strcpy(darray_char_data_get(&tmp), str);
+
+ /* Parse "lower" super shape */
+ tk = strtok_r(darray_char_data_get(&tmp), "#", &tk_ctx);
+ if(!tk) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ res = parse_super_shape(tk, args->geom.data.super_shape.lower);
+ if(res != RES_OK) goto error;
+
+ /* Parse "upper" super shape */
+ tk = strtok_r(NULL, "#", &tk_ctx);
+ if(!tk) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ res = parse_super_shape(tk, args->geom.data.super_shape.upper);
+ if(res != RES_OK) goto error;
+
+ /* Parse the optionnal "nslices" attrib */
+ tk = strtok_r(NULL, "#", &tk_ctx);
+ if(!tk) {
+ args->geom.data.super_shape.nslices = SCHIFF_SUPER_SHAPE_DEFAULT.nslices;
+ } else {
+ res = cstr_to_uint(tk, &args->geom.data.super_shape.nslices);
+ if(res != RES_OK) goto error;
+ }
+
+ args->geom.type = SCHIFF_SUPER_SHAPE;
+
+ sshape = &args->geom.data.super_shape;
+ FOR_EACH(i, 0, 6) {
+ if(sshape->lower[0][i] > sshape->upper[0][i]
+ || sshape->lower[1][i] > sshape->upper[1][i]) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+exit:
+ darray_char_release(&tmp);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
parse_wavelengths(const char* str, struct schiff_args* args)
{
size_t len;
@@ -232,7 +353,7 @@ schiff_args_init
ASSERT(argc && argv && args);
*args = SCHIFF_ARGS_NULL;
- while((opt = getopt(argc, argv, "c:d:g:G:ho:s:w:")) != -1) {
+ while((opt = getopt(argc, argv, "c:d:g:G:ho:s:u:w:")) != -1) {
res_T res = RES_OK;
switch(opt) {
case 'c': res = parse_cylinder_distribution(optarg, args); break;
@@ -245,6 +366,7 @@ schiff_args_init
return RES_OK;
case 'o': args->output_filename = optarg; break;
case 's': res = parse_sphere_distribution(optarg, args); break;
+ case 'u': res = parse_super_shape_distribution(optarg, args); break;
case 'w': res = parse_wavelengths(optarg, args); break;
default: res = RES_BAD_ARG; break;
}
diff --git a/src/schiff_geometry.c b/src/schiff_geometry.c
@@ -46,9 +46,8 @@ struct cylinder_context {
};
/* 3D Context of the super shape geometry */
-enum { A, B, M, N0, N1, N2 };
-struct supershape_context {
- double super_formulas[2][6];
+struct super_shape_context {
+ double formulas[2][6];
};
/* Distribution context of a sphere geometry */
@@ -65,9 +64,9 @@ struct cylinder_distribution_context {
};
/* Distribtion context of a super shape geometry */
-struct supershape_distribution_context {
- double lower[2][6]; /* Lower bound of the super formula parameters */
- double upper[2][6]; /* Upper bound of the super formula parameters */
+struct super_shape_distribution_context {
+ double lower[2][6];
+ double upper[2][6];
};
/* 3D context of a generic geometry */
@@ -77,7 +76,7 @@ struct geometry_context {
union {
struct cylinder_context cylinder;
struct sphere_context sphere;
- struct supershape_context supershape;
+ struct super_shape_context super_shape;
} data;
};
@@ -89,7 +88,7 @@ struct geometry_distribution_context {
union {
struct cylinder_distribution_context cylinder;
struct sphere_distribution_context sphere;
- struct supershape_distribution_context supershape;
+ struct super_shape_distribution_context super_shape;
} data;
};
@@ -136,6 +135,78 @@ sphere_get_position(const unsigned ivert, float vertex[3], void* ctx)
vertex[2] *= geom->data.sphere.radius;
}
+static void
+super_shape_get_position(const unsigned ivert, float vertex[3], void* ctx)
+{
+ struct geometry_context* geom = ctx;
+ float angles[2];
+ double cos_angles[2];
+ double sin_angles[2];
+ double uv[2];
+ int iform;
+ ASSERT(geom && geom->type == SCHIFF_SUPER_SHAPE);
+
+ schiff_mesh_get_polar_position(geom->mesh, ivert, angles);
+
+ FOR_EACH(iform, 0, 2) {
+ double m, k, g;
+ double* form = geom->data.super_shape.formulas[iform];
+
+ m = cos(form[M] * angles[iform] / 4.0);
+ m = fabs(m) / form[A];
+ k = sin(form[M] * angles[iform] / 4.0);
+ k = fabs(k) / form[B];
+ g = pow(m, form[N1]) + pow(k, form[N2]);
+ uv[iform] = pow(g, (-1.0/form[N0]));
+
+ cos_angles[iform] = cos(angles[iform]);
+ sin_angles[iform] = sin(angles[iform]);
+ }
+
+ vertex[0] = (float)(uv[0] * cos_angles[0] * uv[1] * cos_angles[1]);
+ vertex[1] = (float)(uv[0] * sin_angles[0] * uv[1] * cos_angles[1]);
+ vertex[2] = (float)(uv[1] * sin_angles[1]);
+
+}
+
+static res_T
+geometry_sample_cylinder
+ (struct ssp_rng* rng,
+ struct sschiff_material* mtl,
+ struct s3d_shape* shape,
+ void* context)
+{
+ struct geometry_distribution_context* distrib = context;
+ struct geometry_context geom;
+ struct s3d_vertex_data attrib;
+ size_t nverts, nprims;
+ double y;
+ 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);
+ 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);
+
+ attrib.usage = S3D_POSITION;
+ 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) / 3/*#coords*/;
+ nprims = darray_uint_size_get(&distrib->mesh->indices) / 3/*#indices per prim*/;
+
+ return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims,
+ geometry_get_indices, (unsigned)nverts, &attrib, 1, &geom);
+}
+
static res_T
geometry_sample_sphere
(struct ssp_rng* rng,
@@ -170,7 +241,7 @@ geometry_sample_sphere
}
static res_T
-geometry_sample_cylinder
+geometry_sample_super_shape
(struct ssp_rng* rng,
struct sschiff_material* mtl,
struct s3d_shape* shape,
@@ -179,29 +250,32 @@ geometry_sample_cylinder
struct geometry_distribution_context* distrib = context;
struct geometry_context geom;
struct s3d_vertex_data attrib;
+ const struct super_shape_distribution_context* sshape;
size_t nverts, nprims;
- double y;
+ int i;
ASSERT(rng && mtl && shape && context);
- ASSERT(distrib->type == SCHIFF_CYLINDER);
+ ASSERT(distrib->type == SCHIFF_SUPER_SHAPE);
- y = ssp_ran_lognormal
- (rng, distrib->data.cylinder.log_mean_radius, distrib->data.cylinder.log_sigma);
- geom.type = SCHIFF_CYLINDER;
+ sshape = &distrib->data.super_shape;
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);
+ geom.type = SCHIFF_SUPER_SHAPE;
+
+ FOR_EACH(i, 0, 6) {
+ geom.data.super_shape.formulas[0][i] = (float)
+ ssp_rng_uniform_double(rng, sshape->lower[0][i], sshape->upper[0][i]);
+ geom.data.super_shape.formulas[1][i] = (float)
+ ssp_rng_uniform_double(rng, sshape->lower[1][i], sshape->upper[1][i]);
+ }
attrib.usage = S3D_POSITION;
attrib.type = S3D_FLOAT3;
- attrib.get = cylinder_get_position;
+ attrib.get = super_shape_get_position;
mtl->get_property = get_material_property;
mtl->material = distrib->properties;
- nverts = darray_float_size_get(&distrib->mesh->vertices) / 3/*#coords*/;
- nprims = darray_uint_size_get(&distrib->mesh->indices) / 3/*#indices per prim*/;
+ nverts = darray_float_size_get(&distrib->mesh->vertices) / 2/*#theta/phi*/;
+ nprims = darray_uint_size_get(&distrib->mesh->indices) / 3/*#indices*/;
return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims,
geometry_get_indices, (unsigned)nverts, &attrib, 1, &geom);
@@ -217,6 +291,7 @@ schiff_geometry_distribution_init
struct schiff_optical_properties* properties)
{
struct geometry_distribution_context* ctx = NULL;
+ int i;
res_T res = RES_OK;
ASSERT(distrib && geom);
@@ -252,6 +327,17 @@ schiff_geometry_distribution_init
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);
+ FOR_EACH(i, 0, 6) {
+ ctx->data.super_shape.lower[0][i] = geom->data.super_shape.lower[0][i];
+ ctx->data.super_shape.lower[1][i] = geom->data.super_shape.lower[1][i];
+ ctx->data.super_shape.upper[0][i] = geom->data.super_shape.upper[0][i];
+ ctx->data.super_shape.upper[1][i] = geom->data.super_shape.upper[1][i];
+ }
+ distrib->sample = geometry_sample_super_shape;
+ break;
default: FATAL("Unreachable code.\n"); break;
}
if(res != RES_OK) {
diff --git a/src/schiff_geometry.h b/src/schiff_geometry.h
@@ -34,6 +34,7 @@
enum schiff_geometry_type {
SCHIFF_CYLINDER,
SCHIFF_SPHERE,
+ SCHIFF_SUPER_SHAPE,
SCHIFF_NONE
};
@@ -58,11 +59,27 @@ struct schiff_cylinder {
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[2][6]; /* Lower bound of the formula params */
+ double upper[2][6]; /* Upper bound of the formula params */
+ unsigned nslices;
+};
+
+#define SCHIFF_SUPER_SHAPE_DEFAULT__ { \
+ {{ 1, 1, 2, 1, 1, 1 }, { 1, 1, 5, 1, 1, 3 }}, \
+ {{ 1, 1, 7, 1, 1, 2 }, { 1, 1, 7, 1, 1, 7 }}, \
+ 64 \
+}
+static const struct schiff_super_shape SCHIFF_SUPER_SHAPE_DEFAULT =
+ SCHIFF_SUPER_SHAPE_DEFAULT__;
+
struct schiff_geometry {
enum schiff_geometry_type type;
union {
struct schiff_cylinder cylinder;
struct schiff_sphere sphere;
+ struct schiff_super_shape super_shape;
} data;
};
diff --git a/src/schiff_mesh.h b/src/schiff_mesh.h
@@ -84,7 +84,7 @@ static INLINE void
schiff_mesh_get_cartesian_position
(struct schiff_mesh* mesh,
const unsigned ivert,
- float vertex[])
+ float vertex[3])
{
const size_t i = ivert * 3;
ASSERT(mesh && vertex && mesh->coordinates == SCHIFF_CARTESIAN);
@@ -95,5 +95,19 @@ schiff_mesh_get_cartesian_position
vertex[2] = darray_float_data_get(&mesh->vertices)[i + 2];
}
+static INLINE void
+schiff_mesh_get_polar_position
+ (struct schiff_mesh* mesh,
+ const unsigned ivert,
+ float angles[2]) /* theta, phi */
+{
+ const size_t i = ivert * 2;
+ ASSERT(mesh && angles && mesh->coordinates == SCHIFF_POLAR);
+ ASSERT(darray_float_size_get(&mesh->vertices) % 2 == 0);
+ ASSERT(ivert < darray_float_size_get(&mesh->vertices) / 2);
+ angles[0] = darray_float_data_get(&mesh->vertices)[i + 0];
+ angles[1] = darray_float_data_get(&mesh->vertices)[i + 1];
+}
+
#endif /* SBOX_SCHIFF_MESH_H */
