test_ssf_phase_hg.c (3129B)
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 "test_ssf_utils.h" 18 19 #define NG 10 20 #define NSAMPS 10000 21 22 int 23 main(int argc, char** argv) 24 { 25 struct mem_allocator allocator; 26 struct ssf_phase* phase; 27 struct ssf_phase* dummy; 28 struct ssp_rng* rng; 29 double wo[3]; 30 double wi[3]; 31 size_t i; 32 (void)argc, (void)argv; 33 34 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 35 CHK(ssp_rng_create(&allocator, SSP_RNG_MT19937_64, &rng) == RES_OK); 36 37 CHK(ssf_phase_create(&allocator, &ssf_phase_hg, &phase) == RES_OK); 38 CHK(ssf_phase_create(&allocator, &phase_dummy, &dummy) == RES_OK); 39 40 CHK(ssf_phase_hg_setup(NULL, -2) == RES_BAD_ARG); 41 CHK(ssf_phase_hg_setup(phase, -2) == RES_BAD_ARG); 42 CHK(ssf_phase_hg_setup(NULL, -1) == RES_BAD_ARG); 43 CHK(ssf_phase_hg_setup(phase, -1) == RES_OK); 44 CHK(ssf_phase_hg_setup(phase, 0) == RES_OK); 45 CHK(ssf_phase_hg_setup(phase, 1) == RES_OK); 46 CHK(ssf_phase_hg_setup(dummy, 1) == RES_BAD_ARG); 47 48 FOR_EACH(i, 0, NG) { 49 const double g = ssp_rng_uniform_double(rng, -1, +1); 50 double sum_cos = 0; 51 double sum_cos_sqr = 0; 52 double E, V, SE; 53 size_t isamp; 54 55 CHK(ssf_phase_hg_setup(phase, g) == RES_OK); 56 57 ssp_ran_sphere_uniform(rng, wo, NULL); 58 FOR_EACH(isamp, 0, NSAMPS) { 59 double w[3]; 60 double weight; 61 double pdf; 62 double r; 63 double ref; 64 ssf_phase_sample(phase, rng, wo, wi, &pdf); 65 CHK(d3_is_normalized(wi)); 66 CHK(eq_eps(ssf_phase_pdf(phase, wo, wi), pdf, 1.e-6)); 67 d3_minus(w, wo); /* Match the convention */ 68 weight = d3_dot(w, wi); 69 sum_cos += weight; 70 sum_cos_sqr += weight*weight; 71 72 /* HG(theta) = 1/(4*PI) * (1 - g^2) / (1 + g^2 - 2*g*cos(theta))^3/2 */ 73 r = ssf_phase_eval(phase, wo, wi); 74 ref = 1/(4*PI) * (1-g*g) / pow(1+g*g-2*g*d3_dot(w,wi), 1.5); 75 CHK(eq_eps(r, ref, 1.e-6)); 76 CHK(eq_eps(r, pdf, 1.e-6)); 77 } 78 /* On average cos(-wo,wi) should be g */ 79 E = sum_cos / NSAMPS; 80 V = sum_cos_sqr / NSAMPS - E*E; 81 SE = sqrt(V/NSAMPS); 82 CHK(eq_eps(E, g, 3*SE)); 83 } 84 85 ssp_ran_sphere_uniform(rng, wo, NULL); 86 ssf_phase_sample(phase, rng, wo, wi, NULL); 87 CHK(d3_is_normalized(wi)); 88 89 CHK(ssf_phase_ref_put(phase) == RES_OK); 90 CHK(ssf_phase_ref_put(dummy) == RES_OK); 91 CHK(ssp_rng_ref_put(rng) == RES_OK); 92 93 check_memory_allocator(&allocator); 94 mem_shutdown_proxy_allocator(&allocator); 95 CHK(mem_allocated_size() == 0); 96 return 0; 97 } 98