HDR: Better handling of the Sun

Add transmittance to the Sun illuminance value to get correct colors at every time of day.
Also improved the Sun disk rendering.

Compositor/HDR/hdr.xml

@@ -672,6 +672,10 @@
       <unit>11</unit>
       <buffer>aerial-perspective</buffer>
     </binding>
+    <binding>
+      <unit>12</unit>
+      <buffer>transmittance</buffer>
+    </binding>
     <attachment>
       <component>color0</component>
       <buffer>hdr-result</buffer>

Effects/HDR/lighting.eff

@@ -72,6 +72,11 @@
         <type>sampler-2d</type>
         <value type="int">11</value>
       </uniform>
+      <uniform>
+        <name>transmittance_lut</name>
+        <type>sampler-2d</type>
+        <value type="int">12</value>
+      </uniform>
     </pass>
   </technique>
 </PropertyList>

Effects/skydome.eff

@@ -314,6 +314,11 @@
         <vertex-shader>Shaders/HDR/skydome.vert</vertex-shader>
         <fragment-shader>Shaders/HDR/skydome.frag</fragment-shader>
       </program>
+      <uniform>
+        <name>sun_disk</name>
+        <type>bool</type>
+        <value type="bool">true</value>
+      </uniform>
       <uniform>
         <name>sky_view_lut</name>
         <type>sampler-2d</type>
@@ -340,6 +345,11 @@
         <vertex-shader>Shaders/HDR/skydome.vert</vertex-shader>
         <fragment-shader>Shaders/HDR/skydome.frag</fragment-shader>
       </program>
+      <uniform>
+        <name>sun_disk</name>
+        <type>bool</type>
+        <value type="bool">false</value>
+      </uniform>
       <uniform>
         <name>sky_view_lut</name>
         <type>sampler-2d</type>

Shaders/HDR/lighting.frag

@@ -13,11 +13,14 @@ uniform samplerCube prefiltered_envmap;
 uniform sampler2DShadow shadow_tex;
 uniform sampler2D dfg_lut;
 uniform sampler2D aerial_perspective_lut;
+uniform sampler2D transmittance_lut;
 
 uniform mat4 fg_ViewMatrix;
 uniform mat4 fg_ViewMatrixInverse;
 uniform vec3 fg_SunDirection;
+uniform vec3 fg_SunDirectionWorld;
 uniform vec3 fg_CameraPositionCart;
+uniform vec3 fg_CameraPositionGeod;
 
 uniform mat4 fg_LightMatrix_csm0;
 uniform mat4 fg_LightMatrix_csm1;
@@ -48,7 +51,8 @@ const float AERIAL_MAX_DEPTH = 32000.0;
 
 const float MAX_PREFILTERED_LOD = 4.0;
 
-const vec3 SUN_INTENSITY = vec3(20.0);
+const float ATMOSPHERE_RADIUS = 6471e3;
+const vec3 EXTRATERRESTRIAL_SOLAR_ILLUMINANCE = vec3(128.0);
 
 vec3 decodeNormal(vec2 enc);
 vec3 positionFromDepth(vec2 pos, float depth);
@@ -359,6 +363,22 @@ vec4 sampleAerialPerspective(float depth)
     return color;
 }
 
+vec3 getSunIlluminance()
+{
+    float cameraHeight = length(fg_CameraPositionCart);
+    vec3 up = fg_CameraPositionCart / cameraHeight;
+    float cosTheta = dot(fg_SunDirectionWorld, up);
+
+    float earthRadius = cameraHeight - max(fg_CameraPositionGeod.z, 0.0);
+    float normalizedHeight = (cameraHeight - earthRadius)
+        / (ATMOSPHERE_RADIUS - earthRadius);
+
+    vec2 coords = vec2(cosTheta * 0.5 + 0.5, clamp(normalizedHeight, 0.0, 1.0));
+    vec3 transmittance = texture(transmittance_lut, coords).rgb;
+
+    return EXTRATERRESTRIAL_SOLAR_ILLUMINANCE * transmittance;
+}
+
 //------------------------------------------------------------------------------
 
 void main()
