<p>For sampling a uniform volume you use the mean free path of a photon:</p>
<p><code>float dt = -logf(1.0f - Xi) / uT;</code></p>
<p>where:</p>
<ul>
<li>Xi is a random variable in [0, 1]</li>
<li>uT is the extinction coefficient (sum of out-scattering and absorption)</li>
</ul>
<p>The pdf of a raytime (t) in this uniform volume is given by:</p>
<p><code>float pdf = uT * expf(-t * uT);</code></p>
<p>In your case, you do not have a uniform volume, but instead one whose density is exponentially decreasing according to the distance from the surface of a sphere.</p>
<p>Approaches for solving this typically involve raymarching which samples a density at each step and treats the step as a uniform volume. This is approximate and introduces some bias, but there are ways of reducing the bias.</p>
<p>I recommend reading Monte Carlo Methods for Physically Based Volume Rendering:</p>
<ul>
<li><a href="https://jannovak.info/publications/VolumeCourse/index.html" rel="nofollow noreferrer">https://jannovak.info/publications/VolumeCourse/index.html</a></li>
</ul>
<p>Another good resource is pbrtv3, particularly chapters 11 and 15:</p>
<ul>
<li><a href="http://www.pbr-book.org/3ed-2018/Volume_Scattering.html" rel="nofollow noreferrer">http://www.pbr-book.org/3ed-2018/Volume_Scattering.html</a></li>
<li><a href="http://www.pbr-book.org/3ed-2018/Light_Transport_II_Volume_Rendering.html" rel="nofollow noreferrer">http://www.pbr-book.org/3ed-2018/Light_Transport_II_Volume_Rendering.html</a></li>
</ul>

CC BY-SA 4.0
2020-10-17T19:40:33.637Z
2020-10-17T19:40:33.637Z
Kara
<p><strong>Update</strong></p>

<p>I am editing and posting this question in a different way; this time from the point of view of <a href="http://nishitalab.org/user/nis/cdrom/sig93_nis.pdf" rel="nofollow noreferrer">Nishita paper</a>.</p>

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

<p>Sunlight gets scattered at P and attenuated before reaching Pv. Therefore intensity reaching Pv coming from P can be obtained by multiplying the attenuation by intensity at P:</p>

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

<p>The intensity reaching Pv therefore can be obtained by integrating light due to air molecules along the ray segment PaPb:</p>

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

<p>Where </p>

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

<p>If I am to solve this integral by Monte Carlo method and without ray marching approach, What would be the estimator and what are the good choices for pdf and sampling method?</p>

<p>Thanks.</p>

<hr>

<p>I am trying to implement the concept of distributed ray tracing on Rayleigh scatterting to render sky a bit more realistic. The code I use is borrowed from
<a href="https://www.scratchapixel.com/lessons/procedural-generation-virtual-worlds/simulating-sky" rel="nofollow noreferrer">this website</a>. It takes 16 samples along the ray and 8 samples for sunlight.
This works fine and I was able to generate this image using only 1 sample per pixel(not sure why the sun is so small, any idea?):</p>

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

<p>I then tried to use Monte Carlo to do this. 1 sample uniformly selected along the ray and 1 sample for sunlight. The only changes I did are on four lines:</p>

<pre><code>line 7:  uint32_t numSamples = 1; //16
line 8:  uint32_t numSamplesLight = 1; //8 

//used drand48() instead of 0.5
line 18: Vec3f samplePosition = orig + (tCurrent+
                                segmentLength * drand48()) * dir; 
line 32: Vec3f samplePositionLight = samplePosition + (tCurrentLight + 
                                segmentLengthLight * drand48) * sunDirection; 
</code></pre>

<p>The pdf is 1/segmentLength and the sample is taken uniformly in this distance.
Here is my implementation in c# code:</p>

<pre><code>public override void Render(int pxlX, int pxlY, Random rand)
{
    Vector estimatedColor = new Vector();

    Ray ray = Scene.Camera.GenerateRay(pxlX + rand.NextDouble(), pxlY + rand.NextDouble(), rand.NextDouble(), rand.NextDouble());

    ISect isect = Scene.Trace(ray);
    if (ISect.IsNull(isect))
        estimatedColor = Atmosph.ComputeIncidentLight(ray, rand);
    else
        estimatedColor = isect.Thing.Material.Emission +
                ComputeRadiance(isect, rand, 0);

    ImageFilm.AddToPixel(pxlX, pxlY, estimatedColor/numPixelSamples);
}

