Importance Sampling of Environment Maps What is the best currently known and ideally also production-verified approach for sampling environment maps (EM) in a MIS-based uni-directional path tracer and similar types of renderers? I would prefer solutions which are reasonably complicated while reasonably functional to those which provide perfect sampling at cost of super complicated and hard-to-understand implementation.

What I know so far

There are some easy ways of sampling EMs. One can sample the needed hemisphere in a cosine-weighted manner,
   which ignores both the BSDF and the EM function shapes. As a result,
   it doesn't work for dynamic EMs:



To improve
   the sampling to a usable level, one can sample the luminance of the EM
   over the whole sphere. It is relatively easily implemented and the
   results are quite good. However, the sampling strategy is still
   ignoring the hemispherical visibility information and the cosine
   factor (and also the BSDF), resulting in high noise on the surfaces
   which are not directly lit by high-intensity areas of the EM:



Papers

I have found a few papers on the topic, but have not read them yet. Is any of these worth reading and implementing in a forward uni-directional path tracer, or is there something even better?


Structured Importance Sampling of Environment Maps (2003) by
Agarwal et al.
Steerable Importance Sampling (2007) by
Kartic Subr and Jim Arvo. They claim to present “...an algorithm for
efficient stratified importance sampling of environment maps that
generates samples in the positive hemi-sphere defined by local
orientation of arbitrary surfaces while accounting for cosine
weighting. “ The paper “Importance Sampling Spherical Harmonics”
comments on it: ”They create a triangulated representation of the
environment map and store the illumination premultiplied by each of
the first nine spherical harmonic basis functions at every vertex.
This forms a steerable basis where the clamped-cosine can be
efficiently rotated to any orientation.”
Practical Product Importance Sampling for Direct Illumination (2008) by
Petrik Clarberg and Tomas Akenine-Möller. An algorithm for sampling
the product of environment map lighting and surface reflectance. Uses
wavelet-based importance sampling.
Importance Sampling Spherical Harmonics (2009)
by Jarosz, Carr, and Jensenn. The abstract says: “...we present the
first practical method for importance sampling functions represented
as spherical harmonics (SH)...”
Tone-Mapped Mean-shift Based Environment Map Sampling (2015) by Feng
et al. This is pretty new and I have found neither a reference to it
nor the paper itself.
This is not a full answer, I would just like to share the knowledge I obtained by studying two of the papers mentioned in the question: Steerable Importance Sampling and Practical Product Importance Sampling for Direct Illumination.

Steerable Importance Sampling

In this paper they propose a method for sampling the product of the clamped cosine component and environment map lighting:

$$
L_{EM}\left(\omega_{i}\right)\left(\omega_{i}\cdot n\right)^+
$$

They make use of the fact that a piece-wise linear approximation of the product function can be relatively well expressed and partially pre-computed using the first nine spherical harmonic bases. They build this approximation on top of an adaptively triangulated EM and use it as an importance function for sampling.

They pre-compute and store approximation coefficients for each triangle vertex and also coefficients for computation of approximation integral over the triangle for each triangle. These coefficients are called vertex and triangle weights. Then they make use of the fact that is it possible to easily compute coefficients for an integral over a set of triangles just by summing the individual triangle weights without incorporating additional spherical harmonic bases. This allows them to build a balanced binary tree over the triangles where each node contains coefficients for computing approximation integral over the node's sub-tree triangles.

The sampling procedure consists of selecting a triangle and sampling its area:


A triangle is chosen by descending down the pre-built binary tree    with probability proportional to the sub-integral approximations. This costs $O\left(\log N_{\triangle}\right)$ on-the-fly computations of sub-integrals, each consisting of one inner product of clamped-cosine spherical harmonic coordinates with the pre-computed coefficients.   
The chosen triangle surface is then sampled in $O\left(1\right)$ time in a    bi-linear fashion by a novel stratified sampling strategy proposed in    the paper.


To me, this looks like a promising technique, but the classical question with papers is how it will behave in the real life. On the one hand, there may be pathological cases when the EM is hard to approximate with triangulated piece-wise linear function, which can lead to an enormous amount of triangles and/or to poor sample quality. On the other hand, it can instantly provide a relatively good approximation of the whole EM contribution, which can be useful when sampling multiple light sources.

Practical Product Importance Sampling for Direct Illumination

In this paper they propose a method for sampling the product of environment map lighting and cosine-weighted surface reflectance:

$$
L_{EM}\left(\omega_{i}\right)f_{r}\left(\omega_{i},\omega_{o},n\right)\left(\omega_{i}\cdot n\right)^+
$$

The only pre-processing in this method is computation of a hierarchical representation of the EM (either mipmap or wavelet based). The rest is done on the fly during the sampling.

The sampling procedure:


Building an on-the-fly BRDF approximation: They first draw several BRDF importance samples and evaluate $f_{r}\left(\omega_{i},\omega_{o},n\right)\left(\omega_{i}\cdot n\right)^+$. From these values they build a quadtree-based piece-wise constant approximation of the BRDF, where each leaf of the tree contains exactly one sample.
Computing a product of the BRDF approximation and the EM: Multiplication is done at the BRDF quadtree leaves and averaged values are propagated to parents.
Product sampling: uniform samples are fed through the product tree using simple sample warping.


The procedure should generate relatively good samples at the cost of heavy pre-computation &ndash; they show that roughly 100&ndash;200 BRDF samples are needed for BRDF approximation to achieve the best sampling performance. This may make it suitable for purely direct illumination computations, where one generates many samples per shading point, but is most probably too expensive for global illumination algorithms (e.g. uni- or bi-directional path tracers), where you usually generate only a few samples per shading point.