Volumetric atmospheric light scattering

Original author: Alan Zucconi
  • Transfer

If you have lived on planet Earth long enough, you probably wondered why the sky is usually blue, but blushes at sunset. The optical phenomenon that has become the (main) reason for this is called Rayleigh scattering . In this article, I will describe how to simulate atmospheric scattering in order to simulate many of the visual effects that appear on planets. If you want to learn how to render physically accurate images of alien planets, then this tutorial is definitely worth exploring.


The article is divided into the following parts:

  • Part 1. Volumetric atmospheric scattering
  • Part 2. Theory of atmospheric scattering
  • Part 3. Mathematics of Rayleigh scattering
  • Part 4. Traveling through the atmosphere
  • Part 5. Atmospheric Shader
  • Part 6. Intersection of the atmosphere
  • Part 7. Atmospheric dispersion shader

Part 1. Volumetric atmospheric scattering


It is so difficult to recreate atmospheric phenomena because the sky is not an opaque object. Traditional rendering techniques assume that objects are just empty shells. All graphic calculations are performed on the surfaces of materials and do not depend on what is inside. This powerful simplification allows the rendering of opaque objects very efficiently. However, the properties of some materials are determined by the fact that light can pass through them. The final appearance of translucent objects is the result of the interaction of light with their internal structure. In most cases, this interaction can be very effectively simulated, as can be seen in the tutorial " Fast Shader for Subsurface Scattering in Unity". Unfortunately, in our case, if we want to create a convincing sky, this is not so. Instead of rendering only the" outer shell "of the planet, we have to simulate what happens to the rays of light passing through the atmosphere. Performing calculations inside an object is called volume rendering ; this subject we discussed in detail in a series of articles Volumetric Rendering in this article series, I talked about the two techniques (. raymarching and iconic function of the distance ), which can not be effectively used for the simulation of atmospheric dispersion in this article. we will look at more suitable for rendering translucent objects solid technique, often referred to as a single volumetric scattering .

Single scattering

In a room without light, we will not see anything. Objects become visible only when light rays are reflected from them and enter our eye. Most game engines (such as Unity and Unreal) assume that light moves in a "vacuum." This means that only objects can influence light. In fact, light always moves in the environment. In our case, such a medium is the air we inhale. Therefore, the distance traveled by light in air affects the appearance of an object. On the surface of the Earth, the density of air is relatively small, its influence is so insignificant that it should be taken into account only when light travels long distances. Distant mountains merge with the sky, but atmospheric dispersion has almost no effect on nearby objects.

The first step in recreating the optical effect of atmospheric scattering is to analyze how light passes through media such as air. As stated above, we can only see an object when light falls into our eyes. In the context of 3D graphics, our eye is the camera used to render the scene. The molecules that make up the air around us can reflect light rays passing through them. Therefore, they are able to change the way we perceive objects. If greatly simplified, there are two ways that molecules can affect our eyesight.

Outward scattering

The most obvious way molecules interact with light is that they reflect light, changing its direction. If a ray of light directed into the camera is reflected, then we observe a process of scattering outward .

A real light source can emit quadrillion photons every second, and each of them with a certain probability can collide with an air molecule. The denser the medium in which the light moves, the more likely is the reflection of a single photon. The effect of outward scattering also depends on the distance traveled by the light.

Outward scattering leads to the fact that the light gradually becomes more and more dim, and this property depends on the distance traveled and the density of air.

Scattering inward

When light is reflected by a particle, it may happen that it is redirected to the camera. Such an effect is the opposite of outward scattering, and it is logical that it is called inward scattering .

Under certain conditions, scattering inward allows you to see light sources that are not in the direct field of view of the camera. The most obvious result of this is the optical effect - light halo around light sources. They arise due to the fact that the camera receives from one source both direct and indirect rays of light, which de facto increases the number of photons received.

Unit Volumetric Scattering

One ray of light can be reflected an arbitrary number of times. This means that before you get into the camera. the beam can travel very difficult routes. This becomes a serious difficulty for us, because rendering high-quality translucent materials requires simulating the paths of each of the individual rays of light. This technique is called raytracing , and is currently too costly to implement in real time. The unit scattering technique presented in this tutorial takes into account a single light ray scattering event. Later we will see that such a simplification still allows you to get realistic results in just a small fraction of the calculations necessary for real ray tracing.

The basis of rendering a realistic sky is to simulate what happens to the rays of light when passing through the atmosphere of a planet. The diagram below shows a camera looking through a planet. The main idea of ​​this rendering technique is to calculate how light travels from$ A $ in $ B $scattering affects. This means that it is necessary to calculate the effect that the scattering in and out has on the beam passing to the camera. As stated above, due to scattering to the outside, we observe attenuation. The amount of light present at each point$ P \ in \ overline {AB} $ It has a small probability of reflection away from the camera.

To correctly take into account the amount of outward scattering occurring at each point $ P $, we first need to find out how much light was present in $ P $. If we assume that only one star illuminates the planet, then all the light received$ P $must come from the sun. Part of this light will be subject to scattering inward and will be reflected in the camera:

These two steps are enough to approximate most of the effects observed in the atmosphere. However, everything is complicated by the fact that the amount of light received$ P $ from the sun, also subject to scattering outward, when light passes through the atmosphere from $ C $ in $ P $.

To summarize what we need to do:

  • The scope of the camera enters the atmosphere in $ A $ and is located in $ B $;
  • As an approximation, we will take into account the effect of scattering in and out when it occurs at each point $ P \ in \ overline {AB} $;
  • The amount of light received $ P $ from the sun;
  • The amount of light received $ P $ and subject to scattering outward when passing through the atmosphere $ \ overline {CP} $;
  • Part of the Light Received $ P $ and prone to scattering inward, which redirects the rays into the chamber;
  • Part of the world from $ P $directed into the chamber is scattered outward and reflected from the field of view.

This is not the only way that rays of light can enter the camera!
The solution proposed in this tutorial takes into account scattering inward along the scope $ \ overline {AB} $. Light reaching the point$ P $from the sun, with a certain probability can be reflected in the camera.

However, there are many ways in which rays can enter the camera. For example, one of the rays scattered outward on the way to$ P $, in a second collision, it can scatter back towards the camera (see diagram below). And there may be rays reaching the camera after three reflections, or even four.

The technique we are considering is called single scattering , because it takes into account only scattering inward along the field of view. More sophisticated techniques can expand this process and take into account other ways in which the beam can reach the camera. However, the number of possible paths grows with an increase in the number of scattering events considered exponentially. Fortunately, the likelihood of getting into the camera also decreases exponentially.

Part 2. Theory of atmospheric scattering

In this part, we will begin to derive the equations governing this complex but wonderful optical phenomenon.

Transmission function