public Vector ComputeIncidentLight(Ray r, Random rand)
{
    r.Start.Z += earthRadius;

    ISect isect = earth.Intersect(r);
    if (ISect.IsNull(isect)) return Vector.ZERO;

    int numSamples = 1;
    int numSamplesLight = 1;
    double segmentLength = isect.Dist / numSamples;
    double tCurrent = 0;
    Vector sumR = Vector.ZERO;
    Vector sumM = Vector.ZERO;// mie and rayleigh contribution 
    double opticalDepthR = 0, opticalDepthM = 0;
    double mu = r.Dir.Dot(sunDirection); // mu in the paper which is the cosine of the angle between the sun direction and the ray direction 
    double phaseR = 3d / (16d * Math.PI) * (1 + mu * mu);
    double g = 0.76f;
    double phaseM = 3d / (8d * Math.PI) * ((1d - g * g) * (1d + mu * mu)) / ((2d + g * g) * Math.Pow(1d + g * g - 2d * g * mu, 1.5d));
    for (int i = 0; i &lt; numSamples; ++i)
    {
        Vector samplePosition = r.Start + r.Dir * (tCurrent + segmentLength * rand.NextDouble());
        double height = samplePosition.Length() - earthRadius;
        // compute optical depth for light
        double hr = Math.Exp(-height / Hr) * segmentLength;
        double hm = Math.Exp(-height / Hm) * segmentLength;
        opticalDepthR += hr;
        opticalDepthM += hm;
        // light optical depth
        ISect isectSun = earth.Intersect(new Ray(samplePosition, sunDirection));
        double segmentLengthLight = isectSun.Dist/ numSamplesLight, tCurrentLight = 0;
        double opticalDepthLightR = 0, opticalDepthLightM = 0;
        int j;
        for (j = 0; j &lt; numSamplesLight; ++j)
        {
            Vector samplePositionLight = samplePosition + sunDirection * (tCurrentLight + segmentLengthLight * rand.NextDouble());
            double heightLight = samplePositionLight.Length() - earthRadius;
            if (heightLight &lt; 0) break;
            opticalDepthLightR += Math.Exp(-heightLight / Hr) * segmentLengthLight;
            opticalDepthLightM += Math.Exp(-heightLight / Hm) * segmentLengthLight;
            tCurrentLight += segmentLengthLight;
        }
        if (j == numSamplesLight)
        {
            Vector tau = betaR * (opticalDepthR + opticalDepthLightR) + betaM * 1.1f * (opticalDepthM + opticalDepthLightM);
            Vector attenuation = new Vector(Math.Exp(-tau.X), Math.Exp(-tau.Y), Math.Exp(-tau.Z));
            sumR += attenuation * hr;
            sumM += attenuation * hm;
        }
        tCurrent += segmentLength;
    }

    // We use a magic number here for the intensity of the sun (20). We will make it more
    // scientific in a future revision of this lesson/code
    return (sumR * betaR * phaseR + sumM * betaM * phaseM) * sunIrradiance;
}
</code></pre>

<p>The below image rendered with 400 samples per pixel, however is different from the above:</p>

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

<p>Could anyone shed some light where I've got it wrong?
Thanks,</p>

CC BY-SA 3.0
2017-08-23T07:22:37.037Z
Without reading the code and just from the picture, my intuition would be to double check the integration. It looks like too much light is absorbed along the ray. I would first check that in the first implementation the color stays the same when changing the number of steps (just to make sure it's a good reference image), then I would double check the absorption equation used in the MT implementation.
Julien Guertault
CC BY-SA 3.0
2017-08-26T09:21:48.233Z
Thanks for your comment. Still not working and I do not understand why. I sample a point uniformly along the ray, instead of dividing it to segment, and the pdf is simply 1/dist. But for some reason it doesn't converge to same result.
ali
CC BY-SA 3.0
2017-09-01T00:29:41.157Z
The comment was directed toward people willing to propose answers. But thanks, it makes the question more interesting too. :)
Julien Guertault
CC BY-SA 3.0
2018-03-05T08:06:14.887Z
I would be suspicious of this line> Vector samplePositionLight = samplePosition + sunDirection * (tCurrentLight + segmentLengthLight * rand.NextDouble());     >> This should walk the sample-to-light in even steps, why is it doing a random length walk ?
PaulHK
CC BY-SA 3.0
2018-03-05T08:55:51.550Z
I've implemented this algorithm before and also found it unnecessary to generate multiple rays for each 'pixel'. A single ray, which raymarches from the view direction, and for each of those steps marches towards the light was enough to get high quality results. Have you tried just one call to ComputeIncidentLight per pixel? (Also see my previous post, you have modified the sample->light ray march in a way that may have broken it)
PaulHK
CC BY-SA 3.0
2018-03-05T12:16:17.097Z
@PaulHK, well even with one sample ppx you still need 16 x 8 samples for marching algorithm. You may think this isn't costly but if you do, then you may use Monte Carlo method instead of marching the along the ray in steps.
ali
CC BY-SA 3.0
2018-03-06T02:08:50.527Z
But stepping by rand.NextDouble() in numLightSamples steps means you don't end up at isectSun.Dist as the original code does, this i think is the cause of your problem.
PaulHK
CC BY-SA 3.0
2018-03-06T03:56:59.813Z
Yes 16x8 samples/pix is heavy, you can reduce this a little and get away with good results (usually lower samples makes the atmosphere shadow area more coarse). You could also look into rendering the atmosphere at a lower resolution and upscale
PaulHK
<pathtracing><global-illumination><distribution><volumetric-scattering>
Using Monte carlo on Rayleigh scattering
ali
# Using Monte carlo on Rayleigh scattering

