Could we Kiss Multiple Importance Sampling Good-Bye? Implementing "Practical Product Sampling by Fitting and Composing Warps"

The cover image is a bait, its obviously not the MIS-free render I'm advocating. Also forgive the color banding visible in some images, they don't come from our production renderer where tonemapping and blue noise dithering is applied, they come from a single shader unit test.

TL;DR Area Light made of a triangle, 8 spp, MIS disabled. Projected Solid Angle, followed by Solid Angle and finally Area Sampling.

Forgive the white sphere, its a dielectric that produces NaNs when MIS is disabled, as my renderers are actually not designed to be invoked without MIS.

If you already know the intricacies of Practical Product Sampling, then skip the first three sections, all the way past the implementation. However the source code to an optimized implementation may interest you. If you're a user of Unity or Unreal engines and only care about realtime area lights, do not close this tab, as what I'm about to present will tie into shadows from Linearly Transformed Cosine Area Lights. Same applies if you don't care about Triangular Area Lights.

What is it?

My previous article, covers how to validate that your sampling function actually produces samples that match the claimed PDF in an idiot proof way. I have also derived the relationship between the determinant of the Jacobian of the warp, the pdf of the input random samples and the target output pdf. At that time I was unaware of the Practical Product Sampling paper.

One of the immediate reactions to that article by a PhD student was "I wonder if this could be used to derive mappings", and not long after, while familiarising myself with the reasearch topics of a Computer Graphics lab that I intend to apply to for a PhD, I've ran into the magnificient paper from Eurographics 2020, "Practical Product Sampling by Fitting and Composing Warps".

Its amazing that it took 25 years and the industry heavy-weights (D. Hart, M. Pharr, T. Müller, W. Lopes, M. McGuire and P. Shirley) to come up with something as simple and as briliant as Practical Product Sampling, a nice general method that can solve many product sampling problems, but first and foremost the problem of sampling the Projected Solid Angle of a Triangle without sample rejection, within bounded time and memory (which is not what we can claim of the previous Arvo and Urena papers).

The premise is ridiculously simple, you pass your uniformly distributed samples through more than one warp and take advantage of the fact that the Jacobian of the combined warp is all the separate warp's Jacobians multiplied together (chain rule).

This seems a useless tautology because one of your warps would need to be very specific, the authors give a single warp (plain importance sampling) example in supplemental material video presentation to demonstrate one would end up with the same problem as they started with.

However, we can note that any of the warps does not need to correspond exactly to the exact conditional PDF we are after. So if we can reasonably approximate one of the conditional probability functions, or rather, it's warp, then we can end up with a vastly superior sampling strategy by combining it with the original warp[s] than the original alone.

If the approximation is "reasonable enough" then it's almost as close to perfect product sampling as we can get when measuring efficiency of computational resources spent vs. quality. The authors did actually benchmark their approximation against perfect precomputed product sampling and the difference in MRSE was negligible in many cases, which suggests the Visibility Factor (shadowing) now has a far greater influence on variance. At the end of my article I will argue how we could make the MRSE difference between Practical Product Sampling and Exhaustive BSDF sampling negligible, always.

Projected Solid Angle Triangle Sampling

In this case, we approximate the Projected Solid Angle over a Spherical Triangle by bilinear interpolation in the Primary Sample Space (a complicated term for denoting your [0,1]^2 sample domain). The authors also explored other illuminators such as rectangles and bilinear patches.

Given a current state of the art sampling strategy for Triangular Light Sources (Arvo 1995), we can warp its input samples which are in the Primary Sample Space according to a Bilinearly Interpolated PDF which approximates the `cos(theta)` or `dot(N,L)` term.

It makes far more sense to first warp according to the conditional probability of the projected solid angle given a spherical triangle, and then proceed with Spherical Triangle Sampling as usual, rather than first carry out exact or approximate Spherical Triangle Sampling followed by a non trivial warp of directions on a sphere.

If we imagine for a second that our approximation of projected solid angle over the surface of a spherical triangle equals the actual projected solid angle, then the determinant of our Jacobian of the composed warp would match the product of the PDFs.

