ssp.h (33290B)
1 /* Copyright (C) 2015-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 #ifndef SSP_H 17 #define SSP_H 18 19 #include <rsys/double33.h> 20 #include <rsys/float33.h> 21 #include <rsys/math.h> 22 #include <rsys/rsys.h> 23 24 #include <stdint.h> /* SIZE_MAX */ 25 26 /* Library symbol management */ 27 #if defined(SSP_SHARED_BUILD) /* Build shared library */ 28 #define SSP_API extern EXPORT_SYM 29 #elif defined(SSP_STATIC) /* Use/build static library */ 30 #define SSP_API extern LOCAL_SYM 31 #else /* Use shared library */ 32 #define SSP_API extern IMPORT_SYM 33 #endif 34 35 /* Helper macro that asserts if the invocation of the ssp function `Func' 36 * returns an error. One should use this macro on smc function calls for which 37 * no explicit error checking is performed */ 38 #ifndef NDEBUG 39 #define SSP(Func) ASSERT(ssp_ ## Func == RES_OK) 40 #else 41 #define SSP(Func) ssp_ ## Func 42 #endif 43 44 #define SSP_SEQUENCE_ID_NONE SIZE_MAX 45 46 /* Forward declaration of opaque types */ 47 struct ssp_rng; 48 struct ssp_rng_proxy; 49 struct ssp_ranst_discrete; 50 struct ssp_ranst_gaussian; 51 struct ssp_ranst_gaussian_float; 52 struct ssp_ranst_piecewise_linear; 53 struct ssp_ranst_piecewise_linear_float; 54 55 /* Forward declaration of external types */ 56 struct mem_allocator; 57 58 enum ssp_rng_type { 59 /* David Jones's Keep It Simple Stupid builtin PRNG type. Suitable for fast 60 * basic randomness */ 61 SSP_RNG_KISS, 62 /* 64-bits Mersenne Twister builtin PRNG type of Matsumoto and Nishimura */ 63 SSP_RNG_MT19937_64, 64 /* 48-bits RANLUX builtin PRNG type of Lusher and James, 1994 */ 65 SSP_RNG_RANLUX48, 66 /* A random_device RNG is a uniformly-distributed integer random number 67 * generator that produces non-deterministic random numbers. It may be 68 * implemented in terms of an implementation-defined pseudo-random number 69 * engine if a non-deterministic source (e.g. a hardware device) is not 70 * available to the implementation. In this case each random_device object 71 * may generate the same number sequence. */ 72 SSP_RNG_RANDOM_DEVICE, 73 /* Counter Based RNG threefry of Salmon, Moraes, Dror & Shaw */ 74 SSP_RNG_THREEFRY, 75 /* Counter Based RNG aes of Salmon, Moraes, Dror & Shaw */ 76 SSP_RNG_AES, 77 SSP_RNG_TYPES_COUNT__, 78 SSP_RNG_TYPE_NULL = SSP_RNG_TYPES_COUNT__ 79 }; 80 81 /* Arguments of the ssp_rng_proxy_create2 function */ 82 struct ssp_rng_proxy_create2_args { 83 const struct ssp_rng* rng; /* Original RNG state. May be NULL*/ 84 enum ssp_rng_type type; /* Type of the RN if 'rng' is NULL */ 85 size_t sequence_offset; /* #RNs before the 1st valid sequence */ 86 size_t sequence_size; /* #RNs in a sequence */ 87 size_t sequence_pitch; /* #RNs between sequences. Must be >= sequence_size */ 88 size_t nbuckets; /* #buckets of continuous RNs in a sequence */ 89 }; 90 #define SSP_RNG_PROXY_CREATE2_ARGS_NULL__ {NULL, SSP_RNG_TYPE_NULL, 0, 0, 0, 0} 91 static const struct ssp_rng_proxy_create2_args SSP_RNG_PROXY_CREATE2_ARGS_NULL = 92 SSP_RNG_PROXY_CREATE2_ARGS_NULL__; 93 94 BEGIN_DECLS 95 96 /******************************************************************************* 97 * Random Number Generator API 98 ******************************************************************************/ 99 SSP_API res_T 100 ssp_rng_create 101 (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ 102 const enum ssp_rng_type type, 103 struct ssp_rng** rng); 104 105 SSP_API res_T 106 ssp_rng_ref_put 107 (struct ssp_rng* rng); 108 109 SSP_API res_T 110 ssp_rng_ref_get 111 (struct ssp_rng* rng); 112 113 SSP_API res_T 114 ssp_rng_get_type 115 (const struct ssp_rng* rng, 116 enum ssp_rng_type* type); 117 118 SSP_API res_T 119 ssp_rng_set 120 (struct ssp_rng* rng, 121 const uint64_t seed); 122 123 /* Uniform random distribution in [ssp_rng_min(rng), ssp_rng_max(rng)] */ 124 SSP_API uint64_t 125 ssp_rng_get 126 (struct ssp_rng* rng); 127 128 /* Advances the internal state by n times. 129 Equivalent to calling ssp_rng_get() n times and discarding the result */ 130 SSP_API res_T 131 ssp_rng_discard 132 (struct ssp_rng* rng, uint64_t n); 133 134 /* Uniform random integer distribution in [lower, upper] */ 135 SSP_API uint64_t 136 ssp_rng_uniform_uint64 137 (struct ssp_rng* rng, 138 const uint64_t lower, 139 const uint64_t upper); 140 141 /* Uniform random floating point distribution in [lower, upper) */ 142 SSP_API double 143 ssp_rng_uniform_double 144 (struct ssp_rng* rng, 145 const double lower, 146 const double upper); 147 148 SSP_API float 149 ssp_rng_uniform_float 150 (struct ssp_rng* rng, 151 const float lower, 152 const float upper); 153 154 /* Uniform random floating point distribution in [0, 1) */ 155 SSP_API double 156 ssp_rng_canonical 157 (struct ssp_rng* rng); 158 159 /* Uniform random single precision floating point distribution in [0, 1) */ 160 SSP_API float 161 ssp_rng_canonical_float 162 (struct ssp_rng* rng); 163 164 SSP_API uint64_t 165 ssp_rng_min 166 (struct ssp_rng* rng); 167 168 SSP_API uint64_t 169 ssp_rng_max 170 (struct ssp_rng* rng); 171 172 SSP_API res_T 173 ssp_rng_read 174 (struct ssp_rng* rng, 175 FILE* stream); 176 177 SSP_API res_T 178 ssp_rng_read_cstr 179 (struct ssp_rng* rng, 180 const char* cstr); /* Null terminated string */ 181 182 SSP_API res_T 183 ssp_rng_write 184 (const struct ssp_rng* rng, 185 FILE* stream); 186 187 SSP_API res_T 188 ssp_rng_write_cstr 189 (const struct ssp_rng* rng, 190 char* buf, /* May be NULL */ 191 const size_t bufsz, /* buf capacity */ 192 size_t* len); /* May be NULL. #chars to write into buf without null char */ 193 194 SSP_API double 195 ssp_rng_entropy 196 (const struct ssp_rng* rng); 197 198 /******************************************************************************* 199 * Proxy of Random Number Generators - It manages a list of `nbuckets' RNGs 200 * whose each have independant `infinite' random sequences 201 ******************************************************************************/ 202 SSP_API res_T 203 ssp_rng_proxy_create 204 (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ 205 const enum ssp_rng_type type, 206 const size_t nbuckets, 207 struct ssp_rng_proxy** proxy); 208 209 /* Build the proxy RNG from the state of `rng'. Note that `rng' is read only 210 * argument, i.e. it is not used internally. */ 211 SSP_API res_T 212 ssp_rng_proxy_create_from_rng 213 (struct mem_allocator* allocator, 214 const struct ssp_rng* rng, 215 const size_t nbuckets, 216 struct ssp_rng_proxy** proxy); 217 218 SSP_API res_T 219 ssp_rng_proxy_create2 220 (struct mem_allocator* mem_allocator, 221 const struct ssp_rng_proxy_create2_args* args, 222 struct ssp_rng_proxy** out_proxy); 223 224 SSP_API res_T 225 ssp_rng_proxy_read 226 (struct ssp_rng_proxy* proxy, 227 FILE* stream); 228 229 SSP_API res_T 230 ssp_rng_proxy_write 231 (const struct ssp_rng_proxy* proxy, 232 FILE* stream); 233 234 SSP_API res_T 235 ssp_rng_proxy_ref_get 236 (struct ssp_rng_proxy* proxy); 237 238 SSP_API res_T 239 ssp_rng_proxy_ref_put 240 (struct ssp_rng_proxy* proxy); 241 242 /* Create the RNG of `ibucket'. Return a RES_BAD_ARG error if RNG was already 243 * created for `ibucket'. Call ssp_rng_ref_put to release the returned RNG */ 244 SSP_API res_T 245 ssp_rng_proxy_create_rng 246 (struct ssp_rng_proxy* proxy, 247 const size_t ibucket, 248 struct ssp_rng** rng); 249 250 SSP_API res_T 251 ssp_rng_proxy_get_type 252 (const struct ssp_rng_proxy* proxy, 253 enum ssp_rng_type* type); 254 255 /* Returns the index of the current sequence. This index is independent of the 256 * original seed used by the proxy. When creating the proxy, the sequence index 257 * is 0. It is then incremented by one each time a new sequence is required. 258 * This index is used to identify the state of the proxy relative to the 259 * original seed. */ 260 SSP_API res_T 261 ssp_rng_proxy_get_sequence_id 262 (const struct ssp_rng_proxy* proxy, 263 size_t* sequence_id); 264 265 /* Discard the random numbers of `n' sequences. If `n' is 0, nothing happens. 266 * If `n' is greater than or equal to one, the random numbers of the current 267 * sequence are ignored as well as the random numbers of the following 268 * sequences until the sequence identifier has been incremented by `n-1'. */ 269 SSP_API res_T 270 ssp_rng_proxy_flush_sequences 271 (struct ssp_rng_proxy* proxy, 272 const size_t n); 273 274 /******************************************************************************* 275 * Miscellaneous distributions 276 ******************************************************************************/ 277 /* Random variate from the exponential distribution with mean `1/mu': 278 * P(x) dx = mu * exp(-x * mu) dx */ 279 SSP_API double 280 ssp_ran_exp 281 (struct ssp_rng* rng, 282 const double mu); 283 284 SSP_API float 285 ssp_ran_exp_float 286 (struct ssp_rng* rng, 287 const float mu); 288 289 SSP_API double 290 ssp_ran_exp_pdf 291 (const double x, 292 const double mu); 293 294 SSP_API float 295 ssp_ran_exp_float_pdf 296 (const float x, 297 const float mu); 298 299 /* Truncated exponential distribution */ 300 SSP_API double 301 ssp_ran_exp_truncated 302 (struct ssp_rng* rng, 303 const double mu, 304 const double max); 305 306 SSP_API float 307 ssp_ran_exp_truncated_float 308 (struct ssp_rng* rng, 309 const float mu, 310 const float max); 311 312 SSP_API double 313 ssp_ran_exp_truncated_pdf 314 (const double x, 315 const double mu, 316 const double max); 317 318 SSP_API float 319 ssp_ran_exp_truncated_float_pdf 320 (const float x, 321 const float mu, 322 const float max); 323 324 /* Stateless random variate from the gaussian (or normal) distribution 325 * with mean `mu': 326 * P(x) dx = 1 / (sigma*sqrt(2*PI)) * exp(-1/2*((x-mu)/sigma)^2) dx */ 327 SSP_API double 328 ssp_ran_gaussian 329 (struct ssp_rng* rng, 330 const double mu, 331 const double sigma); 332 333 SSP_API float 334 ssp_ran_gaussian_float 335 (struct ssp_rng* rng, 336 const float mu, 337 const float sigma); 338 339 /* Return the probability distribution for a gaussian random variate */ 340 SSP_API double 341 ssp_ran_gaussian_pdf 342 (const double x, 343 const double mu, 344 const double sigma); 345 346 SSP_API float 347 ssp_ran_gaussian_float_pdf 348 (const float x, 349 const float mu, 350 const float sigma); 351 352 /* Random variate from the lognormal distribution: 353 * P(x) dx = 1/(sigma*x*sqrt(2*PI)) * exp(-(ln(x)-zeta)^2 / (2*sigma^2)) dx */ 354 SSP_API double 355 ssp_ran_lognormal 356 (struct ssp_rng* rng, 357 const double zeta, 358 const double sigma); 359 360 SSP_API float 361 ssp_ran_lognormal_float 362 (struct ssp_rng* rng, 363 const float zeta, 364 const float sigma); 365 366 /* Return the probability distribution for a lognormal random variate */ 367 SSP_API double 368 ssp_ran_lognormal_pdf 369 (const double x, 370 const double zeta, 371 const double sigma); 372 373 SSP_API float 374 ssp_ran_lognormal_float_pdf 375 (const float x, 376 const float zeta, 377 const float sigma); 378 379 /******************************************************************************* 380 * Sphere distribution 381 ******************************************************************************/ 382 /* Uniform sampling of an unit sphere. The PDF of the generated sample is 383 * stored in sample[3] */ 384 SSP_API double* 385 ssp_ran_sphere_uniform 386 (struct ssp_rng* rng, 387 double sample[3], 388 double* pdf); /* Can be NULL => pdf not computed */ 389 390 SSP_API float* 391 ssp_ran_sphere_uniform_float 392 (struct ssp_rng* rng, 393 float sample[3], 394 float* pdf); /* Can be NULL => pdf not computed */ 395 396 /* Return the probability distribution for a sphere uniform random variate */ 397 static INLINE double 398 ssp_ran_sphere_uniform_pdf(void) 399 { 400 return 1 / (4*PI); 401 } 402 403 static INLINE float 404 ssp_ran_sphere_uniform_float_pdf(void) 405 { 406 return 1 / (4*(float)PI); 407 } 408 409 /******************************************************************************* 410 * Circle distribution 411 ******************************************************************************/ 412 SSP_API double* 413 ssp_ran_circle_uniform 414 (struct ssp_rng* rng, 415 double sample[2], 416 double* pdf); /* Can be NULL => pdf not computed */ 417 418 SSP_API float* 419 ssp_ran_circle_uniform_float 420 (struct ssp_rng* rng, 421 float sample[2], 422 float* pdf); /* Can be NULL => pdf not computed */ 423 424 static INLINE double 425 ssp_ran_circle_uniform_pdf(void) 426 { 427 return 1/(2*PI); 428 } 429 430 static INLINE float 431 ssp_ran_circle_uniform_float_pdf(void) 432 { 433 return 1/(2*(float)PI); 434 } 435 436 /******************************************************************************* 437 * Triangle distribution 438 ******************************************************************************/ 439 /* Uniform sampling of a triangle */ 440 SSP_API double* 441 ssp_ran_triangle_uniform 442 (struct ssp_rng* rng, 443 const double v0[3], /* Position of the 1st vertex */ 444 const double v1[3], /* Position of the 2nd vertex */ 445 const double v2[3], /* Position of the 3rd vertex */ 446 double sample[3], /* Sampled position */ 447 double* pdf); /* Can be NULL => pdf not computed */ 448 449 SSP_API float* 450 ssp_ran_triangle_uniform_float 451 (struct ssp_rng* rng, 452 const float v0[3], /* Position of the 1st vertex */ 453 const float v1[3], /* Position of the 2nd vertex */ 454 const float v2[3], /* Position of the 3rd vertex */ 455 float sample[3], /* Sampled position */ 456 float* pdf); /* Can be NULL => pdf not computed */ 457 458 /* Return the probability distribution for a uniform point sampling on a triangle */ 459 static INLINE double 460 ssp_ran_triangle_uniform_pdf 461 (const double v0[3], 462 const double v1[3], 463 const double v2[3]) 464 { 465 double vec0[3], vec1[3], tmp[3]; 466 ASSERT(v0 && v1 && v2); 467 468 d3_sub(vec0, v0, v2); 469 d3_sub(vec1, v1, v2); 470 return 2 / d3_len(d3_cross(tmp, vec0, vec1)); 471 } 472 473 static INLINE float 474 ssp_ran_triangle_uniform_float_pdf 475 (const float v0[3], 476 const float v1[3], 477 const float v2[3]) 478 { 479 float vec0[3], vec1[3], tmp[3]; 480 ASSERT(v0 && v1 && v2); 481 482 f3_sub(vec0, v0, v2); 483 f3_sub(vec1, v1, v2); 484 return 2 / f3_len(f3_cross(tmp, vec0, vec1)); 485 } 486 487 /******************************************************************************* 488 * Tetrahedron distribution 489 ******************************************************************************/ 490 /* Uniform sampling of a tetrahedron */ 491 SSP_API double* 492 ssp_ran_tetrahedron_uniform 493 (struct ssp_rng* rng, 494 const double v0[3], /* Position of the 1st vertex */ 495 const double v1[3], /* Position of the 2nd vertex */ 496 const double v2[3], /* Position of the 3rd vertex */ 497 const double v3[3], /* Position of the 4th vertex */ 498 double sample[3], /* Sampled position */ 499 double* pdf); /* Can be NULL => pdf not computed */ 500 501 SSP_API float* 502 ssp_ran_tetrahedron_uniform_float 503 (struct ssp_rng* rng, 504 const float v0[3], /* Position of the 1st vertex */ 505 const float v1[3], /* Position of the 2nd vertex */ 506 const float v2[3], /* Position of the 3rd vertex */ 507 const float v3[3], /* Position of the 4th vertex */ 508 float sample[3], /* Sampled position */ 509 float* pdf); /* Can be NULL => pdf not computed */ 510 511 /* Return the probability distribution for a uniform point sampling in a tetrahedron */ 512 static INLINE double 513 ssp_ran_tetrahedron_uniform_pdf 514 (const double v0[3], 515 const double v1[3], 516 const double v2[3], 517 const double v3[3]) 518 { 519 double vec0[3], vec1[3], vec2[3], tmp[3]; 520 ASSERT(v0 && v1 && v2 && v3); 521 522 d3_sub(vec0, v1, v2); 523 d3_sub(vec1, v3, v2); 524 d3_sub(vec2, v0, v2); 525 return 6 / fabs(d3_dot(d3_cross(tmp, vec0, vec1), vec2)); 526 } 527 528 static INLINE float 529 ssp_ran_tetrahedron_uniform_float_pdf 530 (const float v0[3], 531 const float v1[3], 532 const float v2[3], 533 const float v3[3]) 534 { 535 float vec0[3], vec1[3], vec2[3], tmp[3]; 536 ASSERT(v0 && v1 && v2 && v3); 537 538 f3_sub(vec0, v1, v2); 539 f3_sub(vec1, v3, v2); 540 f3_sub(vec2, v0, v2); 541 return 6.f / absf(f3_dot(f3_cross(tmp, vec0, vec1), vec2)); 542 } 543 544 /******************************************************************************* 545 * Hemisphere distribution 546 ******************************************************************************/ 547 /* Uniform sampling of an unit hemisphere whose up direction is implicitly the 548 * Z axis. */ 549 SSP_API double* 550 ssp_ran_hemisphere_uniform_local 551 (struct ssp_rng* rng, 552 double sample[3], 553 double* pdf); /* Can be NULL => pdf not computed */ 554 555 SSP_API float* 556 ssp_ran_hemisphere_uniform_float_local 557 (struct ssp_rng* rng, 558 float sample[3], 559 float* pdf); /* Can be NULL => pdf not computed */ 560 561 /* Return the probability distribution for an hemispheric uniform random 562 * variate with an implicit up direction in Z */ 563 static INLINE double 564 ssp_ran_hemisphere_uniform_local_pdf(const double sample[3]) 565 { 566 ASSERT(sample); 567 return sample[2] < 0 ? 0 : 1/(2*PI); 568 } 569 570 static INLINE float 571 ssp_ran_hemisphere_uniform_float_local_pdf(const float sample[3]) 572 { 573 ASSERT(sample); 574 return sample[2] < 0 ? 0 : 1/(2*(float)PI); 575 } 576 577 /* Uniform sampling of an unit hemisphere whose up direction is `up'. */ 578 static INLINE double* 579 ssp_ran_hemisphere_uniform 580 (struct ssp_rng* rng, 581 const double up[3], 582 double sample[3], 583 double* pdf) /* Can be NULL => pdf not computed */ 584 { 585 double basis[9]; 586 double sample_local[3]; 587 ASSERT(rng && up && sample && d3_is_normalized(up)); 588 ssp_ran_hemisphere_uniform_local(rng, sample_local, pdf); 589 return d33_muld3(sample, d33_basis(basis, up), sample_local); 590 } 591 592 static INLINE float* 593 ssp_ran_hemisphere_uniform_float 594 (struct ssp_rng* rng, 595 const float up[3], 596 float sample[3], 597 float* pdf) /* Can be NULL => pdf not computed */ 598 { 599 float basis[9]; 600 float sample_local[3]; 601 ASSERT(rng && up && sample && f3_is_normalized(up)); 602 ssp_ran_hemisphere_uniform_float_local(rng, sample_local, pdf); 603 return f33_mulf3(sample, f33_basis(basis, up), sample_local); 604 } 605 606 /* Return the probability distribution for an hemispheric uniform random 607 * variate with an explicit `up' direction */ 608 static INLINE double 609 ssp_ran_hemisphere_uniform_pdf 610 (const double up[3], 611 const double sample[3]) 612 { 613 ASSERT(up && sample && d3_is_normalized(up)); 614 return d3_dot(sample, up) < 0 ? 0 : 1/(2*PI); 615 } 616 617 static INLINE float 618 ssp_ran_hemisphere_uniform_float_pdf 619 (const float up[3], 620 const float sample[3]) 621 { 622 ASSERT(up && sample && f3_is_normalized(up)); 623 return f3_dot(sample, up) < 0 ? 0 : 1/(2*(float)PI); 624 } 625 626 /* Cosine weighted sampling of an unit hemisphere whose up direction is 627 * implicitly the Z axis. The PDF of the generated sample is stored in 628 * sample[3] */ 629 SSP_API double* 630 ssp_ran_hemisphere_cos_local 631 (struct ssp_rng* rng, 632 double sample[3], 633 double* pdf); /* Can be NULL => pdf not computed */ 634 635 SSP_API float* 636 ssp_ran_hemisphere_cos_float_local 637 (struct ssp_rng* rng, 638 float sample[3], 639 float* pdf); /* Can be NULL => pdf not computed */ 640 641 /* Return the probability distribution for an hemispheric cosine weighted 642 * random variate with an implicit up direction in Z */ 643 static INLINE double 644 ssp_ran_hemisphere_cos_local_pdf(const double sample[3]) 645 { 646 ASSERT(sample); 647 return sample[2] < 0 ? 0 : sample[2]/PI; 648 } 649 650 static INLINE float 651 ssp_ran_hemisphere_cos_float_local_pdf(const float sample[3]) 652 { 653 ASSERT(sample); 654 return sample[2] < 0 ? 0 : sample[2]/(float)PI; 655 } 656 657 /* Cosine weighted sampling of an unit hemisphere whose up direction is `up'. */ 658 static INLINE double* 659 ssp_ran_hemisphere_cos 660 (struct ssp_rng* rng, 661 const double up[3], 662 double sample[3], 663 double* pdf) /* Can be NULL => pdf not computed */ 664 { 665 double sample_local[3]; 666 double basis[9]; 667 ASSERT(rng && up && sample && d3_is_normalized(up)); 668 ssp_ran_hemisphere_cos_local(rng, sample_local, pdf); 669 return d33_muld3(sample, d33_basis(basis, up), sample_local); 670 } 671 672 static INLINE float* 673 ssp_ran_hemisphere_cos_float 674 (struct ssp_rng* rng, 675 const float up[3], 676 float sample[3], 677 float* pdf) /* Can be NULL => pdf not computed */ 678 { 679 float sample_local[3]; 680 float basis[9]; 681 ASSERT(rng && up && sample && f3_is_normalized(up)); 682 ssp_ran_hemisphere_cos_float_local(rng, sample_local, pdf); 683 return f33_mulf3(sample, f33_basis(basis, up), sample_local); 684 } 685 686 /* Return the probability distribution for an hemispheric cosine weighted 687 * random variate with an explicit `up' direction */ 688 static INLINE double 689 ssp_ran_hemisphere_cos_pdf 690 (const double up[3], 691 const double sample[3]) 692 { 693 double dot; 694 ASSERT(up && sample); 695 dot = d3_dot(sample, up); 696 return dot < 0 ? 0 : dot/PI; 697 } 698 699 static INLINE float 700 ssp_ran_hemisphere_cos_float_pdf 701 (const float up[3], 702 const float sample[3]) 703 { 704 float dot; 705 ASSERT(up && sample); 706 dot = f3_dot(sample, up); 707 return dot < 0 ? 0 : dot/(float)PI; 708 } 709 710 /******************************************************************************* 711 * Henyey Greenstein distribution 712 ******************************************************************************/ 713 /* Henyey-Greenstein sampling of an unit sphere for an incident direction 714 * that is implicitly the Z axis. */ 715 SSP_API double* 716 ssp_ran_sphere_hg_local 717 (struct ssp_rng* rng, 718 const double g, 719 double sample[3], 720 double* pdf); /* Can be NULL => pdf not computed */ 721 722 SSP_API float* 723 ssp_ran_sphere_hg_float_local 724 (struct ssp_rng* rng, 725 const float g, 726 float sample[3], 727 float* pdf); /* Can be NULL => pdf not computed */ 728 729 /* Return the probability distribution for a Henyey-Greenstein random 730 * variate with an implicit incident direction in Z */ 731 SSP_API double 732 ssp_ran_sphere_hg_local_pdf 733 (const double g, 734 const double sample[3]); 735 736 SSP_API float 737 ssp_ran_sphere_hg_float_local_pdf 738 (const float g, 739 const float sample[3]); 740 741 /* Henyey-Greenstein sampling of an unit sphere for an incident direction 742 * that is `up'. */ 743 static INLINE double* 744 ssp_ran_sphere_hg 745 (struct ssp_rng* rng, 746 const double up[3], 747 const double g, 748 double sample[3], 749 double* pdf) /* Can be NULL => pdf not computed */ 750 { 751 double sample_local[3]; 752 double basis[9]; 753 ASSERT(-1 <= g && g <= +1 && rng && up && sample && d3_is_normalized(up)); 754 if (fabs(g) == 1) { 755 d3_muld(sample, up, g); 756 if(pdf) *pdf=INF; 757 } else { 758 ssp_ran_sphere_hg_local(rng, g, sample_local, pdf); 759 d33_muld3(sample, d33_basis(basis, up), sample_local); 760 } 761 return sample; 762 } 763 764 static INLINE float* 765 ssp_ran_sphere_hg_float 766 (struct ssp_rng* rng, 767 const float up[3], 768 const float g, 769 float sample[3], 770 float* pdf) /* Can be NULL => pdf not computed */ 771 { 772 float sample_local[3]; 773 float basis[9]; 774 ASSERT(-1 <= g && g <= +1 && rng && up && sample && f3_is_normalized(up)); 775 if(absf(g) == 1) { 776 f3_mulf(sample, up, g); 777 if(pdf) *pdf = (float)INF; 778 } else { 779 ssp_ran_sphere_hg_float_local(rng, g, sample_local, pdf); 780 f33_mulf3(sample, f33_basis(basis, up), sample_local); 781 } 782 return sample; 783 } 784 785 /* Return the probability distribution for a Henyey-Greenstein random 786 * variate with an explicit incident direction that is `up' */ 787 SSP_API double 788 ssp_ran_sphere_hg_pdf 789 (const double up[3], 790 const double g, 791 const double sample[3]); 792 793 SSP_API float 794 ssp_ran_sphere_hg_float_pdf 795 (const float up[3], 796 const float g, 797 const float sample[3]); 798 799 /******************************************************************************* 800 * General discrete distribution with state 801 ******************************************************************************/ 802 /* Create a discrete random variate data structure from a list of weights. 803 * `weights' contain the weights of `nweights' discrete events. Its elements 804 * must be positive but they needn't add up to one. */ 805 SSP_API res_T 806 ssp_ranst_discrete_create 807 (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ 808 struct ssp_ranst_discrete** ran); 809 810 SSP_API res_T 811 ssp_ranst_discrete_setup 812 (struct ssp_ranst_discrete* discrete, 813 const double* weights, 814 const size_t nweights); 815 816 SSP_API res_T 817 ssp_ranst_discrete_ref_get 818 (struct ssp_ranst_discrete* ran); 819 820 SSP_API res_T 821 ssp_ranst_discrete_ref_put 822 (struct ssp_ranst_discrete* ran); 823 824 /* Return the index of the sampled discret event. */ 825 SSP_API size_t 826 ssp_ranst_discrete_get 827 (struct ssp_rng* rng, 828 const struct ssp_ranst_discrete* ran); 829 830 /* Return the PDF of the discret event `i'. */ 831 SSP_API double 832 ssp_ranst_discrete_pdf 833 (const size_t i, 834 const struct ssp_ranst_discrete* ran); 835 836 /******************************************************************************* 837 * Gaussian distribution with state 838 ******************************************************************************/ 839 SSP_API res_T 840 ssp_ranst_gaussian_create 841 (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ 842 struct ssp_ranst_gaussian** ran); 843 844 SSP_API res_T 845 ssp_ranst_gaussian_float_create 846 (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ 847 struct ssp_ranst_gaussian_float** ran); 848 849 SSP_API res_T 850 ssp_ranst_gaussian_setup 851 (struct ssp_ranst_gaussian* ran, 852 const double mu, 853 const double sigma); 854 855 SSP_API res_T 856 ssp_ranst_gaussian_float_setup 857 (struct ssp_ranst_gaussian_float* ran, 858 const float mu, 859 const float sigma); 860 861 SSP_API res_T 862 ssp_ranst_gaussian_float_ref_get 863 (struct ssp_ranst_gaussian_float* ran); 864 865 SSP_API res_T 866 ssp_ranst_gaussian_ref_get 867 (struct ssp_ranst_gaussian* ran); 868 869 SSP_API res_T 870 ssp_ranst_gaussian_ref_put 871 (struct ssp_ranst_gaussian* ran); 872 873 SSP_API res_T 874 ssp_ranst_gaussian_float_ref_put 875 (struct ssp_ranst_gaussian_float* ran); 876 877 SSP_API double 878 ssp_ranst_gaussian_get 879 (const struct ssp_ranst_gaussian* ran, 880 struct ssp_rng* rng); 881 882 SSP_API float 883 ssp_ranst_gaussian_float_get 884 (const struct ssp_ranst_gaussian_float* ran, 885 struct ssp_rng* rng); 886 887 SSP_API double 888 ssp_ranst_gaussian_pdf 889 (const struct ssp_ranst_gaussian* ran, 890 const double x); 891 892 SSP_API float 893 ssp_ranst_gaussian_float_pdf 894 (const struct ssp_ranst_gaussian_float* ran, 895 const float x); 896 897 /******************************************************************************* 898 * Piecewise linear distribution with state 899 ******************************************************************************/ 900 SSP_API res_T 901 ssp_ranst_piecewise_linear_create 902 (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ 903 struct ssp_ranst_piecewise_linear **ran); 904 905 SSP_API res_T 906 ssp_ranst_piecewise_linear_float_create 907 (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ 908 struct ssp_ranst_piecewise_linear_float **ran); 909 910 SSP_API res_T 911 ssp_ranst_piecewise_linear_setup 912 (struct ssp_ranst_piecewise_linear *ran, 913 const double* intervals, 914 const double* weights, 915 const size_t size); 916 917 SSP_API res_T 918 ssp_ranst_piecewise_linear_float_setup 919 (struct ssp_ranst_piecewise_linear_float *ran, 920 const float* intervals, 921 const float* weights, 922 const size_t size); 923 924 SSP_API res_T 925 ssp_ranst_piecewise_linear_ref_get 926 (struct ssp_ranst_piecewise_linear* ran); 927 928 SSP_API res_T 929 ssp_ranst_piecewise_linear_float_ref_get 930 (struct ssp_ranst_piecewise_linear_float* ran); 931 932 SSP_API res_T 933 ssp_ranst_piecewise_linear_ref_put 934 (struct ssp_ranst_piecewise_linear* ran); 935 936 SSP_API res_T 937 ssp_ranst_piecewise_linear_float_ref_put 938 (struct ssp_ranst_piecewise_linear_float* ran); 939 940 SSP_API double 941 ssp_ranst_piecewise_linear_get 942 (const struct ssp_ranst_piecewise_linear *ran, 943 struct ssp_rng* rng); 944 945 SSP_API float 946 ssp_ranst_piecewise_linear_float_get 947 (const struct ssp_ranst_piecewise_linear_float *ran, 948 struct ssp_rng* rng); 949 950 SSP_API double 951 ssp_ranst_piecewise_linear_pdf 952 (const struct ssp_ranst_piecewise_linear *ran, 953 double x); 954 955 SSP_API float 956 ssp_ranst_piecewise_linear_float_pdf 957 (const struct ssp_ranst_piecewise_linear_float *ran, 958 float x); 959 960 /******************************************************************************* 961 * Uniform distribution of a point into a disk. 962 ******************************************************************************/ 963 SSP_API double* 964 ssp_ran_uniform_disk_local 965 (struct ssp_rng* rng, 966 const double radius, 967 double pt[3], /* pt[2] <= 0 */ 968 double* pdf); /* Can be NULL => pdf not computed */ 969 970 SSP_API float* 971 ssp_ran_uniform_disk_float_local 972 (struct ssp_rng* rng, 973 const float radius, 974 float pt[3], /* pt[2] <= 0 */ 975 float* pdf); /* Can be NULL => pdf not computed */ 976 977 static INLINE double 978 ssp_ran_uniform_disk_local_pdf(const double radius) 979 { 980 return 1 / (radius * radius); 981 } 982 983 static INLINE float 984 ssp_ran_uniform_disk_float_local_pdf(const float radius) 985 { 986 return 1 / (radius * radius); 987 } 988 989 static INLINE double* 990 ssp_ran_uniform_disk 991 (struct ssp_rng* rng, 992 const double radius, 993 const double up[3], 994 double pt[3], 995 double* pdf) /* Can be NULL => pdf not computed */ 996 { 997 double sample_local[3]; 998 double basis[9]; 999 ASSERT(rng && up && pt && radius > 0); 1000 ssp_ran_uniform_disk_local(rng, radius, sample_local, pdf); 1001 return d33_muld3(pt, d33_basis(basis, up), sample_local); 1002 } 1003 1004 static INLINE float* 1005 ssp_ran_uniform_disk_float 1006 (struct ssp_rng* rng, 1007 const float radius, 1008 const float up[3], 1009 float pt[3], 1010 float* pdf) /* Can be NULL => pdf not computed */ 1011 { 1012 float sample_local[3]; 1013 float basis[9]; 1014 ASSERT(rng && up && pt && radius > 0); 1015 ssp_ran_uniform_disk_float_local(rng, radius, sample_local, pdf); 1016 return f33_mulf3(pt, f33_basis(basis, up), sample_local); 1017 } 1018 1019 /******************************************************************************* 1020 * Uniform distribution of a point onto a sphere cap 1021 * pdf = 1/(2*PI*(1-cos(aperture)) 1022 ******************************************************************************/ 1023 /* Uniform sampling of unit sphere cap centered in zero whose up direction 1024 * is implicity the Z axis. */ 1025 SSP_API double* 1026 ssp_ran_sphere_cap_uniform_local 1027 (struct ssp_rng* rng, 1028 /* In [-1, 1]. Height where the sphere cap begins. It equal cos(theta_max) 1029 * where theta_max is the aperture angle from the up axis */ 1030 const double height, 1031 double sample[3], 1032 double* pdf); /* Can be NULL => pdf not computed */ 1033 1034 SSP_API float* 1035 ssp_ran_sphere_cap_uniform_float_local 1036 (struct ssp_rng* rng, 1037 const float height, /* In [-1, 1]. Height of the sphere cap. */ 1038 float sample[3], 1039 float* pdf); /* Can be NULL => pdf not computed */ 1040 1041 static INLINE double 1042 ssp_ran_sphere_cap_uniform_pdf(const double height) 1043 { 1044 ASSERT(height <= 1.0 && height >= -1.0); 1045 return height == 1.0 ? INF : 1.0/(2.0*PI*(1-height)); 1046 } 1047 1048 static INLINE float 1049 ssp_ran_sphere_cap_uniform_float_pdf(const float height) 1050 { 1051 ASSERT(height <= 1.f && height >= -1.f); 1052 return height == 1.f ? (float)INF : 1.f/(2.f*(float)PI*(1.f-height)); 1053 } 1054 1055 /* Uniform sampling of unit sphere cap centered in zero whose up direction 1056 * is explicitly defined */ 1057 static INLINE double* 1058 ssp_ran_sphere_cap_uniform 1059 (struct ssp_rng* rng, 1060 const double up[3], 1061 const double height, 1062 double sample[3], 1063 double* pdf) /* Can be NULL */ 1064 { 1065 double sample_local[3]; 1066 double basis[9]; 1067 ASSERT(up); 1068 ssp_ran_sphere_cap_uniform_local(rng, height, sample_local, pdf); 1069 return d33_muld3(sample, d33_basis(basis, up), sample_local); 1070 } 1071 1072 static INLINE float* 1073 ssp_ran_sphere_cap_uniform_float 1074 (struct ssp_rng* rng, 1075 const float up[3], 1076 const float height, 1077 float sample[3], 1078 float* pdf) 1079 { 1080 float sample_local[3]; 1081 float basis[9]; 1082 ASSERT(up); 1083 ssp_ran_sphere_cap_uniform_float_local(rng, height, sample_local, pdf); 1084 return f33_mulf3(sample, f33_basis(basis, up), sample_local); 1085 } 1086 1087 /******************************************************************************* 1088 * Uniform distribution of a point onto a truncated sphere cap 1089 * pdf = 1/(2*PI*(cos(aperture_max)-cos(aperture_min)) 1090 ******************************************************************************/ 1091 /* Uniform sampling of unit sphere cap centered in zero whose up direction 1092 * is implicity the Z axis. */ 1093 SSP_API double* 1094 ssp_ran_spherical_zone_uniform_local 1095 (struct ssp_rng* rng, 1096 const double height_range[2], 1097 double pt[3], 1098 double* pdf); 1099 1100 SSP_API float* 1101 ssp_ran_spherical_zone_uniform_float_local 1102 (struct ssp_rng* rng, 1103 const float height_range[2], 1104 float pt[3], 1105 float* pdf); 1106 1107 /* Uniform sampling of unit truncated sphere cap centered in zero whose 1108 * up direction is explicitly defined */ 1109 static INLINE double* 1110 ssp_ran_spherical_zone_uniform 1111 (struct ssp_rng* rng, 1112 const double up[3], 1113 const double height_range[2], 1114 double sample[3], 1115 double* pdf) /* Can be NULL */ 1116 { 1117 double sample_local[3]; 1118 double basis[9]; 1119 ASSERT(up); 1120 ssp_ran_spherical_zone_uniform_local(rng, height_range, sample_local, pdf); 1121 return d33_muld3(sample, d33_basis(basis, up), sample_local); 1122 } 1123 1124 static INLINE float* 1125 ssp_ran_spherical_zone_uniform_float 1126 (struct ssp_rng* rng, 1127 const float up[3], 1128 const float height_range[2], 1129 float sample[3], 1130 float* pdf) /* Can be NULL */ 1131 { 1132 float sample_local[3]; 1133 float basis[9]; 1134 ASSERT(up); 1135 ssp_ran_spherical_zone_uniform_float_local(rng, height_range, sample_local, pdf); 1136 return f33_mulf3(sample, f33_basis(basis, up), sample_local); 1137 } 1138 1139 static INLINE double 1140 ssp_ran_spherical_zone_uniform_pdf(const double height_range[2]) 1141 { 1142 ASSERT(height_range && height_range[0] <= height_range[1] 1143 && -1 <= height_range[0] && height_range[1] <= 1); 1144 if(height_range[0] == height_range[1]) { 1145 return height_range[0] == 1 ? INF : ssp_ran_circle_uniform_pdf(); 1146 } 1147 return 1.0/(2.0*PI*(height_range[1]-height_range[0])); 1148 } 1149 1150 static INLINE float 1151 ssp_ran_spherical_zone_uniform_float_pdf(const float height_range[2]) 1152 { 1153 ASSERT(height_range && height_range[0] <= height_range[1] 1154 && -1 <= height_range[0] && height_range[1] <= 1); 1155 if(height_range[0] == height_range[1]) { 1156 return height_range[0] == 1 ? 1157 (float)INF : ssp_ran_circle_uniform_float_pdf(); 1158 } 1159 return 1.f/(2.f*(float)PI*(height_range[1]-height_range[0])); 1160 } 1161 1162 END_DECLS 1163 1164 #endif /* SSP_H */