1、Computer GraphicsDr. Norman I. BadlerProfessor, Computer & Information ScienceDirector, SIG Center for Computer Graphics University of Pennsylvania CIS 460 / CIS 560Fall 2015Part D Norman I. Badler(except where otherwise noted)792Rendering via Sampling793Combining Radiosity and Ray Tracing? Radi

2、osity (e.g., by finite element method): accounts for global diffuse illumination: Get area lights, soft shadows, color bleeding, participating media (“fog”). View independent. No arbitrary ambient terms. But expensive meshing and form factor pre-process step. And what about non-uniform, non-diffuse

3、surface properties? Ray tracing accounts for local specular reflections and translucency. View dependent.794Disadvantage to Both: Finally, neither approach can easily account for caustics: bright and shaped areas caused by the focusing and concentration of light by curved (lens-like) surfaces. Leads

4、 to Monte Carlo Path Tracing and Photon Mapping.795The Rendering Equation: Area Formulation(e.g., Finite Elements)2()()()cos(,)cos(,) ()()( ,)( , )erxyeryAxyL xL xL xNNL xL yyx fxV x y dAr Form factor &visibility terms796fr is the Bidirectional Reflectance Distribution Function (BRDF)for the sur

5、face at point xFE solution: Integral reduced to a sum over a finite number of patches.The Rendering Equation: Hemisphere Formulation(e.g., General Irradiance Situation)()()() ()()( ,)cos(,)xererxL xLxLxLxL xfxNd 797fr is the Bidirectional Reflectance Distribution Function (BRDF)for the surface at po

6、int xIntegration via Randomization Algorithm In either of these formulations we have an integral over some domain (hemisphere or surface) that looks quite complicated to solve, since it is inherently recursive. Solution: estimate integral over a finite number of samples. But how to sample? We examin

7、ed the Finite Element approach, but it is costly and often results in poor images. Use a randomized algorithm to sample incoming irradiance. But how to randomize “nicely”? We need to look into this question before we can return to the rendering equation itself.798Las Vegas vs. Monte CarloRandom Samp

8、ling Las Vegas Always gives the same answer. E.g., QuickSort (selecting the pivot element randomly gives O(n log n) complexity, but it always converges to the sorted list whatever the pivot selection method). Monte Carlo Gives the right answer on average. Used on functions that are difficult to eval

9、uate or have no closed forms (such as our recursive integrals!). Slow convergence for N samples. But independent of dimensionality! Artifacts manifest as noise.799N1ONave Integration via Discretization800 NxfNabxxfdxxfINiiNiiba1OError11f(x) Can evaluate f(x), but dont know an analytic way to integra

10、te, e.g. Simple discretization xba Can also sample via Trapezoidal or Gaussian quadrature methods: same general idea. Good in 1D Error is O(1/N) or O(1/N2)Monte Carlo Terminology X, Y, random variables (chosen by a random process). Function of a random variable is a random variable , e.g., Y=f(X). i

11、s a canonical uniform random variable which takes values in domain 0,1) with equal probability. x, y, random values. Pr “Probability of”. P(x): Cumulative Distribution Function (CDF): p(x): Probability Density Function (PDF)801 xXxP Pr dxxdxpPrExample for Uniform Random Variable802Uniform p(x) for 0

12、 110P(x)0 110 11;Pr:CDF10properties PDFPxXxPdxxpxPxPxpxpxpxoD PPdxxpxPxXxPdxxpxPxPxpxpxpxoD,Pr11;Pr:CDF10properties PDFExample for Uniform Random Variable803Uniform p(x) for 0 110P(x)c10 P() P()0 1Expected Value The Expected Value is the mean or average value of the distribution(integral of the func

13、tion times its PDF): E.g., the expected value of a cosine distribution between 0 and : which is what we would expect from its graph.804 dxxpxfxfEDp+ 00sinsin1coscoscosso111since1000dxxdxxpxxExpxpoExpected Value Properties Homogeneity (Linear) Additivity (Linear)805 xfaExafEiiiiXfEXfEMonte Carlo Esti

14、mator (Uniform random variable) Goal: Evaluate Let Xi (i=1,N) be a set of uniform random variables over a,b. The estimator FN of f(x) over the interval a,b is: The PDF p(x) must be a constant and integrate to 1 over a,b. i.e., p(x) =1/(b-a) check: (b-a) p(x) = 1 The expected value EFN of the estimat

