# Adding cos(θ) from rendering equation ruins the renderings

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

• dot(intersection.normal, scattered_ray).normalize() - dot should return a scalar and the scalar should not be normalized. May 19, 2022 at 4:16
• What is this: if attenuation.r <= 0.0 && attenuation.g <= 0.0 && attenuation.b <= 0.0 {return attenuation;}? Also I don't get why you're getting different results, the two estimators ought to converge to the same thing. What are the albedos of your surfaces? Are they in [0,1/PI]^3? If they aren't, then the scene is not energy conserving, which may explain the second case, but it doesn't explain the first. May 19, 2022 at 16:27
• If you don't want to change all your albedos you can just drop the PI in both pdfs. May 19, 2022 at 17:06

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).

• @Rahn I don't think the bug is in the code that you posted, or at the very least I cannot find it there, it may be some part that you have not posted. What does your dot function look like? You also should implement what I suggested regarding the without cosine part. May 19, 2022 at 11:19
• update again. Code added (it's in Rust, but not so different from C++).
– Rahn
May 19, 2022 at 12:58

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.