@@ -399,7 +419,7 @@ void main()
                      NdotL, NdotV, NdotH, VdotH,
                      f0);
 
-    vec3 sunIlluminance = SUN_INTENSITY * NdotL;
+    vec3 sunIlluminance = getSunIlluminance() * NdotL;
 
     vec3 f = F_SchlickRoughness(NdotV, f0, roughness);
     vec3 indirectSpecular = IBL_Specular(n, v, NdotV, roughness, f);
@@ -411,7 +431,8 @@ void main()
     vec3 color = ambient + brdf * sunIlluminance * shadowFactor;
 
     vec4 aerialPerspective = sampleAerialPerspective(length(pos));
-    color = color * aerialPerspective.a + aerialPerspective.rgb * SUN_INTENSITY;
+    color = color * aerialPerspective.a + aerialPerspective.rgb
+        * EXTRATERRESTRIAL_SOLAR_ILLUMINANCE;
 
     fragHdrColor = color;
 }

Shaders/HDR/skydome.frag

@@ -5,14 +5,17 @@ out vec4 fragColor;
 in vec3 vRayDir;
 in vec3 vRayDirView;
 
+uniform bool sun_disk;
 uniform sampler2D sky_view_lut;
 uniform sampler2D transmittance_lut;
 uniform vec3 fg_SunDirection;
 
-const float PI = 3.141592653;
-const vec3 SUN_INTENSITY = vec3(20.0);
+uniform mat4 osg_ModelViewMatrix;
 
-const float sun_solid_angle = 0.53*PI/180.0; // ~half a degree
+const float PI = 3.141592653;
+const vec3 EXTRATERRESTRIAL_SOLAR_ILLUMINANCE = vec3(128.0);
+
+const float sun_solid_angle = 0.545*PI/180.0; // ~half a degree
 const float sun_cos_solid_angle = cos(sun_solid_angle);
 
 void main()
@@ -24,13 +27,37 @@ void main()
     float elev = sqrt(abs(l) / (PI * 0.5)) * sign(l) * 0.5 + 0.5;
 
     vec3 color = texture(sky_view_lut, vec2(azimuth, elev)).rgb;
-    color *= SUN_INTENSITY;
+    color *= EXTRATERRESTRIAL_SOLAR_ILLUMINANCE;
 
-    // Render the Sun disk
-    // XXX: Apply transmittance
-    float cosTheta = dot(normalize(vRayDirView), fg_SunDirection);
-    if (cosTheta >= sun_cos_solid_angle)
-        color += SUN_INTENSITY;
+    if (sun_disk) {
+        // Render the Sun disk
+        vec3 rayDirView = normalize(vRayDirView);
+        float cosTheta = dot(rayDirView, fg_SunDirection);
+
+        if (cosTheta >= sun_cos_solid_angle) {
+            vec4 groundPoint = osg_ModelViewMatrix * vec4(0.0, 0.0, 0.0, 1.0);
+            float altitude = length(groundPoint);
+            float scaledAltitude = altitude / 100000.0;
+
+            vec3 up = normalize(vec4(0.0, 0.0, 0.0, 1.0) - groundPoint).xyz;
+            float sunZenithCosTheta = dot(rayDirView, up);
+
+            vec2 coords = vec2(sunZenithCosTheta * 0.5 + 0.5,
+                               clamp(scaledAltitude, 0.0, 1.0));
+            vec3 transmittance = texture(transmittance_lut, coords).rgb;
+
+            // Limb darkening
+            // http://www.physics.hmc.edu/faculty/esin/a101/limbdarkening.pdf
+            vec3 u = vec3(1.0);
+            vec3 a = vec3(0.397, 0.503, 0.652);
+            float centerToEdge = 1.0 - (cosTheta - sun_cos_solid_angle)
+                / (1.0 - sun_cos_solid_angle);
+            float mu = sqrt(max(1.0 - centerToEdge * centerToEdge, 0.0));
+            vec3 factor = vec3(1.0) - u * (vec3(1.0) - pow(vec3(mu), a));
+
+            color += EXTRATERRESTRIAL_SOLAR_ILLUMINANCE * transmittance * factor;
+        }
+    }
 
     fragColor = vec4(color, 1.0);
 }