To calculate the amount of light transmitted to the camera, it is useful to take the same journey that the rays from the sun pass. Looking at the diagram below, it is easy to see that the rays of light reaching$ C $go through empty space. Since there is nothing in their path, all the light reaches$ C $without undergoing any scattering effect. Let's denote the variable$ I_c $ the amount of non-scattered light received by the dot $ C $from the sun. During his journey to$ P $, the ray enters the atmosphere of the planet. Some rays collide with molecules hanging in the air and are reflected in different directions. As a result, part of the light is scattered out of the way. Amount of light reaching$ P $ (designated as $ I_P $) will be less than $ I_c $.

Ratio $ I_c $ and $ I_P $called pass :

$ T \ left (\ overline {CP} \ right) = \ frac {I_P} {I_C} $

We can use it to indicate the percentage of light that was not scattered (i.e. was missed ) when traveling from$ C $ before $ P $. Later we explore the nature of this function. So far, the only thing we need to know is that it corresponds to the effect of atmospheric scattering outside.

Consequently, the amount of light received$ P $is equal to:

$ I_P = I_C \, T \ left (\ overline {CP} \ right) $

Scatter function

Point $ P $receives light directly from the sun. However, not all this light passing through$ P $transmitted back to the camera. To calculate the amount of light actually reaching the camera, we introduce a new concept: the scattering function$ S $. It will denote the amount of light reflected in a certain direction. If we look at the diagram below, we will see that only rays reflected at an angle will be directed to the camera$ \ theta $.

Value $ S \ left (\ lambda, \ theta, h \ right) $ denotes the fraction of light reflected on $ \ theta $radian. This function is our main problem, and we will examine its nature in the next part. So far, the only thing you need to know is that it depends on the color of the light received (determined by the wavelength$ \ lambda $), scattering angle$ \ theta $and heights$ h $ points $ P $. Altitude is important because atmospheric density varies as a function of altitude. And density is one of the factors determining the amount of light scattered.

Now we have all the necessary tools for writing a general equation showing the amount of light transmitted from$ P $ to $ A $:

$ I_ {PA} = \ boxed {I_P} \, S \ left (\ lambda, \ theta, h \ right) \, T \ left (\ overline {PA} \ right) $

Thanks to the previous definition, we can expand $ I_P $:

$ I_ {PA} = \ boxed {I_C \, T \ left (\ overline {CP} \ right)} \, S \ left (\ lambda, \ theta, h \ right) \, T \ left (\ overline { PA} \ right) = $

$ = \ underset {\ text {in-scattering}} {\ underbrace {I_C \, S \ left (\ lambda, \ theta, h \ right)}}, \ underset {\ text {out-scattering}} { \ underbrace {T \ left (\ overline {CP} \ right) \, T \ left (\ overline {PA} \ right)}} $

The equation speaks for itself:

  • Light travels from the sun to $ C $without scattering in the vacuum of space;
  • Light enters the atmosphere and passes from $ C $ to $ P $. In this process, due to the scattering outward only part$ T \ left (\ overline {CP} \ right) $ reaches the destination;
  • Part of the light that got from the sun to $ P $reflected back to the camera. The fraction of light subject to inward scattering is$ S \ left (\ lambda, \ theta, h \ right) $;
  • The remaining light passes from $ P $ before $ A $, and again only part $ T \ left (\ overline {PA} \ right) $.

Numerical integration

If you carefully read the previous paragraphs, you might notice that the brightness was recorded in different ways. Symbol$ I_ {PA} $ denotes the amount of light transmitted from $ P $ to $ A $. However, the value does not take into account all the light that receives$ A $. In these simplified models of atmospheric scattering, we take into account scattering inward from each point along the field of view of the camera in the atmosphere.

The total amount of light received$ A $ ($ I_a $), is calculated by adding up the influence of all points $ P \ in \ overline {AB} $. From the point of view of mathematics, on the segment$ \ overline {AB} $there are an infinite number of points; therefore, cyclic traversal of all of them is impossible. However we can divide$ \ overline {AB} $ for smaller lengths $ ds $ (see diagram below), and summarize the impact of each of them.

Such an approximation process is called numerical integration and leads to the following expression:

$I_A = \sum_{P \in \overline{AB}} {I_{PA}\, ds }$

The more points we take into account, the more accurate the final result will be. In reality, in our atmospheric shader, we simply cycle around several points$P_i$ inside the atmosphere, summing up their influence in the overall result.

Why do we multiply by ds?
Here we use the approximation of a continuous phenomenon. The more points we consider, the closer we will be to the real result. Multiplying each point by$ds$, we give its influence weight according to its length. The more points we have, the less important each of them is.

You can look at it differently: multiplying by$ds$, we “average” the influence of all points.

Directional lighting

If the sun is relatively close, then it is best to model it as a point source of light . In this case, the amount received$C$light depends on the distance to the sun. However, if we are talking about planets, we often assume that the sun is so far away that its rays fall on the planet at one angle. If this is true, then we can model the sun as a directional light source . The light received from a directional light source remains constant, regardless of the distance traveled by it. Therefore every point$C$ receives the same amount of light, and the direction relative to the sun is the same for all points.

We can use this assumption to simplify our equations.

Let's replace$I_C$ constant $I_S$, which determines the brightness of the sun .

$I_A = \sum_{P \in \overline{AB}} {\boxed{I_{PA}}\, ds } =$

$= \sum_{P \in \overline{AB}} {\boxed{I_C  \,S\left(\lambda,\theta,h\right) \,T\left(\overline{CP}\right)\,  T\left(\overline{PA}\right)}\, ds }=$

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

There is another optimization we can perform, it includes a scatter function $S\left(\lambda,\theta,h\right)$. If sunlight always falls from one direction, then the angle$\theta$becomes permanent. Below we will see that directed influence can be excluded from the amount$S\left(\lambda,\theta,h\right)$. But for now, we will leave him.

Absorption coefficient

In describing the possible results of the interaction between light and air molecules, we allowed only two options - either through and through, or reflection. But there is a third possibility. Some chemical compounds absorb light. There are many substances in the Earth’s atmosphere that have this property. Ozone, for example, is in the upper atmosphere and actively interacts with ultraviolet radiation. However, its presence has practically no effect on the color of the sky, because it absorbs light outside the visible spectrum. Here on Earth, the effects of light-absorbing substances are often ignored. But in the case of other planets it is impossible to refuse it. For example, the usual coloration of Neptune and Uranus is caused by the presence of a large amount of methane in their atmosphere. Methane absorbs red light, which gives a blue tint.

Почему Солнце стало красным во время урагана Офелия 2017 года?
Если вы живёте в Великобритании, то могли заметить, что во время урагана Офелия солнце становилось красным. Так происходило потому, что Офелия принесла песок из Сахары. Эти крошечные частицы, висевшие в воздухе, усиливали эффект рассеяния. Как мы увидим в следующей части, синий свет рассеивается сильнее, чем красный.

Если посмотреть на цвета видимого спектра (ниже), то легко увидеть, что если рассеивается достаточное количество синего света, то небо и в самом деле может стать жёлтым или красным.

Есть искушение сказать, что цвет окраски неба связан с оттенком жёлтого сахарского песка. Однако такой же эффект можно увидеть во время больших пожаров из-за частиц дыма, которые обычно близки к чёрному цвету.

