commit 86e82071520b52433c94677b8968df20fc31f616
parent eb05706ceee0a91d2e15c760493e263616745c4f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 5 Feb 2021 11:45:58 +0100
Add the RDGFA phase function
Diffstat:
4 files changed, 417 insertions(+), 3 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -66,6 +66,7 @@ set(SSF_FILES_SRC
ssf_phase.c
ssf_phase_hg.c
ssf_phase_rayleigh.c
+ ssf_phase_rdgfa.c
ssf_pillbox_distribution.c
ssf_specular_dielectric_dielectric_interface.c
ssf_specular_reflection.c
diff --git a/src/ssf.h b/src/ssf.h
@@ -153,8 +153,8 @@ SSF_MICROFACET_DISTRIBUTION_TYPE_NULL = SSF_MICROFACET_DISTRIBUTION_TYPE_NULL__;
* direction `wo' and the incoming direction `wi' point *OUTWARD* the scattering
* point. The scattering angle is thus acos(-wo.wi) */
struct ssf_phase_type {
- res_T (*init)(struct mem_allocator* allocator, void* bsdf); /*Can be NULL*/
- void (*release)(void* bsdf); /* Can be NULL */
+ res_T (*init)(struct mem_allocator* allocator, void* phase); /*Can be NULL*/
+ void (*release)(void* phase); /* Can be NULL */
/* Sample a direction `wi' wrt `wo'. */
void
@@ -185,6 +185,19 @@ struct ssf_phase_type {
#define SSF_PHASE_TYPE_NULL__ {NULL,NULL,NULL,NULL,NULL,0,1}
static const struct ssf_phase_type SSF_PHASE_TYPE_NULL = SSF_PHASE_TYPE_NULL__;
+/* RDGFA phase function input arguments */
+struct ssf_phase_rdgfa_setup_args {
+ double wavelength; /* In nm */
+ double fractal_dimension; /* No unit */
+ double gyration_radius; /* In nm */
+
+ /* #intervals used to discretized the both angular domains */
+ size_t nintervals[2];
+};
+#define SSF_PHASE_RDGFA_SETUP_DEFAULT_ARGS__ {0, 0, 0, {1000, 1000}}
+static const struct ssf_phase_rdgfa_setup_args
+SSF_PHASE_RDGFA_SETUP_DEFAULT_ARGS = SSF_PHASE_RDGFA_SETUP_DEFAULT_ARGS__;
+
BEGIN_DECLS
/*******************************************************************************
@@ -309,6 +322,15 @@ SSF_API const struct ssf_phase_type ssf_phase_hg;
* p(wo, wi) = 3/(16*PI) * (1+(-wo.wi)^2) */
SSF_API const struct ssf_phase_type ssf_phase_rayleigh;
+/* RDGFA phase function normalized to 1:
+ * p(wo, wi) = 3/(16*PI) * f(theta)/g * (1+cos(theta)^2);
+ * theta = acos(wo.wi)
+ *
+ * "Effects of multiple scattering on radiative properties of soot fractal
+ * aggregates" J. Yon et al., Journal of Quantitative Spectroscopy and
+ * Radiative Transfer Vol 133, pp. 374-381, 2014 */
+SSF_API const struct ssf_phase_type ssg_phase_rdgfa;
+
/*******************************************************************************
* BSDF API - Bidirectional Scattering Distribution Function. Describes the way
* the light is scattered by a surface. Note that by convention the outgoing
@@ -442,6 +464,11 @@ ssf_phase_hg_setup
(struct ssf_phase* phase,
const double g); /* Asymmetric parameter in [-1, 1] */
+SSF_API res_T
+ssf_phase_rdgfa_setup
+ (struct ssf_phase* phase,
+ const struct ssf_phase_rdgfa_setup_args* args);
+
/*******************************************************************************
* Fresnel API - Define the equation of the fresnel term
******************************************************************************/
diff --git a/src/ssf_phase_hg.c b/src/ssf_phase_hg.c
@@ -110,4 +110,3 @@ ssf_phase_hg_setup(struct ssf_phase* phase, const double g)
((struct hg*)phase->data)->g = g;
return RES_OK;
}
-
diff --git a/src/ssf_phase_rdgfa.c b/src/ssf_phase_rdgfa.c
@@ -0,0 +1,387 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * 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 "ssf.h"
+#include "ssf_phase_c.h"
+
+#include <star/ssp.h>
+
+#include <rsys/algorithm.h>
+#include <rsys/dynamic_array_double.h>
+
+#define EXP1 2.7182818284590452354
+
+struct rdgfa {
+ double wavelength; /* In nm */
+ double fractal_dimension; /* No unit */
+ double gyration_radius; /* In nm */
+
+ double rcp_normalize_factor; /* Reciprocal normalization factor of the CDF */
+
+ /* Discretized cumulative (#entries = nintervals[0] + nintervals[1]) */
+ struct darray_double cdf;
+
+ /* Precomputed values */
+ double Rg2; /* gyration_radius^2 */
+ double cst_4pi_div_lambda; /* (4*PI)/wavelength */
+ double cst_3Df_div_2E; /* 3*Df/(2*exp(1)) */
+ double g; /* g function */
+
+ /* #intervals used to discretized both the [0,theta_limit[ and
+ * [theta_limit,pi[ intervals */
+ size_t nintervals[2];
+
+ /* Length of an angular interval for the tow angular domains */
+ double dtheta[2]; /* In rad */
+
+ /* Limite angle that split the 2 angular domains */
+ double theta_limit; /* In rad */
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE int
+check_phase_rdgfa_setup_args(const struct ssf_phase_rdgfa_setup_args* args)
+{
+ return args
+ && args->wavelength > 0
+ && args->fractal_dimension > 0
+ && args->gyration_radius > 0
+ && args->nintervals[0] != 0
+ && args->nintervals[1] != 0;
+}
+
+static INLINE int
+cmp_dbl(const void* a, const void* b)
+{
+ const double key = *((double*)a);
+ const double val = *((double*)b);
+ if(key < val) {
+ return -1;
+ } else if(key > val) {
+ return 1;
+ } else {
+ return 0;
+ }
+}
+
+#if 0
+/* Not normalized */
+static INLINE double
+eval__(struct rdgfa* rdgfa, const double theta)
+{
+ /* Input arguments */
+ const double lambda = rdgfa->wavelength;
+ const double Df = rdgfa->fractal_dimension;
+ const double Rg = rdgfa->gyration_radius;
+ const double Rg2 = Rg*Rg;
+
+ /* TODO precompute constants */
+ const double k = 2*PI / lambda; /* [nm^-1] */
+ const double k2 = k*k;
+
+ const double g = pow(1 + 4*k2*Rg2/(3*Df), -Df*0.5);
+
+ const double q = 4*PI / lambda * sin(theta*0.5);
+ const double qRg = q * Rg;
+ const double qRg2 = qRg*qRg;
+ const double f = qRg2 < 1.5*Df
+ ? exp(-1.0/3.0 * qRg2)
+ : pow(3.0*Df/(2.0*EXP1) * 1.0/qRg2, Df*0.5);
+
+ ASSERT(d3_is_normalized(wo) && d3_is_normalized(wi));
+ return f;
+}
+#endif
+
+static INLINE double
+eval_f(struct rdgfa* rdgfa, const double theta)
+{
+ double Df, Rg2, q, q2Rg2, f;
+ ASSERT(rdgfa);
+
+ /* Input arguments */
+ Df = rdgfa->fractal_dimension;
+
+ /* Fech precomputed constants */
+ Rg2 = rdgfa->Rg2;
+
+ /* Precompute values */
+ q = rdgfa->cst_4pi_div_lambda * sin(theta*0.5);
+ q2Rg2 = q*q*Rg2;
+
+ /* Evaluate f(theta) */
+ f = q2Rg2 < 1.5*Df
+ ? exp(-1.0/3.0 * q2Rg2)
+ : pow(rdgfa->cst_3Df_div_2E * 1/q2Rg2, Df*0.5);
+ return f;
+}
+
+/* Not normalized */
+static INLINE double
+eval2
+ (struct rdgfa* rdgfa,
+ const double theta,
+ const double cos_theta)
+{
+ double f, cos2_theta, phase;
+ ASSERT(rdgfa && eq_eps(cos_theta, cos(theta), 1.e-6));
+
+ /* Precompute values */
+ cos2_theta = cos_theta * cos_theta;
+
+ /* Evaluate the phase(theta) */
+ f = eval_f(rdgfa, theta);
+ phase = 3.0/(16*PI) * f / rdgfa->g * (1 + cos2_theta);
+ return phase;
+}
+
+
+/* Not normalized */
+static FINLINE double
+eval(struct rdgfa* rdgfa, const double theta)
+{
+ return eval2(rdgfa, theta, cos(theta));
+}
+
+static res_T
+compute_cumulative(struct rdgfa* rdgfa)
+{
+ double* cdf = NULL;
+ double Df, Rg;
+ double f1;
+ double theta1;
+ size_t i;
+ size_t nintervals_total;
+ res_T res = RES_OK;
+
+ /* Compute the overall number of angular intervals */
+ nintervals_total = rdgfa->nintervals[0] + rdgfa->nintervals[1];
+
+ /* Allocate the cumulative array */
+ res = darray_double_resize(&rdgfa->cdf, nintervals_total);
+ if(res != RES_OK) goto error;
+ cdf = darray_double_data_get(&rdgfa->cdf);
+
+ /* Fetch input arguments */
+ Df = rdgfa->fractal_dimension;
+ Rg = rdgfa->gyration_radius;
+
+ /* Compute the limit angle between the 2 angular domains */
+ rdgfa->theta_limit = 2 * asin(rdgfa->wavelength/(4*PI) * sqrt(1.5*Df) / Rg);
+
+ /* Compute the angular step for the 2 angular domains */
+ rdgfa->dtheta[0] = rdgfa->theta_limit / (double)rdgfa->nintervals[0];
+ rdgfa->dtheta[1] = (PI - rdgfa->theta_limit) / (double)rdgfa->nintervals[1];
+
+ theta1 = 0;
+ f1 = eval(rdgfa, 0);
+ FOR_EACH(i, 0, nintervals_total) {
+ /* Fetch the current angular step */
+ const double dtheta_curr = i < rdgfa->nintervals[0]
+ ? rdgfa->dtheta[0] : rdgfa->dtheta[1];
+
+ /* Compute the upper bound of the currant angular range */
+ const double theta2 = theta1 + dtheta_curr;
+
+ /* Compute the (unormalized) cumulative for current interval */
+ const double delta_omega = 2*PI*sin((theta1+theta2)*0.5)*dtheta_curr;
+ const double f2 = eval(rdgfa, theta2);
+ const double tmp = (f1 + f2) * 0.5 * delta_omega;
+ cdf[i] = (i == 0 ? tmp : tmp + cdf[i-1]);
+
+ /* Go to the next interval */
+ f1 = f2;
+ theta1 = theta2;
+ }
+
+ /* Save the normamlization factor */
+ rdgfa->rcp_normalize_factor = 1.0 / cdf[nintervals_total-1];
+
+ /* Finally normalize the CDF */
+ FOR_EACH(i, 0, nintervals_total) {
+ cdf[i] *= rdgfa->rcp_normalize_factor;
+ }
+
+exit:
+ return res;
+error:
+ darray_double_clear(&rdgfa->cdf);
+ goto exit;
+}
+
+/*******************************************************************************
+ * Private functions
+ ******************************************************************************/
+static res_T
+rdgfa_init(struct mem_allocator* allocator, void* phase)
+{
+ struct rdgfa* rdgfa = phase;
+ ASSERT(phase);
+ memset(rdgfa, 0, sizeof(*rdgfa));
+ darray_double_init(allocator, &rdgfa->cdf);
+ return RES_OK;
+}
+
+static void
+rdgfa_release(void* phase)
+{
+ struct rdgfa* rdgfa = phase;
+ ASSERT(phase);
+ darray_double_release(&rdgfa->cdf);
+}
+
+static double
+rdgfa_eval(void* data, const double wo[3], const double wi[3])
+{
+ const struct rdgfa* rdgfa = data;
+ const double cos_theta = d3_dot(wo, wi);
+ const double theta = acos(cos_theta);
+ return eval2(data, theta, cos_theta) * rdgfa->rcp_normalize_factor;
+}
+
+static void
+rdgfa_sample
+ (void* data,
+ struct ssp_rng* rng,
+ const double wo[3],
+ double wi[3],
+ double* pdf)
+{
+ struct rdgfa* rdgfa = data;
+ const double* cdf = NULL;
+ double* find = NULL;
+ double frame[9];
+ double w[3];
+ double theta1;
+ double theta2;
+ double theta;
+ double sin_theta;
+ double phi;
+ double r;
+ double u;
+ size_t n;
+ size_t i;
+ ASSERT(data && rng && wo && wi);
+
+ /* Fetch the CDF array and its number of entries */
+ cdf = darray_double_cdata_get(&rdgfa->cdf);
+ n = darray_double_size_get(&rdgfa->cdf);
+
+ /* Sample the CDF */
+ r = ssp_rng_canonical(rng);
+ find = search_lower_bound(&r, cdf, n, sizeof(double), cmp_dbl);
+ ASSERT(find && (*find) >= r);
+ i = (size_t)(find - cdf);
+
+ /* Compute the angle interval into which the sample lies */
+ if(i < rdgfa->nintervals[0]) {
+ theta1 = rdgfa->dtheta[0] * (double)(i + 0);
+ theta2 = rdgfa->dtheta[0] * (double)(i + 1);
+ } else {
+ theta1 = rdgfa->dtheta[1] * (double)(i - rdgfa->nintervals[0] + 0);
+ theta2 = rdgfa->dtheta[1] * (double)(i - rdgfa->nintervals[0] + 1);
+ theta1 = theta1 + rdgfa->theta_limit;
+ theta2 = theta2 + rdgfa->theta_limit;
+ }
+
+ /* Compute the sampled theta angle by linearly interpolate it from the
+ * boundaries of its interval */
+ if(i == 0) {
+ u = r / cdf[0];
+ } else {
+ u = (r - cdf[i-1]) / (cdf[i] - cdf[i-1]);
+ }
+ ASSERT(0 <= u && u < 1);
+ theta = u * (theta2 - theta1) + theta1;
+
+ /* Uniformly sample a phi angle in [0, 2PI[ */
+ phi = ssp_rng_uniform_double(rng, 0, 2*PI);
+ sin_theta = sin(theta);
+
+ /* Compute the cartesian coordinates of the sampled direction in the _local_
+ * phase function space */
+ wi[0] = cos(phi) * sin_theta;
+ wi[1] = sin(phi) * sin_theta;
+ wi[2] = cos(theta);
+
+ /* Compute the transformation matrix from local phase function to world
+ * space. Note that by convention, in Star-SF the directions point outward
+ * the scattering position. Revert 'wo' to match the convention used by the
+ * previous computations */
+ d33_basis(frame, d3_minus(w, wo));
+ d33_muld3(wi, frame, wi);
+
+ if(pdf) *pdf = rdgfa_eval(rdgfa, wo, wi);
+}
+
+/*******************************************************************************
+ * Exported symbols
+ ******************************************************************************/
+const struct ssf_phase_type ssf_phase_rdgfa = {
+ rdgfa_init,
+ rdgfa_release,
+ rdgfa_sample,
+ rdgfa_eval,
+ rdgfa_eval,
+ sizeof(struct rdgfa),
+ ALIGNOF(struct rdgfa)
+};
+
+res_T
+ssf_phase_rdgfa_setup
+ (struct ssf_phase* phase,
+ const struct ssf_phase_rdgfa_setup_args* args)
+{
+ struct rdgfa* rdgfa = NULL;
+ double lambda, Df, Rg, k, k2;
+ res_T res = RES_OK;
+
+ if(!phase || !check_phase_rdgfa_setup_args(args)) {
+ res = RES_OK;
+ goto error;
+ }
+
+ rdgfa = phase->data;
+ rdgfa->wavelength = args->wavelength;
+ rdgfa->fractal_dimension = args->fractal_dimension;
+ rdgfa->gyration_radius = args->gyration_radius;
+ rdgfa->nintervals[0] = args->nintervals[0];
+
+ /* Fetch input data */
+ lambda = rdgfa->wavelength;
+ Df = rdgfa->fractal_dimension;
+ Rg = rdgfa->gyration_radius;
+
+ /* Precompute constants */
+ rdgfa->Rg2 = Rg*Rg; /* gyration_radius^2 */
+ rdgfa->cst_4pi_div_lambda = 4*PI / lambda;
+ rdgfa->cst_3Df_div_2E = 3*Df/(2.0*EXP1);
+
+ /* Precompute the function g */
+ k = (2.0*PI) / lambda;
+ k2 = k*k;
+ rdgfa->g = pow(1 + 4*k2*rdgfa->Rg2/(3*Df), -Df*0.5);
+
+ /* Precompute the phase function cumulative */
+ res = compute_cumulative(rdgfa);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
