|
57 | 57 | basic characteristic of simple programs producing noisy but ever-better answers is what MC is all
|
58 | 58 | about, and it is especially good for applications like graphics where great accuracy is not needed.
|
59 | 59 |
|
60 |
| -As an example, let’s estimate $ \pi $. There are many ways to do this, with the Buffon Needle problem |
61 |
| -being a classic case study. We’ll do a variation inspired by that. Suppose you have a circle |
| 60 | +As an example, let’s estimate $\pi$. There are many ways to do this, with the Buffon Needle |
| 61 | +problem being a classic case study. We’ll do a variation inspired by that. Suppose you have a circle |
62 | 62 | inscribed inside a square:
|
63 | 63 |
|
64 | 64 | 
|
|
67 | 67 | up inside the circle should be proportional to the area of the circle. The exact fraction should in
|
68 | 68 | fact be the ratio of the circle area to the square area. Fraction:
|
69 | 69 |
|
70 |
| - $$ \frac{Pi \cdot R^2}{(2R)^2} = \frac{Pi}{4} $$ |
| 70 | + $$ \frac{Pi \cdot R^2}{(2R)^2} = \frac{Pi}{4} $$ |
71 | 71 |
|
72 | 72 | Since the *R* cancels out, we can pick whatever is computationally convenient. Let’s go with R=1
|
73 | 73 | centered at the origin:
|
|
272 | 272 | How do we generate a random number with that pdf $p(r)$? For that we will need some more machinery.
|
273 | 273 | Don’t worry this doesn’t go on forever!
|
274 | 274 |
|
275 |
| -Given a random number from `d = random_double()` that is uniform and between 0 and 1, we should be able |
276 |
| -to find some function $f(d)$ that gives us what we want. Suppose $e = f(d) = d^2$. That is no |
| 275 | +Given a random number from `d = random_double()` that is uniform and between 0 and 1, we should be |
| 276 | +able to find some function $f(d)$ that gives us what we want. Suppose $e = f(d) = d^2$. That is no |
277 | 277 | longer a uniform _pdf_ . The _pdf_ of $e$ will be bigger near 0 than it is near 1 (squaring a
|
278 | 278 | number between 0 and 1 makes it smaller). To take this general observation to a function, we need
|
279 | 279 | the cumulative probability distribution function $P(x)$:
|
|
294 | 294 |
|
295 | 295 | This says _the probability that a random variable with our pdf is less than one is 25%_ . This
|
296 | 296 | gives rise to a clever observation that underlies many methods to generate non-uniform random
|
297 |
| -numbers. We want a function `f()` that when we call it as `f(random_double())` we get a return value with |
298 |
| -a pdf $\frac{x^2}{4}$ . We don’t know what that is, but we do know that 25% of what it returns |
| 297 | +numbers. We want a function `f()` that when we call it as `f(random_double())` we get a return value |
| 298 | +with a pdf $\frac{x^2}{4}$ . We don’t know what that is, but we do know that 25% of what it returns |
299 | 299 | should be less than 1, and 75% should be above one. If $f()$ is increasing, then we would expect
|
300 | 300 | $f(0.25) = 1.0$. This can be generalized to figure out $f()$ for every possible input:
|
301 | 301 |
|
|
324 | 324 |
|
325 | 325 | $$ e = \sqrt{4*random_double()} $$
|
326 | 326 |
|
327 |
| -Note that does range 0 to 2 as hoped, and if we send in $1/4$ for `random_double()` we get 1 as desired. |
| 327 | +Note that does range 0 to 2 as hoped, and if we send in $1/4$ for `random_double()` we get 1 as |
| 328 | +desired. |
328 | 329 |
|
329 | 330 | We can now sample our old integral
|
330 | 331 |
|
|
448 | 449 | unit-sphere. The same methodology as before applies. But now we need to have a pdf defined over 2D.
|
449 | 450 | Suppose we have this integral over all directions:
|
450 | 451 |
|
451 |
| - $$ \int cos^2(\theta) $$ |
| 452 | + $$ \int cos^2(\theta) $$ |
452 | 453 |
|
453 | 454 | By MC integration, we should just be able to sample $cos^2(\theta) / p(direction)$. But what is
|
454 | 455 | direction in that context? We could make it based on polar coordinates, so $p$ would be in terms of
|
|
498 | 499 | }
|
499 | 500 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
500 | 501 |
|
501 |
| -The analytic answer (if you remember enough advanced calc, check me!) is $\frac{4}{3} \pi$, and the code |
502 |
| -above produces that. Next, we are ready to apply that in ray tracing! |
| 502 | +The analytic answer (if you remember enough advanced calc, check me!) is $\frac{4}{3} \pi$, and the |
| 503 | +code above produces that. Next, we are ready to apply that in ray tracing! |
503 | 504 |
|
504 | 505 | The key point here is that all the integrals and probability and all that are over the unit sphere.
|
505 | 506 | The area on the unit sphere is how you measure the directions. Call it direction, solid angle, or
|
506 |
| -area-- it’s all the same thing. Solid angle is the term usually used. If you are comfortable with |
| 507 | +area -- it’s all the same thing. Solid angle is the term usually used. If you are comfortable with |
507 | 508 | that, great! If not, do what I do and imagine the area on the unit sphere that a set of directions
|
508 |
| -goes through. The solid angle $\omega$ and the projected area $A$ on the unit sphere are the same thing. |
| 509 | +goes through. The solid angle $\omega$ and the projected area $A$ on the unit sphere are the same |
| 510 | +thing. |
509 | 511 |
|
510 | 512 | 
|
511 | 513 |
|
|
524 | 526 | First, is the light absorbed?
|
525 | 527 |
|
526 | 528 | Probability of light scattering: $A$
|
| 529 | + |
527 | 530 | Probability of light being absorbed: $1-A$
|
528 | 531 |
|
529 | 532 | Here $A$ stands for _albedo_ (latin for _whiteness_). Albedo is a precise technical term in some
|
|
535 | 538 |
|
536 | 539 | If the light does scatter, it will have a directional distribution that we can describe as a pdf
|
537 | 540 | over solid angle. I will refer to this as its _scattering pdf_: $s(direction)$. The scattering pdf
|
538 |
| -can also vary with incident direction, as you will notice when you look at reflections off a road-- |
| 541 | +can also vary with incident direction, as you will notice when you look at reflections off a road -- |
539 | 542 | they become mirror-like as your viewing angle approaches grazing.
|
540 | 543 |
|
541 | 544 | The color of a surface in terms of these quantities is:
|
542 | 545 |
|
543 |
| - $$ Color = \int A \cdot s(direction) \cdot color(direction) $$ |
| 546 | + $$ Color = \int A \cdot s(direction) \cdot color(direction) $$ |
544 | 547 |
|
545 | 548 | Note that $A$ and $s()$ might depend on the view direction, so of course color can vary with view
|
546 | 549 | direction. Also $A$ and $s()$ may vary with position on the surface or within the volume.
|
547 | 550 |
|
548 | 551 | If we apply the MC basic formula we get the following statistical estimate:
|
549 | 552 |
|
550 |
| - $$ Color = (A \cdot s(direction) \cdot color(direction)) / p(direction) $$ |
| 553 | + $$ Color = (A \cdot s(direction) \cdot color(direction)) / p(direction) $$ |
551 | 554 |
|
552 | 555 | Where $p(direction)$ is the pdf of whatever direction we randomly generate.
|
553 | 556 |
|
|
559 | 562 |
|
560 | 563 | To see that remember that in spherical coordinates remember that:
|
561 | 564 |
|
562 |
| - $$ \delta A = sin(\theta) \delta \theta \delta \phi $$ |
| 565 | + $$ \delta A = sin(\theta) \delta \theta \delta \phi $$ |
563 | 566 |
|
564 | 567 | So:
|
565 | 568 |
|
566 |
| - $$ Area = \int_{0}^{2 \pi} \int_{0}^{\pi / 2} cos(\theta)sin(\theta) \delta \theta \delta \phi = |
| 569 | + $$ Area = \int_{0}^{2 \pi} \int_{0}^{\pi / 2} cos(\theta)sin(\theta) \delta \theta \delta \phi = |
567 | 570 | 2 \pi \cdot \frac{1}{2} = \pi $$
|
568 | 571 |
|
569 | 572 | So for a Lambertian surface the scattering pdf is:
|
570 | 573 |
|
571 |
| - $$ s(direction) = cos(\theta) / \pi $$ |
| 574 | + $$ s(direction) = cos(\theta) / \pi $$ |
572 | 575 |
|
573 | 576 | If we sample using the same pdf, so $p(direction) = cos(\theta) / \pi$, the numerator and
|
574 | 577 | denominator cancel out and we get:
|
575 | 578 |
|
576 |
| - $$ Color = A * s(direction) $$ |
| 579 | + $$ Color = A * s(direction) $$ |
577 | 580 |
|
578 | 581 | This is exactly what we had in our original color() function! But we need to generalize now so we
|
579 | 582 | can send extra rays in important directions such as toward the lights.
|
|
584 | 587 | If you read the literature, you’ll see reflection described by the bidirectional reflectance
|
585 | 588 | distribution function (BRDF). It relates pretty simply to our terms:
|
586 | 589 |
|
587 |
| - $$ BRDF = A * s(direction) / cos(\theta) $$ |
| 590 | + $$ BRDF = A * s(direction) / cos(\theta) $$ |
588 | 591 |
|
589 | 592 | So for a Lambertian surface for example, $BRDF = A / \pi$. Translation between our terms and BRDF is
|
590 | 593 | easy.
|
|
604 | 607 | linear mixtures of them to form mixture densities that are also pdfs. For example, the simplest
|
605 | 608 | would be:
|
606 | 609 |
|
607 |
| - $$ p(direction) = 0.5 \cdot pLight(direction) + 0.5* \cdot pSurface(direction) $$ |
| 610 | + $$ p(direction) = 0.5 \cdot pLight(direction) + 0.5* \cdot pSurface(direction) $$ |
608 | 611 |
|
609 | 612 | As long as the weights are positive and add up to one, any such mixture of pdfs is a pdf. And
|
610 |
| -remember, we can use any pdf-- it doesn’t change the answer we converge to! So, the game is to |
611 |
| -figure out how to make the pdf larger where the product $s(direction)*color(direction)$ is large. For |
612 |
| -diffuse surfaces, this is mainly a matter of guessing where $color(direction)$ is high. |
| 613 | +remember, we can use any pdf -- it doesn’t change the answer we converge to! So, the game is to |
| 614 | +figure out how to make the pdf larger where the product $s(direction)*color(direction)$ is large. |
| 615 | +For diffuse surfaces, this is mainly a matter of guessing where $color(direction)$ is high. |
613 | 616 |
|
614 | 617 | For a mirror, $s()$ is huge only near one direction, so it matters a lot more. Most renderers in
|
615 |
| -fact make mirrors a special case and just make the $s/p$ implicit-- our code currently does that. |
| 618 | +fact make mirrors a special case and just make the $s/p$ implicit -- our code currently does that. |
616 | 619 |
|
617 | 620 | Let’s do simple refactoring and temporarily remove all materials that aren’t Lambertian. We can use
|
618 | 621 | our Cornell Box scene again, and let’s generate the camera in the function that generates the model:
|
|
654 | 657 | light.
|
655 | 658 |
|
656 | 659 | First, let’s instrument the code so that it explicitly samples some pdf and then normalizes for
|
657 |
| -that. Remember MC basics: $\int f(x) \approx f(r)/p(r)$ . For the Lambertian material, |
| 660 | +that. Remember MC basics: $\int f(x) \approx f(r)/p(r)$. For the Lambertian material, |
658 | 661 | let’s sample like we do now: $p(direction) = cos(\theta) / \pi$.
|
659 | 662 |
|
660 | 663 | We modify the base-class _material_ to enable this importance sampling:
|
|
749 | 752 |
|
750 | 753 | It’s pretty close to our old picture, but there are differences that are not noise. The front of the
|
751 | 754 | tall box is much more uniform in color. So I have the most difficult kind of bug to find in a Monte
|
752 |
| -Carlo program-- a bug that produces a reasonable looking image. And I don’t know if the bug is the |
| 755 | +Carlo program -- a bug that produces a reasonable looking image. And I don’t know if the bug is the |
753 | 756 | first version of the program or the second, or even in both!
|
754 | 757 |
|
755 | 758 | Let’s build some infrastructure to address this.
|
|
783 | 786 |
|
784 | 787 | $$ \phi = 2 \pi \cdot r_1 $$
|
785 | 788 |
|
786 |
| -For theta we have: |
| 789 | +For $\theta$ we have: |
787 | 790 |
|
788 | 791 | $$ r_2 = \int_{0}^{\theta} 2 \pi f(t) \sin(t) $$
|
789 | 792 |
|
|
841 | 844 |
|
842 | 845 | On the plot.ly website you can rotate that around and see that it appears uniform.
|
843 | 846 |
|
844 |
| -Now let’s derive uniform on the hemisphere . The density being uniform on the hemisphere means |
| 847 | +Now let’s derive uniform on the hemisphere. The density being uniform on the hemisphere means |
845 | 848 | $p(direction) = 1 / (2 \pi)$ and just changing the constant in the theta equations yields:
|
846 | 849 |
|
847 | 850 | $$ \cos(\theta) = 1 - r_2 $$
|
|
1139 | 1142 | form a mixture density. Any weighted average of _pdf_ s is a _pdf_. For example, we could just
|
1140 | 1143 | average the two densities:
|
1141 | 1144 |
|
1142 |
| - $$ pdf\_mixture(direction) = \frac{1}{2} pdf\_reflection(direction) + |
| 1145 | + $$ pdf\_mixture(direction) = \frac{1}{2} pdf\_reflection(direction) + |
1143 | 1146 | \frac{1}{2}pdf\_light(direction)
|
1144 |
| - $$ |
| 1147 | + $$ |
1145 | 1148 |
|
1146 | 1149 | How would we instrument our code to do that? There is a very important detail that makes this not
|
1147 | 1150 | quite as easy as one might expect. Choosing the random direction is simple:
|
|
1630 | 1633 | really just sampling a cone uniformly (the cone is tangent to the sphere). Let’s say the code has
|
1631 | 1634 | `theta_max`. Recall from Chapter 9, that to sample $\theta$ we have:
|
1632 | 1635 |
|
1633 |
| - $$ r2 = \int_{0}^{\theta} 2 \pi \cdot f(t) \cdot sin(t) dt $$ |
| 1636 | + $$ r2 = \int_{0}^{\theta} 2 \pi \cdot f(t) \cdot sin(t) dt $$ |
1634 | 1637 |
|
1635 | 1638 | Here $f(t)$ is an as yet uncalculated constant $C$, so:
|
1636 | 1639 |
|
1637 |
| - $$ r2 = \int_{0}^{\theta} 2 \pi \cdot C \cdot sin(t) dt $$ |
| 1640 | + $$ r2 = \int_{0}^{\theta} 2 \pi \cdot C \cdot sin(t) dt $$ |
1638 | 1641 |
|
1639 | 1642 | Doing some algebra/calculus this yields:
|
1640 | 1643 |
|
1641 |
| - $$ r2 = 2 \pi \cdot C \cdot (1-cos(\theta)) $$ |
| 1644 | + $$ r2 = 2 \pi \cdot C \cdot (1-cos(\theta)) $$ |
1642 | 1645 |
|
1643 | 1646 | So
|
1644 | 1647 |
|
1645 |
| - $$ cos(\theta) = 1 - r2/(2 \pi \cdot C) $$ |
| 1648 | + $$ cos(\theta) = 1 - r2/(2 \pi \cdot C) $$ |
1646 | 1649 |
|
1647 | 1650 | We know that for $r2 = 1$ we should get $\theta_{max}$, so we can solve for $C$:
|
1648 | 1651 |
|
1649 |
| - $$ cos(\theta) = 1 + r2 \cdot (cos(\theta_{max})-1) $$ |
| 1652 | + $$ cos(\theta) = 1 + r2 \cdot (cos(\theta_{max})-1) $$ |
1650 | 1653 |
|
1651 | 1654 | $\phi$ we sample like before, so:
|
1652 | 1655 |
|
1653 |
| - $$ z = cos(\theta) = 1 + r2 \cdot (cos(\theta_{max})-1) $$ |
1654 |
| - $$ x = cos(\phi) \cdot sin(\theta) = cos(2 \pi \cdot r1) \cdot \sqrt{1-z^2} $$ |
1655 |
| - $$ y = sin(\phi) \cdot sin(\theta) = sin(2 \pi \cdot r1) \cdot \sqrt{1-z^2} $$ |
| 1656 | + $$ z = cos(\theta) = 1 + r2 \cdot (cos(\theta_{max})-1) $$ |
| 1657 | + $$ x = cos(\phi) \cdot sin(\theta) = cos(2 \pi \cdot r1) \cdot \sqrt{1-z^2} $$ |
| 1658 | + $$ y = sin(\phi) \cdot sin(\theta) = sin(2 \pi \cdot r1) \cdot \sqrt{1-z^2} $$ |
1656 | 1659 |
|
1657 | 1660 | Now what is $\theta_{max}$?
|
1658 | 1661 |
|
1659 | 1662 | 
|
1660 | 1663 |
|
1661 | 1664 | We can see from the figure that $sin(\theta_{max}) = R / length( c - p )$. So:
|
1662 | 1665 |
|
1663 |
| - $$ cos(\theta_{max}) = \sqrt{1- (R / length( c - p ))^2} $$ |
| 1666 | + $$ cos(\theta_{max}) = \sqrt{1- (R / length( c - p ))^2} $$ |
1664 | 1667 |
|
1665 | 1668 | We also need to evaluate the _pdf_ of directions. For directions toward the sphere this is
|
1666 | 1669 | $1/solid\_angle$. What is the solid angle of the sphere? It has something to do with the $C$ above.
|
1667 | 1670 | It, by definition, is the area on the unit sphere, so the integral is
|
1668 | 1671 |
|
1669 |
| - $$ solid\_angle = \int_{0}^{2 \pi} \int_{0}^{\theta_{max}} sin(\theta) d \theta d \phi |
| 1672 | + $$ solid\_angle = \int_{0}^{2 \pi} \int_{0}^{\theta_{max}} sin(\theta) d \theta d \phi |
1670 | 1673 | = 2 \pi \cdot (1-cos(\theta_{max})) $$
|
1671 | 1674 |
|
1672 | 1675 | It’s good to check the math on all such calculations. I usually plug in the extreme cases (thank you
|
1673 |
| -for that concept, Mr. Horton--my high school physics teacher). For a zero radius sphere |
| 1676 | +for that concept, Mr. Horton -- my high school physics teacher). For a zero radius sphere |
1674 | 1677 | $cos(\theta_{max}) = 0$, and that works. For a sphere tangent at $p$, $cos(\theta_{max}) = 0$, and
|
1675 | 1678 | $2 \pi$ is the area of a hemisphere, so that works too.
|
1676 | 1679 |
|
|
0 commit comments