star-phor

Radiative transfer solver for photoreactors.
git clone https://www.edstar.cnrs.fr/git/star-phor.git
Log | Files | Refs | README | LICENSE

commit be688d494a1d919d71f9683f1193386d1b9402ca
parent a51410c83da1e972f0c26bb9f7316e566d0bc822
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date:   Wed, 14 Jan 2026 15:13:16 +0100

Continue refactoring of realization function

This commit makes a series of minor modifications in the
compute_MVREA_realization function in order to improve semantic clarity
(sort local variables, explicit ambiguous inputs and outputs in function
calls, simplify error handling...).

Also, the wavelength is no longer a member variable of the ray structure,
but is a local variable in the realization function. This makes the ray a
purely geometry-related structure, used exclusively in the intersection
search procedure.

Diffstat:
Msrc/sphor_compute_mvrea.c | 177++++++++++++++++++++++++++++++++++++++-----------------------------------------
Msrc/sphor_interface.c | 24------------------------
Msrc/sphor_interface.h | 9---------
3 files changed, 86 insertions(+), 124 deletions(-)

diff --git a/src/sphor_compute_mvrea.c b/src/sphor_compute_mvrea.c @@ -508,8 +508,8 @@ error: static res_T setup_bsdf_from_refractive_indices (struct sphor* sphor, - struct ray* ray, struct intersection* intersection, + double wavelength, struct ssf_bsdf** out_bsdf) { double n_real_i = 0, n_imag_i = 0; @@ -519,17 +519,16 @@ setup_bsdf_from_refractive_indices res_T res = RES_OK; ASSERT(NULL != sphor); - ASSERT(NULL != ray); ASSERT(NULL != intersection); /* Get refractive index of the medium in the incident side */ res = primitive_get_in_refraction_index - (sphor, primitive, ray->wavelength, &n_real_i, &n_imag_i); + (sphor, primitive, wavelength, &n_real_i, &n_imag_i); if (RES_OK != res){ goto error; } /* Get refractive index of the medium in the transmitted side */ res = primitive_get_out_refraction_index - (sphor, primitive, ray->wavelength, &n_real_t, &n_imag_t); + (sphor, primitive, wavelength, &n_real_t, &n_imag_t); if (RES_OK != res){ goto error; } /* Sample an interaction type and a new direction using ssf */ @@ -555,6 +554,7 @@ sample_interface_ray_interaction struct ssp_rng* rng, struct ray* ray, struct intersection* intersection, + double wavelength, double dir[3], enum ray_interface_interaction_type* interaction_type) { @@ -605,7 +605,7 @@ sample_interface_ray_interaction * index of each side of the interface to construct the BRDF */ else { res = setup_bsdf_from_refractive_indices - (sphor, ray, intersection, &bsdf); + (sphor, intersection, wavelength, &bsdf); if (RES_OK != res){ goto error; } } @@ -664,56 +664,70 @@ compute_MVREA_realization struct htable_accum2id* accum2id, struct darray_accum* accums) { - double dir[3] = {0}; /* Current direction in path sampling */ - double response_function = 0; - double ka = 0; - double s = 0; /* Random number */ + /* Star-Phor Input */ + struct sphin_sensor_volume* sensor_volume = NULL; + struct sphin_sensor_surface* sensor_surface = NULL; + struct sphin_volume* volume = NULL; + + /* The sampled source */ + struct source_view* source_view = NULL; + + /* Ray */ enum ray_interface_interaction_type ray_interface_interaction_type = RAY_INTERFACE_INTERACTION_NONE__; struct intersection intersection = INTERSECTION_NULL; struct ray ray = RAY_DEFAULT; - struct source_view* source_view = NULL; - struct sphin_sensor_volume* sensor_volume = NULL; - struct sphin_sensor_surface* sensor_surface = NULL; - struct sphin_volume* volume = NULL; + + /* Miscellaneous */ + double dir[3] = {0}; /* Sampled direction during path sampling */ + double response_function = 0; + double ka = 0; + double free_path = 0; + double wavelength = 0; + int stop = 0; + + /* For error management */ res_T res = RES_OK; + #define CALL(Function) { \ + res = Function; \ + if(RES_OK != res) goto error; \ + } (void)0 + /* Pre-conditions */ ASSERT(NULL != sphor); ASSERT(NULL != rng); ASSERT(NULL != accum2id); ASSERT(NULL != accums); /* Sample a random position from a source in the whole scene */ - res = sample_source_position - (sphor, rng, &source_view, &ray.origin_on_prim, ray.origin); - if (RES_OK != res){ goto error; } + CALL(sample_source_position(sphor, + /* in: */ rng, + /* out: */ &source_view, &ray.origin_on_prim, ray.origin)); /* Sample a direction according to the source direction distribution */ - res = source_surface_sample_direction - (sphor, rng, &ray.origin_on_prim, ray.direction); - if (RES_OK != res){ goto error; } + CALL(source_surface_sample_direction(sphor, + /* in: */ rng, &ray.origin_on_prim, + /* out: */ ray.direction)); /* Sample a wavelength according to the source emission spectrum */ - res = source_sample_wavelength - (sphor, rng, source_view, &ray.wavelength); - if (RES_OK != res){ goto error; } + CALL(source_sample_wavelength + (sphor, rng, source_view, &wavelength)); + /* Retrieve the volume defined by the sampled primitive */ volume = volume_get_from_primitive(sphor, &ray.origin_on_prim.primitive); /* Retrieve the ka of the present volume */ - res = volume_compute_total_ka(sphor, volume, ray.wavelength, &ka); - if (RES_OK != res){ goto error; } + CALL(volume_compute_total_ka(sphor, volume, wavelength, &ka)); /* While no absorption takes place and an interface is found */ - for(;;) { + while(!stop) { /* Sample the free path length to absorption from an exponential * distribution with rate parameter ka */ - s = ssp_ran_exp(rng, ka); + free_path = ssp_ran_exp(rng, ka); /* Trace a ray in the scene from the sampled primitive in the sampled * direction and retrive the hit distance */ - res = trace_ray(sphor, &ray, &intersection); - if (RES_OK != res){ goto error; } + CALL(trace_ray(sphor, &ray, &intersection)); /* If there is no intersection in the sampled direction from the sampled * position, the photon is lost and it counts zero in the weight. */ @@ -722,11 +736,10 @@ compute_MVREA_realization /* If the distance between the ray origin and the intersection is bigger * than the absorption free path length, an absorption takes place (check * if the volume is a sensor and return the weight accordingly) */ - if (s < intersection.distance) { - /* Check if the volume is a sensor */ - res = sphin_volume_get_sensor(volume, &sensor_volume); - if (RES_OK != res){ goto error; } + if (free_path < intersection.distance) { + CALL(sphin_volume_get_sensor(volume, &sensor_volume)); + /* Check if the volume is a sensor */ if (NULL == sensor_volume) { /* If the volume is not a sensor, the absorbed photon does not count in * the weight. The weight is implictly 0 */ @@ -734,27 +747,28 @@ compute_MVREA_realization /* If the volume is a sensor, compute the weight using the response * function */ - /* TODO res = update_mc_weights(...); */ + /* TODO res = update_mc_weights_MVREA_volume(...,mc_weight,...); */ struct accum_key accum_key = ACCUM_KEY_NULL; size_t prop_rad_count = 0; size_t volume_id = SIZE_MAX; size_t i_proprad = 0; - double vol = 0; + double vol_size = 0; double tot_pow = 0; double mc_weight = 0; - res = sphin_sensor_volume_get_response_function - (sensor_volume, &response_function); - if (RES_OK != res){ goto error; } - res = sphin_volume_compute_total_volume(volume, &vol); + /* Get the total three-dimensional space occupied by the volume */ + res = sphin_volume_compute_total_size(volume, &vol_size/*[m^3]*/); if (RES_OK != res){ goto error; } res = sphin_config_get_total_surface_power(sphor->config, &tot_pow); if (RES_OK != res){ goto error; } + + mc_weight = tot_pow / vol_size; + res = sphin_sensor_volume_get_response_function + (sensor_volume, &response_function); + if (RES_OK != res){ goto error; } res = sphin_volume_get_prop_rad_count(volume, &prop_rad_count); if (RES_OK != res){ goto error; } - mc_weight = tot_pow / vol; - volume_id = intersection.position.primitive.interface->volumes [intersection.position.primitive.side]; @@ -777,81 +791,62 @@ compute_MVREA_realization struct sphin_prop_rad* prop_rad = NULL; res = sphin_volume_get_prop_rad(volume, i_proprad, &prop_rad); if (RES_OK != res){ goto error; } - res = prop_rad_compute_ka(prop_rad, ray.wavelength, &ka_i); + res = prop_rad_compute_ka(prop_rad, wavelength, &ka_i); if (RES_OK != res){ goto error; } accum_key.component_id = i_proprad; - mc_weight = ka_i / ka * tot_pow / vol; + mc_weight = ka_i / ka * tot_pow / vol_size; update_accum(accum2id, &accum_key, mc_weight, accums); } } /* Absorption => stop the path */ - break; - } - - /* Else, the distance between the ray origin and the intersection is smaller - * than the absorption free path length. The photon intersects a primitive - * in the scene view. Sample the type of interaction that the photon has - * with the interface as well as the new direction sampled */ - res = sample_interface_ray_interaction - (sphor, rng, &ray, &intersection, dir, &ray_interface_interaction_type); - if (RES_OK != res){ goto error; } - - int stop = 0; - switch (ray_interface_interaction_type) { - /* An absorption took place, the weight is zero and the path stops */ - case RAY_INTERFACE_INTERACTION_ABSORPTION: + stop = 1; + + } else { + /* Else, the distance between the ray origin and the intersection is + * smaller than the absorption free path length. The photon intersects a + * primitive in the scene view. Sample the type of interaction that the + * photon has with the interface as well as the new direction sampled */ + CALL(sample_interface_ray_interaction(sphor, + /* in: */ rng, &ray, &intersection, wavelength, + /* out: */ dir, &ray_interface_interaction_type)); + + if(RAY_INTERFACE_INTERACTION_ABSORPTION == ray_interface_interaction_type) { /* TODO: verify if surface is sensor and update its accum - res = primitive_get_sensor_surface(sphor, &intersection.position.primitive, &sensor_surface); - if (RES_OK != res){ goto error; } */ + res = primitive_get_sensor_surface(sphor, + &intersection.position.primitive, &sensor_surface); if (RES_OK != + res){ goto error; } */ /* If the volume is a sensor, compute the weight using the response * function */ - if (NULL != sensor_surface) - { - /* res = update_mc_weights(...= */ - } - stop = 1; - break; + if (NULL != sensor_surface){ + /* res = update_mc_weights_surface_MVREA(...= */ + } - /* Else, if a reflection takes place, sample a new direction according to - * the photon arriving direction, geometry normak and surface brdf - * properties */ - case RAY_INTERFACE_INTERACTION_REFLECTION: - break; + /* Absorption => stop the path */ + stop = 1; - /* If the ray traverses the interface, it simply continue in the same - * direction as before */ - case RAY_INTERFACE_INTERACTION_TRANSMISSION: + } else if(RAY_INTERFACE_INTERACTION_TRANSMISSION == ray_interface_interaction_type) { /* Continue the path on the other side of the primitive.*/ intersection.position.primitive.side = !intersection.position.primitive.side; + CALL(ray_update(/* in/out: */ &ray, /* in: */ &intersection, dir)); + /* Update the current medium */ volume = volume_get_from_primitive (sphor, &intersection.position.primitive); - /* Retrieve the ka of the present volume */ - res = volume_compute_total_ka(sphor, volume, ray.wavelength, &ka); - if (RES_OK != res){ goto error; } - break; - - /* Handle the possible exceptions */ - case RAY_INTERFACE_INTERACTION_NONE__: - res = RES_BAD_ARG; - goto error; - default: - FATAL("Unreachable code\n"); - break; - } - if (stop) { break; } - - /* Update ray with the intersection information */ - res = ray_update(&ray, &intersection, dir); - if (RES_OK != res){ goto error; } + CALL(volume_compute_total_ka(sphor, volume, wavelength, &ka)); + } else if(RAY_INTERFACE_INTERACTION_REFLECTION == ray_interface_interaction_type) { + CALL(ray_update(/* in/out: */ &ray, /* in: */ &intersection, dir)); + } + } } + #undef CALL + exit: return res; error: diff --git a/src/sphor_interface.c b/src/sphor_interface.c @@ -384,27 +384,3 @@ ray_update d3_set(ray->direction, dir); return res; } - -res_T -interface_get_side_id - (const double dir[3], /* Direction of arrival of the hit into the geometry */ - const struct s3d_hit* hit, - enum sphin_side* side) -{ - double normal[3] = {0}; - res_T res = RES_OK; - - ASSERT(NULL != dir); - ASSERT(NULL != side); - - /* hit->normal: unormalized hit normal that uses the left hand convention */ - d3_normalize(normal, d3_set_f3(normal, hit->normal)); - - /* Test the angle between the two vectors dir and normal using the dot - * product. If it is negative, the S3D normal points to the hemisphere of the - * arriving hit, which corresponds to the BACK side in the sphin convention. - * Otherwise, it corresponds to the front side. */ - *side = (d3_dot(dir, normal) < 0 ) ? SPHIN_SIDE_BACK : SPHIN_SIDE_FRONT; - - return res; -} diff --git a/src/sphor_interface.h b/src/sphor_interface.h @@ -86,7 +86,6 @@ static const struct intersection INTERSECTION_NULL = INTERSECTION_NULL__; struct ray { double origin[3]; double direction[3]; - double wavelength; double range[2]; /* Origin of the ray when located on a primitive. @@ -96,7 +95,6 @@ struct ray { #define RAY_DEFAULT__ { \ {0,0,0}, \ {0,0,0}, \ - DBL_MAX, \ {0, DBL_MAX}, \ PRIMITIVE_POS_NULL__ \ } @@ -148,11 +146,4 @@ ray_update const struct intersection* intersection, const double dir[3]); -/* TODO this should disappear */ -extern LOCAL_SYM res_T -interface_get_side_id - (const double dir[3], - const struct s3d_hit* hit, - enum sphin_side* side); - #endif /* SPHOR_INTERFACE_H */