Turns out that bilinear interpolation in Primary Sample Space warped onto a spherical triangle is quite similar to approximating the cosine function in the -PI/2 to PI/2 range with a tent function. So its a "reasonable enough" approximation.

Implementation (for a Triangle)

We need an implementation of Spherical Triangle Sampling to start off with, such as the one from my previous article. The following implementation is made available on the same license as the previous, basically, you can't sue me if it doesn't work and you have to let people know that I've saved you some valuable time.

/**
Copyright 2020 Mateusz "DevSH" Kielan

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    https://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
**/

All that needs to change for light sampling is the sample generation function:

// There are two different modes of sampling, one for BSDF and one for BRDF (depending if we throw away bottom hemisphere or not)
vec3 irr_glsl_sampling_generateProjectedSphericalTriangleSample(out float rcpPdf, in mat3 sphericalVertices, in vec3 receiverNormal, in bool isBSDF, vec2 u)
{
    // pre-warp according to proj solid angle approximation
    u = irr_glsl_sampling_generateBilinearSample(rcpPdf,irr_glsl_sampling_computeBilinearPatchForProjSphericalTriangle(sphericalVertices,receiverNormal,isBSDF),u);

    // now warp the points onto a spherical triangle
    float solidAngle;
    const vec3 L = irr_glsl_sampling_generateSphericalTriangleSample(solidAngle,sphericalVertices,u);
    rcpPdf *= solidAngle;

    return L;
}
// utility function
vec3 irr_glsl_sampling_generateProjectedSphericalTriangleSample(out float rcpPdf, in mat3 vertices, in vec3 origin, in vec3 receiverNormal, in bool isBSDF, in vec2 u)
{
    return irr_glsl_sampling_generateProjectedSphericalTriangleSample(rcpPdf,irr_glsl_shapes_getSphericalTriangle(vertices,origin),receiverNormal,isBSDF,u);
}

You can see three most important differences:

1. The pre-warping of the Primary Sample Space sample coordinate `u`
2. The computation of the PDF
3. Suddenly we need to keep track of the receiver's Surface Normal and if it is a BSDF or BRDF (does it transmit)

We'll work through them backwards.

The Surface Normal is needed to know the Solid Angle Projection plane, before we were only concerned about generating points in the domain of our integration uniformly, now we're concerned about an integral that contains something more than dOmega . We also need to know if we're dealing with a BRDF or a BSDF, i.e. whether we want to throw away the parts of the triangle below the hemisphere (no transmission) or not (transmission). All of this only affects our bilinear approximation of the Projected Solid Angle with respect to the Solid Angle.

The PDF will be inversely proportional to the product of the determinants of the warps' Jacobians, so it follows that the PDF will be a multiplication of the individual warp's PDFs. However pay attention at the coordinates at which the PDFs are evaluated, you want to multiply the PDF of the sample that has been generated by each successive warp. In this case its not immediately apparent because the PDF of a Spherical Triangle is a constant. Also because the pre-warp approximates the differential Projection of the differential Solid Angle we cannot factor out the cosine term from the Light Transport Equation's Integral. The cosine factor will need to be computed as usual and the PDF of the pre-warp will be slightly different.

The pre-warping is relatively straight forward, we take the four corners of our Primary Sample Space, evaluate the conditional probability of the Projected Solid Angle given the Solid Angle, and we generate samples as-if for the PDF resulting from a normalized Bilinear Interpolation function. Because Spherical Triangle Sampling produces a constant differential Solid Angle from a uniform differential Primary Sample Space, the aforementioned conditional probability is just the clamped or absolute value of the cosine between the Surface Normal and direction towards the sample.

