Special topics in three-dimensional engine design-atmospheric scattering effects

When doing GIS, there will be a atmospheric circle around the earth, atmospheric scattering, this aspect of the algorithm is the field of computer graphics in-depth research, but there are several mature scattering algorithms. I borrowed from <<gpu 2. High-performance graphics chips and general-purpose computer Programming Tips >> The 16th chapter of the algorithm, the realization of an atmospheric scattering. Effect.

The blue sky in the picture is the effect of scattering, please check the algorithm of the book by yourself.

Steps:

1: Creates an ellipsoid, generates vertices, and an array of vertex indices. This ellipsoid generation algorithm is subsequently posted.

2: Pass uniform According to algorithm, run shader

3: About the image, the elevation of the processing, follow-up post.

Shader Source code:

Vertex:

Uniform VEC3 U_VIEWERPOSITIONWC;
Uniform VEC3 U_SUNDIRECTIONWC;

const float U_PI = 3.141592653589793;
Uniform mat4 u_modelviewprojection;
Uniform mat4 U_modelview;


#define Sky_from_atmosphere

attribute VEC4 position;

Uniform float fcameraheight;
Uniform float fCameraHeight2;
Uniform float Fouterradius; The outer (atmosphere) Radius
Uniform float fOuterRadius2; Fouterradius^2
Uniform float Finnerradius; The inner (planetary) Radius
Uniform float Fscale; 1/(Fouterradius-finnerradius)
Uniform float fscaledepth; The scale depth (i.e. the altitude at which the atmosphere ' s average density is found)
Uniform float fscaleoverscaledepth; Fscale/fscaledepth

const float KR = 0.0025;
const FLOAT FKR4PI = Kr * 4.0 * U_PI;
Const-Float Km = 0.0015;
const FLOAT FKM4PI = Km * 4.0 * U_PI;
Const float ESun = 15.0;
const FLOAT Fkmesun = Km * ESUN;
const FLOAT Fkresun = Kr * ESUN;
Const VEC3 V3invwavelength = VEC3 (
5.60204474633241,//Red = 1.0/math.pow (0.650, 4.0)
9.473284437923038,//Green = 1.0/math.pow (0.570, 4.0)
19.643802610477206); Blue = 1.0/math.pow (0.475, 4.0)
Const float rayleighscaledepth = 0.25;

const int nsamples = 2;
Const FLOAT Fsamples = 2.0;

Varying VEC3 v_rayleighcolor;
Varying VEC3 v_miecolor;
Varying VEC3 V_tocamera;
Varying VEC3 V_positionec;

Float scale (float fcos)
{
float x = 1.0-fcos;
Return fscaledepth * EXP ( -0.00287 + x* (0.459 + x* (3.83 + x* ( -6.80 + x*5.25)));
}

void Main (void)
{
Get the ray from the camera to the vertex and its length (which are the far point of the ray passing through the ATMOSPH Ere
VEC3 v3pos = position.xyz;
VEC3 V3ray = V3POS-U_VIEWERPOSITIONWC;
float Ffar = Length (V3ray);
V3ray/= Ffar;

#ifdef Sky_from_space
Calculate the closest intersection of the ray with the outer atmosphere (which are the the near point of the ray passing THR Ough the atmosphere)
Float B = 2.0 * DOT (U_VIEWERPOSITIONWC, v3ray);
float C = fcameraheight2-fouterradius2;
float Fdet = max (0.0, b*b-4.0 * C);
float Fnear = 0.5 * (-b-sqrt (fdet));

Calculate the ray ' s starting position, then Calculate its scattering offset
VEC3 V3start = U_VIEWERPOSITIONWC + V3ray * fnear;
Ffar-= fnear;
float Fstartangle = dot (V3ray, v3start)/Fouterradius;
float fstartdepth = exp ( -1.0/fscaledepth);
float Fstartoffset = Fstartdepth*scale (Fstartangle);
#else//Sky_from_atmosphere
Calculate the ray ' s starting position, then Calculate its scattering offset
VEC3 V3start = U_VIEWERPOSITIONWC;
float fheight = Length (V3start);
float fdepth = exp (fscaleoverscaledepth * (finnerradius-fcameraheight));
float Fstartangle = dot (V3ray, v3start)/fheight;
float Fstartoffset = Fdepth*scale (Fstartangle);
#endif

Initialize The scattering loop variables
float fsamplelength = ffar/fsamples;
float fscaledlength = fsamplelength * Fscale;
VEC3 V3sampleray = V3ray * FSAMPLELENGTH;
VEC3 v3samplepoint = v3start + V3sampleray * 0.5;

Now loop through the sample rays
VEC3 V3frontcolor = vec3 (0.0, 0.0, 0.0);
for (int i=0; i<nsamples; i++)
{
float fheight = Length (v3samplepoint);
float fdepth = exp (fscaleoverscaledepth * (finnerradius-fheight));
VEC3 lightposition = normalize (U_VIEWERPOSITIONWC); U_sundirectionwc
float Flightangle = dot (lightposition, v3samplepoint)/fheight;
float Fcameraangle = dot (V3ray, v3samplepoint)/fheight;
float Fscatter = (fstartoffset + fdepth* (Scale (Flightangle)-scale (Fcameraangle)));
VEC3 v3attenuate = exp (-fscatter * (v3invwavelength * fkr4pi + fkm4pi));
V3frontcolor + = V3attenuate * (fdepth * fscaledlength);
V3samplepoint + = V3sampleray;
}

Finally, scale the Mie and Rayleigh colors and set up the varying variables for the pixel shader
V_miecolor = V3frontcolor * FKMESUN;
V_rayleighcolor = V3frontcolor * (v3invwavelength * Fkresun);
V_tocamera = U_viewerpositionwc-v3pos;
V_positionec = (U_modelview * position). xyz;
gl_position = u_modelviewprojection * Position;
}

Fragment:

#ifdef Gl_fragment_precision_high
Precision HIGHP float;
#else
Precision Mediump float;
#endif

const float u_infinity = 5906376272000.0;


struct u_raysegment
{
float start;
float stop;
};

Const U_raysegment u_emptyraysegment = u_raysegment (-u_infinity,-u_infinity);

Const U_raysegment u_fullraysegment = u_raysegment (0.0, u_infinity);

struct U_ray
{
VEC3 origin;
VEC3 direction;
};
Uniform mat4 U_inversemodelview;
struct u_ellipsoid
{
VEC3 Center;
VEC3 radii;
VEC3 inverseradii;
VEC3 inverseradiisquared;
};
Uniform mat4 U_view;
Uniform VEC3 U_SUNDIRECTIONWC;
U_raysegment U_rayellipsoidintersectioninterval (U_ray Ray, u_ellipsoid ellipsoid)
{
Ray and ellipsoid Center in the eye coordinates. Radii in model coordinates.
VEC3 q = ellipsoid.inverseradii * (U_inversemodelview * VEC4 (Ray.origin, 1.0)). xyz;
VEC3 w = ellipsoid.inverseradii * (U_inversemodelview * VEC4 (ray.direction, 0.0)). xyz;

Q = q-ellipsoid.inverseradii * (U_inversemodelview * VEC4 (Ellipsoid.center, 1.0)). xyz;

float q2 = dot (q, q);
float QW = dot (q, W);

if (Q2 > 1.0)//Outside ellipsoid.
{
if (qw >= 0.0)//looking outward or tangent (0 intersections).
{
return u_emptyraysegment;
}
else//QW < 0.0.
{
float QW2 = QW * QW;
float difference = q2-1.0; Positively valued.
float W2 = Dot (w, W);
FLOAT Product = W2 * difference;

if (QW2 < product)//imaginary roots (0 intersections).
{
return u_emptyraysegment;
}
else if (qw2 > Product)//Distinct roots (2 intersections).
{
float discriminant = QW * QW-PRODUCT;
Float temp =-qw + sqrt (discriminant); Avoid cancellation.
float root0 = temp/w2;
float root1 = difference/temp;
if (Root0 < ROOT1)
{
U_raysegment i = u_raysegment (root0, ROOT1);
return i;
}
Else
{
U_raysegment i = u_raysegment (ROOT1, root0);
return i;
}
}
else//Qw2 = = product. Repeated roots (2 intersections).
{
float root = sqrt (DIFFERENCE/W2);
U_raysegment i = u_raysegment (root, root);
return i;
}
}
}
else if (Q2 < 1.0)//Inside ellipsoid (2 intersections).
{
float difference = q2-1.0; negatively valued.
float W2 = Dot (w, W);
FLOAT Product = W2 * difference; negatively valued.
float discriminant = QW * QW-PRODUCT;
Float temp =-qw + sqrt (discriminant); Positively valued.
U_raysegment i = u_raysegment (0.0, TEMP/W2);
return i;
}
else//q2 = = 1.0. On ellipsoid.
{
if (QW < 0.0)//looking inward.
{
float W2 = Dot (w, W);
U_raysegment i = u_raysegment (0.0,-QW/W2);
return i;
}
else//QW >= 0.0. looking outward or tangent.
{
return u_emptyraysegment;
}
}
}

Uniform float u_morphtime;

Float u_luminance (VEC3 RGB)
{
Algorithm from Chapter to Graphics Shaders.
Const VEC3 W = VEC3 (0.2125, 0.7154, 0.0721);
return dot (RGB, W);
}


BOOL U_isempty (u_raysegment interval)
{
Return (Interval.stop < 0.0);
}


U_ellipsoid U_getwgs84ellipsoidec ()
{
VEC3 radii = VEC3 (6378137.0, 6378137.0, 6356752.314245);
VEC3 inverseradii = VEC3 (1.0/radii.x, 1.0/RADII.Y, 1.0/radii.z);
VEC3 inverseradiisquared = inverseradii * INVERSERADII;
U_ellipsoid temp = u_ellipsoid (u_view[3].xyz, radii, inverseradii, inverseradiisquared);
return temp;
}




const FLOAT g =-0.95;
const FLOAT g2 = g * g;

Varying VEC3 v_rayleighcolor;
Varying VEC3 v_miecolor;
Varying VEC3 V_tocamera;
Varying VEC3 V_positionec;

void Main (void)
{
Todo:make arbitrary ellipsoid
U_ellipsoid ellipsoid = U_getwgs84ellipsoidec ();

VEC3 Direction = normalize (V_positionec);
U_ray ray = U_ray (VEC3 (0.0), direction);

U_raysegment intersection = U_rayellipsoidintersectioninterval (ray, ellipsoid);
if (!u_isempty (intersection)) {
Discard
}

Extra normalize added for Android
float Fcos = dot (U_SUNDIRECTIONWC, normalize (V_tocamera))/Length (V_tocamera);
float frayleighphase = 0.75 * (1.0 + Fcos*fcos);
float Fmiephase = 1.5 * ((1.0-G2)/(2.0 + G2)) * (1.0 + Fcos*fcos)/POW (1.0 + g2-2.0*g*fcos, 1.5);

Const FLOAT Fexposure = 2.0;

VEC3 RGB = frayleighphase * v_rayleighcolor + fmiephase * v_miecolor;
RGB = VEC3 (1.0)-exp (-fexposure * RGB);
float L = u_luminance (RGB);
Gl_fragcolor = VEC4 (RGB, MIN (Smoothstep (0.0, 0.1, L), 1.0) * Smoothstep (0.0, 1.0, u_morphtime));
}

Special topics in three-dimensional engine design-atmospheric scattering effects