<p>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: <a href="http://home.eps.hw.ac.uk/~ks400/research.html#SteerableImportance" rel="noreferrer">Steerable Importance Sampling</a> and <a href="http://fileadmin.cs.lth.se/graphics/research/papers/2008/pps/" rel="noreferrer">Practical Product Importance Sampling for Direct Illumination</a>.</p>

<h2>Steerable Importance Sampling</h2>

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

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

<p>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.</p>

<p>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.</p>

<p>The sampling procedure consists of selecting a triangle and sampling its area:</p>

<ul>
<li>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.   </li>
<li>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.</li>
</ul>

<p>To me, this looks like a <strong>promising technique</strong>, 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.</p>

<h2>Practical Product Importance Sampling for Direct Illumination</h2>

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

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

<p>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.</p>

<p>The sampling procedure:</p>

<ul>
<li>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.</li>
<li>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.</li>
<li>Product sampling: uniform samples are fed through the product tree using simple sample warping.</li>
</ul>

<p>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 <strong>too expensive for global illumination algorithms</strong> (e.g. uni- or bi-directional path tracers), where you usually generate only a few samples per shading point.</p>

CC BY-SA 3.0
2016-05-21T18:24:35.543Z
2016-05-21T19:29:57.323Z
2016-05-21T19:29:57.323Z
ivokabel
<p>Disclaimer: I have no idea what is the state of the art in the environmental map sampling. In fact, I have very little knowledge about this topic. So this will not be complete answer but I will formulate the problem mathematically and analyze it. I do this mainly for myself, so I make it clear for my self but I hope that OP and others will find it useful.</p>

<hr>

<p>$$
\newcommand{\w}{\omega}
$$</p>

<p>We want to calculate direct illumination at a point i.e. we want to know the value of <a href="https://en.wikipedia.org/wiki/Rendering_equation#Equation_form" rel="nofollow">the integral</a>
$$
I =\int_{S^2} f(\omega_i,\omega_o,n) \, L(\omega_i) \,(\omega_i\cdot n)^+ d\omega_i
$$
where $f(\omega_i,\omega_o,n)$ is BSDF function(I explicitly state dependance on the normal which will be useful later), $L(\omega_i)$ is radiance of environmental map and $(\omega_i \cdot n)^+$ is the cosine term together with the visibility(<a href="https://en.wikipedia.org/wiki/Positive_and_negative_parts" rel="nofollow">that what the $+$ is for</a>) i.e. $(\omega_i \cdot n)^+=0$ if $(\omega_i \cdot n)&lt;0$</p>

<p>We estimate this integral by generating $N$ samples  $\omega_i^1,\dots,\omega_i^N$ with the respect to the probability density function $p(\omega_i)$, the estimator is
$$
I \approx \frac1N \sum_{k=1}^N \frac{f(\omega_i^k,\omega_o,n) \, L(\omega_i^k) \,(\omega_i^k\cdot n)^+ }{p(\omega_i^k)}
$$</p>

<p>The question is: How do we choose the pdf $p$ such that we are able to generate the samples in acceptable time and the variance of the above estimator is reasonably small.</p>

<hr>

<p><strong><s>Best</s> method</strong> Pick $p$ proportional to the integrand
$$
p(\omega_i) \sim f(\omega_i,\omega_o,n) \, L(\omega_i) \,(\omega_i\cdot n)^+ 
$$
But most of the times it is very expensive to generate a sample according to this pdf, so it is not useful in practice.</p>

<p>Methods suggested by OP:</p>

<p><strong>Method one</strong>: Choose $p$ proportional to the cosine term
$$
p(\omega_i) \sim (\omega_i\cdot n)^+
$$
<strong>Method two</strong>: Choose $p$ proportional to the EM
$$
p(\omega_i) \sim L(\omega_i)
$$</p>

<hr>

<p>Based on the names of mentioned papers I can partially guess what they do(unfortunately I do not have the time and energy to read them right now). But before discussing what they most probably do, let's talk about power series a little bit :D</p>