Часть 3. Математика рэлеевского рассеяния.

In this part, we will get acquainted with the mathematics of Rayleigh scattering - an optical phenomenon due to which the sky looks blue. The equations derived in this part of the equation will be transferred to the shader code of the next part.


In the previous part, we derived an equation that provides a good basis for approximating atmospheric scattering in a shader. However, we have missed the fact that one equation will not give us convincing results. If we need a beautiful looking atmospheric shader, then we need to go a little deeper into mathematics.

The interaction between light and matter is incredibly complex, and we will not be able to fully describe it in an easy way. Actually, modeling atmospheric scattering is very time consuming. Part of the problem is that the atmosphere is not a homogeneous environment. Its density and composition vary significantly as a function of height, which makes it almost impossible to create an “ideal” model.

That is why several scattering models are given in the scientific literature, each of which is intended to describe a subset of the optical phenomena that arise under certain conditions. Most of the optical effects demonstrated by the planets can be recreated, taking into account only two different models: Rayleigh scattering and scattering of light by a spherical particle . These two mathematical tools make it possible to predict light scattering on objects of various sizes. The first models the reflection of light by oxygen and nitrogen molecules, which make up most of the air. The latter models the reflection of light in larger structures present in the lower atmosphere, such as pollen, dust and pollutants.

Rayleigh scattering is the cause of blue sky and red sunsets. Scattering of light by a spherical particle gives the clouds their white color. If you want to know how this happens, then we will have to dive deeper into the mathematics of scattering.

Rayleigh scattering

What is the fate of a photon striking a particle? To answer this question, we need to rephrase it more formally. Imagine that a ray of light passes through an empty space and suddenly collides with a particle. The result of such a collision is highly dependent on the particle size and color of the light beam. If the particle is small enough (the size of atoms or molecules), then the behavior of light is better predicted using Rayleigh scattering .

What is going on? Part of the world continues its journey, “not feeling” any influence. However, a small percentage of this source light interacts with the particle and scatters in all directions. However, not all directions receive the same amount of light. Photons are more likely to pass directly through the particle or bounce back. That is, the photon reflection option is 90 degrees less large. This behavior can be seen in the diagram below. The blue line shows the most likely directions of the scattered light.

This optical phenomenon is mathematically described by the Rayleigh scattering equation$S \left(\lambda, \theta, h \right )$that gives us a share of the source light $I_0$reflected in the direction $\theta$:

$I = I_0 \, S \left(\lambda, \theta, h\right)$

$S \left(\lambda, \theta, h\right ) = \frac{\pi^2 \left(n^2-1 \right )^2}{2} \underset{\text{density}}{\underbrace{\frac{\rho\left(h\right)}{N}}} \overset{\text{wavelength}}{\overbrace{\frac{1}{\lambda^4}}} \underset{\text{geometry}}{\underbrace{\left(1+\cos^2\theta \right )}}$


  • $\lambda$: wavelength of the incoming light;
  • $\theta$: Angle rasseiyaniya (scattering angle) ;
  • $h$: высота (altitude) точки;
  • $n=1.00029$: коэффициент преломления воздуха;
  • $N=2.504 \cdot 10^{25}$: количество молекул на кубический метр стандартной атмосферы;
  • $\rho\left(h\right)$: коэффициент плотности.Это число на уровне моря равно $1$, и экспоненциально уменьшается с увеличением $h$. Об этой функции можно многое сказать, и мы рассмотрим её в следующих частях.

Но это не уравнение рэлеевского рассеяния!
Если вы встречались с рэлеевским рассеянием не в области компьютерной графики, то есть вероятность, что вы видели другое уравнение. Например, представленное в статье «Рэлеевское рассеяни» на Википедии, сильное отличается.

Использованное в этом туториале уравнение взято из научной статьи Display of The Earth Taking into Account Atmospheric Scattering, Нишиты et al.

Откуда взялась эта функция?
Одна из задач этого блога — объяснение вывода всех получаемых величин. К сожалению, это не относится к рэлеевскому рассеянию.

Если вам всё-таки интересно разобраться, почему частицы света так странно отражаются от молекул воздуха, то интуитивное понимание происходящего может дать следующее объяснение.

Причиной рэлеевского рассеяния на самом деле является не «отталкивание» частиц. Свет — это электромагнитная волна, и она может взаимодействовать с неравновесием зарядов, присутствующим в определённых молекулах. Эти заряды насыщаются поглощением получаемого электромагнитного излучения, которое позже испускается снова. Показанная в угловой функции двудольная форма показывает, как молекулы воздуха становятся электрическими диполями, излучающими подобно микроскопическим антеннам.

The first thing that can be seen in Rayleigh scattering is that in some directions more light travels than in others. The second important aspect is that the amount of scattered light is highly dependent on the wavelength.$\lambda$coming light. The pie chart below renders$S \left(\lambda, \theta, 0\right )$for three different wavelengths. Calculation$S$ at $h=0$often called sea ​​level scattering .

The figure below shows the visualization of the scattering coefficients for the continuous wavelength / color range of the visible spectrum (the code is available in ShaderToy ).

The center of the image appears black because the wavelengths in this range are outside the visible spectrum.

Rayleigh scattering coefficient

The Rayleigh scattering equation shows how much light is scattered in a certain direction. However, it does not tell us how much energy has been dissipated. To calculate this, we need to take into account energy dissipation in all directions. The derivation of the equation is not easy; if you haven’t mastered complex matanalysis, here is the result:

$\beta \left(\lambda, h \right )=\frac{8\pi^3 \left(n^2-1 \right )^2}{3} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4}$

Where $\beta \left(\lambda, h \right )$shows the fraction of energy that is lost due to scattering after a collision with a single particle. Often this value is called the Rayleigh scattering coefficient .

If you read the previous part of the tutorial, you can guess that$\beta$is actually the attenuation coefficient that was used in determining the transmittance$T\left(\overline{AB}\right)$ on the segment $\overline{AB}$.

Unfortunately, the calculation$\beta$very expensive. In the shader that we are going to write, I will try to save as many calculations as possible. For this reason, it will be useful to “extract” from the scattering coefficient all factors that are constants. What gives us$\beta \left(\lambda\right )$, which is called the Rayleigh scattering coefficient at sea level ($h=0$):

$\beta \left(\lambda \right )=\frac{8\pi^3 \left(n^2-1 \right )^2}{3} \frac{1}{N} \frac{1}{\lambda^4}$

Show me the calculations!
The Rayleigh scattering equation shows us the fraction of the energy scattered in a certain direction. To calculate the total loss, we need to consider all possible directions. For summation over a continuous interval, integration is required .

You can try to integrate$S \left(\lambda, \theta, h\right )$ by $\theta$ on the interval $\left[0,2\pi\right]$, но это будет ошибкой.

