ssf_phase_discrete.c (10183B)
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 /* nextafter */ 17 18 #include "ssf.h" 19 #include "ssf_phase_c.h" 20 21 #include <star/ssp.h> 22 23 #include <rsys/algorithm.h> 24 #include <rsys/double3.h> 25 #include <rsys/dynamic_array_double.h> 26 27 #include <math.h> /* nextafter */ 28 29 struct discrete_item { 30 double cos_theta; 31 double value; 32 }; 33 34 /* Define the dynamic array of discrete items */ 35 #define DARRAY_NAME discrete_item 36 #define DARRAY_DATA struct discrete_item 37 #include <rsys/dynamic_array.h> 38 39 struct discrete { 40 struct darray_discrete_item items; 41 struct darray_double cumulative; 42 }; 43 44 /******************************************************************************* 45 * Helper functions 46 ******************************************************************************/ 47 static INLINE res_T 48 check_ssf_discrete_setup_args 49 (const struct ssf_discrete_setup_args* args) 50 { 51 struct ssf_discrete_item item; 52 if(!args || args->nitems < 2 || !args->get_item) 53 return RES_BAD_ARG; 54 55 /* Invalid 1st angle */ 56 args->get_item(0, &item, args->context); 57 if(item.theta != 0) 58 return RES_BAD_ARG; 59 60 /* Invalid last angle */ 61 args->get_item(args->nitems-1, &item, args->context); 62 if(!eq_eps(item.theta, PI, 1.e-4)) 63 return RES_BAD_ARG; 64 65 return RES_OK; 66 } 67 68 static res_T 69 setup_discrete_items 70 (struct discrete* discrete, 71 const struct ssf_discrete_setup_args* args) 72 { 73 struct discrete_item* items = NULL; 74 size_t i = 0; 75 res_T res = RES_OK; 76 ASSERT(discrete && args); 77 78 res = darray_discrete_item_resize(&discrete->items, args->nitems); 79 if(res != RES_OK) goto error; 80 items = darray_discrete_item_data_get(&discrete->items); 81 82 FOR_EACH(i, 0, args->nitems) { 83 struct ssf_discrete_item item; 84 args->get_item(i, &item, args->context); 85 items[i].cos_theta = cos(item.theta); 86 items[i].value = item.value; 87 88 /* Angles must be sorted in ascending order */ 89 if(i > 0 && items[i].cos_theta > items[i-1].cos_theta) { 90 res = RES_BAD_ARG; 91 goto error; 92 } 93 } 94 95 exit: 96 return res; 97 error: 98 goto exit; 99 } 100 101 static res_T 102 compute_cumulative(struct discrete* discrete) 103 { 104 struct discrete_item* items = NULL; 105 double* cumul = NULL; 106 double alpha = 0; 107 size_t n; 108 size_t i; 109 res_T res = RES_OK; 110 ASSERT(discrete); 111 112 n = darray_discrete_item_size_get(&discrete->items); 113 ASSERT(n >= 2); 114 115 res = darray_double_resize(&discrete->cumulative, n); 116 if(res != RES_OK) goto error; 117 118 items = darray_discrete_item_data_get(&discrete->items); 119 cumul = darray_double_data_get(&discrete->cumulative); 120 121 alpha = 0; 122 cumul[0] = alpha; 123 FOR_EACH(i, 1, n) { 124 alpha += PI 125 * (items[i-1].value + items[i].value) 126 * (items[i-1].cos_theta - items[i].cos_theta); 127 128 cumul[i] = alpha; 129 } 130 131 /* Ensure that the phase function is normalized */ 132 if(!eq_eps(alpha, 1, 1.e-6)) { 133 const double rcp_alpha = 1.0 / alpha; 134 FOR_EACH(i, 0, n) { 135 items[i].value *= rcp_alpha; 136 cumul[i] *= rcp_alpha; 137 } 138 } 139 140 exit: 141 return res; 142 error: 143 goto exit; 144 } 145 146 static INLINE int 147 cmp_discrete_item(const void* a, const void* b) 148 { 149 const struct discrete_item* item; 150 double key; 151 ASSERT(a && b); 152 153 key = *((double*)a); 154 item = b; 155 156 /* The array is sorted in ascending order of the angles while it stores the 157 * cosine of the angles. In the following, we reverse the test to find the 158 * first entry whose _angle_ is not less than the angle whose cosine is the 159 * key sought. This ensures that the array is sorted in ascending order as 160 * expected by the search_lower_bound routine */ 161 if(key < item->cos_theta) { 162 return +1; 163 } else if(key > item->cos_theta) { 164 return -1; 165 } else { 166 return 0; 167 } 168 } 169 170 static INLINE int 171 cmp_cumul(const void* a, const void* b) 172 { 173 double key; 174 double cumul; 175 ASSERT(a && b); 176 177 key = *((double*)a); 178 cumul = *((double*)b); 179 180 if(key < cumul) { 181 return -1; 182 } else if(key > cumul) { 183 return +1; 184 } else { 185 return 0; 186 } 187 } 188 189 /******************************************************************************* 190 * Private functions 191 ******************************************************************************/ 192 static res_T 193 discrete_init(struct mem_allocator* allocator, void* phase) 194 { 195 struct discrete* discrete = phase; 196 ASSERT(phase); 197 darray_discrete_item_init(allocator, &discrete->items); 198 darray_double_init(allocator, &discrete->cumulative); 199 return RES_OK; 200 } 201 202 static void 203 discrete_release(void* phase) 204 { 205 struct discrete* discrete = phase; 206 ASSERT(phase); 207 darray_discrete_item_release(&discrete->items); 208 darray_double_release(&discrete->cumulative); 209 } 210 211 static double 212 discrete_eval(void* phase, const double wo[3], const double wi[3]) 213 { 214 const struct discrete* discrete = phase; 215 const struct discrete_item* items = NULL; 216 double v[3]; 217 double cos_theta = 0; 218 double value = 0; 219 size_t nitems = 0; 220 ASSERT(phase && wo && wi); 221 ASSERT(d3_is_normalized(wo) && d3_is_normalized(wi)); 222 223 d3_minus(v, wo); 224 cos_theta = d3_dot(v, wi); 225 226 items = darray_discrete_item_cdata_get(&discrete->items); 227 nitems = darray_discrete_item_size_get(&discrete->items); 228 ASSERT(nitems >= 2); 229 230 if(eq_eps(cos_theta, 1, 1.e-6)) { 231 value = items[0].value; 232 } else if(eq_eps(cos_theta, 0, 1.e-6)) { 233 value = items[nitems-1].value; 234 } else { 235 const struct discrete_item* found_item = NULL; 236 size_t iitem = 0; 237 double u = 0; 238 239 /* Search for the discrete item whose angle is greater or equal to the 240 * angle between wo and wi */ 241 found_item = search_lower_bound 242 (&cos_theta, items, nitems, sizeof(*items), cmp_discrete_item); 243 iitem = (size_t)(found_item - items); 244 ASSERT(found_item && iitem > 0 && iitem < nitems); 245 ASSERT(cos_theta < items[iitem-1].cos_theta); 246 ASSERT(cos_theta >= items[iitem].cos_theta); 247 248 /* Linearly interpolate the phase function value */ 249 u = 250 (cos_theta - items[iitem].cos_theta) 251 / (items[iitem-1].cos_theta - items[iitem].cos_theta); 252 value = items[iitem].value + u*(items[iitem-1].value - items[iitem].value); 253 } 254 255 return value; 256 } 257 258 static void 259 discrete_sample 260 (void* phase, 261 struct ssp_rng* rng, 262 const double wo[3], 263 double wi[3], 264 double *pdf) 265 { 266 const struct discrete* discrete = phase; 267 const struct discrete_item* items = NULL; 268 const double* cumul = NULL; 269 double frame[9]; 270 double w[3]; 271 double cos_theta = 0; 272 double sin_theta = 0; 273 double phi = 0; 274 size_t ncumul = 0; 275 size_t i = 0; 276 double r = 0; 277 ASSERT(phase && rng && wo && wi); 278 ASSERT(d3_is_normalized(wo)); 279 280 items = darray_discrete_item_cdata_get(&discrete->items); 281 cumul = darray_double_cdata_get(&discrete->cumulative); 282 ncumul = darray_double_size_get(&discrete->cumulative); 283 ASSERT(ncumul == darray_discrete_item_size_get(&discrete->items)); 284 285 /* Sample r in [0, 1] */ 286 r = ssp_rng_uniform_double(rng, 0, nextafter(1, 2)); 287 if(r == 0) { 288 cos_theta = 1; 289 sin_theta = 0; 290 } else if(r ==1) { 291 cos_theta =-1; 292 sin_theta = 0; 293 } else { 294 double* found_cumul = NULL; 295 double u; 296 297 /* Search for the first entry in the cumulative that is greater than or equal 298 * to the sampled number in [0, 1] */ 299 found_cumul = search_lower_bound(&r, cumul, ncumul, sizeof(*cumul), cmp_cumul); 300 301 i = (size_t)(found_cumul - cumul); 302 ASSERT(found_cumul && i > 0 && i < ncumul); 303 304 /* Linearly interpolate the sampled cos_theta */ 305 u = (r - cumul[i-1]) / (cumul[i] - cumul[i-1]); 306 ASSERT(r >= cumul[i-1] && r < cumul[i]); 307 cos_theta = items[i-1].cos_theta + u*(items[i].cos_theta - items[i-1].cos_theta); 308 309 /* Compute the sinus of theta from its cosine */ 310 sin_theta = sqrt(1 - cos_theta*cos_theta); 311 } 312 313 /* Uniformly sample phi in [0, 2PI[ */ 314 phi = ssp_rng_uniform_double(rng, 0, 2*PI); 315 316 /* Calculate the Cartesian coordinates of the direction in the local 317 * coordinate system of the phase function */ 318 wi[0] = cos(phi) * sin_theta; 319 wi[1] = sin(phi) * sin_theta; 320 wi[2] = cos_theta; 321 322 /* Calculate the transformation matrix from the local coordinate system of 323 * the phase function to the absolute coordinate system. Note that by 324 * convention in Star-SF the directions point outward from the scattering 325 * position. That's why we reverse 'wo' to point inside the scattering 326 * position */ 327 d33_basis(frame, d3_minus(w, wo)); 328 d33_muld3(wi, frame, wi); 329 ASSERT(eq_eps(d3_dot(wi, w), cos_theta, fabs(cos_theta*1.e-6))); 330 331 if(pdf) *pdf = discrete_eval(phase, wo, wi); 332 } 333 334 /******************************************************************************* 335 * Exported symbols 336 ******************************************************************************/ 337 const struct ssf_phase_type ssf_phase_discrete = { 338 discrete_init, 339 discrete_release, 340 discrete_sample, 341 discrete_eval, 342 discrete_eval, 343 sizeof(struct discrete), 344 ALIGNOF(struct discrete) 345 }; 346 347 res_T 348 ssf_phase_discrete_setup 349 (struct ssf_phase* phase, 350 const struct ssf_discrete_setup_args* args) 351 { 352 struct discrete* discrete = NULL; 353 res_T res = RES_OK; 354 355 if(!phase || !PHASE_TYPE_EQ(&phase->type, &ssf_phase_discrete)) { 356 res = RES_BAD_ARG; 357 goto error; 358 } 359 res = check_ssf_discrete_setup_args(args); 360 if(res != RES_OK) goto error; 361 362 discrete = phase->data; 363 364 res = setup_discrete_items(discrete, args); 365 if(res != RES_OK) goto error; 366 res = compute_cumulative(discrete); 367 if(res != RES_OK) goto error; 368 369 exit: 370 return res; 371 error: 372 if(discrete) { 373 darray_discrete_item_purge(&discrete->items); 374 darray_double_purge(&discrete->cumulative); 375 } 376 goto exit; 377 }