star-sf

Set of surface and volume scattering functions
git clone git://git.meso-star.fr/star-sf.git
Log | Files | Refs | README | LICENSE

ssf_phase_hg.c (3134B)


      1 /* Copyright (C) 2016-2018, 2021-2025 |Méso|Star> (contact@meso-star.com)
      2  *
      3  * This program is free software: you can redistribute it and/or modify
      4  * it under the terms of the GNU General Public License as published by
      5  * the Free Software Foundation, either version 3 of the License, or
      6  * (at your option) any later version.
      7  *
      8  * This program is distributed in the hope that it will be useful,
      9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     11  * GNU General Public License for more details.
     12  *
     13  * You should have received a copy of the GNU General Public License
     14  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     15 
     16 #include "ssf.h"
     17 #include "ssf_phase_c.h"
     18 
     19 #include <rsys/double3.h>
     20 
     21 #include <star/ssp.h>
     22 
     23 struct hg { double g; };
     24 
     25 /*******************************************************************************
     26  * Private functions
     27  ******************************************************************************/
     28 static res_T
     29 hg_init(struct mem_allocator* allocator, void* phase)
     30 {
     31   ASSERT(phase);
     32   (void)allocator;
     33   ((struct hg*)phase)->g = 0.0;
     34   return RES_OK;
     35 }
     36 
     37 static double
     38 hg_pdf
     39   (void* data,
     40    const double wo[3],
     41    const double wi[3])
     42 {
     43   double w[3];
     44   const struct hg* hg = data;
     45   ASSERT(data && wo && wi);
     46   return ssp_ran_sphere_hg_pdf(d3_minus(w, wo), hg->g, wi);
     47 }
     48 
     49 /* HG(theta) =  1/(4*PI) * (1 - g^2) / (1 + g^2 - 2*g*cos(theta))^3/2 */
     50 static double
     51 hg_eval(void* data, const double wo[3], const double wi[3])
     52 {
     53   const struct hg* hg = data;
     54   double w[3];
     55   double g;
     56   double cos_theta;
     57   double denom;
     58   double val;
     59   ASSERT(data && wo && wi);
     60   ASSERT(d3_is_normalized(wo) && d3_is_normalized(wi));
     61 
     62   g = hg->g;
     63   /* By convention wo points outward the scattering point. Revert it to point
     64    * inward the scattering point in order to compute the cosine of the
     65    * scattering angle */
     66   d3_minus(w, wo);
     67   cos_theta = d3_dot(w, wi);
     68   denom = 1 + g*g - 2*g*cos_theta;
     69   ASSERT(denom != 0);
     70   val = 1.0/(4.0*PI) * (1 - g*g) / (denom*sqrt(denom));
     71   ASSERT(eq_eps(val, hg_pdf(data, wo, wi), 1.e-6));
     72   return val;
     73 }
     74 
     75 static void
     76 hg_sample
     77   (void* data,
     78    struct ssp_rng* rng,
     79    const double wo[3],
     80    double wi[3],
     81    double* pdf)
     82 {
     83   const struct hg* hg = data;
     84   double w[3];
     85   ASSERT(data && wo && wi);
     86 
     87   /* By convention wo points outward the scattering point. Revert it to point
     88    * inward the scattering point as expected by the SSP library */
     89   ssp_ran_sphere_hg(rng, d3_minus(w, wo), hg->g, wi, pdf);
     90 }
     91 
     92 /*******************************************************************************
     93  * Exported symbols
     94  ******************************************************************************/
     95 const struct ssf_phase_type ssf_phase_hg = {
     96   hg_init,
     97   NULL,
     98   hg_sample,
     99   hg_eval,
    100   hg_pdf,
    101   sizeof(struct hg),
    102   ALIGNOF(struct hg)
    103 };
    104 
    105 res_T
    106 ssf_phase_hg_setup(struct ssf_phase* phase, const double g)
    107 {
    108   if(!phase || g < -1 || g > 1) return RES_BAD_ARG;
    109   if(!PHASE_TYPE_EQ(&phase->type, &ssf_phase_hg)) return RES_BAD_ARG;
    110   ((struct hg*)phase->data)->g = g;
    111   return RES_OK;
    112 }