Atmospheric Scattering

46407cc1bdfbd2db4f6e8876d74f990a
0
Kenneth_Gorking 101 Feb 09, 2005 at 12:23

I have implemented atmospheric scattering in my engine, the one made by Nathaniel Hoffman and Arcot J. Preetham. Click and be amazed:

http://www.ati.com/developer/dx9/ATI-LightScattering.pdf

http://www.ati.com/developer/SIGGRAPH03/Pr…CourseNotes.pdf

http://www.imm.dtu.dk/pubdb/views/edoc_dow…pdf/imm2554.pdf

Not that this is the “old” version. He did another for the book “Graphics Programming Methods”, which is where I got my inspiration. The demo he supplied on the accompanying cd kicks ass.

My problem however is that after I have implemented it into my engine, all I get is a very bright spot on an otherwise pitch black sky. That is clearly not supposed to happen. I am using all of their values to the letter, and it still doesn’t work…

I did convert the source he provided from assembly to Cg, and I fear I made an error somewhere. I compared the compiled output from my Cg version to his assembly version, and they are NOTHING alike. I don’t know what the hell the Cg compiler is doing, but the output sure is mangled. Here is my Cg code:

//
// Shader used for calculating atmospheric scattering
//

// constants
static const float log2_e = 1.0f/log(2.0f);

// Structure used for easier setting of parameters
struct ScatteringParams
{
    float4 eye;
    float4 sunDirection;
    float4 beta1_plus_Beta2;
    float4 terrainReflectance;
    float4 hg;
    float4 beta1;
    float4 beta2;
    float4 betaDash1;
    float4 betaDash2;
    float4 sunColor;
    float4 oneOver_Beta1_plus_Beta2;
    float4 termMultiplier;
};

//
// Calculate atmospheric scattering
void CalcAtmosphericScattering_VP
(
    float4 position,
    float4x4 world,
    float4x4 worldView,
    ScatteringParams params,
    out float4 OutExtinction,
    out float4 OutInscattering
)
{
    if(g_bUseAtmosphericScattering == true)
    {
 // Transform position to world space
 float4 V = mul(position, world);
 V = normalize(params.eye - V);

 // Terms used in the scattering equation.
 // Scattering = [cos(theta), 1+cos^2(theta), s] 
 float4 Scattering;
 Scattering.x = dot(V, params.sunDirection); // cos(theta) = V.L
 Scattering.y = 1+(Scattering.x*Scattering.x);  // 1+cos^2(theta) = Phase1(theta)
 Scattering.z = mul(position, worldView).z; // s

 // e^(-(beta_1 + beta_2) * s)
 float4 Extinction = -params.beta1_plus_Beta2 * Scattering.z * log2_e;
 Extinction.x = exp( Extinction.x );
 Extinction.y = exp( Extinction.y );
 Extinction.z = exp( Extinction.z );

 // Apply Reflectance to E to get total net effective E
 float4 TotalExtinction = Extinction * params.terrainReflectance;

 // Phase2(theta) = (1-g^2)/(1+g^2-2g*cos(theta))^(3/2)
 // theta is 180 - actual theta (this corrects for sign)
 float Phase2;
 Phase2 = 1/((Scattering.x * params.hg.z) + params.hg.y);
 Phase2 = Phase2 * Phase2;
 Phase2 = Phase2 * Phase2;
 Scattering.w = Phase2 * params.hg.x;

 // Inscattering (I) = (Beta'_1 * Phase_1(theta) + Beta'_2 * Phase_2(theta)) * [1-exp(-Beta_1*s).exp(-Beta_2*s)] / (Beta_1 + Beta_2)
 float4 Inscattering;
 Inscattering = (params.betaDash1*Scattering.y);
 Inscattering += (params.betaDash2*Scattering.w);
 Inscattering *= (1-Extinction);
 Inscattering *= params.oneOver_Beta1_plus_Beta2;

 // Apply Inscattering contribution factors.
 Inscattering *= params.termMultiplier.y;

 // Scale with Sun color & intesity.
 Inscattering   *= params.sunColor;
 Inscattering   *= params.sunColor.w;
 TotalExtinction    *= params.sunColor;
 TotalExtinction    *= params.sunColor.w;

 OutExtinction = TotalExtinction;
 OutInscattering = Inscattering;
    }
    else
    {
 // No scattering
 OutExtinction = float4(1,1,1,1);
 OutInscattering = float4(0,0,0,0);
    }
}

