Intersecting The Atmosphere

As discussed before, the only way we can calculate the optical depth of a segment that passes through the atmosphere, is via a numerical integration. This means dividing out interval in smaller segments of length $ds$ , and calculating the optical depth of each one assuming its density is constant.

In the image above, the optical depth of the $\overline{AB}$ is calculated with 4 samples, each one only taking into account the density at the centre of the segment itself.

The first step is, obviously, finding the points $A$ and $B$ . If we operate under the assumption that the objects we are rendering a sphere, Unity will try to render its surface. Each pixel on the screen corresponds to a point on the sphere. In the diagram below, that point is called $O$ for origin. In a surface shader, $O$ corresponds to the worldPos variable inside the Input structure. This is how far the shader goes; the only piece of information available to us are $O$ , the direction $D$ which indicates the direction of the view ray, and the atmospheric sphere centred at $C$ with radius $R$ . The challenge is to calculate $A$ and $B$ . The fastest way is to use a geometrical approach, reducing the problem to finding the intersection between the atmospheric sphere and the view ray from the camera.

First of all, we should notice that $O$ , $A$ and $B$ are all lying on the view ray. This means that we can refer to their position not as a point in 3D space, but as the distance on the view ray from the origin. While $A$ is the actual point (represented in the shader as a float3), $AO$ will be its distance from the origin $O$ (as a float). Both $A$ and $AO$ are two equally valid ways to indicate the same point, and it holds that:

$A = O + \overline{AO}\,D$

$B = O + \overline{BO}\,D$

where the overline notation $\overline{XY}$ indicates the length of the segment between some arbitrary points $X$ and $Y$ .

For efficiency reasons, in the shader code we will use $AO$ and $BO$ , and calculate it from $OT$ :

$\overline{AO} = \overline{OT} - \overline{AT}$

$\overline{BO} = \overline{OT} + \overline{BT}$

We should also notice that the segments $\overline{AT}$ and $\overline{BT}$ have the same length. What we need now to find the intersection points is to calculate $\overline{AO}$ and $\overline{AT}$ .

The segment $\overline{OT}$ the easiest to calculate. If we look at the diagram above, we can see that $\overline{OT}$ is the projection of the vector $CO$ onto the view ray. Mathematically this projection can be done with the dot product. If you are familiar with shaders you might know the dot product as a measure of how “aligned” two directions are. When it is applied to two vectors and the second one has unitary length, it becomes a projection operator:

$\overline{OT} = \left(C-O\right) \cdot D$

One should notice that $\left(C-O\right)$ is a 3D vector, not the length of the segment between $C$ and $O$ .

What we need to calculate next is the length of the segment $\overline{AT}$ . This can be calculated using Pythagoras’ theorem on the triangle $\overset{\triangle}{ACT}$ . In fact, it holds that:

$R^2 = \overline{AT}^2 + \overline{CT}$

which means that:

$\overline{AT} = \sqrt{R^2 - \overline{CT}}$

The length of $\overline{CT}$ is still unknown. However, it can be calculated by applying once again Pythagoras’ theorem on the triangle $\overset{\triangle}{OCT}$ :

$\overline{CO}^2 = \overline{OT}^2 + \overline{CT}^2$

$\overline{CT} = \sqrt{\overline{CO}^2 - \overline{OT}^2}$

We know have all the quantities that we need. To sum it up:

$\overline{OT} = \left(C-O\right) \cdot D$

$\overline{CT} = \sqrt{\overline{CO}^2 - \overline{OT}^2}$

$\overline{AT} = \sqrt{R^2 - \overline{CT}^2}$

$\overline{AO} = \overline{OT} - \overline{AT}$

$\overline{BO} = \overline{OT} + \overline{AT}$

That set of equation contains square roots. They are only defined on non-negative numbers. If $R^2 > \overline{CT}^2$ , then there is no solution, meaning that the view ray does not intersect the sphere.

We can translate this into the following Cg function:

bool rayIntersect
(
	// Ray
	float3 O, // Origin
	float3 D, // Direction

	// Sphere
	float3 C, // Centre
	float R,	// Radius
	out float AO, // First intersection time
	out float BO  // Second intersection time
)
{
	float3 L = C - O;
	float DT = dot (L, D);
	float R2 = R * R;

	float CT2 = dot(L,L) - DT*DT;
	
	// Intersection point outside the circle
	if (CT2 > R2)
		return false;

	float AT = sqrt(R2 - CT2);
	float BT = AT;

	AO = DT - AT;
	BO = DT + BT;
	return true;
}

