K. H. Ko School of Mechatronics Gwangju Institute of Science and Technology

1. Photo-realistic Rendering and Global Illumination in Computer Graphics Spring 2012Stochastic Radiosity K. H. Ko School of Mechatronics Gwangju Institute of Science and Technology

2. Discrete Random Walk Methods for Radiosity • The methods belonging to this class are based on the concept of a random walk in a so-called discrete state space. • It turns out that these algorithms are not better than the stochastic Jacobi algorithm.

3. Random Walks in a Discrete State Space • Consider one experiment, involving a set of n urns, labeled i, i=1,…,n. • One of the urns contains a ball, subject to the following “game of chance”: • They ball is initially inserted in a randomly chosen urn. The probability that the ball is stored in urn i, I = 1, …, n, is πi. These probabilities are properly normalized: Σπi = 1. They are called source or birth probabilities. • The ball is randomly moved from one urn to another. The probability pij of moving the ball from urn i to urn j is called the transition probability. The transition probabilities from a fixed urn i need not sum to one. If the ball is in urn i, then the game will be terminated with probability αi=1-Σj=1pij. αi is called the termination or absorption probability at urn i. The sum of the transition probabilities and the termination probability for any given urn is equal to unity. • The previous step is repeated until termination is sampled.

4. Random Walks in a Discrete State Space • Suppose the game is played N times. • During the games, a tally is kept of how many times each urn i is visited by the ball. • The expected number of times Ci that the ball will be observed in each urn i becomes

5. Random Walks in a Discrete State Space • The first term on the right-hand side of the equation indicates the expected number of times that a ball is initially inserted in urn i. • The second term indicates the expected number of times that a ball is moved to urn i from another urn j.

6. Random Walks in a Discrete State Space • The urns are called states and the ball is called a particle. • The game of chance outlined before is an example of a discrete random walk process. • Discrete: because the set of states is countable. • The expected number of visits Ci per random walk is called the collision density χi. • The collision density of a discrete random walk process with source probabilities πi and transition probabilities pij is the solution of a linear system of equations:

7. Random Walks in a Discrete State Space • Note that χi can be larger than unity. • χi is called the collision density rather than a probability. • In summary, at least a certain class of linear systems of equations can be solved by simulating random walks and keeping count of how often each state is being visited. • The states of the random walk correspond to the unknowns of the system.

8. Shooting Random Walk Methods for Radiosity • The power system is similar to the equation

9. Shooting Random Walk Methods for Radiosity • However, the source terms Pei in the system do not sum to one. • The remedy is to divide both sides of the equations by the total self-emitted power PeT = Σi Pei:

10. Shooting Random Walk Methods for Radiosity • This system of equations suggests a discrete random walk process with: • Birth probabilities πi = Pei/PeT: particles are generated randomly on light sources, with a probability proportional to the self-emitted power of each light source. • Transition probabilities pij = Fijρj: first, a candidate transition is sampled by tracing, for instance, a local line. After candidate transition, the particle is subjected to an acceptance/rejection test with survival probability equal to the reflectivity ρj. If the particle does not survive the test, it is said to be absorbed.

11. Shooting Random Walk Methods for Radiosity • By simulating N random walks in this way, and keeping a count Ci of random walk visits to each patch i, the light power Pi can be estimated as • Because the simulated particles originate at the light sources, this random walk method for radiosity is called a shooting random walk method. • It is called a survival random walk estimator because particles are only counted if they survive the rejection test.

12. Shooting Random Walk Methods for Radiosity

13. Shooting Random Walk Methods for Radiosity • Collision Estimation • Transition sampling is suboptimal. • Candidate transition sampling involves an expensive ray-shooting operation. • If the candidate transition is not accepted, this expensive operation has been performed in vain. • This will often be the case if a dark surface is hit. • It will always be more efficient to count the particles visiting a patch, whether they survive or not. • The collision random walk estimates are C’i denotes the total number of particles hitting patch i. The expected number of particles that survive on i is ρiC’i≈Ci.

