Using Monte carlo on Rayleigh scattering

Question

Update

I am editing and posting this question in a different way; this time from the point of view of Nishita paper.

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:

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

Where

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. 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?):

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:

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

Answer 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