There is not a single value, but three to return: $\overline{AO}$ , $\overline{BO}$ and whether or not there is an intersection. The lengths of those two segments are returned using the out keywords, which makes any change the function does on those parameters persistent after its termination.

⭐ Suggested Unity Assets ⭐

Unity is free, but you can upgrade to Unity Pro or Unity Plus subscriptions plans to get more functionality and training resources to power up your projects.

Colliding with the Planet

There is an additional issue that we have to take into account. Certain view rays hit the planet, hence their journey through the atmosphere reaches an early termination. One approach could be to revise the derivation presented above.

An easier, yet less efficient approach, is to run rayIntersect twice, and then to adjust the ending point if needed.

This translates to the following code:

// Intersections with the atmospheric sphere
float tA;	// Atmosphere entry point (worldPos + V * tA)
float tB;	// Atmosphere exit point  (worldPos + V * tB)
if (!rayIntersect(O, D, _PlanetCentre, _AtmosphereRadius, tA, tB))
	return fixed4(0,0,0,0); // The view rays is looking into deep space

// Is the ray passing through the planet core?
float pA, pB;
if (rayIntersect(O, D, _PlanetCentre, _PlanetRadius, pA, pB))
	tB = pA;

Coming Next…

This post showed how it is possible to find intersections between a sphere and a ray. We will use this in the next post to calculate the entrance and exit points of the view ray with the atmospheric sphere.

You can find all the post in this series here:

Part 1. Volumetric Atmospheric Scattering
Part 2. The Theory Behind Atmospheric Scattering
Part 3. The Mathematics of Rayleigh Scattering
Part 4. A Journey Through the Atmosphere
Part 5. A Shader for the Atmospheric Sphere
Part 6. Intersecting The Atmosphere
Part 7. Atmospheric Scattering Shader
🔒 Part 8. An Introduction to Mie Theory

Download

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

Feature	Standard	Premium
Volumetric Shader	✅	✅
Clouds Support	❌	✅
Day & Night Cycle	❌	✅
Earth 8K	✅	✅
Mars 8K	❌	✅
Venus 8K	❌	✅
Neptune 2K	❌	✅
Postprocessing	❌	✅
Download	Standard	Premium

💖 Support this blog

This website exists thanks to the contribution of patrons on Patreon. If you think these posts have either helped or inspired you, please consider supporting this blog.

Follow @AlanZucconi

📧 Stay updated

You will be notified when a new tutorial is released!

📝 Licensing

You are free to use, adapt and build upon this tutorial for your own projects (even commercially) as long as you credit me.

You are not allowed to redistribute the content of this tutorial on other platforms, especially the parts that are only available on Patreon.

If the knowledge you have gained had a significant impact on your project, a mention in the credit would be very appreciated. ❤️🧔🏻

Philip

April 1, 2018

Hi Alan, this is a great tutorial series. I am getting a page not found for part 8 and also the cheat sheet, fyi. Not sure if they’re just not done yet, or if the links are broken.

Part8 about Mie: https://www.alanzucconi.com/?p=7578
Cheat Sheet: https://www.alanzucconi.com/?p=7766

Reply to Philip
- Alan Zucconi
  
  April 4, 2018
  
  Hey!
  Unfortunately, I haven’t had the chance to finalise them yet! T_T
  
  Reply to Alan
Max

December 3, 2018

Hello, Alan. Thank you for a great tutorial.

I’ve got a question about _PlanetCentre variable. Are you setting it with material.SetVector from a C# script? It’s a bit confusing, because you calculated centre and added it to Input struct in a previous chapter. Are you not using that? And if you are, how do you pass a value from Input struct to a lighting function?

Reply to Max
- Alan Zucconi
  
  December 4, 2018
  
  Hey!
  Most the properties (such as _PlanetCentre) are properties that I initialise manually directly in the Material inspector!
  
  Reply to Alan
Tamas

October 2, 2019

Hi Alan!

I’d like to ask in which function does the last snippet go?
I’m new to shader writing and got lost at this part.

Thank you for your clarification!

Reply to Tamas

Published

October 10, 2017

Alan Zucconi in Maths, Shaders, Tutorial, Unity | October 10, 2017