<p>I am currently experimenting with some atmospheric scattering in glsl in my OpenGL engine.
I am experiencing a strange phenomena where a strange band / separation artifact appears along the middle of the screen:</p>
<p><a href="https://i.sstatic.net/sovGD.png" rel="nofollow noreferrer">Artifact from low altitude</a></p>
<p><a href="https://i.sstatic.net/3FcnL.png" rel="nofollow noreferrer">Artifact from mid altitude</a></p>
<p><a href="https://i.sstatic.net/w7uJD.png" rel="nofollow noreferrer">Artifact from high altitude</a></p>
<p>Here is the code responsible for the actual scattering:</p>
<pre><code>#define PI 3.1415926535897932384626433832795
#define iSteps 16
#define jSteps 8

vec2 rsi(vec3 r0, vec3 rd, float sr) {
    // ray-sphere intersection that assumes
    // the sphere is centered at the origin.
    // No intersection when result.x &gt; result.y
    float a = dot(rd, rd);
    float b = 2.0 * dot(rd, r0);
    float c = dot(r0, r0) - (sr * sr);
    float d = (b*b) - 4.0*a*c;
    if (d &lt; 0.0) return vec2(1e5,-1e5);
    return vec2(
        (-b - sqrt(d))/(2.0*a),
        (-b + sqrt(d))/(2.0*a)
    );
}

vec3 atmosphere(vec3 r, vec3 r0, vec3 pSun, float iSun, float rPlanet, float rAtmos, vec3 kRlh, float kMie, float shRlh, float shMie, float g) {
    // Normalize the sun and view directions.
    pSun = normalize(pSun);
    r = normalize(r);

    // Calculate the step size of the primary ray.
    vec2 p = rsi(r0, r, rAtmos);
    if (p.x &gt; p.y) return vec3(0,0,0);
    p.y = min(p.y, rsi(r0, r, rPlanet).x);
    float iStepSize = (p.y - p.x) / float(iSteps);

    // Initialize the primary ray time.
    float iTime = 0.0;

    // Initialize accumulators for Rayleigh and Mie scattering.
    vec3 totalRlh = vec3(0,0,0);
    vec3 totalMie = vec3(0,0,0);

    // Initialize optical depth accumulators for the primary ray.
    float iOdRlh = 0.0;
    float iOdMie = 0.0;

    // Calculate the Rayleigh and Mie phases.
    float mu = dot(r, pSun);
    float mumu = mu * mu;
    float gg = g * g;
    float pRlh = 3.0 / (16.0 * PI) * (1.0 + mumu);
    float pMie = 3.0 / (8.0 * PI) * ((1.0 - gg) * (mumu + 1.0)) / (pow(1.0 + gg - 2.0 * mu * g, 1.5) * (2.0 + gg));

    // Sample the primary ray.
    for (int i = 0; i &lt; iSteps; i++) {

        // Calculate the primary ray sample position.
        vec3 iPos = r0 + r * (iTime + iStepSize * 0.5);

        // Calculate the height of the sample.
        float iHeight = length(iPos) - rPlanet;

        // Calculate the optical depth of the Rayleigh and Mie scattering for this step.
        float odStepRlh = exp(-iHeight / shRlh) * iStepSize;
        float odStepMie = exp(-iHeight / shMie) * iStepSize;

        // Accumulate optical depth.
        iOdRlh += odStepRlh;
        iOdMie += odStepMie;

        // Calculate the step size of the secondary ray.
        float jStepSize = rsi(iPos, pSun, rAtmos).y / float(jSteps);

        // Initialize the secondary ray time.
        float jTime = 0.0;

        // Initialize optical depth accumulators for the secondary ray.
        float jOdRlh = 0.0;
        float jOdMie = 0.0;

        
        // Sample the secondary ray.
        for (int j = 0; j &lt; jSteps; j++) {

            // Calculate the secondary ray sample position.
            vec3 jPos = iPos + pSun * (jTime + jStepSize * 0.5);

            // Calculate the height of the sample.
            float jHeight = length(jPos) - rPlanet;

            // Accumulate the optical depth.
            jOdRlh += exp(-jHeight / shRlh) * jStepSize;
            jOdMie += exp(-jHeight / shMie) * jStepSize;

            // Increment the secondary ray time.
            jTime += jStepSize;
        }

        // Calculate attenuation.
        vec3 attn = exp(-(kMie * (iOdMie + jOdMie) + kRlh * (iOdRlh + jOdRlh)));

        // Accumulate scattering.
        totalRlh += odStepRlh * attn;
        totalMie += odStepMie * attn;

        // Increment the primary ray time.
        iTime += iStepSize;
    }

    // Calculate and return the final color.
    return iSun * (pRlh * kRlh * totalRlh + pMie * kMie * totalMie);
}
</code></pre>
<p>The function is called as so per fragment in the main function:</p>
<pre><code>vec3 color = atmosphere(
        normalize(position),          
        vec3(0.0, 6372e3 + camPos.y, 0.0), 
        vec3(0.0, 2.0, -1.0),              
        48.0,                           
        6371e3,                         
        6471e3,                         
        vec3(5.5e-6, 13.0e-6, 22.4e-6), 
        21e-6,                          
        8e3,                            
        1.2e3,                          
        0.758                           
    );

    // Apply exposure.
    color = 1.0 - exp(-1.0 * color);
