Path tracing - cannot converge diffuse sampling + oversaturation

Question

I'm in process of creating progressive path tracer using DX12 + DXR. I managed to create working raytracing code + pathtracing code with multiple bounces but I'm experiencing problems stated below (after full RT code):

[shader("raygeneration")]
void RayGen()
{
    // Accumulate for limited amount of frames
    if (g_giCB.maxFrames > 0 && g_giCB.accFrames >= g_giCB.maxFrames)
    {
        return;
    }
    uint2 LaunchIndex = DispatchRaysIndex().xy;
    uint2 LaunchDimensions = DispatchRaysDimensions().xy;
    float4 normalAndDepth = NormalTextureInput.Load(int3(LaunchIndex, 0));

    // No geometry hit - skip pixel and use skybox data
    if (normalAndDepth.w == 0)
    {
        RTOutput[LaunchIndex] = albedoTexture.Load(int3(LaunchIndex, 0));
        return;
    }
    
    // Calculate primary ray direction
    uint seed = initRand(LaunchIndex.x + LaunchIndex.y * LaunchDimensions.x, g_sceneCB.frameCount, 16);
    uint seed2 = 0;
    float2 offset = float2(0, 0);
    if (g_giCB.samplingType == SAMPLE_UNIFORM)
    {
        seed2 = initRand(LaunchIndex.x + LaunchIndex.y * LaunchDimensions.x, g_sceneCB.frameCount, 17);
        offset = HammersleyDistribution(g_giCB.accFrames, g_giCB.maxFrames, uint2(seed, seed2));
    }
    else if (g_giCB.samplingType == SAMPLE_MJ)
    {
        const uint pixelIdx = LaunchIndex.y * LaunchDimensions.x + LaunchIndex.x;
        uint sampleSetIdx = 0;
        offset = SamplePoint(pixelIdx, sampleSetIdx);
        seed = pixelIdx;
        seed2 = sampleSetIdx;
    }
    
    float3 primaryRayOrigin = g_sceneCB.cameraPosition.xyz;
    float3 primaryRayDirection;
    GenerateCameraRay(LaunchIndex, LaunchDimensions, g_sceneCB.projectionToWorld, primaryRayOrigin, primaryRayDirection, offset);
    
    // Prepare payload
    PayloadIndirect indirectPayload;
    indirectPayload.color = float3(0, 0, 0);
    indirectPayload.rndSeed = seed;
    indirectPayload.rndSeed2 = seed2;
    indirectPayload.pathLength = 0;
        
    // Calculate pixel color in current pass and merge with previous frames
    float4 finalColor = float4(shootIndirectRay(primaryRayOrigin, primaryRayDirection, 1e-3f, indirectPayload), 1.0f);
    float4 prevScene = RTOutput[LaunchIndex];
    finalColor = ((float) g_giCB.accFrames * prevScene + finalColor) / ((float) g_giCB.accFrames + 1.0f);
    RTOutput[LaunchIndex] = finalColor;
}

[shader("miss")]
void Miss(inout RayPayload payload : SV_RayPayload)
{
    payload.vis = 1.0f;
}

[shader("closesthit")]
void ClosestHit(inout PayloadIndirect payload, in BuiltInTriangleIntersectionAttributes attribs)
{

}

[shader("miss")]
void MissIndirect(inout PayloadIndirect payload : SV_RayPayload)
{
    // Use skybox as contribution if ray failed to hit geometry (right now, disabled for debug purposes)
    float3 rayDir = WorldRayDirection();
    rayDir.z = -rayDir.z;
    if (g_giCB.useSkybox)
    {
        payload.color += skyboxTexture.SampleLevel(g_sampler, rayDir, 0).rgb;
    }
}

