commit 5e365c2c362833832829289c92a23fa25b41deb9
parent f519ed6cab752075a1749e6656a2a75e59743deb
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 17 Feb 2016 10:02:01 +0100
Big refactor of the estimator code
Add the Monte-Carlo estimation of the diffusion cross section and its
cumulative. Make optical properties constants for a geometry
distribution.
Handle Ray-Tracing self intersection on primitive edges and fix an issue
in the rejection of rays that stopped the integration process.
Diffstat:
5 files changed, 790 insertions(+), 354 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -41,11 +41,13 @@ find_package(RSys 0.3 REQUIRED)
find_package(StarSP 0.3 REQUIRED)
find_package(Star3D 0.3 REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
+find_package(GSL REQUIRED)
include_directories(
${RSys_INCLUDE_DIR}
${StarSP_INCLUDE_DIR}
- ${Star3D_INCLUDE_DIR})
+ ${Star3D_INCLUDE_DIR}
+ ${GSL_INCLUDE_DIRS})
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
@@ -79,7 +81,7 @@ set_target_properties(sschiff PROPERTIES
VERSION ${VERSION}
SOVERSION ${VERSION_MAJOR})
-target_link_libraries(sschiff RSys Star3D StarSP)
+target_link_libraries(sschiff RSys Star3D StarSP GSL::gsl)
if(CMAKE_COMPILER_IS_GNUCC)
target_link_libraries(sschiff m)
endif()
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -58,15 +58,6 @@ struct s3d_device;
struct s3d_shape;
struct ssp_rng;
-/* Type of estimated data */
-enum sschiff_data {
- SSCHIFF_EXTINCTION_CROSS_SECTION,
- SSCHIFF_ABSORPTION_CROSS_SECTION,
- SSCHIFF_SCATTERING_CROSS_SECTION,
- SSCHIFF_AVERAGE_PROJECTED_AREA,
- SSCHIFF_DATA_COUNT__
-};
-
/* Optical properties of a material handled by the Schiff integrator */
struct sschiff_material_properties {
double medium_refractive_index;
@@ -89,25 +80,30 @@ static const struct sschiff_material SSCHIFF_NULL_MATERIAL = { NULL, NULL };
/* User defined distribution of the geometries. The unit of the geometry is the
* micron, i.e. 1.0f == 1 micron */
struct sschiff_geometry_distribution {
+ struct sschiff_material material; /* Material of the geometry distribution */
+ double caracteristic_length;
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* context);
void* context;
};
static const struct sschiff_geometry_distribution
-SSCHIFF_NULL_GEOMETRY_DISTRIBUTION = { NULL, NULL };
+SSCHIFF_NULL_GEOMETRY_DISTRIBUTION = { {NULL, NULL}, -1.0, NULL, NULL };
/* State of the Monte Carlo estimation */
-struct sschiff_estimator_state {
- double wavelength; /* Estimated wavelength in micron */
- struct sschiff_estimator_value { /* Values of the estimation */
- double E; /* Expected value */
- double V; /* Variance */
- double SE; /* Standard error */
- } values[SSCHIFF_DATA_COUNT__];
+struct sschiff_state {
+ double E; /* Expected value */
+ double V; /* Variance */
+ double SE; /* Standard error */
+};
+
+struct sschiff_cross_section {
+ struct sschiff_state extinction;
+ struct sschiff_state absorption;
+ struct sschiff_state scattering;
+ struct sschiff_state average_projected_area;
};
/* Type of the function use to distribute the scattering angles in [0, PI].
@@ -155,7 +151,7 @@ sschiff_integrate
const double* wavelengths, /* list of wavelengths to estimate in micron */
const size_t wavelengths_count, /* # wavelengths to estimate */
const sschiff_scattering_angles_distribution_T angles, /* angles distrib */
- const size_t scattering_angles_count, /* # scattering angles. Must be >= 2 */
+ const size_t scattering_angles_count, /* # scattering angles. Must be >= 3 */
const size_t sampled_geometries_count, /* # geometries to sample */
const size_t sampled_directions_count, /* # directions to sample per geometry */
struct sschiff_estimator** estimator);
@@ -182,19 +178,19 @@ sschiff_estimator_get_realisations_count
/* Retrieve the estimation state of a given wavelength */
SSCHIFF_API res_T
-sschiff_estimator_get_wavelength_state
+sschiff_estimator_get_wavelength_cross_section
(const struct sschiff_estimator* estimator,
const double wavelength, /* In micron */
- struct sschiff_estimator_state* state);
+ struct sschiff_cross_section* cross_section);
-/* Retrieve the estimation state of all wavelengths. The length of `states'
- * must be at least equal to the count of integrated wavelengths. One can use
- * the sschiff_estimator_get_wavelengths_count function to retrieve this
- * information. */
+/* Retrieve the estimated cross sections of all wavelengths. The length of
+ * `cross_sections' must be at least equal to the count of integrated
+ * wavelengths. One can use the sschiff_estimator_get_wavelengths_count
+ * function to retrieve this information. */
SSCHIFF_API res_T
-sschiff_estimator_get_states
+sschiff_estimator_get_cross_sections
(const struct sschiff_estimator* estimator,
- struct sschiff_estimator_state states[]);
+ struct sschiff_cross_section cross_sections[]);
END_DECLS
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -31,6 +31,7 @@
#include "sschiff.h"
#include "sschiff_device.h"
+#include <rsys/algorithm.h>
#include <rsys/float2.h>
#include <rsys/float3.h>
#include <rsys/image.h>
@@ -45,25 +46,75 @@
#include <math.h>
#include <omp.h>
+#include <gsl/gsl_sf.h>
+
#define MAX_FAILURES 10
/* Accumulator of double precision values */
-struct accum { double weights[SSCHIFF_DATA_COUNT__]; };
+struct accum {
+ double extinction_cross_section;
+ double absorption_cross_section;
+ double scattering_cross_section;
+ double avg_projected_area;
+
+ /* Per scattering angles weights */
+ double* differential_cross_section;
+ double* differential_cross_section_cumulative;
+};
/* Monte carlo data */
struct mc_data { double weight; double sqr_weight; };
-
static const struct mc_data MC_DATA_NULL = { 0.0, 0.0 };
/* Accumulator of Monte Carlo data */
-struct mc_accum { struct mc_data mc_data[SSCHIFF_DATA_COUNT__]; };
+struct mc_accum {
+ struct mc_data extinction_cross_section;
+ struct mc_data absorption_cross_section;
+ struct mc_data scattering_cross_section;
+ struct mc_data avg_projected_area;
+
+ /* Per scattering angles weight */
+ struct mc_data* differential_cross_section;
+ struct mc_data* differential_cross_section_cumulative;
+};
+
+/* The geometry context stores data that are constants for a given geometry
+ * distribution. It stores limit scattering angles as well as geometry material
+ * since the particle of a distribution shared the same optical properties.
+ * Precomputed values relying onto these data are also stored in order to
+ * improve the integration efficiency */
+struct geometry_context {
+ const double* angles; /* Pointer toward the estimated scattering angles */
+
+ /* Per wavelength index of the scattering angle from which the MC estimation
+ * of the phase function is no more valid. Must be < nangles. */
+ size_t* limit_angles;
+
+ double* medium_refractive_indices;
+ struct sschiff_material_properties* properties; /* Material properties */
+};
+
+/* The integrator context stores per thread data */
+struct integrator_context {
+ /* Geometry of the micro particle */
+ struct s3d_shape* shape;
+ struct s3d_scene* scene;
+
+ struct ssp_rng* rng;
+ const struct geometry_context* geom_ctx; /* Pointer onto the constant data */
-static const struct mc_accum MC_ACCUM_NULL =
-{{{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}};
+ size_t nwavelengths; /* #wavelengths */
+ size_t nangles; /* #scattering angles */
+
+ /* Per wavelengths data */
+ struct accum* accums; /* Temporary accumulators */
+ struct mc_accum* mc_accums; /* Monte Carlo accumulators */
+};
struct sschiff_estimator {
- double* wavelengths;
- struct mc_accum* accums;
+ double* wavelengths; /* List of wavelengths in vacuum */
+ double* angles; /* List of scattering angles */
+ struct mc_accum* accums; /* Monte Carlo estimation */
struct sschiff_device* dev;
size_t nrealisations;
ref_T ref;
@@ -214,6 +265,38 @@ draw_thumbnail
#undef DEFINITION
}
+static FINLINE int
+hit_on_edge(const struct s3d_hit* hit)
+{
+ const float on_edge_eps = 1.e-4f;
+ float w;
+ ASSERT(hit && !S3D_HIT_NONE(hit));
+ w = 1.f - hit->uv[0] - hit->uv[1];
+ return eq_epsf(hit->uv[0], 0.f, on_edge_eps)
+ || eq_epsf(hit->uv[0], 1.f, on_edge_eps)
+ || eq_epsf(hit->uv[1], 0.f, on_edge_eps)
+ || eq_epsf(hit->uv[1], 1.f, on_edge_eps)
+ || eq_epsf(w, 0.f, on_edge_eps)
+ || eq_epsf(w, 1.f, on_edge_eps);
+}
+
+static FINLINE int
+self_hit
+ (const struct s3d_hit* hit_from,
+ const struct s3d_hit* hit_to,
+ const float epsilon)
+{
+ ASSERT(hit_from && hit_to);
+
+ if(S3D_HIT_NONE(hit_from) || S3D_HIT_NONE(hit_to))
+ return 0;
+ if(S3D_PRIMITIVE_EQ(&hit_from->prim, &hit_to->prim))
+ return 1;
+ if(eq_epsf(hit_to->distance - hit_from->distance, 0.f, epsilon))
+ return hit_on_edge(hit_from) && hit_on_edge(hit_to);
+ return 0;
+}
+
/* Compute the length in micron of the ray part that traverses the scene */
static res_T
compute_path_length
@@ -225,38 +308,38 @@ compute_path_length
double* length)
{
float range[2] = { 0, FLT_MAX };
- struct s3d_primitive prim_from;
- struct s3d_hit hit;
+ struct s3d_hit hit_from, hit;
double len = 0;
float dst;
ASSERT(scn && first_hit && ray_org && ray_dir && !S3D_HIT_NONE(first_hit));
dst = range[0] = first_hit->distance;
- prim_from = first_hit->prim;
+ hit_from = *first_hit;
+
do {
do {
range[0] = nextafterf(range[0], range[1]);
S3D(scene_trace_ray(scn, ray_org, ray_dir, range, &hit));
- } while(S3D_PRIMITIVE_EQ(&hit.prim, &prim_from) /* Self-intersection */
- || hit.distance < range[0]); /* Precision issue */
+ } while(hit.distance < range[0] || self_hit(&hit_from, &hit, 1.e-6f));
if(S3D_HIT_NONE(&hit)) {
- log_error(dev, "Error in computing the radiative path length.\n");
- return RES_BAD_ARG;
+ log_error(dev,
+"Error in computing the radiative path length. This may be due to numerical\n"
+"imprecisions or to a geometry that does not define a closed volume.\n");
+ return RES_BAD_OP;
}
len += hit.distance - dst;
dst = range[0] = hit.distance;
- prim_from = hit.prim;
+ hit_from = hit;
do {
range[0] = nextafterf(range[0], range[1]);
S3D(scene_trace_ray(scn, ray_org, ray_dir, range, &hit));
- } while(S3D_PRIMITIVE_EQ(&hit.prim, &prim_from) /* Self-intersection */
- || hit.distance < range[0]); /* Precision issue */
+ } while(hit.distance < range[0] || self_hit(&hit_from, &hit, 1.e-6f));
range[0] = hit.distance;
- prim_from = hit.prim;
+ hit_from = hit;
} while(!S3D_HIT_NONE(&hit));
@@ -264,63 +347,211 @@ compute_path_length
return RES_OK;
}
-static FINLINE void
+static INLINE void
accum_cross_section
- (const struct sschiff_material_properties* props,
- const double* wavelengths,
- const size_t nwavelengths,
- const double rcp_pdf,
- const double path_length,
- struct accum* accums)
+ (struct integrator_context* ctx,
+ const double plane_area,
+ const double path_length)
{
size_t iwlen;
- ASSERT(props && wavelengths && nwavelengths && rcp_pdf>0 && path_length>0);
+ ASSERT(ctx && plane_area > 0 && path_length > 0);
- FOR_EACH(iwlen, 0, nwavelengths) {
+ FOR_EACH(iwlen, 0, ctx->nwavelengths) {
double n_r, k_r;
- double n_e, k_e, lambda_e;
+ double k_e;
double w_extinction, w_absorption, w_scattering;
- n_e = props[iwlen].medium_refractive_index;
- n_r = props[iwlen].relative_real_refractive_index;
- k_r = props[iwlen].relative_imaginary_refractive_index;
-
- lambda_e = wavelengths[iwlen] / n_e;
- k_e = 2.0 * PI / lambda_e;
+ k_e = ctx->geom_ctx->medium_refractive_indices[iwlen];
+ n_r = ctx->geom_ctx->properties[iwlen].relative_real_refractive_index;
+ k_r = ctx->geom_ctx->properties[iwlen].relative_imaginary_refractive_index;
/* Extinction cross section */
- w_extinction = 2.0 * rcp_pdf *
+ w_extinction = 2.0 * plane_area *
(1.0 - exp(-k_e*k_r*path_length) * cos(k_e*(n_r - 1.0)*path_length));
- accums[iwlen].weights[SSCHIFF_EXTINCTION_CROSS_SECTION] += w_extinction;
+ ctx->accums[iwlen].extinction_cross_section += w_extinction;
/* Absorption cross section */
- w_absorption = rcp_pdf * (1.0 - exp(-2.0*k_e*k_r*path_length));
- accums[iwlen].weights[SSCHIFF_ABSORPTION_CROSS_SECTION] += w_absorption;
+ w_absorption = plane_area * (1.0 - exp(-2.0*k_e*k_r*path_length));
+ ctx->accums[iwlen].absorption_cross_section += w_absorption;
/* Scattering cross section */
w_scattering = w_extinction - w_absorption;
- accums[iwlen].weights[SSCHIFF_SCATTERING_CROSS_SECTION] += w_scattering;
+ ctx->accums[iwlen].scattering_cross_section += w_scattering;
/* Average projected area */
- accums[iwlen].weights[SSCHIFF_AVERAGE_PROJECTED_AREA] += rcp_pdf;
+ ctx->accums[iwlen].avg_projected_area += plane_area;
}
}
+static INLINE void
+accum_differential_cross_section
+ (struct integrator_context* ctx,
+ const double plane_area,
+ const double path_length[2],
+ float ray_org[2][3])
+{
+ double delta_r;
+ size_t iwlen;
+ float vec[3];
+
+ /* Compute the length between the 2 ray starting points in footprint space */
+ delta_r = f3_len(f3_sub(vec, ray_org[0], ray_org[1]));
+
+ FOR_EACH(iwlen, 0, ctx->nwavelengths) {
+ double k_e_delta_r;
+ double n_r, k_r;
+ double k_e;
+ double beta_r[2];
+ double beta_i[2];
+ double tmp;
+ size_t iangle;
+
+ /* Fetch optical properties */
+ k_e = ctx->geom_ctx->medium_refractive_indices[iwlen];
+ n_r = ctx->geom_ctx->properties[iwlen].relative_real_refractive_index;
+ k_r = ctx->geom_ctx->properties[iwlen].relative_imaginary_refractive_index;
+
+ /* Precompute some values. TODO can be stored on the geometry context at
+ * its initialisation. This should improve performances*/
+ k_e_delta_r = k_e * delta_r;
+ beta_r[0] = k_e * (n_r - 1) * path_length[0];
+ beta_r[1] = k_e * (n_r - 1) * path_length[1];
+ beta_i[0] = k_e * k_r * path_length[0];
+ beta_i[1] = k_e * k_r * path_length[1];
+ tmp = (k_e * plane_area) / (2*PI);
+ tmp *= tmp;
+ tmp *= 1
+ + exp(-(beta_i[0] + beta_i[1])) * cos(beta_r[0] - beta_r[1])
+ - exp(-beta_i[0]) * cos(beta_r[0])
+ - exp(-beta_i[1]) * cos(beta_r[1]);
+
+ /* Compute and accumulate the MC weights of the differential cross section
+ * and its cumulative. Note that ctx->limit_angles store the index of the
+ * first scattering angle whos phase function is not estimated by
+ * Monte-Carlo */
+ FOR_EACH(iangle, 0, ctx->geom_ctx->limit_angles[iwlen]) {
+ double bessel_j0, bessel_j1;
+ double weight;
+ double k_e_angle_delta_r;
+
+ k_e_angle_delta_r = ctx->geom_ctx->angles[iangle] * k_e_delta_r;
+
+ bessel_j0 = gsl_sf_bessel_j0(k_e_angle_delta_r);
+ weight = bessel_j0 * tmp;
+ ctx->accums[iwlen].differential_cross_section[iangle] += weight;
+
+ bessel_j1 = gsl_sf_bessel_j1(k_e_angle_delta_r);
+ weight = 2*PI * ctx->geom_ctx->angles[iangle] * bessel_j1 / k_e_delta_r * tmp;
+ ctx->accums[iwlen].differential_cross_section_cumulative[iangle] += weight;
+ }
+ }
+}
+
+#if 0
+struct genevieve_arg {
+ double limit_angle;
+ double scattering_cross_section;
+ double limit_differential_cross_section;
+ double limit_differential_cross_section_cumulative;
+};
+
+int
+genevieve(const gsl_vector* x, void* params, gsl_vector* f)
+{
+ struct genevieve_arg* arg = params;
+ const double x = gsl_vector_get(x, 0);
+ double f0;
+ double square_cos_limit_angle;
+ double square_half_sin_limit_angle;
+
+ if(x==2 || x==4 || x==6)
+ return GSL_EDOM;
+
+ square_cos_limit_angle = cos(limit_angle);
+ square_cos_limit_angle *= square_cos_limit_angle;
+
+ square_half_sin_limit_angle = sin(limit_angle) / 2;
+ square_half_sin_limit_angle *= square_half_sin_limit_angle;
+
+ f0 = arg->limit_differential_cross_section_cumulative
+ - arg->scattering_cross_section
+ + (4*PI * arg->limit_differential_cross_section) / (1+square_cos_limit_angle)
+ * (square_half_sin_limit_angle * ((x*x - 6*x + 8) * cos(2*arg->limit_angle) + 3*x*x - 8 * (x - 2) * cos(arg->limit_angle) - 26*x + 72) - pow(sin(arg->limit_angle/2.0), x) * (4*x*x - 24*x + 40))
+ / ((x-2)*(x-4)*(x-6));
+
+ return GSL_SUCCESS;
+}
+
+static void
+phase_function_post_process()
+{
+ const gsl_multiroot_fsolver* solver = NULL;
+ const gsl_multiroot_function func;
+ double n;
+ double r;
+ struct genevieve_param* arg;
+ struct gsl_vector* x = NULL;
+ struct gsl_vector* root;
+ res_T res = RES_OK;
+
+ solver = gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_dnewton, 1);
+ if(!solver) {
+ log_error(dev,
+ "Not enough memory. Couldn't allocate the GSL Newton-Raphson solver.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ x = gsl_vector_alloc(1);
+ if(!x) {
+ log_error(dev,
+ "Not enough memory. Couldn't allocat the initial GSL vector.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* TODO fill arg */
+ func.f = genevieve;
+ func.n = 1;
+ func.params = &arg;
+
+ gsl_vector_set(x, 0, 3.0);
+ gsl_multiroot_fsolver_set(solver, f, x);
+
+ FOR_EACH(iwlen, 0, nwavelengths) {
+ int status;
+ int iter;
+ do {
+ status = gsl_multiroot_fsolver_iterate(solver);
+ if(states != GSL_SUCCESS) break; /* ERROR */
+
+ status = gsl_multiroot_test_residual(solver->f, 1.e-3);
+ } while(status == GSL_CONTINUE && iter < 1000);
+ }
+
+ root = gsl_multiroot_fsolver_root(solver->x);
+ n = gsl_vector_get(root, 0);
+
+ r = 2 * limit_differential_cross_section * pow(sin(limit_angle/2), n)
+ / (1 + pow(cos(limit_angle), 2));
+
+exit:
+ if(solver) gsl_multiroot_fsolver_free(solver);
+ if(x) gsl_vector_free(x);
+ return res;
+error:
+ goto exit;
+}
+#endif
+
static res_T
accum_monte_carlo_weight
(struct sschiff_device* dev,
- struct s3d_scene* scn,
- struct ssp_rng* rng,
- const struct sschiff_material_properties* props,
- const double* wavelengths, /* In micron */
- const size_t nwavelengths,
- const double* angles, /* Scattering angles */
- const size_t nangles,
+ struct integrator_context* ctx,
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 */
- struct accum* accums)
+ const float upper[3]) /* Upper boundary of the RT volume */
{
struct s3d_hit hit;
double sample[2];
@@ -328,16 +559,15 @@ accum_monte_carlo_weight
size_t nfailures = 0;
float axis_x[3], axis_y[3], axis_z[3];
float plane_size[2];
- float ray_org[3];
+ float ray_org[2][3];
float org[3];
float x[3], y[3];
- float rcp_pdf;
+ float plane_area;
res_T res = RES_OK;
- ASSERT(scn && rng && props && wavelengths && nwavelengths && angles && nangles);
- ASSERT(accums && basis && pos && lower && upper);
+ ASSERT(dev && ctx && basis && pos && lower && upper);
f2_sub(plane_size, upper, lower); /* In micron */
- rcp_pdf = plane_size[0] * plane_size[1];
+ plane_area = plane_size[0] * plane_size[1];
/* Define the projection axis */
f3_mulf(axis_x, basis + 0, plane_size[0] * 0.5f);
@@ -348,36 +578,57 @@ accum_monte_carlo_weight
f3_add(org, pos, f3_mulf(org, axis_z, lower[2]));
do {
- double path_length;
+ double path_length[2];
/* Uniformly sample a position onto the projection plane and use it as the
- * origin of the ray to trace */
- sample[0] = ssp_rng_uniform_double(rng, -1.0, 1.0);
- sample[1] = ssp_rng_uniform_double(rng, -1.0, 1.0);
+ * origin of the 1st ray to trace */
+ sample[0] = ssp_rng_uniform_double(ctx->rng, -1.0, 1.0);
+ sample[1] = ssp_rng_uniform_double(ctx->rng, -1.0, 1.0);
f3_mulf(x, axis_x, (float)sample[0]);
f3_mulf(y, axis_y, (float)sample[1]);
- f3_add(ray_org, f3_add(ray_org, x, y), org);
+ f3_add(ray_org[0], f3_add(ray_org[0], x, y), org);
- S3D(scene_trace_ray(scn, ray_org, axis_z, range, &hit));
-
- if(S3D_HIT_NONE(&hit)) /* NULL cross section and phase function weight */
+ S3D(scene_trace_ray(ctx->scene, ray_org[0], axis_z, range, &hit));
+ /* NULL cross section and differential cross section weight */
+ if(S3D_HIT_NONE(&hit))
break;
-
- res = compute_path_length(dev, scn, &hit, ray_org, axis_z, &path_length);
+ res = compute_path_length
+ (dev, ctx->scene, &hit, ray_org[0], axis_z, &path_length[0]);
if(res != RES_OK) { /* Handle numerical/geometry issues */
++nfailures;
continue;
}
- accum_cross_section
- (props, wavelengths, nwavelengths, rcp_pdf, path_length, accums);
+ /* Compute and accumulate the cross section weight */
+ accum_cross_section(ctx, plane_area, path_length[0]);
- /* TODO phase function integration */
- (void)angles, (void)nangles;
+ /* Uniformly sample a position onto the projection plane and use it as the
+ * origin of the 2nd ray to trace */
+ sample[0] = ssp_rng_uniform_double(ctx->rng, -1.0, 1.0);
+ sample[1] = ssp_rng_uniform_double(ctx->rng, -1.0, 1.0);
+ f3_mulf(x, axis_x, (float)sample[0]);
+ f3_mulf(y, axis_y, (float)sample[1]);
+ f3_add(ray_org[1], f3_add(ray_org[1], x, y), org);
+
+ S3D(scene_trace_ray(ctx->scene, ray_org[1], axis_z, range, &hit));
+ 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]);
+ if(res != RES_OK) { /* Handle numerical/geometry issues */
+ ++nfailures;
+ continue;
+ }
+ /* Compute and accumulate the per scattering angle differential cross
+ * section weight */
+ accum_differential_cross_section(ctx, plane_area, path_length, ray_org);
} while(res != RES_OK && nfailures < MAX_FAILURES);
- if(res != RES_OK) {
+ if(nfailures < MAX_FAILURES) {
+ res = RES_OK;
+ } else {
log_error(dev,
"Too many failures in computing the radiative path length. The sampled geometry\n"
"might not defined a closed volume.\n");
@@ -388,17 +639,9 @@ accum_monte_carlo_weight
static res_T
radiative_properties
(struct sschiff_device* dev,
- struct ssp_rng* rng,
const int istep,
- const double* wavelengths, /* Sorted in ascending order */
- const size_t nwavelengths,
- const double* angles, /* Sorted in ascending order in [0, PI] */
- const size_t nangles,
const size_t ndirs,
- struct s3d_scene* scene,
- const struct sschiff_material_properties* props,
- struct mc_accum* mc_accums,
- struct accum* accums)
+ struct integrator_context* ctx)
{
float lower[3], upper[3];
float aabb_pt[8][3];
@@ -407,11 +650,10 @@ radiative_properties
size_t iwlen;
int i;
res_T res = RES_OK;
- ASSERT(dev && rng && wavelengths && nwavelengths && ndirs && props);
- ASSERT(mc_accums && accums);
+ ASSERT(dev && ndirs && ctx);
(void)istep;
- S3D(scene_get_aabb(scene, lower, upper));
+ S3D(scene_get_aabb(ctx->scene, lower, upper));
/* AABB vertex layout
* 6-------7
@@ -431,14 +673,27 @@ radiative_properties
f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f);
- memset(accums, 0, sizeof(struct accum)*nwavelengths);
+ /* Clear direction accumulators */
+ FOR_EACH(iwlen, 0, ctx->nwavelengths) {
+ ctx->accums[iwlen].extinction_cross_section = 0;
+ ctx->accums[iwlen].absorption_cross_section = 0;
+ ctx->accums[iwlen].scattering_cross_section = 0;
+ ctx->accums[iwlen].avg_projected_area = 0;
+
+ /* Clean up to "limit_angle" since great angles are analytically computed */
+ memset(ctx->accums[iwlen].differential_cross_section, 0,
+ sizeof(double)*ctx->geom_ctx->limit_angles[iwlen]);
+ memset(ctx->accums[iwlen].differential_cross_section_cumulative, 0,
+ sizeof(double)*ctx->geom_ctx->limit_angles[iwlen]);
+ }
+
FOR_EACH(idir, 0, ndirs) {
float dir[4];
float basis[9];
float transform[12];
float pt[8][3];
- ssp_ran_sphere_uniform(rng, dir);
+ ssp_ran_sphere_uniform(ctx->rng, dir);
/* Build the transformation matrix from world to footprint space. Use the
* AABB centroid as the origin of the footprint space. */
@@ -457,8 +712,7 @@ radiative_properties
f3_max(upper, upper, pt[i]);
}
- res = accum_monte_carlo_weight(dev, scene, rng, props, wavelengths,
- nwavelengths, angles, nangles, basis, centroid, lower, upper, accums);
+ res = accum_monte_carlo_weight(dev, ctx, basis, centroid, lower, upper);
if(res != RES_OK) goto error;
#if 0
@@ -470,18 +724,42 @@ radiative_properties
}
#endif
}
- FOR_EACH(iwlen, 0, nwavelengths) {
- int iweight;
- FOR_EACH(iweight, 0, SSCHIFF_DATA_COUNT__) {
- const double w = accums[iwlen].weights[iweight] / (double)ndirs;
- mc_accums[iwlen].mc_data[iweight].weight += w;
- mc_accums[iwlen].mc_data[iweight].sqr_weight += w * w;
+ /* Compute the Monte Carlo weight of the temporary accumulator */
+ FOR_EACH(iwlen, 0, ctx->nwavelengths) {
+ size_t iangle;
+
+ #define MC_ACCUM(Data) { \
+ const double w = ctx->accums[iwlen].Data / (double)ndirs; \
+ ctx->mc_accums[iwlen].Data.weight += w; \
+ ctx->mc_accums[iwlen].Data.sqr_weight += w*w; \
+ } (void)0
+
+ MC_ACCUM(extinction_cross_section);
+ MC_ACCUM(absorption_cross_section);
+ MC_ACCUM(scattering_cross_section);
+ MC_ACCUM(avg_projected_area);
+
+ /* Accum up to "limit angle" since great angles are analitically computed */
+ FOR_EACH(iangle, 0, ctx->geom_ctx->limit_angles[iwlen]) {
+ MC_ACCUM(differential_cross_section[iangle]);
+ MC_ACCUM(differential_cross_section_cumulative[iangle]);
}
+
+ #undef MC_ACCUM
}
exit:
return res;
error:
- memset(mc_accums, 0, sizeof(struct mc_accum)*nwavelengths);
+ FOR_EACH(iwlen, 0, ctx->nwavelengths) {
+ ctx->mc_accums[iwlen].extinction_cross_section = MC_DATA_NULL;
+ ctx->mc_accums[iwlen].absorption_cross_section = MC_DATA_NULL;
+ ctx->mc_accums[iwlen].scattering_cross_section = MC_DATA_NULL;
+ ctx->mc_accums[iwlen].avg_projected_area = MC_DATA_NULL;
+ memset(ctx->mc_accums[iwlen].differential_cross_section, 0,
+ sizeof(double)*ctx->nangles);
+ memset(ctx->mc_accums[iwlen].differential_cross_section_cumulative, 0,
+ sizeof(double)*ctx->nangles);
+ }
goto exit;
}
@@ -489,7 +767,7 @@ static void
get_mc_value
(const struct mc_data* data,
const size_t nrealisations,
- struct sschiff_estimator_value* value)
+ struct sschiff_state* value)
{
ASSERT(data && value);
value->E = data->weight / (double)nrealisations;
@@ -502,7 +780,208 @@ static char
check_distribution(struct sschiff_geometry_distribution* distrib)
{
ASSERT(distrib);
- return distrib->sample != NULL;
+ return distrib->sample != NULL && distrib->caracteristic_length > 0;
+}
+
+static void
+geometry_context_release(struct geometry_context* ctx)
+{
+ ASSERT(ctx);
+ if(ctx->limit_angles) sa_release(ctx->limit_angles);
+ if(ctx->medium_refractive_indices) sa_release(ctx->medium_refractive_indices);
+ if(ctx->properties) sa_release(ctx->properties);
+}
+
+static res_T
+geometry_context_init
+ (struct sschiff_device* dev,
+ struct sschiff_geometry_distribution* distrib,
+ const double* wavelengths, /* In the vacuum. Sorted in ascending order */
+ const size_t nwavelengths,
+ const double* angles, /* In radian. Sorted in ascending order in [0, PI] */
+ const size_t nangles,
+ struct geometry_context* ctx)
+{
+ double rcp_PI_Lc;
+ size_t iwlen;
+ res_T res = RES_OK;
+ ASSERT(dev && distrib && wavelengths && nwavelengths && angles && nangles && ctx);
+ memset(ctx, 0, sizeof(struct geometry_context));
+
+ /* Precomputed valut */
+ rcp_PI_Lc = 1.0 / (PI*distrib->caracteristic_length);
+
+ if(!sa_add(ctx->limit_angles, nwavelengths)) {
+ log_error(dev,
+"Couldn't allocate the per wavelength indices of the limit scattering angle.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ if(!sa_add(ctx->medium_refractive_indices, nwavelengths)) {
+ log_error(dev,
+"Couldn't allocate the per wavelength refractive indices of the medium.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ if(!sa_add(ctx->properties, nwavelengths)) {
+ log_error(dev, "Couldn't allocate the per wavelength optical properties.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ FOR_EACH(iwlen, 0, nwavelengths) {
+ double lambda_e;
+ double theta_l;
+ double* angle;
+
+ /* Fetch the particle optical properties */
+ distrib->material.get_property
+ (distrib->material.material, wavelengths[iwlen], ctx->properties + iwlen);
+
+ /* Precompute the refractive index in the medium */
+ lambda_e = wavelengths[iwlen] / ctx->properties[iwlen].medium_refractive_index;
+ ctx->medium_refractive_indices[iwlen] = 2.0*PI / lambda_e;
+
+ /* Search the limit angle */
+ theta_l = sqrt(lambda_e * rcp_PI_Lc);
+ angle = search_lower_bound
+ (&theta_l, angles, nangles, sizeof(double), cmp_double);
+ if(!angle || eq_eps(*angle, PI, 1.e-6)) {
+ log_error(dev, "Invalid limit angle \"%g\".\n", *angle);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Define the index of the angle that *follows* the limit angle */
+ ctx->limit_angles[iwlen] = (size_t)(angle - angles) + 1;
+ ASSERT(ctx->limit_angles[iwlen] < nangles);
+ }
+
+ ctx->angles = angles;
+
+exit:
+ return res;
+error:
+ geometry_context_release(ctx);
+ goto exit;
+}
+
+static void
+integrator_context_release(struct integrator_context* ctx)
+{
+ size_t iwlen;
+ ASSERT(ctx);
+
+ if(ctx->rng) SSP(rng_ref_put(ctx->rng));
+ if(ctx->shape) S3D(shape_ref_put(ctx->shape));
+ if(ctx->scene) S3D(scene_ref_put(ctx->scene));
+
+ #define CTXFREE(Data) if(Data) sa_release(Data)
+ if(ctx->accums) {
+ FOR_EACH(iwlen, 0, ctx->nwavelengths) {
+ CTXFREE(ctx->accums[iwlen].differential_cross_section);
+ CTXFREE(ctx->accums[iwlen].differential_cross_section_cumulative);
+ }
+ sa_release(ctx->accums);
+ }
+ if(ctx->mc_accums) {
+ FOR_EACH(iwlen, 0, ctx->nwavelengths) {
+ CTXFREE(ctx->mc_accums[iwlen].differential_cross_section);
+ CTXFREE(ctx->mc_accums[iwlen].differential_cross_section_cumulative);
+ }
+ sa_release(ctx->mc_accums);
+ }
+ #undef CTXFREE
+}
+
+static res_T
+integrator_context_init
+ (struct integrator_context* ctx,
+ struct s3d_device* s3d,
+ struct ssp_rng* rng,
+ const struct geometry_context* geom_ctx,
+ const size_t nwavelengths,
+ const size_t nangles)
+{
+ size_t iwlen;
+ res_T res = RES_OK;
+ ASSERT(ctx && s3d && rng && geom_ctx && nwavelengths && nangles >= 3);
+
+ memset(ctx, 0, sizeof(struct integrator_context));
+
+ if(RES_OK!=(res = s3d_shape_create_mesh(s3d, &ctx->shape))) goto error;
+ if(RES_OK!=(res = s3d_scene_create(s3d, &ctx->scene))) goto error;
+ if(RES_OK!=(res = s3d_scene_attach_shape(ctx->scene, ctx->shape))) goto error;
+
+ #define CTXALLOC(Data, N) { \
+ if(!sa_add((Data), (N))) { \
+ res = RES_MEM_ERR; \
+ goto error; \
+ } \
+ memset((Data), 0, sizeof((Data)[0])*(N)); \
+ } (void)0
+ CTXALLOC(ctx->accums, nwavelengths);
+ FOR_EACH(iwlen, 0, nwavelengths) {
+ CTXALLOC(ctx->accums[iwlen].differential_cross_section, nangles);
+ CTXALLOC(ctx->accums[iwlen].differential_cross_section_cumulative, nangles);
+ }
+ CTXALLOC(ctx->mc_accums, nwavelengths);
+ FOR_EACH(iwlen, 0, nwavelengths) {
+ CTXALLOC(ctx->mc_accums[iwlen].differential_cross_section, nangles);
+ CTXALLOC(ctx->mc_accums[iwlen].differential_cross_section_cumulative, nangles);
+ }
+ #undef CTXALLOC
+
+ SSP(rng_ref_get(rng));
+ ctx->rng = rng;
+ ctx->geom_ctx = geom_ctx;
+ ctx->nwavelengths = nwavelengths;
+ ctx->nangles = nangles;
+
+exit:
+ return res;
+error:
+ integrator_context_release(ctx);
+ goto exit;
+}
+
+static res_T
+begin_realisation
+ (struct sschiff_device* dev,
+ struct sschiff_geometry_distribution* distrib,
+ struct integrator_context* ctx)
+{
+ res_T res = RES_OK;
+ ASSERT(dev && distrib && ctx);
+
+ /* Sample a particle */
+ res = distrib->sample(ctx->rng, ctx->shape, distrib->context);
+ if(res != RES_OK) {
+ log_error(dev, "Error in sampling a particle.\n");
+ goto error;
+ }
+
+ /* Build the Star-3D representation of the sampled shape */
+ S3D(scene_begin_session(ctx->scene, S3D_TRACE));
+
+exit:
+ return res;
+error:
+ {
+ /* Disable active session if necessary */
+ int mask;
+ if(S3D(scene_get_session_mask(ctx->scene, &mask)), mask) {
+ S3D(scene_end_session(ctx->scene));
+ }
+ }
+ goto exit;
+}
+
+static FINLINE void
+end_realisation(struct integrator_context* ctx)
+{
+ ASSERT(ctx);
+ S3D(scene_end_session(ctx->scene));
}
static void
@@ -514,7 +993,18 @@ estimator_release(ref_T* ref)
estimator = CONTAINER_OF(ref, struct sschiff_estimator, ref);
dev = estimator->dev;
if(estimator->wavelengths) sa_release(estimator->wavelengths);
- if(estimator->accums) sa_release(estimator->accums);
+ if(estimator->angles) sa_release(estimator->angles);
+ if(estimator->accums) {
+ const size_t nwavelengths = sa_size(estimator->accums);
+ size_t i;
+ FOR_EACH(i, 0, nwavelengths) {
+ if(estimator->accums[i].differential_cross_section)
+ sa_release(estimator->accums[i].differential_cross_section);
+ if(estimator->accums[i].differential_cross_section_cumulative)
+ sa_release(estimator->accums[i].differential_cross_section_cumulative);
+ }
+ sa_release(estimator->accums);
+ }
MEM_RM(dev->allocator, estimator);
SSCHIFF(device_ref_put(dev));
}
@@ -522,17 +1012,22 @@ estimator_release(ref_T* ref)
static res_T
estimator_create
(struct sschiff_device* dev,
+ const double* wavelengths,
const size_t nwavelengths,
+ sschiff_scattering_angles_distribution_T angles_distrib,
+ const size_t nangles,
struct sschiff_estimator** out_estimator)
{
struct sschiff_estimator* estimator = NULL;
+ size_t i;
res_T res = RES_OK;
- ASSERT(dev && out_estimator);
+ ASSERT(dev && wavelengths && nwavelengths && angles_distrib && nangles);
+ ASSERT(out_estimator);
estimator = MEM_CALLOC(dev->allocator, 1, sizeof(struct sschiff_estimator));
if(!estimator) {
log_error
- (dev, "Not enough memory: couldn't allocate the Star-Schiff estimator.\n");
+ (dev, "Couldn't allocate the Star-Schiff estimator.\n");
res = RES_MEM_ERR;
goto error;
}
@@ -540,19 +1035,63 @@ estimator_create
SSCHIFF(device_ref_get(dev));
estimator->dev = dev;
+ /* Setup and sort the wavelengths */
if(!sa_add(estimator->wavelengths, nwavelengths)) {
log_error(dev,
- "Not enough memory: couldn't allocate the list of estimated wavelengths.\n");
+ "Couldn't allocate the list of estimated wavelengths.\n");
res = RES_MEM_ERR;
goto error;
}
- if(!sa_add(estimator->accums, nwavelengths)) {
+ FOR_EACH(i, 0, nwavelengths) estimator->wavelengths[i] = wavelengths[i];
+ qsort(estimator->wavelengths, nwavelengths, sizeof(double), cmp_double);
+
+ /* Generate the scattering angles */
+ if(!sa_add(estimator->angles, nangles)) {
log_error(dev,
- "Not enough memory: couldn't allocate the Monte Carlo accumulator.\n");
+ "Couldn't allocate the list of scattering angles.\n");
res = RES_MEM_ERR;
goto error;
}
+ angles_distrib(estimator->angles, nangles);
+ if(estimator->angles[0] != 0.f
+ || !eq_eps(estimator->angles[nangles-1], PI, 1.e-6)) {
+ log_error(dev, "Invalid scattering angle distribution.\n");
+ log_error(dev, "The first and last angles must be 0 and PI, respectively.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ FOR_EACH(i, 1, nangles) {
+ if(estimator->angles[i] <= estimator->angles[i-1]) {
+ log_error(dev, "Invalid scattering angle distribution.\n");
+ log_error(dev, "Angles must be sorted in ascending order in [0, PI]\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ /* Allocate the estimator accumulators */
+ if(!sa_add(estimator->accums, nwavelengths)) {
+ log_error(dev,
+ "Couldn't allocate the Monte Carlo accumulator of the estimator.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ memset(estimator->accums, 0, sizeof(estimator->accums[0])*nwavelengths);
+ FOR_EACH(i, 0, nwavelengths) {
+ /* FIXME rename the following estimation */
+ if(!sa_add(estimator->accums[i].differential_cross_section, nangles)) {
+ log_error(dev,
+ "Couldn't allocate the differential cross sections to estimate.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ if(!sa_add(estimator->accums[i].differential_cross_section_cumulative, nangles)) {
+ log_error(dev,
+"Couldn't allocate the cumulative differential cross sections to estimate.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ }
exit:
*out_estimator = estimator;
return res;
@@ -580,30 +1119,34 @@ sschiff_integrate
const size_t ndirs,
struct sschiff_estimator** out_estimator)
{
- struct sschiff_material_properties** mtls = NULL;
- struct accum** accums = NULL;
- struct mc_accum** mc_accums = NULL;
+ struct geometry_context geom_ctx;
+ struct integrator_context* ctxs = NULL;
struct sschiff_estimator* estimator = NULL;
- struct s3d_shape** shapes = NULL;
- struct s3d_scene** scenes = NULL;
struct ssp_rng** rngs = NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
- double* angles = NULL;
size_t i;
int igeom;
ATOMIC res = (ATOMIC)RES_OK;
- (void)nangles;
+ memset(&geom_ctx, 0, sizeof(geom_ctx));
if(!dev || !rng || !distrib || !wavelengths || !nwavelengths
- || !angles_distrib || nangles < 2 || !ngeoms || !ndirs || !out_estimator
+ || !angles_distrib || nangles < 3 || !ngeoms || !ndirs || !out_estimator
|| !check_distribution(distrib)) {
return RES_BAD_ARG;
}
- res = estimator_create(dev, nwavelengths, &estimator);
+ /* Create the Schiff estimator */
+ res = estimator_create
+ (dev, wavelengths, nwavelengths, angles_distrib, nangles, &estimator);
if(res != RES_OK) goto error;
estimator->nrealisations = ngeoms;
+ /* Initialise the geometry context */
+ res = geometry_context_init(dev, distrib, estimator->wavelengths,
+ nwavelengths, estimator->angles, nangles, &geom_ctx);
+ if(res != RES_OK) goto error;
+
+ /* Create a RNG proxy from the submitted RNG state */
res = ssp_rng_proxy_create_from_rng
(dev->allocator, rng, dev->nthreads, &rng_proxy);
if(res != RES_OK) {
@@ -611,152 +1154,58 @@ sschiff_integrate
goto error;
}
- /* Setup the scattering angles */
- angles = sa_add(angles, nangles);
- if(!angles) {
- log_error(dev,
- "Not enough memory: couldn't allocate the list of scattering angles.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- angles_distrib(angles, nangles);
- if(angles[0] != 0.f && !eq_eps(angles[1], PI, 1.e-6)) {
- log_error(dev,
-"Invalid scattering angles distribution. The first and the last angles must be\n"
-"0 and PI, respectively\n");
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Create the containers of per thread data structures */
- if(!sa_add(mtls, dev->nthreads)) {
- log_error(dev,
- "No enough memory: couldn't allocate the list of optical properties.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- if(!sa_add(accums, dev->nthreads)) {
- log_error(dev,
- "No enough memory: couldn't allocate the list of temporary accumulators.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- if(!sa_add(mc_accums, dev->nthreads)) {
- log_error(dev,
- "Not enough memory: couldn't allocate the list Monte Carlo accumulators.\n");
+ /* Create per thread data structures */
+ if(!sa_add(ctxs, dev->nthreads)) {
+ log_error(dev, "Couldn't allocate the integrator contexts.\n");
res = RES_MEM_ERR;
goto error;
}
if(!sa_add(rngs, dev->nthreads)) {
- log_error(dev,
- "Not enough memory: couldn't allocate the list of RNGs.\n");
- }
- if(!sa_add(scenes, dev->nthreads)) {
- log_error(dev,
- "Not enough memory: couldn't allocate the list of Star-3D scenes.\n");
+ log_error(dev, "Couldn't allocate the list of RNGs.\n");
res = RES_MEM_ERR;
goto error;
}
- if(!sa_add(shapes, dev->nthreads)) {
- log_error(dev,
- "Not enough memory: couldn't allocate the list of Star-3D shapes.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- memset(mtls, 0, sizeof(struct sschiff_material_properties*)*dev->nthreads);
- memset(accums, 0, sizeof(struct accum*)*dev->nthreads);
- memset(mc_accums, 0, sizeof(struct mc_accum*)*dev->nthreads);
- memset(rngs, 0, sizeof(struct ssp_rng*)*dev->nthreads);
- memset(scenes, 0, sizeof(struct s3d_scene*)*dev->nthreads);
- memset(shapes, 0, sizeof(struct s3d_shape*)*dev->nthreads);
+ memset(ctxs, 0, sizeof(ctxs[0])*dev->nthreads);
+ memset(rngs, 0, sizeof(rngs[0])*dev->nthreads);
- /* Create the per thread data structures */
+ /* Init the per thread data structures */
FOR_EACH(i, 0, dev->nthreads) {
- if(!sa_add(mtls[i], nwavelengths)) {
- log_error(dev,
- "Not enough memory: couldn't allocate the optical properties.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- if(!sa_add(accums[i], nwavelengths)) {
- log_error(dev,
- "Not enough memory: couldn't allocate the temporary accumulator.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- if(!sa_add(mc_accums[i], nwavelengths)) {
- log_error(dev,
- "Not enough memory: couldn't allocate the Monte Carlo accumulator.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- memset(mc_accums[i], 0, sizeof(struct mc_accum)*nwavelengths);
-
res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs + i);
if(res != RES_OK) {
- log_error(dev, "Couldn't create the RNG.\n");
- goto error;
- }
- res = s3d_shape_create_mesh(dev->s3d[i], shapes + i);
- if(res != RES_OK) {
- log_error(dev, "Couldn't create the Star-3D shape.\n");
- goto error;
- }
- res = s3d_scene_create(dev->s3d[i], scenes + i);
- if(res != RES_OK) {
- log_error(dev, "Couldn't create the Star-3D scene.\n");
+ log_error(dev,
+ "Couldn't create the RNG of the thread \"%lu\".\n", (unsigned long)i);
goto error;
}
- res = s3d_scene_attach_shape(scenes[i], shapes[i]);
+ res = integrator_context_init(ctxs + i, dev->s3d[i], rngs[i], &geom_ctx,
+ nwavelengths, nangles);
if(res != RES_OK) {
log_error(dev,
- "Couldn't attach the Star-3D Schiff shape to the Star-3D Schiff scene.\n");
+ "Couldn't initialise the integrator context of the thread \"%lu\".\n",
+ (unsigned long)i);
goto error;
}
}
- /* Clear the accumulators */
- memset(estimator->accums, 0, sizeof(struct mc_accum)*nwavelengths);
-
- /* Setup and sort the wavelengths */
- FOR_EACH(i, 0, nwavelengths) estimator->wavelengths[i] = wavelengths[i];
- qsort(estimator->wavelengths, nwavelengths, sizeof(double), cmp_double);
-
/* Paralell Schiff integration */
#pragma omp parallel for schedule(static)
for(igeom=0; igeom < (int)ngeoms; ++igeom) {
- const double* wlengths = estimator->wavelengths;
- struct sschiff_material material = SSCHIFF_NULL_MATERIAL;
- size_t iwavelength;
const int ithread = omp_get_thread_num();
- ATOMIC res_local;
+ ATOMIC res_local = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue;
- /* Sample a geometry, i.e. a shape and its associated material */
- res_local = distrib->sample
- (rngs[ithread], &material, shapes[ithread], distrib->context);
+ /* Setup the data for the current realisation */
+ res_local = begin_realisation(dev, distrib, ctxs + ithread);
if(res_local != RES_OK) {
- log_error(dev, "Error in sampling a Schiff geometry.\n");
- ATOMIC_SET(&res, res_local);
+ ATOMIC_CAS(&res, res_local, RES_OK);
continue;
}
- /* Fetch the per wavelength material properties */
- FOR_EACH(iwavelength, 0, nwavelengths) {
- material.get_property(material.material, wavelengths[iwavelength],
- mtls[ithread] + iwavelength);
- }
-
- /* Build the Star-3D representation of the sampled shape */
- S3D(scene_begin_session(scenes[ithread], S3D_TRACE));
/* Schiff Estimation */
- res_local = radiative_properties(dev, rngs[ithread], igeom, wlengths,
- nwavelengths, angles, nangles, ndirs, scenes[ithread], mtls[ithread],
- mc_accums[ithread], accums[ithread]);
- if(res != RES_OK) ATOMIC_SET(&res, res_local);
+ res_local = radiative_properties(dev, igeom, ndirs, ctxs+ithread);
+ if(res_local != RES_OK) ATOMIC_CAS(&res, res_local, RES_OK);
- S3D(scene_end_session(scenes[ithread]));
+ end_realisation(ctxs + ithread);
}
if(res != RES_OK) goto error; /* Handle integration error */
@@ -764,44 +1213,36 @@ sschiff_integrate
FOR_EACH(i, 0, dev->nthreads) {
size_t iwlen;
FOR_EACH(iwlen, 0, nwavelengths) {
- size_t idata;
- FOR_EACH(idata, 0, SSCHIFF_DATA_COUNT__) {
- estimator->accums[iwlen].mc_data[idata].weight +=
- mc_accums[i][iwlen].mc_data[idata].weight;
- estimator->accums[iwlen].mc_data[idata].sqr_weight +=
- mc_accums[i][iwlen].mc_data[idata].sqr_weight;
+ size_t iangle;
+ #define MC_ACCUM(Data) { \
+ const struct mc_accum* mc_accum = ctxs[i].mc_accums + iwlen; \
+ estimator->accums[iwlen].Data.weight += mc_accum->Data.weight; \
+ estimator->accums[iwlen].Data.sqr_weight += mc_accum->Data.sqr_weight; \
+ } (void)0
+ MC_ACCUM(extinction_cross_section);
+ MC_ACCUM(absorption_cross_section);
+ MC_ACCUM(scattering_cross_section);
+ MC_ACCUM(avg_projected_area);
+
+ /* Accum up to "limit angle"; Great angles are analitically computed */
+ FOR_EACH(iangle, 0, geom_ctx.limit_angles[iwlen]) {
+ MC_ACCUM(differential_cross_section[iangle]);
+ MC_ACCUM(differential_cross_section_cumulative[iangle]);
}
}
}
exit:
+ geometry_context_release(&geom_ctx);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
- if(angles) sa_release(angles);
- /* Release per thread data */
- if(mtls) {
- FOR_EACH(i, 0, dev->nthreads) if(mtls[i]) sa_release(mtls[i]);
- sa_release(mtls);
- }
- if(accums) {
- FOR_EACH(i, 0, dev->nthreads) if(accums[i]) sa_release(accums[i]);
- sa_release(accums);
- }
- if(mc_accums) {
- FOR_EACH(i, 0, dev->nthreads) if(mc_accums[i]) sa_release(mc_accums[i]);
- sa_release(mc_accums);
+ if(ctxs) {
+ FOR_EACH(i, 0, dev->nthreads) integrator_context_release(ctxs+i);
+ sa_release(ctxs);
}
if(rngs) {
FOR_EACH(i, 0, dev->nthreads) if(rngs[i]) SSP(rng_ref_put(rngs[i]));
sa_release(rngs);
}
- if(shapes) {
- FOR_EACH(i, 0, dev->nthreads) if(shapes[i]) S3D(shape_ref_put(shapes[i]));
- sa_release(shapes);
- }
- if(scenes) {
- FOR_EACH(i, 0, dev->nthreads) if(scenes[i]) S3D(scene_ref_put(scenes[i]));
- sa_release(scenes);
- }
if(out_estimator) *out_estimator = estimator;
return (res_T)res;
error:
@@ -847,19 +1288,19 @@ sschiff_estimator_get_realisations_count
}
res_T
-sschiff_estimator_get_wavelength_state
+sschiff_estimator_get_wavelength_cross_section
(const struct sschiff_estimator* estimator,
const double wavelength,
- struct sschiff_estimator_state* state)
+ struct sschiff_cross_section* cross_section)
{
- const struct mc_accum* accum;
+ const struct mc_accum* acc;
const double* wavelengths;
const double* find;
size_t iwavelength;
size_t nwavelengths;
- int i;
+ size_t N;
- if(!estimator || !state)
+ if(!estimator || !cross_section)
return RES_BAD_ARG;
wavelengths = estimator->wavelengths;
@@ -869,39 +1310,36 @@ sschiff_estimator_get_wavelength_state
if(!find) return RES_BAD_ARG;
iwavelength = (size_t)(find - wavelengths);
- accum = estimator->accums + iwavelength;
-
- FOR_EACH(i, 0, SSCHIFF_DATA_COUNT__) {
- get_mc_value(accum->mc_data+i, estimator->nrealisations, state->values+i);
- }
- state->wavelength = wavelength;
+ acc = estimator->accums + iwavelength;
+
+ N = estimator->nrealisations;
+ get_mc_value(&acc->extinction_cross_section, N, &cross_section->extinction);
+ get_mc_value(&acc->absorption_cross_section, N, &cross_section->absorption);
+ get_mc_value(&acc->scattering_cross_section, N, &cross_section->scattering);
+ get_mc_value(&acc->avg_projected_area, N,
+ &cross_section->average_projected_area);
return RES_OK;
}
res_T
-sschiff_estimator_get_states
+sschiff_estimator_get_cross_sections
(const struct sschiff_estimator* estimator,
- struct sschiff_estimator_state* states)
+ struct sschiff_cross_section* cross_sections)
{
- const double* wavelengths;
- const struct mc_accum* accums;
size_t nwavelengths;
size_t iwlen;
- int i;
- if(!estimator || !states) return RES_BAD_ARG;
+ size_t N;
+ if(!estimator || !cross_sections) return RES_BAD_ARG;
nwavelengths = sa_size(estimator->wavelengths);
- wavelengths = estimator->wavelengths;
- accums = estimator->accums;
-
+ N = estimator->nrealisations;
FOR_EACH(iwlen, 0, nwavelengths) {
- states[iwlen].wavelength = wavelengths[iwlen];
- FOR_EACH(i, 0, SSCHIFF_DATA_COUNT__) {
- get_mc_value
- (accums[iwlen].mc_data + i,
- estimator->nrealisations,
- states[iwlen].values + i);
- }
+ const struct mc_accum* acc = estimator->accums + iwlen;
+ get_mc_value(&acc->extinction_cross_section, N, &cross_sections->extinction);
+ get_mc_value(&acc->absorption_cross_section, N, &cross_sections->absorption);
+ get_mc_value(&acc->scattering_cross_section, N, &cross_sections->scattering);
+ get_mc_value(&acc->avg_projected_area, N,
+ &cross_sections->average_projected_area);
}
return RES_OK;
}
diff --git a/src/test_sschiff_estimator_cylinder.c b/src/test_sschiff_estimator_cylinder.c
@@ -128,10 +128,7 @@ dump_geometry(struct geometry* geom)
static res_T
sample_cylinder
- (struct ssp_rng* rng,
- struct sschiff_material* mtl,
- struct s3d_shape* shape,
- void* sampler_context)
+ (struct ssp_rng* rng, struct s3d_shape* shape, void* sampler_context)
{
struct s3d_vertex_data attrib;
struct sampler_context* sampler_ctx = sampler_context;
@@ -149,9 +146,6 @@ sample_cylinder
attrib.type = S3D_FLOAT3;
attrib.get = get_position;
- mtl->get_property = get_material_property;
- mtl->material = sampler_ctx;
-
nverts = sa_size(cylinder.geometry->vertices) / 3/*#coords*/;
nprims = sa_size(cylinder.geometry->indices) / 3/*#indices per prim*/;
@@ -559,7 +553,7 @@ main(int argc, char** argv)
2.010619438e+2, 2.018997018e+2, 2.027374599e+2
};
const size_t nx = sizeof(x)/sizeof(double);
- const size_t nscatt_angles = 2;
+ const size_t nscatt_angles = 3;
const size_t ngeoms = 100;
const size_t ndirs = 10;
size_t i;
@@ -575,14 +569,17 @@ main(int argc, char** argv)
FOR_EACH(i, 0, nx) {
const double wavelength = 0.6; /* In micron */
- struct sschiff_estimator_state state;
+ struct sschiff_cross_section cross_section;
double interval[2];
- struct sschiff_estimator_value* val;
+ struct sschiff_state* val;
sampler_ctx.aspect_ratio = 0.837;
sampler_ctx.mean_radius = (x[i] * 0.450) / (2*PI);
sampler_ctx.sigma = 1.18;
+ distrib.caracteristic_length = 1;
+ distrib.material.get_property = get_material_property;
+ distrib.material.material = &sampler_ctx;
distrib.sample = sample_cylinder;
distrib.context = &sampler_ctx;
@@ -594,14 +591,13 @@ main(int argc, char** argv)
time_sub(&t0, &t1, &t0);
time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
- CHECK(sschiff_estimator_get_wavelength_state
- (estimator, wavelength, &state), RES_OK);
- CHECK(eq_eps(state.wavelength, wavelength, 1.e-6), 1);
+ CHECK(sschiff_estimator_get_wavelength_cross_section
+ (estimator, wavelength, &cross_section), RES_OK);
printf("%u - x = %g - Wavelength = %g micron - %s\n",
(unsigned)i, x[i], wavelength, buf);
- val = state.values + SSCHIFF_EXTINCTION_CROSS_SECTION;
+ val = &cross_section.extinction;
compute_intersection(interval, results[i].extinction_E,
results[i].extinction_SE, val->E, val->SE);
printf(" Extinction ~ %9.3g +/- %9.3g ~ %9.3g +/- %9.3g (%9.3g)\n",
@@ -609,7 +605,7 @@ main(int argc, char** argv)
val->E, val->SE, interval[1] - interval[0]);
CHECK(interval[0] <= interval[1], 1);
- val = state.values + SSCHIFF_ABSORPTION_CROSS_SECTION;
+ val = &cross_section.absorption;
compute_intersection(interval, results[i].absorption_E,
results[i].absorption_SE, val->E, val->SE);
printf(" Absorption ~ %9.3g +/- %9.3g ~ %9.3g +/- %9.3g (%9.3g)\n",
@@ -617,7 +613,7 @@ main(int argc, char** argv)
val->E, val->SE, interval[1] - interval[0]);
CHECK(interval[0] <= interval[1], 1);
- val = state.values + SSCHIFF_SCATTERING_CROSS_SECTION;
+ val = &cross_section.scattering;
compute_intersection(interval, results[i].scattering_E,
results[i].scattering_SE, val->E, val->SE);
printf(" Scattering ~ %9.3g +/- %9.3g ~ %9.3g +/- %9.3g (%9.3g)\n",
@@ -625,7 +621,7 @@ main(int argc, char** argv)
val->E, val->SE, interval[1] - interval[0]);
CHECK(interval[0] <= interval[1], 1);
- val = state.values + SSCHIFF_AVERAGE_PROJECTED_AREA;
+ val = &cross_section.average_projected_area;
compute_intersection(interval, results[i].avg_proj_area_E,
results[i].avg_proj_area_SE, val->E, val->SE);
printf(" Proj Area ~ %9.3g +/- %9.3g ~ %9.3g +/- %9.3g (%9.3g)\n",
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -111,7 +111,6 @@ get_position(const unsigned ivert, float vertex[3], void* ctx)
static res_T
sample_sphere
(struct ssp_rng* rng,
- struct sschiff_material* mtl,
struct s3d_shape* shape,
void* sampler_context)
{
@@ -121,7 +120,6 @@ sample_sphere
size_t nverts, nprims;
NCHECK(rng, NULL);
- NCHECK(mtl, NULL);
NCHECK(shape, NULL);
sphere.geometry = &sampler_ctx->geometry;
@@ -132,9 +130,6 @@ sample_sphere
attrib.type = S3D_FLOAT3;
attrib.get = get_position;
- mtl->get_property = get_material_property;
- mtl->material = sampler_ctx;
-
nverts = sa_size(sphere.geometry->vertices) / 3/*#coords*/;
nprims = sa_size(sphere.geometry->indices) / 3/*#indices per prim*/;
@@ -242,11 +237,11 @@ check_schiff_estimation
char buf[64];
struct sschiff_estimator* estimator;
struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
- struct sschiff_estimator_state state;
+ struct sschiff_cross_section cross_section;
struct time t0, t1;
size_t i;
- const size_t nscatt_angles = 2;
+ const size_t nscatt_angles = 3;
const size_t ngeoms = 100;
const size_t ndirs = 100;
const double wavelength = sampler_ctx->wavelength;
@@ -255,6 +250,9 @@ check_schiff_estimation
NCHECK(results, NULL);
NCHECK(results, 0);
+ distrib.material.get_property = get_material_property;
+ distrib.material.material = sampler_ctx;
+ distrib.caracteristic_length = 1.0;
distrib.sample = sample_sphere;
distrib.context = sampler_ctx;
@@ -265,44 +263,40 @@ check_schiff_estimation
sampler_ctx->medium_refractive_index);
FOR_EACH(i, 0, nresults) {
- const struct sschiff_estimator_value* val;
+ const struct sschiff_state* val;
double dst;
sampler_ctx->mean_radius = results[i].mean_radius;
time_current(&t0);
- CHECK(sschiff_integrate(dev, rng, &distrib, &wavelength, 1,
+ CHECK(sschiff_integrate(dev, rng, &distrib, &wavelength, 1,
sschiff_uniform_scattering_angles, 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.2f micron - %s: \n", (int)i, results[i].mean_radius, buf);
- CHECK(sschiff_estimator_get_states(estimator, &state), RES_OK);
- CHECK(eq_eps(state.wavelength, wavelength, 1.e-12), 1);
+ CHECK(sschiff_estimator_get_cross_sections(estimator, &cross_section), RES_OK);
- val = state.values + SSCHIFF_EXTINCTION_CROSS_SECTION ;
+ val = &cross_section.extinction;
dst = val->E - results[i].extinction_cross_section;
printf(" Extinction cross section = %7.2f ~ %12.7f +/- %12.7f (%6.2f)\n",
results[i].extinction_cross_section, val->E, val->SE, dst / val->SE);
CHECK(eq_eps(val->E, results[i].extinction_cross_section, 4*val->SE), 1);
- val = state.values + SSCHIFF_ABSORPTION_CROSS_SECTION ;
+ val = &cross_section.absorption;
dst = val->E - results[i].absorption_cross_section;
printf(" Absorption cross section = %7.2f ~ %12.7f +/- %12.7f (%6.2f)\n",
results[i].absorption_cross_section, val->E, val->SE, dst / val->SE);
CHECK(eq_eps(val->E, results[i].absorption_cross_section, 4*val->SE), 1);
- val = state.values + SSCHIFF_SCATTERING_CROSS_SECTION;
+ val = &cross_section.scattering;
dst = val->E - results[i].scattering_cross_section;
printf(" Scattering cross section = %7.2f ~ %12.7f +/- %12.7f (%6.2f)\n",
results[i].scattering_cross_section, val->E, val->SE, dst / val->SE);
CHECK(eq_eps(val->E, results[i].scattering_cross_section, 4*val->SE), 1);
- /*s = log(sampler_ctx->sigma);
- r = sampler_ctx->mean_radius * exp(2.5 * s * s);
- r = PI * r * r * 1.e-12 */ /* Meter^2 to micron^2*/;
- val = state.values + SSCHIFF_AVERAGE_PROJECTED_AREA;
+ val = &cross_section.average_projected_area;
printf(" Averavege projected area = %12.6g +/- %12.7g\n", val->E, val->SE);
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
@@ -317,7 +311,7 @@ main(int argc, char** argv)
struct sschiff_device* dev;
struct sschiff_geometry_distribution distrib = SSCHIFF_NULL_GEOMETRY_DISTRIBUTION;
struct sschiff_estimator* estimator;
- struct sschiff_estimator_state state;
+ struct sschiff_cross_section cross_section;
struct ssp_rng* rng;
const struct result results_n_r_1_01[] = {
{ 1.00, 0.484, 0.159, 0.643 },
@@ -362,6 +356,9 @@ main(int argc, char** argv)
sampler_ctx.mean_radius = 1.0;
sampler_ctx.sigma = 1.18;
+ distrib.material.get_property = get_material_property;
+ distrib.material.material = &sampler_ctx;
+ distrib.caracteristic_length = 1;
distrib.sample = sample_sphere;
distrib.context = &sampler_ctx;
@@ -433,7 +430,7 @@ main(int argc, char** argv)
CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1,
sschiff_uniform_scattering_angles, 1, 1, 1, &estimator), RES_BAD_ARG);
CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1,
- sschiff_uniform_scattering_angles, 2, 1, 1, &estimator), RES_OK);
+ sschiff_uniform_scattering_angles, 3, 1, 1, &estimator), RES_OK);
CHECK(sschiff_estimator_get_wavelengths_count(NULL, NULL), RES_BAD_ARG);
CHECK(sschiff_estimator_get_wavelengths_count(estimator, NULL), RES_BAD_ARG);
@@ -447,15 +444,22 @@ main(int argc, char** argv)
CHECK(sschiff_estimator_get_realisations_count(estimator, &count), RES_OK);
CHECK(count, 1);
- CHECK(sschiff_estimator_get_wavelength_state(NULL, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_state(estimator, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_state(NULL, wlen, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_state(estimator, wlen, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_state(NULL, 0, &state), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_state(estimator, 0, &state), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_state(NULL, wlen, &state), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_state(estimator, wlen, &state), RES_OK);
- CHECK(eq_eps(wlen, state.wavelength, 1.e-12), 1);
+ CHECK(sschiff_estimator_get_wavelength_cross_section
+ (NULL, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_cross_section
+ (estimator, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_cross_section
+ (NULL, wlen, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_cross_section
+ (estimator, wlen, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_cross_section
+ (NULL, 0, &cross_section), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_cross_section
+ (estimator, 0, &cross_section), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_cross_section
+ (NULL, wlen, &cross_section), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelength_cross_section
+ (estimator, wlen, &cross_section), RES_OK);
CHECK(sschiff_estimator_ref_get(NULL), RES_BAD_ARG);
CHECK(sschiff_estimator_ref_get(estimator), RES_OK);
@@ -463,8 +467,8 @@ main(int argc, char** argv)
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
- CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1,
- sschiff_uniform_scattering_angles, 2, 10, 1, &estimator), RES_OK);
+ CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1,
+ sschiff_uniform_scattering_angles, 3, 10, 1, &estimator), RES_OK);
CHECK(sschiff_estimator_get_realisations_count(estimator, &count), RES_OK);
CHECK(count, 10);
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
