commit 191dee04476221ae1f25a4878895af71b5437202
parent 29430c46f57e0db91bf0c09c2506aeb65bb39696
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 19 Oct 2015 14:48:12 +0200
Update the Schiff integration API
Diffstat:
4 files changed, 84 insertions(+), 83 deletions(-)
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -49,21 +49,23 @@
#define SSCHIFF(Func) sschiff_ ## Func
#endif
-/* Forward declaration of external types */
+/* Forward declaration of external types. */
struct mem_allocator;
struct logger;
struct s3d_device;
struct s3d_shape;
struct ssp_rng;
+/* Optical properties of a material handled by the Schiff integrator. */
struct sschiff_material_properties {
double medium_refractive_index;
double relative_real_refractive_index;
double relative_imaginary_refractive_index;
};
+/* Client side material used by the Schiff integrator */
struct sschiff_material {
- /* Retrieve the refractive complex index of the material */
+ /* Retrieve the optical properties of the material */
void (*get_property)
(void* material,
const double wavelength,
@@ -73,23 +75,17 @@ struct sschiff_material {
static const struct sschiff_material SSCHIFF_NULL_MATERIAL = { NULL, NULL };
-/* Descriptor of an integration */
-struct sschiff_integrator_desc {
- res_T (*sample_geometry) /* Sample a geometry i.e. a shape and its material */
+struct sschiff_geometry_distribution {
+ res_T (*sample) /* Sample a geometry i.e. a shape and its material */
(struct ssp_rng* rng,
struct sschiff_material* mtl, /* Sampled material */
struct s3d_shape* shape, /* Sampled shape */
- void* sampler_context);
- void* sampler_context;
- double wavelength; /* Wavelenght to estimate */
- size_t scattering_angles_count; /* # of scattering angles to handle */
- size_t sampled_geometries_count; /* # geometries to sample */
- size_t sampled_directions_count; /* # directions to sample per geometry */
+ void* context);
+ void* context;
};
-static const struct sschiff_integrator_desc SSCHIFF_NULL_INTEGRATOR_DESC = {
- NULL, NULL, 0, 0, 0, 0
-};
+static const struct sschiff_geometry_distribution
+SSCHIFF_NULL_GEOMETRY_DISTRIBUTION = { NULL, NULL };
/* Result of the Monte Carlo estimation */
struct sschiff_estimator_status {
@@ -128,7 +124,11 @@ SSCHIFF_API res_T
sschiff_integrate
(struct sschiff_device* dev,
struct ssp_rng* rng,
- struct sschiff_integrator_desc* desc,
+ struct sschiff_geometry_distribution* distribution,
+ const double wavelength, /* Wavelenght to estimate */
+ const size_t scattering_angles_count, /* # of scattering angles to handle */
+ const size_t sampled_geometries_count, /* # geometries to sample */
+ const size_t sampled_directions_count, /* # directions to sample per geometry */
struct sschiff_estimator** estimator);
SSCHIFF_API res_T
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -26,7 +26,7 @@
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL license and that you accept its terms. */
-#define _POSIX_C_SOURCE 200112L
+#define _POSIX_C_SOURCE 200112L /* nextafterf support */
#include "sschiff.h"
#include "sschiff_device.h"
@@ -44,8 +44,7 @@
#include <math.h>
-#define NDIRS 1000 /* # Sampled directions */
-
+/* Type of computed Monte Carlo weights */
enum weight_type {
WEIGHT_EXTINCTION_CROSS_SECTION,
WEIGHT_ABSORPTION_CROSS_SECTION,
@@ -54,6 +53,7 @@ enum weight_type {
WEIGHTS_COUNT
};
+/* Accumulator of Monte Carlo weights */
struct accum {
double weight;
double sqr_weight;
@@ -82,7 +82,7 @@ morton2D_decode(const uint32_t u32)
return (uint16_t)x;
}
-/* Compute an orthonormal basis where `dir' is the Z axis */
+/* Compute an orthonormal basis where `dir' is the Z axis. */
static void
build_basis(const float dir[3], float basis[9])
{
@@ -93,7 +93,7 @@ build_basis(const float dir[3], float basis[9])
/* Define which axis direction is minimal and use its unit vector to compute
* a vector ortogonal to `dir'. This ensures that the unit vector and `dir'
- * are not colinear and thus that their cross product is not a zero vector */
+ * are not colinear and thus that their cross product is not a zero vector. */
FOR_EACH(i, 0, 3) dir_abs[i] = absf(dir[i]);
if(dir_abs[0] < dir_abs[1]) {
if(dir_abs[0] < dir_abs[2]) f3(axis_min, 1.f, 0.f, 0.f);
@@ -340,8 +340,9 @@ static res_T
radiative_properties
(struct sschiff_device* dev,
struct ssp_rng* rng,
- struct sschiff_integrator_desc* desc,
const int istep,
+ const double wavelength,
+ const size_t ndirs,
struct s3d_scene* scene,
struct sschiff_material_properties* props,
struct accum accums[WEIGHTS_COUNT])
@@ -352,7 +353,7 @@ radiative_properties
float centroid[3];
size_t idir;
int i;
- ASSERT(dev && rng && desc && props && accums);
+ ASSERT(dev && rng && ndirs && props && accums);
(void)istep;
S3D(scene_get_aabb(scene, lower, upper));
@@ -376,7 +377,7 @@ radiative_properties
f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f);
memset(weights_dirs, 0, sizeof(weights_dirs));
- FOR_EACH(idir, 0, desc->sampled_directions_count) {
+ FOR_EACH(idir, 0, ndirs) {
double weights[WEIGHTS_COUNT];
float dir[4];
float basis[9];
@@ -402,8 +403,8 @@ radiative_properties
f3_max(upper, upper, pt[i]);
}
- compute_monte_carlo_weight(dev, scene, rng, props,
- desc->wavelength, basis, centroid, lower, upper, weights);
+ compute_monte_carlo_weight(dev, scene, rng, props, wavelength, basis,
+ centroid, lower, upper, weights);
FOR_EACH(i, 0, WEIGHTS_COUNT)
weights_dirs[i] += weights[i];
@@ -417,7 +418,7 @@ radiative_properties
#endif
}
FOR_EACH(i, 0, WEIGHTS_COUNT) {
- weights_dirs[i] /= (double)desc->sampled_directions_count;
+ weights_dirs[i] /= (double)ndirs;
accums[i].weight += weights_dirs[i];
accums[i].sqr_weight += weights_dirs[i] * weights_dirs[i];
accums[i].naccums += 1;
@@ -437,14 +438,10 @@ get_status(const struct accum* acc, struct sschiff_estimator_status* status)
}
static char
-check_desc(struct sschiff_integrator_desc* desc)
+check_distribution(struct sschiff_geometry_distribution* distrib)
{
- ASSERT(desc);
- return desc->sample_geometry
- && desc->scattering_angles_count
- && desc->sampled_geometries_count
- && desc->sampled_directions_count
- && desc->wavelength > 0.f;
+ ASSERT(distrib);
+ return distrib->sample != NULL;
}
static void
@@ -496,7 +493,11 @@ res_T
sschiff_integrate
(struct sschiff_device* dev,
struct ssp_rng* rng,
- struct sschiff_integrator_desc* desc,
+ struct sschiff_geometry_distribution* distrib,
+ const double wavelength,
+ const size_t scattering_angles_count,
+ const size_t sampled_geometries_count,
+ const size_t sampled_directions_count,
struct sschiff_estimator** out_estimator)
{
struct sschiff_estimator* estimator = NULL;
@@ -504,15 +505,16 @@ sschiff_integrate
struct s3d_scene* scene = NULL;
size_t i;
res_T res = RES_OK;
+ (void)scattering_angles_count;
- if(!dev || !rng || !desc || !check_desc(desc) || !out_estimator) {
+ if(!dev || !rng || !distrib || !check_distribution(distrib) || !out_estimator) {
res = RES_BAD_ARG;
goto error;
}
res = estimator_create(dev, &estimator);
if(res != RES_OK) goto error;
- estimator->wavelength = desc->wavelength;
+ estimator->wavelength = wavelength;
res = s3d_shape_create_mesh(dev->s3d, &shape);
if(res != RES_OK) {
@@ -532,21 +534,21 @@ sschiff_integrate
}
memset(estimator->accums, 0, sizeof(struct accum)*WEIGHTS_COUNT);
- FOR_EACH(i, 0, desc->sampled_geometries_count) {
+ FOR_EACH(i, 0, sampled_geometries_count) {
struct sschiff_material material = SSCHIFF_NULL_MATERIAL;
struct sschiff_material_properties props;
- res = desc->sample_geometry(rng, &material, shape, desc->sampler_context);
+ res = distrib->sample(rng, &material, shape, distrib->context);
if(res != RES_OK) {
log_error(dev, "Error in sampling a Schiff geometry.\n");
goto error;
}
- material.get_property(material.material, desc->wavelength, &props);
+ material.get_property(material.material, wavelength, &props);
S3D(scene_begin_session(scene, S3D_TRACE));
- res = radiative_properties(dev, rng, desc, (int)i, scene, &props,
- estimator->accums);
+ res = radiative_properties(dev, rng, (int)i, wavelength,
+ sampled_directions_count, scene, &props, estimator->accums);
if(res != RES_OK) goto error;
S3D(scene_end_session(scene));
diff --git a/src/test_sschiff_estimator_cylinder.c b/src/test_sschiff_estimator_cylinder.c
@@ -806,7 +806,7 @@ main(int argc, char** argv)
struct mem_allocator allocator;
struct sampler_context sampler_ctx;
struct sschiff_device* dev;
- struct sschiff_integrator_desc desc = SSCHIFF_NULL_INTEGRATOR_DESC;
+ struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
struct sschiff_estimator* estimator;
struct sschiff_estimator_status status;
struct ssp_rng* rng;
@@ -883,24 +883,25 @@ main(int argc, char** argv)
sampler_ctx.aspect_ratio = 0.263;
sampler_ctx.mean_radius = 0.983;
- desc.sample_geometry = sample_cylinder;
- desc.sampler_context = &sampler_ctx;
- desc.scattering_angles_count = 1;
- desc.sampled_geometries_count = 300;
- desc.sampled_directions_count = 30;
+ distrib.sample = sample_cylinder;
+ distrib.context = &sampler_ctx;
FOR_EACH(i, 0, nresults) {
+ const size_t nscatt_angles = 1;
+ const size_t ngeoms = 300;
+ const size_t ndirs = 30;
+ const double wavelength = results[i].wavelength;
double dst;
- desc.wavelength = results[i].wavelength;
time_current(&t0);
- CHECK(sschiff_integrate(dev, rng, &desc, &estimator), RES_OK);
+ CHECK(sschiff_integrate(dev, rng, &distrib, wavelength, nscatt_angles,
+ ngeoms, ndirs, &estimator), RES_OK);
time_current(&t1);
time_sub(&t0, &t1, &t0);
time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
- CHECK(status.nsteps, desc.sampled_geometries_count);
+ CHECK(status.nsteps, ngeoms);
dst = status.E - results[i].extinction_cross_section;
printf(
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -241,21 +241,22 @@ check_schiff_estimation
{
char buf[64];
struct sschiff_estimator* estimator;
- struct sschiff_integrator_desc desc = SSCHIFF_NULL_INTEGRATOR_DESC;
+ struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
struct sschiff_estimator_status status;
struct time t0, t1;
size_t i;
+ const size_t nscatt_angles = 1;
+ const size_t ngeoms = 100;
+ const size_t ndirs = 10;
+ const double wavelength = sampler_ctx->wavelength;
+
NCHECK(sampler_ctx, NULL);
NCHECK(results, NULL);
NCHECK(results, 0);
- desc.sample_geometry = sample_sphere;
- desc.scattering_angles_count = 1;
- desc.sampler_context = sampler_ctx;
- desc.sampled_geometries_count = 100;
- desc.sampled_directions_count = 10;
- desc.wavelength = sampler_ctx->wavelength;
+ distrib.sample = sample_sphere;
+ distrib.context = sampler_ctx;
printf("Schiff sphere estimation - m: %g; n_r: %g; kappa_r: %g; n_e: %g\n",
sampler_ctx->wavelength,
@@ -268,14 +269,15 @@ check_schiff_estimation
sampler_ctx->mean_radius = results[i].mean_radius;
time_current(&t0);
- CHECK(sschiff_integrate(dev, rng, &desc, &estimator), RES_OK);
+ CHECK(sschiff_integrate(dev, rng, &distrib, wavelength, nscatt_angles,
+ ngeoms, ndirs, &estimator), RES_OK);
time_current(&t1);
time_sub(&t0, &t1, &t0);
time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
printf("%d - mean radius: %5.2fe-6 - %s: \n", (int)i, results[i].mean_radius, buf);
CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
- CHECK(status.nsteps, desc.sampled_geometries_count);
+ CHECK(status.nsteps, ngeoms);
dst = status.E - results[i].extinction_cross_section;
printf(" Extinction cross section = %7.2g ~ %12.7g +/- %12.7g (%6.2g)\n",
results[i].extinction_cross_section, status.E, status.SE, dst / status.SE);
@@ -310,7 +312,7 @@ main(int argc, char** argv)
struct mem_allocator allocator;
struct sampler_context sampler_ctx;
struct sschiff_device* dev;
- struct sschiff_integrator_desc desc = SSCHIFF_NULL_INTEGRATOR_DESC;
+ struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
struct sschiff_estimator* estimator;
struct sschiff_estimator_status status;
struct ssp_rng* rng;
@@ -354,29 +356,25 @@ main(int argc, char** argv)
sampler_ctx.mean_radius = 1.0;
sampler_ctx.sigma = 1.18;
- desc.sample_geometry = sample_sphere;
- desc.sampler_context = &sampler_ctx;
- desc.wavelength = sampler_ctx.wavelength;
- desc.scattering_angles_count = 1;
- desc.sampled_geometries_count = 1;
- desc.sampled_directions_count = 1;
-
- CHECK(sschiff_integrate(NULL, NULL, NULL, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, NULL, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, NULL, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, NULL, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, &desc, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, &desc, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, &desc, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, &desc, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, NULL, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, NULL, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, NULL, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, NULL, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, &desc, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, &desc, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, &desc, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, &desc, &estimator), RES_OK);
+ distrib.sample = sample_sphere;
+ distrib.context = &sampler_ctx;
+
+ CHECK(sschiff_integrate(NULL, NULL, NULL, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, 0.6e-6, 1, 1, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, NULL, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, NULL, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, NULL, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, NULL, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, NULL, &distrib, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, NULL, &distrib, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(NULL, rng, &distrib, 0.6e-6, 1, 1, 1, &estimator), RES_BAD_ARG);
+ CHECK(sschiff_integrate(dev, rng, &distrib, 0.6e-6, 1, 1, 1, &estimator), RES_OK);
CHECK(sschiff_estimator_get_extinction_cross_section(NULL, NULL), RES_BAD_ARG);
CHECK(sschiff_estimator_get_extinction_cross_section(estimator, NULL), RES_BAD_ARG);
