Adding cos(θ) from rendering equation ruins the renderings

Question

I'm following Ray Tracing The Rest of Your Life to implement a ray tracer, but the explains on math (mainly pdf part) got me confused so I followed Rendering Equation to understand the math behind. As a result, I sample the scene with plain Monte Carlo method (derived from the Equation) and got some seemly correct result after removing the cos(θ) in Rendering Equation, but got an incorrectly dark image with cos(θ). And worse, the rendering doesn't converge: as number of samples increase, it goes darker and darker.

from left to right: 500 samples without cos(θ), 500 samples with cos(θ), 100 samples with cos(θ):

The code is not complicated. For every sampling, the camera send a ray into function ray_trace(Ray ray, int depth), the function terminates when the ray hit a light or couldn't find one after bouncing enough time (depth <= 0). When it hits something in the scene (the only material in the scene other than light is Lambertian material), a scattered ray is generated from material's scatter(Ray ray, Intersection intersection, Ray &scattered_ray) function. In Lambertian's scatter() function, the scattered ray is generated by uniformly sampling around hemisphere. Then this scattered ray becomes the new ray in function ray_trace, thus depth = depth-1, and start another round of ray tracing.

{
    // initially, this ray is generated from a camera
    if (depth <= 0)
    {
        // return black when no light found
        return Vec3(0.0, 0.0, 0.0);
    }

    auto intersection = scene.intersect(ray);

    if (!intersection.intersected)
    {
        // return black when hits nothing
        return Vec3(0.0, 0.0, 0.0);
    }

    if (intersection.hit_a_light)
    {
        // found a light
        return intersection.material.emit(...);
    }

    Ray scattered_ray;
    Vec3 attenuation = intersection.material.scatter(ray, intersection, &scattered_ray);

    if (with_cosine_theta)
    {
        Vec3 cosine_theta = dot(intersection.normal.normalize(), scattered_ray.direction.normalize());
        return cosine_theta * attenuation * ray_trace(scattered_ray, depth - 1);
    }
    else
    {
        // without cosine(theta)
        return attenuation * ray_trace(scattered_ray, depth - 1);
    }
}

Vec3 Lambertian::scatter(Ray ray, Intersection intersection, Ray &scattered_ray)
{
    double phi = random(0.0, 2.0 * PI);
    double sin_phi = sine(phi);
    double cos_phi = cosine(phi);

    double cos_theta = random(-1.0, 1.0);
    double sin_theta = sqrt(1 - cos_theta * cos_theta);

    Vec3 scattered_direction = Vec3(sin_phi * sin_theta, cos_phi * sin_theta, cos_theta);

    // generate random direction in a hemisphere
    if (dot(scattered_direction, intersection.normal) < 0)
    {
        scattered_direction = -scattered_direction;
    }

    scattered_ray.origin = intersection.hit_point;
    scattered_ray.direction = scattered_direction;

    return this.albedo;
}

Can someone point out what makes my rendering wrong?

---

update:

Following lightxbulb's idea, I tried both approaches (code at https://codeshare.io/bvR8ev):

• keep cos(θ) in ray_trace() and uniformly sample hemisphere in scatter() (WITH_COSINE=true, UNIFORM_SAMPLE=true)
• remove cos(θ) in ray_trace() and sample hemisphere with cos(θ) importance in scatter() (WITH_COSINE=false, UNIFORM_SAMPLE=false)

for the 1st approach I have (100/200/500 samples):

for the 2nd approach I have (100/200/500 samples):

left to right: 100 samples, 200 samples, 500 samples

Answer 1

Your rendering without the cosine is probably wrong since you're using uniform sphere sampling and the cosine must be the there (unless you assume your brdfs contain a term 1/cos(theta)). Also do this:

pdf = 1.0/(2.0*PI);
return cosine_theta * attenuation / pdf * ray_trace(scattered_ray, depth - 1);

and also for the non-cosine part:

pdf = 1.0/(2.0*PI);
return attenuation/pdf * ray_trace(scattered_ray, depth - 1);

If you want to use the non-cosine estimator, then you will have to use cosine distributed sampling for it. It will look like this:

Vec3 Lambertian::scatter_cosine_dist(Ray ray, Intersection intersection, Ray &scattered_ray)
{
    double phi = random(0.0, 2.0 * PI);
    double sin_phi = sine(phi);
    double cos_phi = cosine(phi);

    double cos_theta = random(-1.0, 1.0);
    double sin_theta = sqrt(1 - cos_theta * cos_theta);

    Vec3 scattered_direction = Vec3(sin_phi * sin_theta, cos_phi * sin_theta, cos_theta).normalize();

    scattered_ray.origin = intersection.hit_point;
    scattered_ray.direction = (intersection.normal+scattered_direction).normalize();

    return this.albedo;
}

Then modify the non-cosine estimator's pdf too:

pdf_wo_cos = 1.0/PI;
return attenuation/pdf_wo_cos * ray_trace(scattered_ray, depth - 1);

The pdf in this is actually cos(theta)/PI, but the cosine cancels with the cosine from the rendering equation. If you use the cosine distributed sampling for the estimator without the cosine, and use the uniform sampling for the estimator with the cosine, then you should get equivalent results up to noise (the cosine estimator ought to result in less noise).

Answer 2

What coordinate system did you use for intersection.normal and scattered_ray.direction? A typical way would be sampling the scattering direction in the local space and then transforming it into the world space. Thus the cosine theta will be computed inside the sampling function. I am not sure if it is the coordinate system that caused your issue.