atmosphere and volumetric cloud updates (#264)

WickedEngine/shaders/skyAtmosphere.hlsli

@@ -357,6 +357,15 @@ float3 GetTransmittance(AtmosphereParameters atmosphere, float pHeight, float su
 
 float3 GetAtmosphereTransmittance(float3 worldPosition, float3 worldDirection, AtmosphereParameters atmosphere, Texture2D<float4> transmittanceLutTexture)
 {
+	// If the worldDirection is occluded from this virtual planet, then return.
+	// We do this due to the low resolution LUT, where the stored zenith to horizon never reaches black, to prevent linear interpolation artefacts.
+	// At the most shadowed point of the LUT, pure black with earth shadow is never reached.
+	float2 sol = RaySphereIntersect(worldPosition, worldDirection, float3(0.0f, 0.0f, 0.0f), atmosphere.bottomRadius);
+	if (sol.x > 0.0f || sol.y > 0.0f)
+	{
+		return 0.0f;
+	}
+	
 	float pHeight = length(worldPosition);
 	const float3 UpVector = worldPosition / pHeight;
 	float SunZenithCosAngle = dot(worldDirection, UpVector);

WickedEngine/shaders/volumetricCloud_renderCS.hlsl

@@ -105,6 +105,8 @@ static const float g_LODDistance = 25000.0; // After a certain distance, noises
 static const float g_LODMin = 0.0; // 
 static const float g_BigStepMarch = 3.0; // How long inital rays should be until they hit something. Lower values may ives a better image but may be slower.
 static const float g_TransmittanceThreshold = 0.005; // Default: 0.005. If the clouds transmittance has reached it's desired opacity, there's no need to keep raymarching for performance.
+static const float g_ShadowSampleCount = 5.0f;
+static const float g_GroundContributionSampleCount = 2.0f;
 
 
 float GetHeightFractionForPoint(AtmosphereParameters atmosphere, float3 pos)
@@ -280,19 +282,24 @@ void VolumetricShadow(inout ParticipatingMedia participatingMedia, in Atmosphere
 		extinctionAccumulation[ms] = 0.0f;
 	}
 	
-	const float shadowStepCount = 5.0;
-	const float invShadowStepCount = 1.0 / shadowStepCount;
-
-	float previousNormT = 0.0;
+	const float sampleCount = g_ShadowSampleCount;
+	const float sampleSegmentT = 0.5f;
+	
 	float lodOffset = 0.5;
-	for (float shadowT = invShadowStepCount; shadowT <= 1.0; shadowT += invShadowStepCount)
+	for (float s = 0.0f; s < sampleCount; s += 1.0)
 	{
-		float currentNormT = shadowT * shadowT;
-		float deltaNormT = currentNormT - previousNormT; // 5 samples: 0.04, 0.12, 0.2, 0.28, 0.36
-		float extinctionFactor = deltaNormT;
-		float shadowSampleDistance = g_ShadowStepLength * (previousNormT + deltaNormT * 0.5); // 5 samples: 0.02, 0.1, 0.26, 0.5, 0.82
+		// More expensive but artefact free
+		float t0 = (s) / sampleCount;
+		float t1 = (s + 1.0) / sampleCount;
+		// Non linear distribution of sample within the range.
+		t0 = t0 * t0;
+		t1 = t1 * t1;
 
-		float3 samplePoint = worldPosition + sunDirection * shadowSampleDistance; // Step futher towards the light
+		float delta = t1 - t0; // 5 samples: 0.04, 0.12, 0.2, 0.28, 0.36
+		float t = t0 + delta * sampleSegmentT; // 5 samples: 0.02, 0.1, 0.26, 0.5, 0.82
+		
+		float shadowSampleT = g_ShadowStepLength * t;
+		float3 samplePoint = worldPosition + sunDirection * shadowSampleT; // Step futher towards the light
 
 		float heightFraction = GetHeightFractionForPoint(atmosphere, samplePoint);
 		if (heightFraction < 0.0 || heightFraction > 1.0)
@@ -314,10 +321,9 @@ void VolumetricShadow(inout ParticipatingMedia participatingMedia, in Atmosphere
 		[unroll]
 		for (ms = 0; ms < MS_COUNT; ms++)
 		{
-			extinctionAccumulation[ms] += shadowParticipatingMedia.extinctionCoefficients[ms] * extinctionFactor;
+			extinctionAccumulation[ms] += shadowParticipatingMedia.extinctionCoefficients[ms] * delta;
 		}
 
-		previousNormT = currentNormT;
 		lodOffset += 0.5;
 	}
 
@@ -328,6 +334,74 @@ void VolumetricShadow(inout ParticipatingMedia participatingMedia, in Atmosphere
 	}
 }
 
+void VolumetricGroundContribution(inout float3 environmentLuminance, in AtmosphereParameters atmosphere, float3 worldPosition, float3 sunDirection, float3 sunIlluminance, float3 atmosphereTransmittanceToLight, float3 windOffset, float3 windDirection, float2 coverageWindOffset, float lod)
+{
+	float planetRadius = atmosphere.bottomRadius * SKY_UNIT_TO_M;
+	float3 planetCenterWorld = atmosphere.planetCenter * SKY_UNIT_TO_M;
+
+	float cloudBottomRadius = planetRadius + g_CloudStartHeight;
+
+	float cloudSampleAltitudde = length(worldPosition - planetCenterWorld); // Distance from planet center to tracing sample
+	float cloudSampleHeightToBottom = cloudSampleAltitudde - cloudBottomRadius; // Distance from altitude to bottom of clouds
+	
+	float3 opticalDepth = 0.0;
+	
+	const float contributionStepLength = min(4000.0, cloudSampleHeightToBottom);
+	const float3 groundScatterDirection = float3(0.0, -1.0, 0.0);
+	
+	const float sampleCount = g_GroundContributionSampleCount;
+	const float sampleSegmentT = 0.5f;
+	
+	// Ground Contribution tracing loop, same idea as volumetric shadow
+	float lodOffset = 0.5;
+	for (float s = 0.0f; s < sampleCount; s += 1.0)
+	{
+		// More expensive but artefact free
+		float t0 = (s) / sampleCount;
+		float t1 = (s + 1.0) / sampleCount;
+		// Non linear distribution of sample within the range.
+		t0 = t0 * t0;
+		t1 = t1 * t1;
+
+		float delta = t1 - t0; // 5 samples: 0.04, 0.12, 0.2, 0.28, 0.36		
+		float t = t0 + (t1 - t0) * sampleSegmentT; // 5 samples: 0.02, 0.1, 0.26, 0.5, 0.82
+
+		float contributionSampleT = contributionStepLength * t;
+		float3 samplePoint = worldPosition + groundScatterDirection * contributionSampleT; // Step futher towards the scatter direction
+
+		float heightFraction = GetHeightFractionForPoint(atmosphere, samplePoint);
+		/*if (heightFraction < 0.0 || heightFraction > 1.0) // No impact
+		{
+			break;
+		}*/
+		
+		float3 weatherData = SampleWeather(samplePoint, heightFraction, coverageWindOffset);
+		if (weatherData.r < 0.4)
+		{
+			continue;
+		}
+
+		float contributionCloudDensity = SampleCloudDensity(samplePoint, heightFraction, weatherData, windOffset, windDirection, lod + lodOffset, true);
+
+		float3 contributionExtinction = g_ExtinctionCoefficient * contributionCloudDensity;
+
+		opticalDepth += contributionExtinction * contributionStepLength * delta;
+		
+		lodOffset += 0.5;
+	}
+	
+	const float3 planetSurfaceNormal = float3(0.0, 1.0, 0.0); // Ambient contribution from the clouds is only done on a plane above the planet
+	const float3 groundBrdfNdotL = saturate(dot(sunDirection, planetSurfaceNormal)) * (atmosphere.groundAlbedo / PI); // Lambert BRDF diffuse shading
+
+	const float uniformPhase = UniformPhase();
+	const float groundHemisphereLuminanceIsotropic = (2.0f * PI) * uniformPhase; // Assumes the ground is uniform luminance to the cloud and solid angle is bottom hemisphere 2PI
+	const float3 groundToCloudTransfertIsoScatter = groundBrdfNdotL * groundHemisphereLuminanceIsotropic;
+	
+	const float3 scatteredLuminance = atmosphereTransmittanceToLight * sunIlluminance * groundToCloudTransfertIsoScatter;
+
+	environmentLuminance += scatteredLuminance * exp(-opticalDepth);
+}
+
 struct ParticipatingMediaPhase
 {
 	float phase[MS_COUNT];
@@ -398,8 +472,24 @@ void VolumetricCloudLighting(AtmosphereParameters atmosphere, float3 startPositi
 	ParticipatingMedia participatingMedia = SampleParticipatingMedia(albedo, extinction, g_MultiScatteringScattering, g_MultiScatteringExtinction, atmosphereTransmittanceToLight);
 	
 
-	// Calcualte volumetric shadow
-	VolumetricShadow(participatingMedia, atmosphere, worldPosition, sunDirection, windOffset, windDirection, coverageWindOffset, lod);
+	// Sample environment lighting
+	float3 environmentLuminance = SampleAmbientLight(heightFraction);
+
+
+	// Only render if there is any sign of scattering (albedo * extinction)
+	if (any(participatingMedia.scatteringCoefficients[0] > 0.0))
+	{
+		// Calcualte volumetric shadow
+		VolumetricShadow(participatingMedia, atmosphere, worldPosition, sunDirection, windOffset, windDirection, coverageWindOffset, lod);
+
+
+		// Calculate bounced light from ground onto clouds
+		const float maxTransmittanceToView = max(max(transmittanceToView.x, transmittanceToView.y), transmittanceToView.z);
+		if (maxTransmittanceToView > 0.01f)
+		{
+			VolumetricGroundContribution(environmentLuminance, atmosphere, worldPosition, sunDirection, sunIlluminance, atmosphereTransmittanceToLight, windOffset, windDirection, coverageWindOffset, lod);
+		}
+	}
 
 
 	// Sample dual lob phase with multiple scattering
@@ -407,11 +497,6 @@ void VolumetricCloudLighting(AtmosphereParameters atmosphere, float3 startPositi
 	ParticipatingMediaPhase participatingMediaPhase = SampleParticipatingMediaPhase(phaseFunction, g_MultiScatteringEccentricity);
 
 
-	// Sample environment lighting
-	// Todo: Ground contribution?
-	float3 environmentLuminance = SampleAmbientLight(heightFraction);
-	
-
 	// Update depth sampling
 	float depthWeight = min(transmittanceToView.r, min(transmittanceToView.g, transmittanceToView.b));
 	depthWeightedSum += depthWeight * length(worldPosition - startPosition);

WickedEngine/shaders/volumetricCloud_reprojectCS.hlsl

@@ -25,20 +25,21 @@ inline float HdrWeight4(float3 color, float exposure)
 	return rcp(Luma4(color) * exposure + 4.0f);
 }
 
-float4 clip_aabb(float3 aabb_min, float3 aabb_max, float4 p, float4 q)
+// Different aabb clipping method from eg. SSR temporal, suitable for clouds in this case
+float4 clip_aabb(float4 aabb_min, float4 aabb_max, float4 prev_sample)
 {
-	float3 p_clip = 0.5 * (aabb_max + aabb_min);
-	float3 e_clip = 0.5 * (aabb_max - aabb_min) + 0.00000001f;
+	float4 p_clip = 0.5 * (aabb_max + aabb_min);
+	float4 e_clip = 0.5 * (aabb_max - aabb_min) + 0.00000001f;
 
-	float4 v_clip = q - float4(p_clip, p.w);
-	float3 v_unit = v_clip.xyz / e_clip;
-	float3 a_unit = abs(v_unit);
-	float ma_unit = max(a_unit.x, max(a_unit.y, a_unit.z));
+	float4 v_clip = prev_sample - p_clip;
+	float4 v_unit = v_clip / e_clip;
+	float4 a_unit = abs(v_unit);
+	float ma_unit = max(max(a_unit.x, max(a_unit.y, a_unit.z)), a_unit.w);
 
 	if (ma_unit > 1.0)
-		return float4(p_clip, p.w) + v_clip / ma_unit;
+		return p_clip + v_clip / ma_unit;
 	else
-		return q; // point inside aabb
+		return prev_sample; // point inside aabb
 }
 
 inline void ResolverAABB(Texture2D<float4> currentColor, SamplerState currentSampler, float sharpness, float exposureScale, float AABBScale, float2 uv, float2 texelSize, inout float4 currentMin, inout float4 currentMax, inout float4 currentAverage, inout float4 currentOutput)
@@ -73,6 +74,8 @@ inline void ResolverAABB(Texture2D<float4> currentColor, SamplerState currentSam
                        sampleColors[5] * sampleWeights[5] + sampleColors[6] * sampleWeights[6] + sampleColors[7] * sampleWeights[7] + sampleColors[8] * sampleWeights[8]) / totalWeight;
 #endif
 
+
+#if 0 // Standard clipping
 	
     // Variance Clipping (AABB)
     
@@ -87,7 +90,39 @@ inline void ResolverAABB(Texture2D<float4> currentColor, SamplerState currentSam
 
 	float4 mean = m1 / 9.0;
 	float4 stddev = sqrt((m2 / 9.0) - sqr(mean));
+
+#else // Depth check
+
+	float originalLinearDepth = getLinearDepth(texture_depth.SampleLevel(sampler_point_clamp, uv, 0).r);
+	float validSampleCount = 1.0;
+	
+	float4 m1 = 0.0;
+	float4 m2 = 0.0;
+    [unroll]
+	for (uint x = 0; x < 9; x++)
+	{
+		if (x == 4)
+		{
+			m1 += sampleColors[x];
+			m2 += sampleColors[x] * sampleColors[x];
+		}
+		else
+		{
+			float depth = getLinearDepth(texture_depth.SampleLevel(sampler_point_clamp, uv + (SampleOffset[x] / texelSize), 0).r);
+			if (abs(originalLinearDepth - depth) < 1.5)
+			{
+				m1 += sampleColors[x];
+				m2 += sampleColors[x] * sampleColors[x];
+				validSampleCount += 1.0;
+			}
+		}
+	}
+
+	float4 mean = m1 / validSampleCount;
+	float4 stddev = sqrt((m2 / validSampleCount) - sqr(mean));
         
+#endif
+
 	currentMin = mean - AABBScale * stddev;
 	currentMax = mean + AABBScale * stddev;
 
@@ -169,8 +204,9 @@ void main(uint3 DTid : SV_DispatchThreadID)
 	float4 currentMin, currentMax, currentAverage;
 	ResolverAABB(cloud_current, sampler_point_clamp, 0, temporalExposure, temporalScale, uv, xPPResolution, currentMin, currentMax, currentAverage, current);
 
-	previous = clip_aabb(currentMin.xyz, currentMax.xyz, clamp(currentAverage, currentMin, currentMax), previous);
-    
+	//previous = clip_aabb(currentMin.xyz, currentMax.xyz, clamp(currentAverage, currentMin, currentMax), previous);
+	previous = clip_aabb(currentMin, currentMax, previous);
+
 	float4 result = lerp(previous, current, temporalResponse);
     
 	result = is_saturated(prevUV) ? result : current;