Несмотря на то, что мы визуализировали рэлеевское рассеяние в двух измерениях, на самом деле это трёхмерное явление. Угол рассеяния $\theta$ может принимать любое направление в 3D-пространстве. Вычисления с учётом полного распределения функции, которое зависит от $\theta$ в трёхмерном пространстве (как $S$) называется интегрированием по телесному углу:

$\beta \left(\lambda, h \right ) = \int_{0}^{2\pi} \int_{0}^{\pi} S \left(\lambda, \theta, h \right ) \sin\theta \, d\theta d\phi$

Внутренний интеграл перемещает $\theta$ в плоскости XY, а внешний — поворачивает результат вокруг оси X, чтобы учесть третье измерение. Прибавляемый $\sin \theta$ используется для сферических углов.

Процесс интегрирования интересует только то, что зависит от $\theta$. Несколько членов $S \left(\lambda, \theta, h \right )$ постоянны, поэтому их можно перенести из под знака интеграла:

$\beta \left(\lambda, h \right )=\int_{0}^{2\pi} \int_{0}^{\pi} \underset{\text{constant}}{\underbrace{ \frac{\pi^2 \left(n^2-1 \right )^2}{2} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4} }} \left(1+\cos^2\theta \right ) \sin\theta \, d\theta d\phi=$

$=\frac{\pi^2 \left(n^2-1 \right )^2}{2} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4} \int_{0}^{2\pi} \int_{0}^{\pi} \left(1+\cos^2\theta \right ) \sin\theta \, d\theta d\phi$

Это значительно упрощает внутренний интеграл, который теперь принимает следующий вид:

$\beta \left(\lambda, h \right )=\frac{\pi^2 \left(n^2-1 \right )^2}{2} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4} \int_{0}^{2\pi} \boxed{ \int_{0}^{\pi} \left(1+\cos^2\theta \right ) \sin\theta \, d\theta }d\phi=$

$=\frac{\pi^2 \left(n^2-1 \right )^2}{2} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4} \int_{0}^{2\pi} \boxed{\frac{8}{3}} d\phi$

Now we can perform external integration:

$\beta \left(\lambda, h \right )=\frac{\pi^2 \left(n^2-1 \right )^2}{2} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4} \boxed{ \int_{0}^{2\pi} \frac{8}{3} d\phi}=$

$=\frac{\pi^2 \left(n^2-1 \right )^2}{2} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4} \boxed{ \frac{16 \pi}{3}}$

Which brings us to the final form:

$\beta \left(\lambda, h \right )=\frac{8\pi^3 \left(n^2-1 \right )^2}{3} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4}$

Since this is an integral that takes into account the effect of scattering in all directions, it is logical that the expression no longer depends on $\theta$.

This new equation gives us another way to understand how different colors are scattered. The graph below shows the amount of scattering to which light is exposed, as a function of its wavelength.

This is a strong relationship between the scattering coefficient. $\beta$ and $\lambda$which causes a red sky at sunset. Photons received from the sun are distributed over a wider range of wavelengths. Rayleigh scattering tells us that the atoms and molecules of the Earth’s atmosphere scatter blue more than green or red. If the light can move long enough, then its entire blue component will be lost due to scattering. This is what happens at sunset, when the light travels almost parallel to the surface.

Reasoning the same way, we can understand why the sky looks blue. The light of the sun is falling from one direction. However, its blue component is scattered in all directions. When we look at the sky, blue light comes from all directions.

Rayleigh phase function

The original equation describing Rayleigh scattering, $S \left(\lambda, \theta\right )$, can be decomposed into two components. The first is the scattering coefficient we just derived$\beta \left(\lambda\right )$which modulates its intensity. The second part is related to the scattering geometry, and controls its direction:

$S \left(\lambda, \theta, h\right ) = \beta \left(\lambda, h\right ) \gamma \left(\theta\right)$

This new value $\gamma \left(\theta\right)$ can be obtained by dividing $S \left(\lambda, \theta, h\right )$ on the $\beta \left(\lambda\right )$:

$\gamma \left(\theta\right) = \frac{S \left(\lambda, \theta, h\right )} {\beta \left(\lambda\right )}=$

$= \underset{S \left(\lambda, \theta, h\right )}{\underbrace{ \frac{\pi^2 \left(n^2-1 \right )^2}{2} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4} \left(1+\cos^2\theta \right ) }} \, \underset{\frac{1}{\beta \left(\lambda\right )}}{\underbrace{ \frac{3}{8\pi^3 \left(n^2-1 \right )^2} \frac{N}{\rho\left(h\right)} \lambda^4 }}=$

$= \frac{3}{16\pi} \left(1+\cos^2 \theta\right)$

You can see that this new expression is independent of the wavelength of the incoming light. This seems counterintuitive because we know for sure that Rayleigh scattering affects shortwaves more strongly.

However$\gamma \left(\theta\right)$models the shape of the dipole that we saw above. Member$\frac{3}{16\pi}$ serves as a normalization coefficient so that the integral over all possible angles $\theta$ summed up in $1$. If formulated technically, we can say that the integral over$4\pi$ steradians equal $1$.

In the following parts, we will see how the separation of these two components allows us to derive more efficient equations.

To summarize

  • Rayleigh scattering equation : indicates the fraction of light reflected in the direction$\theta$. The amount of scattering depends on the wavelength.$\lambda$ incoming light.

$S \left(\lambda, \theta, h\right ) = \frac{\pi^2 \left(n^2-1 \right )^2}{2} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4} \left(1+\cos^2\theta \right )$


$S \left(\lambda, \theta, h\right ) = \beta \left(\lambda,h \right ) \gamma\left(\theta\right)$

  • Rayleigh scattering coefficient : indicates the fraction of light lost due to scattering after the first collision.

$\beta \left(\lambda,h \right )=\frac{8\pi^3 \left(n^2-1 \right )^2}{3} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4}$

  • Rayleigh scattering coefficient at sea level : this is an analog$\beta \left(\lambda,0\right )$. Creating this additional coefficient will be very useful for deriving more efficient equations.

$\beta \left(\lambda \right )=\beta \left(\lambda,0 \right )=\frac{8\pi^3 \left(n^2-1 \right )^2}{3} \frac{1}{N} \frac{1}{\lambda^4}$

If we consider wavelengths roughly corresponding to the red, green, and blue colors, we will obtain the following results:

$\beta\left(680nm\right) = 0.00000519673$

$\beta\left(550nm\right) = 0.0000121427$

$\beta\left(440nm\right) = 0.0000296453$

These results are calculated under the assumption that $h=0$ (this implies that $\rho=1$) This only happens at sea level, where the dispersion coefficient has a maximum value. Therefore, it serves as a "guide" for the amount of light scattering.

  • Rayleigh phase function : controls the scattering geometry, which denotes the relative fraction of light lost in a particular direction. Coefficient$\frac{3}{16\pi}$ serves as a normalization coefficient, therefore, the integral over the unit sphere will be $1$.

$\gamma\left(\theta\right)= \frac{3}{16\pi} \left(1+\cos^2 \theta\right)$

  • Density coefficient : This function is used to simulate the density of the atmosphere. Its formal definition will be presented below. If you are not against mathematical spoilers, then it is defined as follows:


