The next evolution of GPU programming

 

 

 

 

Normal

 

Rayleigh Scattering

 

Mie Scattering

Aerial Perspective

Author:

David W. Leonard

Date:

29 May 2003

Details:

This project implements Aerial Perspective as described in Rendering Outdoor Light Scattering in Real Time by Hoffman and Preetham.  My version uses Cg shaders with OpenGL, instead of ATI shaders with DirectX.  In addition, the scattering formulae were partitioned between vertex shader and OpenGL code for performance; minimizing total number of computations by maximizing per-frame variable calculations within OpenGL’s display function.  The revised formulae are as follows.

 

Extinction Factor: Fex(s)= e-K1*s, where K1= βRayleigh+βMie and s is vertex distance.

Total In-scattering: Lin(s, θ)= (βRayleigh(θ)+ βMie(θ))*K2, where K2=ESun/K1, and

                                                                                          θ is the Sun’s zenith angle.

Rayleigh Scattering: βRayleigh(θ) =  K3*(1+cos2θ), where K3=3*βRayleigh/(16*π).

Mie Scattering: βMie(θ) =  K4/(K5.x-K5.y*cos θ)1.5, where K4=(K5.x+K5.y)*βMie/(4*π) and

                                                                                   K5=[1+g2, 2*g].

A fragment shader is necessary because current graphics APIs only support monochromatic fog and Aerial Perspective requires tri-chromatic at the least.  This shader’s functionality is L(s, θ)= L0*Fex(s) + Lin(s, θ)*(1-Fex(s)).


Source
Aerial Perspective was implemented on top of NVIDIA’s Cg Procedural Terrain found in the Cg Toolkit.

Files:

cg_terrain.cg, pAerialPersp.cg


Download Associated Files:
      AerialPersp-bin.zip (3,504 KB)

 

 

/*********************************************************************NVMH3****

Path:  NVSDK\Common\media\programs

File:  cg_terrain.cg

 

Copyright NVIDIA Corporation 2002

TO THE MAXIMUM EXTENT PERMITTED BY APPLICABLE LAW, THIS SOFTWARE IS PROVIDED

*AS IS* AND NVIDIA AND ITS SUPPLIERS DISCLAIM ALL WARRANTIES, EITHER EXPRESS

OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, IMPLIED WARRANTIES OF MERCHANTABILITY

AND FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT SHALL NVIDIA OR ITS SUPPLIERS

BE LIABLE FOR ANY SPECIAL, INCIDENTAL, INDIRECT, OR CONSEQUENTIAL DAMAGES

WHATSOEVER (INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS PROFITS,

BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR ANY OTHER PECUNIARY LOSS)

ARISING OUT OF THE USE OF OR INABILITY TO USE THIS SOFTWARE, EVEN IF NVIDIA HAS

BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.

 

 

Comments:

 

  CG noise implementation for vertex program profile

  sgreen 5/02/02

 

  Compile with: cgc.exe -profile vp20 vnoise.cg

 

  This is based on Perlin's original code:

  http://mrl.nyu.edu/~perlin/doc/oscar.html

 

  It combines the permutation and gradient tables into one array of

  float4's to conserve constant memory.

  The table is duplicated twice to avoid modulo operations.

 

  Notes:

  Should use separate tables for 1, 2 and 3D versions

 

******************************************************************************/

 

#define B  32      // table size

#define B2 66      // B*2 + 2

#define BR 0.03125 // 1 / B

 

// this is the smoothstep function f(t) = 3t^2 - 2t^3, without the normalization

float3 s_curve(float3 t)

{

    return t*t*( float3(3.0f, 3.0f, 3.0f) - float3(2.0f, 2.0f, 2.0f)*t);

}

 

float2 s_curve(float2 t)

{

    return t*t*( float2(3.0f, 3.0f) - float2(2.0f, 2.0f)*t);

}

 

float s_curve(float t)

{

    return t*t*(3.0f-2.0f*t);

}

 

// 3D version

float noise(float3 v, const uniform float4 pg[B2])

{

    v = v + float3(10000.0f, 10000.0f, 10000.0f);   // hack to avoid negative numbers

   

    float3 i = frac(v * BR) * B;   // index between 0 and B-1

    float3 f = frac(v);            // fractional position

   

    // lookup in permutation table

    float2 p;

    p[0] = pg[ i[0]     ].w;

    p[1] = pg[ i[0] + 1 ].w;

    p = p + i[1];

   

    float4 b;

    b[0] = pg[ p[0] ].w;

    b[1] = pg[ p[1] ].w;

    b[2] = pg[ p[0] + 1 ].w;

    b[3] = pg[ p[1] + 1 ].w;

    b = b + i[2];

   

    // compute dot products between gradients and vectors

    float4 r;

    r[0] = dot( pg[ b[0] ].xyz, f );

    r[1] = dot( pg[ b[1] ].xyz, f - float3(1.0f, 0.0f, 0.0f) );

    r[2] = dot( pg[ b[2] ].xyz, f - float3(0.0f, 1.0f, 0.0f) );

    r[3] = dot( pg[ b[3] ].xyz, f - float3(1.0f, 1.0f, 0.0f) );

   

    float4 r1;

    r1[0] = dot( pg[ b[0] + 1 ].xyz, f - float3(0.0f, 0.0f, 1.0f) );

    r1[1] = dot( pg[ b[1] + 1 ].xyz, f - float3(1.0f, 0.0f, 1.0f) );

    r1[2] = dot( pg[ b[2] + 1 ].xyz, f - float3(0.0f, 1.0f, 1.0f) );

    r1[3] = dot( pg[ b[3] + 1 ].xyz, f - float3(1.0f, 1.0f, 1.0f) );

   

    // interpolate

    f = s_curve(f);

    r = lerp( r, r1, f[2] );

    r = lerp( r.xyyy, r.zwww, f[1] );

    return lerp( r.x, r.y, f[0] );

}

 

