ssp_ran.c (19150B)
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 #include "ssp_rng_c.h" 17 18 static FINLINE float 19 sinf2cosf(const float d) 20 { 21 return sqrtf(MMAX(0, 1 - d*d)); 22 } 23 24 static FINLINE float 25 cosf2sinf(const float d) 26 { 27 return sinf2cosf(d); 28 } 29 30 /******************************************************************************* 31 * Exported state free random variates 32 ******************************************************************************/ 33 double 34 ssp_ran_exp(struct ssp_rng* rng, const double mu) 35 { 36 ASSERT(rng && mu >= 0); 37 RAN_NAMESPACE::exponential_distribution<double> distribution(mu); 38 return wrap_ran(*rng, distribution); 39 } 40 41 float 42 ssp_ran_exp_float(struct ssp_rng* rng, const float mu) 43 { 44 ASSERT(rng && mu >= 0); 45 RAN_NAMESPACE::exponential_distribution<float> distribution(mu); 46 return wrap_ran(*rng, distribution); 47 } 48 49 double 50 ssp_ran_exp_pdf(const double x, const double mu) 51 { 52 ASSERT(x >= 0 && mu > 0); 53 return mu * exp(-x * mu); 54 } 55 56 float 57 ssp_ran_exp_float_pdf(const float x, const float mu) 58 { 59 ASSERT(x >= 0 && mu > 0); 60 return mu * expf(-x * mu); 61 } 62 63 double 64 ssp_ran_exp_truncated(struct ssp_rng* rng, const double mu, const double max) 65 { 66 double u, r; 67 ASSERT(rng && mu > 0 && max > 0); 68 u = ssp_rng_canonical(rng); 69 r = -log(1 - u * (1 - exp(-mu * max))) / mu; 70 return r; 71 } 72 73 float 74 ssp_ran_exp_truncated_float(struct ssp_rng* rng, const float mu, const float max) 75 { 76 float u, r; 77 ASSERT(rng && mu > 0 && max > 0); 78 u = ssp_rng_canonical_float(rng); 79 r = -logf(1 - u * (1 - expf(-mu * max))) / mu; 80 return r; 81 } 82 83 double 84 ssp_ran_exp_truncated_pdf(const double x, const double mu, const double max) 85 { 86 ASSERT(x >= 0 && mu > 0 && max > 0); 87 if(x > max) 88 return 0; 89 else { 90 double norm = mu / (1 - exp(-max * mu)); 91 return norm * exp(-x * mu); 92 } 93 } 94 95 float 96 ssp_ran_exp_truncated_float_pdf(const float x, const float mu, const float max) 97 { 98 ASSERT(x >= 0 && mu > 0 && max > 0); 99 if(x > max) 100 return 0; 101 else { 102 float norm = mu / (1 - expf(-max * mu)); 103 return norm * expf(-x * mu); 104 } 105 } 106 107 double 108 ssp_ran_gaussian 109 (struct ssp_rng* rng, 110 const double mu, 111 const double sigma) 112 { 113 ASSERT(rng); 114 RAN_NAMESPACE::normal_distribution<double> distribution(mu, sigma); 115 return wrap_ran(*rng, distribution); 116 } 117 118 float 119 ssp_ran_gaussian_float 120 (struct ssp_rng* rng, 121 const float mu, 122 const float sigma) 123 { 124 ASSERT(rng); 125 RAN_NAMESPACE::normal_distribution<float> distribution(mu, sigma); 126 return wrap_ran(*rng, distribution); 127 } 128 129 double 130 ssp_ran_gaussian_pdf 131 (const double x, 132 const double mu, 133 const double sigma) 134 { 135 const double tmp = (x - mu) / sigma; 136 return 1 / (sigma * SQRT_2_PI) * exp(-0.5 * tmp * tmp); 137 } 138 139 float 140 ssp_ran_gaussian_float_pdf 141 (const float x, 142 const float mu, 143 const float sigma) 144 { 145 const float tmp = (x - mu) / sigma; 146 return 1 / (sigma * (float)SQRT_2_PI) * expf(-0.5f * tmp * tmp); 147 } 148 149 double 150 ssp_ran_lognormal 151 (struct ssp_rng* rng, 152 const double zeta, 153 const double sigma) 154 { 155 ASSERT(rng); 156 RAN_NAMESPACE::lognormal_distribution<double> distribution(zeta, sigma); 157 return wrap_ran(*rng, distribution); 158 } 159 160 float 161 ssp_ran_lognormal_float 162 (struct ssp_rng* rng, 163 const float zeta, 164 const float sigma) 165 { 166 ASSERT(rng); 167 RAN_NAMESPACE::lognormal_distribution<float> distribution(zeta, sigma); 168 return wrap_ran(*rng, distribution);; 169 } 170 171 double 172 ssp_ran_lognormal_pdf 173 (const double x, 174 const double zeta, 175 const double sigma) 176 { 177 const double tmp = log(x) - zeta; 178 return 1 / (sigma * x * SQRT_2_PI) * exp(-(tmp*tmp) / (2*sigma*sigma)); 179 } 180 181 float 182 ssp_ran_lognormal_float_pdf 183 (const float x, 184 const float zeta, 185 const float sigma) 186 { 187 const float tmp = logf(x) - zeta; 188 return 1 / (sigma * x * (float)SQRT_2_PI) * expf(-(tmp*tmp) / (2 * sigma*sigma)); 189 } 190 191 double* 192 ssp_ran_sphere_uniform 193 (struct ssp_rng* rng, 194 double sample[3], 195 double* pdf) 196 { 197 double phi, cos_theta, sin_theta, v; 198 ASSERT(rng && sample); 199 phi = ssp_rng_uniform_double(rng, 0, 2 * PI); 200 v = ssp_rng_canonical(rng); 201 cos_theta = 1 - 2 * v; 202 sin_theta = 2 * sqrt(v * (1 - v)); 203 sample[0] = cos(phi) * sin_theta; 204 sample[1] = sin(phi) * sin_theta; 205 sample[2] = cos_theta; 206 if(pdf) *pdf = 1 / (4*PI); 207 return sample; 208 } 209 210 float* 211 ssp_ran_sphere_uniform_float 212 (struct ssp_rng* rng, 213 float sample[3], 214 float* pdf) 215 { 216 float phi, cos_theta, sin_theta, v; 217 ASSERT(rng && sample); 218 phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI); 219 v = ssp_rng_canonical_float(rng); 220 cos_theta = 1 - 2 * v; 221 sin_theta = 2 * sqrtf(v * (1 - v)); 222 sample[0] = cosf(phi) * sin_theta; 223 sample[1] = sinf(phi) * sin_theta; 224 sample[2] = cos_theta; 225 if(pdf) *pdf = 1 / (4*(float)PI); 226 return sample; 227 } 228 229 double* 230 ssp_ran_circle_uniform 231 (struct ssp_rng* rng, 232 double sample[2], 233 double* pdf) 234 { 235 double theta; 236 ASSERT(rng && sample); 237 theta = ssp_rng_uniform_double(rng, 0, 2*PI); 238 sample[0] = cos(theta); 239 sample[1] = sin(theta); 240 if(pdf) *pdf = 1/(2*PI); 241 return sample; 242 } 243 244 float* 245 ssp_ran_circle_uniform_float 246 (struct ssp_rng* rng, 247 float sample[2], 248 float* pdf) 249 { 250 float theta; 251 ASSERT(rng && sample); 252 theta = ssp_rng_uniform_float(rng, 0, 2*(float)PI); 253 sample[0] = cosf(theta); 254 sample[1] = sinf(theta); 255 if(pdf) *pdf = 1/(2*(float)PI); 256 return sample; 257 } 258 259 double* 260 ssp_ran_triangle_uniform 261 (struct ssp_rng* rng, 262 const double v0[3], 263 const double v1[3], 264 const double v2[3], 265 double sample[3], 266 double* pdf) 267 { 268 double sqrt_u, v, one_minus_u; 269 double vec0[3]; 270 double vec1[3]; 271 double tmp[3]; 272 ASSERT(rng && v0 && v1 && v2 && sample); 273 274 sqrt_u = sqrt(ssp_rng_canonical(rng)); 275 v = ssp_rng_canonical(rng); 276 one_minus_u = 1.0 - sqrt_u; 277 278 d3_sub(vec0, v0, v2); 279 d3_sub(vec1, v1, v2); 280 sample[0] = v2[0] + one_minus_u * vec0[0] + v * sqrt_u * vec1[0]; 281 sample[1] = v2[1] + one_minus_u * vec0[1] + v * sqrt_u * vec1[1]; 282 sample[2] = v2[2] + one_minus_u * vec0[2] + v * sqrt_u * vec1[2]; 283 if(pdf) *pdf = 2 / d3_len(d3_cross(tmp, vec0, vec1)); 284 return sample; 285 } 286 287 float* 288 ssp_ran_triangle_uniform_float 289 (struct ssp_rng* rng, 290 const float v0[3], 291 const float v1[3], 292 const float v2[3], 293 float sample[3], 294 float* pdf) 295 { 296 float sqrt_u, v, one_minus_u; 297 float vec0[3]; 298 float vec1[3]; 299 float tmp[3]; 300 ASSERT(rng && v0 && v1 && v2 && sample); 301 302 sqrt_u = sqrtf(ssp_rng_canonical_float(rng)); 303 v = ssp_rng_canonical_float(rng); 304 one_minus_u = 1 - sqrt_u; 305 306 f3_sub(vec0, v0, v2); 307 f3_sub(vec1, v1, v2); 308 sample[0] = v2[0] + one_minus_u * vec0[0] + v * sqrt_u * vec1[0]; 309 sample[1] = v2[1] + one_minus_u * vec0[1] + v * sqrt_u * vec1[1]; 310 sample[2] = v2[2] + one_minus_u * vec0[2] + v * sqrt_u * vec1[2]; 311 if(pdf) *pdf = 2 / f3_len(f3_cross(tmp, vec0, vec1)); 312 return sample; 313 } 314 315 double* 316 ssp_ran_tetrahedron_uniform 317 (struct ssp_rng* rng, 318 const double v0[3], 319 const double v1[3], 320 const double v2[3], 321 const double v3[3], 322 double sample[3], 323 double* pdf) 324 { 325 double a, s, t, u; /* Barycentric coordinates of the sampled point. */ 326 double tmp0[3], tmp1[3], tmp2[3], tmp3[3], tmp4[3], tmp5[3]; 327 ASSERT(rng && v0 && v1 && v2 && v3 && sample); 328 329 s = ssp_rng_canonical(rng); 330 t = ssp_rng_canonical(rng); 331 u = ssp_rng_canonical(rng); 332 333 if(s + t > 1) { /* Cut and fold the cube into a prism */ 334 s = 1 - s; 335 t = 1 - t; 336 } 337 if(t + u > 1) { /* Cut and fold the prism into a tetrahedron */ 338 double swp = u; 339 u = 1 - s - t; 340 t = 1 - swp; 341 } 342 else if(s + t + u > 1) { 343 double swp = u; 344 u = s + t + u - 1; 345 s = 1 - t - swp; 346 } 347 a = 1 - s - t - u; 348 349 d3_add(sample, 350 d3_add(tmp4, d3_muld(tmp0, v0, a), d3_muld(tmp1, v1, s)), 351 d3_add(tmp5, d3_muld(tmp2, v2, t), d3_muld(tmp3, v3, u))); 352 353 if(pdf) 354 *pdf = ssp_ran_tetrahedron_uniform_pdf(v0, v1, v2, v3); 355 356 return sample; 357 } 358 359 float* 360 ssp_ran_tetrahedron_uniform_float 361 (struct ssp_rng* rng, 362 const float v0[3], 363 const float v1[3], 364 const float v2[3], 365 const float v3[3], 366 float sample[3], 367 float* pdf) 368 { 369 float a, s, t, u; /* Barycentric coordinates of the sampled point. */ 370 float tmp0[3], tmp1[3], tmp2[3], tmp3[3], tmp4[3], tmp5[3]; 371 ASSERT(rng && v0 && v1 && v2 && v3 && sample); 372 373 s = ssp_rng_canonical_float(rng); 374 t = ssp_rng_canonical_float(rng); 375 u = ssp_rng_canonical_float(rng); 376 377 if(s + t > 1) { /* Cut and fold the cube into a prism */ 378 s = 1 - s; 379 t = 1 - t; 380 } 381 if(t + u > 1) { /* Cut and fold the prism into a tetrahedron */ 382 float swp = u; 383 u = 1 - s - t; 384 t = 1 - swp; 385 } 386 else if(s + t + u > 1) { 387 float swp = u; 388 u = s + t + u - 1; 389 s = 1 - t - swp; 390 } 391 a = 1 - s - t - u; 392 393 f3_add(sample, 394 f3_add(tmp4, f3_mulf(tmp0, v0, a), f3_mulf(tmp1, v1, s)), 395 f3_add(tmp5, f3_mulf(tmp2, v2, t), f3_mulf(tmp3, v3, u))); 396 397 if(pdf) 398 *pdf = ssp_ran_tetrahedron_uniform_float_pdf(v0, v1, v2, v3); 399 400 return sample; 401 } 402 403 double* 404 ssp_ran_hemisphere_uniform_local 405 (struct ssp_rng* rng, 406 double sample[3], 407 double* pdf) 408 { 409 double phi, cos_theta, sin_theta; 410 ASSERT(rng && sample); 411 phi = ssp_rng_uniform_double(rng, 0, 2 * PI); 412 cos_theta = ssp_rng_canonical(rng); 413 sin_theta = cos2sin(cos_theta); 414 sample[0] = cos(phi) * sin_theta; 415 sample[1] = sin(phi) * sin_theta; 416 sample[2] = cos_theta; 417 if(pdf) *pdf = 1 / (2*PI); 418 return sample; 419 } 420 421 float* 422 ssp_ran_hemisphere_uniform_float_local 423 (struct ssp_rng* rng, 424 float sample[3], 425 float* pdf) 426 { 427 float phi, cos_theta, sin_theta; 428 ASSERT(rng && sample); 429 phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI); 430 cos_theta = ssp_rng_canonical_float(rng); 431 sin_theta = cosf2sinf(cos_theta); 432 sample[0] = cosf(phi) * sin_theta; 433 sample[1] = sinf(phi) * sin_theta; 434 sample[2] = cos_theta; 435 if(pdf) *pdf = 1 / (2*(float)PI); 436 return sample; 437 } 438 439 double* 440 ssp_ran_hemisphere_cos_local 441 (struct ssp_rng* rng, 442 double sample[3], 443 double* pdf) 444 { 445 double phi, cos_theta, sin_theta, v; 446 ASSERT(rng && sample); 447 phi = ssp_rng_uniform_double(rng, 0, 2 * PI); 448 v = ssp_rng_canonical(rng); 449 cos_theta = sqrt(v); 450 sin_theta = sqrt(1 - v); 451 sample[0] = cos(phi) * sin_theta; 452 sample[1] = sin(phi) * sin_theta; 453 sample[2] = cos_theta; 454 if(pdf) *pdf = cos_theta / PI; 455 return sample; 456 } 457 458 float* 459 ssp_ran_hemisphere_cos_float_local 460 (struct ssp_rng* rng, 461 float sample[3], 462 float* pdf) 463 { 464 float phi, cos_theta, sin_theta, v; 465 ASSERT(rng && sample); 466 phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI); 467 v = ssp_rng_canonical_float(rng); 468 cos_theta = sqrtf(v); 469 sin_theta = sqrtf(1 - v); 470 sample[0] = cosf(phi) * sin_theta; 471 sample[1] = sinf(phi) * sin_theta; 472 sample[2] = cos_theta; 473 if(pdf) *pdf = cos_theta / (float)PI; 474 return sample; 475 } 476 477 double 478 ssp_ran_sphere_hg_local_pdf 479 (const double g, 480 const double sample[3]) 481 { 482 double epsilon_g, epsilon_mu; 483 ASSERT(fabs(g) <= 1 && sample && d3_is_normalized(sample)); 484 if(g>0) { 485 epsilon_g = 1 - g; 486 epsilon_mu = 1 - sample[2]; 487 } else { 488 epsilon_g = 1 + g; 489 epsilon_mu = 1 + sample[2]; 490 } 491 return 1 / (4 * PI) 492 *epsilon_g*(2 - epsilon_g) 493 *pow(epsilon_g*epsilon_g + 2 * epsilon_mu*(1 - epsilon_g), -1.5); 494 } 495 496 float 497 ssp_ran_sphere_hg_float_local_pdf 498 (const float g, 499 const float sample[3]) 500 { 501 float epsilon_g, epsilon_mu; 502 ASSERT(fabsf(g) <= 1 && sample && f3_is_normalized(sample)); 503 if(g>0) { 504 epsilon_g = 1 - g; 505 epsilon_mu = 1 - sample[2]; 506 } else { 507 epsilon_g = 1 + g; 508 epsilon_mu = 1 + sample[2]; 509 } 510 return 1 / (4 * (float)PI) 511 *epsilon_g*(2 - epsilon_g) 512 *powf(epsilon_g*epsilon_g + 2 * epsilon_mu*(1 - epsilon_g), -1.5f); 513 } 514 515 double* 516 ssp_ran_sphere_hg_local 517 (struct ssp_rng* rng, 518 const double g, 519 double sample[3], 520 double* pdf) 521 { 522 double phi, R, cos_theta, sin_theta; 523 ASSERT(fabs(g) <= 1 && rng && sample); 524 if(fabs(g) == 1) { 525 d3(sample, 0, 0, g); 526 if(pdf) *pdf = INF; 527 } else { 528 phi = ssp_rng_uniform_double(rng, 0, 2 * PI); 529 R = ssp_rng_canonical(rng); 530 cos_theta = 2 * R*(1 + g)*(1 + g)*(g*(R - 1) + 1) 531 / ((1 - g * (1 - 2 * R))*(1 - g * (1 - 2 * R))) - 1; 532 sin_theta = cos2sin(cos_theta); 533 sample[0] = cos(phi) * sin_theta; 534 sample[1] = sin(phi) * sin_theta; 535 sample[2] = cos_theta; 536 if(pdf) *pdf = ssp_ran_sphere_hg_local_pdf(g, sample); 537 } 538 return sample; 539 } 540 541 float* 542 ssp_ran_sphere_hg_float_local 543 (struct ssp_rng* rng, 544 const float g, 545 float sample[3], 546 float* pdf) 547 { 548 float phi, R, cos_theta, sin_theta; 549 ASSERT(fabsf(g) <= 1 && rng && sample); 550 if(fabsf(g) == 1) { 551 f3(sample, 0, 0, g); 552 if(pdf) *pdf = (float)INF; 553 } else { 554 phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI); 555 R = ssp_rng_canonical_float(rng); 556 cos_theta = 2 * R*(1 + g)*(1 + g)*(g*(R - 1) + 1) 557 / ((1 - g * (1 - 2 * R))*(1 - g * (1 - 2 * R))) - 1; 558 sin_theta = cosf2sinf(cos_theta); 559 sample[0] = cosf(phi) * sin_theta; 560 sample[1] = sinf(phi) * sin_theta; 561 sample[2] = cos_theta; 562 if(pdf) *pdf = ssp_ran_sphere_hg_float_local_pdf(g, sample); 563 } 564 return sample; 565 } 566 567 double 568 ssp_ran_sphere_hg_pdf 569 (const double up[3], 570 const double g, 571 const double sample[3]) 572 { 573 double epsilon_g, epsilon_mu; 574 ASSERT(-1 <= g && g <= +1 && 575 up && sample && d3_is_normalized(up) && d3_is_normalized(sample)); 576 if(fabs(g) == 1) return INF; 577 if(g>0) { 578 epsilon_g = 1 - g; 579 epsilon_mu = 1 - d3_dot(sample, up); 580 } else { 581 epsilon_g = 1 + g; 582 epsilon_mu = 1 + d3_dot(sample, up); 583 } 584 return 1 / (4 * PI) 585 *epsilon_g*(2 - epsilon_g) 586 *pow(epsilon_g*epsilon_g + 2 * epsilon_mu*(1 - epsilon_g), -1.5); 587 } 588 589 float 590 ssp_ran_sphere_hg_float_pdf 591 (const float up[3], 592 const float g, 593 const float sample[3]) 594 { 595 float epsilon_g, epsilon_mu; 596 ASSERT(-1 <= g && g <= +1 && 597 up && sample && f3_is_normalized(up) && f3_is_normalized(sample)); 598 if(fabsf(g) == 1) return (float)INF; 599 if(g>0) { 600 epsilon_g = 1 - g; 601 epsilon_mu = 1 - f3_dot(sample, up); 602 } else { 603 epsilon_g = 1 + g; 604 epsilon_mu = 1 + f3_dot(sample, up); 605 } 606 return 1 / (4 * (float)PI) 607 *epsilon_g*(2 - epsilon_g) 608 *powf(epsilon_g*epsilon_g + 2 * epsilon_mu*(1 - epsilon_g), -1.5f); 609 } 610 611 double* 612 ssp_ran_uniform_disk_local 613 (struct ssp_rng* rng, 614 const double radius, 615 double pt[3], 616 double* pdf) 617 { 618 double theta, r; 619 ASSERT(rng && pt && radius > 0); 620 theta = ssp_rng_uniform_double(rng, 0, 2 * PI); 621 r = radius * sqrt(ssp_rng_canonical(rng)); 622 pt[0] = r * cos(theta); 623 pt[1] = r * sin(theta); 624 pt[2] = 0; 625 if(pdf) *pdf = 1 / (radius * radius); 626 return pt; 627 } 628 629 float* 630 ssp_ran_uniform_disk_float_local 631 (struct ssp_rng* rng, 632 const float radius, 633 float pt[3], 634 float* pdf) 635 { 636 float theta, r; 637 ASSERT(rng && pt && radius > 0); 638 theta = ssp_rng_uniform_float(rng, 0, 2 * (float)PI); 639 r = radius * sqrtf(ssp_rng_canonical_float(rng)); 640 pt[0] = r * cosf(theta); 641 pt[1] = r * sinf(theta); 642 pt[2] = 0; 643 if(pdf) *pdf = 1 / (radius * radius); 644 return pt; 645 } 646 647 double* 648 ssp_ran_sphere_cap_uniform_local 649 (struct ssp_rng* rng, 650 const double height, 651 double pt[3], 652 double* pdf) 653 { 654 ASSERT(rng && pt && (height >= -1 || height <= 1)); 655 656 if(height == 1) { 657 d3(pt, 0, 0, 1); 658 if(pdf) *pdf = INF; 659 } else if(height == -1) { 660 ssp_ran_sphere_uniform(rng, pt, pdf); 661 } else { 662 const double cos_aperture = height; 663 double phi, cos_theta, sin_theta; 664 phi = ssp_rng_uniform_double(rng, 0, 2.0*PI); 665 cos_theta = ssp_rng_uniform_double(rng, cos_aperture, 1); 666 sin_theta = cos2sin(cos_theta); 667 pt[0] = cos(phi) * sin_theta; 668 pt[1] = sin(phi) * sin_theta; 669 pt[2] = cos_theta; 670 if(pdf) *pdf = 1.0/(2.0*PI*(1.0-cos_aperture)); 671 } 672 return pt; 673 } 674 675 float* 676 ssp_ran_sphere_cap_uniform_float_local 677 (struct ssp_rng* rng, 678 const float height, 679 float pt[3], 680 float* pdf) 681 { 682 ASSERT(rng && pt && (height >= -1 || height <= 1)); 683 if(height == 1) { 684 f3(pt, 0, 0, 1); 685 if(pdf) *pdf = (float)INF; 686 } else if(height == -1) { 687 ssp_ran_sphere_uniform_float(rng, pt, pdf); 688 } else { 689 const float cos_aperture = height; 690 float phi, cos_theta, sin_theta; 691 phi = ssp_rng_uniform_float(rng, 0, 2.f*(float)PI); 692 cos_theta = ssp_rng_uniform_float(rng, cos_aperture, 1); 693 sin_theta = (float)cos2sin(cos_theta); 694 pt[0] = cosf(phi) * sin_theta; 695 pt[1] = sinf(phi) * sin_theta; 696 pt[2] = cos_theta; 697 if(pdf) *pdf = 1.f/(2.f*(float)PI*(1.f-cos_aperture)); 698 } 699 return pt; 700 } 701 702 double* 703 ssp_ran_spherical_zone_uniform_local 704 (struct ssp_rng* rng, 705 const double height_range[2], 706 double pt[3], 707 double* pdf) 708 { 709 ASSERT(rng && pt && height_range 710 && height_range[0] <= height_range[1] 711 && height_range[0] >= -1 && height_range[1] <= 1); 712 713 if(height_range[1] == 1) { 714 ssp_ran_sphere_cap_uniform_local(rng, height_range[0], pt, pdf); 715 } 716 else if(height_range[0] == height_range[1]) { 717 double sin_theta = cos2sin(height_range[0]); 718 ssp_ran_circle_uniform(rng, pt, pdf); 719 pt[0] *= sin_theta; 720 pt[1] *= sin_theta; 721 pt[2] = height_range[0]; 722 } else { 723 const double cos_aperture_min = height_range[0]; 724 const double cos_aperture_max = height_range[1]; 725 double phi, cos_theta, sin_theta; 726 phi = ssp_rng_uniform_double(rng, 0, 2.0*PI); 727 cos_theta = ssp_rng_uniform_double(rng, cos_aperture_min, cos_aperture_max); 728 sin_theta = cos2sin(cos_theta); 729 pt[0] = cos(phi) * sin_theta; 730 pt[1] = sin(phi) * sin_theta; 731 pt[2] = cos_theta; 732 if(pdf) *pdf = 1.0/(2.0*PI*(cos_aperture_max-cos_aperture_min)); 733 } 734 return pt; 735 } 736 737 float* 738 ssp_ran_spherical_zone_uniform_float_local 739 (struct ssp_rng* rng, 740 const float height_range[2], 741 float pt[3], 742 float* pdf) 743 { 744 ASSERT(rng && pt && height_range 745 && height_range[0] <= height_range[1] 746 && height_range[0] >= -1 && height_range[1] <= 1); 747 748 if(height_range[1] == 1) { 749 ssp_ran_sphere_cap_uniform_float_local(rng, height_range[0], pt, pdf); 750 } 751 else if(height_range[0] == height_range[1]) { 752 float sin_theta = (float)cos2sin(height_range[0]); 753 ssp_ran_circle_uniform_float(rng, pt, pdf); 754 pt[0] *= sin_theta; 755 pt[1] *= sin_theta; 756 pt[2] = height_range[0]; 757 } else { 758 const float cos_aperture_min = height_range[0]; 759 const float cos_aperture_max = height_range[1]; 760 float phi, cos_theta, sin_theta; 761 phi = ssp_rng_uniform_float(rng, 0, 2.f*(float)PI); 762 cos_theta = ssp_rng_uniform_float(rng, cos_aperture_min, cos_aperture_max); 763 sin_theta = (float)cos2sin(cos_theta); 764 pt[0] = cosf(phi) * sin_theta; 765 pt[1] = sinf(phi) * sin_theta; 766 pt[2] = cos_theta; 767 if(pdf) *pdf = 1.f/(2.f*(float)PI*(cos_aperture_max-cos_aperture_min)); 768 } 769 return pt; 770 } 771