Files
FlaxEngine/Source/Shaders/AtmosphereFog.hlsl
2024-02-26 19:00:48 +01:00

550 lines
22 KiB
HLSL

// Copyright (c) 2012-2024 Wojciech Figat. All rights reserved.
#ifndef __ATMOSPHERE_FOG__
#define __ATMOSPHERE_FOG__
// Provides functions for atmospheric scattering and aerial perspective.
//
// Explanations:
// Scale Height = the altitude (height above ground) at which the average
// atmospheric density is found.
// Optical Depth = also called optical length, airmass, etc.
//
// References:
// [GPUGems2] GPU Gems 2: Accurate Atmospheric Scattering by Sean O'Neil.
// [GPUPro3] An Approximation to the Chapman Grazing-Incidence Function for
// Atmospheric Scattering, GPU Pro3, pp. 105.
// Papers bei Bruneton, Nishita, etc.
//
// This code contains embedded portions of free sample source code from
// http://www-evasion.imag.fr/Membres/Eric.Bruneton/PrecomputedAtmosphericScattering2.zip, Author: Eric Bruneton,
// 08/16/2011, Copyright (c) 2008 INRIA, All Rights Reserved, which have been altered from their original version.
#include "./Flax/Atmosphere.hlsl"
static const float HeightOffset = 0.01f;
// inscattered light along ray x+tv, when sun in direction s (=S[L]-T(x,x0)S[L]|x0)
float3 GetInscatterColor(float fogDepth, float3 X, float T, float3 V, float3 S, float radius, float Mu, out float3 attenuation, bool isSceneGeometry)
{
float3 result = float3(0.0f, 0.0f, 0.0f);
attenuation = float3(1.0f, 1.0f, 1.0f);
float d = -radius * Mu - sqrt(radius * radius * (Mu * Mu - 1.0) + RadiusAtmosphere * RadiusAtmosphere);
if (d > 0.0f)
{
// if X in space and ray intersects atmosphere
// move X to nearest intersection of ray with top atmosphere boundary
X += d * V;
T -= d;
Mu = (radius * Mu + d) / RadiusAtmosphere;
radius = RadiusAtmosphere;
}
float epsilon = 0.005f;
if (radius < RadiusGround + HeightOffset + epsilon)
{
float diff = (RadiusGround + HeightOffset + epsilon) - radius;
X -= diff * V;
T -= diff;
radius = RadiusGround + HeightOffset + epsilon;
Mu = dot(X, V) / radius;
}
if (radius <= RadiusAtmosphere && fogDepth > 0.0f)
{
float3 X0 = X + T * V;
float R0 = length(X0);
float Nu = dot(V, S);
float MuS = dot(X, S) / radius;
float MuHorizon = -sqrt(1.0 - (RadiusGround / radius) * (RadiusGround / radius));
if (isSceneGeometry)
{
Mu = max(Mu, MuHorizon + epsilon + 0.15f);
}
else
{
Mu = max(Mu, MuHorizon + epsilon);
}
float MuOriginal = Mu;
float blendRatio = 0.0f;
if (isSceneGeometry)
{
blendRatio = saturate(exp(-V.z) - 0.5);
if (blendRatio < 1.0)
{
V.z = max(V.z, 0.15);
V = normalize(V);
float3 x1 = X + T * V;
Mu = dot(x1, V) / length(x1);
}
}
float phaseR = PhaseFunctionR(Nu);
float phaseM = PhaseFunctionM(Nu);
float4 inscatter = max(Texture4DSample(AtmosphereInscatterTexture, radius, Mu, MuS, Nu), 0.0);
if (T > 0.0)
{
attenuation = AnalyticTransmittance(radius, Mu, T);
float Mu0 = dot(X0, V) / R0;
float MuS0 = dot(X0, S) / R0;
if (isSceneGeometry)
{
R0 = max(R0, radius);
}
if (R0 > RadiusGround + HeightOffset)
{
if (blendRatio < 1.0)
{
inscatter = max(inscatter - attenuation.rgbr * Texture4DSample(AtmosphereInscatterTexture, R0, Mu0, MuS0, Nu), 0.0);
if (!isSceneGeometry)
{
if (abs(Mu - MuHorizon) < epsilon)
{
Mu = MuHorizon - epsilon;
R0 = sqrt(radius * radius + T * T + 2.0 * radius * T * Mu);
Mu0 = (radius * Mu + T) / R0;
Mu0 = max(MuHorizon + epsilon, Mu0);
float4 inscatter0 = Texture4DSample(AtmosphereInscatterTexture, radius, Mu, MuS, Nu);
float4 inscatter1 = Texture4DSample(AtmosphereInscatterTexture, R0, Mu0, MuS0, Nu);
float4 inscatterA = max(inscatter0 - attenuation.rgbr * inscatter1, 0.0);
Mu = MuHorizon + epsilon;
R0 = sqrt(radius * radius + T * T + 2.0 * radius * T * Mu);
Mu0 = (radius * Mu + T) / R0;
Mu0 = max(MuHorizon + epsilon, Mu0);
inscatter0 = Texture4DSample(AtmosphereInscatterTexture, radius, Mu, MuS, Nu);
inscatter1 = Texture4DSample(AtmosphereInscatterTexture, R0, Mu0, MuS0, Nu);
float4 inscatterB = max(inscatter1 - attenuation.rgbr * inscatter1, 0.0);
float alpha = ((Mu - MuHorizon) + epsilon) * 0.5f / epsilon;
inscatter = lerp(inscatterA, inscatterB, alpha);
}
}
else if (blendRatio > 0.0)
{
inscatter = lerp(inscatter, (1.0 - attenuation.rgbr) * max(Texture4DSample(AtmosphereInscatterTexture, radius, MuOriginal, MuS, Nu), 0.0), blendRatio);
}
}
else
{
inscatter = (1.0 - attenuation.rgbr) * inscatter;
}
}
}
inscatter.w *= smoothstep(0.00, 0.02, MuS);
result = max(inscatter.rgb * phaseR + GetMie(inscatter) * phaseM, 0.0);
}
return result;
}
// Ground radiance at end of ray x+tv, when sun in direction s attenuated between ground and viewer (=R[L0]+R[L*])
float3 GetGroundColor(float4 sceneColor, float3 X, float T, float3 V, float3 S, float radius, float3 attenuation, bool isSceneGeometry)
{
float3 result = float3(0.0f, 0.0f, 0.0f);
if (T > 0.0f)
{
// if ray hits ground surface
// ground Reflectance at end of ray, X0
float3 X0 = X + T * V;
float R0 = length(X0);
float3 N = X0 / R0;
sceneColor.xyz = saturate(sceneColor.xyz + 0.05);
float4 reflectance = sceneColor * float4(0.2, 0.2, 0.2, 1.0);
// direct sun light (radiance) reaching X0
float MuS = dot(N, S);
float3 sunLight = isSceneGeometry ? float3(0.f, 0.f, 0.f) : TransmittanceWithShadow(R0, MuS);
// precomputed sky light (irradiance) (=E[L*]) at X0
float3 groundSkyLight = Irradiance(AtmosphereIrradianceTexture, R0, MuS);
// light reflected at X0 (=(R[L0]+R[L*])/T(X,X0))
float3 groundColor = reflectance.rgb * ((max(MuS, 0.0) * sunLight + groundSkyLight) / PI);
// water specular color due to SunLight
if (!isSceneGeometry && reflectance.w > 0.0)
{
float3 H = normalize(S - V);
float fresnel = 0.02 + 0.98 * pow(1.0 - dot(-V, H), 5.0);
float waterBrdf = fresnel * pow(max(dot(H, N), 0.0), 150.0);
groundColor += reflectance.w * max(waterBrdf, 0.0) * sunLight;
}
result = attenuation * groundColor; //=R[L0]+R[L*]
}
return result;
}
// Direct sun light for ray x+tv, when sun in direction s (=L0)
float3 GetSunColor(AtmosphericFogData atmosphericFog, float3 X, float T, float3 V, float3 S, float radius, float Mu)
{
if (T > 0.0)
return float3(0.0f, 0.0f, 0.0f);
float3 transmittance = radius <= RadiusAtmosphere ? TransmittanceWithShadow(radius, Mu) : float3(1.0, 1.0, 1.0); // T(X,xo)
float sunIntensity = step(cos(PI * atmosphericFog.AtmosphericFogSunDiscScale / 180.0), dot(V, S)); // Lsun
return transmittance * sunIntensity; // Eq (9)
}
float3 inscatter(inout float3 x, inout float t, float3 v, float3 s, out float r, out float mu, out float3 attenuation)
{
float3 result = 0;
r = length(x);
mu = dot(x, v) / r;
float d = -r * mu - sqrt(r * r * (mu * mu - 1.0) + RadiusAtmosphere * RadiusAtmosphere);
// if x in space and ray intersects atmosphere
if (d > 0.0)
{
// move x to nearest intersection of ray with top atmosphere boundary
x += d * v;
t -= d;
mu = (r * mu + d) / RadiusAtmosphere;
r = RadiusAtmosphere;
}
float epsilon = 0.0045f;
// if ray intersects atmosphere
if (r <= RadiusAtmosphere)
{
float nu = dot(v, s);
float muS = dot(x, s) / r;
float phaseR = PhaseFunctionR(nu);
float phaseM = PhaseFunctionM(nu);
float4 inscatter = max(Texture4DSample(AtmosphereInscatterTexture, r, mu, muS, nu), 0.0);
if (t > 0.0)
{
float3 x0 = x + t * v;
float r0 = length(x0);
float rMu0 = dot(x0, v);
float mu0 = rMu0 / r0;
float muS0 = dot(x0, s) / r0;
attenuation = AnalyticTransmittance(r, mu, t);
if (r0 > RadiusGround + 0.01)
{
// computes S[L]-T(x,x0)S[L]|x0
inscatter = max(inscatter - attenuation.rgbr * Texture4DSample(AtmosphereInscatterTexture, r0, mu0, muS0, nu), 0.0);
float muHoriz = -sqrt(1.0 - (RadiusGround / r) * (RadiusGround / r));
if (abs(mu - muHoriz) < epsilon)
{
mu = muHoriz - epsilon;
r0 = sqrt(r * r + t * t + 2.0 * r * t * mu);
mu0 = (r * mu + t) / r0;
float4 inScatter0 = Texture4DSample(AtmosphereInscatterTexture, r, mu, muS, nu);
float4 inScatter1 = Texture4DSample(AtmosphereInscatterTexture, r0, mu0, muS0, nu);
float4 inScatterA = max(inScatter0 - attenuation.rgbr * inScatter1, 0.0);
mu = muHoriz + epsilon;
r0 = sqrt(r * r + t * t + 2.0 * r * t * mu);
mu0 = (r * mu + t) / r0;
inScatter0 = Texture4DSample(AtmosphereInscatterTexture, r, mu, muS, nu);
inScatter1 = Texture4DSample(AtmosphereInscatterTexture, r0, mu0, muS0, nu);
float4 inScatterB = max(inScatter0 - attenuation.rgbr * inScatter1, 0.0);
float alpha = ((mu - muHoriz) + epsilon) / (2.0 * epsilon);
inscatter = lerp(inScatterA, inScatterB, alpha);
}
}
}
inscatter.w *= smoothstep(0.00, 0.02, muS);
result = max(inscatter.rgb * phaseR + GetMie(inscatter) * phaseM, 0.0);
}
return result;
}
static const float EPSILON_ATMOSPHERE = 0.002f;
static const float EPSILON_INSCATTER = 0.004f;
// input - viewPosition: view position in world space
// input - d: view ray in world space
// output - offset: distance to atmosphere or 0 if within atmosphere
// output - maxPathLength: distance traversed within atmosphere
// output - return value: intersection occurred true/false
bool intersectAtmosphere(in float3 viewPosition, in float3 d, out float offset, out float maxPathLength)
{
offset = 0.0f;
maxPathLength = 0.0f;
// vector from ray origin to center of the sphere
float3 l = -viewPosition;
float l2 = dot(l, l);
float s = dot(l, d);
// adjust top atmosphere boundary by small epsilon to prevent artifacts
float r = Rt - EPSILON_ATMOSPHERE;
float r2 = r * r;
if (l2 <= r2)
{
// ray origin inside sphere, hit is ensured
float m2 = l2 - (s * s);
float q = sqrt(r2 - m2);
maxPathLength = s + q;
return true;
}
if (s >= 0)
{
// ray starts outside in front of sphere, hit is possible
float m2 = l2 - (s * s);
if (m2 <= r2)
{
// ray hits atmosphere definitely
float q = sqrt(r2 - m2);
offset = s - q;
maxPathLength = (s + q) - offset;
return true;
}
}
return false;
}
// input - surfacePos: reconstructed position of current pixel
// input - viewDir: view ray in world space
// input/output - attenuation: extinction factor along view path
// input/output - irradianceFactor: surface hit within atmosphere 1.0f, otherwise 0.0f
// output - return value: total in-scattered light
float3 GetInscatteredLight(AtmosphericFogData atmosphericFog, in float3 viewPosition, in float3 surfacePos, in float3 viewDir, in out float3 attenuation, in out float irradianceFactor)
{
float3 inscatteredLight = float3(0.0f, 0.0f, 0.0f);
float offset;
float maxPathLength;
if (intersectAtmosphere(viewPosition, viewDir, offset, maxPathLength))
{
return float3(offset / 10, 0, 0);
float pathLength = distance(viewPosition, surfacePos);
//return pathLength.xxx;
// check if object occludes atmosphere
if (pathLength > offset)
{
//return float3(1,0,0);
// offsetting camera
float3 startPos = viewPosition + offset * viewDir;
float startPosHeight = length(startPos);
pathLength -= offset;
// starting position of path is now ensured to be inside atmosphere
// was either originally there or has been moved to top boundary
float muStartPos = dot(startPos, viewDir) / startPosHeight;
float nuStartPos = dot(viewDir, atmosphericFog.AtmosphericFogSunDirection);
float musStartPos = dot(startPos, atmosphericFog.AtmosphericFogSunDirection) / startPosHeight;
// in-scattering for infinite ray (light in-scattered when
// no surface hit or object behind atmosphere)
float4 inscatter = max(Texture4DSample(AtmosphereInscatterTexture, startPosHeight, muStartPos, musStartPos, nuStartPos), 0.0f);
float surfacePosHeight = length(surfacePos);
float musEndPos = dot(surfacePos, atmosphericFog.AtmosphericFogSunDirection) / surfacePosHeight;
//return float3(muStartPos, 0, 0);
// check if object hit is inside atmosphere
if (pathLength < maxPathLength)
{
//return float3(1,0,0);
// reduce total in-scattered light when surface hit
// within atmosphere
// fíx described in chapter 5.1.1
attenuation = AnalyticTransmittance(startPosHeight, muStartPos, pathLength); // TODO: optimize this!! here
float muEndPos = dot(surfacePos, viewDir) / surfacePosHeight;
float4 inscatterSurface = Texture4DSample(AtmosphereInscatterTexture, surfacePosHeight, muEndPos, musEndPos, nuStartPos);
inscatter = max(inscatter - attenuation.rgbr * inscatterSurface, 0.0f);
irradianceFactor = 1.0f;
}
else
{
//return float3(1,0,0);
// retrieve extinction factor for inifinte ray
// fíx described in chapter 5.1.1
attenuation = AnalyticTransmittance(startPosHeight, muStartPos, pathLength); // TODO: optimize this! and here
}
/*
// avoids imprecision problems near horizon by interpolating between
// two points above and below horizon
// fíx described in chapter 5.1.2
float muHorizon = -sqrt(1.0 - (g_Rg / startPosHeight) * (g_Rg / startPosHeight));
if (abs(muStartPos - muHorizon) < EPSILON_INSCATTER)
{
float mu = muHorizon - EPSILON_INSCATTER;
float samplePosHeight = sqrt(startPosHeight*startPosHeight + pathLength * pathLength + 2.0f * startPosHeight * pathLength*mu);
float muSamplePos = (startPosHeight * mu + pathLength) / samplePosHeight;
float4 inScatter0 = Texture4DSample(AtmosphereInscatterTexture, startPosHeight, mu, musStartPos, nuStartPos);
float4 inScatter1 = Texture4DSample(AtmosphereInscatterTexture, samplePosHeight, muSamplePos, musEndPos, nuStartPos);
float4 inScatterA = max(inScatter0 - attenuation.rgbr * inScatter1, 0.0);
mu = muHorizon + EPSILON_INSCATTER;
samplePosHeight = sqrt(startPosHeight * startPosHeight + pathLength * pathLength + 2.0f * startPosHeight * pathLength * mu);
muSamplePos = (startPosHeight * mu + pathLength) / samplePosHeight;
inScatter0 = Texture4DSample(AtmosphereInscatterTexture, startPosHeight, mu, musStartPos, nuStartPos);
inScatter1 = Texture4DSample(AtmosphereInscatterTexture, samplePosHeight, muSamplePos, musEndPos, nuStartPos);
float4 inScatterB = max(inScatter0 - attenuation.rgbr * inScatter1, 0.0f);
float t = ((muStartPos - muHorizon) + EPSILON_INSCATTER) / (2.0 * EPSILON_INSCATTER);
inscatter = lerp(inScatterA, inScatterB, t);
}
*/
// avoids imprecision problems in Mie scattering when sun is below
//horizon
// fíx described in chapter 5.1.3
inscatter.w *= smoothstep(0.00, 0.02, musStartPos);
float phaseR = PhaseFunctionR(nuStartPos);
float phaseM = PhaseFunctionM(nuStartPos);
inscatteredLight = max(inscatter.rgb * phaseR + GetMie(inscatter) * phaseM, 0.0f);
float sunIntensity = 10;
inscatteredLight *= sunIntensity;
}
}
return inscatteredLight;
}
float4 GetAtmosphericFog(AtmosphericFogData atmosphericFog, float viewFar, float3 viewPosition, float3 viewVector, float sceneDepth, float3 sceneColor)
{
float4 result = float4(1, 0, 0, 1);
#if 0
float scale = 0.00001f * atmosphericFog.AtmosphericFogDistanceScale;// convert cm to km
viewPosition.y = (viewPosition.y - atmosphericFog.AtmosphericFogGroundOffset) * atmosphericFog.AtmosphericFogAltitudeScale;
//viewPosition *= atmosphericFog.AtmosphericFogDistanceScale;
//viewPosition *= scale;
//viewPosition *= scale;
//if(length(worldPosition) > Rg)
// return float4(0, 0,0 ,0);
float3 attenuation = float3(1.0f, 1.0f, 1.0f);
float irradianceFactor = 0.0f;
float3 inscatteredLight = GetInscatteredLight(atmosphericFog, viewPosition, worldPosition, viewVector, attenuation, irradianceFactor);
//float3 inscatteredLight = GetInscatteredLight(atmosphericFog, worldPosition, viewPosition, viewVector, attenuation, irradianceFactor);
return float4(inscatteredLight, 1);
/*float offset = 0.0f;
float maxPathLength = 0.0f;
if(intersectAtmosphere(worldPosition, viewVector, offset, maxPathLength))
{
return float4(1,0,0,1);
}*/
return float4(0, 0, 0, 1);
#endif
#if 1
// TODO: scale viewPosition from cm to km
float scale = 0.0001f * atmosphericFog.AtmosphericFogDistanceScale;
viewPosition.y = (viewPosition.y - atmosphericFog.AtmosphericFogGroundOffset) * atmosphericFog.AtmosphericFogAltitudeScale;
//viewPosition *= atmosphericFog.AtmosphericFogDistanceScale;
viewPosition *= scale;
//viewPosition.xz *= 0.00001f;
//viewPosition *= scale;
viewPosition.y += RadiusGround + HeightOffset;
//viewPosition.y -= atmosphericFog.AtmosphericFogGroundOffset;
//worldPosition
float Radius = length(viewPosition);
float3 V = normalize(viewVector);
float Mu = dot(viewPosition, V) / Radius;
float T = -Radius * Mu - sqrt(Radius * Radius * (Mu * Mu - 1.0) + RadiusGround * RadiusGround);
//-Radius * Mu > sqrt(Radius * Radius * (Mu * Mu - 1.0) + RadiusGround * RadiusGround)
/*
float3 g = viewPosition - float3(0.0, 0.0, RadiusGround + 10.0);
float a = V.x * V.x + V.y * V.y - V.z * V.z;
float b = 2.0 * (g.x * V.x + g.y * V.y - g.z * V.z);
float c = g.x * g.x + g.y * g.y - g.z * g.z;
float d = -(b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
bool cone = d > 0.0 && abs(viewPosition.z + d * V.z - RadiusGround) <= 10.0;
if (T > 0.0)
{
if (cone && d < T)
T = d;
//return float4(1,0,0,1);
}
else if (cone)
{
T = d;
//return float4(0,0,1,1);
}
*/
//return float4(V, 1);
// TODO: create uniform param for that depth limit
float depthThreshold = min(100.0f * atmosphericFog.AtmosphericFogDistanceScale, (viewFar * 0.997f) * scale); // 100km limit or far plane
sceneDepth *= scale;
float fogDepth = max(0.0f, sceneDepth - atmosphericFog.AtmosphericFogStartDistance);
float shadowFactor = 1.0f;
float DistanceRatio = min(fogDepth * 0.1f / atmosphericFog.AtmosphericFogStartDistance, 1.0f);
bool isSceneGeometry = (sceneDepth < depthThreshold); // Assume as scene geometry
//isSceneGeometry = false;
if (isSceneGeometry)
{
shadowFactor = DistanceRatio * atmosphericFog.AtmosphericFogPower;
T = max(sceneDepth + atmosphericFog.AtmosphericFogDistanceOffset, 1.0f);
// TODO: fog for scene objects
return 0;
//return float4(0, 1, 0, 1);
}
float3 attenuation;
float3 inscatterColor = GetInscatterColor(fogDepth, viewPosition, T, V, atmosphericFog.AtmosphericFogSunDirection, Radius, Mu, attenuation, isSceneGeometry); //S[L]-T(viewPosition,xs)S[l]|xs
//float3 inscatterColor = inscatter(viewPosition, T, V, atmosphericFog.AtmosphericFogSunDirection, Radius, Mu, attenuation); //S[L]-T(viewPosition,xs)S[l]|xs
//float3 groundColor = 0;
float3 groundColor = GetGroundColor(float4(sceneColor.rgb, 1.f), viewPosition, T, V, atmosphericFog.AtmosphericFogSunDirection, Radius, attenuation, isSceneGeometry); //R[L0]+R[L*]
float3 sun = GetSunColor(atmosphericFog, viewPosition, T, V, atmosphericFog.AtmosphericFogSunDirection, Radius, Mu); //L0
float3 color = (atmosphericFog.AtmosphericFogSunPower) * (sun + groundColor + inscatterColor) * atmosphericFog.AtmosphericFogSunColor.rgb;
//return float4(color, lerp(saturate(attenuation.r * atmosphericFog.AtmosphericFogDensityScale - atmosphericFog.AtmosphericFogDensityOffset), 1.0f, (1.0f - DistanceRatio)) );
//return lerp(saturate(attenuation.r * atmosphericFog.AtmosphericFogDensityScale - atmosphericFog.AtmosphericFogDensityOffset), 1.0f, (1.0f - DistanceRatio));
//color *= lerp(saturate(attenuation.r * atmosphericFog.AtmosphericFogDensityScale - atmosphericFog.AtmosphericFogDensityOffset), 1.0f, (1.0f - DistanceRatio));
//return attenuation.bbbb;
return float4(color, 1);
//return float4(sun, 1);
//return float4(inscatterColor + sun, 1);
//return float4(sun + groundColor, 1);
//return float4(sun + groundColor + inscatterColor, 1);
#endif
return result;
}
float4 GetAtmosphericFog(AtmosphericFogData atmosphericFog, float viewFar, float3 worldPosition, float3 cameraPosition)
{
float3 viewVector = worldPosition - cameraPosition;
return GetAtmosphericFog(atmosphericFog, viewFar, cameraPosition, viewVector, length(viewVector), float3(0, 0, 0));
}
#endif