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!
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"
},