test_ssf_phase_discrete.c (7961B)
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 #include <star/ssp.h> 20 21 /* Pre-normalized discrete phase function */ 22 static const struct ssf_discrete_item g_items[] = { 23 {0.0000000000000000, 68.603348420512177}, 24 {1.7453292519943295e-2, 50.735150232030207}, 25 {3.4906585039886591e-2, 24.908999683188814}, 26 {5.2359877559829883e-2, 11.936459532538100}, 27 {6.9813170079773182e-2, 6.5845291580807999}, 28 {8.7266462599716474e-2, 4.1115681738409000}, 29 {0.10471975511965977, 2.7971255874985190}, 30 {0.12217304763960307, 2.0200365258652173}, 31 {0.13962634015954636, 1.5241311066143122}, 32 {0.15707963267948966, 1.1891144970876193}, 33 {0.17453292519943295, 0.95259719392666198}, 34 {0.26179938779914941, 0.40729722471207858}, 35 {0.34906585039886590, 0.22439173045102234}, 36 {0.52359877559829882, 9.8911118269091880e-2}, 37 {0.69813170079773179, 5.7013070175769544e-2}, 38 {0.87266462599716477, 3.8428151500227763e-2}, 39 {1.0471975511965976, 2.8819207088581208e-2}, 40 {1.2217304763960306, 2.3435147759515759e-2}, 41 {1.3962634015954636, 2.0323680045268303e-2}, 42 {1.5707963267948966, 1.8554414090107985e-2}, 43 {1.7453292519943295, 1.7631650380735577e-2}, 44 {1.9198621771937625, 1.7273221501888446e-2}, 45 {2.0943951023931953, 1.7265595355529996e-2}, 46 {2.2689280275926285, 1.7471501307208134e-2}, 47 {2.4434609527920612, 1.7768921015187674e-2}, 48 {2.6179938779914944, 1.8089219162242556e-2}, 49 {2.7925268031909272, 1.8356134284788293e-2}, 50 {2.8797932657906435, 1.8455274187448138e-2}, 51 {2.9670597283903604, 1.8531535651032636e-2}, 52 {2.9845130209103035, 1.8546787943749535e-2}, 53 {3.0019663134302466, 1.8554414090107985e-2}, 54 {3.0194196059501901, 1.8562040236466435e-2}, 55 {3.0368728984701332, 1.8569666382824885e-2}, 56 {3.0543261909900763, 1.8584918675541785e-2}, 57 {3.0717794835100198, 1.8638301700050933e-2}, 58 {3.0892327760299634, 1.9156879652425504e-2}, 59 {3.1066860685499069, 2.1879413902392029e-2}, 60 {3.1241393610698500, 2.8323507575281980e-2}, 61 {3.1415926535897931, 3.2632280267806027e-2} 62 }; 63 static const size_t g_nitems = sizeof(g_items) / sizeof(*g_items); 64 65 struct context { 66 const struct ssf_discrete_item* items; 67 size_t nitems; 68 }; 69 70 /******************************************************************************* 71 * Helper functions 72 ******************************************************************************/ 73 static void 74 get_item(const size_t id, struct ssf_discrete_item* item, void* context) 75 { 76 const struct context* ctx = context; 77 CHK(id < ctx->nitems); 78 *item = ctx->items[id]; 79 } 80 81 static void 82 test_setup(struct ssf_phase* discrete) 83 { 84 struct ssf_discrete_setup_args args = SSF_DISCRETE_SETUP_ARGS_NULL; 85 struct context ctx; 86 struct ssf_phase* dummy = NULL; 87 ASSERT(discrete); 88 89 ctx.items = g_items; 90 ctx.nitems = g_nitems; 91 92 args.get_item = get_item; 93 args.context = &ctx; 94 args.nitems = g_nitems; 95 96 CHK(ssf_phase_discrete_setup(NULL, &args) == RES_BAD_ARG); 97 CHK(ssf_phase_discrete_setup(discrete, NULL) == RES_BAD_ARG); 98 CHK(ssf_phase_discrete_setup(discrete, &args) == RES_OK); 99 100 /* Invalid phase function type */ 101 CHK(ssf_phase_create(&mem_default_allocator, &phase_dummy, &dummy) == RES_OK); 102 CHK(ssf_phase_discrete_setup(dummy, &args) == RES_BAD_ARG); 103 CHK(ssf_phase_ref_put(dummy) == RES_OK); 104 105 /* Invalid last #items */ 106 ctx.nitems = args.nitems = 1; 107 CHK(ssf_phase_discrete_setup(discrete, &args) == RES_BAD_ARG); 108 109 /* Last angle is not PI */ 110 ctx.nitems = args.nitems = g_nitems - 1; 111 CHK(ssf_phase_discrete_setup(discrete, &args) == RES_BAD_ARG); 112 113 /* First angle is not 0 */ 114 ctx.items = g_items + 1; 115 CHK(ssf_phase_discrete_setup(discrete, &args) == RES_BAD_ARG); 116 } 117 118 static void 119 test_eval(struct ssf_phase* discrete, struct ssp_rng* rng) 120 { 121 struct ssf_discrete_setup_args args = SSF_DISCRETE_SETUP_ARGS_NULL; 122 struct context ctx; 123 double wo[3]; 124 double wi[3]; 125 size_t i; 126 127 ctx.items = g_items; 128 ctx.nitems = g_nitems; 129 args.get_item = get_item; 130 args.context = &ctx; 131 args.nitems = g_nitems; 132 CHK(ssf_phase_discrete_setup(discrete, &args) == RES_OK); 133 134 d3(wo, 0, 0,-1); 135 d3(wi, 0, 0, 1); 136 CHK(ssf_phase_eval(discrete, wo, wi) == g_items[0].value); 137 d3(wo, 0, 0, 1); 138 CHK(ssf_phase_eval(discrete, wo, wi) == g_items[g_nitems-1].value); 139 140 FOR_EACH(i, 0, 10) { 141 double cos_theta; 142 double val; 143 double ref; 144 size_t iitem; 145 146 ssp_ran_sphere_uniform(rng, wo, NULL); 147 ssp_ran_sphere_uniform(rng, wi, NULL); 148 val = ssf_phase_eval(discrete, wo, wi); 149 150 cos_theta = d3_dot(d3_minus(wo, wo), wi); 151 if(cos_theta == 0) { 152 ref = g_items[0].value; 153 } else { 154 const double theta = acos(cos_theta); 155 double u; 156 /* Look for the phase function discrete items to consider regarding the 157 * sampled wo and wi directions */ 158 FOR_EACH(iitem, 0, g_nitems) { 159 if(g_items[iitem].theta >= theta) break; 160 } 161 ASSERT(iitem < g_nitems && iitem > 0); 162 163 /* Compute the parameter of the linear interpolation */ 164 u = cos_theta - cos(g_items[iitem-1].theta); 165 u /= cos(g_items[iitem].theta) - cos(g_items[iitem-1].theta); 166 167 ref = g_items[iitem-1].value + u*(g_items[iitem].value - g_items[iitem-1].value); 168 } 169 CHK(eq_eps(val, ref, 1.e-6)); 170 } 171 } 172 173 static void 174 test_sample(struct ssf_phase* discrete, struct ssp_rng* rng) 175 { 176 struct ssf_discrete_setup_args args = SSF_DISCRETE_SETUP_ARGS_NULL; 177 struct context ctx; 178 double ref = 0; 179 size_t iitem; 180 size_t i; 181 182 ctx.items = g_items; 183 ctx.nitems = g_nitems; 184 args.get_item = get_item; 185 args.context = &ctx; 186 args.nitems = g_nitems; 187 CHK(ssf_phase_discrete_setup(discrete, &args) == RES_OK); 188 189 FOR_EACH(iitem, 1, g_nitems) { 190 const double mu0 = cos(g_items[iitem-1].theta); 191 const double mu1 = cos(g_items[iitem-0].theta); 192 const double phi0 = g_items[iitem-1].value; 193 const double phi1 = g_items[iitem-0].value; 194 const double delta_mu = mu0 - mu1; 195 ref += (mu0*phi0 + mu1*phi1)/2 * delta_mu; 196 } 197 ref *= 2*PI; 198 199 FOR_EACH(i, 0, 10) { 200 const size_t N = 10000; 201 double wo[3]; 202 double wi[3]; 203 double sum = 0; 204 double sum2 = 0; 205 double E = 0; 206 double V = 0; 207 double SE = 0; 208 size_t ireal; 209 ssp_ran_sphere_uniform(rng, wo, NULL); 210 211 FOR_EACH(ireal, 0, N) { 212 double w[3]; 213 double mu; 214 double pdf; 215 ssf_phase_sample(discrete, rng, wo, wi, &pdf); 216 217 mu = d3_dot(d3_minus(w, wo), wi); 218 CHK(pdf == ssf_phase_eval(discrete, wo, wi)); 219 220 sum += mu; 221 sum2 += mu*mu; 222 } 223 E = sum / (double)N; 224 V = sum2 / (double)N - E*E; 225 SE = sqrt(V/(double)N); 226 CHK(eq_eps(E, ref, 3*SE)); 227 } 228 } 229 230 /******************************************************************************* 231 * Main function 232 ******************************************************************************/ 233 int 234 main(int argc, char** argv) 235 { 236 struct ssp_rng* rng = NULL; 237 struct ssf_phase* discrete = NULL; 238 (void)argc, (void)argv; 239 240 CHK(ssf_phase_create 241 (&mem_default_allocator, &ssf_phase_discrete, &discrete) == RES_OK); 242 CHK(ssp_rng_create(&mem_default_allocator, SSP_RNG_MT19937_64, &rng) == RES_OK); 243 244 test_setup(discrete); 245 test_eval(discrete, rng); 246 test_sample(discrete, rng); 247 248 CHK(ssf_phase_ref_put(discrete) == RES_OK); 249 CHK(ssp_rng_ref_put(rng) == RES_OK); 250 return 0; 251 }