From 9b9ae5cf38c1e2935e8682908c7e5292eb461205 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fernando=20Garc=C3=ADa=20Li=C3=B1=C3=A1n?= Date: Wed, 28 Jul 2021 09:40:04 +0200 Subject: [PATCH] HDR pipeline: even better atmospheric scattering MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implementation of "A Scalable and Production Ready Sky and Atmosphere Rendering Technique" by Sébastien Hillaire (2020). --- Compositor/HDR/hdr.xml | 103 ++++++++++++---- Effects/HDR/atmos-aerial-perspective.eff | 12 +- Effects/HDR/atmos-multiple-scattering.eff | 18 +++ Effects/HDR/atmos-sky-view.eff | 12 +- Effects/HDR/atmos-transmittance.eff | 13 ++ Effects/HDR/lighting.eff | 7 +- Effects/skydome.eff | 10 ++ Shaders/HDR/atmos-aerial-perspective.frag | 127 +++++++++++++++---- Shaders/HDR/atmos-include.frag | 78 ++++++++++++ Shaders/HDR/atmos-multiple-scattering.frag | 135 +++++++++++++++++++++ Shaders/HDR/atmos-sky-view.frag | 109 +++++++++++++---- Shaders/HDR/atmos-transmittance.frag | 55 +++++++++ Shaders/HDR/atmosphere-include.frag | 133 -------------------- Shaders/HDR/lighting.frag | 48 ++++---- Shaders/HDR/skydome.frag | 18 ++- Shaders/HDR/skydome.vert | 7 +- 16 files changed, 647 insertions(+), 238 deletions(-) create mode 100644 Effects/HDR/atmos-multiple-scattering.eff create mode 100644 Effects/HDR/atmos-transmittance.eff create mode 100644 Shaders/HDR/atmos-include.frag create mode 100644 Shaders/HDR/atmos-multiple-scattering.frag create mode 100644 Shaders/HDR/atmos-transmittance.frag delete mode 100644 Shaders/HDR/atmosphere-include.frag diff --git a/Compositor/HDR/hdr.xml b/Compositor/HDR/hdr.xml index badef3b72..2da27f59b 100644 --- a/Compositor/HDR/hdr.xml +++ b/Compositor/HDR/hdr.xml @@ -65,6 +65,28 @@ + + transmittance + 2d + 256 + 64 + rgb16f + linear + linear + clamp-to-edge + clamp-to-edge + + + multiple-scattering + 2d + 32 + 32 + rgb16f + linear + linear + clamp-to-edge + clamp-to-edge + sky-view 2d @@ -77,18 +99,7 @@ clamp-to-edge - aerial-inscatter - 2d - 512 - 32 - rgba16f - linear - linear - clamp-to-edge - clamp-to-edge - - - aerial-transmittance + aerial-perspective 2d 512 32 @@ -213,28 +224,70 @@ - + + + atmos-transmittance + quad + Effects/HDR/atmos-transmittance + + color + transmittance + + + + atmos-multiple-scattering + quad + Effects/HDR/atmos-multiple-scattering + + 0 + transmittance + + + color + multiple-scattering + + atmos-sky-view quad Effects/HDR/atmos-sky-view + + 0 + transmittance + + + 1 + multiple-scattering + color sky-view - atmos-aerial-perspective quad Effects/HDR/atmos-aerial-perspective + + 0 + transmittance + + + 1 + multiple-scattering + - color0 - aerial-inscatter - - - color1 - aerial-transmittance + color + aerial-perspective @@ -617,11 +670,7 @@ 11 - aerial-inscatter - - - 12 - aerial-transmittance + aerial-perspective color0 @@ -648,6 +697,10 @@ 11 sky-view + + 13 + transmittance + color hdr-result diff --git a/Effects/HDR/atmos-aerial-perspective.eff b/Effects/HDR/atmos-aerial-perspective.eff index 727714dac..e07dd4495 100644 --- a/Effects/HDR/atmos-aerial-perspective.eff +++ b/Effects/HDR/atmos-aerial-perspective.eff @@ -6,9 +6,19 @@ Shaders/HDR/trivial.vert Shaders/HDR/atmos-aerial-perspective.frag - Shaders/HDR/atmosphere-include.frag + Shaders/HDR/atmos-include.frag Shaders/HDR/gbuffer-include.frag + + transmittance_lut + sampler-2d + 0 + + + multiscattering_lut + sampler-2d + 1 + diff --git a/Effects/HDR/atmos-multiple-scattering.eff b/Effects/HDR/atmos-multiple-scattering.eff new file mode 100644 index 000000000..8f55d189c --- /dev/null +++ b/Effects/HDR/atmos-multiple-scattering.eff @@ -0,0 +1,18 @@ + + + Effects/HDR/atmos-multiple-scattering + + + + Shaders/HDR/trivial.vert + Shaders/HDR/atmos-multiple-scattering.frag + Shaders/HDR/atmos-include.frag + + + transmittance_lut + sampler-2d + 0 + + + + diff --git a/Effects/HDR/atmos-sky-view.eff b/Effects/HDR/atmos-sky-view.eff index b7ea7f888..0c9cb3ae5 100644 --- a/Effects/HDR/atmos-sky-view.eff +++ b/Effects/HDR/atmos-sky-view.eff @@ -6,8 +6,18 @@ Shaders/HDR/trivial.vert Shaders/HDR/atmos-sky-view.frag - Shaders/HDR/atmosphere-include.frag + Shaders/HDR/atmos-include.frag + + transmittance_lut + sampler-2d + 0 + + + multiscattering_lut + sampler-2d + 1 + diff --git a/Effects/HDR/atmos-transmittance.eff b/Effects/HDR/atmos-transmittance.eff new file mode 100644 index 000000000..76c8179d8 --- /dev/null +++ b/Effects/HDR/atmos-transmittance.eff @@ -0,0 +1,13 @@ + + + Effects/HDR/atmos-transmittance + + + + Shaders/HDR/trivial.vert + Shaders/HDR/atmos-transmittance.frag + Shaders/HDR/atmos-include.frag + + + + diff --git a/Effects/HDR/lighting.eff b/Effects/HDR/lighting.eff index d62e28ac9..09264917e 100644 --- a/Effects/HDR/lighting.eff +++ b/Effects/HDR/lighting.eff @@ -68,15 +68,10 @@ 10 - aerial_inscatter_lut + aerial_perspective_lut sampler-2d 11 - - aerial_transmittance_lut - sampler-2d - 12 - diff --git a/Effects/skydome.eff b/Effects/skydome.eff index 0e0574c1c..c0237c536 100644 --- a/Effects/skydome.eff +++ b/Effects/skydome.eff @@ -319,6 +319,11 @@ sampler-2d 11 + + transmittance_lut + sampler-2d + 13 + @@ -340,6 +345,11 @@ sampler-2d 11 + + transmittance_lut + sampler-2d + 13 + diff --git a/Shaders/HDR/atmos-aerial-perspective.frag b/Shaders/HDR/atmos-aerial-perspective.frag index 41ab755b5..5ec3ac62b 100644 --- a/Shaders/HDR/atmos-aerial-perspective.frag +++ b/Shaders/HDR/atmos-aerial-perspective.frag @@ -1,7 +1,18 @@ +// An implementation of Sébastien Hillaire's "A Scalable and Production Ready +// Sky and Atmosphere Rendering Technique". +// +// This shader generates the aerial perspective LUT. This LUT is used by opaque +// and transparent objects to apply atmospheric scattering. In-scattering is +// stored in the RGB channels, while transmittance is stored in the alpha +// channel. +// Unlike the paper, we are using a tiled 2D texture instead of a true 3D +// texture. For some reason the overhead of rendering to a texture a lot of +// times (the depth of the 3D texture) seems to be too high, probably because +// OSG is not sharing state between those passes. + #version 330 core -layout(location = 0) out vec3 out_inscatter; -layout(location = 1) out vec3 out_transmittance; +out vec4 fragColor; in vec2 texCoord; @@ -9,39 +20,111 @@ uniform mat4 fg_ViewMatrixInverse; uniform vec3 fg_CameraPositionCart; uniform vec3 fg_CameraPositionGeod; uniform vec3 fg_SunDirectionWorld; -uniform vec2 fg_NearFar; +uniform sampler2D transmittance_lut; +uniform sampler2D multiscattering_lut; + +const float PI = 3.141592653; +const float ATMOSPHERE_RADIUS = 6471e3; const float TOTAL_SLICES = 16.0; const float DEPTH_RANGE = 32000.0; +const int AERIAL_PERSPECTIVE_SAMPLES = 20; +const vec3 ONE_OVER_THREE = vec3(1.0 / 3.0); vec3 positionFromDepth(vec2 pos, float depth); -void calculateScattering(in vec3 rayOrigin, - in vec3 rayDir, - in float tmax, - in vec3 lightDir, - in float earthRadius, - out vec3 inscatter, - out vec3 transmittance); + +float raySphereIntersection(vec3 ro, vec3 rd, float radius); +vec3 sampleMedium(in float height, + out float mieScattering, out float mieAbsorption, + out vec3 rayleighScattering, out vec3 ozoneAbsorption); +float miePhaseFunction(float cosTheta); +float rayleighPhaseFunction(float cosTheta); +vec3 getValueFromLUT(sampler2D lut, float sunCosTheta, float normalizedHeight); void main() { + vec3 up = normalize(fg_CameraPositionCart); + float sunCosTheta = dot(fg_SunDirectionWorld, up); + + // Account for the depth layer we are currently in + // FIXME: We should probably be writing the pixel center float x = texCoord.x * TOTAL_SLICES; vec2 coord = vec2(fract(x), texCoord.y); - - float depth = mix(fg_NearFar.x, DEPTH_RANGE, ceil(x) / TOTAL_SLICES); + // Depth goes from the 0 to DEPTH_RANGE in a squared distribution. + // The first slice is not at 0 since that would waste a slice. + float w = ceil(x) / TOTAL_SLICES; + w *= w; + float depth = w * DEPTH_RANGE; vec3 fragPos = positionFromDepth(coord * 2.0 - 1.0, 0.0); vec3 rayDir = vec4(fg_ViewMatrixInverse * vec4(normalize(fragPos), 0.0)).xyz; - float earthRadius = length(fg_CameraPositionCart) - fg_CameraPositionGeod.z; + float cameraHeight = length(fg_CameraPositionCart); + float earthRadius = cameraHeight - fg_CameraPositionGeod.z; - vec3 inscatter, transmittance; - calculateScattering(fg_CameraPositionCart, - rayDir, - depth, - fg_SunDirectionWorld, - earthRadius, - inscatter, transmittance); - out_inscatter = inscatter; - out_transmittance = transmittance; + vec3 rayOrigin = fg_CameraPositionCart; + + float atmosDist = raySphereIntersection(rayOrigin, rayDir, ATMOSPHERE_RADIUS); + float groundDist = raySphereIntersection(rayOrigin, rayDir, earthRadius); + + float tmax; + if (cameraHeight < ATMOSPHERE_RADIUS) { + // We are inside the atmosphere + if (groundDist < 0.0) { + // No ground collision, use the distance to the outer atmosphere + tmax = atmosDist; + } else { + // Use the distance to the ground + tmax = groundDist; + } + } else { + // We are in outer space, skip + fragColor = vec4(0.0, 0.0, 0.0, 1.0); + return; + } + // Clip the max distance to the depth of this slice + tmax = min(tmax, depth); + + float cosTheta = dot(rayDir, fg_SunDirectionWorld); + float miePhase = miePhaseFunction(cosTheta); + float rayleighPhase = rayleighPhaseFunction(cosTheta); + + vec3 L = vec3(0.0); + vec3 throughput = vec3(1.0); + float t = 0.0; + + for (int i = 0; i < AERIAL_PERSPECTIVE_SAMPLES; ++i) { + float newT = ((float(i) + 0.3) / AERIAL_PERSPECTIVE_SAMPLES) * tmax; + float dt = newT - t; + t = newT; + + vec3 samplePos = rayOrigin + rayDir * t; + float height = length(samplePos) - earthRadius; + float normalizedHeight = height / (ATMOSPHERE_RADIUS - earthRadius); + + float mieScattering, mieAbsorption; + vec3 rayleighScattering, ozoneAbsorption; + vec3 extinction = sampleMedium(height, mieScattering, mieAbsorption, + rayleighScattering, ozoneAbsorption); + + vec3 sampleTransmittance = exp(-dt*extinction); + + vec3 sunTransmittance = getValueFromLUT( + transmittance_lut, sunCosTheta, normalizedHeight); + vec3 multiscattering = getValueFromLUT( + multiscattering_lut, sunCosTheta, normalizedHeight); + + vec3 S = + rayleighScattering * (rayleighPhase * sunTransmittance + multiscattering) + + mieScattering * (miePhase * sunTransmittance + multiscattering); + + vec3 Sint = (S - S * sampleTransmittance) / extinction; + L += throughput * Sint; + throughput *= sampleTransmittance; + } + + // Instead of storing an entire vec3, store the mean of its components + float transmittance = dot(throughput, ONE_OVER_THREE); + + fragColor = vec4(L, transmittance); } diff --git a/Shaders/HDR/atmos-include.frag b/Shaders/HDR/atmos-include.frag new file mode 100644 index 000000000..59558b74d --- /dev/null +++ b/Shaders/HDR/atmos-include.frag @@ -0,0 +1,78 @@ +// An implementation of Sébastien Hillaire's "A Scalable and Production Ready +// Sky and Atmosphere Rendering Technique". + +#version 330 core + +const float PI = 3.141592653; + +// Atmosphere parameters +// Section 2.1 of [Bruneton08], units are inverse meters +const float mie_density_height_scale = 8.33333e-4; // Hm=1.2km +const float rayleigh_density_height_scale = 1.25e-4; // Hr=8km +const float mie_scattering = 3.996e-6; +const float mie_absorption = 4.4e-6; +const vec3 rayleigh_scattering = vec3(5.802, 13.558, 33.1) * vec3(1e-6); +const vec3 ozone_absorption = vec3(0.650, 1.881, 0.085) * vec3(1e-6); + +const float g = 0.8; +const float gg = g*g; +const float mie_phase_scale = 3.0/(8.0*PI); +const float rayleigh_phase_scale = 3.0/(16.0*PI); + +// Returns the distance between ro and the first intersection with the sphere +// or -1.0 if there is no intersection. +// -1.0 is also returned if the ray is pointing away from the sphere. +float raySphereIntersection(vec3 ro, vec3 rd, float radius) +{ + float b = dot(ro, rd); + float c = dot(ro, ro) - radius*radius; + if (c > 0.0 && b > 0.0) return -1.0; + float d = b*b - c; + if (d < 0.0) return -1.0; + if (d > b*b) return (-b+sqrt(d)); + return (-b-sqrt(d)); +} + +// Sample the sky medium properties at a height in meters, 0 being the ground +// and ~100km being the top of the atmosphere. +// Returns the total atmospheric extinction. +vec3 sampleMedium(in float height, + out float mieScattering, out float mieAbsorption, + out vec3 rayleighScattering, out vec3 ozoneAbsorption) +{ + float densityMie = exp(-mie_density_height_scale * height); + float densityRayleigh = exp(-rayleigh_density_height_scale * height); + float densityOzone = max(0.0, 1.0 - abs(height-25.0e3)/15.0e3); + + mieScattering = mie_scattering * densityMie; + mieAbsorption = mie_absorption * densityMie; + + rayleighScattering = rayleigh_scattering * densityRayleigh; + + ozoneAbsorption = ozone_absorption * densityOzone; + + return mieScattering + mieAbsorption + rayleighScattering + ozoneAbsorption; +} + +// Approximation of the Mie phase function with the Cornette-Shanks phase function +float miePhaseFunction(float cosTheta) +{ + float num = (1.0 - gg) * (1.0 + cosTheta*cosTheta); + float den = (2.0 + gg) * pow((1.0 + gg - 2.0 * g * cosTheta), 1.5); + return mie_phase_scale * num / den; +} + +float rayleighPhaseFunction(float cosTheta) +{ + return rayleigh_phase_scale * (1.0 + cosTheta*cosTheta); +} + +// Sample one of the LUTs (transmittance or multiple scattering) for a given +// normalized height inside the atmosphere [0,1] and the cosine of the Sun +// zenith angle. +vec3 getValueFromLUT(sampler2D lut, float sunCosTheta, float normHeight) +{ + float x = clamp(sunCosTheta * 0.5 + 0.5, 0.0, 1.0); + float y = clamp(normHeight, 0.0, 1.0); + return texture(lut, vec2(x, y)).rgb; +} diff --git a/Shaders/HDR/atmos-multiple-scattering.frag b/Shaders/HDR/atmos-multiple-scattering.frag new file mode 100644 index 000000000..671a0b065 --- /dev/null +++ b/Shaders/HDR/atmos-multiple-scattering.frag @@ -0,0 +1,135 @@ +// An implementation of Sébastien Hillaire's "A Scalable and Production Ready +// Sky and Atmosphere Rendering Technique". +// +// This shader generates the multiple scattering LUT. + +#version 330 core + +out vec3 fragColor; + +in vec2 texCoord; + +uniform vec3 fg_CameraPositionCart; +uniform vec3 fg_CameraPositionGeod; + +uniform sampler2D transmittance_lut; + +const float PI = 3.141592653; +const float ATMOSPHERE_RADIUS = 6471e3; +const int SQRT_SAMPLES = 4; +const float INV_SAMPLES = 1.0 / float(SQRT_SAMPLES*SQRT_SAMPLES); +const int MULTIPLE_SCATTERING_SAMPLES = 20; + +const vec3 ground_albedo = vec3(0.3); + +float raySphereIntersection(vec3 ro, vec3 rd, float radius); +vec3 sampleMedium(in float height, + out float mieScattering, out float mieAbsorption, + out vec3 rayleighScattering, out vec3 ozoneAbsorption); +float miePhaseFunction(float cosTheta); +float rayleighPhaseFunction(float cosTheta); +vec3 getValueFromLUT(sampler2D lut, float sunCosTheta, float normalizedHeight); + +vec3 generateRayDir(float theta, float phi) +{ + float cosPhi = cos(phi); + float sinPhi = sin(phi); + float cosTheta = cos(theta); + float sinTheta = sin(theta); + return vec3(cosTheta * sinPhi, sinTheta * sinPhi, cosPhi); +} + +void main() +{ + float sunCosTheta = texCoord.x * 2.0 - 1.0; + vec3 sunDir = vec3(-sqrt(1.0 - sunCosTheta*sunCosTheta), 0.0, sunCosTheta); + + float earthRadius = length(fg_CameraPositionCart) - fg_CameraPositionGeod.z; + float altitude = mix(earthRadius, ATMOSPHERE_RADIUS, texCoord.y); + vec3 rayOrigin = vec3(0.0, 0.0, altitude); + + vec3 Ltotal = vec3(0.0); + vec3 LMStotal = vec3(0.0); + + for (int i = 0; i < SQRT_SAMPLES; ++i) { + for (int j = 0; j < SQRT_SAMPLES; ++j) { + float theta = 2.0 * PI * (float(i) + 0.5) / float(SQRT_SAMPLES); + float phi = PI * (float(j) + 0.5) / float(SQRT_SAMPLES); + vec3 rayDir = generateRayDir(theta, phi); + + float atmosDist = raySphereIntersection(rayOrigin, rayDir, ATMOSPHERE_RADIUS); + float groundDist = raySphereIntersection(rayOrigin, rayDir, earthRadius); + + float tmax; + if (groundDist < 0.0) { + // No ground collision, use the distance to the outer atmosphere + tmax = atmosDist; + } else { + // Use the distance to the ground + tmax = groundDist; + } + + float cosTheta = dot(rayDir, sunDir); + float miePhase = miePhaseFunction(cosTheta); + float rayleighPhase = rayleighPhaseFunction(-cosTheta); + + vec3 L = vec3(0.0); + vec3 LMS = vec3(0.0); + vec3 throughput = vec3(1.0); + float t = 0.0; + + for (int k = 0; k < MULTIPLE_SCATTERING_SAMPLES; ++k) { + float newT = ((float(k) + 0.3) / MULTIPLE_SCATTERING_SAMPLES) * tmax; + float dt = newT - t; + t = newT; + + vec3 samplePos = rayOrigin + rayDir * t; + float height = length(samplePos) - earthRadius; + float normalizedHeight = height / (ATMOSPHERE_RADIUS - earthRadius); + + float mieScattering, mieAbsorption; + vec3 rayleighScattering, ozoneAbsorption; + vec3 extinction = sampleMedium(height, mieScattering, mieAbsorption, + rayleighScattering, ozoneAbsorption); + + vec3 sampleTransmittance = exp(-dt*extinction); + + vec3 sunTransmittance = getValueFromLUT( + transmittance_lut, sunCosTheta, normalizedHeight); + + vec3 S = (rayleighScattering * rayleighPhase + + mieScattering * miePhase) * sunTransmittance; + + // Not using the power serie + vec3 MS = mieScattering + rayleighScattering; + vec3 MSint = (MS - MS * sampleTransmittance) / extinction; + LMS += throughput * MSint; + + vec3 Sint = (S - S * sampleTransmittance) / extinction; + L += throughput * Sint; + throughput *= sampleTransmittance; + } + + if (groundDist >= 0.0) { + // Account for bounced light off the Earth + vec3 p = rayOrigin + rayDir * groundDist; + float pHeight = length(p); + vec3 up = p / pHeight; + + float normHeight = (pHeight - earthRadius) + / (ATMOSPHERE_RADIUS - earthRadius); + float sunZenithCosTheta = dot(sunDir, up); + + vec3 transmittanceFromGround = getValueFromLUT( + transmittance_lut, sunZenithCosTheta, normHeight); + L += transmittanceFromGround * throughput + * clamp(sunZenithCosTheta, 0.0, 1.0) * ground_albedo / PI; + } + + Ltotal += L * INV_SAMPLES; + LMStotal += LMS * INV_SAMPLES; + } + } + + fragColor = Ltotal / (1.0 - LMStotal); +} diff --git a/Shaders/HDR/atmos-sky-view.frag b/Shaders/HDR/atmos-sky-view.frag index d138a540a..e0a5f6a4a 100644 --- a/Shaders/HDR/atmos-sky-view.frag +++ b/Shaders/HDR/atmos-sky-view.frag @@ -1,3 +1,12 @@ +// An implementation of Sébastien Hillaire's "A Scalable and Production Ready +// Sky and Atmosphere Rendering Technique". +// +// This shader generates the sky-view texture. Since the sky generally has low +// frequency detail, it's possible to pre-compute it on a small texture and +// sample it later when rendering the skydome. This effectively bypasses the +// need for raymarching on screen-sized textures, which is specially costly on +// larger resolutions like 4K. + #version 330 core out vec3 fragColor; @@ -8,22 +17,29 @@ uniform vec3 fg_CameraPositionCart; uniform vec3 fg_CameraPositionGeod; uniform vec3 fg_SunDirectionWorld; +uniform sampler2D transmittance_lut; +uniform sampler2D multiscattering_lut; + const float PI = 3.141592653; -void calculateScattering(in vec3 rayOrigin, - in vec3 rayDir, - in float tmax, - in vec3 lightDir, - in float earthRadius, - out vec3 inscatter, - out vec3 transmittance); +const float ATMOSPHERE_RADIUS = 6471e3; +const int SCATTERING_SAMPLES = 32; + +float raySphereIntersection(vec3 ro, vec3 rd, float radius); +vec3 sampleMedium(in float height, + out float mieScattering, out float mieAbsorption, + out vec3 rayleighScattering, out vec3 ozoneAbsorption); +float miePhaseFunction(float cosTheta); +float rayleighPhaseFunction(float cosTheta); +vec3 getValueFromLUT(sampler2D lut, float sunCosTheta, float normalizedHeight); void main() { // Always leave the sun right in the middle of the texture as the skydome // model is already being rotated. - float sunCosTheta = dot(fg_SunDirectionWorld, normalize(fg_CameraPositionCart)); - vec3 lightDir = vec3(-sqrt(1.0 - sunCosTheta*sunCosTheta), 0.0, sunCosTheta); + vec3 up = normalize(fg_CameraPositionCart); + float sunCosTheta = dot(fg_SunDirectionWorld, up); + vec3 sunDir = vec3(-sqrt(1.0 - sunCosTheta*sunCosTheta), 0.0, sunCosTheta); float azimuth = 2.0 * PI * texCoord.x; // [0, 2pi] // Apply a non-linear transformation to the elevation to dedicate more @@ -32,20 +48,67 @@ void main() float elev = l*l * sign(l) * PI * 0.5; // [-pi/2, pi/2] vec3 rayDir = vec3(cos(elev) * cos(azimuth), cos(elev) * sin(azimuth), sin(elev)); - float cameraPosLength = length(fg_CameraPositionCart); - // Since FG internally uses WG84 coordinates, we use the current Earth - // radius under the camera, which varies with the latitude. For practical - // purposes we model the Earth as a perfect sphere with this radius. - float earthRadius = cameraPosLength - fg_CameraPositionGeod.z; + float cameraHeight = length(fg_CameraPositionCart); + float earthRadius = cameraHeight - fg_CameraPositionGeod.z; - vec3 rayOrigin = vec3(0.0, 0.0, cameraPosLength); + vec3 rayOrigin = vec3(0.0, 0.0, cameraHeight); - vec3 inscatter, transmittance; - calculateScattering(rayOrigin, - rayDir, - 9.0e8, - lightDir, - earthRadius, - inscatter, transmittance); - fragColor = inscatter; + float atmosDist = raySphereIntersection(rayOrigin, rayDir, ATMOSPHERE_RADIUS); + float groundDist = raySphereIntersection(rayOrigin, rayDir, earthRadius); + + float tmax; + if (cameraHeight < ATMOSPHERE_RADIUS) { + // We are inside the atmosphere + if (groundDist < 0.0) { + // No ground collision, use the distance to the outer atmosphere + tmax = atmosDist; + } else { + // Use the distance to the ground + tmax = groundDist; + } + } else { + // We are in outer space, skip + fragColor = vec3(0.0); + return; + } + + float cosTheta = dot(rayDir, sunDir); + float miePhase = miePhaseFunction(cosTheta); + float rayleighPhase = rayleighPhaseFunction(-cosTheta); + + vec3 L = vec3(0.0); + vec3 throughput = vec3(1.0); + float t = 0.0; + + for (int i = 0; i < SCATTERING_SAMPLES; ++i) { + float newT = ((float(i) + 0.3) / SCATTERING_SAMPLES) * tmax; + float dt = newT - t; + t = newT; + + vec3 samplePos = rayOrigin + rayDir * t; + float height = length(samplePos) - earthRadius; + float normalizedHeight = height / (ATMOSPHERE_RADIUS - earthRadius); + + float mieScattering, mieAbsorption; + vec3 rayleighScattering, ozoneAbsorption; + vec3 extinction = sampleMedium(height, mieScattering, mieAbsorption, + rayleighScattering, ozoneAbsorption); + + vec3 sampleTransmittance = exp(-dt*extinction); + + vec3 sunTransmittance = getValueFromLUT( + transmittance_lut, sunCosTheta, normalizedHeight); + vec3 multiscattering = getValueFromLUT( + multiscattering_lut, sunCosTheta, normalizedHeight); + + vec3 S = + rayleighScattering * (rayleighPhase * sunTransmittance + multiscattering) + + mieScattering * (miePhase * sunTransmittance + multiscattering); + + vec3 Sint = (S - S * sampleTransmittance) / extinction; + L += throughput * Sint; + throughput *= sampleTransmittance; + } + + fragColor = L; } diff --git a/Shaders/HDR/atmos-transmittance.frag b/Shaders/HDR/atmos-transmittance.frag new file mode 100644 index 000000000..4121475f7 --- /dev/null +++ b/Shaders/HDR/atmos-transmittance.frag @@ -0,0 +1,55 @@ +// An implementation of Sébastien Hillaire's "A Scalable and Production Ready +// Sky and Atmosphere Rendering Technique". +// +// This shader generates the transmittance LUT. It stores the transmittance to +// the Sun through the atmosphere for a given Sun zenith angle and a height +// inside the atmosphere (0 being the ground). + +#version 330 core + +out vec3 fragColor; + +in vec2 texCoord; + +uniform vec3 fg_CameraPositionCart; +uniform vec3 fg_CameraPositionGeod; + +const float ATMOSPHERE_RADIUS = 6471e3; +const int TRANSMITTANCE_STEPS = 40; + +float raySphereIntersection(vec3 ro, vec3 rd, float radius); +vec3 sampleMedium(in float height, + out float mieScattering, out float mieAbsorption, + out vec3 rayleighScattering, out vec3 ozoneAbsorption); + +void main() +{ + float sunCosTheta = texCoord.x * 2.0 - 1.0; + vec3 sunDir = vec3(-sqrt(1.0 - sunCosTheta*sunCosTheta), 0.0, sunCosTheta); + + float earthRadius = length(fg_CameraPositionCart) - fg_CameraPositionGeod.z; + float altitude = mix(earthRadius, ATMOSPHERE_RADIUS, texCoord.y); + vec3 rayOrigin = vec3(0.0, 0.0, altitude); + + float dist = raySphereIntersection(rayOrigin, sunDir, ATMOSPHERE_RADIUS); + float t = 0.0; + vec3 transmittance = vec3(1.0); + + for (int i = 0; i < TRANSMITTANCE_STEPS; ++i) { + float newT = ((float(i) + 0.3) / TRANSMITTANCE_STEPS) * dist; + float dt = newT - t; + t = newT; + + vec3 samplePos = rayOrigin + sunDir * t; + float height = length(samplePos) - earthRadius; + + float mieScattering, mieAbsorption; + vec3 rayleighScattering, ozoneAbsorption; + vec3 extinction = sampleMedium(height, mieScattering, mieAbsorption, + rayleighScattering, ozoneAbsorption); + + transmittance *= exp(-dt * extinction); + } + + fragColor = transmittance; +} diff --git a/Shaders/HDR/atmosphere-include.frag b/Shaders/HDR/atmosphere-include.frag deleted file mode 100644 index 726d2c6b4..000000000 --- a/Shaders/HDR/atmosphere-include.frag +++ /dev/null @@ -1,133 +0,0 @@ -#version 330 core - -uniform vec3 beta_rayleigh = vec3(5.5e-6, 13.0e-6, 22.4e-6); -uniform float beta_mie = 21e-6; -uniform vec3 beta_absortion = vec3(2.04e-5, 4.97e-5, 1.95e-6); -uniform float beta_ambient = 0.0; -uniform float rayleigh_scale_height = 8e3; -uniform float mie_scale_height = 1.2e3; -uniform float absortion_scale_height = 30e3; -uniform float absortion_falloff = 3e3; -uniform int num_samples = 32; -uniform int num_light_samples = 4; - -const float PI = 3.141592653; -const float ATMOSPHERE_RADIUS = 6471e3; -const vec3 SUN_INTENSITY = vec3(20.0); - -vec2 raySphereIntersection(vec3 ro, vec3 rd, float radius) -{ - vec3 tc = -ro; - float b = dot(tc, rd); - float d = b*b - dot(tc, tc) + radius*radius; - if (d < 0.0) return vec2(-1.0); - float s = sqrt(d); - return vec2(b-s, b+s); -} - -void calculateScattering(in vec3 rayOrigin, - in vec3 rayDir, - in float tmax, - in vec3 lightDir, - in float earthRadius, - out vec3 inscatter, - out vec3 transmittance) -{ - vec2 hit = raySphereIntersection(rayOrigin, rayDir, ATMOSPHERE_RADIUS); - vec2 hitEarth = raySphereIntersection(rayOrigin, rayDir, earthRadius - 1.0); - if (hitEarth.y > 0.0) - tmax = max(0.0, hitEarth.x); - - float tmin = max(hit.x, 0.0); - tmax = min(hit.y, tmax); - if (tmax < 0.0) - discard; - - float stepSize = (tmax - tmin) / float(num_samples); - - const float g = 0.758; // Mie scattering direction - const float gg = g*g; - float mu = dot(rayDir, lightDir); - float mumu = mu*mu; - float phaseRayleigh = 3.0 / (50.2654824574 /* 16*PI */) * (1.0 + mumu); - float phaseMie = 3.0 / (25.1327412287 /* 8*PI */) * ((1.0 - gg) * (mumu + 1.0)) / - (pow(1.0 + gg - 2.0 * mu * g, 1.5) * (2.0 + gg)); - - float opticalDepthRayleigh = 0.0; - float opticalDepthMie = 0.0; - float opticalDepthAbsortion = 0.0; - - float primaryTime = tmin; - - vec3 extinctionFactor = vec3(0.0); - vec3 totalRayleigh = vec3(0.0); - vec3 totalMie = vec3(0.0); - - for (int i = 0; i < num_samples; ++i) { - vec3 samplePoint = rayOrigin + rayDir * (primaryTime + stepSize * 0.5); - - float altitude = length(samplePoint) - earthRadius; - - float densityRayleigh = exp(-altitude / rayleigh_scale_height); - float densityMie = exp(-altitude / mie_scale_height); - float densityAbsortion = clamp((1.0 / cosh((absortion_scale_height - altitude) / absortion_falloff)) * densityRayleigh, 0.0, 1.0); - float stepOpticalDepthRayleigh = densityRayleigh * stepSize; - float stepOpticalDepthMie = densityMie * stepSize; - float stepOpticalDepthAbsortion = densityAbsortion * stepSize; - opticalDepthRayleigh += stepOpticalDepthRayleigh; - opticalDepthMie += stepOpticalDepthMie; - opticalDepthAbsortion += stepOpticalDepthAbsortion; - - vec2 pl = raySphereIntersection(samplePoint, lightDir, ATMOSPHERE_RADIUS); - float stepSizeLight = pl.y / float(num_light_samples); - - float opticalDepthLightRayleigh = 0.0; - float opticalDepthLightMie = 0.0; - float opticalDepthLightAbsortion = 0.0; - - float secondaryTime = 0.0; - - int j; - for (j = 0; j < num_light_samples; ++j) { - vec3 samplePointLight = samplePoint + lightDir * - (secondaryTime + stepSizeLight * 0.5); - - float altitudeLight = length(samplePointLight) - earthRadius; - if (altitudeLight < 0.0) - break; - - float densityLightRayleigh = exp(-altitudeLight / rayleigh_scale_height); - float densityLightMie = exp(-altitudeLight / mie_scale_height); - float densityLightAbsortion = clamp((1.0 / cosh((absortion_scale_height - altitudeLight) / absortion_falloff)) * densityLightRayleigh, 0.0, 1.0); - opticalDepthLightRayleigh += densityLightRayleigh * stepSizeLight; - opticalDepthLightMie += densityLightMie * stepSizeLight; - opticalDepthLightAbsortion += densityLightAbsortion * stepSizeLight; - - secondaryTime += stepSizeLight; - } - - if (j == num_light_samples) { - vec3 tau = - beta_rayleigh * (opticalDepthRayleigh + opticalDepthLightRayleigh) + - beta_mie * (opticalDepthMie + opticalDepthLightMie) + - beta_absortion * (opticalDepthAbsortion + opticalDepthLightAbsortion); - vec3 attenuation = exp(-tau); - - extinctionFactor += attenuation; - - totalRayleigh += stepOpticalDepthRayleigh * attenuation; - totalMie += stepOpticalDepthMie * attenuation; - } - - primaryTime += stepSize; - } - - transmittance = exp(-(beta_rayleigh * opticalDepthRayleigh + - beta_mie * opticalDepthMie + - beta_absortion * opticalDepthAbsortion)); - - inscatter = SUN_INTENSITY * - (totalRayleigh * beta_rayleigh * phaseRayleigh + - totalMie * beta_mie * phaseMie + - opticalDepthRayleigh * beta_ambient); -} diff --git a/Shaders/HDR/lighting.frag b/Shaders/HDR/lighting.frag index 186049e53..7365135be 100644 --- a/Shaders/HDR/lighting.frag +++ b/Shaders/HDR/lighting.frag @@ -12,14 +12,12 @@ uniform sampler2D ao_tex; uniform samplerCube prefiltered_envmap; uniform sampler2DShadow shadow_tex; uniform sampler2D dfg_lut; -uniform sampler2D aerial_inscatter_lut; -uniform sampler2D aerial_transmittance_lut; +uniform sampler2D aerial_perspective_lut; uniform mat4 fg_ViewMatrix; uniform mat4 fg_ViewMatrixInverse; uniform vec3 fg_SunDirection; uniform vec3 fg_CameraPositionCart; -uniform vec2 fg_NearFar; uniform mat4 fg_LightMatrix_csm0; uniform mat4 fg_LightMatrix_csm1; @@ -329,30 +327,34 @@ vec3 BRDF(in vec3 albedo, in float metalness, in float roughness, } //------------------------------------------------------------------------------ +// Atmospheric scattering -float map(float value, float min1, float max1, float min2, float max2) { - return min2 + (value - min1) * (max2 - min2) / (max1 - min1); -} - -vec3 sampleAerialPerspectiveSlice(sampler2D tex, int slice) +vec4 sampleAerialPerspectiveSlice(int slice) { + // Sample at the pixel center float offset = slice * AERIAL_LUT_TILE_SIZE + AERIAL_LUT_TEXEL_SIZE * 0.5; float x = texCoord.x * (AERIAL_LUT_TILE_SIZE - AERIAL_LUT_TEXEL_SIZE) + offset; - return texture(tex, vec2(x, texCoord.y)).rgb; + return texture(aerial_perspective_lut, vec2(x, texCoord.y)); } -vec3 sampleAerialPerspective(sampler2D tex, vec3 zero, float depth) +vec4 sampleAerialPerspective(float depth) { - vec3 color; - depth = min(abs(depth), AERIAL_MAX_DEPTH); - float d = map(depth, fg_NearFar.x, AERIAL_MAX_DEPTH, 0.0, AERIAL_SLICES); - if (d <= 1.0) { - color = mix(zero, sampleAerialPerspectiveSlice(tex, 0), d); + vec4 color; + // Map to [0,1] + float w = depth / AERIAL_MAX_DEPTH; + // Squared distribution + w = sqrt(clamp(w, 0.0, 1.0)); + // Remap to [0,16] to sample the right tile + w *= AERIAL_SLICES; + if (w <= 1.0) { + // Handle special case of fragments behind the first slice + color = mix(vec4(0.0, 0.0, 0.0, 1.0), sampleAerialPerspectiveSlice(0), w); } else { - d -= 1.0; - color = mix(sampleAerialPerspectiveSlice(tex, int(floor(d))), - sampleAerialPerspectiveSlice(tex, int(ceil(d))), - fract(d)); + w -= 1.0; // [0,15] + // Manually linearly interpolate between slices + color = mix(sampleAerialPerspectiveSlice(int(floor(w))), + sampleAerialPerspectiveSlice(int(ceil(w))), + fract(w)); } return color; } @@ -403,10 +405,8 @@ void main() float shadowFactor = getShadowing(pos, n, NdotL); vec3 color = ambient + brdf * sunIlluminance * shadowFactor; - vec3 inscatter = sampleAerialPerspective( - aerial_inscatter_lut, vec3(0.0), length(pos)); - vec3 transmittance = sampleAerialPerspective( - aerial_transmittance_lut, vec3(1.0), length(pos)); + vec4 aerialPerspective = sampleAerialPerspective(length(pos)); + color = color * aerialPerspective.a + aerialPerspective.rgb * SUN_INTENSITY; - fragHdrColor = color * transmittance + inscatter; + fragHdrColor = color; } diff --git a/Shaders/HDR/skydome.frag b/Shaders/HDR/skydome.frag index a55fda92b..113880021 100644 --- a/Shaders/HDR/skydome.frag +++ b/Shaders/HDR/skydome.frag @@ -2,19 +2,35 @@ out vec4 fragColor; -in vec3 rayDir; +in vec3 vRayDir; +in vec3 vRayDirView; 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); + +const float sun_solid_angle = 0.53*PI/180.0; // ~half a degree +const float sun_cos_solid_angle = cos(sun_solid_angle); void main() { + vec3 rayDir = normalize(vRayDir); float azimuth = atan(rayDir.y, rayDir.x) / PI * 0.5 + 0.5; // Undo the non-linear transformation from the sky-view LUT float l = asin(rayDir.z); 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; + + // Render the Sun disk + // XXX: Apply transmittance + float cosTheta = dot(normalize(vRayDirView), fg_SunDirection); + if (cosTheta >= sun_cos_solid_angle) + color += SUN_INTENSITY; + fragColor = vec4(color, 1.0); } diff --git a/Shaders/HDR/skydome.vert b/Shaders/HDR/skydome.vert index 1d77a203c..a141fd672 100644 --- a/Shaders/HDR/skydome.vert +++ b/Shaders/HDR/skydome.vert @@ -2,12 +2,15 @@ layout(location = 0) in vec4 pos; -out vec3 rayDir; +out vec3 vRayDir; +out vec3 vRayDirView; +uniform mat4 osg_ModelViewMatrix; uniform mat4 osg_ModelViewProjectionMatrix; void main() { gl_Position = osg_ModelViewProjectionMatrix * pos; - rayDir = normalize(pos.xyz); + vRayDir = normalize(pos.xyz); + vRayDirView = (osg_ModelViewMatrix * vec4(vRayDir, 0.0)).xyz; }