// 2D version

float noise(float2 v, const uniform float4 pg[B2])

{

    v = v + float2(10000.0f, 10000.0f);

   

    float2 i = frac(v * BR) * B;   // index between 0 and B-1

    float2 f = frac(v);            // fractional position

   

    // lookup in permutation table

    float2 p;

    p[0] = pg[ i[0]   ].w;

    p[1] = pg[ i[0]+1 ].w;

    p = p + i[1];

   

    // compute dot products between gradients and vectors

    float4 r;

    r[0] = dot( pg[ p[0] ].xy,   f);

    r[1] = dot( pg[ p[1] ].xy,   f - float2(1.0f, 0.0f) );

    r[2] = dot( pg[ p[0]+1 ].xy, f - float2(0.0f, 1.0f) );

    r[3] = dot( pg[ p[1]+1 ].xy, f - float2(1.0f, 1.0f) );

   

    // interpolate

    f = s_curve(f);

    r = lerp( r.xyyy, r.zwww, f[1] );

    return lerp( r.x, r.y, f[0] );

}

 

// 1D version

float noise(float v, const uniform float4 pg[B2])

{

    v = v + 10000.0f;

   

    float i = frac(v * BR) * B;   // index between 0 and B-1

    float f = frac(v);            // fractional position

   

    // compute dot products between gradients and vectors

    float2 r;

    r[0] = pg[i].x * f;

    r[1] = pg[i + 1].x * (f - 1.0f);

   

    // interpolate

    f = s_curve(f);

    return lerp( r[0], r[1], f);

}

 

// Ridged multifractal terrain

// See "Texturing & Modeling, A Procedural Approach", Chapter 12

 

// ridge function

float2 ridge(float2 h, float offset)

{

    h = abs(h);

    h = offset - h;

    h = h * h;

    return h;

}

 

struct inputs

{

    float4 Position : POSITION;

    /*float4 Normal   : NORMAL;

    float4 Texture0 : TEXCOORD0;*/

};

 

struct outputs

{

    float4 Position : POSITION;

    float4 Color    : COLOR0;

    float3 L_in     : COLOR1;

    float3 TexCoord : TEXCOORD0;

    float3 F_ext    : TEXCOORD1;

};

 

outputs main(inputs In,

             uniform float4x4 modelViewProjection,

             uniform float4x4 noiseMatrix,

             uniform float4 freq_amp,              // noise frequency and amplitude

             uniform float4 tex_offset,

             uniform float4 tex_scale,

             uniform float4 terrain_param,

             uniform float4 pg[B2],                // permutation/gradient table

            

             uniform float4x4 Modelview,      // Modelview Transformation

             uniform float3   SunDir,         // Sun Direction assumed normalized

             uniform float3   K1,

             uniform float3   K2,

             uniform float3   K3,

             uniform float3   K4,

             uniform float2   K5

             )

{

  outputs Out;

 

  float4 noisePos = mul(noiseMatrix, In.Position);

 

  // get noise values

  float2 octave;

  octave[0] = noise(noisePos.xz * freq_amp[0], pg);

  octave[1] = noise(noisePos.xz * freq_amp[1], pg);

 

  // apply ridge function to both octaves

  octave = ridge(octave, terrain_param[0]);

 

  // sum octaves (scale second octave by first)

  float height = (octave[0]*freq_amp[2] + octave[1]*octave[0]*freq_amp[3]);

 

  Out.TexCoord.xy = noisePos.xz * tex_scale.xy;

 

  float y = (height - tex_scale.w) * tex_scale.z;

  Out.TexCoord.z = y;

 

  Out.Color.xyz = y;

 

  float4 P = In.Position;

  P.y = height;

  Out.Position = mul(modelViewProjection, P);

 

  /* Aerial Perspective */

      // Compute vertex distance (z-depth in eye space)

      float3 WorldPos = mul(Modelview, P).xyz;

      float s = length(WorldPos);

 

      // Extinction factor

    Out.F_ext = exp(K1 * -s);

 

      // Cosine of incoming sunlight angle 0 (theta)

      float3 VertexDir = normalize(WorldPos);

      float cos0 = dot(SunDir, VertexDir);

 

      // Rayleigh scattering from angle 0 (theta)

      float3 BetaR0 = 1 + cos0*cos0;

 

      // Mie scattering from angle 0

      float tmp = K5.x + K5.y * cos0;    

      float3 BetaM0 = pow(tmp,-1.5);

 

      // In-scattered color from sun

      Out.L_in = K2 * (K3*BetaR0 + K4*BetaM0);

 

  return Out;

}

 

 

struct inputs

{

    float4 Position : POSITION;   // Transformed Vertex Position

    float4 L0       : COLOR0;     // Object Color

    float3 L_in     : COLOR1;     // In-scattered Color

    float3 TexCoord : TEXCOORD0;  // Texture Coordinate

    float3 F_ext    : TEXCOORD1;  // Extinction Factor

};

 

fragout pAerialPersp (inputs IN, const uniform sampler3D BaseTexture)

{

      fragout OUT;

 

      float4 texColor = tex3D(BaseTexture, IN.TexCoord.xyz);

      OUT.col.xyz = lerp(IN.L_in, texColor.xyz, IN.F_ext);

 

      return OUT;

}