commit 7e84fcd1788b967a6af17d953ed911a10febb949
parent ab5834dfd4dc116343fd0deacbd706fd49149f8a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 2 Mar 2016 11:46:16 +0100
Refactor the estimator
Enhance error message and move code.
Diffstat:
4 files changed, 301 insertions(+), 255 deletions(-)
diff --git a/src/sschiff_device.c b/src/sschiff_device.c
@@ -75,6 +75,20 @@ log_error(struct sschiff_device* dev, const char* msg, ...)
}
}
+void
+log_warning(struct sschiff_device* dev, const char* msg, ...)
+{
+ va_list vargs_list;
+ ASSERT(dev && msg);
+ if(dev->verbose) {
+ res_T res; (void)res;
+ va_start(vargs_list, msg);
+ res = logger_vprint(dev->logger, LOG_WARNING, msg, vargs_list);
+ ASSERT(res == RES_OK);
+ va_end(vargs_list);
+ }
+}
+
/*******************************************************************************
* API functions
******************************************************************************/
diff --git a/src/sschiff_device.h b/src/sschiff_device.h
@@ -54,5 +54,16 @@ log_error
__attribute((format(printf, 2, 3)))
#endif
;
+
+extern LOCAL_SYM void
+log_warning
+ (struct sschiff_device* dev,
+ const char* msg,
+ ...)
+#ifdef COMPILER_GCC
+ __attribute((format(printf, 2, 3)))
+#endif
+;
+
#endif /* SSCHIFF_DEVICE_H */
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -371,8 +371,7 @@ accum_differential_cross_section
n_r = ctx->estimator->properties[iwlen].relative_real_refractive_index;
k_r = ctx->estimator->properties[iwlen].relative_imaginary_refractive_index;
- /* Precompute some values. TODO can be stored on the geometry context at
- * its initialisation. This should improve performances*/
+ /* Precompute some values. TODO can be stored on the context */
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];
@@ -407,6 +406,216 @@ accum_differential_cross_section
}
}
+static res_T
+accum_monte_carlo_weight
+ (struct sschiff_device* dev,
+ 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 s3d_hit hit;
+ double sample[2];
+ const float range[2] = { 0, FLT_MAX };
+ size_t nfailures = 0;
+ float axis_x[3], axis_y[3], axis_z[3];
+ float plane_size[2];
+ float ray_org[2][3];
+ float org[3];
+ float x[3], y[3];
+ float plane_area;
+ res_T res = RES_OK;
+ ASSERT(dev && ctx && basis && pos && lower && upper);
+
+ f2_sub(plane_size, upper, lower); /* In micron */
+ plane_area = plane_size[0] * plane_size[1];
+
+ /* Define the projection axis */
+ f3_mulf(axis_x, basis + 0, plane_size[0] * 0.5f);
+ f3_mulf(axis_y, basis + 3, plane_size[1] * 0.5f);
+ f3_set(axis_z, basis + 6);
+
+ /* Compute the origin of the projection plane */
+ f3_add(org, pos, f3_mulf(org, axis_z, lower[2]));
+
+ do {
+ double path_length[2];
+
+ /* Uniformly sample a position onto the projection plane and use it as the
+ * 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[0], f3_add(ray_org[0], x, y), org);
+
+ 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, ctx->scene, &hit, ray_org[0], axis_z, &path_length[0]);
+ if(res != RES_OK) { /* Handle numerical/geometry issues */
+ ++nfailures;
+ continue;
+ }
+ /* Compute and accumulate the cross section weight */
+ accum_cross_section(ctx, plane_area, path_length[0]);
+
+ /* 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(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");
+ }
+ return res;
+}
+
+static res_T
+radiative_properties
+ (struct sschiff_device* dev,
+ const int istep,
+ const size_t ndirs,
+ struct integrator_context* ctx)
+{
+ float lower[3], upper[3];
+ float aabb_pt[8][3];
+ float centroid[3];
+ size_t idir;
+ size_t iwlen;
+ int i;
+ res_T res = RES_OK;
+ ASSERT(dev && ndirs && ctx);
+ (void)istep;
+
+ S3D(scene_get_aabb(ctx->scene, lower, upper));
+
+ /* AABB vertex layout
+ * 6-------7
+ * /' /|
+ * 2-------3 |
+ * | 4.....|.5
+ * |, |/
+ * 0-------1 */
+ f3(aabb_pt[0], lower[0], lower[1], lower[2]);
+ f3(aabb_pt[1], upper[0], lower[1], lower[2]);
+ f3(aabb_pt[2], lower[0], upper[1], lower[2]);
+ f3(aabb_pt[3], upper[0], upper[1], lower[2]);
+ f3(aabb_pt[4], lower[0], lower[1], upper[2]);
+ f3(aabb_pt[5], upper[0], lower[1], upper[2]);
+ f3(aabb_pt[6], lower[0], upper[1], upper[2]);
+ f3(aabb_pt[7], upper[0], upper[1], upper[2]);
+
+ f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f);
+
+ /* 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->estimator->limit_angles[iwlen]);
+ memset(ctx->accums[iwlen].differential_cross_section_cumulative, 0,
+ sizeof(double)*ctx->estimator->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(ctx->rng, dir);
+
+ /* Build the transformation matrix from world to footprint space. Use the
+ * AABB centroid as the origin of the footprint space. */
+ build_basis(dir, basis);
+ build_transform(centroid, basis, transform);
+
+ /* Transform the AABB from world to footprint space */
+ FOR_EACH(i, 0, 8) {
+ f3_add(pt[i], f33_mulf3(pt[i], transform, aabb_pt[i]), transform + 9);
+ }
+
+ /* Compute the AABB of the footprint */
+ f3_splat(lower, FLT_MAX); f3_splat(upper,-FLT_MAX);
+ FOR_EACH(i, 0, 8) {
+ f3_min(lower, lower, pt[i]);
+ f3_max(upper, upper, pt[i]);
+ }
+
+ res = accum_monte_carlo_weight(dev, ctx, basis, centroid, lower, upper);
+ if(res != RES_OK) goto error;
+ }
+ /* 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->estimator->limit_angles[iwlen]) {
+ MC_ACCUM(differential_cross_section[iangle]);
+ MC_ACCUM(differential_cross_section_cumulative[iangle]);
+ }
+
+ #undef MC_ACCUM
+ }
+exit:
+ return res;
+error:
+ 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;
+}
+
struct function_arg {
double scattering_cross_section;
double angle;
@@ -551,16 +760,25 @@ solve_wide_angles_phase_function_model_parameters
/* Retrieve the estimated "n" value */
root = gsl_multiroot_fdfsolver_root(solver);
*n = gsl_vector_get(root, 0);
- printf("theta_limit = %g, n = %g\n", arg.angle, *n); /* TODO remove this */
if(status != GSL_SUCCESS) {
- /* TODO enhance the error message. */
- log_error(dev, "Error in Newton estimation - GSL error: %s - %d.\n",
+ log_error(dev,
+"Cannot estimate the parameters of the wide scattering phase function model.\n"
+"GSL error: %s - %d.\n",
gsl_strerror(status), iter);
res = RES_BAD_OP;
goto error;
}
+ if(*n < 0) {
+ log_error(dev,
+"The estimated parameter `n' of the wide scattering angles phase function\n"
+"model cannot be negative.\n"
+"n = %g\n", *n);
+ res = RES_BAD_OP;
+ goto error;
+ }
+
/* Compute r from the estimated n */
*r = 2 * arg.differential_cross_section * pow(arg.sin_half_angle, *n)
/ (1 + arg.sqr_cos_angle);
@@ -601,8 +819,9 @@ compute_phase_function
double cos_2_limit_angle;
double coef_limit;
double n, r; /* Connector values */
- double scattering_cross_section;
double rcp_scattering_cross_section;
+ double rcp_sqr_scattering_cross_section;
+ struct sschiff_state scattering_cross_section;
struct sschiff_state limit_differential_cross_section;
struct sschiff_state limit_differential_cross_section_cumulative;
size_t ilimit_angle; /* Index of the limit angle of the current wavelength */
@@ -618,9 +837,10 @@ compute_phase_function
cos_2_limit_angle = cos(2.0*limit_angle);
/* Find the connector values between small and whide angles */
- scattering_cross_section = expected_value
+ get_mc_value
(&estimator->accums[iwlen].scattering_cross_section,
- estimator->nrealisations);
+ estimator->nrealisations,
+ &scattering_cross_section);
get_mc_value
(&estimator->accums[iwlen].differential_cross_section[ilimit_angle],
estimator->nrealisations,
@@ -629,12 +849,6 @@ compute_phase_function
(&estimator->accums[iwlen].differential_cross_section_cumulative[ilimit_angle],
estimator->nrealisations,
&limit_differential_cross_section_cumulative);
- printf("sigma_s = %g, Ws = %g +/- %g, Wc = %g +/- %g\n",
- scattering_cross_section,
- limit_differential_cross_section.E,
- limit_differential_cross_section.SE,
- limit_differential_cross_section_cumulative.E,
- limit_differential_cross_section_cumulative.SE);
res = solve_wide_angles_phase_function_model_parameters
(dev,
@@ -642,11 +856,35 @@ compute_phase_function
sin_half_limit_angle,
cos_limit_angle,
cos_2_limit_angle,
- scattering_cross_section,
+ scattering_cross_section.E,
limit_differential_cross_section.E,
limit_differential_cross_section_cumulative.E,
&n, &r);
- if(res != RES_OK) goto error;
+ if(res != RES_OK) {
+ log_error(estimator->dev,
+"Couldn't estimate the parameters of the phase function for wide scattering\n"
+"angles. The phase function is thus not computed.\n"
+"Wavelength = %g microns\n",
+ estimator->wavelengths[iwlen]);
+ goto error;
+ }
+
+ if(n < 2 || n > 4) {
+ log_warning(estimator->dev,
+"The wide scattering angles phase function model parameter `n' is not in\n"
+"[0, 4]. One might increase the number of realisations.\n"
+" n = %g\n"
+" scattering cross section = %g +/- %g\n"
+" differential cross section = %g +/- %g\n"
+" differential cross section cumulative = %g +/- %g\n",
+ n,
+ scattering_cross_section.E,
+ scattering_cross_section.SE,
+ limit_differential_cross_section.E,
+ limit_differential_cross_section.SE,
+ limit_differential_cross_section_cumulative.E,
+ limit_differential_cross_section_cumulative.SE);
+ }
/* Precomputed vlaue */
coef_limit = compute_differential_cross_section_cumulative_term
@@ -683,22 +921,24 @@ compute_phase_function
(estimator->accums[iwlen].differential_cross_section_cumulative[nangles-1].weight,
scattering_cross_section, 1.e-3));*/
- /* Compute the [cumulative] phase function for short angles */
- rcp_scattering_cross_section = 1.0/scattering_cross_section;
- FOR_EACH(iangle, 0, ilimit_angle) {
- get_mc_value
- (&estimator->accums[iwlen].differential_cross_section[iangle],
- estimator->nrealisations,
+ /* Compute the [cumulative] phase function for small angles */
+ rcp_scattering_cross_section = 1.0/scattering_cross_section.E;
+ rcp_sqr_scattering_cross_section =
+ rcp_scattering_cross_section * rcp_scattering_cross_section;
+ FOR_EACH(iangle, 0, ilimit_angle+1) {
+ struct mc_data mc_data;
+
+ mc_data = estimator->accums[iwlen].differential_cross_section[iangle];
+ mc_data.weight *= rcp_scattering_cross_section;
+ mc_data.sqr_weight *= rcp_sqr_scattering_cross_section;
+ get_mc_value(&mc_data, estimator->nrealisations,
&estimator->phase_functions[iwlen].values[iangle]);
- estimator->phase_functions[iwlen].values[iangle].E *=
- rcp_scattering_cross_section;
- get_mc_value
- (&estimator->accums[iwlen].differential_cross_section_cumulative[iangle],
- estimator->nrealisations,
+ mc_data = estimator->accums[iwlen].differential_cross_section_cumulative[iangle];
+ mc_data.weight *= rcp_scattering_cross_section;
+ mc_data.sqr_weight *= rcp_sqr_scattering_cross_section;
+ get_mc_value(&mc_data, estimator->nrealisations,
&estimator->phase_functions[iwlen].cumulative[iangle]);
- estimator->phase_functions[iwlen].cumulative[iangle].E *=
- rcp_scattering_cross_section;
}
/* Compute the [cumulative] phase function for wide angles */
FOR_EACH(iangle, ilimit_angle + 1, nangles) {
@@ -710,239 +950,20 @@ compute_phase_function
estimator->accums[iwlen].differential_cross_section_cumulative[iangle].weight
* rcp_scattering_cross_section;
- /* THe phase function for wide angles is analitically computed, i.e. there
+ /* The phase function for wide angles is analitically computed, i.e. there
* is no variance or standard error */
- estimator->phase_functions[iwlen].values[iangle].V = -1;
- estimator->phase_functions[iwlen].values[iangle].SE = -1;
- estimator->phase_functions[iwlen].cumulative[iangle].V = -1;
- estimator->phase_functions[iwlen].cumulative[iangle].SE = -1;
+ estimator->phase_functions[iwlen].values[iangle].V = 0;
+ estimator->phase_functions[iwlen].values[iangle].SE = 0;
+ estimator->phase_functions[iwlen].cumulative[iangle].V = 0;
+ estimator->phase_functions[iwlen].cumulative[iangle].SE = 0;
}
- /* Setup the [cumulative] phase function for the limit angle */
- estimator->phase_functions[iwlen].values[ilimit_angle] =
- limit_differential_cross_section;
- estimator->phase_functions[iwlen].values[ilimit_angle].E *=
- rcp_scattering_cross_section;
- estimator->phase_functions[iwlen].cumulative[ilimit_angle] =
- limit_differential_cross_section_cumulative;
- estimator->phase_functions[iwlen].cumulative[ilimit_angle].E *=
- rcp_scattering_cross_section;
-
exit:
return res;
error:
goto exit;
}
-static res_T
-accum_monte_carlo_weight
- (struct sschiff_device* dev,
- 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 s3d_hit hit;
- double sample[2];
- const float range[2] = { 0, FLT_MAX };
- size_t nfailures = 0;
- float axis_x[3], axis_y[3], axis_z[3];
- float plane_size[2];
- float ray_org[2][3];
- float org[3];
- float x[3], y[3];
- float plane_area;
- res_T res = RES_OK;
- ASSERT(dev && ctx && basis && pos && lower && upper);
-
- f2_sub(plane_size, upper, lower); /* In micron */
- plane_area = plane_size[0] * plane_size[1];
-
- /* Define the projection axis */
- f3_mulf(axis_x, basis + 0, plane_size[0] * 0.5f);
- f3_mulf(axis_y, basis + 3, plane_size[1] * 0.5f);
- f3_set(axis_z, basis + 6);
-
- /* Compute the origin of the projection plane */
- f3_add(org, pos, f3_mulf(org, axis_z, lower[2]));
-
- do {
- double path_length[2];
-
- /* Uniformly sample a position onto the projection plane and use it as the
- * 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[0], f3_add(ray_org[0], x, y), org);
-
- 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, ctx->scene, &hit, ray_org[0], axis_z, &path_length[0]);
- if(res != RES_OK) { /* Handle numerical/geometry issues */
- ++nfailures;
- continue;
- }
- /* Compute and accumulate the cross section weight */
- accum_cross_section(ctx, plane_area, path_length[0]);
-
- /* 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(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");
- }
- return res;
-}
-
-static res_T
-radiative_properties
- (struct sschiff_device* dev,
- const int istep,
- const size_t ndirs,
- struct integrator_context* ctx)
-{
- float lower[3], upper[3];
- float aabb_pt[8][3];
- float centroid[3];
- size_t idir;
- size_t iwlen;
- int i;
- res_T res = RES_OK;
- ASSERT(dev && ndirs && ctx);
- (void)istep;
-
- S3D(scene_get_aabb(ctx->scene, lower, upper));
-
- /* AABB vertex layout
- * 6-------7
- * /' /|
- * 2-------3 |
- * | 4.....|.5
- * |, |/
- * 0-------1 */
- f3(aabb_pt[0], lower[0], lower[1], lower[2]);
- f3(aabb_pt[1], upper[0], lower[1], lower[2]);
- f3(aabb_pt[2], lower[0], upper[1], lower[2]);
- f3(aabb_pt[3], upper[0], upper[1], lower[2]);
- f3(aabb_pt[4], lower[0], lower[1], upper[2]);
- f3(aabb_pt[5], upper[0], lower[1], upper[2]);
- f3(aabb_pt[6], lower[0], upper[1], upper[2]);
- f3(aabb_pt[7], upper[0], upper[1], upper[2]);
-
- f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f);
-
- /* 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->estimator->limit_angles[iwlen]);
- memset(ctx->accums[iwlen].differential_cross_section_cumulative, 0,
- sizeof(double)*ctx->estimator->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(ctx->rng, dir);
-
- /* Build the transformation matrix from world to footprint space. Use the
- * AABB centroid as the origin of the footprint space. */
- build_basis(dir, basis);
- build_transform(centroid, basis, transform);
-
- /* Transform the AABB from world to footprint space */
- FOR_EACH(i, 0, 8) {
- f3_add(pt[i], f33_mulf3(pt[i], transform, aabb_pt[i]), transform + 9);
- }
-
- /* Compute the AABB of the footprint */
- f3_splat(lower, FLT_MAX); f3_splat(upper,-FLT_MAX);
- FOR_EACH(i, 0, 8) {
- f3_min(lower, lower, pt[i]);
- f3_max(upper, upper, pt[i]);
- }
-
- res = accum_monte_carlo_weight(dev, ctx, basis, centroid, lower, upper);
- if(res != RES_OK) goto error;
- }
- /* 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->estimator->limit_angles[iwlen]) {
- MC_ACCUM(differential_cross_section[iangle]);
- MC_ACCUM(differential_cross_section_cumulative[iangle]);
- }
-
- #undef MC_ACCUM
- }
-exit:
- return res;
-error:
- 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;
-}
static res_T
inverse_cumulative_phase_function
@@ -1440,7 +1461,7 @@ sschiff_integrate
}
}
-#if 0
+#if 1
/* TODO use parallelism */
FOR_EACH(i, 0, nwavelengths) {
res = compute_phase_function(dev, i, nangles, estimator);
diff --git a/src/test_sschiff_estimator_rhodo.c b/src/test_sschiff_estimator_rhodo.c
@@ -89,7 +89,7 @@ main(int argc, char** argv)
const double wavelength = 0.4; /* In micron */
const size_t nscatt_angles = 1000;
const size_t ngeoms = 10000;
- const size_t ndirs = 100;
+ const size_t ndirs = 10;
(void)argc, (void)argv;
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
