We have not considered any radiance contribution coming from the indirect sources to the surface of interactions. The model we have used so far couldn't simulate the global illumination and we had to fake it by adding an ambient term. Even if the ambient term can be adjusted to a reasonable level, we cannot simulate phenomena, such as color bleeding, caustics, etc., that can only be simulated correctly by global illumination algorithms. Therefore, it is a must to implement a global illumination algorithm to produce more realistic images.
In this assignment, we are required to solve
the rendering equation by using monte carlo techniques. Producing an unbiased image with a low variance is not a simple task and has been a research topic since the early days of this field. Thankfully, there are great resources offering several variance reduction techniques to reduce the noise and for faster convergence. Here are some of the resources I have been using:
After reading and learning a lot from these resources, I started to implement my path tracer. Although my path tracer is very simple and using a few of the techniques that I've learned, I am planning to implement them in time. For now, my path tracer has the following features:
 Progressive path tracing: It is great to see your images converge in realtime and to move around the scene you created. For most of the scenes, I take interactive frame rates.
 Cosine weighted importance sampling: It is the only variance reduction strategy I used for the indirect lighting calculation. However, in future, an importance sampling for the bsdf term seems inevitable.
 Explicit light sampling: At every bounce, I take samples from each light source. If many light sources are involved, we may choose one of them uniformly or by a pdf that we defined according to some criteria, for example, power of lights.
 Smooth dielectrics: Perfectly reflective and refractive surfaces are easy to implement if you already implemented them in your ray tracer. Since they have delta distributions, only a single path is meaningful to follow. Therefore, everything is the same as in the case of our old ray tracer. However, we do not send both reflected and refracted rays for our refractive surfaces since the recursion would be an overhead in general and it is already deadly for the performance of the CUDA programs. Instead, we choose one of them with probability determined by the fresnel equations, which I talked about in the previous posts.
Also, I am planning to
 read this paper to implement both refraction and reflection for the rough surfaces.
 add multiple importance sampling for direct lighting to reduce the noise related to glossy reflections.
 add importance sampling for the bsdf term of the indirect lighting.
Finally, let me sum the basic things up for a path tracer in some brief paragraphs. We separate incoming radiance sources into two parts as direct and indirect. We sample lights explicitly for the direct part and while sampling this part, it is reasonable to importance sample the lights as well as BSDFs. For the indirect part, however, we have no idea about which direction may potentially give us higher radiance values since it is the problem itself that we want to solve. However, we can say which direction can contribute more according to the cosine term and BSDF.
Furthermore, since we explicitly sample the lights, if we hit any area light source with the indirect rays we have sent, we cannot take any contributions from the light sources we already sampled in the previous bounce. That is, if we hit a light that we already sampled, we consider as if it does not contribute. Also, beware that we do not sample the lights explicitly for perfectly reflective and refractive surfaces because that would be meaningless due to their delta distributions. Therefore, if the rays that have been sent after the interaction with smooth surfaces hit a light source, we take the contribution coming from these rays.
As before, I used jittered samples to send the rays not only to the center of the pixels but to the entire square owned by the pixel.
I tried several russian roulette strategies and fixed depth values for ray terminations. While russian roulette can make a substantial decrease in fps, fixed depth values result in better performance. However, unlike fixed depth values, russian roulette generates unbiased images, which is more important for me now. The outputs shown below are generated using russian roulette that considers the radiance contribution weight of the ray.
Let us see the outputs.

Cornellbox rendered with 3800 samples 

Sponza considering only direct illumination 

Sponza rendered with 650 samples. 

Color Bleeding Effect 

Caustic Effect 

Nothing Special 