Rendering equation - brdf / pdf > 1?

Question

I am creating a path tracer and I got some strange results when calculating the strength of light.

I am using a rendering equation inspired by this diagram:

1. When I use the cosine everything gets really dark and the edges of my spheres are almost black. That makes sense because I am almost multiplying by cos(PI/2) = 0. I have read that it suppose to cancel out with something so I have removed it which help. But what does it cancel out with?

2. I am using Lambertian diffuse BRDF which means that the "Fr" part of the equation should be "albedo / PI". I am also doing uniform hemisphere sampling which gave me a PDF of 1 / PI. Now when I calculate BRDF / PDF I get "2 * albedo". That does not seem right. This means that the more and more bounces I make the bigger influence the light has which is nonsense. The light should be dimmer from the 3rd bounce than from the 1st bounce. Where did I make a mistake?

EDIT 1: I am including pseudo-code. The cosine is the dot product. And after boun There are some things I have not included because they do not contribute to the equation.

<!-- language: lang-cpp -->
//Initialization
Ray ray;
float3 passthrought = (1,1,1);
float3 color = (0,0,0);

//Bouncing
for(int I = 0; i < 4; i++){
    Object obj = ray.intersect(...);

    color += obj.emission * passthrought;

    float3 newDirection = sample_hemisphere();
    float pdf = 1 / (2 * PI);
    float3 brdf = obj.color / PI; 

    if(pdf > EPSILON)
        passthrought *= dot(normal, newDirection) * brdf / pdf;
    else
        passthrought *= 0;

    //Prepare new bounce
    ray.origin = ray.direction * ray.distance + ray.origin;
    ray.direction = newDirection;
}


float3 sample_hemisphere(){
    //Tranformation to world space
    float3 w = ray->worldNormal;
    float3 axis = fabs(w.x) > 0.1f ? (float3)(0.0f, 1.0f, 0.0f) : (float3)(1.0f, 0.0f, 0.0f);
    float3 u = normalize(cross(axis, w));
    float3 v = cross(w, u);

    //Sampling from 2 random variables <0, 1>
    float e1 = rand();
    float e2 = rand();

    float s = sqrt(1.0 - e1 * e1);
    float phi = 2 * M_PI_F * e2;

    return normalize(u * cos(2 * M_PI_F * e2) * s + v * sin(2 * M_PI_F * e2) * s + w * e1);
}

Edit 2: I just wanted to include an example render with both variants and somehow it works now. The black border around objects is gone. I think I had to calculate the cosine with normal and the ray towards the camera which was wrong. Anyway, I am including 2 renders one with a cosine of normal and outgoing ray and one without it. Keep in mind that Reinhard tone mapping with white point set to max luminance and gamma correction of 1/2.2 was applied. Which one of them is more right?

With cosine:

without cosine (is smoother because it has like 10 times more samples):

Answer 1

When I use the cosine everything gets really dark and the edges of my spheres are almost black. That makes sense because I am almost multiplying by cos(PI/2) = 0. I have read that it suppose to cancel out with something so I have removed it which help. But what does it cancel out with?

That cosine has a geometrical meaning, keep it, or if you take it off, be sure you have a mathematical reason to do so. Grazing angles distribute much less of their radiance to the surface, and at an angle of 90°, the received radiance is zero.

I am using Lambertian diffuse BRDF which means that the "Fr" part of the equation should be "albedo / PI". I am also doing uniform hemisphere sampling which gave me a PDF of 1 / PI. Now when I calculate BRDF / PDF I get "2 * albedo". That does not seem right. This means that the more and more bounces I make the bigger influence the light has which is nonsense. The light should be dimmer from the 3rd bounce than from the 1st bounce. Where did I make a mistake?

The probability density function $p$, by its definition, must have an integral of 1 on the whole definition set, thus on the hemisphere. If you do the math, you'll get:

$\int_\Omega p \, \text{d}\Omega = \int_0^{2\pi} \int_0^{\pi/2} p \, \sin \theta \, \text{d}\theta \, \text{d}\phi = 2\pi \Rightarrow p = \frac{1}{2\pi}$

Sampling your PDF will give you $\theta = \arccos \xi_1 $, $\phi = 2\pi \xi_2$, where $\xi_1, \xi_2$ are independant random real numbers between 0 and 1. Sampling a PDF is not always easy and is another topic in itself, but for more details you can read this excerpt from the excellent PBR.

But as you saw it, you can do better : the cosine term in the rendering equation that bothers you from the beginning tells you that you should favor the directions that maximize this term. So let's choose a pdf such that

$p = \alpha \cos \theta_i$

If you want to know the value of $\alpha$, you have to remember that the integral of the pdf must be equal to one. If you do the math, this will give you $p = \frac{\cos \theta_i}{\pi}$.

Sampling this new PDF with the method shown in PBR will give you $\theta = \arccos \sqrt{\xi_1}$, $\phi = 2\pi \xi_2$.

Now, and only now, you can notice that the cosine of the rendering equation simplifies with the one in the PDF. You can then remove them.