star-phor

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

commit a3f743987d61709209ed359eba7b2c1cf7d55208
parent 934e3783f95d2e5b0f6aa5a3bb05f645f4404554
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date:   Fri, 27 Feb 2026 08:40:54 +0100

Prepare internal API for extended interface optical handling

This commit refactors the way interface optical properties are handled
during ray–primitive intersections in path propagation, in order to
comply with the new star-phor-input configuration specification and the
updated sphin library, which now support extended BSDF features,
including more detailed BRDF models and explicit BTDF definitions.

The configuration file specification now explicitly enforces the
conceptual separation of the BSDF into BRDF and BTDF components. This
differs from the current Star Scattering Functions (ssf) library, which
is used to sample interface scattering in star-phor and does not fully
reflect this separation at the API level. As a consequence, an
intermediate abstraction layer has been introduced. This layer:
 1. Samples the interaction type (reflection, transmission, absorption)
    according to the interface optical properties.
 2. Instantiates the appropriate ssf BSDF object.
 3. Samples the outgoing direction.

Since ssf does not yet cover the complete set of BRDF and BTDF models
supported by star-phor, some models will need to be implemented
internally in future work.

The BSDF sampling code, which was previously embedded in
sphor_compute_mvrea.c, has been extracted into a dedicated compilation
unit (sphor_ran_bsdf.c). This change reinforces the star-phor
architectural principle of clearly separated abstraction layers.

New functions have also been introduced to query reflectivity and
transmissivity from primitive properties. These quantities are handled
as spectral properties and evaluated at the current wavelength. The
analytical test has been updated to consider the reflectivity a spectral
property.

Diffstat:
M.gitignore | 1+
MMakefile | 4+++-
Msrc/sphor_compute_mvrea.c | 237++++---------------------------------------------------------------------------
Msrc/sphor_interface.c | 286+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/sphor_interface.h | 22++++++++++++++++++++++
Asrc/sphor_ran_bsdf.c | 315+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/sphor_ran_bsdf.h | 52++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/test_sphor_MVREA_analytical1.c | 19+++++++++++++++++--
8 files changed, 706 insertions(+), 230 deletions(-)

