commit fbe60e6000694c1cf95afed6521719449e74736f
parent 46e0ec243da62bcdefee4e080fdc107808433ec0
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date: Wed, 9 Jul 2025 18:00:04 +0200
Implement first Monte Carlo algorithm
This commit implements a Monte Carlo algorithm that computes the
fraction of photons emitted in a scene that are absorbed in a volume
acting as a sensor. This type of algorithm is particularly useful for
computing the Mean Volumetric Rate of Energy Absorption (MVREA). This
quantity is related to absorptivity via its product with a factor that
accounts for the source power and the specific illuminated surface.
Several useful abstractions are also introduced, such as functions for
sampling path departure directions according to a directional
distribution, and handling reflections on BRDF surfaces.
The realization loop is parallelized. As a result, this commit
introduces a dependency on the OpenMPI library.
Diffstat:
16 files changed, 951 insertions(+), 113 deletions(-)
diff --git a/Makefile b/Makefile
@@ -38,7 +38,12 @@ all: executable library tests
################################################################################
SRC_LIB =\
src/sphor.c\
+ 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
OBJ_LIB = $(SRC_LIB:.c=.o)
DEP_LIB = $(SRC_LIB:.c=.d)
@@ -58,6 +63,7 @@ $(LIBNAME): $(OBJ_LIB)
$(CC) $(CFLAGS_LIB) -o $@ $(OBJ_LIB) $(LDFLAGS_LIB)
.config: config.mk
+ $(PKG_CONFIG) --atleast-version $(MPI_VERSION) $(MPI_PC)
$(PKG_CONFIG) --atleast-version $(RSYS_VERSION) rsys
$(PKG_CONFIG) --atleast-version $(S3D_VERSION) s3d
$(PKG_CONFIG) --atleast-version $(SPHIN_VERSION) sphin
diff --git a/README.md b/README.md
@@ -12,6 +12,7 @@ Solver for radiative transfer in photoreactors.
- star-3d
- star-sp
- star-uniq
+- Open MPI
## Installation
diff --git a/config.mk b/config.mk
@@ -18,18 +18,22 @@ LD = ld
OBJCOPY = objcopy
PKG_CONFIG = pkg-config
RANLIB = ranlib
+MPI_PC = ompi
################################################################################
# Dependencies
################################################################################
+MPI_VERSION = 2
RSYS_VERSION=0.14
S3D_VERSION=0.10
SPHIN_VERSION=0.0
SSP_VERSION=0.14
SUNIQ_VERSION=0.0
-INCS = $$($(PKG_CONFIG) --cflags rsys s3d sphin star-sp suniq)
-LIBS = $$($(PKG_CONFIG) --libs rsys s3d sphin star-sp suniq)
+INCS = $$($(PKG_CONFIG) --cflags $(MPI_PC) rsys s3d sphin star-sp suniq)\
+ -fopenmp -lm
+LIBS = $$($(PKG_CONFIG) --libs $(MPI_PC) rsys s3d sphin star-sp suniq)\
+ -fopenmp -lm
################################################################################
# Compilation options
diff --git a/src/sphor.c b/src/sphor.c
@@ -24,6 +24,7 @@
#include "sphor.h"
#include "sphor_c.h"
+#include "sphor_compute_mvrea.h"
#include <star/sphin.h>
#include <star/s3d.h>
@@ -147,11 +148,14 @@ sphor_ref_put(struct sphor* sphor)
res_T
sphor_run(struct sphor* sphor)
{
+ double estim = 0;
+ double std = 0;
res_T res = RES_OK;
if (NULL == sphor) { goto error; }
- /* TODO Run Monte Carlo algorithm */
+ res = sphor_compute_MVREA(sphor, &estim, &std);
+ if (RES_OK != res) { goto error; }
exit:
return res;
diff --git a/src/sphor_compute_mvrea.c b/src/sphor_compute_mvrea.c
@@ -0,0 +1,334 @@
+/* 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.h"
+#include "sphor_c.h"
+#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>
+#include <star/ssp.h>
+#include <omp.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+compute_MVREA_realization
+(struct sphor* sphor,
+ struct ssp_rng* rng,
+ double* weight/* TODO make this a list in sphor, since we can have multiple
+ sensors. Each sensor should have its own weight */)
+{
+ double ka = 0; /* Absorption coefficient of current volume */
+ double normal[3] = {0}; /* Geometry normal */
+ double dir[3] = {0}; /* Current direction in path sampling */
+ double pos[3] = {0}; /* Current position in path sampling */
+ double sample[3] = {0};
+ double reflectivity = 0; /* Surface reflectivity */
+ double response_function = 0;
+ double s = 0; /* Random number */
+ enum sphin_brdf_type brdf_type = SPHIN_BRDF_NONE__;
+ enum sphin_side side_id = SPHIN_SIDE_NONE__;
+ /* Position and direction using when calling s3d functions */
+ float r[3] = {0}; /* Position */
+ float w[3] = {0}; /* Direction */
+ /* Allowed range for rays during path sampling */
+ float range[2] = {0.f, FLT_MAX};
+ float st[2] = {0}; /* Parametric coordinates of a point in a s3d_primitive */
+ int i = 0; /* Iterator */
+ size_t triangle_id;
+ struct interface interface = INTERFACE_NULL;
+ struct s3d_attrib attrib;
+ struct s3d_hit hit = S3D_HIT_NULL;
+ struct s3d_primitive prim = S3D_PRIMITIVE_NULL;
+ struct source_view* source_view = NULL;
+ struct sphin_brdf* brdf = NULL;
+ struct sphin_sensor_volume* sensor_volume = NULL;
+ struct sphin_source_surface* source = NULL;
+ struct sphin_surface* surface = NULL;
+ struct sphin_volume* volume = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != rng);
+
+ /* Sample a random position from a source in the whole scene */
+ res = sample_source_position(sphor, rng, &source_view, &prim, st);
+ if (RES_OK != res){ goto error; }
+ res = s3d_primitive_get_attrib (&prim, S3D_POSITION, st, &attrib);
+ if (RES_OK != res){ goto error; }
+ d3_set_f3(pos, attrib.value);
+
+ /* Retrieve geometry normal of the sampled primitive */
+ res = s3d_primitive_get_attrib (&prim, S3D_GEOMETRY_NORMAL, st, &attrib);
+ if (RES_OK != res){ goto error; }
+ d3_normalize(normal, d3_minus(normal, d3_set_f3(normal, attrib.value)));
+
+ /* TODO Sample a wavelength when spectral information available in sphin */
+
+ res = sphin_config_get_surface
+ (sphor->config, source_view->sphin_id, &surface);
+ if (RES_OK != res){ goto error; }
+ res = sphin_surface_get_source(surface, &source);
+ if (RES_OK != res){ goto error; }
+ if (NULL == source){ res = RES_BAD_ARG; goto error; }
+
+ /* Sample a direction according to the source direction distribution */
+ res = source_surface_sample_direction(source, rng, &prim, dir);
+ if (RES_OK != res){ goto error; }
+
+ /* While no absorption takes place and an interface is found */
+ for(;;) {
+ /* Trace a ray in the scene from the sampled primitive in the sampled
+ * direction and retrive the hit distance */
+ f3_set_d3(r, pos);
+ f3_set_d3(w, dir);
+ res = s3d_scene_view_trace_ray
+ (sphor->scene_view, r, w, range, &prim, &hit);
+ if (RES_OK != res){ goto error; }
+
+ /* If there is no intersection in the sampled direction from the sampled
+ * position, the photon is lost and it counts zero in the weight. */
+ if (S3D_HIT_NONE(&hit)) {
+ *weight = 0;
+ break;
+ }
+
+ /* Retrieve the intercepted triangle and the corresponding interface in the
+ * sphor structure */
+ triangle_id = hit.prim.prim_id;
+ interface = darray_interface_data_get(&sphor->interfaces)[triangle_id];
+
+ /* Retrieve the position of the hit and set the new position */
+ res = s3d_primitive_get_attrib
+ (&hit.prim, S3D_POSITION, st, &attrib);
+ if (RES_OK != res){ goto error; }
+ d3_set_f3(pos, attrib.value);
+ /* Update current primitive */
+ prim = hit.prim;
+
+ /* Retrieve the side of the intercepted interface corresponding to the
+ * income direction of the hit */
+ res = interface_get_side_id(dir, &hit, &side_id);
+ if (RES_OK != res){ goto error; }
+
+ /* Get the absorption coefficient of the current volume to compute an
+ * absorption free length. If the intersected interface does not correspond
+ * to any volume in the config, the absorption coefficient is supposed to be
+ * zero and the absorption free length is set to FLT_MAX */
+ if (INVALID_ID != interface.volumes[side_id]) {
+ /* Retrieve the volume we are currently in */
+ res = sphin_config_get_volume
+ (sphor->config, interface.volumes[side_id], &volume);
+ if (RES_OK != res){ goto error; }
+
+ res = sphin_volume_get_ka(volume, &ka);
+ if (RES_OK != res){ goto error; }
+
+ /* Sample free path length according to the volume absorption coefficent */
+ s = ssp_ran_exp(rng, ka);}
+ else {
+ s = FLT_MAX;
+ }
+
+ /* If the distance between the ray origin and the intersection is bigger
+ * than the absorption free path length, an absorption takes place (check
+ * if the volume is a sensor and return the weight accordingly) */
+ if (s < hit.distance) {
+ /* Check if the volume is a sensor */
+ res = sphin_volume_get_sensor(volume, &sensor_volume);
+ /* If the volume is a sensor, compute the weight using the response
+ * function */
+ if (NULL != sensor_volume) {
+ res = sphin_sensor_volume_get_response_function
+ (sensor_volume, &response_function);
+ if (RES_OK != res){ goto error; }
+ ASSERT(response_function == 1);
+ *weight = response_function;
+ }
+ /* If the volume is not a sensor, the absorbed photon does not count in
+ * the weight */
+ else {
+ *weight = 0;
+ }
+ break;
+ }
+
+ /* Else, the distance between the ray origin and the intersection is smaller
+ * than the absorption free path length. The photon intersects a primitive
+ * in the scene view. */
+
+ /* Use sphin to retrieve the properties of the intersected surface */
+ FOR_EACH(i, 0, interface.surface_count[side_id]) {
+ res = sphin_config_get_surface
+ (sphor->config, interface.surfaces[side_id][i], &surface);
+ if (RES_OK != res){ goto error; }
+ res = sphin_surface_get_brdf(surface, &brdf);
+ if (RES_OK != res){ goto error; }
+ if (NULL != brdf){ break; }
+ }
+
+ /* 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 */
+ if((NULL != brdf) & (interface.surface_count[side_id] > 0)){
+ 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; }
+ 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;
+ }
+
+ /* Test if the photon was absorbed or reflected */
+ s = ssp_rng_canonical(rng);
+ /* An absorption took place, the weight is zero and the path stops */
+ if (s > reflectivity) {
+ *weight = 0;
+ break;
+ }
+
+ /* Else, if a reflection takes place, sample a new direction according to
+ * the photon arriving direction, geometry normak and surface brdf
+ * properties */
+ res = ran_brdf_reflection_direction(brdf, rng, &hit.prim, dir, sample);
+ if (RES_OK != res){ goto error; }
+ d3_set(dir, sample);
+ }
+
+ /* If the interface side does not have a BRDF (i.e., it defines only a
+ * sphin_volume or a sphin_surface without a BRDF), the photon
+ * crosses the interface the path continues */
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+sphor_compute_MVREA
+ (struct sphor* sphor,
+ double* out_estim,
+ double* out_std)
+{
+ res_T res = RES_OK;
+ size_t nthreads, samples;
+ size_t i; /* iterator */
+ struct ssp_rng_proxy *rng_proxy = NULL;
+ struct ssp_rng **rngs = NULL;
+ double* sum_thread = NULL;
+ double* sum2_thread = NULL;
+ double sum = 0, sum2 = 0;
+ double estim = 0, std = 0;
+
+ samples = sphor->samples;
+ nthreads = sphor->nthreads;
+ if (UINT_MAX == nthreads) { /* use all threads available if not in args */
+ nthreads = (size_t)omp_get_num_procs();
+ }
+
+ /* Création du générateur mandataire RNG_MT19937_64 (Mersenne Twister) */
+ res = ssp_rng_proxy_create(NULL, SSP_RNG_MT19937_64, nthreads, &rng_proxy);
+ if (res != RES_OK) { goto error; }
+
+ /* Allocation des generateurs aleatoire pour chaque processus */
+ rngs = mem_calloc(nthreads, sizeof(*rngs));
+ if (NULL == rngs) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ /* Le generateur mandataire attribue un generateur par sequence aleatoire
+ unique que nous stockons dans le tableau rngs (un par processus) */
+ FOR_EACH(i, 0, nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, &rngs[i]);
+ if (RES_OK != res) { goto error; }
+ }
+ omp_set_num_threads((int)nthreads);
+ sum_thread = mem_calloc(nthreads, sizeof(*sum_thread));
+ sum2_thread = mem_calloc(nthreads, sizeof(*sum2_thread));
+ /* Realizations loop */
+#pragma omp parallel for schedule(static)
+ for(i=0; i<samples; i++) {
+ const int ithread = omp_get_thread_num();
+ double w = 0;
+ res_T res_local = RES_OK;
+
+ if (RES_OK != res) continue;
+
+ res_local = compute_MVREA_realization(sphor, rngs[ithread], &w);
+ if (RES_OK != res_local) {
+ res = res_local;
+ } else {
+ sum_thread[ithread] += w;
+ sum2_thread[ithread] += w*w;
+ }
+ }
+ if(res != RES_OK) goto error;
+
+ FOR_EACH(i, 0, nthreads) {
+ sum += sum_thread[i];
+ sum2 += sum2_thread[i];
+ }
+ estim = sum/(double)samples;
+ std = sqrt
+ ((sum2/(double)samples - estim*estim)/((double)samples-1));
+
+ printf("%f +/- %f\n", estim, std);
+
+exit:
+ if(sum_thread) mem_rm(sum_thread);
+ if(sum2_thread) mem_rm(sum2_thread);
+ /* Liberation de l'ensemble des generateurs */
+ if(rng_proxy) ssp_rng_proxy_ref_put(rng_proxy);
+ if(rngs) {
+ FOR_EACH(i, 0, nthreads) {
+ if(rngs[i]) ssp_rng_ref_put(rngs[i]);
+ }
+ mem_rm(rngs);
+ }
+ *out_estim = estim;
+ *out_std = std;
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sphor_compute_mvrea.h b/src/sphor_compute_mvrea.h
@@ -0,0 +1,39 @@
+/* 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_COMPUTE_MVREA_H
+#define SPHOR_COMPUTE_MVREA_H
+
+#include <rsys/rsys.h>
+
+/* Forward declarations */
+struct sphor;
+struct sphor_create_args;
+
+extern LOCAL_SYM res_T
+sphor_compute_MVREA
+ (struct sphor* sphor,
+ double *out_estim,
+ double *out_std);
+
+#endif /* SPHOR_COMPUTE_MVREA_H */
diff --git a/src/sphor_interface.c b/src/sphor_interface.c
@@ -0,0 +1,61 @@
+/* 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 <rsys/double3.h>
+#include <star/s3d.h>
+#include <star/sphin.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+interface_get_side_id
+ (const double dir[3], /* Direction of arrival of the hit into the geometry */
+ const struct s3d_hit* hit,
+ enum sphin_side* side)
+{
+ double normal[3] = {0};
+ res_T res = RES_OK;
+
+ ASSERT(NULL != dir);
+ ASSERT(NULL != side);
+
+ /* hit->normal: unormalized hit normal that uses the left hand convention */
+ d3_normalize(normal, d3_set_f3(normal, hit->normal));
+
+ /* Test the angle between the two vectors dir and normal using the dot
+ * product. If it is negative, the S3D normal points to the hemisphere of the
+ * arriving hit, which corresponds to the BACK side in the sphin convention.
+ * Otherwise, it corresponds to the front side. */
+ *side = (d3_dot(dir, normal) < 0 ) ? SPHIN_SIDE_BACK : SPHIN_SIDE_FRONT;
+
+ return res;
+}
diff --git a/src/sphor_interface.h b/src/sphor_interface.h
@@ -25,7 +25,10 @@
#ifndef SPHOR_INTERFACE_H
#define SPHOR_INTERFACE_H
+#include "sphor.h"
+
#include <rsys/dynamic_array.h>
+#include <star/sphin.h>
/* Maximum number of surfaces a single triangle can be part of (per side).
* We use a constant to avoid dynamic memory allocation.
@@ -34,6 +37,9 @@
#define MAX_SURFACE_COUNT 4
#define INVALID_ID SIZE_MAX
+/* Forward declarations */
+struct s3d_hit;
+
struct interface {
/* Interface side; FRONT and BACK */
size_t volumes[2/*#sides*/];
@@ -51,4 +57,10 @@ static const struct interface INTERFACE_NULL = INTERFACE_NULL__;
#define DARRAY_DATA struct interface
#include <rsys/dynamic_array.h>
+extern LOCAL_SYM res_T
+interface_get_side_id
+ (const double dir[3],
+ const struct s3d_hit* hit,
+ enum sphin_side* side);
+
#endif /* SPHOR_INTERFACE_H */
diff --git a/src/sphor_ran_brdf.c b/src/sphor_ran_brdf.c
@@ -0,0 +1,99 @@
+/* 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_sources.h"
+#include "sphor_ran_brdf.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 s3d_primitive* prim,
+ double initial_dir[3],
+ double final_dir[3])
+{
+ struct s3d_attrib attrib;
+ enum sphin_brdf_type brdf_type = SPHIN_BRDF_NONE__;
+ float st[2] = {0};
+ double normal[3] = {0};
+ double new_dir[3] = {0};
+ double dot = 0;
+ double invert_normal = 1.f;
+ res_T res = RES_OK;
+
+ res = s3d_primitive_get_attrib
+ (prim, S3D_GEOMETRY_NORMAL, st, &attrib);
+ if (RES_OK != res){ goto error; }
+
+ d3_set_f3(normal, attrib.value);
+ /* 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.f : 1.f;
+
+ 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
@@ -0,0 +1,42 @@
+/* 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 s3d_primitive* prim,
+ double initial_dir[3],
+ double final_dir[3]);
+
+#endif /* SPHOR_RAN_BRDF */
diff --git a/src/sphor_ran_geometry.c b/src/sphor_ran_geometry.c
@@ -0,0 +1,102 @@
+/* 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_ran_geometry.h"
+
+#include <star/s3d.h>
+#include <star/ssp.h>
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+/* Compute the rotation of the vector 'vec' by an angle of 'theta' radians in
+ * the axis 'axis' using the Rodrigues' rotation formula:
+ * https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula */
+res_T
+rotate_vector
+ (double* vrot,
+ double* vec,
+ double* axis,
+ double theta)
+{
+ res_T res = RES_OK;
+
+ double cos_theta = cos(theta);
+ double sin_theta = sin(theta);
+ double dot;
+ double cross[3];
+ double tmp1[3];
+ double tmp2[3];
+ double tmp3[3];
+
+ dot = d3_dot(axis, vec);
+ d3_cross(cross, axis, vec);
+
+ /* Compute the three terms of the Rodrigues' rotation formula */
+ d3_muld(tmp1, axis, (1-cos_theta) * dot);
+ d3_muld(tmp2, cross, sin_theta);
+ d3_muld(tmp3, vec, cos_theta);
+
+ /* Sum the three terms of the formula */
+ d3_add(vrot, tmp1, tmp2);
+ d3_add(vrot, vrot, tmp3);
+
+ return res;
+}
+
+res_T
+sample_surface_position
+ (struct sphor* sphor,
+ struct ssp_rng* rng,
+ struct s3d_scene_view* view,
+ struct s3d_primitive* primitive,
+ float st[2])
+{
+ float u = 0;
+ float v = 0;
+ float w = 0;
+ float uv[2];
+ res_T res = RES_OK;
+
+ ASSERT(NULL != rng);
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != view);
+
+ (void)sphor;
+
+ u = (float)ssp_rng_canonical(rng);
+ v = (float)ssp_rng_canonical(rng);
+ w = (float)ssp_rng_canonical(rng);
+
+ res = s3d_scene_view_sample(view, u, v, w, primitive, uv);
+ if (RES_OK != res) { goto error; }
+
+ st[0] = uv[0];
+ st[1] = uv[1];
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sphor_ran_geometry.h b/src/sphor_ran_geometry.h
@@ -0,0 +1,51 @@
+/* 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_GEOMETRY_H
+#define SPHOR_RAN_GEOMETRY_H
+
+#include <rsys/rsys.h>
+
+/* Forward declarations */
+struct sphor;
+struct ssp_rng;
+struct s3d_scene_view;
+struct s3d_primitive;
+
+extern LOCAL_SYM res_T
+rotate_vector
+ (double* vrot,
+ double* vec,
+ double* axis,
+ double theta);
+
+extern LOCAL_SYM res_T
+sample_surface_position
+ (struct sphor* sphor,
+ struct ssp_rng* rng,
+ struct s3d_scene_view* view,
+ struct s3d_primitive* primitive,
+ float st[2]);
+
+#endif /* SPHOR_RAN_GEOMETRY_H */
diff --git a/src/sphor_ran_source.c b/src/sphor_ran_source.c
@@ -0,0 +1,133 @@
+/* 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_c.h"
+#include "sphor_sources.h"
+#include "sphor_ran_source.h"
+#include "sphor_ran_geometry.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
+source_surface_sample_direction
+ (struct sphin_source_surface* source,
+ struct ssp_rng* rng,
+ struct s3d_primitive* prim,
+ double dir[3])
+{
+ struct s3d_attrib attrib;
+ struct sphin_source_surface_direction_distribution src_dir_dist
+ = SPHIN_SOURCE_SURFACE_DIRECTION_DISTRIBUTION_NULL;
+ double normal[3] = {0};
+ double sample[3] = {0};
+ float st[2] = {0};
+ res_T res = RES_OK;
+
+ /* Sample a direction according to the source direction distribution */
+ res = sphin_source_surface_get_direction_distribution(source, &src_dir_dist);
+ if (RES_OK != res){ goto error; }
+
+ res = s3d_primitive_get_attrib
+ (prim, S3D_GEOMETRY_NORMAL, st, &attrib);
+ if (RES_OK != res){ goto error; }
+ d3_minus(normal, d3_set_f3(normal, attrib.value));
+
+ switch (src_dir_dist.type) {
+ case SPHIN_SOURCE_DIRECTION_COLLIM:
+ d3_set(sample, normal);
+ break;
+ case SPHIN_SOURCE_DIRECTION_COS_POW_N:
+ /* TODO Sample direction with collimation degree */
+ break;
+ case SPHIN_SOURCE_DIRECTION_ISOTROPIC:
+ ssp_ran_hemisphere_cos(rng, normal, sample, NULL);
+ break;
+ case SPHIN_SOURCE_DIRECTION_NONE__:
+ res = RES_BAD_ARG;
+ goto error;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ d3_normalize(dir, sample);
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+sample_source_view
+ (struct sphor* sphor,
+ struct ssp_rng* rng,
+ struct source_view** source_view)
+{
+ struct source_view* source_views = NULL;
+ size_t isource_view = 0;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != rng);
+
+ isource_view = ssp_ranst_discrete_get(rng, sphor->source_distrib_power);
+
+ source_views = darray_source_view_data_get(&sphor->source_views);
+ *source_view = &source_views[isource_view];
+
+ return RES_OK;
+}
+
+res_T
+sample_source_position
+ (struct sphor* sphor,
+ struct ssp_rng* rng,
+ struct source_view** source,
+ struct s3d_primitive* primitive,
+ float st[2])
+{
+ struct source_view* source_to_sample = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != rng);
+
+ res = sample_source_view(sphor, rng, &source_to_sample);
+ if (RES_OK != res) { goto error; }
+
+ res = sample_surface_position(sphor, rng, source_to_sample->view, primitive, st);
+ if (RES_OK != res) { goto error; }
+
+ *source = source_to_sample;
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sphor_ran_source.h b/src/sphor_ran_source.h
@@ -0,0 +1,59 @@
+/* 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_SOURCES_H
+#define SPHOR_RAN_SOURCES_H
+
+#include <rsys/rsys.h>
+
+/* Forward declarations */
+struct sphin_source_surface;
+struct ssp_rng;
+struct s3d_primitive;
+
+extern LOCAL_SYM res_T
+source_surface_sample_direction
+ (struct sphin_source_surface* source,
+ struct ssp_rng* rng,
+ struct s3d_primitive* prim,
+ double dir[3]);
+
+/* Sample a source view considering the discrete distribution of power of all
+ * sources in the config */
+extern LOCAL_SYM res_T
+sample_source_view
+ (struct sphor* sphor,
+ struct ssp_rng* rng,
+ struct source_view** source_view);
+
+/* Sample a source considering the discrete distribution of power and then
+ * sample a position uniformly in this source. */
+extern LOCAL_SYM res_T
+sample_source_position
+ (struct sphor* sphor,
+ struct ssp_rng* rng,
+ struct source_view** source,
+ struct s3d_primitive* primitive,
+ float st[2]);
+
+#endif /* SPHOR_RAN_SOURCES_H */
diff --git a/src/sphor_sources.c b/src/sphor_sources.c
@@ -160,7 +160,6 @@ register_source_view
ASSERT(NULL != sphor);
ASSERT(NULL != args);
- ASSERT(NULL != surface);
source_init(sphor->allocator, &source);
@@ -279,85 +278,3 @@ exit:
error:
goto exit;
}
-
-res_T
-sample_source_view
- (struct sphor* sphor,
- struct ssp_rng* rng,
- const struct source_view** source_view)
-{
- const struct source_view* source_views = NULL;
- size_t isource_view = 0;
-
- ASSERT(NULL != sphor);
- ASSERT(NULL != rng);
-
- isource_view = ssp_ranst_discrete_get(rng, sphor->source_distrib_power);
-
- source_views = darray_source_view_cdata_get(&sphor->source_views);
- *source_view = &source_views[isource_view];
-
- return RES_OK;
-}
-
-res_T
-sample_source_position
- (struct sphor* sphor,
- struct ssp_rng* rng,
- const struct source_view* source_view,
- struct s3d_primitive* primitive,
- double st[2])
-{
- float u = 0;
- float v = 0;
- float w = 0;
- float uv[2];
- res_T res = RES_OK;
-
- ASSERT(NULL != rng);
- ASSERT(NULL != sphor);
- ASSERT(NULL != source_view);
-
- (void)sphor;
-
- u = (float)ssp_rng_canonical(rng);
- v = (float)ssp_rng_canonical(rng);
- w = (float)ssp_rng_canonical(rng);
-
- res = s3d_scene_view_sample(source_view->view, u, v, w, primitive, uv);
- if (RES_OK != res) { goto error; }
-
- st[0] = uv[0];
- st[1] = uv[1];
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-res_T
-sample_source_view_and_position
- (struct sphor* sphor,
- struct ssp_rng* rng,
- struct s3d_primitive* primitive,
- double st[2])
-{
- const struct source_view* source_to_sample = NULL;
- res_T res = RES_OK;
-
- ASSERT(NULL != sphor);
- ASSERT(NULL != rng);
-
- res = sample_source_view(sphor, rng, &source_to_sample);
- if (RES_OK != res) { goto error; }
-
- res = sample_source_position(sphor, rng, source_to_sample, primitive, st);
- if (RES_OK != res) { goto error; }
-
-exit:
- return res;
-error:
- goto exit;
-
-}
diff --git a/src/sphor_sources.h b/src/sphor_sources.h
@@ -35,7 +35,7 @@ struct ssp_rng;
struct source_view {
struct s3d_scene_view* view; /* View of the source */
- size_t sphin_id; /* Source identifier in star-phor-input */
+ size_t sphin_id; /* Surface identifier in star-phor-input */
};
#define SOURCE_VIEW_NULL__ {NULL, 0}
@@ -116,30 +116,4 @@ setup_source_distrib_power
(struct sphor* sphor,
const struct sphor_create_args* args);
-/* Sample a source view considering the discrete distribution of power of all
- * sources in the config */
-extern LOCAL_SYM res_T
-sample_source_view
- (struct sphor* sphor,
- struct ssp_rng* rng,
- const struct source_view** source_view);
-
-/* Sample a position unifomly in a given source */
-extern LOCAL_SYM res_T
-sample_source_position
- (struct sphor* sphor,
- struct ssp_rng* rng,
- const struct source_view* source_view,
- struct s3d_primitive* primitive,
- double st[2]);
-
-/* Sample a source considering the discrete distribution of power and then
- * sample a position uniformly in this source. */
-extern LOCAL_SYM res_T
-sample_source_view_and_position
- (struct sphor* sphor,
- struct ssp_rng* rng,
- struct s3d_primitive* primitive,
- double st[2]);
-
#endif /* SPHOR_SOURCES_H */
