commit 45429811b5b7c03fc09eaf307671877e8be78d85
parent fc4dfff17129da72a2d3408571d6d18418c74141
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date: Fri, 6 Mar 2026 18:01:41 +0100
Merge branch 'feature_brdf_btdf'
Diffstat:
11 files changed, 1975 insertions(+), 231 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -3,11 +3,13 @@
.config
.test
file.txt
+rho.txt
spec.txt
*.pc
*.so
*.swp
*.stl
+*.dat
tags
input
output
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
@@ -173,7 +174,10 @@ uninstall:
################################################################################
TEST_SRC =\
src/test_sphor_lib.c\
- src/test_sphor_MVREA_analytical1.c
+ src/test_sphor_MVREA_analytical1.c\
+ src/test_sphor_MVREA_btdf_keep_current_dir.c\
+ src/test_sphor_MVREA_btdf_snell_dielectric.c\
+ src/test_sphor_MVREA_btdf_lambertian.c
TEST_OBJ = $(TEST_SRC:.c=.o)
TEST_DEP = $(TEST_SRC:.c=.d)
TEST_TGT = $(TEST_SRC:.c=.t)
@@ -221,10 +225,14 @@ $(TEST_OBJ): config.mk sphor-local.pc
test_sphor_lib\
test_sphor_MVREA_analytical1\
+test_sphor_MVREA_btdf_keep_current_dir\
+test_sphor_MVREA_btdf_snell_dielectric\
+test_sphor_MVREA_btdf_lambertian\
: config.mk sphor-local.pc $(LIBNAME)
$(CC) $(CFLAGS_TEST) -o $@ src/$@.o $(LDFLAGS_TEST)
clean_test:
rm -f $(TEST_DEP) $(TEST_OBJ) $(TEST_TGT) \
- input output cube.stl cube_back.stl cube_front.stl spec.txt
+ input output *.stl *.dat \
+ 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,703 @@
+/* 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>
+
+/*******************************************************************************
+ * Internal BSDFs
+ ******************************************************************************/
+
+/* BTDF LAMBERTIAN */
+struct lambertian_transmission {
+ double transmissivity;
+};
+
+static double
+lambertian_transmission_eval
+ (void* data,
+ const double wo[3],
+ const double N[3],
+ const double wi[3])
+{
+ struct lambertian_transmission* btdf = data;
+
+ ASSERT(NULL != data);
+ ASSERT(NULL != N);
+ ASSERT(NULL != wi);
+ ASSERT(d3_is_normalized(N) && d3_is_normalized(wi));
+ ASSERT(d3_dot(wi, N) < 0 && d3_dot(wo, N) > 0);
+ (void)wo, (void)N, (void)wi;
+
+ return btdf->transmissivity/ PI;
+}
+
+static double
+lambertian_transmission_sample
+ (void* data,
+ struct ssp_rng* rng,
+ const double wo[3],
+ const double N[3],
+ double wi[3],
+ int* type,
+ double* pdf)
+{
+ double sample[3];
+
+ ASSERT(NULL != data);
+ ASSERT(NULL != rng);
+ ASSERT(NULL != N);
+ ASSERT(NULL != wi);
+ ASSERT(d3_is_normalized(wo) && d3_is_normalized(N) && d3_dot(wo, N) > 0);
+ (void)wo;
+
+ ssp_ran_hemisphere_cos(rng, N, sample, pdf);
+ d3_muld(wi, sample, -1);
+ if (type) *type = SSF_TRANSMISSION | SSF_DIFFUSE;
+ return ((struct lambertian_transmission*)data)->transmissivity;
+}
+
+static double
+lambertian_transmission_pdf
+ (void* data,
+ const double wo[3],
+ const double N[3],
+ const double wi[3])
+{
+ double cos_wi_N;
+ ASSERT(NULL != data);
+ ASSERT(NULL != N);
+ ASSERT(NULL != wi);
+ ASSERT(d3_is_normalized(N) && d3_is_normalized(wi));
+ (void)data, (void)wo;
+
+ cos_wi_N = d3_dot(wi, N);
+ return cos_wi_N < 0.0 ? -cos_wi_N / PI : 0.0;
+}
+
+static res_T
+lambertian_transmission_setup
+ (struct ssf_bsdf* bsdf,
+ const double transmissivity)
+{
+ void* ptr = NULL;
+
+ ASSERT(NULL != bsdf);
+ ASSERT(0 <= transmissivity && transmissivity <=1);
+
+ ssf_bsdf_get_data(bsdf, &ptr);
+ struct lambertian_transmission* data = ptr;
+ data->transmissivity = transmissivity;
+
+ return RES_OK;
+}
+
+const struct ssf_bsdf_type lambertian_transmission = {
+ NULL,
+ NULL,
+ lambertian_transmission_sample,
+ lambertian_transmission_eval,
+ lambertian_transmission_pdf,
+ sizeof(struct lambertian_transmission),
+ ALIGNOF(struct lambertian_transmission)
+};
+
+/* BTDF KEEP_CURRENT_DIR */
+struct keep_current_dir_transmission {
+ struct ssf_fresnel* fresnel;
+};
+
+static void
+keep_current_dir_transmission_release
+ (void* data)
+{
+ struct keep_current_dir_transmission* btdf = data;
+ ASSERT(data);
+ if (btdf->fresnel) SSF(fresnel_ref_put(btdf->fresnel));
+}
+
+static double
+keep_current_dir_transmission_sample
+ (void* data,
+ struct ssp_rng* rng,
+ const double wo[3],
+ const double N[3],
+ double wi[3],
+ int* type,
+ double* pdf)
+{
+ struct keep_current_dir_transmission* btdf = data;
+ double cos_wo_N;
+
+ ASSERT(NULL != data);
+ ASSERT(NULL != rng);
+ ASSERT(NULL != N);
+ ASSERT(NULL != wi);
+ ASSERT(d3_is_normalized(wo) && d3_is_normalized(N) && d3_dot(wo, N) > 0);
+ (void)rng;
+
+ /* In ssf convention, wo points outwards the surface */
+ d3_minus(wi, wo);
+ if (pdf) *pdf = INF;
+ if (type) *type = SSF_TRANSMISSION;
+
+ cos_wo_N = d3_dot(wo, N);
+ return 1 - ssf_fresnel_eval(btdf->fresnel, cos_wo_N);
+}
+
+static double
+keep_current_dir_transmission_eval
+ (void* data,
+ const double wo[3],
+ const double N[3],
+ const double wi[3])
+{
+ (void)data, (void)wi, (void)N, (void)wo;
+ return 0.0;
+}
+
+static double
+keep_current_dir_transmission_pdf
+ (void* data,
+ const double wo[3],
+ const double N[3],
+ const double wi[3])
+{
+ (void)data, (void)wi, (void)N, (void)wo;
+ return 0.0;
+}
+
+static res_T
+keep_current_dir_transmission_setup
+ (struct ssf_bsdf* bsdf,
+ struct ssf_fresnel* fresnel)
+{
+ void* ptr = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != bsdf);
+ ASSERT(NULL != fresnel);
+
+ ssf_bsdf_get_data(bsdf, &ptr);
+ struct keep_current_dir_transmission* data = ptr;
+
+ if (data->fresnel != fresnel) {
+ if (NULL != data->fresnel) {
+ res = ssf_fresnel_ref_put(data->fresnel);
+ if (RES_OK != res) { return res; }
+ }
+ res = ssf_fresnel_ref_get(fresnel);
+ if (RES_OK != res) { return res; }
+
+ data->fresnel = fresnel;
+ }
+
+ return RES_OK;
+}
+
+const struct ssf_bsdf_type keep_current_dir_transmission = {
+ NULL,
+ keep_current_dir_transmission_release,
+ keep_current_dir_transmission_sample,
+ keep_current_dir_transmission_eval,
+ keep_current_dir_transmission_pdf,
+ sizeof(struct keep_current_dir_transmission),
+ ALIGNOF(struct keep_current_dir_transmission)
+};
+
+/* BTDF SNELL DIELECTRIC */
+struct snell_dielectric_transmission {
+ struct ssf_fresnel* fresnel;
+ double eta_i; /* Refractive index of the incoming medium */
+ double eta_t; /* Refractive index of the transmissive medium */
+};
+
+/* Refract the vect V wrt the normal N using the relative refractive index eta.
+ * Eta is the refraction index of the outside medium (where N points into)
+ * devided by the refraction index of the inside medium. By convention N and V
+ * points on the same side of the surface. */
+static INLINE void
+refract(double res[3], const double V[3], const double N[3], const double eta)
+{
+ double tmp0[3];
+ double tmp1[3];
+ double cos_theta_i;
+ double cos_theta_t;
+ double sin2_theta_i;
+ double sin2_theta_t;
+
+ ASSERT(res && V && N);
+ ASSERT(d3_is_normalized(V) && d3_is_normalized(N));
+ cos_theta_i = d3_dot(V, N);
+ sin2_theta_i = MMAX(0, 1.0 - cos_theta_i*cos_theta_i);
+ sin2_theta_t = eta * eta * sin2_theta_i;
+ cos_theta_t = sqrt(1 - sin2_theta_t);
+
+ d3_muld(tmp0, V, eta);
+ d3_muld(tmp1, N, eta * cos_theta_i - cos_theta_t);
+ d3_sub(res, tmp1, tmp0);
+}
+
+static void
+snell_dielectric_transmission_release
+ (void* data)
+{
+ struct keep_current_dir_transmission* btdf = data;
+ ASSERT(data);
+ if (btdf->fresnel) SSF(fresnel_ref_put(btdf->fresnel));
+}
+
+static double
+snell_dielectric_transmission_sample
+ (void* data,
+ struct ssp_rng* rng,
+ const double wo[3],
+ const double N[3],
+ double wi[3],
+ int* type,
+ double* pdf)
+{
+ struct snell_dielectric_transmission* btdf = data;
+ double wt[3];
+ double cos_wo_N;
+ double eta; /* Ratio of eta_i / eta_t */
+
+ ASSERT(NULL != btdf);
+ ASSERT(NULL != data);
+ ASSERT(NULL != rng);
+ ASSERT(NULL != N);
+ ASSERT(NULL != wi);
+ ASSERT(d3_is_normalized(wo) && d3_is_normalized(N));
+ ASSERT(d3_dot(wo, N) > -1.e-6);
+ (void)rng;
+
+ eta = btdf->eta_i / btdf->eta_t;
+ refract(wt, wo, N, eta);
+
+ cos_wo_N = MMAX(0.0, d3_dot(wo, N));
+
+ if(pdf) *pdf = INF;
+ d3_set(wi, wt);
+ if(type) *type = SSF_SPECULAR | SSF_TRANSMISSION;
+ return 1 - ssf_fresnel_eval(btdf->fresnel, cos_wo_N);
+}
+
+static double
+snell_dielectric_transmission_eval
+ (void* bsdf,
+ const double wo[3],
+ const double N[3],
+ const double wi[3])
+{
+ (void)bsdf, (void)wo, (void)N, (void)wi;
+ return 0.0;
+}
+
+static double
+snell_dielectric_transmission_pdf
+ (void* bsdf,
+ const double wo[3],
+ const double N[3],
+ const double wi[3])
+{
+ (void)bsdf, (void)wo, (void)N, (void)wi;
+ return 0.0;
+}
+
+static res_T
+snell_dielectric_transmission_setup
+ (struct ssf_bsdf* bsdf,
+ struct ssf_fresnel* fresnel,
+ double eta_i,
+ double eta_t)
+{
+ void* ptr = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != bsdf);
+
+ ssf_bsdf_get_data(bsdf, &ptr);
+ struct snell_dielectric_transmission* data = ptr;
+
+ if (data->fresnel != fresnel) {
+ if (NULL != data->fresnel) {
+ res = ssf_fresnel_ref_put(data->fresnel);
+ if (RES_OK != res) { return res; }
+ }
+ res = ssf_fresnel_ref_get(fresnel);
+ if (RES_OK != res) { return res; }
+
+ data->fresnel = fresnel;
+ }
+ data->eta_i = eta_i;
+ data->eta_t = eta_t;
+
+ return RES_OK;
+}
+
+const struct ssf_bsdf_type snell_dielectric_transmission = {
+ NULL,
+ snell_dielectric_transmission_release,
+ snell_dielectric_transmission_sample,
+ snell_dielectric_transmission_eval,
+ snell_dielectric_transmission_pdf,
+ sizeof(struct snell_dielectric_transmission),
+ ALIGNOF(struct snell_dielectric_transmission)
+};
+
+/*******************************************************************************
+ * 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)
+{
+ double eta_i = 0, k_i = 0;
+ double eta_t = 0, k_t = 0;
+ 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)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:
+ res = ssf_bsdf_create
+ (sphor->allocator, &lambertian_transmission, &bsdf);
+ if (RES_OK != res) { goto error; }
+
+ res = lambertian_transmission_setup(bsdf, transmissivity);
+ if (RES_OK != res) { goto error; }
+
+ break;
+
+ case SPHIN_BTDF_DIRECTION_KEEP_CURRENT_DIR:
+ res = ssf_bsdf_create
+ (sphor->allocator, &keep_current_dir_transmission, &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, 1 - transmissivity);
+ if (RES_OK != res) { goto error; }
+
+ res = keep_current_dir_transmission_setup(bsdf, fresnel);
+ if (RES_OK != res) { goto error; }
+
+ break;
+
+ case SPHIN_BTDF_DIRECTION_SNELL_DIELECTRIC:
+
+ res = ssf_bsdf_create
+ (sphor->allocator, &snell_dielectric_transmission, &bsdf);
+ if (RES_OK != res) { goto error; }
+
+ /* 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);
+
+ res = ssf_fresnel_create
+ (sphor->allocator, &ssf_fresnel_constant, &fresnel);
+ if (RES_OK != res) { goto error; }
+
+ res = ssf_fresnel_constant_setup(fresnel, 1 - transmissivity);
+ if (RES_OK != res) { goto error; }
+
+ res = snell_dielectric_transmission_setup(bsdf, fresnel, eta_i, eta_t);
+ if (RES_OK != res) { goto error; }
+
+ 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);
diff --git a/src/test_sphor_MVREA_btdf_keep_current_dir.c b/src/test_sphor_MVREA_btdf_keep_current_dir.c
@@ -0,0 +1,272 @@
+/* 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/>. */
+#define _POSIX_C_SOURCE 200112L
+
+#include "sphor.h"
+
+#include <rsys/cstr.h>
+#include <rsys/logger.h>
+#include <rsys/str.h>
+#include <rsys/text_reader.h>
+
+#include <math.h>
+#include <stdio.h>
+
+static void
+write_emitting_surface()
+{
+ FILE* fp = fopen("emitting_surface.stl", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp,
+ "solid emitting_surface\n"
+ " facet normal 0.5 0.8660254037844386 0\n"
+ " outer loop\n"
+ " vertex -1.4433756729740643 -1.5 -1\n"
+ " vertex -0.5773502691896257 -2 1\n"
+ " vertex -0.5773502691896257 -2 -1\n"
+ " endloop\n"
+ " endfacet\n"
+ " facet normal 0.5 0.8660254037844386 0\n"
+ " outer loop\n"
+ " vertex -0.5773502691896257 -2 1\n"
+ " vertex -1.4433756729740643 -1.5 -1\n"
+ " vertex -1.4433756729740643 -1.5 1\n"
+ " endloop\n"
+ " endfacet\n"
+ "endsolid emitting_surface\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_btdf_surface()
+{
+ FILE* fp = fopen("btdf_surface.stl", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp,
+ "solid btdf_surface\n"
+ " facet normal 1 0 0\n"
+ " outer loop\n"
+ " vertex 0 -1 -1\n"
+ " vertex 0 1 1\n"
+ " vertex 0 -1 1\n"
+ " endloop\n"
+ " endfacet\n"
+ " facet normal 1 0 0\n"
+ " outer loop\n"
+ " vertex 0 -1 -1\n"
+ " vertex 0 1 -1\n"
+ " vertex 0 1 1\n"
+ " endloop\n"
+ " endfacet\n"
+ "endsolid btdf_surface\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_absorbing_surface()
+{
+ FILE* fp = fopen("absorbing_surface.stl", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp,
+ "solid absorbing_surface\n"
+ " facet normal -0.5 -0.8660254037844386 0\n"
+ " outer loop\n"
+ " vertex 1.4433756729740643 1.5 -1\n"
+ " vertex 0.5773502691896257 2 1\n"
+ " vertex 0.5773502691896257 2 -1\n"
+ " endloop\n"
+ " endfacet\n"
+ " facet normal -0.5 -0.8660254037844386 0\n"
+ " outer loop\n"
+ " vertex 0.5773502691896257 2 1\n"
+ " vertex 1.4433756729740643 1.5 -1\n"
+ " vertex 1.4433756729740643 1.5 1\n"
+ " endloop\n"
+ " endfacet\n"
+ "endsolid absorbing_surface\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_spectrum()
+{
+ FILE* fp = fopen("emission_spectrum.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 1\n");
+ fprintf(fp, "401 1\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_btdf()
+{
+ FILE* fp = fopen("btdf.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 1\n");
+ fprintf(fp, "401 1\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_bsdf_null()
+{
+ FILE* fp = fopen("bsdf_null.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 0\n");
+ fprintf(fp, "401 0\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_input_file(FILE* fp)
+{
+ write_emitting_surface();
+ write_btdf_surface();
+ write_absorbing_surface();
+ write_spectrum();
+ write_btdf();
+ write_bsdf_null();
+
+ fprintf(fp, "surface: \"source\"\n");
+ fprintf(fp, "\tgeometry: FRONT emitting_surface.stl\n");
+ fprintf(fp, "\tsource:\n");
+ fprintf(fp, "\tflux_density: 1 umol/m^2/s emission_spectrum.dat nm nm^-1\n");
+ fprintf(fp, "\t\tdirection: COLLIM NORMAL\n");
+ fprintf(fp, "\n");
+
+ fprintf(fp, "surface: \"transmitter\"\n");
+ fprintf(fp, "\tgeometry: BACK btdf_surface.stl\n");
+ fprintf(fp, "\tbtdf: KEEP_CURRENT_DIR btdf.dat\n");
+ fprintf(fp, "\n");
+
+ fprintf(fp, "surface: \"absorber\"\n");
+ fprintf(fp, "\tgeometry: FRONT absorbing_surface.stl\n");
+ fprintf(fp, "\tbrdf: SPECULAR bsdf_null.dat\n");
+ fprintf(fp, "\tsensor:\n");
+ fprintf(fp, "\t\tresponse_function: 1\n");
+
+ CHK(fflush(fp) == 0);
+}
+
+int main(int argc, char** argv)
+{
+ struct sphor* sphor = NULL;
+ struct sphor_create_args args = SPHOR_CREATE_ARGS_DEFAULT;
+ struct txtrdr* txtrdr = NULL;
+ struct str line;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ char input_filename[] = "input";
+ char output_filename[] = "output";
+ double avg = 0;
+ double std = 0;
+ double nsamples = 0;
+
+ FILE* fp = NULL;
+
+ fp = fopen(input_filename, "w+");
+ CHK(NULL != fp);
+
+ write_input_file(fp);
+ CHK(fclose(fp) == 0);
+
+ args.input_filename = input_filename;
+ args.output_filename = output_filename;
+ args.force = 1;
+
+ (void)argc;
+ (void)argv;
+
+ /* Test sphor_create with NULL logger and allocator */
+ CHK(sphor_create(&args, &sphor) == RES_OK);
+
+ /* Run sphor with the config file */
+ CHK(sphor_run(sphor) == RES_OK);
+
+ CHK(sphor_ref_put(sphor) == RES_OK);
+
+ fp = fopen(output_filename, "r");
+ CHK(NULL != fp);
+
+ str_init(NULL, &line);
+ CHK(txtrdr_stream(NULL, fp, output_filename, '#', &txtrdr) == RES_OK);
+
+ nsamples = (double)args.samples;
+
+ /* Verify first level: total scene MVREA */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "MVREA") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 0, 2 * std / sqrt(nsamples)));
+
+ /* Verify first level: total scene losses */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "LOSSES") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 1, 2 * std / sqrt(nsamples)));
+
+ /* Verify second level: per surface losses */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "LOSSES") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(strcmp(token, "absorber") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 1, 2 * std / sqrt(nsamples)));
+
+ CHK(fclose(fp) == 0);
+
+ str_release(&line);
+ txtrdr_ref_put(txtrdr);
+
+ CHK(mem_allocated_size() == 0);
+
+ return 0;
+}
diff --git a/src/test_sphor_MVREA_btdf_lambertian.c b/src/test_sphor_MVREA_btdf_lambertian.c
@@ -0,0 +1,292 @@
+/* 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/>. */
+#define _POSIX_C_SOURCE 200112L
+
+#include "sphor.h"
+
+#include <rsys/cstr.h>
+#include <rsys/logger.h>
+#include <rsys/str.h>
+#include <rsys/text_reader.h>
+
+#include <math.h>
+#include <stdio.h>
+
+static void
+write_emitting_surface()
+{
+ FILE* fp = fopen("emitting_surface.stl", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp,
+ "solid emitting_surface\n"
+ " facet normal 0.5 0.8660254037844386 0\n"
+ " outer loop\n"
+ " vertex -1.4433756729740643 -1.5 -1\n"
+ " vertex -0.5773502691896257 -2 1\n"
+ " vertex -0.5773502691896257 -2 -1\n"
+ " endloop\n"
+ " endfacet\n"
+ " facet normal 0.5 0.8660254037844386 0\n"
+ " outer loop\n"
+ " vertex -0.5773502691896257 -2 1\n"
+ " vertex -1.4433756729740643 -1.5 -1\n"
+ " vertex -1.4433756729740643 -1.5 1\n"
+ " endloop\n"
+ " endfacet\n"
+ "endsolid emitting_surface\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_btdf_surface()
+{
+ FILE* fp = fopen("btdf_surface.stl", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp,
+ "solid btdf_surface\n"
+ " facet normal 1 0 0\n"
+ " outer loop\n"
+ " vertex 0 -1 -1\n"
+ " vertex 0 1 1\n"
+ " vertex 0 -1 1\n"
+ " endloop\n"
+ " endfacet\n"
+ " facet normal 1 0 0\n"
+ " outer loop\n"
+ " vertex 0 -1 -1\n"
+ " vertex 0 1 -1\n"
+ " vertex 0 1 1\n"
+ " endloop\n"
+ " endfacet\n"
+ "endsolid btdf_surface\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_absorbing_surface()
+{
+ FILE* fp = fopen("absorbing_surface.stl", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp,
+ "solid absorbing_surface\n"
+ " facet normal 1 0 0\n"
+ " outer loop\n"
+ " vertex 1 -1 -1\n"
+ " vertex 1 1 1\n"
+ " vertex 1 -1 1\n"
+ " endloop\n"
+ " endfacet\n"
+ " facet normal 1 0 0\n"
+ " outer loop\n"
+ " vertex 1 -1 -1\n"
+ " vertex 1 1 -1\n"
+ " vertex 1 1 1\n"
+ " endloop\n"
+ " endfacet\n"
+ "endsolid absorbing_surface\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_spectrum()
+{
+ FILE* fp = fopen("emission_spectrum.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 1\n");
+ fprintf(fp, "401 1\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_btdf()
+{
+ FILE* fp = fopen("btdf.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 1\n");
+ fprintf(fp, "401 1\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_bsdf_null()
+{
+ FILE* fp = fopen("bsdf_null.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 0\n");
+ fprintf(fp, "401 0\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_input_file(FILE* fp)
+{
+ write_emitting_surface();
+ write_btdf_surface();
+ write_absorbing_surface();
+ write_spectrum();
+ write_btdf();
+ write_bsdf_null();
+
+ fprintf(fp, "surface: \"source\"\n");
+ fprintf(fp, "\tgeometry: FRONT emitting_surface.stl\n");
+ fprintf(fp, "\tsource:\n");
+ fprintf(fp, "\tflux_density: 2 umol/m^2/s emission_spectrum.dat nm nm^-1\n");
+ fprintf(fp, "\t\tdirection: COLLIM NORMAL\n");
+ fprintf(fp, "\n");
+
+ fprintf(fp, "surface: \"transmitter\"\n");
+ fprintf(fp, "\tgeometry: BACK btdf_surface.stl\n");
+ fprintf(fp, "\tbtdf: LAMBERT btdf.dat\n");
+ fprintf(fp, "\n");
+
+ fprintf(fp, "surface: \"absorber\"\n");
+ fprintf(fp, "\tgeometry: BACK absorbing_surface.stl\n");
+ fprintf(fp, "\tbrdf: SPECULAR bsdf_null.dat\n");
+ fprintf(fp, "\tsensor:\n");
+ fprintf(fp, "\t\tresponse_function: 1\n");
+
+ CHK(fflush(fp) == 0);
+}
+
+int main(int argc, char** argv)
+{
+ struct sphor* sphor = NULL;
+ struct sphor_create_args args = SPHOR_CREATE_ARGS_DEFAULT;
+ struct txtrdr* txtrdr = NULL;
+ struct str line;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ char input_filename[] = "input";
+ char output_filename[] = "output";
+ double avg = 0;
+ double std = 0;
+ double nsamples = 0;
+
+ FILE* fp = NULL;
+
+ fp = fopen(input_filename, "w+");
+ CHK(NULL != fp);
+
+ write_input_file(fp);
+ CHK(fclose(fp) == 0);
+
+ args.input_filename = input_filename;
+ args.output_filename = output_filename;
+ args.force = 1;
+
+ (void)argc;
+ (void)argv;
+
+ /* Test sphor_create with NULL logger and allocator */
+ CHK(sphor_create(&args, &sphor) == RES_OK);
+
+ /* Run sphor with the config file */
+ CHK(sphor_run(sphor) == RES_OK);
+
+ CHK(sphor_ref_put(sphor) == RES_OK);
+
+ fp = fopen(output_filename, "r");
+ CHK(NULL != fp);
+
+ str_init(NULL, &line);
+ CHK(txtrdr_stream(NULL, fp, output_filename, '#', &txtrdr) == RES_OK);
+
+ nsamples = (double)args.samples;
+
+ /* Verify first level: total scene MVREA */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "MVREA") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 0, 2 * std / sqrt(nsamples)));
+
+ /* Verify first level: total scene losses */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "LOSSES") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare against the analytical view factor.
+ *
+ * This computes the exact configuration (view) factor F12 between two
+ * coaxial, parallel square surfaces of side length 2 separated by a
+ * distance l = 1.
+ *
+ * The expression evaluated is:
+ *
+ * F12 = (1 / A1) int_{A1} int_{A2} ( l^2 / (pi r^4) ) dA2 dA1
+ *
+ * where:
+ * - A1 and A2 are the two square surfaces,
+ * - l is the axial separation,
+ * - r = ||p2 - p1|| is the distance between differential elements,
+ * - the integrand results from the general view factor formula
+ * (cos theta_1 cos theta_2) / (pi r^2),
+ * using cos theta_1 = cos theta_2 = l / r for parallel planes.
+ *
+ * The result is the density of diffuse power leaving the btdf surface
+ * that reaches the absorbing surface, weighted by the surface area ratio of
+ * the emitting surface and the absorbing surface.
+ */
+ CHK(eq_eps(avg, 0.4152532835771472, 2 * std / sqrt(nsamples)));
+
+ /* Verify second level: per surface losses */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "LOSSES") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(strcmp(token, "absorber") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ CHK(eq_eps(avg, 0.4152532835771472, 2 * std / sqrt(nsamples)));
+
+ CHK(fclose(fp) == 0);
+
+ str_release(&line);
+ txtrdr_ref_put(txtrdr);
+
+ CHK(mem_allocated_size() == 0);
+
+ return 0;
+}
diff --git a/src/test_sphor_MVREA_btdf_snell_dielectric.c b/src/test_sphor_MVREA_btdf_snell_dielectric.c
@@ -0,0 +1,309 @@
+/* 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/>. */
+#define _POSIX_C_SOURCE 200112L
+
+#include "sphor.h"
+
+#include <rsys/cstr.h>
+#include <rsys/logger.h>
+#include <rsys/str.h>
+#include <rsys/text_reader.h>
+
+#include <math.h>
+#include <stdio.h>
+
+static void
+write_emitting_surface()
+{
+ FILE* fp = fopen("emitting_surface.stl", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp,
+ "solid emitting_surface\n"
+ " facet normal 0.5 0.8660254037844386 0\n"
+ " outer loop\n"
+ " vertex -1.4433756729740643 -1.5 -1\n"
+ " vertex -0.5773502691896257 -2 1\n"
+ " vertex -0.5773502691896257 -2 -1\n"
+ " endloop\n"
+ " endfacet\n"
+ " facet normal 0.5 0.8660254037844386 0\n"
+ " outer loop\n"
+ " vertex -0.5773502691896257 -2 1\n"
+ " vertex -1.4433756729740643 -1.5 -1\n"
+ " vertex -1.4433756729740643 -1.5 1\n"
+ " endloop\n"
+ " endfacet\n"
+ "endsolid emitting_surface\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_btdf_surface()
+{
+ FILE* fp = fopen("btdf_surface.stl", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp,
+ "solid btdf_surface\n"
+ " facet normal 1 0 0\n"
+ " outer loop\n"
+ " vertex 0 -1 -1\n"
+ " vertex 0 1 1\n"
+ " vertex 0 -1 1\n"
+ " endloop\n"
+ " endfacet\n"
+ " facet normal 1 0 0\n"
+ " outer loop\n"
+ " vertex 0 -1 -1\n"
+ " vertex 0 1 -1\n"
+ " vertex 0 1 1\n"
+ " endloop\n"
+ " endfacet\n"
+ "endsolid btdf_surface\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_absorbing_surface()
+{
+ FILE* fp = fopen("absorbing_surface.stl", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp,
+ "solid absorbing_surface\n"
+ "facet normal 0 0 0\n"
+ "outer loop\n"
+ "vertex 1.7320508075688772 0 -1\n"
+ "vertex 1.7320508075688772 0 1\n"
+ "vertex 0.8660254037844386 1.5 1\n"
+ "endloop\n"
+ "endfacet\n"
+ "facet normal 0 0 0\n"
+ "outer loop\n"
+ "vertex 0.8660254037844386 1.5 1\n"
+ "vertex 0.8660254037844386 1.5 -1\n"
+ "vertex 1.7320508075688772 0 -1\n"
+ "endloop\n"
+ "endfacet\n"
+ "endsolid absorbing_surface\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_spectrum()
+{
+ FILE* fp = fopen("emission_spectrum.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 1\n");
+ fprintf(fp, "401 1\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_btdf()
+{
+ FILE* fp = fopen("btdf.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 1\n");
+ fprintf(fp, "401 1\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_bsdf_null()
+{
+ FILE* fp = fopen("bsdf_null.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 0\n");
+ fprintf(fp, "401 0\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_n_back()
+{
+ FILE* fp = fopen("n_back.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 1\n");
+ fprintf(fp, "401 1\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_n_front()
+{
+ FILE* fp = fopen("n_front.dat", "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "400 1.7320508075688772\n");
+ fprintf(fp, "401 1.7320508075688772\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_input_file(FILE* fp)
+{
+ write_emitting_surface();
+ write_btdf_surface();
+ write_absorbing_surface();
+ write_spectrum();
+ write_btdf();
+ write_bsdf_null();
+ write_n_front();
+ write_n_back();
+
+ fprintf(fp, "surface: \"source\"\n");
+ fprintf(fp, "\tgeometry: FRONT emitting_surface.stl\n");
+ fprintf(fp, "\tsource:\n");
+ fprintf(fp, "\tflux_density: 1.7320508075688772 umol/m^2/s");
+ fprintf(fp, " emission_spectrum.dat nm nm^-1\n");
+ fprintf(fp, "\t\tdirection: COLLIM NORMAL\n");
+ fprintf(fp, "\n");
+
+ fprintf(fp, "surface: \"transmitter\"\n");
+ fprintf(fp, "\tgeometry: BACK btdf_surface.stl\n");
+ fprintf(fp, "\tbtdf: SNELL_DIELECTRIC btdf.dat\n");
+ fprintf(fp, "\n");
+
+ fprintf(fp, "surface: \"absorber\"\n");
+ fprintf(fp, "\tgeometry: FRONT absorbing_surface.stl\n");
+ fprintf(fp, "\tbrdf: SPECULAR bsdf_null.dat\n");
+ fprintf(fp, "\tsensor:\n");
+ fprintf(fp, "\t\tresponse_function: 1\n");
+
+ fprintf(fp, "volume: \"n_front\"\n");
+ fprintf(fp, "\tgeometry: FRONT btdf_surface.stl\n");
+ fprintf(fp, "\trefractive_index:\n");
+ fprintf(fp, "\t\tn_real: n_front.dat\n");
+
+ fprintf(fp, "volume: \"n_back\"\n");
+ fprintf(fp, "\tgeometry: BACK btdf_surface.stl\n");
+ fprintf(fp, "\trefractive_index:\n");
+ fprintf(fp, "\t\tn_real: n_back.dat\n");
+
+ CHK(fflush(fp) == 0);
+}
+
+int main(int argc, char** argv)
+{
+ struct sphor* sphor = NULL;
+ struct sphor_create_args args = SPHOR_CREATE_ARGS_DEFAULT;
+ struct txtrdr* txtrdr = NULL;
+ struct str line;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ char input_filename[] = "input";
+ char output_filename[] = "output";
+ double avg = 0;
+ double std = 0;
+ double nsamples = 0;
+
+ FILE* fp = NULL;
+
+ fp = fopen(input_filename, "w+");
+ CHK(NULL != fp);
+
+ write_input_file(fp);
+ CHK(fclose(fp) == 0);
+
+ args.input_filename = input_filename;
+ args.output_filename = output_filename;
+ args.force = 1;
+
+ (void)argc;
+ (void)argv;
+
+ /* Test sphor_create with NULL logger and allocator */
+ CHK(sphor_create(&args, &sphor) == RES_OK);
+
+ /* Run sphor with the config file */
+ CHK(sphor_run(sphor) == RES_OK);
+
+ CHK(sphor_ref_put(sphor) == RES_OK);
+
+ fp = fopen(output_filename, "r");
+ CHK(NULL != fp);
+
+ str_init(NULL, &line);
+ CHK(txtrdr_stream(NULL, fp, output_filename, '#', &txtrdr) == RES_OK);
+
+ nsamples = (double)args.samples;
+
+ /* Verify first level: total scene MVREA */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "MVREA") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 0, 2 * std / sqrt(nsamples)));
+
+ /* Verify first level: total scene losses */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "LOSSES") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 1, 2 * std / sqrt(nsamples)));
+
+ /* Verify second level: per surface losses */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "LOSSES") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(strcmp(token, "absorber") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 1, 2 * std / sqrt(nsamples)));
+
+ CHK(fclose(fp) == 0);
+
+ str_release(&line);
+ txtrdr_ref_put(txtrdr);
+
+ CHK(mem_allocated_size() == 0);
+
+ return 0;
+}