float4 CalcAtmosphericScattering_FP( float4 SurfColor, float4 Extinction, float4 Inscattering)
{
    return (SurfColor * Extinction) + Inscattering;
}

And the sky code that uses it:

// Set this to true if scattering is wanted
static const bool g_bUseAtmosphericScattering = true;

#include "AtmosphericScattering.cg"


// Inputs to sky shader
struct SkyInput
{
    float3 Position : POSITION;
    float2 TexCoord0    : TEXCOORD0;
};

// Outputs from sky shader
struct SkyOutput
{
    float4 HPosition    : POSITION;
    float4 Extinction   : COLOR0;
    float4 Inscattering : COLOR1;
};

SkyOutput vertmain( const SkyInput IN,
    uniform float4x4 mvp,
    uniform float4x4 world,
    uniform float4x4 worldView,
    uniform ScatteringParams scatteringParams
  )
{
    SkyOutput OUT;

    float4 pos = float4(IN.Position,1);

    // Transform vertex position into homogenous clip-space.
    OUT.HPosition = mul( pos, mvp );

    CalcAtmosphericScattering_VP
    (
 pos,
 world,
 worldView,
 scatteringParams,
 OUT.Extinction,
 OUT.Inscattering
    );

    return OUT;
}

// Sky is just the inscattering term
float4 fragmain(const SkyOutput IN) : COLOR
{
    return IN.Inscattering;
}

And finally the output from cgc:

