commit 50f608d49311598fe30ebd099457e50eb1a7d791
parent 3ffacc05990b35757f550cbfd556de89af1ecd6d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 1 Apr 2016 16:02:37 +0200
Add the volume scaling data to the geometry distribution
The sampled geometry provide a scale factor of its volume that is taken
into account by the integrator to adjust the length of the sampled
radiative path and the projective area of the particle. This feature is
well suited for a distribution of geometries that have the same shape
but a different volume.
Diffstat:
6 files changed, 81 insertions(+), 28 deletions(-)
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -72,6 +72,7 @@ struct sschiff_geometry_distribution {
res_T (*sample) /* Sample a geometry i.e. a shape and its material */
(struct ssp_rng* rng,
struct s3d_shape* shape, /* Sampled shape */
+ double* scale_factor, /* Scale factor to apply to the shape volume */
void* context);
void* context;
};
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -249,6 +249,7 @@ compute_path_length
const struct s3d_hit* first_hit,
const float ray_org[3],
const float ray_dir[3],
+ const double scale,
double* length)
{
float range[2] = { 0, FLT_MAX };
@@ -256,6 +257,7 @@ compute_path_length
double len = 0;
float dst;
ASSERT(scn && first_hit && ray_org && ray_dir && !S3D_HIT_NONE(first_hit));
+ ASSERT(scale > 0);
(void)dev;
dst = range[0] = first_hit->distance;
@@ -284,7 +286,7 @@ compute_path_length
} while(!S3D_HIT_NONE(&hit));
- *length = len;
+ *length = len * scale;
return RES_OK;
}
@@ -398,7 +400,9 @@ accum_monte_carlo_weight
const float basis[9], /* World space basis of the RT volume */
const float pos[3], /* World space centroid of the RT volume */
const float lower[3], /* Lower boundary of the RT volume */
- const float upper[3]) /* Upper boundary of the RT volume */
+ const float upper[3], /* Upper boundary of the RT volume */
+ const double area_scaling, /* Scale factor to apply to the areas */
+ const double dist_scaling) /* Scale factor to apply to the distances */
{
struct s3d_hit hit;
double sample[2];
@@ -412,9 +416,10 @@ accum_monte_carlo_weight
float plane_area;
res_T res = RES_OK;
ASSERT(dev && ctx && basis && pos && lower && upper);
+ ASSERT(area_scaling > 0 && dist_scaling > 0);
f2_sub(plane_size, upper, lower); /* In micron */
- plane_area = plane_size[0] * plane_size[1];
+ plane_area = plane_size[0] * plane_size[1] * (float)area_scaling;
/* Define the projection axis */
f3_mulf(axis_x, basis + 0, plane_size[0] * 0.5f);
@@ -440,8 +445,8 @@ accum_monte_carlo_weight
if(S3D_HIT_NONE(&hit))
break;
- res = compute_path_length
- (dev, ctx->scene, &hit, ray_org[0], axis_z, &path_length[0]);
+ res = compute_path_length(dev, ctx->scene, &hit, ray_org[0], axis_z,
+ dist_scaling, &path_length[0]);
if(res != RES_OK) { /* Handle numerical/geometry issues */
++nfailures;
continue;
@@ -464,8 +469,8 @@ accum_monte_carlo_weight
if(S3D_HIT_NONE(&hit)) /* NULL differential cross section weight */
break;
- res = compute_path_length
- (dev, ctx->scene, &hit, ray_org[1], axis_z, &path_length[1]);
+ res = compute_path_length(dev, ctx->scene, &hit, ray_org[1], axis_z,
+ dist_scaling, &path_length[1]);
if(res != RES_OK) { /* Handle numerical/geometry issues */
++nfailures;
continue;
@@ -491,6 +496,8 @@ radiative_properties
(struct sschiff_device* dev,
const int istep,
const size_t ndirs,
+ const double area_scaling, /* Scale factor to apply to the areas */
+ const double dist_scaling, /* Scale factor to apply to the distances */
struct integrator_context* ctx)
{
float lower[3], upper[3];
@@ -500,7 +507,7 @@ radiative_properties
size_t iwlen;
int i;
res_T res = RES_OK;
- ASSERT(dev && ndirs && ctx);
+ ASSERT(dev && ndirs && ctx && area_scaling > 0 && dist_scaling > 0);
(void)istep;
S3D(scene_get_aabb(ctx->scene, lower, upper));
@@ -565,7 +572,8 @@ radiative_properties
f3_max(upper, upper, pt[i]);
}
- res = accum_monte_carlo_weight(dev, ctx, basis, centroid, lower, upper);
+ res = accum_monte_carlo_weight
+ (dev, ctx, basis, centroid, lower, upper, area_scaling, dist_scaling);
if(res != RES_OK) goto error;
}
/* Compute the Monte Carlo weight of the temporary accumulator */
@@ -1156,17 +1164,24 @@ static res_T
begin_realisation
(struct sschiff_device* dev,
struct sschiff_geometry_distribution* distrib,
+ double* volume_scaling,
struct integrator_context* ctx)
{
res_T res = RES_OK;
- ASSERT(dev && distrib && ctx);
+ ASSERT(dev && distrib && volume_scaling && ctx);
/* Sample a particle */
- res = distrib->sample(ctx->rng, ctx->shape, distrib->context);
+ res = distrib->sample
+ (ctx->rng, ctx->shape, volume_scaling, distrib->context);
if(res != RES_OK) {
log_error(dev, "Error in sampling a particle.\n");
goto error;
}
+ if(*volume_scaling <= 0) {
+ log_error(dev, "Invalid volume scale factor `%g'.\n", *volume_scaling);
+ res = RES_BAD_ARG;
+ goto error;
+ }
/* Build the Star-3D representation of the sampled shape */
S3D(scene_begin_session(ctx->scene, S3D_TRACE));
@@ -1505,19 +1520,28 @@ sschiff_integrate
#pragma omp parallel for schedule(static)
for(igeom=0; igeom < (int)ngeoms; ++igeom) {
const int ithread = omp_get_thread_num();
+ double volume_scaling;
+ double area_scaling;
+ double dist_scaling;
ATOMIC res_local = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue;
/* Setup the data for the current realisation */
- res_local = begin_realisation(dev, distrib, ctxs + ithread);
+ res_local = begin_realisation(dev, distrib, &volume_scaling, ctxs+ithread);
if(res_local != RES_OK) {
ATOMIC_CAS(&res, res_local, RES_OK);
continue;
}
+ /* Define the scale factor to apply the the area and the distance from the
+ * scale factor of the volume */
+ area_scaling = pow(volume_scaling, 2.0/3.0);
+ dist_scaling = pow(volume_scaling, 1.0/3.0);
+
/* Schiff Estimation */
- res_local = radiative_properties(dev, igeom, ndirs, ctxs+ithread);
+ res_local = radiative_properties
+ (dev, igeom, ndirs, area_scaling, dist_scaling, ctxs+ithread);
if(res_local != RES_OK) ATOMIC_CAS(&res, res_local, RES_OK);
end_realisation(ctxs + ithread);
diff --git a/src/test_sschiff_estimator_cylinder.c b/src/test_sschiff_estimator_cylinder.c
@@ -35,8 +35,8 @@ struct result {
};
struct sampler_context {
- struct geometry geometry;
- double aspect_ratio;
+ struct s3d_shape* shape;
+ double cylinder_volume;
double mean_radius;
double sigma;
};
@@ -56,19 +56,20 @@ get_material_property
static res_T
sample_cylinder
- (struct ssp_rng* rng, struct s3d_shape* shape, void* sampler_context)
+ (struct ssp_rng* rng,
+ struct s3d_shape* shape,
+ double* volume_scaling,
+ void* sampler_context)
{
struct sampler_context* sampler_ctx = sampler_context;
- struct cylinder cylinder;
+ double sphere_volume;
double sample;
(void)rng, (void)sampler_context;
sample = ssp_ran_lognormal(rng, log(sampler_ctx->mean_radius), log(sampler_ctx->sigma));
- cylinder.geometry = &sampler_ctx->geometry;
- cylinder.radius = (float)(sample / pow(3.0 / (2.0*sampler_ctx->aspect_ratio), 1.0/3.0));
- cylinder.height = (float)(2.f * cylinder.radius / sampler_ctx->aspect_ratio);
-
- return cylinder_setup_s3d_shape(&cylinder, shape);
+ sphere_volume = 4.0*PI*sample*sample*sample / 3.0;
+ *volume_scaling = sphere_volume / sampler_ctx->cylinder_volume;
+ return s3d_mesh_copy(sampler_ctx->shape, shape);
}
static INLINE void
@@ -403,20 +404,33 @@ main(int argc, char** argv)
1.977109116e+2, 1.985486697e+2, 1.993864277e+2, 2.002241858e+2,
2.010619438e+2, 2.018997018e+2, 2.027374599e+2
};
+ struct s3d_device* s3d = NULL;
+ struct s3d_shape* shape = NULL;
+ struct geometry geometry;
+ struct cylinder cylinder;
const size_t nx = sizeof(x)/sizeof(double);
const size_t nscatt_angles = 1000;
const size_t ngeoms = 100;
const size_t ndirs = 100;
+ double cylinder_volume = 0;
size_t i;
(void)argc, (void)argv;
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+ CHECK(s3d_device_create(NULL, NULL, 0, &s3d), RES_OK);
+ CHECK(s3d_shape_create_mesh(s3d, &shape), RES_OK);
+
CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
CHECK(sschiff_device_create
- (NULL, &allocator, SSCHIFF_NTHREADS_DEFAULT, 1, NULL, &dev), RES_OK);
+ (NULL, &allocator, SSCHIFF_NTHREADS_DEFAULT, 1, s3d, &dev), RES_OK);
- geometry_init_cylinder(&sampler_ctx.geometry, 64);
+ geometry_init_cylinder(&geometry, 64);
+ cylinder.geometry = &geometry;
+ cylinder.radius = 1.0;
+ cylinder.height = 1.0;
+ cylinder_setup_s3d_shape(&cylinder, shape);
+ cylinder_volume = PI * cylinder.radius * cylinder.radius * cylinder.height;
FOR_EACH(i, 0, nx) {
const double wavelength = 0.6; /* In micron */
@@ -425,9 +439,11 @@ main(int argc, char** argv)
struct sschiff_state* val;
struct sschiff_state result;
- sampler_ctx.aspect_ratio = 0.837;
+ /*sampler_ctx.aspect_ratio = 0.837;*/
sampler_ctx.mean_radius = (x[i] * 0.450) / (2*PI);
sampler_ctx.sigma = 1.18;
+ sampler_ctx.shape = shape;
+ sampler_ctx.cylinder_volume = cylinder_volume;
distrib.characteristic_length = sampler_ctx.mean_radius;
distrib.material.get_property = get_material_property;
@@ -491,7 +507,10 @@ main(int argc, char** argv)
CHECK(sschiff_device_ref_put(dev), RES_OK);
CHECK(ssp_rng_ref_put(rng), RES_OK);
- geometry_release(&sampler_ctx.geometry);
+ CHECK(s3d_device_ref_put(s3d), RES_OK);
+ CHECK(s3d_shape_ref_put(shape), RES_OK);
+
+ geometry_release(&geometry);
check_memory_allocator(&allocator);
mem_shutdown_proxy_allocator(&allocator);
CHECK(mem_allocated_size(), 0);
diff --git a/src/test_sschiff_estimator_rhodo.c b/src/test_sschiff_estimator_rhodo.c
@@ -2551,17 +2551,23 @@ get_material_property
static res_T
sample_cylinder
- (struct ssp_rng* rng, struct s3d_shape* shape, void* sampler_context)
+ (struct ssp_rng* rng,
+ struct s3d_shape* shape,
+ double* volume_scaling,
+ void* sampler_context)
{
struct sampler_context* sampler_ctx = sampler_context;
struct cylinder cylinder;
double sample;
(void)rng, (void)sampler_context;
+ NCHECK(volume_scaling, NULL);
+
sample = ssp_ran_lognormal(rng, log(sampler_ctx->mean_radius), log(sampler_ctx->sigma));
cylinder.geometry = &sampler_ctx->geometry;
cylinder.radius = (float)(sample / pow(3.0 / (2.0*sampler_ctx->aspect_ratio), 1.0/3.0));
cylinder.height = (float)(2.f * cylinder.radius / sampler_ctx->aspect_ratio);
+ *volume_scaling = 1.0;
return cylinder_setup_s3d_shape(&cylinder, shape);
}
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -63,16 +63,19 @@ static res_T
sample_sphere
(struct ssp_rng* rng,
struct s3d_shape* shape,
+ double* volume_scaling,
void* sampler_context)
{
struct sampler_context* sampler_ctx = sampler_context;
struct sphere sphere;
NCHECK(rng, NULL);
+ NCHECK(volume_scaling, NULL);
sphere.geometry = &sampler_ctx->geometry;
sphere.radius = (float)ssp_ran_lognormal
(rng, log(sampler_ctx->mean_radius), log(sampler_ctx->sigma));
+ *volume_scaling = 1.0;
return sphere_setup_s3d_shape(&sphere, shape);
}
diff --git a/src/test_sschiff_utils.h b/src/test_sschiff_utils.h
@@ -318,7 +318,7 @@ sphere_setup_s3d_shape(struct sphere* sphere, struct s3d_shape* shape)
static INLINE void
compute_estimation_intersection
(double intersection[2],
- const double scale,
+ const double scale,
const struct sschiff_state* state0,
const struct sschiff_state* state1)
{