Where $H=8500$meters is called the reduced height .

Part 4. Travel through the atmosphere.

In this part, we will consider modeling atmospheric density at different heights. This is a necessary step, because the density of the atmosphere is one of the parameters necessary for the correct calculation of Rayleigh scattering.

Atmospheric density coefficient

So far, we have not considered the role of the atmospheric density coefficient $\rho$. It is logical that the scattering force is proportional to the density of the atmosphere. The more molecules per cubic meter, the greater the probability of scattering of photons. The difficulty is that the structure of the atmosphere is very complex, it consists of several layers with different pressure, density and temperature. Fortunately for us, the bulk of Rayleigh scattering occurs in the first 60 km of the atmosphere. In the troposphere, temperature decreases linearly, and pressure decreases exponentially.

The diagram below shows the relationship between density and height in the lower atmosphere.

Value $\rho \left(h\right)$ is a measurement of the atmosphere at altitude $h$meters, normalized in a way that starts from zero. In many scientific articles$\rho$also called the density coefficient , because it can also be defined as follows:

$\rho\left(h\right) = \frac{density\left(h\right)}{density\left(0\right)}$

Dividing the true density by $density\left(0\right)$we get that $\rho\left(h\right)$ at sea level is equal $1$. However, as emphasized above, the calculation$density\left(h\right)$far from nontrivial. We can approximate it as an exponential curve; some of you could already understand that the density in the lower atmosphere is decreasing exponentially .

If we want to approximate the density coefficient through an exponential curve, then we can do this as follows:

$\rho\left(h\right) = exp\left\{-\frac{h}{H}\right\}$

Where $H_0$Is a scale factor called the reduced height . For Rayleigh scattering in the lower atmosphere of the Earth, it is often believed that$H=8500$meters (see diagram below). For scattering light by a spherical particle, the value is often approximately equal$1200$ meters.

The value used for $H$does not give the best approximation $\rho\left(h\right)$. However, this is actually not a problem. Most of the values ​​presented in this tutorial undergo serious approximations. To create visually pleasing results, it will be much more efficient to adjust the available parameters to the reference images.

Exponential decrease

In the previous part of the tutorial, we derived an equation showing us how to account for the outward scattering to which a ray of light is exposed after interacting with an individual particle. The value used to model this phenomenon was called the scattering coefficient$\beta$. To take it into account, we introduced the coefficients$\beta$.

In the case of Rayleigh scattering, we also derived a closed form for calculating the amount of light subject to atmospheric scattering in a single interaction:

$\beta \left(\lambda, h \right )=\frac{8\pi^3 \left(n^2-1 \right )^2}{3} \frac{\rho\left(h\right)}{N} \frac{1}{\lambda^4}$

When calculating at sea level, i.e. when $h=0$, the equation gives the following results:

$\beta\left(680nm\right) = 0.00000519673$

$\beta\left(550nm\right) = 0.0000121427$

$\beta\left(440nm\right) = 0.0000296453$

Where $680$, $550$ and $440$- wavelengths roughly corresponding to red, green and blue.

What is the meaning of these numbers? They represent the fraction of light that is lost in a single interaction with a particle. Assuming that a ray of light initially has brightness$I_0$ and passes a segment of the atmosphere with a (general) scattering coefficient $\beta$, then the amount of light remaining after scattering is:

$I_1=\underset{\text{initial energy}}{\underbrace{I_0}} - \underset{\text{energy lost}}{\underbrace{I_0 \beta}}=I_0 \left(1-\beta\right)$

This is true for a single collision, but we are interested in what kind of energy is dissipated during a certain distance. This means that the remaining light undergoes this process at every point.

When light passes through a homogeneous medium with a scattering coefficient$\beta$, how can we calculate how much of it will remain when moving a certain distance?

For those of you who have studied matanalysis, this may seem familiar. When a multiplicative process is repeated on a continuous segment$\left(1-\beta\right)$then the Euler number enters the scene . The amount of light remaining after scattering after passage$x$ meters equals:

$I = I_0 exp \left\{-\beta x \right\}$

And again, we are faced with an exponential function. It has nothing to do with the exponential function describing the density coefficient$\rho$. Both phenomena are described by exponential functions because they are both subject to exponential reduction. Other than that, there is no connection between them.

Where did exp come from?
If you are unfamiliar with matanalysis, you may not understand why such a simple process as $I_1 = I_0 \left(1-\beta\right)$ converted to $I = I_0 exp \left\{-\beta x \right\}$.

We can consider the original expression as the first approximation to the solution. If we want to get a closer approximation, then we need to take into account two scattering events, halving the scattering coefficient for each of them:

$I_1 = I_0 \left(1-\frac{\beta}{2}\right)$

$I_2 = \boxed{I_1} \left(1-\frac{\beta}{2}\right)$

Which leads us to the following:

$I_2 = \boxed{I_0 \left(1-\frac{\beta}{2}\right)} \left(1-\frac{\beta}{2}\right)=$

$=I_0 \left(1-\frac{\beta}{2}\right)^2$

This is a new expression for $I_2$shows the amount of energy saved after two collisions. What about the three collisions? Or four? Or five? In general terms, this can be written as follows:

$I=\lim_{n\rightarrow \infty }I_0 \left(1-\frac{\beta}{n}\right)^n$

Where $\lim_{n\rightarrow \infty }$- This is a mathematical construction that allows you to repeat the process an infinite number of times. We need this because$\infty$ technically is not a chill, as in this context it does not make sense to write something like $\frac{\beta}{\infty}$.

This expression is the very definition of an exponential function:

$\lim_{n\rightarrow \infty } \left(1-\frac{\beta}{n}\right)^n = e^{-\beta}=exp\left\{-\beta\right\}$

which describes the amount of energy stored after the multiplicative process on a continuous interval.

Uniform transmittance

In the second part of the tutorial we introduced the concept of transmission$T$as the fraction of light remaining after the scattering process when passing through the atmosphere. Now we have all the elements to deduce the equation describing it.

Let's look at the diagram below and think about how to calculate the transmittance for the segment$\overline{CP}$. You can easily see that the rays of light reaching$C$pass through empty space. Therefore, they are not scattered. Consequently, the amount of light in$C$actually equal to the brightness of the sun$I_S$. During this journey to$P$part of the light is scattered out of the way; therefore the amount of light reaching$P$ ($I_P$) will be less than $I_S$.

The amount of scattered light depends on the distance traveled. The longer the path, the greater the attenuation. According to the law of exponential reduction, the amount of light at a point$I_P$ can be calculated as follows:

$I_P = I_S  \exp{\left\{-\beta \overline{CP}\right\}}$

Where $\overline{CP}$ - the length of the segment from $C$ before $P$, a $\exp{\left\{x\right\}}$- exponential function$e^{x}$.

Atmospheric transmission