<hr>

<p>If we have a function of one real variable e.g. $f(x)$. Then if it is well behaved then it can be expanded into a power series
$$
f(x) = \sum_{k=0}^\infty a_k x^k
$$
Where $a_k$ are constants. This can be used to approximate $f$ by truncating the sum at some step $n$
$$
f(x) \approx \sum_{k=0}^n a_k x^k
$$
If $n$ is sufficiently high then the error is really small.</p>

<p>Now if we have function in two variables e.g. $f(x,y)$ we can expand it only in the first argument
$$
f(x,y) = \sum_{k=0}^\infty b_k(y) \, x^k
$$
where  $b_k(y)$ are functions only in $y$. It can be also expanded in both arguments
$$
f(x,y) = \sum_{k,l=0}^\infty c_{kl} x^k y^l
$$
where $c_{kl}$ are constants. So function with real arguments can be expanded as sum of powers of that argument. Something similar can be done for functions defined on sphere.</p>

<p>Now, let's have a function defined on sphere e.g. $f(\omega)$. Such a function can be also expanded in similar fashion as function of one real parameter
$$
f(\omega) =\sum_{k=0}^\infty \alpha_k S_k(\omega)
$$
where $\alpha_k$ are constants and $S_k(\omega)$ are <a href="https://en.wikipedia.org/wiki/Spherical_harmonics" rel="nofollow">spherical harmonics</a>. Spherical harmonics are normally indexed by two indices and are written as function in spherical coordinates but that is not important here. The important thing is that $f$ can be written as a sum of some known functions.</p>

