Atmospheric Scattering Shader

This tutorial finally concludes our journey to simulate Rayleigh Scattering for a planet’s atmosphere.

The next (and final) part will show how to change the shader to also include an additional type of scattering, known as Mie Scattering.

You can find all the post in this series here:

You can download the Unity package for this tutorial at the bottom of the page.

Sampling the View Ray

Let’s recall the equation for the atmospheric scattering that we have recently derived:

    \[I= I_S  \sum_{P \in \overline{AB}} {S\left(\lambda, \theta, h\right)  T\left(\overline{CP}\right)  T\left(\overline{PA}\right) ds }\]

The amount of light that we receive is equal to the amount of light emitted from the sun, I_S, multiplied by the sum of the individual contributions from each point P in the segment \overline{AB}.

We could go one and implement that function directly in our shader. However, there are few optimisations that can be done. It was hinted, in a previous tutorial, that the expression could be simplified even further. The first step we can take is to decompose the scattering function into its two basic components:

    \[S \left(\lambda, \theta, h\right ) = \beta \left(\lambda, h \right ) \gamma\left(\theta\right) = \beta \left(\lambda\right )\rho\left(h\right) \gamma\left(\theta\right)\]

The phase function \gamma\left(\theta\right) and the scattering coefficient at sea level \beta \left(\lambda\right ) are constant with respects to the summation, since the angle \theta and the wavelength \lambda do not depend on the sampled point. Hence, they can be factored out:

    \[I = I_S  \, \beta \left(\lambda\right )  \gamma\left(\theta\right) \sum_{P \in \overline{AB}} {   T\left(\overline{CP}\right)  T\left(\overline{PA}\right) \rho\left(h\right) ds }\]

This new expression is mathematically equivalent to the previous one, but is more efficient to calculate since some of the most heavy parts have been taken out of the summation.

We are not ready to start implementing it. There are infinitely many points P we should take into account. A reasonable approximation to I is to split \overline{AB} into several smaller segments of length ds, and to accumulate the contribution of each individual segment. In doing so, we assume that each segment is small enough to have a constant density. Generally speaking, that is not the case, but if the ds is small enough, we can still reach a reasonably good approximation.

The number of segments in  \overline{AB} is referred to as the view samples, since all segments lie on the view ray. In the shader, this will be the _ViewSamples property. By having it as a property, it is accessible from the material inspector. This allows us to reduce the precision of the shader, in favour of its performance.

The following piece of code allows looping through all the segments in the atmosphere.

// Numerical integration to calculate
// the light contribution of each point P in AB
float3 totalViewSamples = 0;
float time = tA;
float ds = (tB-tA) / (float)(_ViewSamples);
for (int i = 0; i < _ViewSamples; i ++)
{
	// Point position
	// (sampling in the middle of the view sample segment)
	float3 P = O + D * (time + ds * 0.5);

	// T(CP) * T(PA) * ρ(h) * ds
	totalViewSamples += viewSampling(P, ds);

	time += ds;
}

// I = I_S * β(λ) * γ(θ) * totalViewSamples
float3 I = _SunIntensity *  _ScatteringCoefficient * phase * totalViewSamples;

The variable time is used to keep track of how far we are from the origin point O, and is incremented by ds after each iteration.

Optical Depth PA

Each point along the view ray \overline{AB} brings its own contribution to the final colour of the pixel that we are drawing. That contribution is, mathematically speaking, the quantity inside the summation:

    \[I = I_S  \, \beta \left(\lambda\right )  \gamma\left(\theta\right) \sum_{P \in \overline{AB}}   \underset{\text{light contribution of}\,L\left(P\right)}{\underbrace{T\left(\overline{CP}\right)  T\left(\overline{PA}\right) \rho\left(h\right) ds}}\]

Like we did in the previous paragraph, let’s try to simplify it. We can expand the expression above even further, by replacing T with its actual definition:

    \[T\left(\overline{XY}\right) =\exp\left\{- \beta\left(\lambda\right)D\left(\overline{XY}\right)\right\}\]

