star-phor

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

commit a51410c83da1e972f0c26bb9f7316e566d0bc822
parent 4e881934e3b6593e2ff6219bc34133f78548dfbb
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date:   Fri,  9 Jan 2026 16:21:40 +0100

Refactor ray–interface interactions to rely exclusively on ssf

This commit unifies the handling of all ray-interface interactions by
sampling them entirely with the star scattering functions (ssf) library.

Previously, ssf was used only for dielectric-dielectric interfaces,
while interfaces involving user-defined sphin_brdf s implemented their
own logic for reflectivity tests, interaction type selection, and
outgoing direction sampling directly in star-phor. This resulted in
duplicated implementations of Lambertian and specular reflection models,
since ssf is also capable of doing these computations.

The code has been refactored to construct appropriate ssf BSDFs  for
sphin_brdf based interfaces as well. All interactions are now sampled
through ssf_bsdf_sample, regardless of the interface type
(dielectric-dielectric, specular, Lambertian).

Diffstat:
MMakefile | 1-
Msrc/sphor_compute_mvrea.c | 184+++++++++++++++++++++++++++++++++++--------------------------------------------
Dsrc/sphor_ran_brdf.c | 95-------------------------------------------------------------------------------
Dsrc/sphor_ran_brdf.h | 42------------------------------------------
4 files changed, 82 insertions(+), 240 deletions(-)

