# Using Monte carlo on Rayleigh scattering

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

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 +

}

public Vector ComputeIncidentLight(Ray r, Random rand)
{

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,

• 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. Aug 23 '17 at 7:22
• 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
Aug 26 '17 at 9:21
• The comment was directed toward people willing to propose answers. But thanks, it makes the question more interesting too. :) Sep 1 '17 at 0:29
• 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 ? Mar 5 '18 at 8:06
• 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) Mar 5 '18 at 8:55

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:

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