<p>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 <code>uniform</code> and <code>cosine-weight</code> 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:</p>
<ul>
<li>Direct illumination estimation:</li>
</ul>
<pre class="lang-cpp prettyprint-override"><code>sample_color += radiance * throughput / random_direction_pdf;
</code></pre>
<p>I don't think this is correct. You see, the <code>radiance</code> 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:</p>
<pre class="lang-cpp prettyprint-override"><code>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
sample_color += radiance * throughput;
// 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 = ...;
</code></pre>
<p>Why? See the following figure:</p>
<p><a href="https://i.sstatic.net/6tHKs.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/6tHKs.png" alt="enter image description here" /></a></p>
<p>The red line segment is shadow ray, and the sampling decision is only related to how you:</p>
<ol>
<li>sample an emitter (if there are multiple)</li>
<li>sample a point from that emitter.</li>
<li>get PDF (area product measure) for sampling the very point on that emitter</li>
<li>Transform to solid angle measure (falloff, cosine attenuation, etc. I think you did it right)</li>
</ol>
<p>Essentially, after <code>if (!in_shadow)</code> branch being finished, the direct illumination already has its unbiased estimation. So you just add it to <code>radiance</code> and account for the path throughput from the previous bounces.</p>
<ul>
<li>Throughput update:</li>
</ul>
<pre class="lang-cpp prettyprint-override"><code>throughput *= mat.diffuse;
</code></pre>
<p>Note that all the subsequent bounces are <strong>sampled</strong> (ray direction sampling and intersection test to obtain new ray origin and direction), you should account for the sampling event here by:</p>
<pre class="lang-cpp prettyprint-override"><code>throughput *= mat.diffuse * sycl::max(dot(closest_hit_info.normal_at_inter, random_dir), 0.0f) / random_direction_pdf;
</code></pre>
<p>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.</p>
<p><a href="https://i.sstatic.net/JXo9R.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/JXo9R.png" alt="enter image description here" /></a></p>
<p>Since our light incident direction (corresponds to <code>random_dir</code>) is sampled, we have the following estimator and we can easily verify its expectation:</p>
<p><span class="math-container">$$
\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
$$</span>
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.</p>
<hr />
<p>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:</p>
<ul>
<li>Direct illumination: use <code>sample_color += radiance * throughput;</code> instead.</li>
<li>Path through update: use <code>throughput *= mat.diffuse * sycl::max(dot(closest_hit_info.normal_at_inter, random_dir), 0.0f) / random_direction_pdf;</code> instead.</li>
</ul>
<p>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 <a href="https://github.com/mmp/pbrt-v3/blob/master/src/integrators/path.cpp" rel="nofollow noreferrer">pbrt-v3: src/integrators/path.cpp:Li()</a>. If you find anything wrong or difficult to understand in this answer, comment to let me know.</p>
<p>By the way, can you point me to the repo for your code? I'd like to try this when I 'm available.</p>
<hr />
<p>The output of my renderer (surface material setting might be different, they are all Lambertian though), which matches with pbrt-v3 and mitsuba0.6:</p>
<p><a href="https://i.sstatic.net/swMGs.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/swMGs.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2023-10-18T19:56:03.057Z
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: https://github.com/TomClabault/SYCL-ray-tracing
Tom Clabault
CC BY-SA 4.0
2023-10-22T20:20:58.590Z
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: https://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?
Tom Clabault
CC BY-SA 4.0
2023-10-23T01:20:27.510Z
(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.
Enigmatisms
CC BY-SA 4.0
2023-10-23T01:23:43.073Z
You can refer to pbrt-v3 [pbrt-v3/src/integrators/path.cpp: Li() line 80](https://github.com/mmp/pbrt-v3/blob/master/src/integrators/path.cpp) and also my renderer [AdaPT/renderer/vanilla_renderer.py: render() line 36](https://github.com/Enigmatisms/AdaPT/blob/master/renderer/vanilla_renderer.py) for the complete implementation logic.
Enigmatisms
CC BY-SA 4.0
2023-10-18T04:49:46.063Z
2023-10-23T01:08:29.997Z
2023-10-23T01:08:29.997Z
Enigmatisms
<p>I'm writing a small path tracer that currently:</p>
<ul>
<li>Samples a random light source at each bounce (direct lighting)</li>
<li>Bounces rays around multiple times (indirect lighting)</li>
<li>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.</li>
</ul>
<p><strong>The issue:</strong> 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 <code>rand_2</code> gets close to 0, see the code below).</p>
<p><strong>Results:</strong></p>
<p>Cosine weighted on the left, uniform distribution over the hemisphere on the right:</p>
<p><a href="https://i.sstatic.net/ilfmr.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ilfmr.jpg" alt="enter image description here" /></a></p>
<p>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).</p>
<p><strong>The code</strong>:</p>
<p><strong>Direction generation (uniform and cosine):</strong></p>
<pre><code>void branchlessONB(const Vector&amp; n, Vector&amp; b1, Vector&amp; 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&amp; normal, const Vector&amp; 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&amp; normal, float&amp; pdf, xorshift32_generator&amp; 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&amp; normal, float&amp; pdf, xorshift32_generator&amp; 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);
}
</code></pre>
<p><strong>Ray bounce part:</strong></p>
<pre><code>for (int bounce = 0; bounce &lt; 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;
                    Vector shadow_ray_direction = random_light_point - shadow_ray_origin;
                    float distance_to_light = length(shadow_ray_direction);
                    Vector shadow_ray_direction_normalized = normalize(shadow_ray_direction);

                    Ray shadow_ray(shadow_ray_origin, shadow_ray_direction_normalized);

                    bool in_shadow = evaluate_shadow_ray(shadow_ray, distance_to_light);

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

                        radiance = emissive_triangle_material.emission;
                        //Cosine angle on the illuminated surface
                        radiance *= sycl::max(dot(closest_hit_info.normal_at_inter, shadow_ray_direction_normalized), 0.0f);
                        //Cosine angle on the light surface
                        radiance *= sycl::max(dot(light_source_info.light_source_normal, -shadow_ray_direction_normalized), 0.0f);
                        //Falloff of the light intensity with the distance squared
                        radiance /= distance_to_light * distance_to_light;
                        //PDF: Probability of having chosen this point on this exact light source
                        radiance /= pdf;
                        //The illuminated surface is Lambertian
                        radiance /= M_PI;
                    }

                    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;
        }
</code></pre>
<p><strong>Light source sampling part:</strong></p>
<pre><code>Point RenderKernel::sample_random_point_on_lights(xorshift32_generator&amp; random_number_generator, float&amp; pdf, LightSourceInformation&amp; 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 &lt;&lt; triangle_area &lt;&lt; sycl::endl;
    pdf = 1 / (nb_emissive_triangles * triangle_area);

    return random_point_on_triangle;
}

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

<pathtracing><importance-sampling>
Cosine weighted hemisphere sampling is a little bit darker and arguably noisier than reference
Tom Clabault
# Cosine weighted hemisphere sampling is a little bit darker and arguably noisier than reference

(Asked by: username_0 on 2023-10-17)

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:

[![enter image description here](https://i.sstatic.net/ilfmr.jpg)](https://i.sstatic.net/ilfmr.jpg)

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;
                    Vector shadow_ray_direction = random_light_point - shadow_ray_origin;
                    float distance_to_light = length(shadow_ray_direction);
                    Vector shadow_ray_direction_normalized = normalize(shadow_ray_direction);

                    Ray shadow_ray(shadow_ray_origin, shadow_ray_direction_normalized);

                    bool in_shadow = evaluate_shadow_ray(shadow_ray, distance_to_light);

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

                        radiance = emissive_triangle_material.emission;
                        //Cosine angle on the illuminated surface
                        radiance *= sycl::max(dot(closest_hit_info.normal_at_inter, shadow_ray_direction_normalized), 0.0f);
                        //Cosine angle on the light surface
                        radiance *= sycl::max(dot(light_source_info.light_source_normal, -shadow_ray_direction_normalized), 0.0f);
                        //Falloff of the light intensity with the distance squared
                        radiance /= distance_to_light * distance_to_light;
                        //PDF: Probability of having chosen this point on this exact light source
                        radiance /= pdf;
                        //The illuminated surface is Lambertian
                        radiance /= M_PI;
                    }

                    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;
}

```

**Tags:** pathtracing, importance-sampling


---

## 1 answers

### Answer by username_1 on 2023-10-18. Score: 2 | Accepted answer

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:

```cpp
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:

```cpp
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
sample_color += radiance * throughput;
// 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:

[![enter image description here](https://i.sstatic.net/6tHKs.png)](https://i.sstatic.net/6tHKs.png)

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:

```cpp
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:

```cpp
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.

[![enter image description here](https://i.sstatic.net/JXo9R.png)](https://i.sstatic.net/JXo9R.png)

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()](https://github.com/mmp/pbrt-v3/blob/master/src/integrators/path.cpp). 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:

[![enter image description here](https://i.sstatic.net/swMGs.png)](https://i.sstatic.net/swMGs.png)

> **Comments:**
> **username_0**: 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: https://github.com/username_0/SYCL-ray-tracing  
> **username_0**: 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: https://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?  
> **username_1**: (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.  
> **username_1**: You can refer to pbrt-v3 [pbrt-v3/src/integrators/path.cpp: Li() line 80](https://github.com/mmp/pbrt-v3/blob/master/src/integrators/path.cpp) and also my renderer [AdaPT/renderer/vanilla_renderer.py: render() line 36](https://github.com/username_1/AdaPT/blob/master/renderer/vanilla_renderer.py) for the complete implementation logic.  

---