15、or FN is indeed the integral:806 dxxfdxfNdxxpxfNabXfENabXfNabEFEbaNibaNibaNiiNiiN x11111NiiNXfNabF1 badxxfVisualizing MC Estimator from 4 Samples807/wjarosz/publications/dissertation/appendixA.pdfMonte Carlo Estimator (Arbitrary random variable) If X is a random variable dr

16、awn from an arbitrary PDF p(x), then the estimator is: (and p(Xi) 0 when |f(Xi)|0) Choosing a “good” PDF is important. (Imagine p(xi) are “weights” that integrate to 1.)808NiiiNXpXfNF11 baNibaNibaNiiiNiiiNdxxfdxxfNdxxpxpxfNXpXfENXpXfNEFE11111111Monte Carlo Estimator (Multiple Dimensions) Suppose tha

17、t N samples Xi are drawn from an arbitrary multidimensional (“joint”) PDF p(X), and the 3D integral is: If samples are chosen uniformly from the box (x0, y0, z0) to (x1, y1, z1), then PDF p(X) is the constant And the estimator is809dzdydxzyxfxxyyzz101010,010101111zzyyxxiiXfNzzyyxx010101Monte Carlo P

18、roperties The number of samples can be chosen independently of the dimensionality of the integrand. (Standard numerical quadrature is exponential in the dimension!) The integral domain need not be a Cartesian box, but any domain (e.g., a hemisphere) as long as we can establish a suitable PDF. Monte

19、Carlo convergence rate is . (Standard quadrature converges faster in 1D, but exponentially worse in higher dimensions.) Monte Carlo is therefore the only practical numerical integration technique for high dimensional integrals! Moreover, in computer graphics rendering using path tracing, the integra

20、nds are infinite dimensional!810N1OSo what do we need for Monte Carlo Integration? N samples Xi and an “appropriate” PDF p(Xi): We need to be able to evaluate f at any Xi. Generating Xi samples: methods include Discrete bins Inversion Rejection Metropolis (not done here) Assume uniform random variab

21、le 0,1) generated by a pseudo-uniform random number generator.811NiiiNXpXfNF11 1Toy Example: Discrete BinsGenerate the following PDF using (note ):(each pi is uniform)812p1p2p3p4 1 2 3 4 20% 30% 40% 10%p1p1p1p1p2p2p2p3p3p4 141xpiiStack PDF to get CDF: ijiixpxP1 1Toy Example: Discrete Bins813Suppose

