commit ab5834dfd4dc116343fd0deacbd706fd49149f8a
parent 1da84d8ec665aeff9b43f4e20f6d4d2d0708fa2d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 2 Mar 2016 10:32:35 +0100
Refactor the schiff estimator
Remove the geometry context. Move its data to the estimator.
Diffstat:
2 files changed, 160 insertions(+), 330 deletions(-)
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -81,55 +81,47 @@ struct mc_accum {
/* Estimated phase function */
struct phase_function {
- size_t ilimit_angle; /* Index of the firs wide angle */
-
/* Per angle values */
struct sschiff_state* values;
struct sschiff_state* 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* wavenumbers;
- struct sschiff_material_properties* properties; /* Material properties */
-};
-
-static struct geometry_context GEOMETRY_CONTEXT_NULL = {NULL, NULL, NULL, NULL};
-
/* 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 */
-
+ struct ssp_rng* rng; /* Random Number Generator */
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 */
+
+ /* Pointer toward the estimator */
+ struct sschiff_estimator* estimator;
};
struct sschiff_estimator {
double* wavelengths; /* List of wavelengths in vacuum */
double* angles; /* List of scattering angles */
- struct phase_function* phase_functions; /* Per wavelength phase function */
- struct mc_accum* accums; /* Per wavelength Monte Carlo estimation */
+
+ /* 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* wavenumbers;
+ struct sschiff_material_properties* properties; /* Material properties */
+
+ /* Per wavelength results */
+ struct phase_function* phase_functions;
+ struct mc_accum* accums; /* Monte Carlo estimation */
+
+ size_t nrealisations; /* # realisation used by MC estimations */
+
struct sschiff_device* dev;
- size_t nrealisations;
ref_T ref;
};
@@ -232,80 +224,6 @@ get_mc_value
value->SE = sqrt(value->V / (double)nrealisations);
}
-/* Dump a thumbnail of the sampled geometry in given sampled direction */
-static INLINE void
-draw_thumbnail
- (struct s3d_scene* scn,
- 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 char* name) /* Filename of the thumbnail */
-{
-#define DEFINITION 128
- STATIC_ASSERT(IS_POW2(DEFINITION), Invalid_thumbnail_definition);
- const float range[2] = { 0, FLT_MAX };
- const float pixel_size = 1.f/(float)DEFINITION;
- float axis_x[3], axis_y[3], axis_z[3];
- float plane_size[2];
- float ray_org[3];
- float org[3];
- uint32_t npixels = (uint32_t)(DEFINITION * DEFINITION);
- uint32_t mcode;
- unsigned char* img = NULL;
- ASSERT(basis && pos && lower && upper);
-
- /* Allocate the thumbnail buffer */
- img = sa_add(img, DEFINITION*DEFINITION);
- ASSERT(img != NULL);
-
- /* Define the projection axis */
- f2_sub(plane_size, upper, lower);
- 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]));
-
- FOR_EACH(mcode, 0, npixels) {
- struct s3d_hit hit;
- size_t ipixel;
- float sample[2], x[3], y[3];
- uint16_t ipixel_x, ipixel_y;
-
- /* Compute the sample position in [0, 1)^2 onto the projection plane */
- ipixel_x = morton2D_decode(mcode);
- ipixel_y = morton2D_decode(mcode>>1);
- sample[0] = ((float)ipixel_x + 0.5f) * pixel_size;
- sample[1] = ((float)ipixel_y + 0.5f) * pixel_size;
-
- /* Compute the ray origin */
- f3_mulf(x, axis_x, sample[0]*2.f - 1.f);
- f3_mulf(y, axis_y, sample[1]*2.f - 1.f);
- f3_add(ray_org, f3_add(ray_org, x, y), org);
-
- S3D(scene_trace_ray(scn, ray_org, axis_z, range, &hit));
-
- /* Simple shading */
- ipixel = (size_t)(ipixel_y * DEFINITION + ipixel_x);
- if(S3D_HIT_NONE(&hit)) {
- img[ipixel] = 0;
- } else {
- float N[3] = {0.f,0.f,0.f};
- float cosine;
- f3_normalize(N, hit.normal);
- if(0 > (cosine = f3_dot(N, axis_z))) {
- cosine = f3_dot(f3_minus(N, N), axis_z);
- }
- img[ipixel] = (unsigned char)(cosine*255.f + 0.5f/*Round*/);
- }
- }
- image_ppm_write(name, DEFINITION, DEFINITION, 1, img);
- sa_release(img);
-#undef DEFINITION
-}
-
static FINLINE int
hit_on_edge(const struct s3d_hit* hit)
{
@@ -403,9 +321,9 @@ accum_cross_section
double k_e;
double w_extinction, w_absorption, w_scattering;
- k_e = ctx->geom_ctx->wavenumbers[iwlen];
- n_r = ctx->geom_ctx->properties[iwlen].relative_real_refractive_index;
- k_r = ctx->geom_ctx->properties[iwlen].relative_imaginary_refractive_index;
+ k_e = ctx->estimator->wavenumbers[iwlen];
+ n_r = ctx->estimator->properties[iwlen].relative_real_refractive_index;
+ k_r = ctx->estimator->properties[iwlen].relative_imaginary_refractive_index;
/* Extinction cross section */
w_extinction = 2.0 * plane_area *
@@ -449,9 +367,9 @@ accum_differential_cross_section
size_t iangle;
/* Fetch optical properties */
- k_e = ctx->geom_ctx->wavenumbers[iwlen];
- n_r = ctx->geom_ctx->properties[iwlen].relative_real_refractive_index;
- k_r = ctx->geom_ctx->properties[iwlen].relative_imaginary_refractive_index;
+ k_e = ctx->estimator->wavenumbers[iwlen];
+ 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*/
@@ -471,19 +389,19 @@ accum_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]) {
+ FOR_EACH(iangle, 0, ctx->estimator->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;
+ k_e_angle_delta_r = ctx->estimator->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;
+ weight = 2*PI * ctx->estimator->angles[iangle] * bessel_j1 / k_e_delta_r * tmp;
ctx->accums[iwlen].differential_cross_section_cumulative[iangle] += weight;
}
}
@@ -524,13 +442,13 @@ function_fdf(const gsl_vector* vec, void* params, gsl_vector* res, gsl_matrix* J
+ 3*n*n
- 8*(n - 2) * arg->cos_angle
- 26*n + 72 )
- - pow(arg->sin_half_angle, n) * (4*n*n - 24*n + 40) );
+ - pow(arg->sin_half_angle, n) * (4*n*n - 24*n + 64) );
v = ((n - 2) * (n - 4) * (n - 6));
u_prime = arg->sqr_sin_half_angle
* ((2*n - 6) * arg->cos_2_angle + 6*n - 8*arg->cos_angle - 26)
- pow(arg->sin_half_angle, n)
- * (log(arg->sin_half_angle) * (4*n*n - 24*n + 40) + 8*n-24);
+ * (log(arg->sin_half_angle) * (4*n*n - 24*n + 64) + 8*n-24);
v_prime = (n-4)*(n-6) + (n-2)*(n-6) + (n-2)*(n-4);
f = arg->differential_cross_section_cumulative
@@ -675,7 +593,6 @@ compute_phase_function
(struct sschiff_device* dev,
const size_t iwlen,
const size_t nangles,
- const struct geometry_context* ctx,
struct sschiff_estimator* estimator)
{
double limit_angle;
@@ -691,10 +608,10 @@ compute_phase_function
size_t ilimit_angle; /* Index of the limit angle of the current wavelength */
size_t iangle;
res_T res = RES_OK;
- ASSERT(dev && ctx && estimator);
+ ASSERT(dev && estimator);
/* Fetch the limit angle and precompute some values */
- ilimit_angle = ctx->limit_angles[iwlen]-1;
+ ilimit_angle = estimator->limit_angles[iwlen]-1;
limit_angle = estimator->angles[ilimit_angle];
sin_half_limit_angle = sin(0.5*limit_angle);
cos_limit_angle = cos(limit_angle);
@@ -955,9 +872,9 @@ radiative_properties
/* 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]);
+ sizeof(double)*ctx->estimator->limit_angles[iwlen]);
memset(ctx->accums[iwlen].differential_cross_section_cumulative, 0,
- sizeof(double)*ctx->geom_ctx->limit_angles[iwlen]);
+ sizeof(double)*ctx->estimator->limit_angles[iwlen]);
}
FOR_EACH(idir, 0, ndirs) {
@@ -987,15 +904,6 @@ radiative_properties
res = accum_monte_carlo_weight(dev, ctx, basis, centroid, lower, upper);
if(res != RES_OK) goto error;
-
-#if 0
- {
- /* Compute an image of the shape transformed in the footprint space */
- char thumb_name[64];
- sprintf(thumb_name, "%d_%lu.ppm", istep, idir);
- draw_thumbnail(scene, basis, centroid, lower, upper, thumb_name);
- }
-#endif
}
/* Compute the Monte Carlo weight of the temporary accumulator */
FOR_EACH(iwlen, 0, ctx->nwavelengths) {
@@ -1013,7 +921,7 @@ radiative_properties
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]) {
+ FOR_EACH(iangle, 0, ctx->estimator->limit_angles[iwlen]) {
MC_ACCUM(differential_cross_section[iangle]);
MC_ACCUM(differential_cross_section_cumulative[iangle]);
}
@@ -1050,7 +958,7 @@ inverse_cumulative_phase_function
ASSERT(cumulative >= 0.0 && cumulative <= 1.0);
ASSERT(iwlen < sa_size(estimator->phase_functions));
- nsmall_angles = estimator->phase_functions[iwlen].ilimit_angle;
+ nsmall_angles = estimator->limit_angles[iwlen];
if(cumulative < cumulative_small_angles[nsmall_angles-1]) {
/* Look for the the cumulative in the filtered small angles array */
double* found = search_lower_bound(&cumulative, cumulative_small_angles,
@@ -1102,138 +1010,32 @@ check_distribution(struct sschiff_geometry_distribution* distrib)
}
static void
-geometry_context_release(struct geometry_context* ctx)
-{
- ASSERT(ctx);
- if(ctx->limit_angles) sa_release(ctx->limit_angles);
- if(ctx->wavenumbers) sa_release(ctx->wavenumbers);
- if(ctx->properties) sa_release(ctx->properties);
-}
-
-/* Return RES_BAD_OP if the limit angle of the phase function is invalid */
-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;
- int skip_differential_cross_section = 0;
- 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->wavenumbers, nwavelengths)) {
- log_error(dev, "Couldn't allocate the list of wavenumbers.\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->wavenumbers[iwlen] = 2.0*PI / lambda_e;
-
- if(skip_differential_cross_section)
- continue;
-
- /* Search the limit angle */
- theta_l = sqrt(lambda_e * rcp_PI_Lc);
- if(theta_l <= 0 || theta_l >= PI) {
- log_error(dev,
-"Invalid theta limit, i.e. angle between small and wide scattering angles\n"
-"`%g'. The phase function for the wavelengths `%g' and its [inverse]\n"
-"cumulative will be not computed.\n",
- theta_l, wavelengths[iwlen]);
- res = RES_BAD_ARG;
- goto error;
- }
- angle = search_lower_bound
- (&theta_l, angles, nangles, sizeof(double), cmp_double);
- if(!angle) {
- log_error(dev,
-"Can't find a limit scattering angle for theta limit `%g'.\n", theta_l);
- res = RES_BAD_ARG;
- goto error;
- }
- if(eq_eps(*angle, PI, 1.e-6)) {
- log_error(dev, "Invalid limit scattering 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);
- }
-
- if(skip_differential_cross_section) {
-
- }
-
- ctx->angles = angles;
-
-exit:
- return res;
-error:
- geometry_context_release(ctx);
- *ctx = GEOMETRY_CONTEXT_NULL;
- 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->estimator) SSCHIFF(estimator_ref_put(ctx->estimator));
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)
+ #define RELEASE(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);
+ RELEASE(ctx->accums[iwlen].differential_cross_section);
+ RELEASE(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);
+ RELEASE(ctx->mc_accums[iwlen].differential_cross_section);
+ RELEASE(ctx->mc_accums[iwlen].differential_cross_section_cumulative);
}
sa_release(ctx->mc_accums);
}
- #undef CTXFREE
+ #undef RELEASE
}
static res_T
@@ -1241,13 +1043,13 @@ integrator_context_init
(struct integrator_context* ctx,
struct s3d_device* s3d,
struct ssp_rng* rng,
- const struct geometry_context* geom_ctx,
+ struct sschiff_estimator* estimator,
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);
+ ASSERT(ctx && s3d && rng && nwavelengths && nangles >= 3);
memset(ctx, 0, sizeof(struct integrator_context));
@@ -1255,28 +1057,29 @@ integrator_context_init
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) { \
+ #define RESIZE(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);
+ RESIZE(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);
+ RESIZE(ctx->accums[iwlen].differential_cross_section, nangles);
+ RESIZE(ctx->accums[iwlen].differential_cross_section_cumulative, nangles);
}
- CTXALLOC(ctx->mc_accums, nwavelengths);
+ RESIZE(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);
+ RESIZE(ctx->mc_accums[iwlen].differential_cross_section, nangles);
+ RESIZE(ctx->mc_accums[iwlen].differential_cross_section_cumulative, nangles);
}
- #undef CTXALLOC
+ #undef RESIZE
SSP(rng_ref_get(rng));
ctx->rng = rng;
- ctx->geom_ctx = geom_ctx;
+ SSCHIFF(estimator_ref_get(estimator));
+ ctx->estimator = estimator;
ctx->nwavelengths = nwavelengths;
ctx->nangles = nangles;
@@ -1331,33 +1134,36 @@ estimator_release(ref_T* ref)
{
struct sschiff_estimator* estimator;
struct sschiff_device* dev;
+ size_t nwavelengths;
+ size_t i;
ASSERT(ref);
+
estimator = CONTAINER_OF(ref, struct sschiff_estimator, ref);
dev = estimator->dev;
- if(estimator->wavelengths) sa_release(estimator->wavelengths);
- if(estimator->angles) sa_release(estimator->angles);
+ nwavelengths = sa_size(estimator->accums);
+
+ #define RELEASE(Data) if(Data) sa_release(Data)
+ RELEASE(estimator->wavelengths);
+ RELEASE(estimator->angles);
+ RELEASE(estimator->limit_angles);
+ RELEASE(estimator->wavenumbers);
+ RELEASE(estimator->properties);
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);
+ RELEASE(estimator->accums[i].differential_cross_section);
+ RELEASE(estimator->accums[i].differential_cross_section_cumulative);
}
sa_release(estimator->accums);
}
if(estimator->phase_functions) {
- const size_t nwavelengths = sa_size(estimator->phase_functions);
- size_t i;
FOR_EACH(i, 0, nwavelengths) {
- if(estimator->phase_functions[i].values)
- sa_release(estimator->phase_functions[i].values);
- if(estimator->phase_functions[i].cumulative)
- sa_release(estimator->phase_functions[i].cumulative);
+ RELEASE(estimator->phase_functions[i].values);
+ RELEASE(estimator->phase_functions[i].cumulative);
}
sa_release(estimator->phase_functions);
}
+ #undef RELEASE
+
MEM_RM(dev->allocator, estimator);
SSCHIFF(device_ref_put(dev));
}
@@ -1365,6 +1171,7 @@ estimator_release(ref_T* ref)
static res_T
estimator_create
(struct sschiff_device* dev,
+ struct sschiff_geometry_distribution* distrib,
const double* wavelengths,
const size_t nwavelengths,
sschiff_scattering_angles_distribution_T angles_distrib,
@@ -1372,10 +1179,11 @@ estimator_create
struct sschiff_estimator** out_estimator)
{
struct sschiff_estimator* estimator = NULL;
+ double rcp_PI_Lc;
size_t i;
res_T res = RES_OK;
- ASSERT(dev && wavelengths && nwavelengths && angles_distrib && nangles);
- ASSERT(out_estimator);
+ ASSERT(dev && distrib && wavelengths && nwavelengths);
+ ASSERT(angles_distrib && nangles && out_estimator);
estimator = MEM_CALLOC(dev->allocator, 1, sizeof(struct sschiff_estimator));
if(!estimator) {
@@ -1388,23 +1196,24 @@ estimator_create
SSCHIFF(device_ref_get(dev));
estimator->dev = dev;
+ #define RESIZE(Array, Count, ErrMsg) { \
+ if(!sa_add(Array, Count)) { \
+ log_error(dev, ErrMsg); \
+ res = RES_MEM_ERR; \
+ goto error; \
+ } \
+ memset(Array, 0, sizeof(Array[0])*Count); \
+ } (void)0
+
/* Setup and sort the wavelengths */
- if(!sa_add(estimator->wavelengths, nwavelengths)) {
- log_error(dev,
- "Couldn't allocate the list of estimated wavelengths.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- FOR_EACH(i, 0, nwavelengths) estimator->wavelengths[i] = wavelengths[i];
+ RESIZE(estimator->wavelengths, nwavelengths,
+ "Couldn't allocate the list of estimated wavelengths.\n");
+ memcpy(estimator->wavelengths, wavelengths, nwavelengths*sizeof(double));
qsort(estimator->wavelengths, nwavelengths, sizeof(double), cmp_double);
/* Generate the scattering angles */
- if(!sa_add(estimator->angles, nangles)) {
- log_error(dev,
- "Couldn't allocate the list of scattering angles.\n");
- res = RES_MEM_ERR;
- goto error;
- }
+ RESIZE(estimator->angles, nangles,
+ "Couldn't allocate the list of scattering angles.\n");
angles_distrib(estimator->angles, nangles);
if(estimator->angles[0] != 0.f
|| !eq_eps(estimator->angles[nangles-1], PI, 1.e-6)) {
@@ -1422,52 +1231,81 @@ estimator_create
}
}
- /* 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);
+ /* Allocate miscellaneous array of temporary/precomputed data */
+ RESIZE(estimator->limit_angles, nwavelengths,
+ "Couldn't allocate the per wavelength indices of the limit scattering angle.\n");
+ RESIZE(estimator->wavenumbers, nwavelengths,
+ "Couldn't allocate the list of wavenumbers.\n");
+ RESIZE(estimator->properties, nwavelengths,
+ "Couldn't allocate the per wavelength optical properties.\n");
+
+ rcp_PI_Lc = 1.0 / (PI*distrib->caracteristic_length);
FOR_EACH(i, 0, nwavelengths) {
- if(!sa_add(estimator->accums[i].differential_cross_section, nangles)) {
+ double lambda_e, theta_l;
+ double* angle;
+
+ /* Fetch the particle optical properties */
+ distrib->material.get_property
+ (distrib->material.material,
+ estimator->wavelengths[i],
+ &estimator->properties[i]);
+
+ /* Precompute the wavenumbers */
+ lambda_e =
+ estimator->wavelengths[i]
+ / estimator->properties[i].medium_refractive_index;
+ estimator->wavenumbers[i] = 2.0*PI/lambda_e;
+
+ /* Search for limit scattering angle */
+ theta_l = sqrt(lambda_e * rcp_PI_Lc);
+ if(theta_l <= 0 || theta_l >= PI) {
log_error(dev,
- "Couldn't allocate the differential cross sections to estimate.\n");
- res = RES_MEM_ERR;
+"Invalid theta limit, i.e. angle between small and wide scattering angles\n"
+"`%g'. The phase function for the wavelengths `%g' and its [inverse]\n"
+"cumulative will be not computed.\n",
+ theta_l, estimator->wavelengths[i]);
+ res = RES_BAD_ARG;
goto error;
}
- if(!sa_add(estimator->accums[i].differential_cross_section_cumulative, nangles)) {
+ angle = search_lower_bound(&theta_l, estimator->angles, nangles,
+ sizeof(double), cmp_double);
+ if(!angle) {
log_error(dev,
-"Couldn't allocate the cumulative differential cross sections to estimate.\n");
- res = RES_MEM_ERR;
+"Can't find a limit scattering angle for theta limit `%g'.\n", theta_l);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(eq_eps(*angle, PI, 1.e-6)) {
+ log_error(dev, "Invalid limit scattering angle `%g'.\n", *angle);
+ res = RES_BAD_ARG;
goto error;
}
+ /* Define the index of the first "wide scattering angle" */
+ estimator->limit_angles[i] = (size_t)(angle - estimator->angles);
+ ASSERT(estimator->limit_angles[i] < nangles);
}
- /* Allocate the phase function data */
- if(!sa_add(estimator->phase_functions, nwavelengths)) {
- log_error(dev,
- "Couldn't allocate the per wavelength phase functions data.\n");
- res = RES_MEM_ERR;
- goto error;
+ /* Allocate the estimator accumulators */
+ RESIZE(estimator->accums, nwavelengths,
+ "Couldn't allocate the Monte Carlo accumulator of the estimator.\n");
+ memset(estimator->accums, 0, sizeof(estimator->accums[0])*nwavelengths);
+ FOR_EACH(i, 0, nwavelengths) {
+ RESIZE(estimator->accums[i].differential_cross_section, nangles,
+ "Couldn't allocate the differential cross sections to estimate.\n");
+ RESIZE(estimator->accums[i].differential_cross_section_cumulative, nangles,
+ "Couldn't allocate the cumulative differential cross sections to estimate.\n");
}
- memset(estimator->phase_functions, 0,
- sizeof(estimator->phase_functions[0])*nwavelengths);
+
+ /* Allocate the phase function data */
+ RESIZE(estimator->phase_functions, nwavelengths,
+ "Couldn't allocate the per wavelength phase functions data.\n");
FOR_EACH(i, 0, nwavelengths) {
- if(!sa_add(estimator->phase_functions[i].values, nangles)) {
- log_error(dev,
- "Couldn't allocate the per angle phase function values.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- if(!sa_add(estimator->phase_functions[i].cumulative, nangles)) {
- log_error(dev,
- "Couldn't allocate the per angle cumulative phase function.\n");
- res = RES_MEM_ERR;
- goto error;
- }
+ RESIZE(estimator->phase_functions[i].values, nangles,
+ "Couldn't allocate the per angle phase function values.\n");
+ RESIZE(estimator->phase_functions[i].cumulative, nangles,
+ "Couldn't allocate the per angle cumulative phase function.\n");
}
+ #undef RESIZE
exit:
*out_estimator = estimator;
@@ -1496,7 +1334,6 @@ sschiff_integrate
const size_t ndirs,
struct sschiff_estimator** out_estimator)
{
- struct geometry_context geom_ctx = GEOMETRY_CONTEXT_NULL;
struct integrator_context* ctxs = NULL;
struct sschiff_estimator* estimator = NULL;
struct ssp_rng** rngs = NULL;
@@ -1505,7 +1342,6 @@ sschiff_integrate
size_t iwlen;
int igeom;
ATOMIC res = (ATOMIC)RES_OK;
- memset(&geom_ctx, 0, sizeof(geom_ctx));
if(!dev || !rng || !distrib || !wavelengths || !nwavelengths
|| !angles_distrib || nangles < 3 || !ngeoms || !ndirs || !out_estimator
@@ -1514,16 +1350,11 @@ sschiff_integrate
}
/* Create the Schiff estimator */
- res = estimator_create
- (dev, wavelengths, nwavelengths, angles_distrib, nangles, &estimator);
+ res = estimator_create(dev, distrib, 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);
@@ -1554,7 +1385,7 @@ sschiff_integrate
"Couldn't create the RNG of the thread \"%lu\".\n", (unsigned long)i);
goto error;
}
- res = integrator_context_init(ctxs + i, dev->s3d[i], rngs[i], &geom_ctx,
+ res = integrator_context_init(ctxs + i, dev->s3d[i], rngs[i], estimator,
nwavelengths, nangles);
if(res != RES_OK) {
log_error(dev,
@@ -1602,23 +1433,22 @@ sschiff_integrate
MC_ACCUM(avg_projected_area);
/* Accum up to "limit angle"; Great angles are analitically computed */
- FOR_EACH(iangle, 0, geom_ctx.limit_angles[iwlen]) {
+ FOR_EACH(iangle, 0, estimator->limit_angles[iwlen]) {
MC_ACCUM(differential_cross_section[iangle]);
MC_ACCUM(differential_cross_section_cumulative[iangle]);
}
}
}
-#if 1
+#if 0
/* TODO use parallelism */
FOR_EACH(i, 0, nwavelengths) {
- res = compute_phase_function(dev, i, nangles, &geom_ctx, estimator);
+ res = compute_phase_function(dev, i, nangles, estimator);
if(res != RES_OK) goto error;
}
#endif
exit:
- geometry_context_release(&geom_ctx);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
if(ctxs) {
FOR_EACH(i, 0, dev->nthreads) integrator_context_release(ctxs+i);
@@ -1751,7 +1581,7 @@ sschiff_estimator_inverse_cumulative_phase_function
/* Allocate the "filtered" phase function cumulative of small angles, i.e.
* the increasing cumulative defined from the Monte-Carlo estimation. */
- nsmall_angles = estimator->phase_functions[iwlen].ilimit_angle;
+ nsmall_angles = estimator->limit_angles[iwlen];
if(!sa_add(cumulative_small_angles, nsmall_angles)) {
log_error(estimator->dev,
"Couldn't allocate the filtered phase function cumulative of small angles.\n");
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -106,7 +106,7 @@ check_schiff_estimation
size_t i;
const size_t nscatt_angles = 1000;
- const size_t ngeoms = 1000;
+ const size_t ngeoms = 100;
const size_t ndirs = 100;
const double wavelength = sampler_ctx->wavelength;