The product of the transmittance over \overline{CP} and \overline{PA} becomes:

    \[T\left(\overline{CP}\right) T\left(\overline{PA}\right)=\]

    \[=\underset{T\left(\overline{CP}\right) }{\underbrace{\exp\left\{-\beta\left(\lambda\right)D\left(\overline{CP}\right)\right \}}} \,\underset{T\left(\overline{PA}\right) }{\underbrace{\exp\left\{-\beta\left(\lambda\right)D\left(\overline{PA}\right)\right \}}}=\]

    \[=\exp\left\{-\beta\left(\lambda\right)\left(D\left(\overline{CP}\right)+D\left(\overline{PA}\right)\right)\right \}\]

The combined transmittance is modelled as an exponential decay with which coefficient is the sum of the optical depths over the path travelled by the light (\overline{CP} and \overline{PA}), multiplied by the scattering coefficient at sea level (\beta with h=0).

The first quantity that we start calculating is the optical depth for the segment \overline{PA}, which goes from the point of entry through the atmosphere, up to the point we are currently sampling in the for loop. Let’s recall the definition of optical depth:

    \[D\left( \overline{PA}\right)=\sum_{Q \in \overline{PA}}{\exp\left\{-\frac{h_Q}{H}\right\}}\, ds\]

If one had to implement this naively, we would create a function called opticalDepth that samples points between P and A in a loop. This is possible, but very inefficient. As a matter of fact, D\left( \overline{PA}\right) is the optical depth of the segment that we are already analysing the outermost for loop. We can save a lot of calculations if we calculate the optical depth the current segmenet centred at P (opticalDepthSegment), and keep accumulating it in the for loop (opticalDepthPA).

// Accumulator for the optical depth
float opticalDepthPA = 0;

// Numerical integration to calculate
// the light contribution of each point P in AB
float time = tA;
float ds = (tB-tA) / (float)(_ViewSamples);
for (int i = 0; i < _ViewSamples; i ++)
{
	// Point position
	// (sampling in the middle of the view sample segment)
	float3 P = O + D * (time + viewSampleSize*0.5);

	// Optical depth of current segment
	// ρ(h) * ds
	float height = distance(C, P) - _PlanetRadius;
	float opticalDepthSegment = exp(-height / _ScaleHeight) * ds;

	// Accumulates the optical depths
	// D(PA)
	opticalDepthPA += opticalDepthSegment;

	...	

	time += ds;
}

Light Sampling

If we look back at the expression for the light contribution of P, we see that the only quantity needed is the optical depth of the segment \overline{CP},

    \[L\left(P\right) =\underset{\text{combined transmittance}}{\underbrace{\exp\left\{-\beta\left(\lambda\right)\left(D\left(\overline{CP}\right)+D\left(\overline{PA}\right)\right)\right \}}}\, \underset{\text{optical depth of}\,ds}{\underbrace{\rho\left(h\right) ds}}\]

We will move the code that calculates the optical depth of the segment \overline{CP} into a function called lightSampling. The name comes from the light ray, which is the segment that starts at P and is directed towards the sun. We have called the point where it exits the atmosphere is C.

The lightSampling function, however, will not just calculate the optical depth of \overline{CP}. So far we have only considered the contribution of the atmosphere, ignoring the role of the actual planet. Our equations do not account for the possibility that the light ray that goes from P towards the sun might hit the planet. If this happens, all calculations done so far must be discarded, as no light would actually reach the camera.

In the diagram above, it’s easy to see that the light contribution of P_0 should be ignored, since no sun light is reaching P_0. While looping through the points in between P and C, the lightSampling function will also check if the planet has been hit. This can be done by checking if the altitude of a point is negative.

bool lightSampling
(	float3 P,	// Current point within the atmospheric sphere
	float3 S,	// Direction towards the sun
	out float opticalDepthCA
)
{
	float _; // don't care about this one
	float C;
	rayInstersect(P, S, _PlanetCentre, _AtmosphereRadius, _, C);

	// Samples on the segment PC
	float time = 0;
	float ds = distance(P, P + S * C) / (float)(_LightSamples);
	for (int i = 0; i < _LightSamples; i ++)
	{
		float3 Q = P + S * (time + lightSampleSize*0.5);
		float height = distance(_PlanetCentre, Q) - _PlanetRadius;
		// Inside the planet
		if (height < 0)
			return false;

		// Optical depth for the light ray
		opticalDepthCA += exp(-height / _RayScaleHeight) * ds;

		time += ds;
	}

	return true;
}