22、0.9,1), then are in p4 Suppose 0.5,0.9), then are in p3 Suppose 0.2,0.5), then are in p2 Suppose 0,0.2), then are in p1 p1p1p1p1p2p2p2p3p3p4 =0.96 =0.71 =0.32 =0.11Generate a value for ; then use selected pi to generate X: X is generated uniformly within the pi interval.CDFXXXX 1Inversion Method (If

23、 CDF is continuous and invertible) Get PDF p(x) Compute CDF Compute inverse P-1(x) Get Compute X=P-1()814 =0.31XPDFCDF xdxpxPx0p(x) P(x) Example: Power Distribution (arises from microfacet model)Sampling from the power distribution815 nxxp 111110110Thus)(and:Hence1so1011, 1SincePDFnnnnnnXxxPxdxxpxPn

24、cnxcdxcxcxxpThe Rejection Method (Throw darts!) Cant integrate f(x) to find its PDF? Cant analytically invert its CDF? p(x) is “ugly” or “weird” or just complicated? Then throw darts at it! Suppose f(x) (c p(x) for some constant scalar c and we know how to sample from p(x).loop forever (i.e., until

25、“happy”): sample X from ps distribution generate if f(X) / (c p(X) thenreturn X Works for any number of dimensions.816c*p(x) f(x) x221x1Examples of Rejection Sampling Uniformly sample a disk (embedded in a unit square) by rejection sampling. (Result IS uniform but “clumpy” in undesirable ways. Recal

26、l stratified sampling discussion.)817Direct Approach to Disk Sampling Compute samples uniformly over a unit disk in polar coordinates. An obvious sampling strategy is but this causes samples to “clump” toward the disk center (a). Without proof, a good uniform sampling function for the unit disk is (

27、b): (The counteracts the disk “clumping”.)818212 ,r1212 ,r212 ,r212 ,r(a)(b)Monte Carlo Integration: Sampling over a Hemisphere819 We need to sample irradiance due to a light source over a hemisphere: How to do this? What might be good choices for p? NiipLNIEddLdLI1200sincos1sincoscos2Sampling over

28、Spheres and Hemispheres820 Uniform direction samples over a complete sphere: And uniform samples over a hemisphere:1112221122221122sin12sin122cos12coszzyzx1212212cos12sinsinsin12coscossinzyxCosine-Weighted Hemisphere Sampling Recall that for irradiance, incident energy is “more important” than grazi

29、ng irradiance: i.e., energy importance follows a cosine distribution, just like Lamberts Law. So if we have to sample a hemisphere around a surface point, it makes sense that we should bias the samples in a cosine-weighted fashion: more samples near the surface normal. We can find these sample direc

30、tions using “Malleys method” (no, not Adams ).821Malleys Method The PDF we want is proportional to the cosine: Normalizing So Finally in terms of solid angle: Note that this is non-negative over 0,/2, and we insured it integrates to 1.822 cosp 10coscos102cos22sin21sincos21sincos121212121221202020202

31、 ccccdcdcddcdpc1c cos1pProjection into Hemisphere Uniformly sample points on the unit disk and just project them up onto the unit hemisphere. Vertical direction: Circular direction:823H2(N)212 ,rsamples21; 00sinrrrdiskon Forming the Sample Ray To form sample ray: Have surface point P, Hemisphere sam

32、ple point is S=(x,y,z) where (x,y) comes from (r,) sample: And 824H2(N), r221yxzsin;cosryrxSPBack to Monte Carlo Path TracingSummary1. Choose sample paths according to an appropriate Probability Density Function.2. Evaluate function at sample.3. Average these appropriately weighted values using a su

33、itable estimator.Error slow convergence (twice as good need 4 times more samples).Variance in estimator noise in image.Good for high dimensional functions.Theoretically best we can do (good estimator and more samples better image).825N1OMonte Carlo Path Tracing For each pixel, estimate radiance valu

34、e as (appropriately chosen and weighted average of (lots of) random paths along rays from eye through the pixel. Compute radiance value for each visible point from Rendering Equation: Self Emission Direct Illumination (shadow rays) Indirect Illumination (Recursive) Keep lengths of recursive paths ma

35、nageable Russian Roulette.826Cutler & Durand, MITMore Samples = Higher Cost but Less Noise!8278 samples/pixel 40 samples/pixel 200 samples/pixel13 sec 63 sec 5 min1K samples/pixel 5K samples/pixel 25K samples/pixel25 min 124 min 10.3 hrs2.4 GHz Intel Core 2 Quad CPU using 4 threads. See the Radi

36、ance Arriving at the Eye828pixel-wise sampling filter dppheyexLLimageplanepixelreflected lightemittedxLxLxLeyexLreestimator for reflected light NiixiirirpNxfxLNxLE1,cos,1Recursive step: form ray from x in direction i and find irradiance from that directioniiixrLxL,How to Stop Recursion? Typically us

37、e “Russian Roulette” Unbiased Absorption probability P: =1P is proportional to hemispherical reflectance Check: black non-reflective surface, = 0, P = 1, so recursion must stop here; mirror (specular) surface, = 1, P = 0, so recursion cannot stop here. Threshold (fixed or adaptive)829Sampling Le(x )

38、Direct Illumination: First Option1.Sample each light source separatelyG(x, yi) is microfacet term in direction yix.V(x, yi) is visibility (“shadow ray”) between x and yi .Possible PDFs:Uniform light area samplingUniform sampling of solid angle subtended by light source830 sNiiiiiriiesdirectypyxVyxGx

39、yxfxyyLNxLE1,1Sampling Le(x )Direct Illumination: Second Option2.Combine light sources as one integral domainG(x, yi) is microfacet term in direction yix.V(x, yi) is visibility (“shadow ray”) between x and yi .PDFs:Uniform source selection; uniform sampling of light area.Power proportional source se

40、lection; uniform sampling of light area.Works well for “natural” light scenes (non-discrete sources), e.g., an (HDR) light field map or “baked-on” light “textures”.831 NiiiiLiiiriiedirectkypkpyxVyxGxyxfxyyLNxLE1,1Sampling 2D Piecewise-Constant Functions(Used for texture maps and environment maps)(a)

41、 The original HDR environment map u,v.(b) A low-resolution version of the row-wise v density function.(c) The piece-wise constant distributions for rows of the image.First (b) is used to select a row v of the image to sample. Rows with bright pixels are more likely to be sampled.Then, given a row, a

42、 value u is sampled from that rows 1D distribution. (HDR map from Paul Debevec.)832(a)(b) (c)uvSampling Lr(x ) Indirect Illumination Uniform hemisphere sampling Importance sampling (PDF choice) of any of the following: Cosine BRDF Incident radiance Combinations of any of the above: Multiple Importan

43、ce Sampling833Multiple Importance Sampling Variance dependent upon the choice of PDF of the estimator. Importance sampling the integral means taking samples that are likely to contribute most significantly to the result. (Reminds one of progressive radiosity method.) IrradianceBRDFcosine834dNxfxrLxL

44、xrrindirectx,cos,Importance Sampling: Light Source(s) Uniform Light Source Sampling Power based Light Source Sampling 1 path per pixel 1 shadow ray per point Direct Illumination Only835Importance Sampling: Cosine Uniform Sampling of Hemisphere Importance Sampling (cosine) of Hemisphere 49 paths per

45、pixel 1 shadow ray per point Max eye path depth = 10 Direct + Indirect Illumination836Importance Sampling: BRDFUniform Sampling of BRDF Importance Sampling of BRDF49 paths per pixel 1 shadow ray per point Max eye path depth = 10 Direct + Indirect IlluminationUniformly sampling a directionalBRDF prod

46、uces noise837MC Path Tracing Difficult Case Monte Carlo Path Tracing is a general technique, but it can do poorly in some scenes, especially with small areas of direct illumination and large areas lit by indirect illumination.838Rarely sampled lit area on ceiling.Shadow feelers almost always miss hi

47、tting light.Image will have high variance: mostly dark except for the few “l(fā)ucky” paths that hit the ceiling area.Bi-Directional Path Tracing Trace paths forward from the camera (as in MC Path Tracing) BUT ALSO backward from the lights. Path from camera: p1, p2, , pi Path from light: q1, q2, , qj Fo

48、r suitable i and j. Then, by concatenation, a path from camera to light is p1, , pi , qj, , q1 A shadow feeler ray is traced from pi to qj to see if they are mutually visible. Can re-use these paths by using combinations of prefixes: paths range from 2 to i+j vertices long, e.g.:p1, q1, p1, p2 , q3,

49、 q2, q1, p1, p2 , q1, etc.839Bi-Directional Path Tracing Showing a couple of possible paths and one blocked path. Various other optimizations possible.840p1q1q2p2p1, q2, q1p1, p2, q2, q1q3p1, p2, q3, q2, q1 Temporal Anti-Aliasing - Motion Blur Blur image of moving object to combat high frequency jer

50、kiness or strobing caused by frame-rate sampling. Earliest work in 1984: Korein and Badler did spheres and polygons. Lucasfilm introduced spatial and temporal “noise” into ray-tracing: Monte Carlo or distribution path tracing. Also applied to camera depth-of-field841http:/ Carlo Path Tracing and Tem

51、poral Anti-Aliasing Extend Monte Carlo path trace into time domain. Because a small amount of noise is added, it is best to supersample and filter. Can also do camera depth-of-field.842Monte Carlo Path Tracing in Time/Camera SpacelightsurfacelensimageplaneFor each primary ray:1. Jitter spatial locat

52、ion on image plane within the pixel.samplepoint843Monte Carlo Path Tracing in Time/Camera SpacelightsurfacelensimageplaneFor each primary ray:2. Jitter time T within inter-frame interval.samplepoint844Monte Carlo Path Tracing in Time/Camera SpacelightsurfacelensimageplaneFor each primary ray:3. Move

53、 camera and objects to position at that time T.samplepoint845Monte Carlo Path Tracing in Time/Camera SpacelightsurfacefocalpointlensimageplaneFor each primary ray:4. Determine focal point.samplepoint846Monte Carlo Path Tracing in Time/Camera SpacelightsurfacefocalpointlensimageplaneFor each primary ray:5. Jitter lens loca