[shader("closesthit")]
void ClosestHitIndirect(inout PayloadIndirect payload, in BuiltInTriangleIntersectionAttributes attribs)
{
    // Load hit data
    float3 hitPos = WorldRayOrigin() + WorldRayDirection() * RayTCurrent();
    float3 triangleNormal, triangleTangent, triangleBitangent;
    loadHitData(triangleNormal, triangleTangent, triangleBitangent, attribs);

    // Use white albedo for all textures (DEBUG version)
    float4 albedo = albedoTexture.Load(int3(DispatchRaysIndex().xy, 0));
    albedo = float4(1, 1, 1, 1);
    
    // Iterate over all lights
    float lightsCount = g_lightCB.lightPositionAndType[15].w;
    for (int i = 0; i < lightsCount; i++)
    {
        // Calculate each light data
        float3 lightColor = g_lightCB.lightDiffuseColor[i].rgb;
        float3 toLight = g_lightCB.lightPositionAndType[i].xyz - hitPos;
        float distToLight = length(toLight);
        toLight = normalize(toLight);
        
        // Check visibility
        float NoL = saturate(dot(triangleNormal.xyz, toLight));
        float visibility = shadowRayVisibility(hitPos, toLight, 1e-3f, distToLight);

        // Calculate light contribution to point in world (diffuse lambertian term)
        payload.color += visibility * NoL * albedo.rgb * INV_PI;
    }
    
    if (g_giCB.useIndirect == 1)
    {
        // Continue spawning rays if path left has not reached maximum
        if (payload.pathLength < g_giCB.bounceCount)
        {
            // Find next direction
            float3 rayDirWS = float3(0, 0, 0);
            if (g_giCB.samplingType == SAMPLE_UNIFORM)
            {
                float3x3 tangentToWorld = float3x3(triangleTangent, triangleBitangent, triangleNormal);
                float2 hammersley = HammersleyDistribution(payload.pathLength, g_giCB.bounceCount, uint2(payload.rndSeed, payload.rndSeed2));
                float3 rayDirTS = UniformSampleHemisphere(hammersley.x, hammersley.y);
                rayDirWS = normalize(mul(rayDirTS, tangentToWorld));
            }
            else if (g_giCB.samplingType == SAMPLE_MJ)
            {
                float3x3 tangentToWorld = float3x3(triangleTangent, triangleBitangent, triangleNormal);
                float2 brdfSample = SamplePoint(payload.rndSeed, payload.rndSeed2);
                float3 rayDirTS = SampleDirectionCosineHemisphere(brdfSample.x, brdfSample.y);
                rayDirWS = normalize(mul(rayDirTS, tangentToWorld));
            }
            else if (g_giCB.samplingType == SAMPLE_RANDOM)
            {
                rayDirWS = getCosHemisphereSample(payload.rndSeed, triangleNormal, triangleTangent, triangleBitangent);
                nextRand(payload.rndSeed);
            }
            
            // Prepare payload
            PayloadIndirect newPayload;
            newPayload.pathLength = payload.pathLength + 1;
            newPayload.rndSeed = payload.rndSeed;
            newPayload.rndSeed2 = payload.rndSeed2;
            newPayload.color = float3(0, 0, 0);
            
            // Calculate next ray bounce color contribution
            float3 bounceColor = shootIndirectRay(hitPos, rayDirWS, 1e-3f, newPayload);
            payload.color += bounceColor * albedo.rgb;
        }
    }
}

1. Image doesn't converge to ground truth over time

I am using only Lambertian BRDF term in my code, while implementing multi sampling schemes (Multi-Jittered sampling, Uniform sampling and random direction sampling). All of them provides similar results, hardly possible to distinguish. Here is an image generated with almost 5000 frames. It doesn't differ from image generated with about 100 frames. I am using albedo = (1,1,1) for all textures, for debug purposes in this image:

It looks even much more noisy when source of light is outside the window and the you can can never access light directly except primary ray (situation similar to described here - http://www.pbr-book.org/3ed-2018/Light_Transport_III_Bidirectional_Methods/Bidirectional_Path_Tracing.html ). Therefore I placed point light in the middle of room for simpler example.

I'm not suspecting my equation to be incorrect, because: $$ L_o(\textbf{p}, \mathbf{w_o}) = L_e + \int_{\Omega} L_i(\textbf{p}, \mathbf{w_i}) fr(\mathbf{w_o}, \mathbf{w_i}) \cos \theta d\omega $$

$$ \frac{1}{N} \sum_{k=1}^{N} \frac{ L_i(\textbf{p}, \mathbf{w_k}) fr(\mathbf{w_k}, w_o) \cos \theta }{p(\mathbf{w_k})} $$

Monte Carlo equation above will be simplified, because I have single light with intensity = 1.0, so $L_i$ term will be terminated. BRDF for Lambertian diffuse is equal to NdotL (I'll be using notation $NoL$). PDF for sampling cosine hemisphere is $NoL / \pi$. Let's simplify equation for now:

$$ \frac{1}{N} \sum_{k=1}^{N} \frac{ \cos \theta }{\pi} $$

Which is basically: $$ \frac{1}{N} \sum_{k=1}^{N} \frac{ NoL }{\pi} $$

And translates to code as:

color = NoL * INV_PI;

However, we need to check if light $L$ that we're referring to is visible, therefore full equation, which could be seen in code is:

// Calculate light contribution to point in world (diffuse lambertian term)
payload.color += visibility * NoL * albedo.rgb * INV_PI;

Regarding sampling, I've tried multiple method, especially I was hopeful about Multi-Jittered sampling [Kensler2013]. I am using native implementation, copied from MJP's Path Tracer.

2. Oversaturation of albedo textures

Primary ray + single bounce image with albedo textures

Primary ray + 4 ray bounces image; Notice oversaturation of albedo textures

If I simplify code to the maximum, then ray bounces boils down to very few lines of code used in practice. The problem might be adding albedo too many times but I checked it for sure, and both MJP and cwyman are using same code for calculating path tracing (multiple bounces of rays).

[shader("closesthit")]
void ClosestHitIndirect(inout PayloadIndirect payload, in BuiltInTriangleIntersectionAttributes attribs)
{
    float4 albedo = // ...
    
    // Use single light
    float3 toLight = g_lightCB.lightPositionAndType[0].xyz - hitPos;
    float distToLight = length(toLight);
    toLight = normalize(toLight);
        
    // Check visibility
    float NoL = saturate(dot(triangleNormal.xyz, toLight));
    float visibility = shadowRayVisibility(hitPos, toLight, 1e-3f, distToLight);

    // Calculate light contribution to point in world (diffuse lambertian term)
    payload.color += visibility * NoL * albedo.rgb * INV_PI;
    
    if (g_giCB.useIndirect == 1)
    {
        // Continue spawning rays if path left has not reached maximum
        if (payload.pathLength < g_giCB.bounceCount)
        {
            // Find next direction
            float3 rayDirWS = // ... multiple lighting schemes, check full code above
            
            // Prepare payload
            PayloadIndirect newPayload;
            newPayload.pathLength = payload.pathLength + 1;
            newPayload.rndSeed = payload.rndSeed;
            newPayload.rndSeed2 = payload.rndSeed2;
            newPayload.color = float3(0, 0, 0);
            
            // Calculate next ray bounce color contribution
            float3 bounceColor = shootIndirectRay(hitPos, rayDirWS, 1e-3f, newPayload);
            payload.color += bounceColor * albedo.rgb;
        }
    }
}

For now, I'd like to creating converging diffuse path tracer. I'll continue with Specular GGX in order to have a ground truth reference for Bidirectional Path Tracing and other methods, which can speed up converging rate of an image. Thanks in advance for all your advices.

Also, for future readers - I'd recommend checking questions below, which are related to the topic, but didn't manage to solve my problems:

• Is it expected that a naive path tracer takes many, many samples to converge?
• Non Real Time Raytracing
• Path tracing the Cook-Torrance BRDF

---

Edit after Nathan's Reed answer:

The oversaturated image looks better now. It looks a bit dull after moving from linear to sRGB space (following Nathan's order of exposure -> tone mapping -> linear_To_sRGB):

Primary ray + 4 ray bounces image; No skybox miss shader applied to fair comparison to oversaturated image before

However, I checked case with light source from the outside and I really like the result. There is a room for improvement, for sure, like adding some AA, moving to PBR (GGX). I doubt I'll have a time for that, but using BDPT and implementing soft shadows would definatelly improve case shown below:

To sum up for future reference: adding tone mapper (on top of existing exposure settings) and moving from linear to sRGB space helped with oversaturation. Using R16G16B16A16_FLOAT instead of R8B8G8A8_UNORM format for backbuffer and target output textures, solved problem with not converging images. In the future, R32G32B32A32_FLOAT might be used, but with my current DX12 setup, this format causes errors when used in backbuffer and I cannot change that for now. 4x16 format proved to be enough for now, so I'll stay with that.

Answer 1

One point: the Lambert BRDF is not $N\cdot L$, it's just the albedo divided by pi. The $N \cdot L$ factor comes from the $\cos \theta$ in the rendering equation. So, when sampling with a cosine-weighted distribution the $N \cdot L$s and pis will cancel out and you should just be accumulating $\frac{1}{N} \sum L_i * \text{albedo}$.

It looks like in your code you're doing this correctly for the cosine-weighted distribution. But then in the case of the uniform hemispherical distribution you're not putting in the needed factor of $2(N \cdot L)$. (The factor of 2 comes from dividing by the uniform hemispherical pdf of 1/2π, and canceling the 1/π in the BRDF.) So I think that your uniform sampling method does not give quite correct results here.

I'm not quite sure what you meant by "image doesn't converge to ground truth". Do you mean that it's still noisy even after quite a lot of samples? The problem could be in your accumulation strategy. Depending on the bit depth of your framebuffer, after some number of frames the contribution of one more frame will be so small as to be smaller than the least-significant bit of the accumulated value. If you're using 16-bit float format for instance this will happen after about 1000 frames; if using R11G11B10_FLOAT then it will take only 30-60 frames before new accumulated values will no longer have any effect. The accumulation framebuffer should be 32-bit float at a minimum.

With multiple bounces, it doesn't look like you're accounting for throughput along a path correctly. You have payload.color += bounceColor * albedo.rgb, but note that this only takes account of the albedo at the current surface. If this is the Nth bounce, the color should be multiplied by the albedos of all previous surfaces in the path—since that's how this light is getting to the camera, by bouncing through all of them. The path payload structure needs to include not just an accumulated color, but a value typically called "throughput", that contains the product of all the (BRDF * cos(theta) / pdf) factors along the path so far. On each bounce you update the throughput by multiplying in the factor for the latest bounce, then update the color by multiplying the sampled radiance with the accumulated throughput.

Also, do you have any exposure control / tonemapping on the final image? An image with multiple bounces is expected to be brighter overall than the same scene with 1 bounce, as more light is being accumulated. If you don't adjust the exposure and apply a tone curve of some sort (and gamma correction), you can end up with things looking bad/wrong on screen, even if the internal HDR framebuffer is correct.