vec4 irr_glsl_sampling_computeBilinearPatchForProjSphericalTriangle(in mat3 sphericalVertices, in vec3 receiverNormal, in bool isBSDF)
{
    // a positive value would prevent us from a scenario where `irr_glsl_sampling_rcpProbBilinearSample` will return NAN
    const float minimumProjSolidAngle = 0.0;

    vec3 bxdfPdfAtVertex = transpose(sphericalVertices)*receiverNormal;
    // take abs of the value is we have a BSDF
    bxdfPdfAtVertex = uintBitsToFloat(floatBitsToUint(bxdfPdfAtVertex)&uvec3(isBSDF ? 0x7fFFffFFu:0xffFFffFFu));
    bxdfPdfAtVertex = max(bxdfPdfAtVertex,vec3(minimumProjSolidAngle));

    // the swizzle needs to match the mapping of the [0,1]^2 square to the triangle vertices
    return bxdfPdfAtVertex.yyxz;
}
// The square's vertex values are defined in Z-order, so indices 0,1,2,3 (xyzw) correspond to (0,0),(1,0),(0,1),(1,1)
vec2 irr_glsl_sampling_generateBilinearSample(out float rcpPdf, in vec4 bilinearCoeffs, vec2 u)
{
    const vec2 twiceAreasUnderXCurve = vec2(bilinearCoeffs[0]+bilinearCoeffs[1],bilinearCoeffs[2]+bilinearCoeffs[3]);
    u.y = irr_glsl_sampling_generateLinearSample(twiceAreasUnderXCurve,u.y);

    const vec2 ySliceEndPoints = vec2(mix(bilinearCoeffs[0],bilinearCoeffs[2],u.y),mix(bilinearCoeffs[1],bilinearCoeffs[3],u.y));
    u.x = irr_glsl_sampling_generateLinearSample(ySliceEndPoints,u.x);

    rcpPdf = (twiceAreasUnderXCurve[0]+twiceAreasUnderXCurve[1])/(4.0*mix(ySliceEndPoints[0],ySliceEndPoints[1],u.x));

    return u;
}

This is where you might trip up in the implementation, if you didnt use the same notation Arvo uses or just ripped your implementation off the internet, since depending which arcs of the Spherical Triangle you choose to split (and in what direction) or how to label them, the corners of the Primary Sample Space will correspond to different triangle vertices.

I do not concern myself at all with the fact that not clamping the cosine term to a positive minimum, choosing instead to clamp to 0, will result in NaN PDFs for the warp, this is because the only time that will happen is when the receiver is on the triangle's horizon or vice versa or if the receiver has a BRDF material and the triangle is below the receiver's horizon. All of these cases are when we absolutely do not want to trace a ray towards the light because the light contribution will be null. All you need to do is to keep you NEE ray trace guard as `if (rcpPdf<FLT_MAX)`.

You can already notice that this approximation will (sometimes grossly) underestimate the ideal PDF for a case where the triangle straddles the zenith of the clamped cosine function, and it will overestimate when the triangle straddles the receiver's horizon, furthermore degenerating into plain Solid Angle sampling for receivers with transmissive BSDFs. The authors of the paper suggest using a Biquadratic fit instead of Bilinear which however requires solving a cubic instead of a quadratic equation to generate samples.

float irr_glsl_sampling_generateLinearSample(in vec2 linearCoeffs, in float u)
{
    const float rcpDiff = 1.0/(linearCoeffs[0]-linearCoeffs[1]);
    const vec2 squaredCoeffs = linearCoeffs*linearCoeffs;
    // 1.0/0.0 gives inf
    return abs(rcpDiff)<FLT_MAX ? (linearCoeffs[0]-sqrt(mix(squaredCoeffs[0],squaredCoeffs[1],u)))*rcpDiff:u;
}

I believe a potentially better solution would be to split the Bilinear fit into a Piecewise Bilinear fit with 4 separate patches and at least one adaptively placed control point (i.e. the horizon or the zenith if the triangle straddles it). Maybe at some point I'll look into it.

Let's say I want/need Multiple Importance Sampling

Now you didn't think I've implemented this sampling without MIS, did you?

Here's a GIF comparison at a much lower spp of 4 with a worst case split of 3 BSDF samples to 1 NEE.

