 What a wonderful day we have today and what a wonderful time it is to write a path tracer. So why don't we get started? What we are going to be looking at is a program called SmallPaint which is a small path tracer in effectively 250 lines of code and it contains everything that we have learned during this course. We're going to be able to compute soft shadows, anti-aliasing, Monte Carlo integration and even quasi Monte Carlo sampling which basically means low discrepancy sampling. This version of the program is going to be able to compute refraction, color breathing and caustics and in the end the binary you compile from the code can be compressed into less than 4 kilobytes. This is how the end result looks like it has a beautiful painterly look which actually comes from a bug and you can also see that the light source up there, the whitish looking sphere is you could say perfectly anti-aliased. In order to achieve this with a recursive ray tracer and no global illumination algorithm you would need to compute the very same image on a much larger resolution and then scale it down to a smaller image like this. This anti-aliasing effect you get for free if you compute Monte Carlo path tracing. The question is how is this exactly done and now is the best time to put everything into use that we have learned so far. Let's get started. We have a completely standard vector class. It is a three-dimensional vector. It has its own constructor, the classical operators that you would expect and we also have a dot product for the vector and a cross product for the vector. It is also obviously possible to compute the length of this vector so nothing too exciting or important here but we definitely need to build on a solid vector class. Now the representation of a ray, a ray has an origin and the direction and if you take a close look at the constructor you can see that when you initialize a ray with a direction then this direction is normed and the reason is that when we compute the dot products between these vectors most of these information needs to be directional information so we are not interested in the magnitude of the vector only interested in the direction of the vector. A good way to get rid of problems where you would initialize your ray with directions that are not vectors of unit length but you can do is to normalize this input in the constructor so you will never have headaches about incorrect results where you have no idea what is really happening. What is the representation of an object? Well an object has a color. This is denoted by CL. This is actually very sloppy notation and this you can imagine as Elbedo but it is not a double so it's not a number it is actually a vector and the reason for this is the fact that we need to define the Elbedo's how much light in different wavelengths is being reflected and how much is being absorbed by the given object. Now object may also have a mission if they have some non-zero emission then they are light sources and by type we have an integer that would specify what kind of BRDF we have. It is also important to have an intersection routine and some other function that can compute the normal of the object. Now these are of course virtual functions we don't define them or an abstract object but they would have to be implemented in other classes that inherit from the object. Let's take a look at the sphere. So sphere has this C and R C is the center of the objects and R is the radius. The constructor is trivial. We have the intersection function. Now if you hopefully remember all three of the elements of the quadratic function that we need to solve but if you take a good look then you will see that A is missing from here and the question is why is that? The answer is we are multiplying D with D the direction vector of array with itself and if it's a vector that is normed so it is of length 1 then this A will always be 1. After that we can deal with the discriminant. The discriminant is normally not B squared minus 4 AC. Remember A is 1 here so it's omitted but it is behind the square root and this square root is completely omitted here which looks like a mistake but it is not. It is essentially an optimization step because you will see that we are testing the discriminant if it's less than 0. If it's less than 0 then we don't need to deal with this equation because the solutions for the quadratic equation are going to exist in the plane of complex numbers and that's not a valid T. It's not a valid distance where we would intersect the sphere. If this is not happening then we can compute the square root for the discriminant. Why after that? Because if the discriminant is bigger than 0 then the square root is not going to make a difference. So what we can essentially do is to postpone the square root calculation after the test. Note that square roots are really really expensive so we can save up lots of computational time if you omit this calculation. There is a really nice square root hack in the source code of quick 3 which is by the way open source. Take a look at how people are trying to hack together functions in order to work better and faster than they should because square roots are super expensive and there are some really interesting hacks in order to speed them up. We have the plus and minus term and the division by 2a is again postponed. And that's also another interesting question. Why is this postponed? So you can see that the sol 2 is divided by 2 and the sol 1 is also divided by 2 but only after the test. So it is possible that if we have the solution 2 if it is bigger than epsilon then we have the first expression after the question mark but if it's not then we are looking for the second expression after it which is another test and if the answer is no for that as well then we are going to return 0. This would mean that we don't have any hits or the hits are behind us and we are not interested in intersections that are behind our array. There is a possibility that we encounter this and in this case I don't want to waste my time by dividing these solutions by 2 because I'm not going to use them. Why am I splitting hairs here? That's an important question. Why do we need to optimize so much? Because if you grab a profiler a program that is able to show you in what functions are you spending most of your time in then this profiler would show you that more than 90% of the execution time is spent in the intersection routines. So you have to have a really well optimized intersection routine. Some programs have replaced these expressions with assembly in order to speed it up even more. So how do we compute the normal of a sphere? Well very simple. What we have here is p minus c. Now what does it mean? So if I have a minus p then this means a vector that points from b to a. So what this would mean look at the figure here. If I would have a circle then this would mean that from the center it is pointing towards the given points of the sphere. Now this is also divided by r because you could imagine you have a sphere that is not of unit radius and if it's not of unit radius then this normal vector would have a length. You could compute a normalization we have a normalized function in our vector implementation but it also has a square root so it would be much better to have something that's faster. Well if you just divide by the radius of the sphere then you would immediately get a vector of unit length. So in the end we can get the correct result by just simply dividing by r. Excellent. Now we have a perspective camera here. Hopefully you remember from the first lecture we are basically just copy-pasting this expression here we have derived them rigorously and analyzed how this exactly works. A simple intuition basically what we are doing we have as an input an x and a y. Basically this means give me a pixel with this displacement and what it would give you back is the world's base coordinates of these pixels. Uniform sampling of a hemisphere what is this for if we are encountering a diffuse object what we would like to do is to send a ray out on the hemisphere of this object. This we would need to uniform example this means that the diffuse PRDF is one over pi or row over pi if you take into consideration the elbows and you need transforms in order to do it there is a reading behind this link how it works exactly is drawing uniform samples on a plane which is simple and then we are projecting it on the hemisphere that's basically all there is. What about the trace function? As you can see here in the first line this code says that there is a maximum depth. Now clamping after a maximum depth value is not really optimal because whatever number you put in there the higher order bounces are going to be completely omitted. Now the real solution would be Russian-Rulet past termination which we fortunately also have after depth of an arbitrary number like five you start the Russian-Rulet routine which basically says there is a probability for stopping the light path right there and we generate a random number and compare to this probability. If we don't hit this probability then we will continue our path but we will multiply the output and the contribution of this light path by this given number that we have specified in one of the previous lectures. So this was implemented by Christian Mahacek and kind thanks. And you can see that later we are going to use this RR factor in order to multiply the contribution of a ray later. Now what about the intersection routine? This is definitely not the best way to do it but it's sure as hell the easiest way to do it. We specify this T which is going to be the intersection distance. How far we are from the start of the ray and how far is this intersection exactly? ID is basically which object did we hit and then we iterate through all of the objects in the scene and what we are interested in is calling the intersection routine. This will return a T how far is the intersection and what I am interested in an intersection that is the smallest number. This means the closest intersection and also something that is larger than epsilon because if I would tolerate zero then this would mean that self intersections are accepted. Therefore every single object that I am on is going to be the first intersection. I'm not interested in this. I know where I am. I just want to know where I am continuing my path. If there is no intersection then we return. There is zero contribution. Where is the intersection in world space? We denote this by HP. This means hit point and where we have started a ray we traveled along the direction of the ray with this T amount where the intersection is. So this is where we ended up and the origin of the new ray is going to be this hit point. What is the normal going to be? Well we just grabbed the object that we intersected and we are taking the normal with the given function. What is the return radiance? We simply add the emission term. This is the emission term on all three wavelengths. There is a magic multiplier disregard that and then we continue evaluating the second part of the rendering equation. The first part is the emission and the second is the reflected amount of light. Let's continue with the inside of the trace function. If we hit an object with a type of one then this is a diffuse BRDF. The next functions compute the next random number for the low discrepancy halton sampler and the direction is going to be a random sample in the hemisphere, a completely uniform random sample in the hemisphere of this object. What we have here is this N plus the hemisphere function. This is intuition. This is not exactly what is happening. I have just shortened the code slightly in order to simplify what is going on here. The code that you will download will have the real deal in there. Now then we compute the cosine term very straightforward and on the TMP we are going to instantiate a new vector and this is going to hold the recursion. So the subsequent samples that we shoot out from the hemisphere are going to be added and added to this TMP. Now is the time for recursion. We pass the ray and the scene to the trace function. The ray is actually not the current one, it's the new one. So basically we set up the new hit point and the new direction of the ray and this is what we are going to pass to the trace function. We increment the depth variable because we have computed the bounds. The TMP is going to be a variable where we gather all these radians and we put every other parameter that is needed to compute one more bounds. Now the color is going to contain the cosine term and all these radians that is collected from the recursion and we multiply it with the CL.xyz which is basically the BRDF. So this is the right side of the rendering equation for a diffuse BRDF. This is multiplied by 0.1. This is just a magic constant. Now what about a specular BRDF? What if we hit a mirror? Well very simple. We compute the perfect reflection direction. You can see the ray dot D and we again set up this variable to collect the radians in there and we are not doing anything. We are just going to add the radians as we get reflected off of this mirror. Then we are going to compute subsequent bounces and this is going to be stored on this TMP. So this is what we are going to add to this radians. What about a refractive material? Well we have every bit of knowledge that we need for this because essentially this is the vector version of Snals Law. What does it mean? Well the original Snals Law that we have computed is in 1D. So it only gives you one angle. But if you are in 3D you are interested in angles in two different dimensions. This is nothing but the extension of the very same law into higher dimension. Now where is this implemented exactly? You can see the cosine of theta 2. Note that n1 and n2 is considered differently because one of these media is always going to be air. Therefore one of the indices of refraction is always going to be 1. The rest is just copy paste. And again you can see that the square root is missing and we are going to postpone this after the test of cosine t2 because if it is actually not larger than 0 then we are not going to need this variable at all. Therefore we can postpone this after the test again. What about the direction of the outgoing gray? Well this is just copy paste from the formula that we have derived before. So simple as that. Obviously we again need the recursion because if we go inside a glass sphere then we are going to compute the refraction. So we are going to be inside of the sphere. What does it mean? One that we have to invert the normal because we are inside so the normals are flipped. And again we need to compute the trace function which is the recursion. So we are also interested in the higher order bounces. Onwards to Fresnel's law. What is the probability of reflection and refraction when rays are bouncing off in different directions in different angles of refractive surfaces? Implemented by Christian Hathner. So a big thanks for him. It is very simple. You can see that it is exactly the same as what we have learned in mathematics. So this is the R0 term. This is the probability of reflection in normal incidents. And we are interested in the square of that. And note that you don't see the n1 and then 2. This is because one of them is always going to be air or vacuum. So it is going to have the index of refraction of one. Now what about the final probability of reflection? It is also coming from the formula. We have every bit of information we need. So we just put in there this term with the cosine attenuation. How does the main function look like? Well we have some wizardry with C++11 lambda functions. But basically this is just a shortcut in order to be able to add a new sphere or a new plane to the scene in one line of code. Spheres are given by their radius position color by color. We obviously mean albedos emission and type. Type means what kind of BRDF we have? A diffuse a specular or a refractive BRDF. Now for planes we have position normal color emission and obviously type. So what kind of material we have? So using just one line of code you can add a new object and specify everything every information that you would need to it. Now we also add the light source and we specify the index of refraction for the refractive BRDFs. And we also specify how many samples per pixel would we like to compute? Onwards to the main loop we have two for loops that iterate through the width and the height of the image plane. Now vector C is color. It's again very sloppy. What it means is actually the radiance that we compute. We instantiate array. What is going to be the origin of the ray? This is going to be a 0 0 0. So this is where the camera is placed. What is going to be the direction of the ray? Well we connect this ray to the camera plane. And we specify which pixel we are computing with i and j and then we add this weird random number to it. Now what this means is actually filtering. In recursive ray tracing what you would do is you would only send the ray through the midpoint of a pixel and that's it. You would compute one sample per pixel. In Monte Carlo part tracing you're computing many samples per pixel and they don't have to go through the midpoint of the pixel. You would sample the area of the pixel. And this gives you anti aliasing effects for free if you use it correctly. What is going to be the direction of the ray? Well this is again the same a minus b. The b is the origin of the ray and a is the camera coordinate. So what does it mean? That it is pointing from 0 to the camera plane. And we normalize this expression to have a ray of unit length. Now we obviously call the trace function. The number of bounces is 0 and we pass every information that needs to be known in order to compute these bounces. So we provide this initial ray and the scene and everything else. Obviously we also pass the C and this is going to collect all the radiance there is in the subsequent bounces. And then after this recursion is done we deposit all this energy, all these radiance to the individual pixels. And then we divide by the number of samples because if we wouldn't do this then you remember the one over n multiplier everywhere in Monte Carlo integration. If you wouldn't do this then the more samples you compute the brighter image you would get. And this is obviously not what we're looking for. At the very end we create a file. This is the ppm file format where you can easily write all your contributions in there. We also start a stopwatch in order to measure how long we have been tracing all these rays. So very simple, very trivial and when we are done we close the file. It has the image in there and we done write how long the rendering algorithm has been running for. And basically that's it. That's it. This is effectively 250 lines of code that can compute indirect illumination, caustics and every global illumination effect. And it can compute images like this. This is one student submission from previous years. Absolutely gorgeous. This is the fixed version of small paint where there are no errors in the sampling. Another one from Michal Kama. This actually looks like I don't know if you are into the music band Boards of Canada but this looks exactly like one of their album covers. Love it. Really cool. And also Syrpinski Triangles from Christiane Kusla. You can find the link for the code in there and take a crack at it. Just try it, build different scenes, try to understand what is going on in there, try to mess the code up. I wonder what happens if I would not normalize this vector. Play with it. It's a really small, concise and really understandable path tracer. So take your time and play with it. It's lots of fun and you can create lots of beautiful, beautiful images with global illumination. Thank you.