14. Shooting Random Walk Methods for Radiosity • Absorption Estimation • A third random walk estimator only counts particles if they are absorbed. • The resulting absorption random walks estimates are • C’’i : the number of particles that are absorbed on i. • C’i = Ci + C’’i • The expected number of particles being absorbed on i is (1-ρi)C’i ≈C”i.

15. Gathering Random Walk Methods for Radiosity • It is possible to estimate the radiosity on a given patch i, by means of particles that originate at i and that are counted when they hit a light source. • Gathering random walk estimators • A gathering random walk estimator corresponds to a shooting random walk estimator for solving an adjoint system of equations.

16. Gathering Random Walk Methods for Radiosity • Adjoint Systems of Equations • Consider a linear system of equations Cx=e. • C is the coefficient matrix of the system with elements cij. • e is the source vector. • x is the vector of unknowns. • A well-known result from algebra states that each scalar product <x,w>=Σni=1xiwi of the solution x of the linear system, with an arbitrary weight vector w can also be obtained as a scalar product <e,y> of the source term e with the solution of the adjoint system of linear equations CTy = w. • <w,x> = <CTy,x> = <y,Cx> = <y,e>

17. Gathering Random Walk Methods for Radiosity • Adjoint of the radiosity system • Adjoint systems corresponding to the radiosity system of equations look like • Consider the power Pk emitted by a patch k. • Pk can be written as a scalar product Pk = AkBk = <B,W> with Wi=Aiδik. • All components of the direct importance vector W are 0, except the kth component, which is equal to Wk = Ak. • Then, Pk can also be obtained as Pk = <Y,E> = ΣiYiBei. • It is a weighted sum of the self-emitted radiosities at the light sources in the scene. • The solution Y of the adjoint system indicates to what extent each light source contributes to the radiosity at k. • Y is called the importance or potential in the literature.

18. Gathering Random Walk Methods for Radiosity • Gathering random walk estimators • The adjoints of the radiosity system also have the indices of the form factors in the right order. • They can be solved using a random walk simulation with transitions sampled with local or global lines. • The particles are now shot from the patch of interest, instead of from the light sources. • The Transition probabilities are pji = ρjFji. • First, an absorption/survival test is performed. • If the particle survives, it is propagated to a new patch with probabilities corresponding to the form factors. • A nonzero contribution to the radiosity of patch k results whenever the imaginary particle hits a light source. • Its physical interpretation is that of gathering. • The gathering random walk estimator is a collision estimator.

19. Discussion • Discrete random walk estimators for radiosity can be classified according to the following criteria: • Whether they are shooting or gathering. • According to where they generate a contribution: at absorption, survival, at every collision. • Comparison of the variants can be made by considering the variance of each method.

20. Discussion • Shooting versus Gathering • The shooting estimators have lower variance, except on small patches which have low probability of being hit by rays shot from light sources. • Unlike shooting estimators, the variance of gathering estimators does not depend on the patch area Ak. • For sufficiently small patches, gathering will be more efficient. • Gathering could be used in order to “clean” noisy artifacts on small patches after shooting.

21. Discussion • Absorption, Survival or Collision • The survival estimators are always worse than the corresponding collision estimators. • The reflectivity ρk (shooting) or ρs (gathering) is always smaller than 1. • As a rule, the collision estimators also have lower variance than the absorption estimators. • When the transition probabilities are modulated, to shoot more rays into important directions, an absorption estimation can be better than a collision estimator. • It can be shown that a collision estimator can never be perfect, because random walks can contribute a variable number of scores. • An absorption estimator always yields a single score, so it does not suffer from this source of variance. • In theory, absorption estimators can be made perfect.

22. Discussion • Discrete Collision Shooting Random Walks versus Stochastic Jacobi Relaxation • For the same number of rays, discrete collision shooting random walks and incremental power-shooting Jacobi iterations are approximately equally efficient. • However, this conclusion is no longer true when higher-order approximations are used, or with low-discrepancy sampling or in combination with variance reduction techniques. • Many variance-reduction techniques and low-discrepancy samplings are easier to implement and appear more effective for stochastic relaxation than with random walks.

