commit 66cb2687eda5448174d0af13f3ccd3a20cf4eb3d
parent f17035cbb3f42d9abda8a89d1cae3ff99a0ffc46
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 15 Sep 2021 10:22:33 +0200
Add an exemple that computes the translation sensitivity
Diffstat:
6 files changed, 277 insertions(+), 17 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -29,7 +29,7 @@ option(NO_TEST "Do not build tests" OFF)
################################################################################
find_package(RCMake 0.4 REQUIRED)
find_package(RSys 0.12 REQUIRED)
-find_package(Star3D 0.7 REQUIRED)
+find_package(Star3D 0.8 REQUIRED)
find_package(StarMC 0.4 REQUIRED)
find_package(StarSP 0.9 REQUIRED)
@@ -55,6 +55,7 @@ set(SGS_FILES_SRC
sgs_args.c
sgs.c
sgs_compute_4v_s.c
+ sgs_compute_trans_sensib.c
sgs_geometry.c
sgs_geometry_box.c
sgs_geometry_slope.c
diff --git a/src/sgs.c b/src/sgs.c
@@ -144,7 +144,11 @@ sgs_run(struct sgs* sgs)
res = sgs_geometry_dump_vtk(sgs->geom, stdout);
if(res != RES_OK) goto error;
} else {
+#if 1
+ res = compute_trans_sensib(sgs);
+#else
res = compute_4v_s(sgs);
+#endif
if(res != RES_OK) goto error;
}
diff --git a/src/sgs_c.h b/src/sgs_c.h
@@ -20,6 +20,7 @@
#ifndef SGS_C_H
#define SGS_C_H
+#include <rsys/double3.h>
#include <rsys/logger.h>
#include <rsys/rsys.h>
@@ -37,6 +38,21 @@ struct sgs {
struct mem_allocator* allocator;
};
+/* Reflect the vector V wrt the normal N. By convention V points outward the
+ * surface. */
+static INLINE double*
+reflect(double res[3], const double V[3], const double N[3])
+{
+ double tmp[3];
+ double cos_V_N;
+ ASSERT(res && V && N);
+ ASSERT(d3_is_normalized(V) && d3_is_normalized(N));
+ cos_V_N = d3_dot(V, N);
+ d3_muld(tmp, N, 2*cos_V_N);
+ d3_sub(res, tmp, V);
+ return res;
+}
+
extern LOCAL_SYM void
setup_logger
(struct sgs* sgs);
@@ -45,4 +61,8 @@ extern LOCAL_SYM res_T
compute_4v_s
(struct sgs* sgs);
+extern LOCAL_SYM res_T
+compute_trans_sensib
+ (struct sgs* sgs);
+
#endif /* SGS_C_H */
diff --git a/src/sgs_compute_4v_s.c b/src/sgs_compute_4v_s.c
@@ -31,21 +31,6 @@
/*******************************************************************************
* Helper function
******************************************************************************/
-/* Reflect the vector V wrt the normal N. By convention V points outward the
- * surface. */
-static INLINE double*
-reflect(double res[3], const double V[3], const double N[3])
-{
- double tmp[3];
- double cos_V_N;
- ASSERT(res && V && N);
- ASSERT(d3_is_normalized(V) && d3_is_normalized(N));
- cos_V_N = d3_dot(V, N);
- d3_muld(tmp, N, 2*cos_V_N);
- d3_sub(res, tmp, V);
- return res;
-}
-
static res_T
realisation(void* weight, struct ssp_rng* rng, const unsigned ithread, void* ctx)
{
diff --git a/src/sgs_compute_trans_sensib.c b/src/sgs_compute_trans_sensib.c
@@ -0,0 +1,249 @@
+/* Copyright (C) 2021
+ * CNRS/RAPSODEE,
+ * CNRS/LMAP,
+ * |Meso|Star> (contact@meso-star.com),
+ * UPS/Laplace.
+ *
+ * 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 "sgs_c.h"
+#include "sgs_geometry.h"
+#include "sgs_log.h"
+
+#include <rsys/cstr.h>
+
+#include <star/smc.h>
+#include <star/ssp.h>
+
+enum {X, Y, Z};
+
+static const double RECV_MIN[2] = {0.5, 1.5};
+static const double RECV_MAX[2] = {1.5, 2.5};
+static const double EMIT_THRESHOLD = 2;
+static const double V[3] = {0, 0, 1};
+
+/* FIXME this should be a variable */
+static const double EMIT_SZ[3] = {4, 4, 0};
+static const double EMIT_E_SZ[3] = {0, 4, 2};
+
+#if 0
+static const double RECV_MIN[2] = {0.125, 0.375};
+static const double RECV_MAX[2] = {0.375, 0.625};
+static const double EMIT_THRESHOLD = 0.5;
+static const double V[3] = {0, 0, 1};
+#endif
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+realisation(void* weight, struct ssp_rng* rng, const unsigned ithread, void* ctx)
+{
+ struct sgs_hit hit = SGS_HIT_NULL;
+
+ struct sgs_fragment frag = SGS_FRAGMENT_NULL;
+ struct sgs* sgs = ctx;
+
+ double dir_emit[3] = {0, 0, 0};
+ double dir_spec[3] = {0, 0, 0};
+ double pos_emit[3] = {0, 0, 0};
+ double pos_recv[3] = {0, 0, 0};
+ double pos_emit_e[3] = {0, 0, 0};
+ double dir_spec_e[3] = {0, 0, 0};
+
+ double range[2] = {0, DBL_MAX};
+
+ double u_emit[3];
+ double u_spec[3];
+ double u_emit_e[3];
+ double u_spec_e[3];
+ double tmp[3];
+ double beta_emit;
+ double beta_spec;
+ double alpha_emit;
+ double alpha_spec;
+ double beta_emit_e;
+ double beta_spec_e;
+ double alpha_emit_e;
+ double alpha_spec_e;
+
+ double rho;
+ double grad_rho[3];
+ double I;
+ double grad_I[3];
+
+ double S;
+
+ double w = 0;
+
+ enum sgs_surface_type surf_emit = SGS_SURFACE_NONE;
+
+ res_T res = RES_OK;
+ (void)ithread;
+
+ /* Sample the sensitivity emissive surface */
+ sgs_geometry_sample(sgs->geom, rng, &frag);
+ surf_emit = frag.surface;
+
+ /* Sample the cosine weighted sampling of the emissive direction */
+ ssp_ran_hemisphere_cos(rng, frag.normal, dir_emit, NULL);
+
+ pos_emit[X] = frag.position[X];
+ pos_emit[Y] = frag.position[Y];
+ pos_emit[Z] = frag.position[Z];
+
+ range[0] = 0;
+ range[1] = INF;
+
+ sgs_geometry_trace_ray(sgs->geom, pos_emit, dir_emit, range, surf_emit, &hit);
+ if(SGS_HIT_NONE(&hit)) {
+ res = RES_OK;
+ goto error;
+ }
+ pos_recv[X] = pos_emit[X] + hit.distance*dir_emit[X];
+ pos_recv[Y] = pos_emit[Y] + hit.distance*dir_emit[Y];
+ pos_recv[Z] = pos_emit[Z] + hit.distance*dir_emit[Z];
+
+ /* Does not reach the receiver */
+ if(hit.surface != SGS_SURFACE_Z_NEG
+ || RECV_MIN[X] > pos_recv[X] || pos_recv[X] > RECV_MAX[X]
+ || RECV_MIN[Y] > pos_recv[Y] || pos_recv[Y] > RECV_MAX[Y]) {
+ goto exit;
+ }
+
+ reflect(dir_spec, dir_emit, frag.normal);
+ sgs_geometry_trace_ray(sgs->geom, pos_emit, dir_spec, range, surf_emit, &hit);
+ if(SGS_HIT_NONE(&hit)) {
+ res = RES_OK;
+ goto error;
+ }
+ pos_emit_e[X] = pos_emit[X] + hit.distance*dir_spec[X];
+ pos_emit_e[Y] = pos_emit[Y] + hit.distance*dir_spec[Y];
+ pos_emit_e[Z] = pos_emit[Z] + hit.distance*dir_spec[Z];
+
+ if(hit.surface != SGS_SURFACE_X_POS || pos_emit_e[Z] > EMIT_THRESHOLD) {
+ goto exit;
+ }
+
+ alpha_emit = d3_dot(V, frag.normal) / d3_dot(dir_emit, frag.normal);
+ alpha_spec = d3_dot(V, frag.normal) / d3_dot(dir_spec, frag.normal);
+ d3_sub(u_emit, V, d3_muld(tmp, dir_emit, alpha_emit));
+ d3_sub(u_spec, V, d3_muld(tmp, dir_spec, alpha_spec));
+ beta_emit = d3_normalize(u_emit, u_emit);
+ beta_spec = d3_normalize(u_spec, u_spec);
+
+ d3_minus(dir_spec_e, dir_spec);
+
+ alpha_emit_e = d3_dot(u_emit, hit.normal) / d3_dot(dir_spec_e, hit.normal);
+ alpha_spec_e = d3_dot(u_spec, hit.normal) / d3_dot(dir_spec_e, hit.normal);
+ d3_sub(u_emit_e, u_emit, d3_muld(tmp, dir_spec_e, alpha_emit_e));
+ d3_sub(u_spec_e, u_spec, d3_muld(tmp, dir_spec_e, alpha_spec_e));
+ beta_emit_e = d3_normalize(u_emit_e, u_emit_e);
+ beta_spec_e = d3_normalize(u_spec_e, u_spec_e);
+
+ rho = 0.25
+ * (1 - cos(2*PI*pos_emit[X]/EMIT_SZ[X]))
+ * (1 - cos(2*PI*pos_emit[Y]/EMIT_SZ[Y]));
+ grad_rho[X] = 0.25
+ * (((2*PI)/EMIT_SZ[X])*sin(2*PI*pos_emit[X]/EMIT_SZ[X]))
+ * (1 - cos(2*PI*pos_emit[Y]/EMIT_SZ[Y]));
+ grad_rho[Y] = 0.25
+ * (((2*PI)/EMIT_SZ[Y])*sin(2*PI*pos_emit[Y]/EMIT_SZ[Y]))
+ * (1 - cos(2*PI*pos_emit[X]/EMIT_SZ[X]));
+ grad_rho[Z] = 0;
+
+ I =
+ (1 - cos(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]))
+ * (1 - cos(2*PI*pos_emit_e[Z]/EMIT_E_SZ[Z]));
+ grad_I[X] = 0;
+ grad_I[Y] =
+ (((2*PI)/EMIT_E_SZ[Y])*sin(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]))
+ * (1 - cos(2*PI*pos_emit_e[Z]/EMIT_E_SZ[Z]));
+ grad_I[Z] =
+ (((2*PI)/EMIT_E_SZ[Z])*sin(2*PI*pos_emit_e[Z]/EMIT_E_SZ[Z]))
+ * (1 - cos(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]));
+
+ S =
+ - beta_emit * d3_dot(grad_rho, u_emit) * I
+ - rho * beta_emit * beta_emit_e * d3_dot(grad_I, u_emit_e)
+ + rho * beta_spec * beta_spec_e * d3_dot(grad_I, u_spec_e);
+
+ /* w = I*EMIT_SZ[X]*EMIT_SZ[Y]*PI;*/ /* Weight */
+ w = S*EMIT_SZ[X]*EMIT_SZ[Y]*PI; /* Sensib */
+
+exit:
+ SMC_DOUBLE(weight) = w;
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+res_T
+compute_trans_sensib(struct sgs* sgs)
+{
+ struct smc_integrator integrator;
+ struct smc_estimator_status status;
+ struct smc_device* smc = NULL;
+ struct smc_estimator* estimator = NULL;
+ res_T res = RES_OK;
+ ASSERT(sgs);
+
+ /* Create the Star-MonteCarlo device */
+ res = smc_device_create
+ (&sgs->logger, sgs->allocator, sgs->nthreads, &ssp_rng_mt19937_64, &smc);
+ if(res != RES_OK) {
+ sgs_log_err(sgs, "Could not create the Star-MonteCarlo device -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ /* Setup the integrator */
+ integrator.integrand = realisation;
+ integrator.type = &smc_double;
+ integrator.max_realisations = sgs->nrealisations;
+ integrator.max_failures = (size_t)((double)sgs->nrealisations*0.01);
+
+ /* Run Monte-Carlo integration */
+ res = smc_solve(smc, &integrator, sgs, &estimator);
+ if(res != RES_OK) {
+ sgs_log_err(sgs, "Error while solving 4V/S -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ /* Retrieve the estimated value */
+ SMC(estimator_get_status(estimator, &status));
+ if(status.NF >= integrator.max_failures) {
+ sgs_log_err(sgs, "Too many failures: %lu\n", (unsigned long)status.NF);
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ /* Print the result */
+ sgs_log(sgs, "%g +/- %g; #failures = %lu/%lu\n",
+ SMC_DOUBLE(status.E),
+ SMC_DOUBLE(status.SE),
+ (unsigned long)status.NF,
+ (unsigned long)sgs->nrealisations);
+
+exit:
+ if(smc) SMC(device_ref_put(smc));
+ if(estimator) SMC(estimator_ref_put(estimator));
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sgs_geometry.c b/src/sgs_geometry.c
@@ -54,12 +54,13 @@ hit_filter
(const struct s3d_hit* hit,
const float ray_org[3],
const float ray_dir[3],
+ const float ray_range[2],
void* ray_data,
void* filter_data)
{
const struct hit_filter_context* ctx = ray_data;
enum sgs_surface_type surface;
- (void)ray_org, (void)ray_dir, (void)filter_data;
+ (void)ray_org, (void)ray_dir, (void)ray_range, (void)filter_data;
if(!ctx) /* Nothing to do */
return 0;