To get the PDF of a given ray/sample, we need to retrace our steps in the opposite direction. All we need to do is apply the steps in reverse. Rember that you have no reason to call the following function unless the closest hit ray you've just obtained the intersection for, intersected an actual triangular light that could have been sampled with the methods outlined above.

If you are paying attention you will notice the memory cost of this technique, which is that we need to store some information about the receiver in the ray's payload, something we've never had to do before. This is because when we want to know an already generated light sample's PDF its usually exclusively for MIS after we've hit some light with our BSDF or Path Guiding generated ray.

float irr_glsl_sampling_probProjectedSphericalTriangleSample(in mat3 sphericalVertices, in vec3 receiverNormal, in bool receiverWasBSDF, in vec3 L)
{
    float pdf;
    const vec2 u = irr_glsl_sampling_generateSphericalTriangleSampleInverse(pdf,sphericalVertices,L);

    return pdf*irr_glsl_sampling_probBilinearSample(irr_glsl_sampling_computeBilinearPatchForProjSphericalTriangle(sphericalVertices,receiverNormal,receiverWasBSDF),u);
}

Obtaining the PDF of the of the bilinear function is almost trivial, its simply evaluating its normalized version, but I'll list it anyway. This is where paying attention to the fact that you can multiply the warp's associated PDFs but evaluated at warped locations, really pays its dividends (like not having to worry about the inverse warp of the bilinear).

float irr_glsl_sampling_probBilinearSample(in vec4 bilinearCoeffs, vec2 u)
{
    return 4.0*mix(mix(bilinearCoeffs[0],bilinearCoeffs[1],u.x),mix(bilinearCoeffs[2],bilinearCoeffs[3],u.x),u.y)/(bilinearCoeffs[0]+bilinearCoeffs[1]+bilinearCoeffs[2]+bilinearCoeffs[3]);
}

The only problem you'll encounter is the inverse of the Spherical Triangle Sampling function, its not something you can just easily find on the web outside of PBRT v4, and you definitely don't want as many `cross` calls, which is okay for PBRT because its a educational renderer, made and intended to be legible.

vec2 irr_glsl_sampling_generateSphericalTriangleSampleInverse(out float pdf, in mat3 sphericalVertices, in vec3 L)
{
    // for angles between view-to-vertex vectors
    float cos_a,cos_c,csc_b,csc_c;
    // Both vertices and angles at the vertices are denoted by the same upper case letters A, B, and C. The angles A, B, C of the triangle are equal to the angles between the planes that intersect the surface of the sphere or, equivalently, the angles between the tangent vectors of the great circle arcs where they meet at the vertices. Angles are in radians. The angles of proper spherical triangles are (by convention) less than PI
    vec3 cos_vertices,sin_vertices;
    // get solid angle, which is also the reciprocal of the probability
    pdf = 1.0/irr_glsl_shapes_SolidAngleOfTriangle(sphericalVertices,cos_vertices,sin_vertices,cos_a,cos_c,csc_b,csc_c);

    // get the modified B angle of the first subtriangle by getting it from the triangle formed by vertices A,B and the light sample L
    const float cosAngleAlongBC_s = dot(L,sphericalVertices[1]);
    const float csc_a_ = inversesqrt(1.0-cosAngleAlongBC_s*cosAngleAlongBC_s); // only NaN if L is close to B which implies v=0
    const float cos_b_ = dot(L,sphericalVertices[0]);

    const float cosB_ = (cos_b_-cosAngleAlongBC_s*cos_c)*csc_a_*csc_c; // only NaN if `csc_a_` (L close to B) is NaN OR if `csc_c` is NaN (which would mean zero solid angle triangle to begin with, so uv can be whatever)
    const float sinB_ = sqrt(1.0-cosB_*cosB_);

    // now all that remains is to obtain the modified C angle, which is the angle at the unknown vertex `C_s`
    const float cosC_ = sin_vertices[0]*sinB_*cos_c-cos_vertices[0]*cosB_; // if cosB_ is NaN then cosC_ doesn't matter because the subtriangle has zero Solid Angle (we could pretend its `-cos_vertices[0]`)
    const float sinC_ = sqrt(1.0-cosC_*cosC_);

    const float subTriSolidAngleRatio = irr_glsl_getArccosSumofABC_minus_PI(cos_vertices[0],cosB_,cosC_,sin_vertices[0],sinB_,sinC_)*pdf; // will only be NaN if either the original triangle has zero solid angle or the subtriangle has zero solid angle (all can be satisfied with u=0)
    const float u = subTriSolidAngleRatio>FLT_MIN ? subTriSolidAngleRatio:0.0; // tiny overruns of u>1.0 will not affect the PDF much because a bilinear warp is used and the gradient has a bound (won't be true if LTC will get used)

    // INF if any angle is 0 degrees, which implies L lays along BA arc, if the angle at A is PI minus the angle at either B_ or C_ while the other of C_ or B_ has a zero angle, we get a NaN (which is also a zero solid angle subtriangle, implying L along AB arc)
    const float cosBC_s = (cos_vertices[0]+cosB_*cosC_)/(sinB_*sinC_);
    // if cosBC_s is really large then we have numerical issues (should be 1.0 which means the arc is really short), if its NaN then either the original or sub-triangle has zero solid angle, in both cases we can consider that the BC_s arc is actually the BA arc and substitute
    const float v = (1.0-cosAngleAlongBC_s)/(1.0-(cosBC_s<uintBitsToFloat(0x3f7fffff) ? cosBC_s:cos_c));

    return vec2(u,v);
}