23. Photon Density Estimation Methods • The radiosity system of linear equations can be solved stochastically. • By sampling according to the form factors, the numerical value of the form factors was never needed. • Instead of solving the radiosity system of equations, the general rendering equation or the radiosity integral equation can be solved. • Random walks in a continuous state space can be used. • The random walks are nothing but simulated trajectories of photons emitted by light sources and bouncing throughout a scene. • The surface hit points of these photons are recorded in a data structure. • An essential property of such particle hit points is that their density at any given location is proportional to the radiosity at that location.

24. Photon Density Estimation Methods • The main benefit of this approach is that nondiffuse light emission and scattering can be taken into account to a certain extent. • The rendering equation does not have to be solved exactly at every surface point. • The photon density estimation methods open the way to more sophisticated and pleasing representations of the illumination than the average radiosity on surface patches. • In the photon mapping method, the representation of illumination is independent of scene geometry. • This allows us to use nonpolygonized geometry representations, procedural geometry such as fractals, and object instantiation in a straightforward manner.

25. Photon Transport Simulation and Radiosity • Photon trajectory simulation is called analog photon trajectory simulation. • We start out by selecting an initial particle location x0 on a light source. • We do that with a probability proportional to the self-emitted radiosity: • x0 is the starting point for a random walk. • S(x0) is called the birth or source density.

26. Photon Transport Simulation and Radiosity • Next, an initial particle direction Θ0 is selected using the directional light emission distribution of the light source at x0, times the outgoing cosine. • For a diffuse light source,

27. Photon Transport Simulation and Radiosity • Consider now a ray shot from the sampled location x0 into the selected direction Θ0. The density of hit point x1 of such rays with object surfaces depends on surface orientation, distance, and visibility with regard to x0. The transparent surface shows this density on the bottom surface of the shown model in the right figure.

28. Photon Transport Simulation and Radiosity • The density of incoming hits is shown, taking into account these geometric factors as well as the (diffuse) light emission characteristics at x0.

29. Photon Transport Simulation and Radiosity • A survival test is carried out at the obtained surface hit point x1: • A random decision is made whether or not to sample absorption (and path termination) or reflection. • We take the probability σ(x1) of sampling reflection equal to the albedo ρ(x1,-Θ0), the fraction of power coming in from x0 that gets reflected at x1. • The full transition density from x0 to x1 is thus

30. Photon Transport Simulation and Radiosity • A survival test is carried out at the obtained surface hit point x1:

31. Photon Transport Simulation and Radiosity • If survival is sampled, a reflected particle direction is chosen according to the BRDF times the outgoing cosine. • For a diffuse surface, again only the outgoing cosine remains.

32. Photon Transport Simulation and Radiosity • Subsequent transitions are sampled in the same way, by shooting a ray, performing a survival test, and sampling reflection if the particle is not absorbed. The influence of surface orientation, distance, and visibility on the left surface of the scene with regard to x1.

33. Photon Transport Simulation and Radiosity • Subsequent transitions are sampled in the same way, by shooting a ray, performing a survival test, and sampling reflection if the particle is not absorbed. The combined effect of the cosine distribution at x1 and the former.

34. Photon Transport Simulation and Radiosity • The following figures indicate that the process has been performed twice more, this time for nondiffuse reflection.

35. Photon Transport Simulation and Radiosity • Now, consider the expected number χ(x) of particle hits resulting from such a simulation, per unit of area near a surface location x. • This particle hit density consists of two contributions: • The density of particles being born near x as given by S(x) • The density of particles visiting x after visiting some other surface location y. • The density of particles coming from elsewhere depends on the density χ(y) elsewhere and the transition density T(x|y) to x:

36. Photon Transport Simulation and Radiosity • For a diffuse environment, • In other words, the number of particle hits per unit area expected near a surface location x is proportional to the radiosity B(x).

37. Photon Transport Simulation and Radiosity • The problem of computing radiosity has thus been reduced to the problem of estimating particle hit densities. • To estimating the number of particle hits per unit area at a given surface location. • The problem of estimating a density like that, given nothing more than a set of sample point locations, has been studied intensively in statistics. • An alternative, equivalent point of view is to regard the problem at hand as a Monte Carlo integration problem. • We want to compute integrals of the unknown radiosity function B(x) with a given measurement, or response, function M(x)