diff --git a/Makefile b/Makefile @@ -42,7 +42,6 @@ SRC_LIB =\ src/sphor_compute_mvrea.c\ src/sphor_config.c\ src/sphor_interface.c\ - src/sphor_ran_brdf.c\ src/sphor_ran_geometry.c\ src/sphor_ran_source.c\ src/sphor_sources.c diff --git a/src/sphor_compute_mvrea.c b/src/sphor_compute_mvrea.c @@ -28,7 +28,6 @@ #include "sphor_compute_mvrea.h" #include "sphor_interface.h" #include "sphor_ran_source.h" -#include "sphor_ran_brdf.h" #include <star/sphin.h> #include <star/s3d.h> @@ -426,25 +425,24 @@ volume_get_from_primitive } static res_T -interface_sample_ray_interaction_type_brdf - (struct sphor* sphor, - struct ssp_rng* rng, - const struct ray* ray, - const struct intersection* intersection, - enum ray_interface_interaction_type* interaction_type) +setup_bsdf_from_sphin_brdf + (struct sphor* sphor, + struct ray* ray, + struct intersection* intersection, + struct ssf_fresnel** out_fresnel, + struct ssf_bsdf** out_bsdf) { - enum sphin_brdf_type brdf_type = SPHIN_BRDF_NONE__; - struct sphin_brdf* brdf = NULL; - struct primitive* primitive = NULL; double reflectivity = 0; - double s = 0; /* random number */ + 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 != rng); ASSERT(NULL != ray); ASSERT(NULL != intersection); - ASSERT(NULL != interaction_type); (void)ray; @@ -452,7 +450,7 @@ interface_sample_ray_interaction_type_brdf res = primitive_get_brdf(sphor, primitive, &brdf); if (RES_OK != res){ goto error; } - if (NULL != brdf) { res = RES_BAD_ARG; 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; } @@ -461,70 +459,35 @@ interface_sample_ray_interaction_type_brdf 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; } - break; - case SPHIN_BRDF_NONE__: - res = RES_BAD_ARG; - goto error; - default: - FATAL("Unreachable code\n"); - break; - } - /* Bernoulli test if the photon was absorbed or reflected */ - s = ssp_rng_canonical(rng); - if (s < reflectivity) { - *interaction_type = RAY_INTERFACE_INTERACTION_REFLECTION; - } - else { *interaction_type = RAY_INTERFACE_INTERACTION_ABSORPTION; } - -exit: - return res; -error: - goto exit; -} - -static res_T -sample_interface_ray_interaction_from_brdf -(struct sphor* sphor, - struct ssp_rng* rng, - struct ray* ray, - struct intersection* intersection, - double dir[3], - enum ray_interface_interaction_type* interaction_type) -{ - struct sphin_brdf* brdf = NULL; - struct primitive* primitive = NULL; - res_T res = RES_OK; - - ASSERT(NULL != sphor); - ASSERT(NULL != rng); - ASSERT(NULL != ray); - ASSERT(NULL != intersection); - primitive = &intersection->position.primitive; + res = ssf_bsdf_create + (sphor->allocator, &ssf_specular_reflection, &bsdf); + if (RES_OK != res){ goto error; } - res = primitive_get_brdf(sphor, primitive, &brdf); - if (RES_OK != res){ goto error; } - if (NULL != brdf) { res = RES_BAD_ARG; goto error; } + res = ssf_fresnel_create + (sphor->allocator, &ssf_fresnel_constant, &fresnel); + if (RES_OK != res){ goto error; } - res = interface_sample_ray_interaction_type_brdf - (sphor, rng, ray, intersection, interaction_type); - if (RES_OK != res){ goto error; } + res = ssf_fresnel_constant_setup(fresnel, reflectivity); + if (RES_OK != res){ goto error; } - switch (*interaction_type) { - case RAY_INTERFACE_INTERACTION_ABSORPTION: - break; - case RAY_INTERFACE_INTERACTION_REFLECTION: - res = ran_brdf_reflection_direction - (brdf, rng, primitive, ray->direction, dir); - break; - case RAY_INTERFACE_INTERACTION_TRANSMISSION: - d3_set(dir, ray->direction); + res = ssf_specular_reflection_setup(bsdf, fresnel); + if (RES_OK != res){ goto error; } break; - case RAY_INTERFACE_INTERACTION_NONE__: + + case SPHIN_BRDF_NONE__: res = RES_BAD_ARG; goto error; default: @@ -533,37 +496,31 @@ sample_interface_ray_interaction_from_brdf } 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 -sample_interface_ray_interaction_from_refractive_index +setup_bsdf_from_refractive_indices (struct sphor* sphor, - struct ssp_rng* rng, struct ray* ray, struct intersection* intersection, - double dir[3], - enum ray_interface_interaction_type* interaction_type) + struct ssf_bsdf** out_bsdf) { - double invert_normal = 1.; - double normal[3] = {0}; double n_real_i = 0, n_imag_i = 0; double n_real_t = 0, n_imag_t = 0; - double pdf = 0; - double sample[3] = {0}; - int flag = 0; struct primitive* primitive = &intersection->position.primitive; struct ssf_bsdf* bsdf = NULL; res_T res = RES_OK; ASSERT(NULL != sphor); - ASSERT(NULL != rng); ASSERT(NULL != ray); ASSERT(NULL != intersection); - ASSERT(NULL != dir); - ASSERT(NULL != interaction_type); /* Get refractive index of the medium in the incident side */ res = primitive_get_in_refraction_index @@ -575,11 +532,6 @@ sample_interface_ray_interaction_from_refractive_index (sphor, primitive, ray->wavelength, &n_real_t, &n_imag_t); if (RES_OK != res){ goto error; } - d3_set(normal, primitive->normal); - /* Ensure normal and incoming direction point to the same hemisphere */ - invert_normal = (d3_dot(normal, ray->direction) > 0 ) ?1. : -1.; - d3_normalize(normal, d3_muld(normal, normal, invert_normal)); - /* Sample an interaction type and a new direction using ssf */ res = ssf_bsdf_create (sphor->allocator, &ssf_specular_dielectric_dielectric_interface, &bsdf); @@ -588,20 +540,12 @@ sample_interface_ray_interaction_from_refractive_index res = ssf_specular_dielectric_dielectric_interface_setup (bsdf, n_real_i, n_real_t); if (RES_OK != res){ goto error; } - ssf_bsdf_sample(bsdf, rng, ray->direction, normal, sample, &flag, &pdf); - - 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 != bsdf) { SSF(bsdf_ref_put(bsdf)); } + *out_bsdf = bsdf; return res; error: + if (NULL != bsdf) { SSF(bsdf_ref_put(bsdf)); } goto exit; } @@ -614,8 +558,16 @@ sample_interface_ray_interaction double dir[3], enum ray_interface_interaction_type* interaction_type) { - struct sphin_brdf* brdf = NULL; + double invert_normal = 0; + double normal[3] = {0}; + double wo[3] = {0}; + double pdf = 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); @@ -625,6 +577,14 @@ sample_interface_ray_interaction primitive = &intersection->position.primitive; + /* Ensure that the ray direction and the geometry normal are compatible with + * th 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); @@ -634,20 +594,39 @@ sample_interface_ray_interaction /* 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 = sample_interface_ray_interaction_from_brdf - (sphor, rng, ray, intersection, dir, interaction_type); + 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 = sample_interface_ray_interaction_from_refractive_index - (sphor, rng, ray, intersection, dir, interaction_type); + else { + res = setup_bsdf_from_refractive_indices + (sphor, ray, intersection, &bsdf); if (RES_OK != res){ goto error; } } + + ssf_bsdf_sample(bsdf, rng, wo, normal, sample, &flag, &pdf); + + if (flag & SSF_TRANSMISSION) { + if (NULL != brdf) { + *interaction_type = RAY_INTERFACE_INTERACTION_ABSORPTION; + } + else { + *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; @@ -870,6 +849,7 @@ compute_MVREA_realization /* Update ray with the intersection information */ res = ray_update(&ray, &intersection, dir); if (RES_OK != res){ goto error; } + } exit: diff --git a/src/sphor_ran_brdf.c b/src/sphor_ran_brdf.c @@ -1,95 +0,0 @@ -/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique - * Copyright (C) 2024-2025 Clermont Auvergne INP - * Copyright (C) 2024-2025 INSA Lyon - * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux - * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse - * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com) - * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com) - * Copyright (C) 2024-2025 Université de Lorraine - * Copyright (C) 2024-2025 Université Paul Sabatier - * Copyright (C) 2024-2025 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_interface.h" -#include "sphor_ran_brdf.h" -#include "sphor_sources.h" - -#include <star/sphin.h> -#include <star/ssp.h> -#include <star/suniq.h> -#include <rsys/dynamic_array_double.h> - -/******************************************************************************* - * Helper functions - ******************************************************************************/ - -/******************************************************************************* - * Local functions - ******************************************************************************/ -res_T -ran_brdf_reflection_direction - (struct sphin_brdf* brdf, - struct ssp_rng* rng, - struct primitive* primitive, - double initial_dir[3], - double final_dir[3]) -{ - enum sphin_brdf_type brdf_type = SPHIN_BRDF_NONE__; - double normal[3] = {0}; - double new_dir[3] = {0}; - double dot = 0; - double invert_normal = 1.; - res_T res = RES_OK; - - d3_set(normal, primitive->normal); - - /* Ensure the geometry normal points into the same hemisphere as the incoming - * direction. If the dot product between the normal and the incoming (initial) - * direction is positive, then the normal points in the same direction as the - * incoming ray, which is incorrect for reflection. Invert the normal in that - * case. */ - invert_normal = (d3_dot(normal, initial_dir) > 0 ) ?-1. : 1.; - - res = sphin_brdf_get_type(brdf, &brdf_type); - if (RES_OK != res){ goto error; } - - d3_normalize(normal, d3_muld(normal, normal, invert_normal)); - - switch (brdf_type) { - case SPHIN_BRDF_LAMBERT: - /* Sample a cosine-weighted direction over the hemisphere oriented - * by the surface normal. */ - ssp_ran_hemisphere_cos(rng, normal, new_dir, NULL); - break; - case SPHIN_BRDF_SPECULAR: - /* Use the formula to compute the reflection: - * reflected = initial_dir − 2 (initial_dir dot normal) normal */ - dot = d3_dot(normal, initial_dir); - d3_add(new_dir, initial_dir, d3_muld(new_dir, normal, -2*dot)); - break; - case SPHIN_BRDF_NONE__: - res = RES_BAD_ARG; - goto error; - default: - FATAL("Unreachable code\n"); - break; - } - d3_normalize(final_dir, new_dir); - -exit: - return res; -error: - goto exit; -} diff --git a/src/sphor_ran_brdf.h b/src/sphor_ran_brdf.h @@ -1,42 +0,0 @@ -/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique - * Copyright (C) 2024-2025 Clermont Auvergne INP - * Copyright (C) 2024-2025 INSA Lyon - * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux - * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse - * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com) - * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com) - * Copyright (C) 2024-2025 Université de Lorraine - * Copyright (C) 2024-2025 Université Paul Sabatier - * Copyright (C) 2024-2025 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_BRDF_H -#define SPHOR_RAN_BRDF_H - -#include <rsys/rsys.h> - -/* Forward declarations */ -struct sphin_brdf; -struct ssp_rng; -struct s3d_primitive; - -extern LOCAL_SYM res_T -ran_brdf_reflection_direction - (struct sphin_brdf* brdf, - struct ssp_rng* rng, - struct primitive* primitive, - double initial_dir[3], - double final_dir[3]); - -#endif /* SPHOR_RAN_BRDF */