First part proceeds as usual, because you need to know the original Solid Angle of the triangle. Then the juicy meat of the code starts, we need to find the Solid Angle of the sub-triangle produced by the first subdivision along the arc A to C to obtain our first unwarped coordinate which should be a ratio of the subtriangle's Solid Angle to the original triangle's Solid Angle.

Now this is where people would usually try to reconstruct the original `C_s` vertex along the arc A to C, I've skipped that step altogether by applying the Spherical Triangle Law of Cosines (the form that uses vertex angles instead of side arclengths) to find the subtriangle's angle at vertex B. It is then relatively simple to deduce the angle at vertex `C_s`. The whole process you can find on Wikipedia under "Solving Spherical Triangles", its the Angle-Side-Angle case (ASA). No cross products used! `C_s` remains unknown.

To obtain our second coordinate we perform a pretty easy reverse computation of the arguments to the last slerp in the forward warp function, since we know the cosine of the angle of the sampled direction to the B vertex of the original triangle and can apply the Cosine Law to deduce the arclength between B and C_s the unknown vertex of the first subtriangle.

Of course this process can produce NaN but this will only happen in the following cases:

1. The Triangle was below the Receiver's Horizon
2. The Receiver was on the Triangle's Horizon
3. The Triangle's Solid Angle was incredibly small

So all you need to check is basically if the pdf is greater than a threshold before deciding to weigh a non-NEE sample with MIS that takes the NEE generator probability into account.

As for the NaNs and INFs produced during the inverse warp computation, we protect the `u` and `v` coordinates against them, see the analysis in the code's comments.

Caveat - It looks like you still need MIS?

Ofcourse we can only generate samples on the surface of an emissive triangle tagged as a light, to get any indirect lighting and direct lighting from surfaces not tagged as lights we will still need to cast BSDF rays and NEE rays. There is however a choice of whether we apply MIS or not.

We're often told that MIS is essentially free, its a partial myth, computing the PDFs of samples generated with other other technique's generators is not free. Even if we consider the "free" to mean "no extra rays cast", we either need to remember the samples we've generated with the different generators which leads to needing a stack-like mechanism to share path prefixes, or if we choose (or are forced by our sampler not tolerating correlated path prefixes) to go stackless, then we need to stochastically choose between a BSDF or NEE sample which is no longer "free". Furthermore the "indirect" lighting part of the light transport does not benefit from the NEE rays at all, because the probability of generating an NEE sample on the surface of a "non-light" is 0 and allowing the NEE rays to account for hits of other lights would make the NEE PDF impractical to compute (you'd need to know all the lights that overlap along a ray, which is an exhaustive search instead of a closest or any hit).

