# Cosine weighted hemisphere sampling is a little bit darker and arguably noisier than reference

I'm writing a small path tracer that currently:

• Samples a random light source at each bounce (direct lighting)
• Bounces rays around multiple times (indirect lighting)
• The scene only contains lambertian materials so the directions for the ray bounces are generated according either to a uniform distribution over the hemisphere or a cosine weighted distribution.

The issue: Generating the directions with the cosine weighted distribution gives a slightly darker results in some part of the image than when using a uniform distribution and the variance doesn't seem to be improved as much as it should. It also introduces fireflies (supposedly because of the PDF that can get very close to 0 when rand_2 gets close to 0, see the code below).

Results:

Cosine weighted on the left, uniform distribution over the hemisphere on the right:

The fireflies are obvious when zooming in a little bit. The front face (visible) of the small cube is darker than in the right image (uniform distribution).

The code:

Direction generation (uniform and cosine):

void branchlessONB(const Vector& n, Vector& b1, Vector& b2)
{
float sign = sycl::copysign(1.0f, n.z);
const float a = -1.0f / (sign + n.z);
const float b = n.x * n.y * a;
b1 = Vector(1.0f + sign * n.x * n.x * a, sign * b, -sign * n.x);
b2 = Vector(b, sign + n.y * n.y * a, -n.y);
}

Vector RenderKernel::rotate_vector_around_normal(const Vector& normal, const Vector& random_dir_local_space) const
{
Vector tangent, bitangent;
branchlessONB(normal, tangent, bitangent);

//Transforming from the random_direction in its local space to the space around the normal
//given in parameter (the space with the given normal as the Z up vector)
return random_dir_local_space.x * tangent + random_dir_local_space.y * bitangent + random_dir_local_space.z * normal;
}

Vector RenderKernel::uniform_direction_around_normal(const Vector& normal, float& pdf, xorshift32_generator& random_number_generator) const
{
float rand_1 = random_number_generator();
float rand_2 = random_number_generator();

float phi = 2.0f * M_PI * rand_1;
float root = sycl::sqrt(1 - rand_2 * rand_2);

pdf = 1.0f / (2.0f * M_PI);

//Generating a random direction in a local space with Z as the Up vector
Vector random_dir_local_space(sycl::cos(phi) * root, sycl::sin(phi) * root, rand_2);
return rotate_vector_around_normal(normal, random_dir_local_space);
}

Vector RenderKernel::cosine_weighted_direction_around_normal(const Vector& normal, float& pdf, xorshift32_generator& random_number_generator) const
{
float rand_1 = random_number_generator();
float rand_2 = random_number_generator();

float sqrt_rand_2 = sycl::sqrt(rand_2);
float phi = 2.0f * M_PI * rand_1;
float cos_theta = sqrt_rand_2;
float sin_theta = sycl::sqrt(sycl::max(0.0f, 1.0f - cos_theta * cos_theta));

pdf = sqrt_rand_2 / M_PI;

//Generating a random direction in a local space with Z as the Up vector
Vector random_dir_local_space = Vector(sycl::cos(phi) * sin_theta, sycl::sin(phi) * sin_theta, sqrt_rand_2);
return rotate_vector_around_normal(normal, random_dir_local_space);
}


Ray bounce part:

for (int bounce = 0; bounce < MAX_BOUNCES; bounce++)
{
if (next_ray_state == BOUNCE)
{
HitInfo closest_hit_info;
bool intersection_found = intersect_scene(ray, closest_hit_info);

if (intersection_found)
{
//Indirect lighting
int material_index = m_materials_indices_buffer[closest_hit_info.triangle_index];
SimpleMaterial mat = m_materials_buffer_access[material_index];

throughput *= mat.diffuse;
if (bounce == 0)
sample_color += mat.emission;

// Direct lighting
float pdf;
LightSourceInformation light_source_info;
Point random_light_point = sample_random_point_on_lights(random_number_generator, pdf, light_source_info);
Point shadow_ray_origin = closest_hit_info.inter_point + closest_hit_info.normal_at_inter * 1.0e-4f;

Color radiance = Color(0.0f, 0.0f, 0.0f);
{
const SimpleMaterial& emissive_triangle_material = m_materials_buffer_access[m_materials_indices_buffer[light_source_info.emissive_triangle_index]];

//Cosine angle on the illuminated surface
//Cosine angle on the light surface
//Falloff of the light intensity with the distance squared
//PDF: Probability of having chosen this point on this exact light source
//The illuminated surface is Lambertian
}

Vector random_dir = cosine_weighted_direction_around_normal(closest_hit_info.normal_at_inter, random_direction_pdf, random_number_generator);
Point new_ray_origin = closest_hit_info.inter_point + closest_hit_info.normal_at_inter * 1.0e-4f;

ray = Ray(new_ray_origin, normalize(random_dir));
next_ray_state = RayState::BOUNCE;

sample_color += radiance * throughput / random_direction_pdf;
}
else
next_ray_state = RayState::MISSED;
}
else if (next_ray_state == MISSED)
{
//Handle skysphere here
break;
}
else if (next_ray_state == TERMINATED)
break;
}


Light source sampling part:

Point RenderKernel::sample_random_point_on_lights(xorshift32_generator& random_number_generator, float& pdf, LightSourceInformation& light_info) const
{
light_info.emissive_triangle_index = random_number_generator() * m_emissive_triangle_indices_buffer.size();
light_info.emissive_triangle_index = m_emissive_triangle_indices_buffer[light_info.emissive_triangle_index];
Triangle random_emissive_triangle = m_triangle_buffer_access[light_info.emissive_triangle_index];

float rand_1 = random_number_generator();
float rand_2 = random_number_generator();

float sqrt_r1 = sycl::sqrt(rand_1);
float u = 1.0f - sqrt_r1;
float v = (1.0f - rand_2) * sqrt_r1;

Vector AB = random_emissive_triangle.m_b - random_emissive_triangle.m_a;
Vector AC = random_emissive_triangle.m_c - random_emissive_triangle.m_a;

Point random_point_on_triangle = random_emissive_triangle.m_a + AB * u + AC * v;

Vector normal = cross(AB, AC);
float length_normal = length(normal);
light_info.light_source_normal = normal / length_normal; //Normalized
float triangle_area = length_normal * 0.5f;
float nb_emissive_triangles = m_emissive_triangle_indices_buffer.size();
//m_out_stream << triangle_area << sycl::endl;
pdf = 1 / (nb_emissive_triangles * triangle_area);

return random_point_on_triangle;
}

bool RenderKernel::evaluate_shadow_ray(Ray& ray, float t_max) const
{
HitInfo hit_info;
intersect_scene(ray, hit_info);
if (hit_info.t + 1.0e-4f < t_max)
//There is something in between the light and the origin of the ray
return true;
else
return false;
}


One thing confuses me in your code. Since you say that you are sampling BRDF and light source together, I figured it might be possible that you did change the sampling methods between uniform and cosine-weight but you forgot to have the corresponding PDF for the MIS weight. Then I found that you aren't using any MIS, and instead, the evaluation for direct and indirect components seem weird:

• Direct illumination estimation:
sample_color += radiance * throughput / random_direction_pdf;


I don't think this is correct. You see, the radiance here is the direct component and has nothing to do with ray direction sampling (unless you are using MIS with BRDF PDF evaluation). That means:

Vector random_dir = cosine_weighted_direction_around_normal(..., random_direction_pdf, ...);
Point new_ray_origin = ...;

ray = Ray(new_ray_origin, normalize(random_dir));
next_ray_state = ...;
sample_color += radiance * throughput / random_direction_pdf;
/// ======== should actually just be like the following ==========
// accumulate direct illunination
// update ray direction to account for indirect subsequent bounces
Vector random_dir = cosine_weighted_direction_around_normal(..., random_direction_pdf, ...);
Point new_ray_origin = ...;

ray = Ray(new_ray_origin, normalize(random_dir));
next_ray_state = ...;


Why? See the following figure:

The red line segment is shadow ray, and the sampling decision is only related to how you:

1. sample an emitter (if there are multiple)
2. sample a point from that emitter.
3. get PDF (area product measure) for sampling the very point on that emitter
4. Transform to solid angle measure (falloff, cosine attenuation, etc. I think you did it right)

Essentially, after if (!in_shadow) branch being finished, the direct illumination already has its unbiased estimation. So you just add it to radiance and account for the path throughput from the previous bounces.

• Throughput update:
throughput *= mat.diffuse;


Note that all the subsequent bounces are sampled (ray direction sampling and intersection test to obtain new ray origin and direction), you should account for the sampling event here by:

throughput *= mat.diffuse * sycl::max(dot(closest_hit_info.normal_at_inter, random_dir), 0.0f) / random_direction_pdf;


This should be correct for throughput update, the derivation is simple: the following figure shows how the incident ray gets attenuated by BRDF and cosine term.

Since our light incident direction (corresponds to random_dir) is sampled, we have the following estimator and we can easily verify its expectation:

$$\hat{I}=L_i(x, w_i)\frac{f\cdot \cos\theta}{\text{PDF}(w_i)}\\ \mathbb{E}(\hat{I}) = \int_{\Omega_{\text{hemisphere}}}\frac{L_i(x, w_i)f\cdot \cos\theta}{\text{PDF}(w_i)}\text{PDF}(w_i) dw_i = \int_{\Omega_{\text{hemisphere}}}L_i(x, w_i)f\cdot \cos\theta dw_i$$ and obviously the expectation is what we are looking for: to integrate the incident radiance from all directions on the hemisphere and apply surface interaction terms.

So, in all, I think you should solve the following problems first before you can dive deeper into what the rendering results should actually be:

• Direct illumination: use sample_color += radiance * throughput; instead.
• Path through update: use throughput *= mat.diffuse * sycl::max(dot(closest_hit_info.normal_at_inter, random_dir), 0.0f) / random_direction_pdf; instead.

and since you are using BRDF sampling and emitter sampling at the same time, you can try to use MIS to combine both sampling methods to converge even faster. For more information on the logic for simple path tracer, refer to pbrt-v3: src/integrators/path.cpp:Li(). If you find anything wrong or difficult to understand in this answer, comment to let me know.

By the way, can you point me to the repo for your code? I'd like to try this when I 'm available.

The output of my renderer (surface material setting might be different, they are all Lambertian though), which matches with pbrt-v3 and mitsuba0.6:

• Thanks for the very detailed answer! I haven't had the time to properly integrate the changes you proposed yet so I don't have a feedback on my progress yet. Here's the repo in the meantime though: github.com/TomClabault/SYCL-ray-tracing Commented Oct 18, 2023 at 19:56
• If I understood correctly, the only changes I had to make were the cosine factor for the throughput as well as dividing by the PDF of the bounce direction chosen (and not when accumulating the direct lighting). Basically what you wrote as a summary at the end of your answer. This is the result I get for 128 samples: imgur.com/a/OIsIQen. It does look darker than before but maybe that's how it's supposed to look? Also one detail is still unclear: what's the PDF for the direction chosen of the very first ray? Since it's not sampled, I assume it's 1? Commented Oct 22, 2023 at 20:20
• (1) Yes, you only divide the PDF corresponds to the sampling event: ray direction sampling ---> ray direction sampling PDF, and emitter sampling (shadow ray next event estimation) ---> emitter PDF. (2) The very first ray is normally sampled (if you are using some anti-aliasing methods, like stratified sampling within a pixel) but the PDF is usually uniform so it degrades to a constant scaling factor. btw, mine output (updated in the answer) is much darker and the constrast and color bleeding are obvious, your output... seems quite weird to me. Commented Oct 23, 2023 at 1:20
• You can refer to pbrt-v3 pbrt-v3/src/integrators/path.cpp: Li() line 80 and also my renderer AdaPT/renderer/vanilla_renderer.py: render() line 36 for the complete implementation logic. Commented Oct 23, 2023 at 1:23