commit 47de63c3c9a62ef4b0c5a605154c811393d9cdd1
parent e22a1821b16b63a81a32d4b1c33b30d9da6eee71
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 15 Oct 2015 14:50:50 +0200
Fix the Schiff estimator whit several sampled directions
The Monte Carlo weight of the sampled directions were incorrectly
accumulated. Let N directions sampled for a given geometry, the
resulting MC weight is W = 1/N Sum_i^N(w).
Diffstat:
2 files changed, 42 insertions(+), 36 deletions(-)
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -52,6 +52,8 @@ struct accum {
size_t naccums;
};
+static const struct accum ACCUM_NULL = { 0, 0, 0 };
+
struct sschiff_estimator {
double wavelength;
@@ -180,7 +182,7 @@ compute_monte_carlo_weight
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 */
- struct accum* accum)
+ double* weight)
{
struct s3d_hit hit;
double sample[2];
@@ -192,7 +194,7 @@ compute_monte_carlo_weight
float x[3], y[3];
float rcp_pdf;
res_T res = RES_OK;
- ASSERT(scn && rng && mtl && basis && pos && lower && upper && accum);
+ ASSERT(scn && rng && mtl && basis && pos && lower && upper && weight);
f2_sub(plane_size, upper, lower);
rcp_pdf = plane_size[0] * plane_size[1] * 1.e-12f/*Metter^2 to micron^2*/;
@@ -215,12 +217,12 @@ compute_monte_carlo_weight
S3D(scene_trace_ray(scn, ray_org, axis_z, range, &hit));
+ *weight = 0;
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(dev, scn, &hit, ray_org, axis_z, &path_length);
if(res != RES_OK) goto error;
@@ -234,13 +236,9 @@ compute_monte_carlo_weight
lambda_e = wavelength / n_e;
k_e = 2.0 * PI / lambda_e;
- weight = 2.0 * rcp_pdf
+ *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->naccums += 1;
exit:
return res;
@@ -332,6 +330,7 @@ radiative_properties
struct sschiff_material* material,
struct accum* accum)
{
+ double weight_dirs = 0.0;
float lower[3], upper[3];
float aabb_pt[8][3];
float centroid[3];
@@ -360,6 +359,7 @@ radiative_properties
f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f);
FOR_EACH(idir, 0, desc->sampled_directions_count) {
+ double weight;
float dir[4];
float basis[9];
float transform[12];
@@ -386,7 +386,9 @@ radiative_properties
}
compute_monte_carlo_weight(dev, scene, rng, material,
- desc->wavelength, basis, centroid, lower, upper, accum);
+ desc->wavelength, basis, centroid, lower, upper, &weight);
+
+ weight_dirs += weight;
#if 0
/* Compute an image of the shape transformed in the footprint space */
{
@@ -397,6 +399,10 @@ radiative_properties
}
#endif
}
+ weight_dirs /= (double)desc->sampled_directions_count;
+ accum->weight += weight_dirs;
+ accum->sqr_weight += weight_dirs * weight_dirs;
+ accum->naccums += 1;
return RES_OK;
}
@@ -495,7 +501,7 @@ sschiff_integrate
goto error;
}
- memset(&estimator->extinction_cross_section, 0, sizeof(struct accum));
+ estimator->extinction_cross_section = ACCUM_NULL;
FOR_EACH(i, 0, desc->sampled_geometries_count) {
struct sschiff_material material = SSCHIFF_NULL_MATERIAL;
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -41,13 +41,13 @@ struct result {
double extinction_cross_section;
};
-struct sphere {
+struct geometry {
float* vertices; /* List of float[3] */
unsigned* indices; /* List of unsigned[3] */
};
struct sampler_context {
- struct sphere unit_sphere;
+ struct geometry geometry;
double wavelength;
double medium_refractive_index;
double relative_real_refractive_index;
@@ -55,8 +55,8 @@ struct sampler_context {
double mean_radius;
};
-struct sphere_context {
- struct sphere* unit_sphere;
+struct sphere {
+ struct geometry* geometry;
float radius;
};
@@ -82,27 +82,27 @@ get_material_property
static void
get_indices(const unsigned itri, unsigned ids[3], void* ctx)
{
- struct sphere_context* sphere = ctx;
+ struct sphere* sphere = ctx;
const size_t i = itri * 3;
- CHECK(sa_size(sphere->unit_sphere->indices) % 3, 0);
- CHECK(itri < sa_size(sphere->unit_sphere->indices) / 3, 1);
- ids[0] = sphere->unit_sphere->indices[i + 0];
- ids[1] = sphere->unit_sphere->indices[i + 1];
- ids[2] = sphere->unit_sphere->indices[i + 2];
+ CHECK(sa_size(sphere->geometry->indices) % 3, 0);
+ CHECK(itri < sa_size(sphere->geometry->indices) / 3, 1);
+ ids[0] = sphere->geometry->indices[i + 0];
+ ids[1] = sphere->geometry->indices[i + 1];
+ ids[2] = sphere->geometry->indices[i + 2];
}
static INLINE void
get_position(const unsigned ivert, float vertex[3], void* ctx)
{
- struct sphere_context* sphere = ctx;
+ struct sphere* sphere = ctx;
const size_t i = ivert * 3;
- CHECK(sa_size(sphere->unit_sphere->vertices) % 3, 0);
- CHECK(ivert < sa_size(sphere->unit_sphere->vertices) / 3, 1);
- vertex[0] = sphere->unit_sphere->vertices[i + 0] * sphere->radius;
- vertex[1] = sphere->unit_sphere->vertices[i + 1] * sphere->radius;
- vertex[2] = sphere->unit_sphere->vertices[i + 2] * sphere->radius;
+ CHECK(sa_size(sphere->geometry->vertices) % 3, 0);
+ CHECK(ivert < sa_size(sphere->geometry->vertices) / 3, 1);
+ vertex[0] = sphere->geometry->vertices[i + 0] * sphere->radius;
+ vertex[1] = sphere->geometry->vertices[i + 1] * sphere->radius;
+ vertex[2] = sphere->geometry->vertices[i + 2] * sphere->radius;
}
static res_T
@@ -114,14 +114,14 @@ sample_sphere
{
struct s3d_vertex_data attrib;
struct sampler_context* sampler_ctx = sampler_context;
- struct sphere_context sphere;
+ struct sphere sphere;
size_t nverts, nprims;
NCHECK(rng, NULL);
NCHECK(mtl, NULL);
NCHECK(shape, NULL);
- sphere.unit_sphere = &sampler_ctx->unit_sphere;
+ sphere.geometry = &sampler_ctx->geometry;
sphere.radius = (float)ssp_ran_lognormal
(rng, log(sampler_ctx->mean_radius), log(1.18));
@@ -132,15 +132,15 @@ sample_sphere
mtl->get_property = get_material_property;
mtl->material = sampler_ctx;
- nverts = sa_size(sphere.unit_sphere->vertices) / 3/*#coords*/;
- nprims = sa_size(sphere.unit_sphere->indices) / 3/*#indices per prim*/;
+ nverts = sa_size(sphere.geometry->vertices) / 3/*#coords*/;
+ nprims = sa_size(sphere.geometry->indices) / 3/*#indices per prim*/;
return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims,
get_indices, (unsigned)nverts, &attrib, 1, &sphere);
}
static void
-unit_sphere_init(struct sphere* sphere, const unsigned ntheta)
+sphere_init(struct geometry* sphere, const unsigned ntheta)
{
const unsigned nphi = (unsigned)(((double)ntheta + 0.5) / 2.0);
const double step_theta = 2*PI / (double)ntheta;
@@ -220,7 +220,7 @@ unit_sphere_init(struct sphere* sphere, const unsigned ntheta)
}
static void
-unit_sphere_release(struct sphere* sphere)
+sphere_release(struct geometry* sphere)
{
sa_release(sphere->vertices);
sa_release(sphere->indices);
@@ -250,8 +250,8 @@ check_schiff_estimation
desc.sample_geometry = sample_sphere;
desc.scattering_angles_count = 1;
desc.sampler_context = sampler_ctx;
- desc.sampled_geometries_count = 1000;
- desc.sampled_directions_count = 1;
+ desc.sampled_geometries_count = 100;
+ desc.sampled_directions_count = 100;
desc.wavelength = 6.e-7;
printf("Schiff sphere estimation - m: %g; n_r: %g; kappa_r: %g; n_e: %g\n",
@@ -327,7 +327,7 @@ main(int argc, char** argv)
CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
CHECK(sschiff_device_create(NULL, &allocator, 0, NULL, &dev), RES_OK);
- unit_sphere_init(&sampler_ctx.unit_sphere, 64);
+ sphere_init(&sampler_ctx.geometry, 64);
sampler_ctx.wavelength = 0.6e-6;
sampler_ctx.relative_real_refractive_index = 1.1;
@@ -381,7 +381,7 @@ main(int argc, char** argv)
CHECK(sschiff_device_ref_put(dev), RES_OK);
CHECK(ssp_rng_ref_put(rng), RES_OK);
- unit_sphere_release(&sampler_ctx.unit_sphere);
+ sphere_release(&sampler_ctx.geometry);
check_memory_allocator(&allocator);
mem_shutdown_proxy_allocator(&allocator);
CHECK(mem_allocated_size(), 0);