The situation gets even more cumbersome if we throw a third strategy such as Path Guiding into the mix, so what we'd really want to do is collapse this back into just two strategies, for example the Product of BSDF and Light Distribution plus the Path Guiding. I'll actually get back to it at the very end.

Obviously due to the bilinear or biquadtratic interpolation, especially looking at the quality of the highlights on the relatively smooth spheres in my render, it seems like we cannot do Practical Sampling of the product of the Light Distribution and the BSDF.

You couldn't be more wrong! Linearly Transformed Cosines [Eric Heitz, Jonathan Dupuy, Stephen Hill and David Neubelt] to the rescue!

We can fit one or more (two necessary for BSDFs) Clamped Cosine Distributions to any popular BSDF with relatively little error. The fit is usually off-line and for isotropic BSDFs it can be stored in a small texture indexed by the angles of Incoming and Outgoing directions to the Surface Normal, analytical fits could also be possible.

This means I don't need MIS or specialized BSDF sampling, as I can just fit an LTC to any BSDF lobe no matter the roughness!

Not only does this speed everything up, because a distribution given by a 3x3 matrix multiply and a dot product are insanely faster to evaluate and sample than all of the non-zero roughness BSDFs, but also because this distribution is probably the only currently known non-uniform distribution on a sphere which has trivial integrals over arbitrary domains (thanks to Lambert developing Irradiance Tensors in the 19th century and Stokes'/Green's theorems about Contour Integrals), especially over polygons where they're a simply sum of arccosines (but not the same ones as to obtain the Solid Angle). Spherical Caps also have a nice solution, which is equivalent to sampling Projected Spherical Caps. This property is so remarkable it allows for non-sampled real-time implementations of Area Lighting, albeit without correct shadows, from arbitrary polygons, spheres and thin lines have been implemented in Unity or even WebGL. The only implementation detail is that you need to clip or split your polygons or other domains such as spherical caps at the horizon.

However during their presentation's Q&A one of the authors acknowledged that this technique does not yet bring anything groundbreaking to the table for Monte Carlo rendering in terms of variance reduction because it could not be sampled trivially for polygons, so one would need to revert to either just sampling the BSDF (which, granted, is a lot faster) or the spherical domain the light covers (at best Solid Angle sampling, at worst area sampling) or both combined with MIS.

With "Practical Product Sampling" we have the fast technique capable of sampling the Linearly Trasformed Cosine over an arbitrary polygon, and with "Sampling Projected Spherical Caps" we have the technique capable of sampling the LTC over a sphere. We'd simply add another warp which is the LTC's associated matrix to transform the triangle's vertices before the usual sampling routine.

This is my whole argument for my initial claim that the MRSE difference between Practical Product Sampling and Exhaustive BSDF sampling could always be negligible, if the sampling was ameliorated with an application of LTC.

It's so powerful that combined with technologies like KHR_raytracing we could have LTC area lighting without disturbingly hacked and abused shadows in real-time, failing that we could at least finally use all LTC's benefits for offline rendering.

This now allows us to abandon MIS completely if we chose to, which we in reality would never want to do because MIS helps with the Visibility (shadows) term variance, especially if we are able to cast BSDF and NEE rays with shared path-prefixes where the tradeoff is negligible.

The real deal is when you have Path Guiding thrown in as the customary third estimator, now you can reduce this to just two estimators which would be

1. The Product of LTC BSDF and Light Distribution
2. The Product of LTC BSDF and Path Guiding PDF

The second has recently been done by Jakob et al. coming out of EPFL (the creators of Mitsuba 2), and unless I've been reading the Practical Path Guiding paper carelessly the first approach was only ever so slightly hinted at.

Kissing good-bye to the one extra estimator introduced by Path Guiding into the MIS formulation while improving direct light sampling is a non-trivial gain, and sounds like a great opportunity for improvement which I must eventually investigate!

If you manage to get it running before I do, please drop me a message!