ssf_phase_rayleigh.c (3027B)
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 #define _POSIX_C_SOURCE 200112L /* cbrt support */ 17 18 #include "ssf.h" 19 #include "ssf_phase_c.h" 20 21 #include <rsys/double3.h> 22 #include <star/ssp.h> 23 #include <math.h> 24 25 /******************************************************************************* 26 * Private functions 27 ******************************************************************************/ 28 static res_T 29 rayleigh_init(struct mem_allocator* allocator, void* phase) 30 { 31 (void)allocator, (void)phase; 32 return RES_OK; 33 } 34 35 static void 36 rayleigh_release(void* phase) 37 { (void)phase; } 38 39 40 /* Rayleigh(theta) = 3/(16*PI)*(1+cos(theta)^2) */ 41 static double 42 rayleigh_eval(void* data, const double wo[3], const double wi[3]) 43 { 44 double cos_theta; 45 double w[3]; 46 ASSERT(wo && wi); 47 ASSERT(d3_is_normalized(wo) && d3_is_normalized(wi)); 48 (void)data; 49 /* By convention wo points outward the scattering point. Revert it to point 50 * inward the scattering point in order to compute the cosine of the 51 * scattering angle */ 52 d3_minus(w, wo); 53 cos_theta = d3_dot(w, wi); 54 return 3.0/(16.0*PI)*(1.0+cos_theta*cos_theta); 55 } 56 57 static void 58 rayleigh_sample 59 (void* data, 60 struct ssp_rng* rng, 61 const double wo[3], 62 double wi[3], 63 double* pdf) 64 { 65 double frame[9]; 66 double sample[3]; 67 double w[3]; 68 double phi, cos_theta, sin_theta; 69 double u; 70 double s; 71 ASSERT(rng && wo && wi); 72 (void)data; 73 74 phi = ssp_rng_uniform_double(rng, 0, 2*PI); 75 76 u = 4*ssp_rng_canonical(rng)-2; 77 s = cbrt(u + sqrt(1+u*u)); 78 cos_theta = s - 1.0/s; 79 sin_theta = cos2sin(cos_theta); 80 81 sample[0] = cos(phi) * sin_theta; 82 sample[1] = sin(phi) * sin_theta; 83 sample[2] = cos_theta; 84 85 /* By convention wo points outward the scattering point. Revert it to point 86 * inward the scattering point */ 87 d3_minus(w, wo); 88 d33_muld3(wi, d33_basis(frame, w), sample); 89 90 if(pdf) *pdf = rayleigh_eval(data, wo, wi); 91 } 92 93 static double 94 rayleigh_pdf(void* data, const double wo[3], const double wi[3]) 95 { 96 return rayleigh_eval(data, wo, wi); 97 } 98 99 /******************************************************************************* 100 * Exported symbols 101 ******************************************************************************/ 102 const struct ssf_phase_type ssf_phase_rayleigh = { 103 rayleigh_init, 104 rayleigh_release, 105 rayleigh_sample, 106 rayleigh_eval, 107 rayleigh_pdf, 108 0, 109 1 110 }; 111