</code></pre>
<p>What is causing these artifacts? Or is this intended behaviour?</p>

CC BY-SA 4.0
2020-12-04T18:15:04.970Z
If you mean the band near the horizon that is intended. The ray of light coming from the horizon travels a longer distance and so more of its blue light is scattered and the red light is left.
ali
CC BY-SA 4.0
2020-12-04T18:15:16.137Z
Have a look at here https://www.scratchapixel.com/lessons/procedural-generation-virtual-worlds/simulating-sky/simulating-colors-of-the-sky
ali
CC BY-SA 4.0
2020-12-04T19:15:43.980Z
From the screenshots it looks like the artifact may be located at an altitude equal to the camera. There may be something in the math that is depending on the sign of of the ray direction (whether it's traveling up or down relative to the planet) or something that's numerically unstable there. I would try to isolate where the artifact is coming from by commenting out different parts of the code / different lighting components, and check whether the artifact is still produced. That should help you narrow down what part of the math is creating it.
Nathan Reed
CC BY-SA 4.0
2020-12-04T22:28:26.153Z
@NathanReed that is actually a very good observation, thanks! I will start debugging the shader. If you get anymore insights please let me know. Thanks again :P
user10037795
<opengl><3d><glsl><c++><atmospherics>
Atmospheric scattering band/seperation artifact
user10037795
# Atmospheric scattering band/seperation artifact

(Asked by: username_0 on 2020-12-04)

I am currently experimenting with some atmospheric scattering in glsl in my OpenGL engine.
I am experiencing a strange phenomena where a strange band / separation artifact appears along the middle of the screen:

[Artifact from low altitude](https://i.sstatic.net/sovGD.png)

[Artifact from mid altitude](https://i.sstatic.net/3FcnL.png)

[Artifact from high altitude](https://i.sstatic.net/w7uJD.png)

Here is the code responsible for the actual scattering:

```
#define PI 3.1415926535897932384626433832795
#define iSteps 16
#define jSteps 8

vec2 rsi(vec3 r0, vec3 rd, float sr) {
    // ray-sphere intersection that assumes
    // the sphere is centered at the origin.
    // No intersection when result.x > result.y
    float a = dot(rd, rd);
    float b = 2.0 * dot(rd, r0);
    float c = dot(r0, r0) - (sr * sr);
    float d = (b*b) - 4.0*a*c;
    if (d < 0.0) return vec2(1e5,-1e5);
    return vec2(
        (-b - sqrt(d))/(2.0*a),
        (-b + sqrt(d))/(2.0*a)
    );
}

vec3 atmosphere(vec3 r, vec3 r0, vec3 pSun, float iSun, float rPlanet, float rAtmos, vec3 kRlh, float kMie, float shRlh, float shMie, float g) {
    // Normalize the sun and view directions.
    pSun = normalize(pSun);
    r = normalize(r);

    // Calculate the step size of the primary ray.
    vec2 p = rsi(r0, r, rAtmos);
    if (p.x > p.y) return vec3(0,0,0);
    p.y = min(p.y, rsi(r0, r, rPlanet).x);
    float iStepSize = (p.y - p.x) / float(iSteps);

    // Initialize the primary ray time.
    float iTime = 0.0;

    // Initialize accumulators for Rayleigh and Mie scattering.
    vec3 totalRlh = vec3(0,0,0);
    vec3 totalMie = vec3(0,0,0);

    // Initialize optical depth accumulators for the primary ray.
    float iOdRlh = 0.0;
    float iOdMie = 0.0;

    // Calculate the Rayleigh and Mie phases.
    float mu = dot(r, pSun);
    float mumu = mu * mu;
    float gg = g * g;
    float pRlh = 3.0 / (16.0 * PI) * (1.0 + mumu);
    float pMie = 3.0 / (8.0 * PI) * ((1.0 - gg) * (mumu + 1.0)) / (pow(1.0 + gg - 2.0 * mu * g, 1.5) * (2.0 + gg));

    // Sample the primary ray.
    for (int i = 0; i < iSteps; i++) {

        // Calculate the primary ray sample position.
        vec3 iPos = r0 + r * (iTime + iStepSize * 0.5);

        // Calculate the height of the sample.
        float iHeight = length(iPos) - rPlanet;

        // Calculate the optical depth of the Rayleigh and Mie scattering for this step.
        float odStepRlh = exp(-iHeight / shRlh) * iStepSize;
        float odStepMie = exp(-iHeight / shMie) * iStepSize;

        // Accumulate optical depth.
        iOdRlh += odStepRlh;
        iOdMie += odStepMie;

        // Calculate the step size of the secondary ray.
        float jStepSize = rsi(iPos, pSun, rAtmos).y / float(jSteps);

        // Initialize the secondary ray time.
        float jTime = 0.0;

        // Initialize optical depth accumulators for the secondary ray.
        float jOdRlh = 0.0;
        float jOdMie = 0.0;

        
        // Sample the secondary ray.
        for (int j = 0; j < jSteps; j++) {

            // Calculate the secondary ray sample position.
            vec3 jPos = iPos + pSun * (jTime + jStepSize * 0.5);

            // Calculate the height of the sample.
            float jHeight = length(jPos) - rPlanet;

            // Accumulate the optical depth.
            jOdRlh += exp(-jHeight / shRlh) * jStepSize;
            jOdMie += exp(-jHeight / shMie) * jStepSize;

            // Increment the secondary ray time.
            jTime += jStepSize;
        }

        // Calculate attenuation.
        vec3 attn = exp(-(kMie * (iOdMie + jOdMie) + kRlh * (iOdRlh + jOdRlh)));

        // Accumulate scattering.
        totalRlh += odStepRlh * attn;
        totalMie += odStepMie * attn;

        // Increment the primary ray time.
        iTime += iStepSize;
    }

    // Calculate and return the final color.
    return iSun * (pRlh * kRlh * totalRlh + pMie * kMie * totalMie);
}

```

The function is called as so per fragment in the main function:

```
vec3 color = atmosphere(
        normalize(position),          
        vec3(0.0, 6372e3 + camPos.y, 0.0), 
        vec3(0.0, 2.0, -1.0),              
        48.0,                           
        6371e3,                         
        6471e3,                         
        vec3(5.5e-6, 13.0e-6, 22.4e-6), 
        21e-6,                          
        8e3,                            
        1.2e3,                          
        0.758                           
    );

    // Apply exposure.
    color = 1.0 - exp(-1.0 * color);

```

What is causing these artifacts? Or is this intended behaviour?

**Tags:** opengl, 3d, glsl, c++, atmospherics

> **Comments:**
> **username_1**: If you mean the band near the horizon that is intended. The ray of light coming from the horizon travels a longer distance and so more of its blue light is scattered and the red light is left.  
> **username_1**: Have a look at here https://www.scratchapixel.com/lessons/procedural-generation-virtual-worlds/simulating-sky/simulating-colors-of-the-sky  
> **username_2**: From the screenshots it looks like the artifact may be located at an altitude equal to the camera. There may be something in the math that is depending on the sign of of the ray direction (whether it's traveling up or down relative to the planet) or something that's numerically unstable there. I would try to isolate where the artifact is coming from by commenting out different parts of the code / different lighting components, and check whether the artifact is still produced. That should help you narrow down what part of the math is creating it.  
> **username_0**: @username_2 that is actually a very good observation, thanks! I will start debugging the shader. If you get anymore insights please let me know. Thanks again :P  

---

## 0 answers

