Light scattering FX: fixed bug with extinction not taking into account intersection with the Earth

1 files changed, 37 insertions, 12 deletions

diff --git a/Shaders/PostProcess/EpipolarLightScattering/private/Extinction.fxh b/Shaders/PostProcess/EpipolarLightScattering/private/Extinction.fxh
index ae2e9c6..6b09fae 100644
--- a/Shaders/PostProcess/EpipolarLightScattering/private/Extinction.fxh
+++ b/Shaders/PostProcess/EpipolarLightScattering/private/Extinction.fxh
@@ -7,18 +7,24 @@
 // (mu=cos(view zenith angle)), intersections with ground ignored
 float2 GetDensityIntegralAnalytic(float r, float mu, float d, float EarthRadius, float4 ParticleScaleHeight) 
 {
-    float2 f2A = sqrt( (float2(0.5, 0.5) * ParticleScaleHeight.zw) * r );
-    float4 f4A01 = f2A.xxyy * float2(mu, mu + d / r).xyxy;
-    float4 f4A01s = sign(f4A01);
-    float4 f4A01sq = f4A01*f4A01;
+    float2 f2A     = sqrt(ParticleScaleHeight.zw * 0.5 * r);
+    float4 f4A01   = f2A.xxyy * float2(mu, mu + d / r).xyxy;
+    float4 f4A01s  = sign(f4A01);
+    float4 f4A01sq = f4A01 * f4A01;
     
     float2 f2X;
     f2X.x = f4A01s.y > f4A01s.x ? exp(f4A01sq.x) : 0.0;
     f2X.y = f4A01s.w > f4A01s.z ? exp(f4A01sq.z) : 0.0;
     
-    float4 f4Y = f4A01s / (2.3193*abs(f4A01) + sqrt(1.52*f4A01sq + 4.0)) * float3(1.0, exp(-float2(d,d)*ParticleScaleHeight.zw*(d/(2.0*r)+mu))).xyxz;
+    float4 f4Y =
+        f4A01s 
+        / (2.3193 * abs(f4A01) + sqrt(1.52 * f4A01sq + 4.0)) 
+        * float3(1.0, exp(-ParticleScaleHeight.zw * d * (d / (2.0 * r) + mu)) ).xyxz;
 
-    return sqrt((6.2831*ParticleScaleHeight.xy)*r) * exp((EarthRadius-r) * ParticleScaleHeight.zw) * (f2X + float2( dot(f4Y.xy, float2(1.0, -1.0)), dot(f4Y.zw, float2(1.0, -1.0)) ));
+    return 
+        sqrt((6.2831 * ParticleScaleHeight.xy) * r) 
+        * exp((EarthRadius - r) * ParticleScaleHeight.zw) 
+        * (f2X + float2(f4Y.x - f4Y.y, f4Y.z - f4Y.w));
 }
 
 
@@ -57,14 +63,33 @@ float3 GetExtinction(float3 f3StartPos,
     float fRayLength = length(f3EyeDir);
     f3EyeDir /= fRayLength;
 
-    float2 f2RayAtmTopIsecs = float2(0.0, 0.0); 
-    // Compute intersections of the view ray with the atmosphere
-    GetRaySphereIntersection(f3StartPos, f3EyeDir, f3EarthCentre, fAtmTopRadius, f2RayAtmTopIsecs);
+    // Compute intersections of the view ray with the atmosphere and the Earth
+    float4 f4Isecs;
+    GetRaySphereIntersection2( f3StartPos, f3EyeDir, f3EarthCentre,
+                               float2(fEarthRadius, fAtmTopRadius), 
+                               f4Isecs);
+    float2 f2RayEarthIsecs  = f4Isecs.xy;
+    float2 f2RayAtmTopIsecs = f4Isecs.zw;
+
     // If the ray misses the atmosphere, there is no extinction
-    if( f2RayAtmTopIsecs.y < 0.0 )return float3(1.0, 1.0, 1.0);
+    if (f2RayAtmTopIsecs.y < 0.0)
+    {
+        return float3(1.0, 1.0, 1.0);
+    }
+
+    // Limit the ray length by the distance to the top of the atmosphere
+    fRayLength = min(fRayLength, f2RayAtmTopIsecs.y);
+
+    if (f2RayEarthIsecs.x > 0.0)
+    {
+        // The ray hits the Earth surface
+        fRayLength = min(fRayLength, f2RayEarthIsecs.x);
+    }
+
+    // Update the end point
+    f3EndPos = f3StartPos + f3EyeDir * fRayLength;
 
-    // Do not let the start and end point be outside the atmosphere
-    f3EndPos = f3StartPos + f3EyeDir * min(f2RayAtmTopIsecs.y, fRayLength);
+    // If the ray starts outside of the atmosphere, move the starting point to the intersection
     f3StartPos += f3EyeDir * max(f2RayAtmTopIsecs.x, 0.0);
 
     return GetExtinctionUnverified(f3StartPos, f3EndPos, f3EyeDir, f3EarthCentre, fEarthRadius, f4ParticleScaleHeight);