<p>Now function which takes two points on sphere e.g. $f(\omega,\omega')$ can be expanded only in its first arguments
$$
f(\omega,\omega') = \sum_{k=0}^\infty \beta_k(\omega') \, S_k(\omega)
$$
or in both its arguments
$$
f(\omega,\omega') = \sum_{k,l=0}^\infty \gamma_{kl} \, S_k(\omega)S_l(\omega')
$$</p>

<hr>

<p>So how is this all useful?</p>

<p>I propose the  <strong>CMUNSM</strong>(Crazy mental useless no sampling method):
Lets assume that we have expansions for all the function i.e.
\begin{align}
f(\omega_i,\omega_o,n) &amp;= \sum_{k,l,m=0}^\infty \alpha_{klm} S_k(\omega_i)S_l(\omega_o) S_m(n) \\
L(\omega_i ) &amp;= \sum_{n=0}^\infty \beta_n S_n(\omega) \\
(\omega_i\cdot n)^+ &amp;= \sum_{p,q=0}^\infty \gamma_{pq} S_p(\omega_i)S_q(n)
\end{align}
If we plug this into the integral we get
$$
I = \sum_{k,l,m,n,p,q=0}^\infty \alpha_{klm} \beta_n  \gamma_{pq} S_l(\omega_o) S_m(n) S_q(n) \int_{S^2} S_k(\omega_i) S_n(\omega)  S_p(\omega_i) d\omega_i
$$</p>

<p>Actually we no longer need Monte Carlo because we can calculate values of the integrals $\int_{S^2} S_k(\omega_i) S_n(\omega)  S_p(\omega_i) d\omega_i$ beforehand and then evaluate the sum(actually approximate the sum, we would sum only first few terms) and we get desired result.</p>

<p>This is all nice <strong>but</strong> we might not know the expansions of BSDF or environmental map or the expansions converge very slowly therefore we would have to take a lots of terms in the sum to get reasonably accurate answer.</p>

<hr>

<p>So the idea is not to expand in all arguments. One method which might be worth investigating would be to ignore BSDF and expand only the environmental map i.e.
$$
L(\omega_i) \approx \sum_{n=0}^K \beta_n S_n(\omega_i)
$$
this would lead to pdf:
$$
p(\omega_i) \sim \sum_{n=0}^K \beta_n S_n(\omega_i) (\omega \cdot n)^+
$$</p>

<p>We already know how to do this for $K=0$, this is nothing but the <strong>method one</strong>. My guess is, it is done in one of the papers for higher $K$.</p>

<hr>

<p>Further extensions. You can expand different functions in different arguments and do similar stuff as above. Another thing is, that you can expand in different basis, i.e. do not use spherical harmonics but different functions.</p>

<p>So this is my take on the topic, I hope you have found it at least a little bit useful and now I'm off to GoT and bed.</p>

CC BY-SA 3.0
2016-05-16T20:33:39.170Z
Haha, when I posted the answer, SE asked me if I'm a human or a robot, the site wasn't sure :D I hope it is not because of the length of the answer, It got a little bit out of hand.
tom
CC BY-SA 3.0
2016-05-16T20:39:10.357Z
you want to make my brain melt, don't you. ;-) BTW: I already managed to read two of the papers/presentations so I'll hopefully extend the question or write a superficial answer at the end of this week. And now, GoT FTW!
ivokabel
CC BY-SA 3.0
2016-05-16T20:30:13.717Z
2016-05-16T20:30:13.717Z
tom
<p>While the product sampling methods provides better (perfect) distribution for rays I would say that using MIS (multiple importance sampling) is a method verified in production. Since shadowing information is unknown product sampling doesn't become perfect anyway and it is quite hard to implmenet. Shooting more rays might be worth more! Depends on your situation and ray budgets of course!</p>

<p>Short description of MIS: In essence you trace both a BSDF-ray (as you would anyway for doing indirect lighting) and an explicit ray towards the EM. MIS give you weights so that you can combine them in a way that removes a lot of the noise. MIS is especially good at choosing "technique" (implicit or explicit sampling) based on the situation that arises. This happens naturally without the user having to make hard choices based on roughness etc.</p>

<p>Chapter 9 of <a href="http://graphics.stanford.edu/papers/veach_thesis/" rel="nofollow">http://graphics.stanford.edu/papers/veach_thesis/</a> covers this in detail. Also see <a href="https://www.shadertoy.com/view/4sSXWt" rel="nofollow">https://www.shadertoy.com/view/4sSXWt</a> for a demo of MIS in action with area lights.</p>

CC BY-SA 3.0
2016-05-14T12:28:28.730Z
Yes, MIS is an important production-verified technique, which helps a lot and I employ it in my solution (I guess, I should have stated that more clearly in the question).

However, the overall performance of a MIS-based estimator depends on the quality of its partial sampling strategies. What I am trying to do here is to improve one of the the sub-strategies to improve the overall performance of the estimator. In my experience, it is usually more efficient to use less high-quality samples than may be more expensive to generate than more easily-generated low-quality ones.
ivokabel
CC BY-SA 3.0
2016-05-13T22:34:31.890Z
2016-05-13T22:34:31.890Z
anders
<p>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.</p>

<h2>What I know so far</h2>

<p>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:</p>

<p><a href="https://i.sstatic.net/TBYYn.png" rel="noreferrer"><img src="https://i.sstatic.net/TBYYn.png" alt="enter image description here"></a></p>

<p>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:</p>

<p><a href="https://i.sstatic.net/1DOSr.png" rel="noreferrer"><img src="https://i.sstatic.net/1DOSr.png" alt="enter image description here"></a></p>

<h2>Papers</h2>

<p>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?</p>

<ul>
<li><p><a href="https://cseweb.ucsd.edu/~ravir/papers/structured/" rel="noreferrer">Structured Importance Sampling of Environment Maps</a> (2003) by
Agarwal et al.</p></li>
<li><p><a href="http://home.eps.hw.ac.uk/~ks400/research.html#SteerableImportance" rel="noreferrer">Steerable Importance Sampling</a> (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.”</p></li>
<li><p><a href="http://fileadmin.cs.lth.se/graphics/research/papers/2008/pps/" rel="noreferrer">Practical Product Importance Sampling for Direct Illumination</a> (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.</p></li>
<li><p><a href="https://www.cs.dartmouth.edu/~wjarosz/publications/jarosz09importance.html" rel="noreferrer">Importance Sampling Spherical Harmonics</a> (2009)
by Jarosz, Carr, and Jensenn. The abstract says: “...we present the
first practical method for importance sampling functions represented
as spherical harmonics (SH)...”</p></li>
<li><p><a href="http://www.ncbi.nlm.nih.gov/pubmed/26584494" rel="noreferrer">Tone-Mapped Mean-shift Based Environment Map Sampling</a> (2015) by Feng
et al. This is pretty new and I have found neither a reference to it
nor the paper itself.</p></li>
</ul>

CC BY-SA 3.0
2016-05-17T20:42:02.350Z
I have one question. Is the second picture generated only by sampling the EM? Or is it MISed version of sampling cosine and sampling EM? I really hope that it is the MISed version, because if so, then I might have a remedy for the high noise in the shadowy part.
tom
CC BY-SA 3.0
2016-05-17T22:17:39.087Z
No @tom, it uses sperical EM sampling only, ignoring both the (Lambert) BRDF and the cosine factor. 64 samples were used and no image-space filtering applied, just averaging over pixel area. When MIS is applied to combine the EM sampling with the cosine sampling, the noise in the shadow decreases a lot, but slightly increases in the sunlit part.
ivokabel
<rendering><sampling><pathtracing><importance-sampling><ibl>
Importance Sampling of Environment Maps
ivokabel
# Importance Sampling of Environment Maps

(Asked by: username_0 on 2016-05-04)

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:

[![enter image description here](https://i.sstatic.net/TBYYn.png)](https://i.sstatic.net/TBYYn.png)

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:

[![enter image description here](https://i.sstatic.net/1DOSr.png)](https://i.sstatic.net/1DOSr.png)

## 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](https://cseweb.ucsd.edu/~ravir/papers/structured/) (2003) by
  Agarwal et al.
* [Steerable Importance Sampling](http://home.eps.hw.ac.uk/~ks400/research.html#SteerableImportance) (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](http://fileadmin.cs.lth.se/graphics/research/papers/2008/pps/) (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](https://www.cs.dartmouth.edu/~wjarosz/publications/jarosz09importance.html) (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](http://www.ncbi.nlm.nih.gov/pubmed/26584494) (2015) by Feng
  et al. This is pretty new and I have found neither a reference to it
  nor the paper itself.

**Tags:** rendering, sampling, pathtracing, importance-sampling, ibl

> **Comments:**
> **username_1**: I have one question. Is the second picture generated only by sampling the EM? Or is it MISed version of sampling cosine and sampling EM? I really hope that it is the MISed version, because if so, then I might have a remedy for the high noise in the shadowy part.  
> **username_0**: No @username_1, it uses sperical EM sampling only, ignoring both the (Lambert) BRDF and the cosine factor. 64 samples were used and no image-space filtering applied, just averaging over pixel area. When MIS is applied to combine the EM sampling with the cosine sampling, the noise in the shadow decreases a lot, but slightly increases in the sunlit part.  

---

## 3 answers

### Answer by username_0 on 2016-05-21. Score: 7

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](http://home.eps.hw.ac.uk/~ks400/research.html#SteerableImportance) and [Practical Product Importance Sampling for Direct Illumination](http://fileadmin.cs.lth.se/graphics/research/papers/2008/pps/).

## 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 – they show that roughly 100–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.


---

### Answer by username_1 on 2016-05-16. Score: 4

Disclaimer: I have no idea what is the state of the art in the environmental map sampling. In fact, I have very little knowledge about this topic. So this will not be complete answer but I will formulate the problem mathematically and analyze it. I do this mainly for myself, so I make it clear for my self but I hope that OP and others will find it useful.

---

$$
\newcommand{\w}{\omega}
$$

We want to calculate direct illumination at a point i.e. we want to know the value of [the integral](https://en.wikipedia.org/wiki/Rendering_equation#Equation_form)
$$
I =\int\_{S^2} f(\omega\_i,\omega\_o,n) \, L(\omega\_i) \,(\omega\_i\cdot n)^+ d\omega\_i
$$
where $f(\omega\_i,\omega\_o,n)$ is BSDF function(I explicitly state dependance on the normal which will be useful later), $L(\omega\_i)$ is radiance of environmental map and $(\omega\_i \cdot n)^+$ is the cosine term together with the visibility([that what the $+$ is for](https://en.wikipedia.org/wiki/Positive_and_negative_parts)) i.e. $(\omega\_i \cdot n)^+=0$ if $(\omega\_i \cdot n)<0$

We estimate this integral by generating $N$ samples $\omega\_i^1,\dots,\omega\_i^N$ with the respect to the probability density function $p(\omega\_i)$, the estimator is
$$
I \approx \frac1N \sum\_{k=1}^N \frac{f(\omega\_i^k,\omega\_o,n) \, L(\omega\_i^k) \,(\omega\_i^k\cdot n)^+ }{p(\omega\_i^k)}
$$

The question is: How do we choose the pdf $p$ such that we are able to generate the samples in acceptable time and the variance of the above estimator is reasonably small.

---

**~~Best~~ method** Pick $p$ proportional to the integrand
$$
p(\omega\_i) \sim f(\omega\_i,\omega\_o,n) \, L(\omega\_i) \,(\omega\_i\cdot n)^+
$$
But most of the times it is very expensive to generate a sample according to this pdf, so it is not useful in practice.

Methods suggested by OP:

**Method one**: Choose $p$ proportional to the cosine term
$$
p(\omega\_i) \sim (\omega\_i\cdot n)^+
$$
**Method two**: Choose $p$ proportional to the EM
$$
p(\omega\_i) \sim L(\omega\_i)
$$

---

Based on the names of mentioned papers I can partially guess what they do(unfortunately I do not have the time and energy to read them right now). But before discussing what they most probably do, let's talk about power series a little bit :D

---

If we have a function of one real variable e.g. $f(x)$. Then if it is well behaved then it can be expanded into a power series
$$
f(x) = \sum\_{k=0}^\infty a\_k x^k
$$
Where $a\_k$ are constants. This can be used to approximate $f$ by truncating the sum at some step $n$
$$
f(x) \approx \sum\_{k=0}^n a\_k x^k
$$
If $n$ is sufficiently high then the error is really small.

Now if we have function in two variables e.g. $f(x,y)$ we can expand it only in the first argument
$$
f(x,y) = \sum\_{k=0}^\infty b\_k(y) \, x^k
$$
where $b\_k(y)$ are functions only in $y$. It can be also expanded in both arguments
$$
f(x,y) = \sum\_{k,l=0}^\infty c\_{kl} x^k y^l
$$
where $c\_{kl}$ are constants. So function with real arguments can be expanded as sum of powers of that argument. Something similar can be done for functions defined on sphere.

Now, let's have a function defined on sphere e.g. $f(\omega)$. Such a function can be also expanded in similar fashion as function of one real parameter
$$
f(\omega) =\sum\_{k=0}^\infty \alpha\_k S\_k(\omega)
$$
where $\alpha\_k$ are constants and $S\_k(\omega)$ are [spherical harmonics](https://en.wikipedia.org/wiki/Spherical_harmonics). Spherical harmonics are normally indexed by two indices and are written as function in spherical coordinates but that is not important here. The important thing is that $f$ can be written as a sum of some known functions.

Now function which takes two points on sphere e.g. $f(\omega,\omega')$ can be expanded only in its first arguments
$$
f(\omega,\omega') = \sum\_{k=0}^\infty \beta\_k(\omega') \, S\_k(\omega)
$$
or in both its arguments
$$
f(\omega,\omega') = \sum\_{k,l=0}^\infty \gamma\_{kl} \, S\_k(\omega)S\_l(\omega')
$$

---

So how is this all useful?

I propose the **CMUNSM**(Crazy mental useless no sampling method):
Lets assume that we have expansions for all the function i.e.
\begin{align}
f(\omega\_i,\omega\_o,n) &= \sum\_{k,l,m=0}^\infty \alpha\_{klm} S\_k(\omega\_i)S\_l(\omega\_o) S\_m(n) \\
L(\omega\_i ) &= \sum\_{n=0}^\infty \beta\_n S\_n(\omega) \\
(\omega\_i\cdot n)^+ &= \sum\_{p,q=0}^\infty \gamma\_{pq} S\_p(\omega\_i)S\_q(n)
\end{align}
If we plug this into the integral we get
$$
I = \sum\_{k,l,m,n,p,q=0}^\infty \alpha\_{klm} \beta\_n \gamma\_{pq} S\_l(\omega\_o) S\_m(n) S\_q(n) \int\_{S^2} S\_k(\omega\_i) S\_n(\omega) S\_p(\omega\_i) d\omega\_i
$$

Actually we no longer need Monte Carlo because we can calculate values of the integrals $\int\_{S^2} S\_k(\omega\_i) S\_n(\omega) S\_p(\omega\_i) d\omega\_i$ beforehand and then evaluate the sum(actually approximate the sum, we would sum only first few terms) and we get desired result.

This is all nice **but** we might not know the expansions of BSDF or environmental map or the expansions converge very slowly therefore we would have to take a lots of terms in the sum to get reasonably accurate answer.

---

So the idea is not to expand in all arguments. One method which might be worth investigating would be to ignore BSDF and expand only the environmental map i.e.
$$
L(\omega\_i) \approx \sum\_{n=0}^K \beta\_n S\_n(\omega\_i)
$$
this would lead to pdf:
$$
p(\omega\_i) \sim \sum\_{n=0}^K \beta\_n S\_n(\omega\_i) (\omega \cdot n)^+
$$

We already know how to do this for $K=0$, this is nothing but the **method one**. My guess is, it is done in one of the papers for higher $K$.

---

Further extensions. You can expand different functions in different arguments and do similar stuff as above. Another thing is, that you can expand in different basis, i.e. do not use spherical harmonics but different functions.

So this is my take on the topic, I hope you have found it at least a little bit useful and now I'm off to GoT and bed.

> **Comments:**
> **username_1**: Haha, when I posted the answer, SE asked me if I'm a human or a robot, the site wasn't sure :D I hope it is not because of the length of the answer, It got a little bit out of hand.  
> **username_0**: you want to make my brain melt, don't you. ;-) BTW: I already managed to read two of the papers/presentations so I'll hopefully extend the question or write a superficial answer at the end of this week. And now, GoT FTW!  

---

### Answer by username_2 on 2016-05-13. Score: 0

While the product sampling methods provides better (perfect) distribution for rays I would say that using MIS (multiple importance sampling) is a method verified in production. Since shadowing information is unknown product sampling doesn't become perfect anyway and it is quite hard to implmenet. Shooting more rays might be worth more! Depends on your situation and ray budgets of course!

Short description of MIS: In essence you trace both a BSDF-ray (as you would anyway for doing indirect lighting) and an explicit ray towards the EM. MIS give you weights so that you can combine them in a way that removes a lot of the noise. MIS is especially good at choosing "technique" (implicit or explicit sampling) based on the situation that arises. This happens naturally without the user having to make hard choices based on roughness etc.

Chapter 9 of <http://graphics.stanford.edu/papers/veach_thesis/> covers this in detail. Also see <https://www.shadertoy.com/view/4sSXWt> for a demo of MIS in action with area lights.

> **Comments:**
> **username_0**: Yes, MIS is an important production-verified technique, which helps a lot and I employ it in my solution (I guess, I should have stated that more clearly in the question).  However, the overall performance of a MIS-based estimator depends on the quality of its partial sampling strategies. What I am trying to do here is to improve one of the the sub-strategies to improve the overall performance of the estimator. In my experience, it is usually more efficient to use less high-quality samples than may be more expensive to generate than more easily-generated low-quality ones.  

---