The function above first calculates the point C using rayInstersect. It then divides the segment \overline{PA} into _LightSamples segments of length ds. The calculation for the optical depth is the same used in the outermost loop.

The function returns false if the planet has been hit. We can use this to update the missing code the outermost loop, replacing the ....

	// D(CP)
	float opticalDepthCP = 0;
	bool overground = lightSampling(P, S);

	if (overground)
	{
		// Combined transmittance
		// T(CP) * T(PA) = T(CPA) = exp{ -β(λ) [D(CP) + D(PA)]}
		float transmittance = exp
		(
			-_ScatteringCoefficient *
			(opticalDepthCP + opticalDepthPA)
		);

		// Light contribution
		// T(CPA) * ρ(h) * ds
		totalViewSamples += transmittance * opticalDepthSegment;
	}

Now that we have taken into account all the elements, our shader is complete.

Coming Next…

This post (finally!) completes the volumetric shader that simulates atmospheric scattering. So far, we have only taken into account the contribution of the Rayleigh scattering. There are many optical phenomena that cannot be explained with Rayleigh scattering alone. The next post will introduce a second type of scattering, which is known as Mie scattering.

You can find all the post in this series here:

Download

Become a Patron!
You can download all the assets necessary to reproduce the volumetric atmospheric scattering presented in this tutorial.

FeatureStandardPremium
Volumetric Shader
Clouds Support
Day & Night Cycle
Earth 8K
Mars 8K
Venus 8K
Neptune 2K
Postprocessing
DownloadStandardPremium

Comments

10 responses to “Atmospheric Scattering Shader”

  1. Could you do tutorials on generating planet landscapes using noise generation? I don’t mean for the purpose of landing on them, but like for generating them as a backdrop. Don’t know if you ever saw this asset before it became deprecated for unity but was top notch in its day:

    https://www.youtube.com/watch?v=CbVuHbQ6PAk

    He uses noise gen in run time, and it seems like an art in an of itself to mastur noise gen but most people just do ” here’s a couple perlin noise layers and now you have a planet ” but they never go into how you can control the noise to get things like islands/continents etc and not just some random noise.

    1. Funny thing, I’m literally working on something noise-related RIGHT now!
      I would suggest having a look at RedBlobGames’ amazing blog post:
      https://www.redblobgames.com/maps/terrain-from-noise/

      It might help you with what you need!

  2. Hi, I just want to say big THANK to you for this detailed tutorial.
    I followed it step by step (1-7 in sereis, Mie scattering not yet) and wrote a Unity [URP] version shader.
    I prepare a open sourced repo on Github so if anyone want to check it out.
    https://github.com/riveranb/AtmosphericScatteringWeekend

    At first I was confused about if I’m doing something wrong so the result image is quite dark. And I tried to accumulate optical depth in the case of ‘overground == false’, then I found that dark-side hemisphere of planet also got brightened. So I guess the basic assumption is ‘SINGLE scattering’ therefore if light source is blocked by planet there should be no optical integration at that sample point.
    This is my understanding from your tutorial series. If I have enough time I definitely will access your patreon contents and give Mie scattering and other cool stuff a try.

    1. Also if there are any violations to your license, please let me know 🙂

  3. Hey Alan, wanted to thank you for the great tutorial. I was able to create my own atmospheric renderer thanks in large part to this series (you can see it in action here: http://davidson16807.github.io/tectonics.js/).

    I wanted to point out something I found: you can dramatically improve performance by using a closed form approximation for the “lightSampling()” function in your tutorial. I wrote up a blog post about this, which I hope will make sense: http://davidson16807.github.io/tectonics.js//2019/03/24/fast-atmospheric-scattering.html

    1. Thank you, that is really interesting!
      I hope other people will be able to use your implementation!

  4. Steffen avatar

    How does this shader look like when inside the atmosphere?

    1. Hey!
      This particular version doesn’t really work inside the atmosphere.
      But it can be changed very easily. The first step is to have a material that does not cull back faces.

  5. […] Atmospheric Scattering Shader […]

Leave a Reply

Your email address will not be published. Required fields are marked *