commit f383438b2f5de3c9e2fe98be37ed691c6fe452eb
parent eeee76acf3d587adfdcadd36c08444def80b684a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 13 Oct 2015 15:25:35 +0200
Update the schiff_integrate API
One can control the number of sampled geometries and the number of
sampled direction per geometry.
Diffstat:
3 files changed, 66 insertions(+), 56 deletions(-)
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -73,6 +73,7 @@ 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 ssp_rng* rng,
@@ -80,11 +81,14 @@ struct sschiff_integrator_desc {
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 */
};
static const struct sschiff_integrator_desc SSCHIFF_NULL_INTEGRATOR_DESC = {
- NULL, NULL, 0
+ NULL, NULL, 0, 0, 0, 0
};
/* Result of the Monte Carlo estimation */
@@ -92,7 +96,7 @@ struct sschiff_estimator_status {
double E; /* Expected value */
double V; /* Variance */
double SE; /* Standard error */
- size_t nsteps;
+ size_t nsteps; /* # Monte Carlo weights */
};
/* Opaque types */
@@ -125,8 +129,6 @@ sschiff_integrate
(struct sschiff_device* dev,
struct ssp_rng* rng,
struct sschiff_integrator_desc* desc,
- const double wavelength, /* Wavelength to integrate */
- const size_t steps_count, /* #realisations */
struct sschiff_estimator** estimator);
SSCHIFF_API res_T
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -44,7 +44,7 @@
#include <math.h>
-#define NDIRS 1 /* # Sampled directions */
+#define NDIRS 1000 /* # Sampled directions */
struct accum {
double weight;
@@ -53,7 +53,6 @@ struct accum {
};
struct sschiff_estimator {
- size_t nsteps;
double wavelength;
struct accum extinction_cross_section;
@@ -125,7 +124,8 @@ build_transform
static res_T
compute_path_length
- (struct s3d_scene* scn,
+ (struct sschiff_device* dev,
+ struct s3d_scene* scn,
const struct s3d_hit* first_hit,
const float ray_org[3],
const float ray_dir[3],
@@ -145,7 +145,10 @@ compute_path_length
S3D(scene_trace_ray(scn, ray_org, ray_dir, range, &hit));
} while(S3D_PRIMITIVE_EQ(&hit.prim, &prim_from)); /* Self-intersection */
- if(S3D_HIT_NONE(&hit)) return RES_BAD_ARG;
+ if(S3D_HIT_NONE(&hit)) {
+ log_error(dev, "Error in computing the radiative path length.\n");
+ return RES_BAD_ARG;
+ }
len += hit.distance - range[0];
range[0] = hit.distance;
@@ -167,7 +170,8 @@ compute_path_length
static res_T
compute_monte_carlo_weight
- (struct s3d_scene* scn,
+ (struct sschiff_device* dev,
+ struct s3d_scene* scn,
struct ssp_rng* rng,
struct sschiff_material* mtl,
const double wavelength,
@@ -179,7 +183,6 @@ compute_monte_carlo_weight
{
struct s3d_hit hit;
double sample[2];
- double weight;
const float range[2] = { 0, FLT_MAX };
float axis_x[3], axis_y[3], axis_z[3];
float plane_size[2];
@@ -210,15 +213,14 @@ compute_monte_carlo_weight
f3_add(ray_org, f3_add(ray_org, x, y), org);
S3D(scene_trace_ray(scn, ray_org, axis_z, range, &hit));
- if(S3D_HIT_NONE(&hit)) {
- weight = 0.0;
- } else {
+ if(!S3D_HIT_NONE(&hit)) {
struct sschiff_material_properties props;
double n_r, k_r;
double n_e, k_e, lambda_e;
double path_length;
+ double weight;
- res = compute_path_length(scn, &hit, ray_org, axis_z, &path_length);
+ res = compute_path_length(dev, scn, &hit, ray_org, axis_z, &path_length);
if(res != RES_OK) goto error;
mtl->get_property(mtl->material, wavelength, &props);
@@ -232,10 +234,10 @@ compute_monte_carlo_weight
weight = 2.0 * rcp_pdf
* (1.0 - exp(-k_e*k_r*path_length) * cos(k_e*(n_r - 1.0)*path_length));
- }
- accum->weight += weight;
- accum->sqr_weight += weight * weight;
+ accum->weight += weight;
+ accum->sqr_weight += weight * weight;
+ }
accum->naccums += 1;
exit:
@@ -323,7 +325,6 @@ radiative_properties
(struct sschiff_device* dev,
struct ssp_rng* rng,
struct sschiff_integrator_desc* desc,
- const double wavelength,
const int istep,
struct s3d_scene* scene,
struct s3d_shape* shape,
@@ -333,7 +334,8 @@ radiative_properties
float lower[3], upper[3];
float aabb_pt[8][3];
float centroid[3];
- int idir;
+ size_t idir;
+ int mask;
res_T res = RES_OK;
ASSERT(dev && rng && desc && accum);
@@ -364,7 +366,7 @@ radiative_properties
f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f);
- FOR_EACH(idir, 0, NDIRS) {
+ FOR_EACH(idir, 0, desc->sampled_directions_count) {
float dir[4];
float basis[9];
float transform[12];
@@ -392,14 +394,13 @@ radiative_properties
#if 1
(void)istep;
- res = compute_monte_carlo_weight
- (scene, rng, &material, wavelength, basis, centroid, lower, upper, accum);
- if(res != RES_OK) goto error;
+ res = compute_monte_carlo_weight(dev, scene, rng, &material,
+ desc->wavelength, basis, centroid, lower, upper, accum);
+ /*if(res != RES_OK) goto error;*/
#else
/* Compute an image of the shape transformed in the footprint space */
{
char thumb_name[64];
- (void)wavelength;
sprintf(thumb_name, "%d_%d.ppm", istep, idir);
draw_thumbnail(scene, basis, centroid, lower, upper, thumb_name);
}
@@ -408,6 +409,8 @@ radiative_properties
S3D(scene_end_session(scene));
exit:
+ S3D(scene_get_session_mask(scene, &mask));
+ if(mask) S3D(scene_end_session(scene));
return res;
error:
goto exit;
@@ -418,7 +421,10 @@ check_desc(struct sschiff_integrator_desc* desc)
{
ASSERT(desc);
return desc->sample_geometry
- && desc->scattering_angles_count;
+ && desc->scattering_angles_count
+ && desc->sampled_geometries_count
+ && desc->sampled_directions_count
+ && desc->wavelength > 0.f;
}
static void
@@ -471,26 +477,22 @@ sschiff_integrate
(struct sschiff_device* dev,
struct ssp_rng* rng,
struct sschiff_integrator_desc* desc,
- const double wavelength, /* Wavelength to integrate */
- const size_t steps_count, /* #realisations */
struct sschiff_estimator** out_estimator)
{
struct sschiff_estimator* estimator = NULL;
struct s3d_shape* shape = NULL;
struct s3d_scene* scene = NULL;
- int i;
+ size_t i;
res_T res = RES_OK;
- if(!dev || !rng || !desc || !check_desc(desc) || wavelength < 0.0
- || !steps_count || !out_estimator) {
+ if(!dev || !rng || !desc || !check_desc(desc) || !out_estimator) {
res = RES_BAD_ARG;
goto error;
}
res = estimator_create(dev, &estimator);
if(res != RES_OK) goto error;
- estimator->nsteps = steps_count;
- estimator->wavelength = wavelength;
+ estimator->wavelength = desc->wavelength;
res = s3d_shape_create_mesh(dev->s3d, &shape);
if(res != RES_OK) {
@@ -510,8 +512,8 @@ sschiff_integrate
}
memset(&estimator->extinction_cross_section, 0, sizeof(struct accum));
- FOR_EACH(i, 0, (int)steps_count) {
- res = radiative_properties(dev, rng, desc, wavelength, i, scene, shape,
+ FOR_EACH(i, 0, desc->sampled_geometries_count) {
+ res = radiative_properties(dev, rng, desc, (int)i, scene, shape,
&estimator->extinction_cross_section);
if(res != RES_OK) goto error;
}
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -36,8 +36,6 @@
#include <star/s3d.h>
#include <star/ssp.h>
-#define NSTEPS 1000
-
struct result {
double mean_radius;
double extinction_cross_section;
@@ -252,6 +250,9 @@ check_schiff_estimation
desc.sample_geometry = sample_sphere;
desc.scattering_angles_count = 1;
desc.sampler_context = sampler_ctx;
+ desc.sampled_geometries_count = 500;
+ desc.sampled_directions_count = 5;
+ desc.wavelength = 6.e-7;
printf("Schiff sphere estimation - m: %g; n_r: %g; kappa_r: %g; n_e: %g\n",
sampler_ctx->wavelength,
@@ -263,11 +264,13 @@ check_schiff_estimation
sampler_ctx->mean_radius = results[i].mean_radius;
time_current(&t0);
- CHECK(sschiff_integrate(dev, rng, &desc, 6.e-7, NSTEPS, &estimator), RES_OK);
+ CHECK(sschiff_integrate(dev, rng, &desc, &estimator), RES_OK);
time_current(&t1);
CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
+
+
time_sub(&t0, &t1, &t0);
time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
printf(
@@ -280,7 +283,7 @@ check_schiff_estimation
status.E, status.SE,
buf);
- CHECK(eq_eps(status.E, results[i].extinction_cross_section, 2*status.SE), 1);
+ CHECK(eq_eps(status.E, results[i].extinction_cross_section, 3*status.SE), 1);
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
}
}
@@ -335,25 +338,28 @@ main(int argc, char** argv)
sampler_ctx.mean_radius = 1.0;
desc.sample_geometry = sample_sphere;
- desc.scattering_angles_count = 1;
desc.sampler_context = &sampler_ctx;
-
- CHECK(sschiff_integrate(NULL, NULL, NULL, 6e-7, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, NULL, 6e-7, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, NULL, 6e-7, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, NULL, 6e-7, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, &desc, 6e-7, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, &desc, 6e-7, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, &desc, 6e-7, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, &desc, 6e-7, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, NULL, 6e-7, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, NULL, 6e-7, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, NULL, 6e-7, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, NULL, 6e-7, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, NULL, &desc, 6e-7, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, NULL, &desc, 6e-7, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(NULL, rng, &desc, 6e-7, 1, &estimator), RES_BAD_ARG);
- CHECK(sschiff_integrate(dev, rng, &desc, 6e-7, 1, &estimator), RES_OK);
+ 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);
CHECK(sschiff_estimator_get_extinction_cross_section(NULL, NULL), RES_BAD_ARG);
CHECK(sschiff_estimator_get_extinction_cross_section(estimator, NULL), RES_BAD_ARG);
