GLSL Atmospheric Scattering Not Scaling With Transformations

Question

I am trying to implement atmospheric scatting in GLSL version 4.10. I am adapting the shaders from the this Shadertoy shader https://www.shadertoy.com/view/lslXDr. The atmosphere in my program is created from a scaled version of the planet sphere.

I have the actual scattering equations working, but the inner radius of the atmosphere does not line up with the outer radius of the sphere for most camera positions. I know this is from the radius of the atmosphere being bigger than the planet sphere, but I cannot seem to get it to scale right.

My problem is best illustrated here. The model is scaled up in these pictures. As can be seen, the atmosphere inner radius does not match the radius of the planet (the dark blue sphere).

Here the model is scaled and translated. The atmosphere is off center from the camera and the inner atmosphere is still not lined up with the planet.

Here is the vertex shader, which is essentially a pass through shader

#version 410

in vec4 vPosition;
in vec3 vNormal;

out vec3 fPosition;
out mat3 m;

uniform mat4 model;
uniform mat4 view;
uniform mat4 projection;

void main()
{   
    fPosition = vec3(vPosition);
    m = mat3(model);

    gl_Position = projection*view*model*vPosition;
}

And the fragment shader.

#version 410

uniform float time;
uniform vec3 camPosition;
uniform float fInnerRadius;
uniform float fOuterRadius;

in vec3 fPosition;
in mat3 m;

out vec4 FragColor;

const float PI = 3.14159265359;
const float degToRad = PI / 180.0;
const float MAX = 10000.0;

float K_R = 0.166;
const float K_M = 0.0025;
const float E = 14.3;
const vec3 C_R = vec3(0.3, 0.7, 1.0);
const float G_M = -0.85;

float SCALE_H = 4.0 / (fOuterRadius - fInnerRadius);
float SCALE_L = 1.0 / (fOuterRadius - fInnerRadius);

const int numOutScatter = 10;
const float fNumOutScatter = 10.0;
const int numInScatter = 10;
const float fNumInScatter = 10.0;

vec3 rayDirection(vec3 camPosition) {
    vec3 ray = m*fPosition - camPosition;
    float far = length(ray);
    return ray /= far;
}

vec2 rayIntersection(vec3 p, vec3 dir, float radius ) {
    float b = dot( p, dir );
    float c = dot( p, p ) - radius * radius;

    float d = b * b - c;
    if ( d < 0.0 ) {
        return vec2( MAX, -MAX );
    }
    d = sqrt( d );

    float near = -b - d;
    float far = -b + d;

    return vec2(near, far);
 }

// Mie
// g : ( -0.75, -0.999 )
//      3 * ( 1 - g^2 )               1 + c^2
// F = ----------------- * -------------------------------
//      2 * ( 2 + g^2 )     ( 1 + g^2 - 2 * g * c )^(3/2)
float miePhase( float g, float c, float cc ) {
    float gg = g * g;

    float a = ( 1.0 - gg ) * ( 1.0 + cc );

    float b = 1.0 + gg - 2.0 * g * c;
    b *= sqrt( b );
    b *= 2.0 + gg;  

    return 1.5 * a / b;
}

// Reyleigh
// g : 0
// F = 3/4 * ( 1 + c^2 )
float rayleighPhase( float cc ) {
    return 0.75 * ( 1.0 + cc );
}


float density(vec3 p) {
    return exp(-(length(p) - fInnerRadius) * SCALE_H);
}

float optic(vec3 p, vec3 q) {
    vec3 step = (q - p) / fNumOutScatter;
    vec3 v = p + step * 0.5;

    float sum = 0.0;
    for(int i = 0; i < numOutScatter; i++) {
        sum += density(v);
        v += step;
    }
    sum *= length(step)*SCALE_L;
    return sum;
}

vec3 inScatter(vec3 o, vec3 dir, vec2 e, vec3 l) {
    float len = (e.y - e.x) / fNumInScatter;
    vec3 step = dir * len;
    vec3 p = o + dir * e.x;
    vec3 v = p + dir * (len * 0.5);

    vec3 sum = vec3(0.0);
    for(int i = 0; i < numInScatter; i++) {
        vec2 f = rayIntersection(v, l, fOuterRadius);
        vec3 u = v + l * f.y;
        float n = (optic(p, v) + optic(v, u))*(PI * 4.0);
        sum += density(v)* exp(-n * ( K_R * C_R + K_M ));
        v += step;
    }
    sum *= len * SCALE_L;
    float c = dot(dir, -l);
    float cc = c * c;
    return sum * ( K_R * C_R * rayleighPhase( cc ) + K_M * miePhase( G_M, c, cc ) ) * E;
}

void main (void)
{
    vec3 dir = rayDirection(vec3(camPosition.x, 0.0, camPosition.z));
    vec3 eye = vec3(camPosition.x, 0.0, camPosition.z);

    vec3 l = normalize(vec3(0.0, 0.0, 1.0));

    vec2 e = rayIntersection(eye, dir, fOuterRadius);
    if ( e.x > e.y ) {
        discard;
    }
    vec2 f = rayIntersection(eye, dir, fInnerRadius);
    e.y = min(e.y, f.x);

    vec3 I = inScatter(eye, dir, e, l);

    FragColor = vec4(I, 1.0);
}

If needed here is the code that draws the atmosphere. The code that draws the planet has essentially the same transformations sans the scaleFactor.

void drawAtmosphere()
{
    glUseProgram(atmosphereShader);

    v = getViewMatrix();
    vec3 Position = getCameraPosition();

    float scaleFactor = 1.25;
    m = multiplymat4(translate(0.0, 0.0, -10), scale(fScale*scaleFactor));
    float fOuter = (fScale*scaleFactor);
    float fInner = fScale;

    glUniform1f(glGetUniformLocation(atmosphereShader, "fInnerRadius"), fInner);
    glUniform1f(glGetUniformLocation(atmosphereShader, "fOuterRadius"), fOuter);
    glUniform3f(glGetUniformLocation(atmosphereShader, "camPosition"), Position.x, Position.y, Position.z);
    glUniform1f(glGetUniformLocation(atmosphereShader, "time"), glfwGetTime());

    initMVP(atmosphereShader, m, v);

    glBindVertexArray (atmosphereVAO);
    glDrawArrays( GL_TRIANGLES, 0, planet.vertexNumber);
    glBindVertexArray(0);

}

Any help, or anything that can point me in the right direction is appreciated.

Answer 1

Found the problem was caused by incorrect calculation of the camera position and not taking into account the model space of the object. I uploaded a stripped down version of the code here.

Hopefully this will help anyone trying to implement Sean O'Neil's atmosphere code.