Warped Area Sampling Harmonic Weight Computation

This may be long and a little involved. I will try my best to explain the possible issue.

It's an issue of implementation detail of harmonic weight used by warped-area sampling (WAS) differentiable path tracer. Here is the paper

I am talking about this function: float computeHarmonicWeight() in Source\RenderPasses\WARDiffPathTracer\WarpedAreaReparam.slang This is the current implementation:

float computeHarmonicWeight(
    no_diff IntersectionAD isect,
    no_diff float3 origin,
    no_diff float3 auxDirection,
    no_diff float auxSampleY,
    no_diff float kappa,
    float3 direction
)
{
    float boundaryTerm = no_diff computeBoundaryTerm(isect.normalW, auxDirection);

    float sy = max(1.f - auxSampleY, 1e-6f);
    float invVMFDensity = 1.f / ((1.f - sy) * exp(-2.f * kappa) + sy);

    float wDenom = invVMFDensity - 1.f + boundaryTerm;
    float wDenomRcp = select(wDenom > 1e-4f, 1.f / wDenom, 0.f);
    float w = wDenomRcp * wDenomRcp * wDenomRcp * invVMFDensity;
    return w;
}

This function basically implements the computation of weight_i, a single line in Algorithm 2 in the paper. It seems like the code doesn't match the formula given in the paper. Here is some derivation & explanation. (I will use 'w' for omega, unit direction and 'weight' for w in the paper for convenience)

dot(w, w'). We take advantage of a fact that dot product doesn't change after we transform both vectors with same transformation. Thus, it equals the dot product in local space, where w_local is simply (0,0,1) and the sampled w'_local has z component cos(theta) = 1.f + log(exp(-2.f kappa) (1.f - sy) + sy) / kappa; See sampleVonMisesFisher() in Source\Falcor\Utils\Math\MathHelpers.slang for details.
Following that, after some simple derivation, we have the exponential term in the denominator becomes (1.f - sy) exp(-2.f kappa) + sy, which should be "VMF density of w' ".
Thus, weight_i = 1 / (VMF_density - 1 + BoundaryTerm) / VMFDensity
Further, I couldn't understand why there is a "1 / wDenom ^ 3". I didn't find anything in the calculation related to a cubic term.

And here is my modified implementation:

float computeHarmonicWeight(
    no_diff IntersectionAD isect,
    no_diff float3 origin,
    no_diff float3 auxDirection,
    no_diff float auxSampleY,
    no_diff float kappa,
    float3 direction // not used in computation; allows autodiff to find grad w.r.t. this direction
)
{
    float boundaryTerm = no_diff computeBoundaryTerm(isect.normalW, auxDirection);

    float sy = max(1.f - auxSampleY, 1e-6f);
    // see WAS paper Algorithm 2, computation of weight_i
    /* float invVMFDensity = 1.f / ((1.f - sy) * exp(-2.f * kappa) + sy);

    float wDenom = invVMFDensity - 1.f + boundaryTerm;
    float wDenomRcp = select(wDenom > 1e-4f, 1.f / wDenom, 0.f);
    float w = wDenomRcp * wDenomRcp * wDenomRcp * invVMFDensity; */

    float VMFDensity = (1.f - sy) * exp(-2.f * kappa) + sy;
    float w = 1.f / (VMFDensity - 1.f + boundaryTerm) / VMFDensity;
    return w;
}

After the change, here is the bunny_ref test scene. (First is before change, second is after change.) I did see some very subtle difference. But admittedly, I am not an expert on Falcor, differentiable rendering, or WAS, so it's very possible I am wrong. But I do want to point things out in case there is something wrong.

Thank you!

BeforeChange

AfterChange

Finally, here is the launch.json entry for this test scene for your convenience.

        {
            // Launch configuration for Mogwai (Warped Area Sampling Differentiable Rendering)
            "name": "Mogwai WAS",
            "type": "cppvsdbg",
            "request": "launch",
            "program": "${command:cmake.launchTargetPath}",
            "windows": {
                "program": "${command:cmake.launchTargetDirectory}/Mogwai.exe"
            },
            "args": [
                "--script=${workspaceFolder}/scripts/WARDiffPathTracer.py",
                "--scene=${workspaceFolder}/media/inv_rendering_scenes/bunny_ref.pyscene"
            ],
            "stopAtEntry": false,
            "cwd": "${command:cmake.launchTargetDirectory}",
            "environment": [
                {
                    "name": "FALCOR_DEVMODE",
                    "value": "1"
                }
            ],
            "visualizerFile": "${workspaceFolder}/Source/Falcor/Falcor.natvis"
        },

NVIDIAGameWorks / Falcor

Warped Area Sampling Harmonic Weight Computation #426