(Asked by: username_0 on 2017-08-22)

**Update**

I am editing and posting this question in a different way; this time from the point of view of [Nishita paper](http://nishitalab.org/user/nis/cdrom/sig93_nis.pdf).

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

Sunlight gets scattered at P and attenuated before reaching Pv. Therefore intensity reaching Pv coming from P can be obtained by multiplying the attenuation by intensity at P:

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

The intensity reaching Pv therefore can be obtained by integrating light due to air molecules along the ray segment PaPb:

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

Where

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

If I am to solve this integral by Monte Carlo method and without ray marching approach, What would be the estimator and what are the good choices for pdf and sampling method?

Thanks.

---

I am trying to implement the concept of distributed ray tracing on Rayleigh scatterting to render sky a bit more realistic. The code I use is borrowed from
[this website](https://www.scratchapixel.com/lessons/procedural-generation-virtual-worlds/simulating-sky). It takes 16 samples along the ray and 8 samples for sunlight.
This works fine and I was able to generate this image using only 1 sample per pixel(not sure why the sun is so small, any idea?):

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

I then tried to use Monte Carlo to do this. 1 sample uniformly selected along the ray and 1 sample for sunlight. The only changes I did are on four lines:

```
line 7:  uint32_t numSamples = 1; //16
line 8:  uint32_t numSamplesLight = 1; //8 

//used drand48() instead of 0.5
line 18: Vec3f samplePosition = orig + (tCurrent+
                                segmentLength * drand48()) * dir; 
line 32: Vec3f samplePositionLight = samplePosition + (tCurrentLight + 
                                segmentLengthLight * drand48) * sunDirection; 

```

The pdf is 1/segmentLength and the sample is taken uniformly in this distance.
Here is my implementation in c# code:

```
public override void Render(int pxlX, int pxlY, Random rand)
{
    Vector estimatedColor = new Vector();

    Ray ray = Scene.Camera.GenerateRay(pxlX + rand.NextDouble(), pxlY + rand.NextDouble(), rand.NextDouble(), rand.NextDouble());

    ISect isect = Scene.Trace(ray);
    if (ISect.IsNull(isect))
        estimatedColor = Atmosph.ComputeIncidentLight(ray, rand);
    else
        estimatedColor = isect.Thing.Material.Emission +
                ComputeRadiance(isect, rand, 0);

    ImageFilm.AddToPixel(pxlX, pxlY, estimatedColor/numPixelSamples);
}

public Vector ComputeIncidentLight(Ray r, Random rand)
{
    r.Start.Z += earthRadius;

    ISect isect = earth.Intersect(r);
    if (ISect.IsNull(isect)) return Vector.ZERO;

    int numSamples = 1;
    int numSamplesLight = 1;
    double segmentLength = isect.Dist / numSamples;
    double tCurrent = 0;
    Vector sumR = Vector.ZERO;
    Vector sumM = Vector.ZERO;// mie and rayleigh contribution 
    double opticalDepthR = 0, opticalDepthM = 0;
    double mu = r.Dir.Dot(sunDirection); // mu in the paper which is the cosine of the angle between the sun direction and the ray direction 
    double phaseR = 3d / (16d * Math.PI) * (1 + mu * mu);
    double g = 0.76f;
    double phaseM = 3d / (8d * Math.PI) * ((1d - g * g) * (1d + mu * mu)) / ((2d + g * g) * Math.Pow(1d + g * g - 2d * g * mu, 1.5d));
    for (int i = 0; i < numSamples; ++i)
    {
        Vector samplePosition = r.Start + r.Dir * (tCurrent + segmentLength * rand.NextDouble());
        double height = samplePosition.Length() - earthRadius;
        // compute optical depth for light
        double hr = Math.Exp(-height / Hr) * segmentLength;
        double hm = Math.Exp(-height / Hm) * segmentLength;
        opticalDepthR += hr;
        opticalDepthM += hm;
        // light optical depth
        ISect isectSun = earth.Intersect(new Ray(samplePosition, sunDirection));
        double segmentLengthLight = isectSun.Dist/ numSamplesLight, tCurrentLight = 0;
        double opticalDepthLightR = 0, opticalDepthLightM = 0;
        int j;
        for (j = 0; j < numSamplesLight; ++j)
        {
            Vector samplePositionLight = samplePosition + sunDirection * (tCurrentLight + segmentLengthLight * rand.NextDouble());
            double heightLight = samplePositionLight.Length() - earthRadius;
            if (heightLight < 0) break;
            opticalDepthLightR += Math.Exp(-heightLight / Hr) * segmentLengthLight;
            opticalDepthLightM += Math.Exp(-heightLight / Hm) * segmentLengthLight;
            tCurrentLight += segmentLengthLight;
        }
        if (j == numSamplesLight)
        {
            Vector tau = betaR * (opticalDepthR + opticalDepthLightR) + betaM * 1.1f * (opticalDepthM + opticalDepthLightM);
            Vector attenuation = new Vector(Math.Exp(-tau.X), Math.Exp(-tau.Y), Math.Exp(-tau.Z));
            sumR += attenuation * hr;
            sumM += attenuation * hm;
        }
        tCurrent += segmentLength;
    }

    // We use a magic number here for the intensity of the sun (20). We will make it more
    // scientific in a future revision of this lesson/code
    return (sumR * betaR * phaseR + sumM * betaM * phaseM) * sunIrradiance;
}

```

The below image rendered with 400 samples per pixel, however is different from the above:

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

Could anyone shed some light where I've got it wrong?
Thanks,

**Tags:** pathtracing, global-illumination, distribution, volumetric-scattering

> **Comments:**
> **username_1**: Without reading the code and just from the picture, my intuition would be to double check the integration. It looks like too much light is absorbed along the ray. I would first check that in the first implementation the color stays the same when changing the number of steps (just to make sure it's a good reference image), then I would double check the absorption equation used in the MT implementation.  
> **username_0**: Thanks for your comment. Still not working and I do not understand why. I sample a point uniformly along the ray, instead of dividing it to segment, and the pdf is simply 1/dist. But for some reason it doesn't converge to same result.  
> **username_1**: The comment was directed toward people willing to propose answers. But thanks, it makes the question more interesting too. :)  
> **username_2**: I would be suspicious of this line> Vector samplePositionLight = samplePosition + sunDirection * (tCurrentLight + segmentLengthLight * rand.NextDouble());     >> This should walk the sample-to-light in even steps, why is it doing a random length walk ?  
> **username_2**: I've implemented this algorithm before and also found it unnecessary to generate multiple rays for each 'pixel'. A single ray, which raymarches from the view direction, and for each of those steps marches towards the light was enough to get high quality results. Have you tried just one call to ComputeIncidentLight per pixel? (Also see my previous post, you have modified the sample->light ray march in a way that may have broken it)  
> **username_0**: @username_2, well even with one sample ppx you still need 16 x 8 samples for marching algorithm. You may think this isn't costly but if you do, then you may use Monte Carlo method instead of marching the along the ray in steps.  
> **username_2**: But stepping by rand.NextDouble() in numLightSamples steps means you don't end up at isectSun.Dist as the original code does, this i think is the cause of your problem.  
> **username_2**: Yes 16x8 samples/pix is heavy, you can reduce this a little and get away with good results (usually lower samples makes the atmosphere shadow area more coarse). You could also look into rendering the atmosphere at a lower resolution and upscale  

---

## 1 answers

### Answer by username_3 on 2020-10-17. Score: 1

For sampling a uniform volume you use the mean free path of a photon:

`float dt = -logf(1.0f - Xi) / uT;`

where:

* Xi is a random variable in [0, 1]
* uT is the extinction coefficient (sum of out-scattering and absorption)

The pdf of a raytime (t) in this uniform volume is given by:

`float pdf = uT * expf(-t * uT);`

In your case, you do not have a uniform volume, but instead one whose density is exponentially decreasing according to the distance from the surface of a sphere.

Approaches for solving this typically involve raymarching which samples a density at each step and treats the step as a uniform volume. This is approximate and introduces some bias, but there are ways of reducing the bias.

I recommend reading Monte Carlo Methods for Physically Based Volume Rendering:

* <https://jannovak.info/publications/VolumeCourse/index.html>

Another good resource is pbrtv3, particularly chapters 11 and 15:

* <http://www.pbr-book.org/3ed-2018/Volume_Scattering.html>
* <http://www.pbr-book.org/3ed-2018/Light_Transport_II_Volume_Rendering.html>


---