diff --git a/.gitignore b/.gitignore @@ -3,6 +3,7 @@ .config .test file.txt +rho.txt spec.txt *.pc *.so diff --git a/Makefile b/Makefile @@ -42,6 +42,7 @@ SRC_LIB =\ src/sphor_compute_mvrea.c\ src/sphor_config.c\ src/sphor_interface.c\ + src/sphor_ran_bsdf.c\ src/sphor_ran_geometry.c\ src/sphor_ran_source.c\ src/sphor_sources.c @@ -226,5 +227,6 @@ test_sphor_MVREA_analytical1\ clean_test: rm -f $(TEST_DEP) $(TEST_OBJ) $(TEST_TGT) \ - input output cube.stl cube_back.stl cube_front.stl spec.txt + input output cube.stl cube_back.stl cube_front.stl \ + spec.txt rho.txt for i in $(TEST_SRC); do rm -f "$$(basename "$${i}" ".c")"; done diff --git a/src/sphor_compute_mvrea.c b/src/sphor_compute_mvrea.c @@ -26,6 +26,7 @@ #include "sphor_accum.h" #include "sphor_compute_mvrea.h" #include "sphor_interface.h" +#include "sphor_ran_bsdf.h" #include "sphor_ran_source.h" #include <rsys/double3.h> @@ -41,20 +42,12 @@ struct sphin_prop_rad; struct sphin_volume; struct ssf_bsdf; -struct ssf_fresnel; struct ssp_rng; /* Syntactic sugar */ #define ALL_COMPONENTS INVALID_ID #define ALL_SENSORS INVALID_ID -enum ray_interface_interaction_type { - RAY_INTERFACE_INTERACTION_REFLECTION, /* 0 */ - RAY_INTERFACE_INTERACTION_TRANSMISSION, /* 1 */ - RAY_INTERFACE_INTERACTION_ABSORPTION, /* 2 */ - RAY_INTERFACE_INTERACTION_NONE__ -}; - /******************************************************************************* * Helper functions ******************************************************************************/ @@ -433,216 +426,6 @@ volume_get_from_primitive } static res_T -setup_bsdf_from_sphin_brdf - (struct sphor* sphor, - struct ray* ray, - struct intersection* intersection, - struct ssf_fresnel** out_fresnel, - struct ssf_bsdf** out_bsdf) -{ - double reflectivity = 0; - struct primitive* primitive = &intersection->position.primitive; - struct sphin_brdf* brdf = NULL; - struct ssf_bsdf* bsdf = NULL; - struct ssf_fresnel* fresnel = NULL; - enum sphin_brdf_type brdf_type = SPHIN_BRDF_NONE__; - res_T res = RES_OK; - - ASSERT(NULL != sphor); - ASSERT(NULL != ray); - ASSERT(NULL != intersection); - - (void)ray; - - primitive = (struct primitive*)&intersection->position.primitive; - - res = primitive_get_brdf(sphor, primitive, &brdf); - if (RES_OK != res) { goto error; } - if (NULL == brdf) { res = RES_BAD_ARG; goto error; } - - res = sphin_brdf_get_type(brdf, &brdf_type); - if (RES_OK != res) { goto error; } - - switch (brdf_type) { - case SPHIN_BRDF_LAMBERT: - res = sphin_brdf_lambertian_get_reflectivity(brdf, &reflectivity); - if (RES_OK != res) { goto error; } - - res = ssf_bsdf_create - (sphor->allocator, &ssf_lambertian_reflection, &bsdf); - if (RES_OK != res) { goto error; } - - res = ssf_lambertian_reflection_setup(bsdf, reflectivity); - if (RES_OK != res) { goto error; } - - break; - case SPHIN_BRDF_SPECULAR: - res = sphin_brdf_specular_get_reflectivity(brdf, &reflectivity); - if (RES_OK != res) { goto error; } - - res = ssf_bsdf_create - (sphor->allocator, &ssf_specular_reflection, &bsdf); - if (RES_OK != res) { goto error; } - - res = ssf_fresnel_create - (sphor->allocator, &ssf_fresnel_constant, &fresnel); - if (RES_OK != res) { goto error; } - - res = ssf_fresnel_constant_setup(fresnel, reflectivity); - if (RES_OK != res) { goto error; } - - res = ssf_specular_reflection_setup(bsdf, fresnel); - if (RES_OK != res) { goto error; } - break; - - case SPHIN_BRDF_NONE__: - res = RES_BAD_ARG; - goto error; - default: - FATAL("Unreachable code\n"); - break; - } - -exit: - *out_fresnel = fresnel; - *out_bsdf = bsdf; - return res; -error: - if (NULL != bsdf) { SSF(bsdf_ref_put(bsdf)); } - if (NULL != fresnel) { SSF(fresnel_ref_put(fresnel)); } - goto exit; -} - -static res_T -setup_bsdf_from_refractive_indices - (struct sphor* sphor, - struct intersection* intersection, - double wavelength, - struct ssf_bsdf** out_bsdf) -{ - double n_real_i = 0, n_imag_i = 0; - double n_real_t = 0, n_imag_t = 0; - struct primitive* primitive = &intersection->position.primitive; - struct ssf_bsdf* bsdf = NULL; - res_T res = RES_OK; - - ASSERT(NULL != sphor); - ASSERT(NULL != intersection); - - /* Get refractive index of the medium in the incident side */ - res = primitive_get_in_refraction_index - (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, wavelength, &n_real_t, &n_imag_t); - if (RES_OK != res) { goto error; } - - /* Sample an interaction type and a new direction using ssf */ - res = ssf_bsdf_create - (sphor->allocator, &ssf_specular_dielectric_dielectric_interface, &bsdf); - if (RES_OK != res) { goto error; } - - res = ssf_specular_dielectric_dielectric_interface_setup - (bsdf, n_real_i, n_real_t); - if (RES_OK != res) { goto error; } - -exit: - *out_bsdf = bsdf; - return res; -error: - if (NULL != bsdf) { SSF(bsdf_ref_put(bsdf)); } - goto exit; -} - -static res_T -sample_interface_ray_interaction - (struct sphor* sphor, - struct ssp_rng* rng, - struct ray* ray, - struct intersection* intersection, - double wavelength, - double dir[3], - enum ray_interface_interaction_type* interaction_type) -{ - double invert_normal = 0; - double normal[3] = {0}; - double wo[3] = {0}; - double pdf = 0; - double reflectivity = 0; - double sample[3] = {0}; - int flag = 0; - struct primitive* primitive = NULL; - struct sphin_brdf* brdf = NULL; - struct ssf_bsdf* bsdf = NULL; - struct ssf_fresnel* fresnel = NULL; - res_T res = RES_OK; - - ASSERT(NULL != sphor); - ASSERT(NULL != rng); - ASSERT(NULL != ray); - ASSERT(NULL != intersection); - - primitive = &intersection->position.primitive; - - /* Ensure that the ray direction and the geometry normal are compatible with - * the star-sf library. */ - d3_muld(wo, ray->direction, -1); - d3_set(normal, primitive->normal); - /* Ensure normal and direction point to the same hemisphere */ - invert_normal = (d3_dot(normal, wo) > 0 ) ? 1. : -1.; - d3_normalize(normal, d3_muld(normal, normal, invert_normal)); - - /* Verify if one of the surfaces defined in the intersected interface - * side has a BRDF defined */ - res = primitive_get_brdf(sphor, primitive, &brdf); - if (RES_OK != res) { goto error; } - - if (NULL != brdf) { - /* If the interface side is a sphin_surface with a BRDF declared, - * verify if a reflection takes place based on the reflectivity of the - * surface */ - res = setup_bsdf_from_sphin_brdf - (sphor, ray, intersection, &fresnel, &bsdf); - if (RES_OK != res) { goto error; } - - } - - /* If the interface side does not have a BRDF (i.e., it defines only a - * sphin_volume or a sphin_surface without a BRDF), we use the refractive - * index of each side of the interface to construct the BRDF */ - else { - res = setup_bsdf_from_refractive_indices - (sphor, intersection, wavelength, &bsdf); - if (RES_OK != res) { goto error; } - } - - reflectivity = ssf_bsdf_sample(bsdf, rng, wo, normal, sample, &flag, &pdf); - - if (NULL != brdf) { - if (ssp_rng_canonical(rng) < reflectivity) { - *interaction_type = RAY_INTERFACE_INTERACTION_REFLECTION; } - else { - *interaction_type = RAY_INTERFACE_INTERACTION_ABSORPTION; } - } - else if (flag & SSF_TRANSMISSION) { - *interaction_type = RAY_INTERFACE_INTERACTION_TRANSMISSION; } - else if (flag & SSF_REFLECTION) { - *interaction_type = RAY_INTERFACE_INTERACTION_REFLECTION; } - else { res = RES_BAD_ARG; } - - d3_set(dir, sample); - -exit: - if (NULL != fresnel) { SSF(fresnel_ref_put(fresnel)); } - if (NULL != bsdf) { SSF(bsdf_ref_put(bsdf)); } - return res; -error: - goto exit; -} - -static res_T MVREA_update_volume_weights (struct sphor* sphor, struct intersection* intersection, @@ -834,8 +617,8 @@ compute_MVREA_realization struct source_view* source_view = NULL; /* Ray */ - enum ray_interface_interaction_type ray_interface_interaction_type = - RAY_INTERFACE_INTERACTION_NONE__; + enum intersection_ray_interaction_type intersection_ray_interaction_type = + INTERSECTION_RAY_INTERACTION_NONE__; struct intersection intersection = INTERSECTION_NULL; struct ray ray = RAY_DEFAULT; @@ -920,10 +703,10 @@ compute_MVREA_realization * 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)); + /* out: */ dir, &intersection_ray_interaction_type)); - if (RAY_INTERFACE_INTERACTION_ABSORPTION == - ray_interface_interaction_type) { + if (INTERSECTION_RAY_INTERACTION_ABSORPTION == + intersection_ray_interaction_type) { CALL(primitive_get_sensor_surface (sphor, &intersection.position.primitive, &sensor_surface)); @@ -938,8 +721,8 @@ compute_MVREA_realization /* Absorption by a surface => stop the path */ stop = 1; } - else if (RAY_INTERFACE_INTERACTION_TRANSMISSION == - ray_interface_interaction_type) { + else if (INTERSECTION_RAY_INTERACTION_TRANSMISSION == + intersection_ray_interaction_type) { /* Continue the path on the other side of the primitive.*/ intersection.position.primitive.side = !intersection.position.primitive.side; @@ -952,8 +735,8 @@ compute_MVREA_realization CALL(volume_compute_total_ka(sphor, volume, wavelength, &ka)); } - else if (RAY_INTERFACE_INTERACTION_REFLECTION == - ray_interface_interaction_type) { + else if (INTERSECTION_RAY_INTERACTION_REFLECTION == + intersection_ray_interaction_type) { CALL(ray_update(/* in/out: */ &ray, /* in: */ &intersection, dir)); } } diff --git a/src/sphor_interface.c b/src/sphor_interface.c @@ -33,6 +33,7 @@ #include <rsys/rsys.h> #include <star/s3d.h> #include <star/sphin.h> +#include <star/ssf.h> struct sphin_brdf; struct sphin_sensor_surface; @@ -101,6 +102,61 @@ error: } res_T +primitive_get_btdf + (struct sphor* sphor, + const struct primitive* primitive, + struct sphin_btdf** out_btdf) +{ + char* surface_name = NULL; + size_t i = 0; + const struct interface* interface = NULL; + struct sphin_btdf* btdf = NULL; + struct sphin_btdf* found_btdf = NULL; + struct sphin_surface* surface = NULL; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != primitive); + ASSERT(NULL != out_btdf); + + interface = primitive->interface; + + FOR_EACH(i, 0, (size_t)interface->surface_count[(size_t)primitive->side]) { + res = sphin_config_get_surface + (sphor->config, + interface->surfaces[(size_t)primitive->side][i], + &surface); + if (RES_OK != res) { goto error; } + + res = sphin_surface_get_btdf(surface, &btdf); + if (RES_OK != res) { goto error; } + if (NULL != btdf){ + if (NULL == found_btdf) { + sphin_surface_get_name(surface, &surface_name); + found_btdf = btdf; /* Keep first non-NULL */ + } else { + /* Found second non-NULL BTDF - this is an error/conflict */ + char* last_surface_name = NULL; + sphin_surface_get_name(surface, &last_surface_name); + ERROR(sphor, + "Multiple BTDFs defined for same interface side. " + "Conflicting surfaces: '%s' and '%s'.\n", + surface_name, last_surface_name); + res = RES_BAD_ARG; + goto error; + } + } + } + +exit: + *out_btdf = btdf; + return res; +error: + found_btdf = NULL; + goto exit; +} + +res_T primitive_get_source (struct sphor* sphor, const struct primitive* primitive, @@ -318,6 +374,236 @@ primitive_get_out_refraction_index } res_T +primitive_get_reflectivity + (struct sphor* sphor, + struct primitive* primitive, + double wavelength, + double cos_theta, + double* reflectivity) +{ + double eta_i = 0, k_i = 0; + double eta_t = 0, k_t = 0; + double refrel = 0; + enum sphin_brdf_reflectivity_type reflectivity_type = + SPHIN_BRDF_REFLECTIVITY_NONE__; + struct sphin_brdf* brdf = NULL; + struct sphin_spectral_property* reflectivity_spectrum = NULL; + struct ssf_fresnel* fresnel = NULL; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != primitive); + ASSERT(NULL != reflectivity); + + res = primitive_get_brdf(sphor, primitive, &brdf); + if (RES_OK != res) { goto error; } + + if (NULL == brdf) { + res = RES_BAD_ARG; + goto error; + } + + res = sphin_brdf_get_reflectivity_type(brdf, &reflectivity_type); + if (RES_OK != res) { goto error; } + + switch (reflectivity_type) { + case SPHIN_BRDF_REFLECTIVITY_TABULATED: + + res = sphin_brdf_get_reflectivity_value(brdf, &reflectivity_spectrum); + if (RES_OK != res) { goto error; } + + res = sphin_spectral_property_interpolate_at_wavelength + (reflectivity_spectrum, wavelength, + SPHIN_INTERPOLATION_LINEAR, &refrel); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_create + (sphor->allocator, &ssf_fresnel_constant, &fresnel); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_constant_setup(fresnel, refrel); + if (RES_OK != res) { goto error; } + + break; + + case SPHIN_BRDF_REFLECTIVITY_FRESNEL_DIELECTRIC: + + /* Get refractive index of the medium in the incident side */ + res = primitive_get_in_refraction_index + (sphor, primitive, wavelength, &eta_i, &k_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, wavelength, &eta_t, &k_t); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_create + (sphor->allocator, &ssf_fresnel_dielectric_dielectric, &fresnel); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_dielectric_dielectric_setup(fresnel, eta_i, eta_t); + if (RES_OK != res) { goto error; } + + break; + + case SPHIN_BRDF_REFLECTIVITY_FRESNEL_DIELECTRIC_CONDUCTOR: + + /* Get refractive index of the medium in the incident side */ + res = primitive_get_in_refraction_index + (sphor, primitive, wavelength, &eta_i, &k_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, wavelength, &eta_t, &k_t); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_create + (sphor->allocator, &ssf_fresnel_dielectric_conductor, &fresnel); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_dielectric_conductor_setup(fresnel, eta_i, eta_t, k_t); + if (RES_OK != res) { goto error; } + + break; + + case SPHIN_BRDF_REFLECTIVITY_NONE__: + res = RES_BAD_ARG; + goto error; + + default: + FATAL("Unreachable code\n"); + break; + } + + *reflectivity = ssf_fresnel_eval(fresnel, cos_theta); + +exit: + if (NULL != fresnel) { + ssf_fresnel_ref_put(fresnel); + } + return res; +error: + goto exit; +} + +res_T +primitive_get_transmissivity + (struct sphor* sphor, + struct primitive* primitive, + double wavelength, + double cos_theta, + double* transmissivity) +{ + double eta_i = 0, k_i = 0; + double eta_t = 0, k_t = 0; + double transm = 0; + enum sphin_btdf_transmissivity_type transmissivity_type = + SPHIN_BTDF_TRANSMISSIVITY_NONE__; + struct sphin_btdf* btdf = NULL; + struct sphin_spectral_property* transmissivity_spectrum = NULL; + struct ssf_fresnel* fresnel = NULL; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != primitive); + ASSERT(NULL != transmissivity); + + res = primitive_get_btdf(sphor, primitive, &btdf); + if (RES_OK != res) { goto error; } + + if (NULL == btdf) { + res = RES_BAD_ARG; + goto error; + } + + res = sphin_btdf_get_transmissivity_type(btdf, &transmissivity_type); + if (RES_OK != res) { goto error; } + + switch (transmissivity_type) { + case SPHIN_BTDF_TRANSMISSIVITY_TABULATED: + + res = sphin_btdf_get_transmissivity_value(btdf, &transmissivity_spectrum); + if (RES_OK != res) { goto error; } + + res = sphin_spectral_property_interpolate_at_wavelength + (transmissivity_spectrum, wavelength, + SPHIN_INTERPOLATION_LINEAR, &transm); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_create + (sphor->allocator, &ssf_fresnel_constant, &fresnel); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_constant_setup(fresnel, 1. - transm); + if (RES_OK != res) { goto error; } + + break; + + case SPHIN_BTDF_TRANSMISSIVITY_FRESNEL_DIELECTRIC: + + /* Get refractive index of the medium in the incident side */ + res = primitive_get_in_refraction_index + (sphor, primitive, wavelength, &eta_i, &k_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, wavelength, &eta_t, &k_t); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_create + (sphor->allocator, &ssf_fresnel_dielectric_dielectric, &fresnel); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_dielectric_dielectric_setup(fresnel, eta_i, eta_t); + if (RES_OK != res) { goto error; } + + break; + + case SPHIN_BRDF_REFLECTIVITY_FRESNEL_DIELECTRIC_CONDUCTOR: + + /* Get refractive index of the medium in the incident side */ + res = primitive_get_in_refraction_index + (sphor, primitive, wavelength, &eta_i, &k_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, wavelength, &eta_t, &k_t); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_create + (sphor->allocator, &ssf_fresnel_dielectric_conductor, &fresnel); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_dielectric_conductor_setup(fresnel, eta_i, eta_t, k_t); + if (RES_OK != res) { goto error; } + + break; + + case SPHIN_BRDF_REFLECTIVITY_NONE__: + res = RES_BAD_ARG; + goto error; + + default: + FATAL("Unreachable code\n"); + break; + } + + *transmissivity = 1 - ssf_fresnel_eval(fresnel, cos_theta); + +exit: + if (NULL != fresnel) { + ssf_fresnel_ref_put(fresnel); + } + return res; +error: + goto exit; +} + +res_T trace_ray (const struct sphor* sphor, const struct ray* ray, diff --git a/src/sphor_interface.h b/src/sphor_interface.h @@ -120,6 +120,12 @@ primitive_get_brdf struct sphin_brdf** out_brdf); extern LOCAL_SYM res_T +primitive_get_btdf + (struct sphor* sphor, + const struct primitive* primitive, + struct sphin_btdf** out_btdf); + +extern LOCAL_SYM res_T primitive_get_in_refraction_index (struct sphor* sphor, struct primitive* primitive, @@ -136,6 +142,22 @@ primitive_get_out_refraction_index double* n_imag); extern LOCAL_SYM res_T +primitive_get_reflectivity + (struct sphor* sphor, + struct primitive* primitive, + double wavelength, + double cos_theta, + double* reflectivity); + +extern LOCAL_SYM res_T +primitive_get_transmissivity + (struct sphor* sphor, + struct primitive* primitive, + double wavelength, + double cos_theta, + double* transmissivity); + +extern LOCAL_SYM res_T trace_ray (const struct sphor* sphor, const struct ray* ray, diff --git a/src/sphor_ran_bsdf.c b/src/sphor_ran_bsdf.c @@ -0,0 +1,315 @@ +/* Copyright (C) 2024-2026 Centre National de la Recherche Scientifique + * Copyright (C) 2024-2026 Clermont Auvergne INP + * Copyright (C) 2024-2026 INSA Lyon + * Copyright (C) 2024-2026 Institut Mines Télécom Albi-Carmaux + * Copyright (C) 2024-2026 Institut National Polytechnique de Toulouse + * Copyright (C) 2024-2026 |Méso|Star> (contact@meso-star.com) + * Copyright (C) 2024-2026 PhotonLyX (info@photonlyx.com) + * Copyright (C) 2024-2026 Université de Lorraine + * Copyright (C) 2024-2026 Université Paul Sabatier + * Copyright (C) 2024-2026 Université Toulouse - Jean Jaurès + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include "sphor_c.h" +#include "sphor_ran_bsdf.h" +#include "sphor_interface.h" + +#include <rsys/rsys.h> +#include <star/ssf.h> + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static res_T +setup_bsdf_reflection + (struct sphor* sphor, + struct ray* ray, + struct intersection* intersection, + double wavelength, + double reflectivity, + struct ssf_bsdf** out_bsdf) +{ + struct primitive* primitive = &intersection->position.primitive; + struct sphin_brdf* brdf = NULL; + struct ssf_bsdf* bsdf = NULL; + struct ssf_fresnel* fresnel = NULL; + enum sphin_brdf_direction_distribution brdf_dir = SPHIN_BRDF_DIRECTION_NONE__; + + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != ray); + ASSERT(NULL != intersection); + + (void)ray; + (void)wavelength; + + primitive = (struct primitive*)&intersection->position.primitive; + + res = primitive_get_brdf(sphor, primitive, &brdf); + if (RES_OK != res) { goto error; } + if (NULL == brdf) { res = RES_BAD_ARG; goto error; } + + res = sphin_brdf_get_direction_distribution(brdf, &brdf_dir); + if (RES_OK != res) { goto error; } + + switch (brdf_dir) { + case SPHIN_BRDF_DIRECTION_LAMBERT: + + res = ssf_bsdf_create + (sphor->allocator, &ssf_lambertian_reflection, &bsdf); + if (RES_OK != res) { goto error; } + + res = ssf_lambertian_reflection_setup(bsdf, reflectivity); + if (RES_OK != res) { goto error; } + + break; + + case SPHIN_BRDF_DIRECTION_SPECULAR: + res = ssf_bsdf_create + (sphor->allocator, &ssf_specular_reflection, &bsdf); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_create + (sphor->allocator, &ssf_fresnel_constant, &fresnel); + if (RES_OK != res) { goto error; } + + res = ssf_fresnel_constant_setup(fresnel, reflectivity); + if (RES_OK != res) { goto error; } + + res = ssf_specular_reflection_setup(bsdf, fresnel); + if (RES_OK != res) { goto error; } + + break; + + case SPHIN_BRDF_DIRECTION_NONE__: + res = RES_BAD_ARG; + goto error; + + default: + FATAL("Unreachable code\n"); + break; + } + +exit: + if (NULL != fresnel) { SSF(fresnel_ref_put(fresnel)); } + *out_bsdf = bsdf; +error: + return res; + if (NULL != bsdf) { SSF(bsdf_ref_put(bsdf)); } + goto exit; +} + +static res_T +setup_bsdf_transmission + (struct sphor* sphor, + struct ray* ray, + struct intersection* intersection, + double wavelength, + double transmissivity, + struct ssf_bsdf** out_bsdf) +{ + struct primitive* primitive = &intersection->position.primitive; + struct sphin_btdf* btdf = NULL; + struct ssf_bsdf* bsdf = NULL; + struct ssf_fresnel* fresnel = NULL; + enum sphin_btdf_direction_distribution btdf_dir = SPHIN_BTDF_DIRECTION_NONE__; + + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != ray); + ASSERT(NULL != intersection); + + (void)ray; + (void)transmissivity; + (void)wavelength; + + primitive = (struct primitive*)&intersection->position.primitive; + + res = primitive_get_btdf(sphor, primitive, &btdf); + if (RES_OK != res) { goto error; } + if (NULL == btdf) { res = RES_BAD_ARG; goto error; } + + res = sphin_btdf_get_direction_distribution(btdf, &btdf_dir); + if (RES_OK != res) { goto error; } + + switch (btdf_dir) { + case SPHIN_BTDF_DIRECTION_LAMBERT: + break; + + case SPHIN_BTDF_DIRECTION_KEEP_CURRENT_DIR: + break; + + case SPHIN_BTDF_DIRECTION_SNELL_DIELECTRIC: + break; + + case SPHIN_BTDF_DIRECTION_NONE__: + res = RES_BAD_ARG; + goto error; + + default: + FATAL("Unreachable code\n"); + break; + } + +exit: + if (NULL != fresnel) { SSF(fresnel_ref_put(fresnel)); } + *out_bsdf = bsdf; + return res; +error: + if (NULL != bsdf) { SSF(bsdf_ref_put(bsdf)); } + goto exit; +} + +static res_T +sample_interaction_type + (struct sphor* sphor, + struct ssp_rng* rng, + double reflectivity, + double transmissivity, + enum intersection_ray_interaction_type* interaction_type) +{ + double s = 0; /* Random number */ + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != rng); + + (void)sphor; + + if (reflectivity + transmissivity > 1) { + res = RES_BAD_ARG; goto error; } + + if (reflectivity < 0) { + res = RES_BAD_ARG; goto error; } + + if (transmissivity < 0) { + res = RES_BAD_ARG; goto error; } + + s = ssp_rng_canonical(rng); + + if (s >= 0 && s <= reflectivity) { + *interaction_type = INTERSECTION_RAY_INTERACTION_REFLECTION; + } else if (s > reflectivity && s <= reflectivity + transmissivity) { + *interaction_type = INTERSECTION_RAY_INTERACTION_TRANSMISSION; + } else if (s > reflectivity + transmissivity && s <= 1) { + *interaction_type = INTERSECTION_RAY_INTERACTION_ABSORPTION; + } else { res = RES_BAD_ARG; goto error; } + +exit: + return res; +error: + goto exit; +} +/******************************************************************************* + * Local functions + ******************************************************************************/ +res_T +sample_interface_ray_interaction + (struct sphor* sphor, + struct ssp_rng* rng, + struct ray* ray, + struct intersection* intersection, + double wavelength, + double dir[3], + enum intersection_ray_interaction_type* interaction_type) +{ + double invert_normal = 0; + double normal[3] = {0}; + double wo[3] = {0}; + double reflectivity = 0; + double transmissivity = 0; + double sample[3] = {0}; + struct primitive* primitive = NULL; + struct sphin_brdf* brdf = NULL; + struct sphin_btdf* btdf = NULL; + struct ssf_bsdf* bsdf = NULL; + res_T res = RES_OK; + + ASSERT(NULL != sphor); + ASSERT(NULL != rng); + ASSERT(NULL != ray); + ASSERT(NULL != intersection); + + primitive = &intersection->position.primitive; + + /* Ensure that the ray direction and the geometry normal are compatible with + * the star-sf library. */ + d3_minus(wo, ray->direction); + d3_set(normal, primitive->normal); + /* Ensure normal and direction point to the same hemisphere */ + invert_normal = (d3_dot(normal, wo) > 0 ) ? 1. : -1.; + d3_normalize(normal, d3_muld(normal, normal, invert_normal)); + + /* Verify if one of the surfaces defined in the intersected interface + * side has a BRDF defined */ + res = primitive_get_brdf(sphor, primitive, &brdf); + if (RES_OK != res) { goto error; } + + /* Verify if one of the surfaces defined in the intersected interface + * side has a BTDF defined */ + res = primitive_get_btdf(sphor, primitive, &btdf); + if (RES_OK != res) { goto error; } + + if (NULL == brdf && NULL == btdf) { /* Nothing to do */ + transmissivity = 1; + *interaction_type = INTERSECTION_RAY_INTERACTION_TRANSMISSION; + d3_set(dir, ray->direction); + } else { + /* Setup reflectivity */ + if (NULL != brdf) { + res = primitive_get_reflectivity + (sphor, primitive, wavelength, + d3_dot(normal, ray->direction), &reflectivity); + if (RES_OK != res) { goto error; } + } + + /* Setup transmissivity */ + if (NULL != btdf) { + res = primitive_get_transmissivity + (sphor, primitive, wavelength, + d3_dot(normal, ray->direction), &transmissivity); + if (RES_OK != res) { goto error; } + } + + /* Sample interaction type */ + res = sample_interaction_type + (sphor, rng, reflectivity, transmissivity, interaction_type); + if (RES_OK != res) { goto error; } + + if (*interaction_type == INTERSECTION_RAY_INTERACTION_ABSORPTION) { + goto exit; /* Nothing else to do */ + } + + if (*interaction_type == INTERSECTION_RAY_INTERACTION_REFLECTION) { + res = setup_bsdf_reflection + (sphor, ray, intersection, wavelength, reflectivity, &bsdf); + } + if (*interaction_type == INTERSECTION_RAY_INTERACTION_TRANSMISSION) { + res = setup_bsdf_transmission + (sphor, ray, intersection, wavelength, transmissivity, &bsdf); + } + ssf_bsdf_sample(bsdf, rng, wo, normal, sample, NULL, NULL); + + d3_set(dir, sample); + } + +exit: + if (NULL != bsdf) { SSF(bsdf_ref_put(bsdf)); } + return res; +error: + goto exit; +} diff --git a/src/sphor_ran_bsdf.h b/src/sphor_ran_bsdf.h @@ -0,0 +1,52 @@ +/* Copyright (C) 2024-2026 Centre National de la Recherche Scientifique + * Copyright (C) 2024-2026 Clermont Auvergne INP + * Copyright (C) 2024-2026 INSA Lyon + * Copyright (C) 2024-2026 Institut Mines Télécom Albi-Carmaux + * Copyright (C) 2024-2026 Institut National Polytechnique de Toulouse + * Copyright (C) 2024-2026 |Méso|Star> (contact@meso-star.com) + * Copyright (C) 2024-2026 PhotonLyX (info@photonlyx.com) + * Copyright (C) 2024-2026 Université de Lorraine + * Copyright (C) 2024-2026 Université Paul Sabatier + * Copyright (C) 2024-2026 Université Toulouse - Jean Jaurès + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ +#ifndef SPHOR_RAN_BSDF_H +#define SPHOR_RAN_BSDF_H + +#include <rsys/rsys.h> + +/* Forward declarations */ +struct sphor; +struct intersection; +struct ssp_rng; +struct ssf_bsdf; + +enum intersection_ray_interaction_type { + INTERSECTION_RAY_INTERACTION_REFLECTION, /* 0 */ + INTERSECTION_RAY_INTERACTION_TRANSMISSION, /* 1 */ + INTERSECTION_RAY_INTERACTION_ABSORPTION, /* 2 */ + INTERSECTION_RAY_INTERACTION_NONE__ +}; + +extern LOCAL_SYM res_T +sample_interface_ray_interaction + (struct sphor* sphor, + struct ssp_rng* rng, + struct ray* ray, + struct intersection* intersection, + double wavelength, + double dir[3], + enum intersection_ray_interaction_type* interaction_type); + +#endif /* SPHOR_RAN_BSDF_H */ diff --git a/src/test_sphor_MVREA_analytical1.c b/src/test_sphor_MVREA_analytical1.c @@ -115,14 +115,28 @@ write_spectral_file(const char* filename) CHK(fclose(fp) == 0); } +static void +write_reflectivity_file(const char* filename) +{ + FILE* fp = NULL; + + fp = fopen(filename, "w"); + CHK(NULL != fp); + + fprintf(fp, "450 0.5\n"); + fprintf(fp, "550 0.5\n"); + CHK(fclose(fp) == 0); +} static void write_input_file(FILE* fp) { const char* spectral_filename = "spec.txt"; + const char* reflectivity_filename = "rho.txt"; write_cube(); write_spectral_file(spectral_filename); + write_reflectivity_file(reflectivity_filename); fprintf(fp, "surface: \"source\"\n"); fprintf(fp, "\tgeometry: BACK cube_front.stl\n"); @@ -134,7 +148,8 @@ write_input_file(FILE* fp) fprintf(fp, "\tgeometry: BACK cube.stl\n"); fprintf(fp, "\tgeometry: BACK cube_front.stl\n"); fprintf(fp, "\tgeometry: BACK cube_back.stl\n"); - fprintf(fp, "\tprop_rad: \"sigma_a\" SCATTERER\n"); + fprintf(fp, "\tprop_rad: \"sigma_a\" \n"); + fprintf(fp, "\tscatterer:\n"); fprintf(fp, "\t\tconcentration: 1 mol/m^3\n"); fprintf(fp, "\t\tcross_sections:\n"); fprintf(fp, "\t\t\tabs_cross_sec: spec.txt nm m^2/mol\n"); @@ -142,7 +157,7 @@ write_input_file(FILE* fp) fprintf(fp, "\t\tresponse_function: 1\n"); fprintf(fp, "surface: \"mirror\"\n"); fprintf(fp, "\tgeometry: BACK cube_back.stl\n"); - fprintf(fp, "\tbrdf: SPECULAR 0.5\n"); + fprintf(fp, "\tbrdf: SPECULAR rho.txt\n"); fprintf(fp, "\tsensor:\n"); fprintf(fp, "\t\tresponse_function: 1\n"); CHK(fflush(fp) == 0);