38. Photon Transport Simulation and Radiosity • The radiosity B(x), and thus the integrand, cannot be evaluated a priori, but since analog simulation yields surface points xs with density χ(xs)=B(xs)/PeT, M can be estimated as

39. Photon Transport Simulation and Radiosity • Basically all we have to do is simulate a number of photon trajectories and accumulate the value of measurement function M(xs) at the photon surface hit points xs. • The procedure corresponds closely with the survival estimator. • Particles are only taken into account after surviving impact on a surface. • Absorption and collision estimators can be defined and source term estimation can be suppressed. • In practice, collision estimation, that is, counting all particles that land on a surface, is preferred. • This “over-counting” shall be compensated by multiplying all resulting expressions by the reflectivity.

40. Photon Transport Simulation and Radiosity • Histogram Methods • The easiest, and probably most often used, way of estimating density functions is by subdividing the domain of the samples into bins-surface patches in our case- and to count the number of samples Ni in each bin. • The ratio Ni/Ai yields an approximation for the particle density in each bin.

41. Photon Transport Simulation and Radiosity • Histogram Methods

42. Photon Transport Simulation and Radiosity • Orthogonal Series Estimation • Histogram methods yield a single average radiosity value for each surface patch. • It is possible to obtain linear, quadratic, cubic, or other higher order approximations for the radiosity function B(x). • The problem of computing such higher-order approximations comes down to computing the coefficients Bi,α in the following decomposition of B(x)

43. Photon Transport Simulation and Radiosity • Orthogonal Series Estimation • A constant approximation is obtained when using just one basis function Ψi(x) per patch. • The value of the function is 1 one the patch and 0 outside. • In this case, the method corresponds to the histogram method. • The main advantage of orthogonal series estimation over the histogram method is that a smoother approximation of radiosity is possible on a fixed mesh. • The main disadvantage is the cost. • The cost of computing a higher-order approximation with K basis functions to fixed statistical error is about K times the cost of computing a constant approximation.

44. Photon Transport Simulation and Radiosity • Kernel Methods • The radiosity B(z) at a point z could also be written as an integral involving a Dirac impulse function. • Estimation of the integral is not an easy thing. • An approximation for the radiosity at z can be obtained by using a different, normalized density kernel or footprint function F(x,z) with nonzero width centered around z:

45. Photon Transport Simulation and Radiosity • Kernel Methods • Usually, a symmetric kernel is chosen, which depends only on the distance between x and z: F(x,z) = F(z,x) = F(rxz). • With N photon trajectories, the integral can be estimated as • Advantage • They allow a representation of illumination that is not necessarily mesh-based, or that allows a suitable mesh to be constructed a posteriori. • This allows us to get rid of edge discontinuities. • Disadvantage • The choice of the kernel bandwidth is a difficult problem. • A lot of effort is required in order to avoid underestimation oat surface dges, called boundary bias.

46. Photon Transport Simulation and Radiosity • Nearest Neighbor Methods • Rather than fixing a certain area A and counting the number of particles N on that area like in the histogram method, one fixes a number of particles N and looks for an area that contains this number of particles. • If N particles can be found on a small area, the density (N/A) will be high. • If they are only found on a large area, the density will be low. • The main advantage of nearest neighbor estimation is that one can entirely get rid of surface meshes. • All one needs to store is the set of particle hit-point locations.

47. Photon Transport Simulation and Radiosity • Nearest Neighbor Methods

48. Photon Transport Simulation and Radiosity • Instant Radiosity • The right-hand side of the radiosity integral equation is also an integral containing the radiosity function. • This integral can be estimated by tracing a number of photon trajectories. • This results in estimates for B(z) of the following form:

49. Photon Transport Simulation and Radiosity • Instant Radiosity • The sum in the right-hand side can be interpreted as the direct illumination at z due to point light sources with strength PeT/N placed at every photon hit point xs.