sky.cg
./AtmosphericScattering.cg
151 lines, 0 errors.
vs.1.1
// DX8 Vertex shader generated by NVIDIA Cg compiler
// cgc version 1.3.0001, build date Jan 7 2005 14:01:35
// command line args: -profile vs_1_1 -entry vertmain
// source file: sky.cg
// source file: ./AtmosphericScattering.cg
// nv30vp backend compiling 'vertmain' program
def c24, 1, 2, 0.69314718, 1.442695
def c25, 2.71828, 0, 0, 0
//vendor NVIDIA Corporation
//version 1.0.02
//profile vs_1_1
//program vertmain
//semantic vertmain.mvp
//semantic vertmain.world
//semantic vertmain.worldView
//semantic vertmain.scatteringParams
//var float4x4 mvp : : c[0], 4 : 1 : 1
//var float4x4 world : : c[4], 4 : 2 : 1
//var float4x4 worldView : : c[8], 4 : 3 : 1
//var float4 scatteringParams.eye : : c[12] : 4 : 1
//var float4 scatteringParams.sunDirection : : c[13] : 4 : 1
//var float4 scatteringParams.beta1_plus_Beta2 : : c[14] : 4 : 1
//var float4 scatteringParams.terrainReflectance : : c[15] : 4 : 1
//var float4 scatteringParams.hg : : c[16] : 4 : 1
//var float4 scatteringParams.beta1 : : c[17] : 4 : 1
//var float4 scatteringParams.beta2 : : c[18] : 4 : 1
//var float4 scatteringParams.betaDash1 : : c[19] : 4 : 1
//var float4 scatteringParams.betaDash2 : : c[20] : 4 : 1
//var float4 scatteringParams.sunColor : : c[21] : 4 : 1
//var float4 scatteringParams.oneOver_Beta1_plus_Beta2 : : c[22] : 4 : 1
//var float4 scatteringParams.termMultiplier : : c[23] : 4 : 1
//var float3 IN.Position : $vin.POSITION : ATTR0 : 0 : 1
//var float2 IN.TexCoord0 : $vin.TEXCOORD0 : ATTR7 : 0 : 0
//var float4 vertmain.HPosition : $vout.POSITION : HPOS : -1 : 1
//var float4 vertmain.Extinction : $vout.COLOR0 : COL0 : -1 : 1
//var float4 vertmain.Inscattering : $vout.COLOR1 : COL1 : -1 : 1
//const c[24] = 1 2 0.6931472 1.442695
//const c[25] = 2.718282 0 0 0
    mov r0.xyz, v0.xyzz
    mov r0.w, c24.x
    mul r1, r0.y, c1
    mad r1, r0.x, c0, r1
    mad r1, r0.z, c2, r1
    mad oPos, r0.w, c3, r1
    mul r1, r0.y, c5
    mad r1, r0.x, c4, r1
    mad r1, r0.z, c6, r1
    mad r1, r0.w, c7, r1
    add r2, c12, -r1
    dp4 r1.x, r2, r2
    rsq r1.x, r1.x
    mul r1, r1.x, r2
    dp4 r2.xyw, r1, c13
    mad r2.y, r2.x, r2.x, c24.x
    mul r1, r0.y, c9
    mad r1, r0.x, c8, r1
    mad r1, r0.z, c10, r1
    mad r1.z, r0.w, c11, r1
    mad r0.x, r2.x, c16.z, c16.y
    rcp r0.x, r0.x
    mul r0.x, r0.x, r0.x
    mul r0.x, r0.x, r0.x
    mul r0.xyz, r0.x, c16.x
    mov r2.w, r0.xyzx
    mul r0, c20, r2.w
    mad r2, c19, r2.y, r0
    mul r0, -c14, r1.z
    mul r1, r0, c24.w
    mov r0.xy, c25.x
    mov r0.zw, r1.x
    lit r0.z, r0
    mov r1.x, r0.z
    mov r0.xy, c25.x
    mov r0.zw, r1.y
    lit r0.z, r0
    mov r1.y, r0.z
    mov r0.xy, c25.x
    mov r0.zw, r1.z
    lit r1.z, r0
    add r0, c24.x, -r1
    mul r0, r2, r0
    mul r0, r0, c22
    mul r0, r0, c23.y
    mul r0, r0, c21
    mul oD1, r0, c21.w
    mul r0, r1, c15
    mul r0, r0, c21
    mul oD0, r0, c21.w
// 51 instructions
// 3 temp registers

Notice how the tree exp() have dissapered…

I hope that someone else have any experience with this, and can help me out.

5 Replies

Please log in or register to post a reply.

22b3033832c5c699c856814b0cf80cb1
0
bladder 101 Feb 10, 2005 at 08:51

You get the black atmosphere with bright white in the sun effect when the raleigh coefficient is low right? I dont know why the exp functions are not in the assembly output. They definetly should be there. I think your best bet would be to email nvidia and show them the hlsl file and the assembly output and ask them why the exp is not in the assembly output, cause it definetly should be.

Why dont you try setting the extinction values to some constants and see if that works out?

46407cc1bdfbd2db4f6e8876d74f990a
0
Kenneth_Gorking 101 Feb 10, 2005 at 09:23

Well, after playing around with the code last night, I found that the problem is twofold.

The first error was my own. :blush:
In the code that creates the sun, the intensity of the sun defaulted to 1, when it should be 100. Stupid.

This however was a “good” error, because it got me thinking about the massive size of the sun-disc (about 20 degrees!). It would seem that some value was too large, so I went over Cg code, and found an error in Hoffmans Henney-Greenstein phase function. Here is the formula for the function:

Phase2(theta) = (1-g\^2)/(1+g\^2-2g*cos(theta))\^(3/2)

Hoffman translated this into:

// params.hg = [1-g*g, 1+g*g, 2*g, 1]
float Phase2;
Phase2 = 1/((Scattering.x * params.hg.z) + params.hg.y);
Phase2 = Phase2 * Phase2;
Phase2 = Phase2 * Phase2;
Scattering.w = Phase2 * params.hg.x;

Which is clearly wrong, because it becomes:

Phase2(theta) = (1-g\^2)/(1+g\^2-2g*cos(theta))\^3

Changing it to the following, yields a much better result:

// Correct formula. Looks much better!
float Phase2;
Phase2 = params.hg.x / (pow(params.hg.y + params.hg.z * Scattering.x, (3.0f/2.0f)));
Scattering.w = Phase2;
22b3033832c5c699c856814b0cf80cb1
0
bladder 101 Feb 10, 2005 at 09:59

so is the assembly output wrong? or is there some there some other way that exp is being calculated (since its not in the output)

025594d3a6db748ec1f8be6842441d2c
0
robocrop 101 May 24, 2005 at 22:58

@Kenneth Gorking

This however was a “good” error, because it got me thinking about the massive size of the sun-disc (about 20 degrees!). It would seem that some value was too large, so I went over Cg code, and found an error in Hoffmans Henney-Greenstein phase function. Here is the formula for the function:

Phase2(theta) = (1-g\^2)/(1+g\^2-2g*cos(theta))\^(3/2)

Hoffman translated this into:

// params.hg = [1-g*g, 1+g*g, 2*g, 1]
float Phase2;
Phase2 = 1/((Scattering.x * params.hg.z) + params.hg.y);
Phase2 = Phase2 * Phase2;
Phase2 = Phase2 * Phase2;
Scattering.w = Phase2 * params.hg.x;

Which is clearly wrong, because it becomes:

Phase2(theta) = (1-g\^2)/(1+g\^2-2g*cos(theta))\^3

Changing it to the following, yields a much better result:

// Correct formula. Looks much better!
float Phase2;
Phase2 = params.hg.x / (pow(params.hg.y + params.hg.z * Scattering.x, (3.0f/2.0f)));
Scattering.w = Phase2;

[snapback]15958[/snapback]

You are reading Hoffman’s code incorrectly, perhaps because you transcribed it into Cg.

In the original shader code, here is Hoffman’s algorithm:

; Phase2(theta) = (1-g^2)/(1+g^2-2g*cos(theta))^(3/2)
; theta is 180 - actual theta (this corrects for sign)
; c[28] = [1-g^2, 1+g, 2g]
mad r4.x, c[28].z, r0.x, c[28].y; 


rsq r4.x, r4.x   
mul r4.y, r4.x, r4.x    
mul r4.x, r4.y, r4.x;
mul r0.w, r4.x, c[28].x      ; r0.w = Phase2(theta)

Your translation of the code somehow lost the rsq.

Your alternate method actually provides the same results as Hoffman’s code, so I would hope that it looks good :tongue:

However I find that this sample doesn’t really look as good in practice as it does at first glance (or on the fly-through). The sun tends to disperse greatly forming a huge disc, and attempts to make it more directional lose so much of the ambient lighting that it is virtually useless. Does anyone know of other scattering examples which might deal with these issues?

46407cc1bdfbd2db4f6e8876d74f990a
0
Kenneth_Gorking 101 May 28, 2005 at 12:04

Well that sucks! I guess i read it as a rcp…

Anyways, there are some more errors in the code I have discovered.

In the setup on the CPU for their shader, (1-g)\^2 is calculated as 1-g*g, which distorts the phase function.

Also, instead of the pow(a\^(3/2)) they use the sqrt(a\^3), which is fine, except that they calculate it wrong, and end up with sqrt(a\^4), which makes the sun disc too large.

The first phase function, 1+cos\^2(theta), is missing a calculation, (3/(16*3.141592)) needs to be multiplied with the result, or the result is bi-directional, meaning that sunsets are visible on the opposite side of the dome.