We based our equation on the assumption that the probability of reflection ( scattering coefficient$\beta$ same at every point along $\overline{CP}$. Unfortunately, this is not the case.

The scattering coefficient is highly dependent on the density of the atmosphere. The more air molecules per cubic meter, the higher the chance of a collision. The density of the planet’s atmosphere is heterogeneous, and varies with altitude. It also means that we cannot calculate outward scatterings from$\overline{CP}$in one step. To solve this problem, we need to calculate the outward dispersion at each point using the scattering coefficient.

To understand how this works, let's start with the approximation. Section$\overline{CP}$ divided into two $\overline{CQ}$ and $\overline{QP}$.

First we calculate the amount of light from $C$which reaches $Q$:

$I_Q = I_S \exp{\left\{-\beta{\left(\lambda, h_0\right)} \overline{CQ} \right\}}$

Then we use the same approach to calculate the amount of light reaching $P$ of $Q$:

$I_P = \boxed{I_Q} \exp{\left\{-\beta{\left(\lambda, h_1\right)} \overline{QP} \right\}}$

If we substitute $I_Q$ into the second equation and simplify, we get:

If $\overline{CQ}$ and $\overline{QP}$ have the same length $ds$, then we can simplify the expression even more:

In the case of two segments of the same length with different scattering coefficients, the outward scattering is calculated by summing the scattering coefficients of the individual segments and multiplied by the lengths of the segments.

We can repeat this process with an arbitrary number of segments, getting closer and closer to the true value. This will lead us to the following equation:

$I_P = I_S\exp\left\{ -\boxed{ \sum_{Q \in \overline{CP}} { \beta\left( \lambda, h_Q \right) } \, ds } \right \} $

Where $h_Q$ - point height $Q$.

The method of dividing a straight line into many segments, which we just used, is called numerical integration .

If we assume that the initial amount of light received is$1$, then we obtain the equation of atmospheric transmission through an arbitrary segment:

$T\left(\overline{CP}\right) =\exp\left\{ -\sum_{Q \in \overline{CP}} { \beta\left( \lambda, h_Q \right) } \, ds \right \} $

We can expand this expression even further. replacing the total$\beta$ the real value used in Rayleigh scattering, $\beta$:

$T\left(\overline{CP}\right) =\exp\left\{ -\sum_{Q \in \overline{CP}} { \boxed{\frac{8\pi^3 \left(n^2-1 \right )^2}{3} \frac{\rho\left(h_Q\right)}{N} \frac{1}{\lambda^4}} } \, ds \right \} $

Many factors $\beta$ constant, so that they can be carried out for the amount of:

$T\left(\overline{CP}\right) =\exp\left\{ - \underset{\beta\left(\lambda\right)}{ \underset{\text{constant}}{ \underbrace{ \frac{8\pi^3 \left(n^2-1 \right )^2}{3} \frac{1}{N} \frac{1}{\lambda^4} }}} \overset{\text{optical depth}\,D\left(\overline{CP}\right)}{ \overbrace{ \sum_{Q \in \overline{CP}} { \rho\left(h_Q\right) } \, ds}} \right \}$

The value expressed by summation is called the optical thickness$D\left(\overline{CP}\right)$which we will calculate in the shader. The remainder is a magnification factor that can only be calculated once, and it corresponds to a scattering coefficient at sea level . In the finished shader, we will only calculate the optical thickness, and the scattering coefficients at sea level$\beta$We will transmit as input.


$T\left(\overline{CP}\right) =\exp\left\{ - \beta\left(\lambda\right) D\left(\overline{CP}\right) \right\}$

Part 5. Atmospheric shader.


Writing shader

We can start writing a shader for this effect in an infinite number of ways. Since we want to render atmospheric scattering on planets, it is logical to assume that it will be used on spheres.

If you use this tutorial to create a game, chances are that you will apply a shader to an existing planet. Adding atmospheric scatter calculations on top of the swatch is possible, but usually produces poor results. The reason is that the atmosphere is larger than the radius of the planet, so it needs to be rendered on a transparent sphere of a slightly larger size. The figure below shows how the atmosphere extends over the surface of the planet, mixing with the empty space behind it.

The application of dispersible material to a separate sphere is possible, but redundant. In this tutorial, I suggest expanding Standard Surface Shader Unity by adding a shader pass that renders the atmosphere in a slightly larger sphere. We will call it the atmospheric sphere .

Two-pass shader

If you used to work with Unity surface shaders , you may have noticed that they do not support a block Pass, with which several passes are set in the vertex and fragment shaders .

Creating a shader with two passes is possible, just add two separate sections of the CG code in one block SubShader:

Shader "Custom/NewSurfaceShader" {
 Properties {
 _Color ("Color", Color) = (1,1,1,1)
 _MainTex ("Albedo (RGB)", 2D) = "white" {}
 _Glossiness ("Smoothness", Range(0,1)) = 0.5
 _Metallic ("Metallic", Range(0,1)) = 0.0
 SubShader {
 // --- Первый проход ---
 Tags { "RenderType"="Opaque" }
 LOD 200
 // Здесь код Cg
 // ------------------
 // --- Второй проход ---
 Tags { "RenderType"="Opaque" }
 LOD 200
 // Здесь код Cg
 // -------------------
 FallBack "Diffuse"

You can change the first pass to render the planet. From now on, we will focus on the second pass, realizing atmospheric scattering in it.

Normal Extrusion

The atmospheric sphere is slightly larger than the planet. This means that the second pass should extrude the sphere out. If your model uses smooth normals, then we can achieve this effect using a technique called extrusion of normals .

Normal extrusion is one of the oldest shader tricks, and is usually one of the first to be studied. My blog has many references to it; A good start might be a Surface Shader post from the A Gentle Introduction to Shaders series .

If you are unfamiliar with how the extrusion of normals works, then I will explain: all the vertices are processed by the shader using its vertex function. We can use this function to change the position of each vertex, making the sphere bigger.

The first step is to change the directive by pragmaadding it vertex:vert; this will force Unity to execute a function for each vertex vert.

#pragma surface surf StandardScattering vertex:vert
void vert (inout appdata_full v, out Input o)
 v.vertex.xyz += v.normal * (_AtmosphereRadius - _PlanetRadius);

This code fragment shows a vertex function extruding a sphere along its normals. The magnitude of the extrusion of a sphere depends on the size of the atmosphere and the size of the planet. Both of these values ​​need to be passed to the shader as properties that can be accessed through the material inspector .

Our shader also needs to know where the center of the planet is. We can add its calculation to the vertex function too. We discussed the finding of the center point of an object in the space of the world in the article Vertex and Fragment Shaders .

struct Input
 float2 uv_MainTex;
 float3 worldPos; // Автоматически инициализируется Unity
 float3 centre;   // Инициализируется в вершинной функции
void vert (inout appdata_full v, out Input o)
 v.vertex.xyz += v.normal * (_AtmosphereRadius - _PlanetRadius);
 o.centre = mul(unity_ObjectToWorld, half4(0,0,0,1));

Если посмотреть на вершинную функцию, то можно заметить, что она всегда содержит таинственный вызов UNITY_INITIALIZE_OUTPUT. Шейдер получает позицию вершин в пространстве объекта, и ему нужно спроецировать их в координаты мира с помощью предоставляемых Unity положения, масштаба и поворота.

И это одна из операций, выполняемых UNITY_INITIALIZE_OUTPUT. Без него нам бы пришлось писать необходимый для этих вычислений код самостоятельно.

Аддитивное смешивание

Another interesting feature that we need to deal with is transparency. Usually transparent material allows you to see what is behind the object. This solution in our case will not work, because the atmosphere is not just a transparent piece of plastic. Light passes through it, so we need to use the additive mixing mode to increase the luminosity of the planet.

In the standard Unity Surface Shader, no blending mode is enabled by default. To change the situation, we will have to replace the labels in the second pass with the following:

Tags { "RenderType"="Transparent"
LOD 200
Cull Back
Blend One One

The expression is Blend One Oneused by the shader to access the additive blending mode.

Native lighting function

Most often, when writing a surface shader, programmers change its function surf, which is used to provide “physical” properties, such as albedo, surface smoothness, metal properties, etc. All these properties are then used by the shader to calculate realistic shading.

In our case, these calculations are not required. To get rid of them, we need to remove the lighting model used by the shader . I examined this topic in detail; You can study the following posts to figure out how to do this:

The new lighting model will be called StandardScattering; we must create a function for lighting in real time and for global illumination, that is, respectively, LightingStandardScatteringand LightingStandardScattering_GI.

The code that we need to write will also use properties such as the direction of lighting and the direction of view . They can be obtained using the following code snippet.

#pragma surface surf StandardScattering vertex:vert
#include "UnityPBSLighting.cginc"
inline fixed4 LightingStandardScattering(SurfaceOutputStandard s, fixed3 viewDir, UnityGI gi)
 float3 L = gi.light.dir;
 float3 V = viewDir;
 float3 N = s.Normal;
 float3 S = L; // Направление освещения от солнца
 float3 D = -V;  // Направление луча видимости, проходящего через атмосферу
void LightingStandardScattering_GI(SurfaceOutputStandard s, UnityGIInput data, inout UnityGI gi)
 LightingStandard_GI(s, data, gi); 

It ...will contain the shader code itself, which is needed to implement this effect.

Floating point precision

In this tutorial we will assume that all calculations are performed in meters. This means that if we need to simulate the Earth, then we need a sphere with a radius of 6371000 meters. In fact, in Unity this is not possible due to floating point errors that occur when you have to work with very large and very small numbers at the same time.

To circumvent these limitations, you can scale the scattering coefficient to compensate. For example, if the planet has a radius of only 6.371 meters, then the scattering coefficient
$\beta\left(\lambda\right)$ should be more than 1,000,000 times, and the reduced height $H$- 1,000,000 times less.

In my Unity project, all properties and calculations are expressed in meters. This allows us to use real, physical values ​​for the scattering coefficients and the reduced height. However, the shader also gets the size of the sphere in meters, so that it can perform scale conversion from Unity units to real-scale meters.

Part 6. Intersection of the atmosphere.

As mentioned above, the only way to calculate the optical thickness of a segment passing through the atmosphere is through numerical integration . This means that you need to divide the interval into smaller lengths$ds$ and calculate the optical thickness of each, considering its density constant.

Above is the optical thickness. $\overline{AB}$It is calculated by four samples, each of which takes into account the density in the center of the segment itself.

The first step will obviously be finding the points$A$ and $B$. If we assume that the objects we render are a sphere, then Unity will try to render its surface. Each pixel on the screen corresponds to a point on the sphere. In the figure below, this point is called$O$(from origin ). In the surface shader$O$corresponds to a variable worldPosin the structure Input. This is the amount of work the shader does; the only information available to us is$O$direction $D$, which indicates the direction of the line of sight , and the atmospheric sphere centered$C$ and radius $R$. The difficulty lies in computing$A$ and $B$. The geometry will be used the fastest, reducing the problem to finding the intersection of the atmospheric sphere and the camera visibility beam .

Firstly, it is worth noting that $O$, $A$ and $B$lie on the line of sight. This means that we can handle their position not as a point in 3D space, but as with the distance of the line of sight to origin.$A$Is the real point (represented as in the shader float3), and$AO$ - distance to the starting point $O$(stored as float).$A$ and $AO$ - two equally correct ways of designating the same point, that is, it is true:

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

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

where is the entry with a bar above $\overline{XY}$ denotes the length of the segment between two arbitrary points $X$ and $Y$.

In the shader code for efficiency we will use$AO$ and $BO$, and compute them from $OT$:

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

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

It is also worth noting that the segments $\overline{AT}$ and $\overline{BT}$have the same length. Now, to find the intersection points, we need to calculate$\overline{AO}$ and $\overline{AT}$.

Section$\overline{OT}$The easiest way to calculate. If you look at the diagram above, you can see that$\overline{OT}$ is a projection of a vector $CO$on the line of sight. Mathematically, this projection can be done using a scalar product . If you are familiar with shaders, then you may know that a scalar product is a measure of how “aligned” two directions are. When it is applied to two vectors, and one of them has a unit length, it becomes the projection operator:

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

It is also worth noting that
$\left(C-O\right)$ Is a three-dimensional vector, not the length of the segment between $C$ and $O$.

Next we need to calculate the length of the segment$\overline{AT}$. It can be calculated using the Pythagon theorem in a triangle$\overset{\triangle}{ACT}$. She claims that:

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

And that means that:

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

Length $\overline{CT}$still unknown. However, it can be calculated by again applying the Pythagorean theorem to the triangle$\overset{\triangle}{OCT}$:

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

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

Now we have all the necessary quantities. Summarize:

$\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}$

There are square roots in this set of equations. They are defined only for non-negative numbers. If$R^2 > \overline{CT}^2$, then there is no solution, which means that the line of sight does not intersect with the sphere.

We can convert this to the following Cg function:

bool rayIntersect
 // Ray
 float3 O, // Origin
 float3 D, // Направление
 // Сфера
 float3 C, // Центр
 float R, // Радиус
 out float AO, // Время первого пересечения
 out float BO  // Время второго пересечения
 float3 L = C - O;
 float DT = dot (L, D);
 float R2 = R * R;
 float CT2 = dot(L,L) - DT*DT;
 // Точка пересечения за пределами круга
 if (CT2 > R2)
 return false;
 float AT = sqrt(R2 - CT2);
 float BT = AT;
 AO = DT - AT;
 BO = DT + BT;
 return true;

It returns not one, but three values $\overline{AO}$, $\overline{BO}$as well as the binary value of the intersection. The lengths of these two segments are returned using keywords outthat save after the function completes any changes that it made to these parameters.

Collision with the planet

There is another problem to consider. Some rays of visibility collide with the planet, so their journey through the atmosphere ends early. One solution is to revise the above.

A simpler but less effective approach is to perform rayIntersecttwice and change the endpoint if necessary.

This translates to the following code:

// Пересечения с атмосферной сферой
float tA; // Точка входа в атмосферу (worldPos + V * tA)
float tB; // Точка выхода из атмосферы  (worldPos + V * tB)
if (!rayIntersect(O, D, _PlanetCentre, _AtmosphereRadius, tA, tB))
 return fixed4(0,0,0,0); // Лучи видимости смотрят в глубокий космос
// Проходит ли луч через ядро планеты?
float pA, pB;
if (rayIntersect(O, D, _PlanetCentre, _PlanetRadius, pA, pB))
 tB = pA;

Part 7. Atmospheric dispersion shader.

In this part, we will finally complete the work on simulating Rayleigh scattering in the atmosphere of the planet.


Visibility Beam Sampling

Let's recall the equation of atmospheric scattering that we 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 we receive is equal to the amount of light emitted by the sun, $I_S$multiplied by the sum of the individual influences of each point $P$ on the segment $\overline{AB}$.

We can implement this function directly in the shader. However, some optimizations need to be made here. In the previous part, I said that the expression can be simplified even more. The first thing to do is to break up the scattering function into two main 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)$

Phase function$\gamma\left(\theta\right)$and sea ​​level dispersion coefficient$\beta \left(\lambda\right )$ are constant with respect to the sum, because the angle $\theta$ and wavelength $\lambda$independent of the sampled point. Therefore, they can be carried out for the amount of:

$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 similar to the previous one, but more efficient in calculating, because the most "heavy" parts are derived from the sum.

We are not yet ready for its implementation. There are an infinite number of points$P$which we must consider. Reasonable approximation$I$ will be a split $\overline{AB}$ into several shorter segments $ds$and the addition of the influence of each individual segment. Having done so, we can assume that each segment is small enough to have a constant density. In general, this is not so, but if$ds$ small enough, then we can achieve a fairly good approximation.

The number of segments in $\overline{AB}$we will call the samples of visibility because all segments lie on the line of sight. In the shader, this will be a property _ViewSamples. Since this property, we can access it from the material inspector. This allows us to reduce the accuracy of the shader for its performance.

The following code snippet allows you to bypass all segments in the atmosphere.

// Численное интегрирование для вычисления
// влияния света в каждой точке P отрезка AB
float3 totalViewSamples = 0;
float time = tA;
float ds = (tB-tA) / (float)(_ViewSamples);
for (int i = 0; i < _ViewSamples; i ++)
 // Позиция точки
 // (сэмплирование в середине отрезка сэмпла видимости)
 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 is timeused to track how far we are from the starting point.$O$. After each iteration, it increases by ds.

Optical Thickness PA

Every point on the line of sight $\overline{AB}$contributes to the final color of the pixel we draw. This contribution, expressed mathematically, is a value inside the sum:

$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}}$

Let's try to simplify this expression, as in the last paragraph. We can further expand the above expression by replacing$T$ its definition:

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

Result of skipping over $\overline{CP}$ and $\overline{PA}$ turns into:

$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 transmission is modeled as an exponential decrease, the coefficient of which is the sum of the optical thicknesses along the path traveled by light ($\overline{CP}$ and $\overline{PA}$) times the dispersion coefficient at sea level ($\beta$ at $h=0$)
The first value that we begin to calculate is the optical thickness of the segment$\overline{PA}$, which goes from the point of entry into the atmosphere to the point at which we are currently sampling in the for loop. Let's recall the definition of optical thickness:

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

If we had to implement it “forehead”, then we would create a function opticalDepththat samples the points between$P$ and $A$. It is possible, but very inefficient. In fact,$D\left( \overline{PA}\right)$Is the optical thickness of the line segment that we are already analyzing in the outermost for loop. We can save a lot of calculations if we calculate the optical thickness of the current segment centered in$P$( opticalDepthSegment), and continue to summarize it in a for ( opticalDepthPA) loop .

// Накопитель оптической толщины
float opticalDepthPA = 0;
// Численное интегрирование для вычисления
// влияния света каждой точки P в AB
float time = tA;
float ds = (tB-tA) / (float)(_ViewSamples);
for (int i = 0; i < _ViewSamples; i ++)
 // Позиция точки
 // (сэмплируется посередине сэмпла отрезка видимости)
 float3 P = O + D * (time + viewSampleSize*0.5);
 // Оптическая толщина текущего отрезка
 // ρ(h) * ds
 float height = distance(C, P) - _PlanetRadius;
 float opticalDepthSegment = exp(-height / _ScaleHeight) * ds;
 // Накопление оптических толщин
 // D(PA)
 opticalDepthPA += opticalDepthSegment;
 time += ds;

Light sampling

If we return to the expression for the influence of light $P$, then the only necessary value is the optical thickness 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 thickness of the segment $\overline{CP}$into a function called lightSampling. The name is taken from a ray of light , which is a segment starting at$P$and directed towards the sun. We called the point at which it comes out of the atmosphere$C$.

However, the function lightSamplingdoes not just calculate the optical thickness$\overline{CP}$. So far we have considered only the influence of the atmosphere and ignored the role of the planet itself. Our equations do not take into account the possibility that a ray of light moving from$P$to the sun, may collide with the planet. If this happens, then all calculations performed up to this point are not applicable, because the light does not actually reach the camera.

In the diagram below it is easy to see that the influence of light $P_0$ must be ignored because no light of the sun reaches $P_0$. When cycling through points between$P$ and $C$The function lightSamplingalso checks for collisions with the planet. This can be done by checking the height of the point for negativity.

bool lightSampling
( float3 P, // текущая точка внутри атмосферной сферы
 float3 S, // Направление к солнцу
 out float opticalDepthCA
 float _; // это нас не интересует
 float C;
 rayInstersect(P, S, _PlanetCentre, _AtmosphereRadius, _, C);
 // Сэмплы на отрезке 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;
 // Внутри планеты
 if (height < 0)
 return false;
 // Оптическая толщина для луча света
 opticalDepthCA += exp(-height / _RayScaleHeight) * ds;
 time += ds;
 return true;

The function above calculates the point first $C$with the help of rayInstersect. Then she divides the segment$\overline{PA}$into segments _LightSamplesof length ds. The calculation of optical thickness is the same as that used in the outermost loop.

The function returns false if a collision with a planet occurs. We can use this to update the missing code of the outermost loop by replacing ....

 // D(CP)
 float opticalDepthCP = 0;
 bool overground = lightSampling(P, S);
 if (overground)
 // Комбинированное пропускание
 // T(CP) * T(PA) = T(CPA) = exp{ -β(λ) [D(CP) + D(PA)]}
 float transmittance = exp
 -_ScatteringCoefficient *
 (opticalDepthCP + opticalDepthPA)
 // Влияния света
 // T(CPA) * ρ(h) * ds
 totalViewSamples += transmittance * opticalDepthSegment;

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

[Note per .: on the author’s Patreon page, you can buy access to Standard and Premium versions of the finished shader.]

Also popular now: