solstice-solver

Solver library of the solstice app
git clone git://git.meso-star.com/solstice-solver.git
Log | Files | Refs | README | LICENSE

commit 0ee3b390d45d05c990c2d53406df787cd9a1f933
parent 87a16b4cdda7ea6e0f3c77043241e458b24ca843
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 16 Feb 2017 16:35:10 +0100

Merge remote-tracking branch 'origin/primaries' into rcv_primitive

Diffstat:
Msrc/ssol.h | 30++++++++++++++++++++----------
Msrc/ssol_estimator.c | 88++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------
Msrc/ssol_estimator_c.h | 154++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---
Msrc/ssol_instance.c | 12++++++++++++
Msrc/ssol_instance_c.h | 1+
Msrc/ssol_object.c | 15+++++++++++++++
Msrc/ssol_object_c.h | 1+
Msrc/ssol_scene.c | 25+++++++------------------
Msrc/ssol_scene_c.h | 6+++---
Msrc/ssol_shape.c | 69+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------
Msrc/ssol_shape_c.h | 1+
Msrc/ssol_solver.c | 205++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------
Msrc/test_ssol_by_receiver_integration.c | 8+++++---
13 files changed, 504 insertions(+), 111 deletions(-)

diff --git a/src/ssol.h b/src/ssol.h @@ -572,6 +572,12 @@ SSOL_API res_T ssol_object_clear (struct ssol_object* object); +/* Retrieve the area of the object */ +SSOL_API res_T +ssol_object_get_area + (const struct ssol_object* object, + double* area); + /******************************************************************************* * Object Instance API - Clone of an object with a set of per instance data as * world transformation, material parameters, etc. Note that the object @@ -616,6 +622,12 @@ ssol_instance_get_id (const struct ssol_instance* instance, uint32_t* id); +/* Retrieve the area of the instance */ +SSOL_API res_T +ssol_instance_get_area + (const struct ssol_instance* instance, + double* area); + /******************************************************************************* * Param buffer API ******************************************************************************/ @@ -780,6 +792,14 @@ ssol_estimator_get_mc_global struct ssol_mc_global* mc_global); SSOL_API res_T +ssol_estimator_get_primary_entity_x_receiver_status + (struct ssol_estimator* estimator, + const struct ssol_instance* prim_instance, + const struct ssol_instance* recv_instance, + const enum ssol_side_flag side, + struct ssol_mc_receiver* rcv); + +SSOL_API res_T ssol_estimator_get_count (const struct ssol_estimator* estimator, size_t* count); @@ -789,16 +809,6 @@ ssol_estimator_get_failed_count (const struct ssol_estimator* estimator, size_t* count); -SSOL_API res_T -ssol_estimator_get_sampled_area - (const struct ssol_estimator* estimator, - double* area); - -SSOL_API res_T -ssol_estimator_get_primary_area - (const struct ssol_estimator* estimator, - double* area); - /******************************************************************************* * Per receiver MC estimations ******************************************************************************/ diff --git a/src/ssol_estimator.c b/src/ssol_estimator.c @@ -36,6 +36,7 @@ create_mc_receivers { struct htable_instance_iterator it, end; struct mc_receiver mc_rcv_null; + struct mc_per_primary_data mc_samp_null; res_T res = RES_OK; ASSERT(scene && estimator); @@ -43,20 +44,28 @@ create_mc_receivers htable_instance_end(&scene->instances_rt, &end); mc_receiver_init(estimator->dev->allocator, &mc_rcv_null); + init_mc_per_prim_data(estimator->dev->allocator, &mc_samp_null); while(!htable_instance_iterator_eq(&it, &end)) { const struct ssol_instance* inst = *htable_instance_iterator_data_get(&it); htable_instance_iterator_next(&it); - if(!inst->receiver_mask) continue; - - res = htable_receiver_set(&estimator->mc_receivers, &inst, &mc_rcv_null); - if(res != RES_OK) goto error; + if(inst->receiver_mask) { + res = htable_receiver_set(&estimator->mc_receivers, &inst, &mc_rcv_null); + if(res != RES_OK) goto error; + } + if(inst->sample) { + res = htable_primary_set(&estimator->global_primaries, &inst, &mc_samp_null); + if(res != RES_OK) goto error; + } } exit: + mc_receiver_release(&mc_rcv_null); + release_mc_per_prim_data(&mc_samp_null); return res; error: htable_receiver_clear(&estimator->mc_receivers); + htable_primary_clear(&estimator->global_primaries); goto exit; } @@ -69,6 +78,7 @@ estimator_release(ref_T* ref) ASSERT(ref); dev = estimator->dev; htable_receiver_release(&estimator->mc_receivers); + htable_primary_release(&estimator->global_primaries); ASSERT(dev && dev->allocator); MEM_RM(dev->allocator, estimator); SSOL(device_ref_put(dev)); @@ -79,7 +89,7 @@ estimator_release(ref_T* ref) ******************************************************************************/ res_T ssol_estimator_ref_get -(struct ssol_estimator* estimator) + (struct ssol_estimator* estimator) { if (!estimator) return RES_BAD_ARG; ref_get(&estimator->ref); @@ -87,10 +97,9 @@ ssol_estimator_ref_get } res_T -ssol_estimator_ref_put -(struct ssol_estimator* estimator) +ssol_estimator_ref_put(struct ssol_estimator* estimator) { - if (!estimator) return RES_BAD_ARG; + if(!estimator) return RES_BAD_ARG; ref_put(&estimator->ref, estimator_release); return RES_OK; } @@ -116,6 +125,48 @@ ssol_estimator_get_mc_global } res_T +ssol_estimator_get_primary_entity_x_receiver_status + (struct ssol_estimator* estimator, + const struct ssol_instance* prim_instance, + const struct ssol_instance* recv_instance, + const enum ssol_side_flag side, + struct ssol_mc_receiver* rcv) +{ + struct mc_per_primary_data* prim_data = NULL; + struct mc_receiver_1side* mc_rcv1 = NULL; + if(!estimator || !prim_instance || !recv_instance || !rcv + || (side != SSOL_BACK && side != SSOL_FRONT) + || !prim_instance->sample + || !(recv_instance->receiver_mask & (int)side)) + return RES_BAD_ARG; + + /* Check if prim_instance is a primary entity */ + prim_data = estimator_get_primary_entity_data + (&estimator->global_primaries, prim_instance); + if(prim_data == NULL) return RES_BAD_ARG; + + /* realisation count for this primary */ + mc_rcv1 = estimator_get_prim_recv_data(prim_data, recv_instance, side); + if(!prim_data->nb_samples || !mc_rcv1) { + memset(rcv, 0, sizeof(rcv[0])); + } else { + #define SETUP_MC_RESULT(Name) { \ + const double N = (double)prim_data->nb_samples; \ + const struct mc_data* data = &mc_rcv1->Name; \ + rcv->Name.E = data->weight / N; \ + rcv->Name.V = data->sqr_weight / N - rcv->Name.E*rcv->Name.E; \ + rcv->Name.SE = rcv->Name.V > 0 ? sqrt(rcv->Name.V / N) : 0; \ + } (void)0 + SETUP_MC_RESULT(integrated_irradiance); + SETUP_MC_RESULT(absorptivity_loss); + SETUP_MC_RESULT(reflectivity_loss); + SETUP_MC_RESULT(cos_loss); + #undef SETUP_MC_RESULT + } + return RES_OK; +} + +res_T ssol_estimator_get_count (const struct ssol_estimator* estimator, size_t* count) { @@ -133,26 +184,6 @@ ssol_estimator_get_failed_count return RES_OK; } -res_T -ssol_estimator_get_sampled_area - (const struct ssol_estimator* estimator, - double* area) -{ - if (!estimator || !area) return RES_BAD_ARG; - *area = estimator->sampled_area; - return RES_OK; -} - -res_T -ssol_estimator_get_primary_area - (const struct ssol_estimator* estimator, - double* area) -{ - if (!estimator || !area) return RES_BAD_ARG; - *area = estimator->primary_area; - return RES_OK; -} - /******************************************************************************* * Local function ******************************************************************************/ @@ -177,6 +208,7 @@ estimator_create } htable_receiver_init(dev->allocator, &estimator->mc_receivers); + htable_primary_init(dev->allocator, &estimator->global_primaries); SSOL(device_ref_get(dev)); estimator->dev = dev; ref_init(&estimator->ref); diff --git a/src/ssol_estimator_c.h b/src/ssol_estimator_c.h @@ -16,6 +16,7 @@ #ifndef SSOL_ESTIMATOR_C_H #define SSOL_ESTIMATOR_C_H +#include "ssol_device_c.h" #include "ssol_instance_c.h" #include <limits.h> @@ -228,21 +229,142 @@ mc_receiver_copy_and_release #include <rsys/hash_table.h> /******************************************************************************* + * Per primitive MC data + ******************************************************************************/ +struct mc_per_primary_data { + /* global data for this entity */ + struct mc_data cos_loss; + struct mc_data shadow_loss; + double area; + double sun_cos; + size_t nb_samples; + size_t nb_failed; + /* by-receptor data for this entity */ + struct htable_receiver by_receiver; +}; + +static INLINE void +init_mc_per_prim_data + (struct mem_allocator* alloc, + struct mc_per_primary_data* data) +{ + ASSERT(alloc && data); + data->area = 0; + data->sun_cos = 0; + data->cos_loss = MC_DATA_NULL; + data->shadow_loss = MC_DATA_NULL; + data->nb_samples = 0; + data->nb_failed = 0; + htable_receiver_init(alloc, &data->by_receiver); +} + +static INLINE void +release_mc_per_prim_data(struct mc_per_primary_data* data) +{ + ASSERT(data); + htable_receiver_release(&data->by_receiver); +} + +#define PRIM_COPY(Dst, Src) {\ + Dst->area = Src->area;\ + Dst->sun_cos = Src->sun_cos;\ + Dst->nb_failed = Src->nb_failed;\ + Dst->nb_samples = Src->nb_samples;\ + Dst->shadow_loss = Src->shadow_loss;\ + Dst->cos_loss = Src->cos_loss;\ +} (void)0 + +static INLINE res_T +copy_mc_per_prim_data + (struct mc_per_primary_data* dst, + const struct mc_per_primary_data* src) +{ + ASSERT(dst && src); + PRIM_COPY(dst, src); + return htable_receiver_copy(&dst->by_receiver, &src->by_receiver); +} + +static INLINE res_T +copy_and_release_mc_per_prim_data + (struct mc_per_primary_data* dst, + struct mc_per_primary_data* src) +{ + ASSERT(dst && src); + PRIM_COPY(dst, src); + return htable_receiver_copy_and_release(&dst->by_receiver, &src->by_receiver); +} + +static INLINE res_T +copy_and_clear_mc_per_prim_data + (struct mc_per_primary_data* dst, + struct mc_per_primary_data* src) +{ + ASSERT(dst && src); + PRIM_COPY(dst, src); + return htable_receiver_copy_and_clear(&dst->by_receiver, &src->by_receiver); +} + + +static INLINE struct mc_receiver_1side* +mc_per_primary_get_mc_receiver + (struct mc_per_primary_data* primary, + const struct ssol_instance* instance, + const enum ssol_side_flag side) +{ + struct mc_receiver* mc_rcv = NULL; + struct mc_receiver_1side* mc_rcv1 = NULL; + struct mc_receiver mc_rcv_null; + res_T res = RES_OK; + ASSERT(primary && instance); + ASSERT(instance->receiver_mask & (int)side); + + mc_receiver_init(instance->dev->allocator, &mc_rcv_null); + + mc_rcv = htable_receiver_find(&primary->by_receiver, &instance); + if(!mc_rcv) { + res = htable_receiver_set(&primary->by_receiver, &instance, &mc_rcv_null); + if(res != RES_OK) goto error; + mc_rcv = htable_receiver_find(&primary->by_receiver, &instance); + } + mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back; + +exit: + mc_receiver_release(&mc_rcv_null); + return mc_rcv1; +error: + mc_rcv1 = NULL; + goto exit; +} + +#undef PRIM_COPY + +/* Define the htable_primary data structure */ +#define HTABLE_NAME primary +#define HTABLE_KEY const struct ssol_instance* +#define HTABLE_DATA struct mc_per_primary_data +#define HTABLE_DATA_FUNCTOR_INIT init_mc_per_prim_data +#define HTABLE_DATA_FUNCTOR_RELEASE release_mc_per_prim_data +#define HTABLE_DATA_FUNCTOR_COPY copy_mc_per_prim_data +#define HTABLE_DATA_FUNCTOR_COPY_AND_RELEASE copy_and_release_mc_per_prim_data +#define HTABLE_DATA_FUNCTOR_COPY_AND_CLEAR copy_and_clear_mc_per_prim_data +#include <rsys/hash_table.h> + +/******************************************************************************* * Estimator data structure ******************************************************************************/ struct ssol_estimator { size_t realisation_count; size_t failed_count; - /* The implicit MC computations */ + + /* Implicit MC computations */ struct mc_data shadowed; struct mc_data missing; struct mc_data cos_loss; /* TODO compute it */ - /* Per receiver MC */ - struct htable_receiver mc_receivers; + struct htable_receiver mc_receivers; /* Per receiver MC */ + struct htable_primary global_primaries; /* Per sampled MC */ - /* areas */ - double sampled_area, primary_area; + double primary_area; struct ssol_device* dev; ref_T ref; @@ -268,5 +390,27 @@ estimator_get_mc_receiver return side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back; } +static FINLINE struct mc_per_primary_data* +estimator_get_primary_entity_data + (struct htable_primary* primaries, + const struct ssol_instance* instance) +{ + struct mc_per_primary_data* data; + ASSERT(primaries && instance); + if (!instance->sample) return NULL; + data = htable_primary_find(primaries, &instance); + return data; +} + +static FINLINE struct mc_receiver_1side* +estimator_get_prim_recv_data + (struct mc_per_primary_data* primary_data, + const struct ssol_instance* instance, + const enum ssol_side_flag side) +{ + ASSERT(primary_data && instance); + return estimator_get_mc_receiver(&primary_data->by_receiver, instance, side); +} + #endif /* SSOL_ESTIMATOR_C_H */ diff --git a/src/ssol_instance.c b/src/ssol_instance.c @@ -84,10 +84,12 @@ ssol_object_instantiate /* Create the Star-3D instance to ray-trace */ res = s3d_scene_instantiate(object->scn_rt, &instance->shape_rt); if(res != RES_OK) goto error; + instance->shape_rt_area = object->scn_rt_area; /* Create the Star-3D instance to sample */ res = s3d_scene_instantiate(object->scn_samp, &instance->shape_samp); if(res != RES_OK) goto error; + instance->shape_samp_area = object->scn_samp_area; exit: if(out_instance) *out_instance = instance; @@ -185,3 +187,13 @@ ssol_instance_get_id(const struct ssol_instance* instance, uint32_t* id) return RES_OK; } +res_T +ssol_instance_get_area + (const struct ssol_instance* instance, + double* area) +{ + if (!instance || !area) return RES_BAD_ARG;; + /* the area of the 3D surface */ + *area = instance->shape_rt_area; + return RES_OK; +} diff --git a/src/ssol_instance_c.h b/src/ssol_instance_c.h @@ -24,6 +24,7 @@ struct ssol_instance { struct ssol_object* object; /* Instantiated object */ struct s3d_shape* shape_rt; /* Instantiated Star-3D shape to ray-trace */ struct s3d_shape* shape_samp; /* Instantiated Star-3D shape to sample */ + double shape_rt_area, shape_samp_area; double transform[12]; /* Column major 4x3 affine transformation */ int receiver_mask; /* Combination of ssol_side_flag */ int receiver_per_primitive; /* Enable the per primitive receiver */ diff --git a/src/ssol_object.c b/src/ssol_object.c @@ -137,11 +137,13 @@ ssol_object_add_shaded_shape res = s3d_scene_attach_shape(object->scn_rt, shape->shape_rt); if(res != RES_OK) goto error; mask |= BIT(ATTACH_S3D_RT); + object->scn_rt_area += shape->shape_rt_area; /* Add the shape samp to the sampling scene of the object */ res = s3d_scene_attach_shape(object->scn_samp, shape->shape_samp); if(res != RES_OK) goto error; mask |= BIT(ATTACH_S3D_SAMP); + object->scn_samp_area += shape->shape_samp_area; /* Ask for a shaded shape identifier */ i = darray_shaded_shape_size_get(&object->shaded_shapes); @@ -207,6 +209,8 @@ ssol_object_clear(struct ssol_object* obj) darray_shaded_shape_clear(&obj->shaded_shapes); htable_shaded_shape_clear(&obj->shaded_shapes_rt); htable_shaded_shape_clear(&obj->shaded_shapes_samp); + + obj->scn_rt_area = 0; S3D(scene_clear(obj->scn_rt)); S3D(scene_clear(obj->scn_samp)); @@ -214,3 +218,13 @@ ssol_object_clear(struct ssol_object* obj) return RES_OK; } +res_T +ssol_object_get_area +(const struct ssol_object* object, + double* area) +{ + if (!object || !area) return RES_BAD_ARG;; + /* the area of the 3D surface */ + *area = object->scn_rt_area; + return RES_OK; +} +\ No newline at end of file diff --git a/src/ssol_object_c.h b/src/ssol_object_c.h @@ -47,6 +47,7 @@ struct ssol_object { struct s3d_scene* scn_rt; /* RT scene to instantiate */ struct s3d_scene* scn_samp; /* Sampling scene to instantiate */ + double scn_rt_area, scn_samp_area; struct ssol_device* dev; ref_T ref; diff --git a/src/ssol_scene.c b/src/ssol_scene.c @@ -45,7 +45,6 @@ scene_release(ref_T* ref) SSOL(scene_clear(scene)); if(scene->scn_rt) S3D(scene_ref_put(scene->scn_rt)); if(scene->scn_samp) S3D(scene_ref_put(scene->scn_samp)); - if(scene->scn_prim) S3D(scene_ref_put(scene->scn_prim)); if(scene->sun) SSOL(sun_ref_put(scene->sun)); if(scene->atmosphere) SSOL(atmosphere_ref_put(scene->atmosphere)); htable_instance_release(&scene->instances_rt); @@ -84,8 +83,6 @@ ssol_scene_create if(res != RES_OK) goto error; res = s3d_scene_create(dev->s3d, &scene->scn_samp); if(res != RES_OK) goto error; - res = s3d_scene_create(dev->s3d, &scene->scn_prim); - if(res != RES_OK) goto error; exit: if(out_scene) *out_scene = scene; @@ -195,7 +192,6 @@ ssol_scene_clear(struct ssol_scene* scene) htable_instance_clear(&scene->instances_samp); S3D(scene_clear(scene->scn_rt)); S3D(scene_clear(scene->scn_samp)); - S3D(scene_clear(scene->scn_prim)); if (scene->sun) ssol_scene_detach_sun(scene, scene->sun); if (scene->atmosphere) ssol_scene_detach_atmosphere(scene, scene->atmosphere); @@ -280,25 +276,25 @@ res_T scene_create_s3d_views (struct ssol_scene* scn, struct s3d_scene_view** out_view_rt, - struct s3d_scene_view** out_view_samp, - struct s3d_scene_view** out_view_prim) + struct s3d_scene_view** out_view_samp) { struct htable_instance_iterator it, end; struct s3d_scene_view* view_rt = NULL; struct s3d_scene_view* view_samp = NULL; - struct s3d_scene_view* view_prim = NULL; int has_sampled = 0; int has_receiver = 0; res_T res = RES_OK; ASSERT(scn && out_view_rt && out_view_samp); S3D(scene_clear(scn->scn_samp)); - S3D(scene_clear(scn->scn_prim)); htable_instance_clear(&scn->instances_samp); htable_instance_begin(&scn->instances_rt, &it); htable_instance_end(&scn->instances_rt, &end); + scn->sampled_area = 0; + scn->primary_area = 0; + while (!htable_instance_iterator_eq(&it, &end)) { struct ssol_instance* inst = *htable_instance_iterator_data_get(&it); unsigned id; @@ -310,6 +306,9 @@ scene_create_s3d_views if(!inst->sample) continue; + scn->sampled_area += inst->shape_samp_area; + scn->primary_area += inst->shape_rt_area; + /* Note that geometries with virtual material can be sampled without risk * since the solver avoid to shade them and simply pursue the primary ray */ has_sampled = 1; @@ -318,9 +317,6 @@ scene_create_s3d_views res = s3d_scene_attach_shape(scn->scn_samp, inst->shape_samp); if(res != RES_OK) goto error; - /* Attach the instantiated s3d raytraced shape to the s3d primary scene */ - res = s3d_scene_attach_shape(scn->scn_prim, inst->shape_rt); - if(res != RES_OK) goto error; /* Register the instantiated s3d sampling shape */ S3D(shape_get_id(inst->shape_samp, &id)); @@ -346,13 +342,10 @@ scene_create_s3d_views if(res != RES_OK) goto error; res = s3d_scene_view_create(scn->scn_samp, S3D_SAMPLE, &view_samp); if(res != RES_OK) goto error; - res = s3d_scene_view_create(scn->scn_prim, S3D_SAMPLE, &view_prim); - if (res != RES_OK) goto error; exit: *out_view_rt = view_rt; *out_view_samp = view_samp; - *out_view_prim = view_prim; return res; error: S3D(scene_clear(scn->scn_samp)); @@ -365,10 +358,6 @@ error: S3D(scene_view_ref_put(view_samp)); view_samp = NULL; } - if (view_prim) { - S3D(scene_view_ref_put(view_prim)); - view_prim = NULL; - } goto exit; } diff --git a/src/ssol_scene_c.h b/src/ssol_scene_c.h @@ -42,11 +42,12 @@ struct ssol_scene { struct s3d_scene* scn_rt; /* S3D scene to ray trace */ struct s3d_scene* scn_samp; /* S3D scene to sample */ - struct s3d_scene* scn_prim; /* S3D scene of primary objects */ struct ssol_sun* sun; /* Sun of the scene */ struct ssol_atmosphere* atmosphere; /* Atmosphere of the scene */ + double sampled_area, primary_area; /* scene areas */ + struct ssol_device* dev; ref_T ref; }; @@ -57,8 +58,7 @@ extern LOCAL_SYM res_T scene_create_s3d_views (struct ssol_scene* scn, struct s3d_scene_view** view_rt, - struct s3d_scene_view** view_samp, - struct s3d_scene_view** out_view_prim); + struct s3d_scene_view** view_samp); #endif /* SSOL_SCENE_C_H */ diff --git a/src/ssol_shape.c b/src/ssol_shape.c @@ -415,6 +415,43 @@ error: goto exit; } +static double +mesh_compute_area + (const unsigned ntris, + void (*get_indices)(const unsigned itri, unsigned ids[3], void* data), + const unsigned nverts, + void (*get_position)(const unsigned ivert, float position[3], void* data), + void* ctx) +{ + unsigned itri; + double area = 0; + (void)nverts; + + FOR_EACH(itri, 0, ntris) { + float v0[3], v1[3], v2[3]; + double E0[3], E1[3], N[3]; + double V0[3], V1[3], V2[3]; + unsigned IDS[3]; + + get_indices(itri, IDS, ctx); + ASSERT(IDS[0] < nverts); + ASSERT(IDS[1] < nverts); + ASSERT(IDS[2] < nverts); + + get_position(IDS[0], v0, ctx); + get_position(IDS[1], v1, ctx); + get_position(IDS[2], v2, ctx); + d3_set_f3(V0, v0); + d3_set_f3(V1, v1); + d3_set_f3(V2, v2); + d3_sub(E0, V1, V0); + d3_sub(E1, V2, V0); + + area += d3_len(d3_cross(N, E0, E1)); + } + return area * 0.5; +} + /* Setup the Star-3D shape of the quadric to ray-trace, i.e. the clipped 2D * profile of the quadric whose vertices are displaced with respect to the * quadric equation */ @@ -423,12 +460,14 @@ quadric_setup_s3d_shape_rt (const struct ssol_quadric* quadric, const struct darray_double* coords, const struct darray_size_t* ids, - struct s3d_shape* shape) + struct s3d_shape* shape, + double * rt_area) { struct quadric_mesh_context ctx; struct s3d_vertex_data vdata; unsigned nverts; unsigned ntris; + res_T res; ASSERT(quadric && coords && ids && shape); ASSERT(darray_double_size_get(coords)%2 == 0); ASSERT(darray_size_t_size_get(ids)%3 == 0); @@ -458,8 +497,12 @@ quadric_setup_s3d_shape_rt default: FATAL("Unreachable code.\n"); break; } - return s3d_mesh_setup_indexed_vertices + res = s3d_mesh_setup_indexed_vertices (shape, ntris, quadric_mesh_get_ids, nverts, &vdata, 1, &ctx); + if (res != RES_OK) return res; + *rt_area = mesh_compute_area + (ntris, quadric_mesh_get_ids, nverts, quadric_mesh_parabol_get_pos, &ctx); + return RES_OK; } /* Setup the Star-3D shape of the quadric to sample, i.e. the clipped 2D @@ -469,12 +512,14 @@ quadric_setup_s3d_shape_samp (const struct ssol_quadric* quadric, const struct darray_double* coords, const struct darray_size_t* ids, - struct s3d_shape* shape) + struct s3d_shape* shape, + double *samp_area) { struct quadric_mesh_context ctx; struct s3d_vertex_data vdata; unsigned nverts; unsigned ntris; + res_T res; ASSERT(coords && ids && shape); ASSERT(darray_double_size_get(coords)%2 == 0); ASSERT(darray_size_t_size_get(ids)%3 == 0); @@ -490,8 +535,12 @@ quadric_setup_s3d_shape_samp vdata.usage = S3D_POSITION; vdata.type = S3D_FLOAT3; vdata.get = quadric_mesh_plane_get_pos; - return s3d_mesh_setup_indexed_vertices + res = s3d_mesh_setup_indexed_vertices (shape, ntris, quadric_mesh_get_ids, nverts, &vdata, 1, &ctx); + if (res != RES_OK) return res; + *samp_area = mesh_compute_area + (ntris, quadric_mesh_get_ids, nverts, quadric_mesh_plane_get_pos, &ctx); + return RES_OK; } static res_T @@ -979,11 +1028,12 @@ ssol_punched_surface_setup /* Setup the Star-3D shape to ray-trace */ res = quadric_setup_s3d_shape_rt - (psurf->quadric, &coords, &ids, shape->shape_rt); + (psurf->quadric, &coords, &ids, shape->shape_rt, &shape->shape_rt_area); if(res != RES_OK) goto error; /* Setup the Star-3D shape to sample */ - res = quadric_setup_s3d_shape_samp(psurf->quadric, &coords, &ids, shape->shape_samp); + res = quadric_setup_s3d_shape_samp + (psurf->quadric, &coords, &ids, shape->shape_samp, &shape->shape_samp_area); if(res != RES_OK) goto error; exit: @@ -1005,6 +1055,7 @@ ssol_mesh_setup void* data) { struct s3d_vertex_data attrs[SSOL_ATTRIBS_COUNT__]; + void (*get_position)(const unsigned ivert, float position[3], void* data) = NULL; res_T res = RES_OK; unsigned i; @@ -1030,6 +1081,8 @@ ssol_mesh_setup case SSOL_POSITION: attrs[i].usage = SSOL_TO_S3D_POSITION; attrs[i].type = S3D_FLOAT3; + ASSERT(!get_position); + get_position = attrs[i].get; break; case SSOL_NORMAL: attrs[i].usage = SSOL_TO_S3D_NORMAL; @@ -1042,13 +1095,17 @@ ssol_mesh_setup default: FATAL("Unreachable code.\n"); break; } } + ASSERT(get_position); + res = s3d_mesh_setup_indexed_vertices (shape->shape_rt, ntris, get_indices, nverts, attrs, nattribs, data); if(res != RES_OK) goto error; + shape->shape_rt_area = mesh_compute_area(ntris, get_indices, nverts, get_position, data); /* The Star-3D shape to sample is the same of the one to ray-traced */ res = s3d_mesh_copy(shape->shape_rt, shape->shape_samp); if(res != RES_OK) goto error; + shape->shape_samp_area = shape->shape_rt_area; exit: return res; diff --git a/src/ssol_shape_c.h b/src/ssol_shape_c.h @@ -32,6 +32,7 @@ struct ssol_shape { struct s3d_shape* shape_rt; /* Star-3D shape to ray-trace */ struct s3d_shape* shape_samp; /* Star-3D shape to sample */ struct ssol_quadric quadric; + double shape_rt_area, shape_samp_area; struct ssol_device* dev; ref_T ref; diff --git a/src/ssol_solver.c b/src/ssol_solver.c @@ -47,6 +47,7 @@ struct point { const struct ssol_instance* inst; const struct shaded_shape* sshape; + struct mc_per_primary_data* primary; struct s3d_primitive prim; double N[3]; double pos[3]; @@ -63,6 +64,7 @@ struct point { #define POINT_NULL__ { \ NULL, /* Instance */ \ NULL, /* Shaded shape */ \ + NULL, /* Primary data */ \ S3D_PRIMITIVE_NULL__, /* Primitive */ \ {0, 0, 0}, /* Normal */ \ {0, 0, 0}, /* Position */ \ @@ -82,7 +84,8 @@ static int point_init (struct point* pt, struct ssol_scene* scn, - const double sampled_area, + struct htable_primary* prim_data, + struct ssol_estimator* estimator, struct s3d_scene_view* view_samp, struct s3d_scene_view* view_rt, struct ranst_sun_dir* ran_sun_dir, @@ -118,8 +121,8 @@ point_init /* Initialise the Monte Carlo weight */ cos_sun = fabs(d3_dot(pt->N, pt->dir)); - pt->weight = scn->sun->dni * sampled_area * cos_sun; - pt->cos_loss = scn->sun->dni * sampled_area * (1 - cos_sun); + pt->weight = scn->sun->dni * scn->sampled_area * cos_sun; + pt->cos_loss = scn->sun->dni * scn->sampled_area * (1 - cos_sun); pt->absorptivity_loss = pt->reflectivity_loss = 0; /* Retrieve the sampled instance and shaded shape */ @@ -129,6 +132,19 @@ point_init pt->sshape = darray_shaded_shape_cdata_get (&pt->inst->object->shaded_shapes) + id; + /* store primary entity related weights */ + pt->primary = estimator_get_primary_entity_data(prim_data, pt->inst); + if (!pt->primary) { + struct mc_per_primary_data data; + init_mc_per_prim_data(estimator->dev->allocator, &data); + htable_primary_set(prim_data, &pt->inst, &data); + pt->primary = estimator_get_primary_entity_data(prim_data, pt->inst); + ASSERT(pt->primary); + } + pt->primary->cos_loss.weight += pt->cos_loss; + pt->primary->cos_loss.sqr_weight += pt->cos_loss * pt->cos_loss; + pt->primary->nb_samples++; + /* For punched surface, retrieve the sampled position and normal onto the * quadric surface */ if(pt->sshape->shape->type == SHAPE_PUNCHED) { @@ -300,27 +316,27 @@ check_scene(const struct ssol_scene* scene, const char* caller) { ASSERT(scene && caller); - if (!scene->sun) { + if(!scene->sun) { log_error(scene->dev, "%s: no sun attached.\n", caller); return RES_BAD_ARG; } - if (!scene->sun->spectrum) { + if(!scene->sun->spectrum) { log_error(scene->dev, "%s: sun's spectrum undefined.\n", caller); return RES_BAD_ARG; } - if (scene->sun->dni <= 0) { + if(scene->sun->dni <= 0) { log_error(scene->dev, "%s: sun's DNI undefined.\n", caller); return RES_BAD_ARG; } - if (scene->atmosphere) { + if(scene->atmosphere) { int i; ASSERT(scene->atmosphere->type == ATMOS_UNIFORM); i = spectrum_includes (scene->atmosphere->data.uniform.spectrum, scene->sun->spectrum); - if (!i) { + if(!i) { log_error(scene->dev, "%s: sun/atmosphere spectra mismatch.\n", caller); return RES_BAD_ARG; } @@ -328,7 +344,7 @@ check_scene(const struct ssol_scene* scene, const char* caller) return RES_OK; } -static void +static res_T accum_mc_receivers_1side (struct mc_receiver_1side* dst, struct mc_receiver_1side* src) @@ -336,6 +352,7 @@ accum_mc_receivers_1side const struct mc_primitive_1side* src_mc_prim; struct mc_primitive_1side* dst_mc_prim; size_t i; + res_T res = RES_OK; ASSERT(dst && src); #define ACCUM_WEIGHT(Name) { \ @@ -352,6 +369,7 @@ accum_mc_receivers_1side FOR_EACH(i, 0, darray_mc_prim_size_get(&src->mc_prims)) { src_mc_prim = darray_mc_prim_cdata_get(&src->mc_prims) + i; dst_mc_prim = mc_receiver_1side_get_mc_primitive(dst, src_mc_prim->index); + if(!dst_mc_prim) goto error; #define ACCUM_WEIGHT(Name) { \ dst_mc_prim->Name.weight += src_mc_prim->Name.weight; \ dst_mc_prim->Name.sqr_weight += src_mc_prim->Name.sqr_weight; \ @@ -362,6 +380,64 @@ accum_mc_receivers_1side ACCUM_WEIGHT(cos_loss); #undef ACCUM_WEIGHT } +exit: + return res; +error: + goto exit; +} + +static res_T +accum_mc_per_primary_data + (struct mc_per_primary_data* dst, + struct mc_per_primary_data* src) +{ + struct htable_receiver_iterator it, end; + struct mc_receiver mc_rcv_null; + res_T res = RES_OK; + ASSERT(dst && src); + + mc_receiver_init(NULL, &mc_rcv_null); + + #define ACCUM_WEIGHT(Name) { \ + dst->Name.weight += src->Name.weight; \ + dst->Name.sqr_weight += src->Name.sqr_weight; \ + } (void)0 + ACCUM_WEIGHT(cos_loss); + ACCUM_WEIGHT(shadow_loss); + #undef ACCUM_WEIGHT + + dst->nb_samples += src->nb_samples; + dst->nb_failed += src->nb_failed; + + /* dst->by_receiver += src->by_receiver; */ + htable_receiver_begin(&src->by_receiver, &it); + htable_receiver_end(&src->by_receiver, &end); + while(!htable_receiver_iterator_eq(&it, &end)) { + struct mc_receiver* src_mc_rcv = htable_receiver_iterator_data_get(&it); + const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&it); + struct mc_receiver* dst_mc_rcv = htable_receiver_find(&dst->by_receiver, &inst); + htable_receiver_iterator_next(&it); + + if(!dst_mc_rcv) { + res = htable_receiver_set(&dst->by_receiver, &inst, &mc_rcv_null); + if(res != RES_OK) goto error; + dst_mc_rcv = htable_receiver_find(&dst->by_receiver, &inst); + } + + if(inst->receiver_mask & (int)SSOL_FRONT) { + res = accum_mc_receivers_1side(&dst_mc_rcv->front, &src_mc_rcv->front); + if(res != RES_OK) goto error; + } + if(inst->receiver_mask & (int)SSOL_BACK) { + res = accum_mc_receivers_1side(&dst_mc_rcv->back, &src_mc_rcv->back); + if(res != RES_OK) goto error; + } + } +exit: + mc_receiver_release(&mc_rcv_null); + return res; +error: + goto exit; } /******************************************************************************* @@ -375,10 +451,10 @@ ssol_solve FILE* output, struct ssol_estimator** out_estimator) { - struct htable_receiver_iterator it, end; + struct htable_receiver_iterator r_it, r_end; + struct htable_primary_iterator p_it, p_end; struct s3d_scene_view* view_rt = NULL; struct s3d_scene_view* view_samp = NULL; - struct s3d_scene_view* view_prim = NULL; struct ranst_sun_dir* ran_sun_dir = NULL; struct ranst_sun_wl* ran_sun_wl = NULL; struct ssf_bsdf** bsdfs = NULL; @@ -387,11 +463,12 @@ ssol_solve struct mc_data* mc_shadows = NULL; struct mc_data* mc_missings = NULL; struct mc_data* mc_cos_losses = NULL; /* TODO compute it */ + size_t* failed_counts = NULL; /* TODO compute it */ struct htable_receiver* mc_rcvs = NULL; + struct htable_primary* mc_prims = NULL; struct ssol_estimator* estimator = NULL; - float area; int nthreads = 0; - int nrealisations = 0; + int64_t nrealisations = 0; int i = 0; ATOMIC res = RES_OK; ASSERT(nrealisations <= INT_MAX); @@ -404,7 +481,7 @@ ssol_solve /* CL compiler supports OpenMP parallel loop whose indices are signed. The * following line ensures that the unsigned number of realisations does not * overflow the realisation index. */ - if(realisations_count > INT_MAX) { + if(realisations_count > INT64_MAX) { res = RES_BAD_ARG; goto error; } @@ -415,24 +492,20 @@ ssol_solve if(res != RES_OK) goto error; /* Create data structures shared by all threads */ - res = scene_create_s3d_views(scn, &view_rt, &view_samp, &view_prim); + res = scene_create_s3d_views(scn, &view_rt, &view_samp); if(res != RES_OK) goto error; res = sun_create_distributions(scn->sun, &ran_sun_dir, &ran_sun_wl); if(res != RES_OK) goto error; - /* Create the estimator */ - res = estimator_create(scn->dev, scn, &estimator); - if (res != RES_OK) goto error; - S3D(scene_view_compute_area(view_samp, &area)); - estimator->sampled_area = area; - S3D(scene_view_compute_area(view_prim, &area)); - estimator->primary_area = area; - /* Create a RNG proxy from the submitted RNG state */ res = ssp_rng_proxy_create_from_rng (scn->dev->allocator, rng_state, scn->dev->nthreads, &rng_proxy); if(res != RES_OK) goto error; + /* Create the estimator */ + res = estimator_create(scn->dev, scn, &estimator); + if (res != RES_OK) goto error; + /* Create per thread data structures */ #define CREATE(Data) { \ ASSERT(!(Data)); \ @@ -448,6 +521,8 @@ ssol_solve CREATE(mc_missings); CREATE(mc_cos_losses); CREATE(mc_rcvs); + CREATE(mc_prims); + CREATE(failed_counts); #undef CREATE /* Setup per thread data structures */ @@ -459,6 +534,9 @@ ssol_solve htable_receiver_init(scn->dev->allocator, mc_rcvs + i); res = htable_receiver_copy(mc_rcvs + i, &estimator->mc_receivers); if(res != RES_OK) goto error; + htable_primary_init(scn->dev->allocator, mc_prims + i); + res = htable_primary_copy(mc_prims + i, &estimator->global_primaries); + if(res != RES_OK) goto error; } #pragma omp parallel for schedule(static) @@ -470,7 +548,9 @@ ssol_solve struct mc_data* shadow; struct mc_data* missing; struct mc_data* cos_loss; + size_t* nfails; struct htable_receiver* receiver; + struct htable_primary* primary; float org[3], dir[3], range[2] = { 0, FLT_MAX }; const int ithread = omp_get_thread_num(); size_t depth = 0; @@ -486,12 +566,16 @@ ssol_solve cos_loss = mc_cos_losses + ithread; missing = mc_missings + ithread; receiver = mc_rcvs + ithread; + primary = mc_prims + ithread; + nfails = failed_counts + ithread; /* Find a new starting point of the radiative random walk */ - is_lit = point_init(&pt, scn, estimator->sampled_area, view_samp, view_rt, + is_lit = point_init(&pt, scn, primary, estimator, view_samp, view_rt, ran_sun_dir, ran_sun_wl, rng); if(!is_lit) { /* The starting point is not lit */ + pt.primary->shadow_loss.weight += pt.weight; + pt.primary->shadow_loss.sqr_weight += pt.weight; shadow->weight += pt.weight; shadow->sqr_weight += pt.weight * pt.weight; continue; @@ -514,6 +598,8 @@ ssol_solve ATOMIC_SET(&res, res_local); break; } + + /* Per receiver MC accumulation */ mc_rcv1 = estimator_get_mc_receiver(receiver, pt.inst, pt.side); #define ACCUM_WEIGHT(Name, W) { \ mc_rcv1->Name.weight += (W); \ @@ -524,8 +610,8 @@ ssol_solve ACCUM_WEIGHT(reflectivity_loss, pt.reflectivity_loss); ACCUM_WEIGHT(cos_loss, pt.cos_loss); #undef ACCUM_WEIGHT - hit_a_receiver = 1; + /* Per primitive receiver MC accumulation */ if(pt.inst->receiver_per_primitive) { struct mc_primitive_1side* mc_prim; mc_prim = mc_receiver_1side_get_mc_primitive(mc_rcv1, pt.prim.prim_id); @@ -543,6 +629,24 @@ ssol_solve ACCUM_WEIGHT(cos_loss, pt.cos_loss); #undef ACCUM_WEIGHT } + + /* Per-primary/receiver MC accumulation */ + mc_rcv1 = mc_per_primary_get_mc_receiver(pt.primary, pt.inst, pt.side); + #define ACCUM_WEIGHT(Name, W) { \ + mc_rcv1->Name.weight += (W); \ + mc_rcv1->Name.sqr_weight += (W)*(W); \ + } (void)0 + if(!mc_rcv1) { + ATOMIC_SET(&res, RES_MEM_ERR); + break; + } + ACCUM_WEIGHT(integrated_irradiance, pt.weight); + ACCUM_WEIGHT(absorptivity_loss, pt.absorptivity_loss); + ACCUM_WEIGHT(reflectivity_loss, pt.reflectivity_loss); + ACCUM_WEIGHT(cos_loss, pt.cos_loss); + #undef ACCUM_WEIGHT + + hit_a_receiver = 1; } mtl = point_get_material(&pt); @@ -600,7 +704,7 @@ ssol_solve } } - /* Merge per thread estimations */ + /* Merge per thread global MC estimations */ FOR_EACH(i, 0, nthreads) { estimator->shadowed.weight += mc_shadows[i].weight; estimator->shadowed.sqr_weight += mc_shadows[i].sqr_weight; @@ -608,32 +712,52 @@ ssol_solve estimator->missing.sqr_weight += mc_missings[i].sqr_weight; estimator->cos_loss.weight += mc_cos_losses[i].weight; estimator->cos_loss.sqr_weight += mc_cos_losses[i].sqr_weight; + estimator->failed_count += failed_counts[i]; } - /* Merge per thread receiver estimations */ - htable_receiver_begin(&estimator->mc_receivers, &it); - htable_receiver_end(&estimator->mc_receivers, &end); - while(!htable_receiver_iterator_eq(&it, &end)) { - struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&it); - const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&it); - htable_receiver_iterator_next(&it); + /* Merge per thread receiver MC estimations */ + htable_receiver_begin(&estimator->mc_receivers, &r_it); + htable_receiver_end(&estimator->mc_receivers, &r_end); + while(!htable_receiver_iterator_eq(&r_it, &r_end)) { + struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it); + const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it); + htable_receiver_iterator_next(&r_it); FOR_EACH(i, 0, nthreads) { struct mc_receiver* mc_rcv_thread; mc_rcv_thread = htable_receiver_find(mc_rcvs + i, &inst); - if(inst->receiver_mask & (int)SSOL_FRONT) - accum_mc_receivers_1side(&mc_rcv->front, &mc_rcv_thread->front); - if(inst->receiver_mask & (int)SSOL_BACK) - accum_mc_receivers_1side(&mc_rcv->back, &mc_rcv_thread->back); + if(inst->receiver_mask & (int)SSOL_FRONT) { + res = accum_mc_receivers_1side(&mc_rcv->front, &mc_rcv_thread->front); + if(res != RES_OK) goto error; + } + if(inst->receiver_mask & (int)SSOL_BACK) { + res = accum_mc_receivers_1side(&mc_rcv->back, &mc_rcv_thread->back); + if(res != RES_OK) goto error; + } + } + } + + /* Merge per thread sampled instance MC estimations */ + htable_primary_begin(&estimator->global_primaries, &p_it); + htable_primary_end(&estimator->global_primaries, &p_end); + while(!htable_primary_iterator_eq(&p_it, &p_end)) { + struct mc_per_primary_data* pri_data = htable_primary_iterator_data_get(&p_it); + const struct ssol_instance* key = *htable_primary_iterator_key_get(&p_it); + htable_primary_iterator_next(&p_it); + + FOR_EACH(i, 0, nthreads) { + struct mc_per_primary_data* th_data = htable_primary_find(mc_prims + i, &key); + res = accum_mc_per_primary_data(pri_data, th_data); + if(res != RES_OK) goto error; } } + estimator->realisation_count += realisations_count; exit: if(view_rt) S3D(scene_view_ref_put(view_rt)); if(view_samp) S3D(scene_view_ref_put(view_samp)); - if(view_prim) S3D(scene_view_ref_put(view_prim)); if(ran_sun_dir) ranst_sun_dir_ref_put(ran_sun_dir); if(ran_sun_wl) ranst_sun_wl_ref_put(ran_sun_wl); if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); @@ -649,9 +773,14 @@ exit: FOR_EACH(i, 0, nthreads) htable_receiver_release(mc_rcvs + i); sa_release(mc_rcvs); } + if (mc_prims) { + FOR_EACH(i, 0, nthreads) htable_primary_release(mc_prims + i); + sa_release(mc_prims); + } sa_release(mc_shadows); sa_release(mc_missings); sa_release(mc_cos_losses); + sa_release(failed_counts); if(out_estimator) *out_estimator = estimator; return (res_T)res; error: diff --git a/src/test_ssol_by_receiver_integration.c b/src/test_ssol_by_receiver_integration.c @@ -133,6 +133,7 @@ main(int argc, char** argv) #define N__ 10000 #define S_DNI_cos (4 * 1000 * cos(PI / 4)) #define GET_RCV_STATUS ssol_estimator_get_mc_receiver +#define GET_PRIM_X_RCV_STATUS ssol_estimator_get_primary_entity_x_receiver_status CHECK(ssol_solve(scene, rng, N__, NULL, &estimator1), RES_OK); CHECK(GET_RCV_STATUS(estimator1, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", @@ -150,9 +151,10 @@ main(int argc, char** argv) printf("Ir(target) = %g +/- %g\n", mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); CHECK(eq_eps(mc_rcv.integrated_irradiance.E, S_DNI_cos, S_DNI_cos * 1e-1), 1); -#undef N__ -#undef S_DNI_cos -#undef GET_RCV_STATUS + CHECK(GET_PRIM_X_RCV_STATUS(estimator1, heliostat, target, SSOL_FRONT, &mc_rcv), RES_OK); + printf("Ir(heliostat=>target) = %g +/- %g", + mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE); + CHECK(eq_eps(mc_rcv.integrated_irradiance.E, S_DNI_cos, S_DNI_cos * 1e-1), 1); /* Free data */ CHECK(ssol_instance_ref_put(heliostat), RES_OK);