Ambient Occlusion

try and correct it again. For the volume in scattering, we have, we have n photons in the circle, we have, I mean sphere, delta, phi, j, the flux of the photon, the albedo, the phase function with the two angles in the numerator, and in the denominator we have four-thirds pi r cubed, which is the radius of the sphere, and I put a one over n here for the average, and that was the error on the board. We don't find the average, we add up all the flux, and then we divide by the volume of the sphere to find out the flux that's in scattered per unit volume, so it's the total flux, I shouldn't have had the one over n that I wrote on the board. It wasn't in my notes, but somehow I added it in the lecture thinking I needed an average here, but I don't, and the reason you need the albedo in here is because the phase function is a probability. I'm not sure in the book it might have been an F instead of a P, and it integrates to one, and so alpha is the albedo, which says, you know, what fraction of the light that actually is represented by all those photon scattering events actually does scatter rather than be absorbed, because you store a photon at every extinction event where it's either absorbed or scattered. And now this thing could, in black and white, it's just the flux per photon, but in color, instead of doing a complete separate photon propagation pass and a separate photon map for the three colors, what you could do is think of alpha as a vector, alpha R, alpha G, alpha B, and you multiply this flux separately as it's propagating by these three numbers in case they're different. That's assuming that the extinction is independent of the wavelength, it's only the scattering that is. And what I proposed earlier is to use this alpha factor as a Russian roulette probability of continuing. But what you could do at each scattering event is you could let alpha max is the max of the albedo for red, the albedo for green and the albedo for blue. And then what you could do is multiply the flux by the ratio alpha R to alpha max. If the maximum was R, this ratio would be one. So for one of these, the ratio is is one and the other two at some number less than or equal to one. And you could still use the alpha max to determine the continuing probability in the Russian roulette. And then if you didn't divide by it, you're you're basically dividing by it here. You know, you didn't divide by it in the denominator here in the scattering computation. OK, so now I want to talk about precomputed radiance transfer. But I noticed on the page right before it, there was ambient occlusion, which is worth discussing. That could also be determined by ray tracing. And what it is, is the fraction of the distant light. I know I say ambient light from all directions, not occluded by the object. So we can we can compute the ambient occlusion for each object in the scene separately as a precomputation, maybe store the texture and then get the ambient component from it. So, for example, if you have like a fold like this and a piece of paper, it's going to be darker inside the fold because, you know, less of the ambient light from all over the room gets in the fold, especially if this paper wasn't reflecting so much light because it's white. So basically what this is, is if you integrate over the whole hemisphere of the cosine, you're going to get pi. So to turn into a fraction between zero and one, we're going to integrate over this hemisphere of the visibility of distant illumination. That means trace a ray from a point and see if it hits the object. And if it does, this will be zero. Otherwise it'll be one. And then you assume that you have that ambient illumination coming in from that direction, omega d omega. So this is the unoccluded fraction of the hemisphere weighted by this cosine factor. I should have put cosine or I should have just put dot. Let me follow the book by saying cosine at the point in question. All right. So the cosine weighted average visibility of some sort of distant environment representing ambient light coming in from all directions. So if you multiply that ambient light by this factor, you'll get the unoccluded light and you can use that in shading. OK, so now on to precomputed radiance transfer, which is the real topic for today. So I'm following the book's notation, if I can remember to. But it's basically taken from the original paper by Peter Pike Sloan et al. That's listed in the index of the book. And the idea is we're going to represent the two factors in our integral. I don't know if I have that integral yet. I'll put it in a minute. It's basically the BRDF psi, I guess, in the book times cosine and psi. And let's see, I guess I'm going to put that in into the I'm going to put visibility also. This is visibility to a distant environment. We're going to look again at, say, an environment on a sphere at infinity or on the six faces of a cube of possible directions. And it's take a small object and include it in a given environmental illumination. And we're trying to find out the environmental illumination, which means we have to get this integral. The two factors in the integral are the local factor, depending on the geometry and the environmental illumination at infinity, given, I guess, if I'm going to follow the notation by in the book, I should call it psi. So these are two factors and the integral we get. Why don't I call it L zero? I guess it's in the book X going to get with coming from direction. Theta is the integral of. Why don't I call this factor one? So I don't have to write it all again. Times the environment. This one depends on psi. And so does that one. The omega of psi. So each of these factors are supposed to be represented using orthogonal basis functions. Why don't I say using orthonormal basis functions or though normal Y sub I of angle. So what does that mean? Orthonormal means if we take the dot product of two of them and in our context, the dot product Y J means the integral of over the hemisphere of Y I psi times Y J psi the omega psi. That's our definition of inner product in this Hilbert space of functions. And what this should be is the same as delta I J where delta I J equals one if I equals J and zero otherwise. So they're orthogonal, which means if you integrate any two different ones, you'll get zero. And they're orthonormal, which means that if we integrate the product of one with itself, namely the thing squared, then we'll get one. So we're looking at a function space of continuous functions over the unit sphere. And we're turning it into a space with a dot product by integrating the product of the two functions. We're basically making the it's called a L2 function space because the norm of a function in that function space is its square of its integral. The integral of the square of the function. OK, so what we're going to do then is we're going to let factor one. I just want to say face writing it again. I wrote it out in my notes of psi equals the sum of coefficients. T.I. times orthonormal basis functions of the angle and factor two. Maybe I can afford to and write out L.I.Y.I of psi. So basically, if we have a complete set of orthonormal functions in our basis, that means that as you sum. I probably should have done I equals one to end here. You take more and more terms. This is an approximation. But in the limit, the error will go to zero in the L2 sense, because the sum of the infinite series should converge to the true function. So if I have time at the end of the period, I may say something about the wavelet basis function. How many people have heard of wavelets? Har wavelets? Only one. OK, these are alternatives to what I'm mainly going to discuss today, which is spherical harmonic basis functions. So spherical harmonic basis functions are good for representing smooth functions on the sphere without sharp discontinuities. And with sort of low spatial frequency is measured on the area of the sphere or basic angular frequency, whereas har wavelets might be better for, you know, isolated small light sources in the environment or sharp boundaries of an area light source or something. How many people have heard of spherical harmonic? Nobody. How many people have heard of Fourier series? OK, well, let's think of it in terms of the context of Fourier series. So in a Fourier series for a function of theta on the units circle. So if we think of F of theta, let me use lowercase theta. We can represent it as a zero plus a one cosine theta plus B one sine theta plus a two cosine two theta plus B two sine two theta and so on. For higher order terms. Maybe I have a more concise representation here. I guess I don't know. We'll just write it as this is approximated by a finite series or we can think of the infinite series converges some and equals. Well, let's add a zero because that's N equals zero. That's the constant term sort of in terms of functions of time. That's sort of the DC part plus some and equals one to infinity of a and cosine and theta plus B and sine and theta. And another way of writing it is a zero plus the sum and equals one to infinity of some complex constant, which has got real part here and imaginary part here times e to the I and theta. Does this look familiar to. OK, so we want to do the same thing for functions on a sphere instead of a circle. We want this is cosine and sine and just a constant. One say are an orthonormal basis for functions, periodic functions on the circle in the sense that if you integrate sine squared or cosine squared, you'll get one over the unit circle. And if you take sine times cosine or any number of sines times sine and theta times cosine and theta or something that will integrate out to zero. And you can represent any function in the limit as infinitely many terms. So that's Fourier series. So spherical harmonics are the analogous basis functions on the sphere. So one thing I should say about Fourier series is if you want a representation with a small number of terms, then your function has to be fairly smooth. It can't have high frequencies beyond the number of terms you want to use. And it's the same in spherical harmonics. So I'll give you equivalent formulas, both in complex form and real form. So basically when the unit sphere and we'll parameterize it by theta, which is the angle of unit vector on the unit sphere compared to the Z axis. Right. And C is the angle that it's projection on the X Y equatorial plane makes with the X axis. These are the spherical coordinates we've been using all along. And then Y lowercase L lower subscript L upper subscript M of theta and phi is some normalization constant K. L M, which you need to make the square integral of one of these functions, which is with itself add up to one instead of some other complicated formula. In fact, the formula for KLM, I'll write in a minute, e to the I M phi. So in terms of their variation as a function of phi, they just look like Fourier series. But they also vary as a function of theta by P L M. Actually, you only need absolute value. And they're the same for plus and minus of cosine theta. So in terms of the X, Y, Z coordinates, cosine theta will be just the Z coordinate. And it turns out that these polynomials in phi or gives sort of give you the X and Y coordinates of the point. It turns out that. Let's see. Actually, these are complex. I'll write the real form in a minute. So let me just tell you what these normalization factors are. First of all, for L 0 1 up to some number where we're going to stop. But these could go on to infinity. And for each of those else M can go from minus L up to plus L. So for, say, zero, there's only one term for one. There's three terms minus one zero and one. So that's four. And for two, there'll be five more terms. Minus two, minus one, zero, one, two. So in general and square terms. Right. Because you're adding up the first and odd numbers to get the total number of terms when L goes up to N and copy the normalization factor. K L M equals the square root of two L plus one over four pi L minus the absolute value of M factorial L plus the absolute value of M factorial. Now, I didn't tell you what these are, except they're they're polynomials and they're called associated Legendre polynomials. And there is a recursion formula for the accuracy of the recursion formula for this whole expression. So let me give you the analogous real real version, which is easier to think of. And that's usually denoted with both subscripts down. And that's going to be twice the square root of two times the real part of Y L M, which is going to be the square root of two times this K L M times cosine M C, because this is the E to the IMC is a complex number. And then this polynomial P L M cosine theta. And that's for M greater than zero. And it's the square root of two times the imaginary part of Y L M L M. And then we have sine M C P L M cosine theta. And that's for M less than zero. And for M zero, it's just Y L zero, which is K L zero P L zero cosine theta for M equals zero. Because if we have M equals zero, this is e to the zero power, which is one. And we're not getting any imaginary part anyway. Just like in the Fourier series where I made a special term for that. So it turns out these represent all polynomials. Polly spelled it wrong in X, Y, Z on the unit sphere of degree less than or equal to the maximum. That's total degree. When you add up the powers of X and Y and Z in the formula, you think they're more because, you know, you can have both the coefficients X and Y and powers of X, Y and Z. But basically with the relation, say, Z squared equals X squared plus one minus X squared plus Y squared. So you can use this to eliminate all powers of Z higher than one. And once you do that, then if you start counting, they're actually N squared polynomials like that. And so, you know, they're they're simple polynomials in X, Y and Z. But this is the way they're defined, because an arbitrary, you know, the power basis of polynomials one X, X squared. They're not orthogonal, but by reorganizing them, you know, maybe doing Gram-Schmidt or something, you can turn it into orthonormal set of polynomials. And this is one way of doing it. So the best way to think of them is this. They're analogous to Fourier series. And these are just particular polynomials in cosine theta, which is going to be your Z, which again, they're orthonormal. So I think if you end up doing Gram-Schmidt on the polynomials in Z, assuming that, you know, the X and Y variation goes away, you'll get exactly these polynomials for the zero case. And then these are just generalizations of them. So if you look on the Web under spherical harmonics, you'll probably find some Web sites with pictures with a red or green, red when they're positive and green where they're negative. And they just sort of blobs on the sphere, low frequency blobs. OK, so if we do that, I just I just erased my integral, unfortunately. But the integral we wanted to compute, integral factor one times P, environment lighting of C, the omega C over the hemisphere or it could be over the whole sphere. Right. Because that cosine term in factor one was really supposed to be zero. You didn't take the dot product of the normal with the direction C and continue to let it be negative values. It's clamped to zero. And if you clamp it to zero, you can think about integrating over the whole sphere because these basis functions are over the whole sphere, not the hemisphere. But that's OK. Factor one will be zero if your C points with a negative dot product with the normal N. So it's integral over the whole sphere now. That's the integral over the sphere of factor one, which was going to be the sum of T. I. Y. I. C. Times the sum. Say I equals. They will approximate it to end term. Well, to the end squared. Right. If we want to make a single index instead of two indices here, then somewhere I wrote it in the notes. I hope what what the formula is that I equals. I just thought L times L plus one plus M plus one. And this is just a way of turning. You know, you count zero and then a one with a minus one, a one with one and a one with two and then two for L and go from minus two minus one zero one two. It's just a way of ordering that into a single sequence so we can have only one subscript instead of two. And so we'll do the same thing. J equals one to N squared. What I call for the light in the environment I called the coefficients of the basis functions. L.I. I better use J, though, because. Y. J. Upside deep side, which is. Let's take the sums out of the integral and the constants out of the integral. T.I. Maybe I have to make a double sum first. J equals one to N. T.I. L.J. And then the integral of Y. I. Upside Y. J. Upside D. I guess in the book they call this the omega psi. OK, but we just saw those were delta I.J. Right. Because they're orthogonal. That's what it means to have orthogonal basis functions. That's on this. So this integral is one only if I equals J and zero otherwise. So, you know, that's a property you can prove about Fourier series. You can also prove about these specific basis functions for spherical harmonics. And so only if I and J are equal, does this integral equal to anything? In that case, it's one. So what this really is, is the sum I equals one to. I should have N squared here. T.I. And then the J has to be equal to I. So it's T.I. L.I. Right. So this is an approximation to the extent that this many terms are a good representation. Actually, suppose this many terms are a good representation of the slowly varying environmental illumination. Some say it's band limited, so it has only these frequencies in it. Then I could make these terms go up to infinity. All of the terms have an exact representation of the visibility, which has like hard edged shadows in terms of how the visibility of the environment depends on like my eye looking past my nose or something. You know, it's going to be a hard edged angle when the viewing ray hits the tip of my nose. It's going to be then and any sudden changes, discontinuities represent high frequencies. But nevertheless, if the environmental illumination has low frequencies, I can still take this sum up to infinity. But the only terms I'm going to get are when they're non-zero terms here. Right. And so, you know, basically the way I wrote it there, I would have to put J's there. It's still only N squared terms. So if any one of these two functions is band limited, this gives you the right answer, even if the others aren't. OK, so let's see, where am I now? How can we estimate these? Because we want to precompute these. If we can precompute both these coefficients and those coefficients, then in the renderer, it's just something we can put in a vertex program. Say we precompute this one. Did I have a definition of it somewhere? This one, I erased my factor, but it's got the visibility in it and the cosine and the BRDF. That depends only on the position on the surface and the geometry of the local object you're trying to draw in the environment. It doesn't depend on the environmental illumination. So you could precompute that. And in the initial paper, they precomputed it per vertex. But they also suggest that you could precompute it per texel on a parameterized surface and store it in a texture map instead of per vertex. And so that precomputation is done. And this one doesn't depend on the geometry of the object. You put multiple objects in the same environment if you expand the environment in terms of these. And in both cases, you're trying to compute this representation. And because you have an orthonormal basis function, you can say that if you really had this expansion, then if we take, let's take this one. If we take the integral of L environment psi times y j of psi d psi. That's the integral of this series. L i y i of psi y j of psi d omega psi. And this this L I, we can take outside of the integral. And I forgot my summation over the sphere. And this is only true if J equals I, in which case it's one. So this is just L I, which means we can find these coefficients by integrating the function. I'm sorry, it's just the L J. Right. Because J L I was in this in the sum, even even infinite sum. And it's the J which is the particular one that we're trying to find. And so we can find this one by integrating the function we're trying to represent in the orthogonal series by each basis function to find the coefficient. It probably looks familiar from Fourier series, too. So to do that integral, we can just sample rays. So let's choose choose a large collection of, say, random directions on the unit sphere and do Monte Carlo integration. And estimate the integral of the unit sphere of this particular function L environment. Psi times Y J of psi as this average, say, large collection of, say, capital M one over M times the sum K equals one to M of L environment. Let's call our random direction psi I and the probability of and multiplied by Y. I think I was calling him K for my samples and I'm trying to estimate. Was that a J there? I can't even see why J of psi K and divided by the probability, which is just four. Let's see, the probability is one over four pi, right? Because the total unit sphere has four pi on it. So we have to divide by the probability one over four pi. Right. And that's our Monte Carlo integration. If you didn't want to choose them randomly because you thought they might be clumped up, which might give you bad statistics, you might start out, say, or choose one of my directions. So I think I called him Psyche using four times subdivision of an octahedron or an icosahedron. Why don't I say icosahedron? So an icosahedron is a particular collection of five triangles at a vertex that, you know, I haven't drawn it very well, but I've got six on this or seven on this vertex. One. Yeah. Something like that. Right. It's this this collection of 20 triangles that cover that make a sphere. It's one of the platonic solids. You could start by taking 20 things or the center of each of these triangles. And then you could say, well, let's divide this triangle by its midpoint into four and then each of these by midpoints into four, as many times as you want to get smaller and smaller triangles and then choose the center of all of those. That's sort of a quasi equally spaced thing, but it's not really equally spaced because if you think about it, if you take this center triangle here and in terms of directions, you have to project it out to the unit sphere. And when you project it out to the unit sphere, these angles are going to be one fifth of a sphere. Right. Seventy two degrees. But these angles, the six of them have to come together, which means these ones are smaller. That means this length is larger than those two lengths when thought of as an arc on a sphere. So what you really have to do is take the sum of L environment of Psy K Y J of Psy K times the solid angle. Right. Where you're computing the solid angles of spherical triangle. It turns out the solid angle of the spherical triangle. Does anybody know the formula? You have a planar triangle. The sum of the angles adds up to one hundred and eighty degrees. But when you bulge it out to make a spherical triangle, the sum of the angles adds up to more than. Right. Because this would be 60 plus 60 plus 72 here. And so you take the excess over one hundred and eighty degrees for the sum of the angles. And that's a formula. Some of the angles between the plane from this arc and the plane from that arc of those two great circle arcs. Those are the interior angles of a spherical triangle. And you take the excess over pi and that represents the area. So which is by definition solid angle. So you could do it that way, too. Instead of a random one, you could do a deterministic one. OK, so that's for the environment. And similarly for this factor here. So this was my LJ, I guess, my coefficient that I had at the far end of that board. And similarly. DJ is. The integral is, if I remember, was in my factors, my bidirectional reflection function at X. If it's diffuse, this is really just row diffuse over pi is the bidirectional reflection function. So let's do that case first. And then we have cosine between an X and psi. And then we have the visibility in direction psi of this infinite sphere. In other words, this visibility is zero if the ray hits an object and one if it gets out to the environment. The omega psi. And we'll estimate that as a similar sum. One over M of. Well, I guess I can take my row D over pi outside it of cosine and X. Psy K times the visibility from that point on the surface to the environment at infinity in that direction. Let's do this this sum. But you could, you know, do that some too. And so that just means tracing array from the point X, either a vertex or a texel position on the parameterized surface out to infinity to see if you hit another part of the same object. And if you don't, then you include this cosine term into this estimate. And so that's how you precompute from the object geometry, these terms. And then you can put an object in the environment. This you could compute, say, by the six faces of a cube. Right. And then if you're doing it that way, you know, you might want to take directions to the pixels in the cube. And this represents the solid angle of each cube. And then this times that solid angle would represent that hemi cube projection factor we talked about. If you put the cosine in that we talk, I gave you formulas for when we were talking about the hemi cube method. And so you could use the pixels themselves to compute the visibility by, you know, putting a cube at each vertex or each texel instead of two doing random rays, you could use graphics hardware. And so these ones depend on the geometry and these ones depend on the environment. So what that means is you pick an average position. Actually, the original paper said you could do these at multiple points and interpolate if you really wanted the environment not to be at infinity. But say the walls of this room are finite or you say you had finite light sources close to the object. But the usual way is just to assume it's an environment infinity independent of position. So you can put these precomputed objects anywhere you want. You could put multiples in them of them. They just wouldn't affect each other in the sense that it's only the distant environment that's giving you this diffuse shading. And it turns out you can also rotate the object. So what you do there is instead of rotating the object, which would mean rotating, you know, redoing this visibility can do thing for each point on the surface in question. You think about rotating the environment the other way around. And it's sort of like if you if you rotate a circle, you're basically changing the phase of these four coefficients, which are maybe still on the board here. Right. It's just like adding a constant angle to theta, which is just going to change. Well, in this form, it's like adding a constant in this term. It's getting the cosine could change to the sine and vice versa if you turn by 90 degrees. But that's easy to do. It just involves relative changes between these two terms for the same end. And it turns out for the spherical harmonics whose equations are erased. What it does is for the same L, it mixes around the. Let me see. What is which way is it? See if I wrote on my notes. I think it just permutes the L. Right. The L. The L terms. I mean, for a constant L, it's just shuffling those terms around. And so it's the transformation is a bunch of two L plus one by two L plus one matrices, as I recall. In any case, there are analytic formulas to tell you for any rotation angle. They're a bit complicated, but you could precompute them. There are many fewer of them for a small number of terms than there are actual directions on your Hemi cube. If you're going to have high resolution environmental models so you can analytically find a formula to reorganize the basis elements. Basically, it's just a formula which expresses each LJ in the new coordinate system with a big matrix. Multiply by from the. Coefficients in the other coordinate system. And so you can rotate the environment, which is equivalent to rotating your object. And if your environment is really at infinity, you can translate the object also. Or you could replace it by a completely different environment map. And so it gives you a way in real time to compute the illumination from a distant environment just by that product. And you can do that per vertex. Now, the original paper also included glossy reflections where the BRDF did I. Where's my equation here? You could actually put X, Omega, Omega. Well, this is in our notation. This would be Theta for the viewing direction and Psi. And then this thing would depend on Theta. And it turns out that you can use a matrix. Which transforms the distant environment into the environment that's viewed from the object. Now, the one thing I didn't have time to talk about is the further iterations. You can actually do multiple bounces of the object of the light on the object and get the diffuse inner reflection. And so I'll talk about that next time. It's in the book. Hopefully it'll only take me 10 minutes next time, because I think there was something else I wanted to cover next time. At this point, I don't remember what it was. I want to cover hierarchical radiosity next time. So I'll try and talk about including multiple reflections rather than just the occlusion in the environment into this form, into this pre-computed radiance transport. And then I'll talk about the hierarchical radiosity. And remember, everybody, Projects Friday, I think I put a notice on the website, but I want you to bring either laptops or USB drives to be able to show me. Well, we'll use that projector and we'll pull down the screen and we'll all share each other's images. OK, that's it for today.