|
|
@@ -0,0 +1,405 @@
|
|
|
+// Copyright (C) 2009-2023, Panagiotis Christopoulos Charitos and contributors.
|
|
|
+// All rights reserved.
|
|
|
+// Code licensed under the BSD License.
|
|
|
+// http://www.anki3d.org/LICENSE
|
|
|
+
|
|
|
+// Shamelessly copied from https://www.shadertoy.com/view/slSXRW
|
|
|
+
|
|
|
+#include <AnKi/Shaders/Sky.hlsl>
|
|
|
+
|
|
|
+constexpr F32 kAtmosphereRadiusMM = 6.460f;
|
|
|
+
|
|
|
+// These are per megameter.
|
|
|
+constexpr Vec3 kRayleighScatteringBase = Vec3(5.802f, 13.558f, 33.1f);
|
|
|
+constexpr F32 kRayleighAbsorptionBase = 0.0f;
|
|
|
+
|
|
|
+constexpr F32 kMieScatteringBase = 3.996f;
|
|
|
+constexpr F32 kMieAbsorptionBase = 4.4f;
|
|
|
+
|
|
|
+constexpr Vec3 kOzoneAbsorptionBase = Vec3(0.650f, 1.881f, 0.085f);
|
|
|
+
|
|
|
+constexpr F32 kMiltipleScatteringSteps = 20.0f;
|
|
|
+constexpr F32 kSqrtSamples = 8.0f;
|
|
|
+
|
|
|
+constexpr F32 kSunTransmittanceSteps = 40.0f;
|
|
|
+
|
|
|
+constexpr F32 kScatteringSteps = 32.0f;
|
|
|
+
|
|
|
+constexpr Vec3 kGroundAlbedo = 0.3f;
|
|
|
+
|
|
|
+// From https://gamedev.stackexchange.com/questions/96459/fast-ray-sphere-collision-code.
|
|
|
+F32 rayIntersectSphere(Vec3 ro, Vec3 rd, F32 rad)
|
|
|
+{
|
|
|
+ const F32 b = dot(ro, rd);
|
|
|
+ const F32 c = dot(ro, ro) - rad * rad;
|
|
|
+
|
|
|
+ if(c > 0.0f && b > 0.0)
|
|
|
+ {
|
|
|
+ return -1.0;
|
|
|
+ }
|
|
|
+
|
|
|
+ const F32 discr = b * b - c;
|
|
|
+ if(discr < 0.0)
|
|
|
+ {
|
|
|
+ return -1.0;
|
|
|
+ }
|
|
|
+
|
|
|
+ // Special case: inside sphere, use far discriminant
|
|
|
+ if(discr > b * b)
|
|
|
+ {
|
|
|
+ return (-b + sqrt(discr));
|
|
|
+ }
|
|
|
+
|
|
|
+ return -b - sqrt(discr);
|
|
|
+}
|
|
|
+
|
|
|
+void getScatteringValues(Vec3 pos, out Vec3 rayleighScattering, out F32 mieScattering, out Vec3 extinction)
|
|
|
+{
|
|
|
+ const F32 altitudeKM = (length(pos) - kGroundRadiusMM) * 1000.0f;
|
|
|
+ // Note: Paper gets these switched up.
|
|
|
+ const F32 rayleighDensity = exp(-altitudeKM / 8.0f);
|
|
|
+ const F32 mieDensity = exp(-altitudeKM / 1.2f);
|
|
|
+
|
|
|
+ rayleighScattering = kRayleighScatteringBase * rayleighDensity;
|
|
|
+ const F32 rayleighAbsorption = kRayleighAbsorptionBase * rayleighDensity;
|
|
|
+
|
|
|
+ mieScattering = kMieScatteringBase * mieDensity;
|
|
|
+ F32 mieAbsorption = kMieAbsorptionBase * mieDensity;
|
|
|
+
|
|
|
+ const Vec3 ozoneAbsorption = kOzoneAbsorptionBase * max(0.0f, 1.0f - abs(altitudeKM - 25.0f) / 15.0f);
|
|
|
+
|
|
|
+ extinction = rayleighScattering + rayleighAbsorption + mieScattering + mieAbsorption + ozoneAbsorption;
|
|
|
+}
|
|
|
+
|
|
|
+/// Get value from transmittance LUT.
|
|
|
+Vec3 getValFromTLut(Texture2D<Vec4> tex, SamplerState linearAnyClampSampler, Vec3 pos, Vec3 dirToSun)
|
|
|
+{
|
|
|
+ const F32 height = length(pos);
|
|
|
+ const Vec3 up = pos / height;
|
|
|
+ const F32 sunCosZenithAngle = dot(dirToSun, up);
|
|
|
+ const Vec2 uv = Vec2(saturate(0.5f + 0.5f * sunCosZenithAngle), saturate((height - kGroundRadiusMM) / (kAtmosphereRadiusMM - kGroundRadiusMM)));
|
|
|
+ return tex.SampleLevel(linearAnyClampSampler, uv, 0.0f).xyz;
|
|
|
+}
|
|
|
+
|
|
|
+Vec3 getValFromMultiScattLut(Texture2D<Vec4> tex, SamplerState linearAnyClampSampler, Vec3 pos, Vec3 dirToSun)
|
|
|
+{
|
|
|
+ const F32 height = length(pos);
|
|
|
+ const Vec3 up = pos / height;
|
|
|
+ const F32 sunCosZenithAngle = dot(dirToSun, up);
|
|
|
+ const Vec2 uv = Vec2(saturate(0.5 + 0.5 * sunCosZenithAngle), saturate((height - kGroundRadiusMM) / (kAtmosphereRadiusMM - kGroundRadiusMM)));
|
|
|
+ return tex.SampleLevel(linearAnyClampSampler, uv, 0.0f).xyz;
|
|
|
+}
|
|
|
+
|
|
|
+F32 getSunAltitude(F32 hour24)
|
|
|
+{
|
|
|
+ const F32 periodSec = 24.0f;
|
|
|
+ const F32 halfPeriod = periodSec / 2.0f;
|
|
|
+ const F32 sunriseShift = 0.8f; // The closer it is to 1 the more time the sun will stay hidden
|
|
|
+ F32 cyclePoint = (1.0f - abs((fmod(hour24, periodSec) - halfPeriod) / halfPeriod));
|
|
|
+ cyclePoint = (cyclePoint * (1.0f + sunriseShift)) - sunriseShift;
|
|
|
+ return (0.5f * kPi) * cyclePoint;
|
|
|
+}
|
|
|
+
|
|
|
+Vec3 getSunDir(F32 time)
|
|
|
+{
|
|
|
+ const F32 altitude = getSunAltitude(time);
|
|
|
+ return normalize(Vec3(0.0f, sin(altitude), -cos(altitude)));
|
|
|
+}
|
|
|
+
|
|
|
+F32 getMiePhase(F32 cosTheta)
|
|
|
+{
|
|
|
+ const F32 g = 0.8f;
|
|
|
+ const F32 scale = 3.0f / (8.0f * kPi);
|
|
|
+
|
|
|
+ const F32 num = (1.0f - g * g) * (1.0f + cosTheta * cosTheta);
|
|
|
+ const F32 denom = (2.0f + g * g) * pow((1.0f + g * g - 2.0f * g * cosTheta), 1.5f);
|
|
|
+
|
|
|
+ return scale * num / denom;
|
|
|
+}
|
|
|
+
|
|
|
+F32 getRayleighPhase(F32 cosTheta)
|
|
|
+{
|
|
|
+ const F32 k = 3.0f / (16.0f * kPi);
|
|
|
+ return k * (1.0f + cosTheta * cosTheta);
|
|
|
+}
|
|
|
+
|
|
|
+// ===========================================================================
|
|
|
+// Transmitance LUT =
|
|
|
+// ===========================================================================
|
|
|
+#pragma anki technique_start comp SkyTransmittanceLut
|
|
|
+
|
|
|
+[[vk::binding(0)]] RWTexture2D<Vec4> g_tLutUav;
|
|
|
+
|
|
|
+Vec3 getSunTransmittance(Vec3 pos, Vec3 dirToSun)
|
|
|
+{
|
|
|
+ if(rayIntersectSphere(pos, dirToSun, kGroundRadiusMM) > 0.0f)
|
|
|
+ {
|
|
|
+ return 0.0f;
|
|
|
+ }
|
|
|
+
|
|
|
+ const F32 atmoDist = rayIntersectSphere(pos, dirToSun, kAtmosphereRadiusMM);
|
|
|
+ F32 t = 0.0f;
|
|
|
+
|
|
|
+ Vec3 transmittance = 1.0f;
|
|
|
+ for(F32 i = 0.0f; i < kSunTransmittanceSteps; i += 1.0f)
|
|
|
+ {
|
|
|
+ const F32 newT = ((i + 0.3f) / kSunTransmittanceSteps) * atmoDist;
|
|
|
+ const F32 dt = newT - t;
|
|
|
+ t = newT;
|
|
|
+
|
|
|
+ const Vec3 newPos = pos + t * dirToSun;
|
|
|
+
|
|
|
+ Vec3 rayleighScattering, extinction;
|
|
|
+ F32 mieScattering;
|
|
|
+ getScatteringValues(newPos, rayleighScattering, mieScattering, extinction);
|
|
|
+
|
|
|
+ transmittance *= exp(-dt * extinction);
|
|
|
+ }
|
|
|
+ return transmittance;
|
|
|
+}
|
|
|
+
|
|
|
+[numthreads(8, 8, 1)] void main(UVec2 svDispatchThreadId : SV_DISPATCHTHREADID)
|
|
|
+{
|
|
|
+ Vec2 lutSize;
|
|
|
+ g_tLutUav.GetDimensions(lutSize.x, lutSize.y);
|
|
|
+
|
|
|
+ const Vec2 uv = (Vec2(svDispatchThreadId) + 0.5f) / lutSize;
|
|
|
+
|
|
|
+ const F32 sunCosTheta = 2.0f * uv.x - 1.0f;
|
|
|
+ const F32 sunTheta = safeacos(sunCosTheta);
|
|
|
+ const F32 height = lerp(kGroundRadiusMM, kAtmosphereRadiusMM, uv.y);
|
|
|
+
|
|
|
+ const Vec3 pos = Vec3(0.0f, height, 0.0f);
|
|
|
+ const Vec3 dirToSun = normalize(Vec3(0.0f, sunCosTheta, -sin(sunTheta)));
|
|
|
+
|
|
|
+ g_tLutUav[svDispatchThreadId] = Vec4(getSunTransmittance(pos, dirToSun), 1.0f);
|
|
|
+}
|
|
|
+
|
|
|
+#pragma anki technique_end comp SkyTransmittanceLut
|
|
|
+
|
|
|
+// ===========================================================================
|
|
|
+// Multiple scattering LUT =
|
|
|
+// ===========================================================================
|
|
|
+#pragma anki technique_start comp SkyMultipleScatteringLut
|
|
|
+
|
|
|
+[[vk::binding(0)]] Texture2D<Vec4> g_tLutTex;
|
|
|
+[[vk::binding(1)]] SamplerState g_linearAnyClampSampler;
|
|
|
+[[vk::binding(2)]] RWTexture2D<Vec4> g_mLutUav;
|
|
|
+
|
|
|
+Vec3 getSphericalDir(F32 theta, F32 phi)
|
|
|
+{
|
|
|
+ const F32 cosPhi = cos(phi);
|
|
|
+ const F32 sinPhi = sin(phi);
|
|
|
+ const F32 cosTheta = cos(theta);
|
|
|
+ const F32 sinTheta = sin(theta);
|
|
|
+ return Vec3(sinPhi * sinTheta, cosPhi, sinPhi * cosTheta);
|
|
|
+}
|
|
|
+
|
|
|
+// Calculates Equation (5) and (7) from the paper.
|
|
|
+void getMulScattValues(Vec3 pos, Vec3 dirToSun, out Vec3 lumTotal, out Vec3 fms)
|
|
|
+{
|
|
|
+ lumTotal = 0.0f;
|
|
|
+ fms = 0.0f;
|
|
|
+
|
|
|
+ const F32 invSamples = 1.0f / F32(kSqrtSamples * kSqrtSamples);
|
|
|
+ for(F32 i = 0.0; i < kSqrtSamples; i += 1.0f)
|
|
|
+ {
|
|
|
+ for(F32 j = 0.0f; j < kSqrtSamples; j += 1.0f)
|
|
|
+ {
|
|
|
+ // This integral is symmetric about theta = 0 (or theta = PI), so we only need to integrate from zero to PI, not zero to 2*PI.
|
|
|
+ const F32 theta = kPi * (i + 0.5f) / kSqrtSamples;
|
|
|
+ const F32 phi = safeacos(1.0f - 2.0f * (j + 0.5f) / kSqrtSamples);
|
|
|
+ const Vec3 rayDir = getSphericalDir(theta, phi);
|
|
|
+
|
|
|
+ const F32 atmoDist = rayIntersectSphere(pos, rayDir, kAtmosphereRadiusMM);
|
|
|
+ const F32 groundDist = rayIntersectSphere(pos, rayDir, kGroundRadiusMM);
|
|
|
+ F32 tMax = atmoDist;
|
|
|
+ if(groundDist > 0.0f)
|
|
|
+ {
|
|
|
+ tMax = groundDist;
|
|
|
+ }
|
|
|
+
|
|
|
+ const F32 cosTheta = dot(rayDir, dirToSun);
|
|
|
+
|
|
|
+ const F32 miePhaseValue = getMiePhase(cosTheta);
|
|
|
+ const F32 rayleighPhaseValue = getRayleighPhase(-cosTheta);
|
|
|
+
|
|
|
+ Vec3 lum = 0.0f;
|
|
|
+ Vec3 lumFactor = 0.0f;
|
|
|
+ Vec3 transmittance = 1.0f;
|
|
|
+ F32 t = 0.0f;
|
|
|
+ for(F32 stepI = 0.0f; stepI < kMiltipleScatteringSteps; stepI += 1.0f)
|
|
|
+ {
|
|
|
+ const F32 newT = ((stepI + 0.3f) / kMiltipleScatteringSteps) * tMax;
|
|
|
+ const F32 dt = newT - t;
|
|
|
+ t = newT;
|
|
|
+
|
|
|
+ const Vec3 newPos = pos + t * rayDir;
|
|
|
+
|
|
|
+ Vec3 rayleighScattering, extinction;
|
|
|
+ F32 mieScattering;
|
|
|
+ getScatteringValues(newPos, rayleighScattering, mieScattering, extinction);
|
|
|
+
|
|
|
+ const Vec3 sampleTransmittance = exp(-dt * extinction);
|
|
|
+
|
|
|
+ // Integrate within each segment.
|
|
|
+ const Vec3 scatteringNoPhase = rayleighScattering + mieScattering;
|
|
|
+ const Vec3 scatteringF = (scatteringNoPhase - scatteringNoPhase * sampleTransmittance) / extinction;
|
|
|
+ lumFactor += transmittance * scatteringF;
|
|
|
+
|
|
|
+ // This is slightly different from the paper, but I think the paper has a mistake? In equation (6), I think S(x,w_s) should be
|
|
|
+ // S(x-tv,w_s).
|
|
|
+ const Vec3 sunTransmittance = getValFromTLut(g_tLutTex, g_linearAnyClampSampler, newPos, dirToSun);
|
|
|
+
|
|
|
+ const Vec3 rayleighInScattering = rayleighScattering * rayleighPhaseValue;
|
|
|
+ const F32 mieInScattering = mieScattering * miePhaseValue;
|
|
|
+ const Vec3 inScattering = (rayleighInScattering + mieInScattering) * sunTransmittance;
|
|
|
+
|
|
|
+ // Integrated scattering within path segment.
|
|
|
+ const Vec3 scatteringIntegral = (inScattering - inScattering * sampleTransmittance) / extinction;
|
|
|
+
|
|
|
+ lum += scatteringIntegral * transmittance;
|
|
|
+ transmittance *= sampleTransmittance;
|
|
|
+ }
|
|
|
+
|
|
|
+ if(groundDist > 0.0f)
|
|
|
+ {
|
|
|
+ Vec3 hitPos = pos + groundDist * rayDir;
|
|
|
+ if(dot(pos, dirToSun) > 0.0f)
|
|
|
+ {
|
|
|
+ hitPos = normalize(hitPos) * kGroundRadiusMM;
|
|
|
+ lum += transmittance * kGroundAlbedo * getValFromTLut(g_tLutTex, g_linearAnyClampSampler, hitPos, dirToSun);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ fms += lumFactor * invSamples;
|
|
|
+ lumTotal += lum * invSamples;
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+[numthreads(8, 8, 1)] void main(UVec2 svDispatchThreadId : SV_DISPATCHTHREADID)
|
|
|
+{
|
|
|
+ Vec2 lutSize;
|
|
|
+ g_mLutUav.GetDimensions(lutSize.x, lutSize.y);
|
|
|
+
|
|
|
+ const Vec2 uv = (Vec2(svDispatchThreadId) + 0.5f) / lutSize;
|
|
|
+
|
|
|
+ const F32 sunCosTheta = 2.0f * uv.x - 1.0f;
|
|
|
+ const F32 sunTheta = safeacos(sunCosTheta);
|
|
|
+ const F32 height = lerp(kGroundRadiusMM, kAtmosphereRadiusMM, uv.y);
|
|
|
+
|
|
|
+ const Vec3 pos = Vec3(0.0f, height, 0.0f);
|
|
|
+ const Vec3 dirToSun = normalize(Vec3(0.0f, sunCosTheta, -sin(sunTheta)));
|
|
|
+
|
|
|
+ Vec3 lum, f_ms;
|
|
|
+ getMulScattValues(pos, dirToSun, lum, f_ms);
|
|
|
+
|
|
|
+ // Equation 10 from the paper.
|
|
|
+ const Vec3 psi = lum / (1.0f - f_ms);
|
|
|
+
|
|
|
+ g_mLutUav[svDispatchThreadId] = Vec4(psi, 0.0f);
|
|
|
+}
|
|
|
+
|
|
|
+#pragma anki technique_end comp SkyMultipleScatteringLut
|
|
|
+
|
|
|
+// ===========================================================================
|
|
|
+// Sky LUT =
|
|
|
+// ===========================================================================
|
|
|
+#pragma anki technique_start comp SkyLut
|
|
|
+
|
|
|
+[[vk::binding(0)]] Texture2D<Vec4> g_tLutTex;
|
|
|
+[[vk::binding(1)]] Texture2D<Vec4> g_mLutTex;
|
|
|
+[[vk::binding(2)]] SamplerState g_linearAnyClampSampler;
|
|
|
+[[vk::binding(3)]] RWTexture2D<Vec4> g_skyLutUav;
|
|
|
+
|
|
|
+struct Consts
|
|
|
+{
|
|
|
+ Vec3 m_dirLightDirection;
|
|
|
+ F32 m_padding;
|
|
|
+};
|
|
|
+
|
|
|
+[[vk::push_constant]] ConstantBuffer<Consts> g_consts;
|
|
|
+
|
|
|
+Vec3 raymarchScattering(Vec3 pos, Vec3 rayDir, Vec3 dirToSun, F32 tMax, F32 numSteps)
|
|
|
+{
|
|
|
+ const F32 cosTheta = dot(rayDir, dirToSun);
|
|
|
+
|
|
|
+ const F32 miePhaseValue = getMiePhase(cosTheta);
|
|
|
+ const F32 rayleighPhaseValue = getRayleighPhase(-cosTheta);
|
|
|
+
|
|
|
+ Vec3 lum = 0.0f;
|
|
|
+ Vec3 transmittance = 1.0f;
|
|
|
+ F32 t = 0.0f;
|
|
|
+ for(F32 i = 0.0; i < numSteps; i += 1.0)
|
|
|
+ {
|
|
|
+ const F32 newT = ((i + 0.3f) / numSteps) * tMax;
|
|
|
+ const F32 dt = newT - t;
|
|
|
+ t = newT;
|
|
|
+
|
|
|
+ const Vec3 newPos = pos + t * rayDir;
|
|
|
+
|
|
|
+ Vec3 rayleighScattering, extinction;
|
|
|
+ F32 mieScattering;
|
|
|
+ getScatteringValues(newPos, rayleighScattering, mieScattering, extinction);
|
|
|
+
|
|
|
+ const Vec3 sampleTransmittance = exp(-dt * extinction);
|
|
|
+
|
|
|
+ const Vec3 sunTransmittance = getValFromTLut(g_tLutTex, g_linearAnyClampSampler, newPos, dirToSun);
|
|
|
+ const Vec3 psiMS = getValFromMultiScattLut(g_mLutTex, g_linearAnyClampSampler, newPos, dirToSun);
|
|
|
+
|
|
|
+ const Vec3 rayleighInScattering = rayleighScattering * (rayleighPhaseValue * sunTransmittance + psiMS);
|
|
|
+ const Vec3 mieInScattering = mieScattering * (miePhaseValue * sunTransmittance + psiMS);
|
|
|
+ const Vec3 inScattering = (rayleighInScattering + mieInScattering);
|
|
|
+
|
|
|
+ // Integrated scattering within path segment.
|
|
|
+ const Vec3 scatteringIntegral = (inScattering - inScattering * sampleTransmittance) / extinction;
|
|
|
+
|
|
|
+ lum += scatteringIntegral * transmittance;
|
|
|
+ transmittance *= sampleTransmittance;
|
|
|
+ }
|
|
|
+
|
|
|
+ return lum;
|
|
|
+}
|
|
|
+
|
|
|
+[numthreads(8, 8, 1)] void main(UVec2 svDispatchThreadId : SV_DISPATCHTHREADID)
|
|
|
+{
|
|
|
+ Vec2 lutSize;
|
|
|
+ g_skyLutUav.GetDimensions(lutSize.x, lutSize.y);
|
|
|
+
|
|
|
+ const Vec2 uv = (Vec2(svDispatchThreadId) + 0.5f) / lutSize;
|
|
|
+
|
|
|
+ const F32 azimuthAngle = (uv.x - 0.5f) * 2.0f * kPi;
|
|
|
+ // Non-linear mapping of altitude. See Section 5.3 of the paper.
|
|
|
+ F32 adjV;
|
|
|
+ if(uv.y < 0.5f)
|
|
|
+ {
|
|
|
+ const F32 coord = 1.0f - 2.0f * uv.y;
|
|
|
+ adjV = -coord * coord;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ {
|
|
|
+ const F32 coord = uv.y * 2.0f - 1.0f;
|
|
|
+ adjV = coord * coord;
|
|
|
+ }
|
|
|
+
|
|
|
+ const F32 height = length(kViewPos);
|
|
|
+ const Vec3 up = kViewPos / height;
|
|
|
+ const F32 horizonAngle = safeacos(sqrt(height * height - kGroundRadiusMM * kGroundRadiusMM) / height) - 0.5f * kPi;
|
|
|
+ const F32 altitudeAngle = adjV * 0.5f * kPi - horizonAngle;
|
|
|
+
|
|
|
+ const F32 cosAltitude = cos(altitudeAngle);
|
|
|
+ const Vec3 rayDir = Vec3(cosAltitude * sin(azimuthAngle), sin(altitudeAngle), -cosAltitude * cos(azimuthAngle));
|
|
|
+
|
|
|
+ const F32 sunAltitude = (0.5f * kPi) - acos(dot(g_consts.m_dirLightDirection, up));
|
|
|
+ const Vec3 dirToSun = Vec3(0.0f, sin(sunAltitude), -cos(sunAltitude));
|
|
|
+
|
|
|
+ const F32 atmoDist = rayIntersectSphere(kViewPos, rayDir, kAtmosphereRadiusMM);
|
|
|
+ const F32 groundDist = rayIntersectSphere(kViewPos, rayDir, kGroundRadiusMM);
|
|
|
+ const F32 tMax = (groundDist < 0.0f) ? atmoDist : groundDist;
|
|
|
+ const Vec3 lum = raymarchScattering(kViewPos, rayDir, dirToSun, tMax, kScatteringSteps);
|
|
|
+
|
|
|
+ g_skyLutUav[svDispatchThreadId] = Vec4(lum, 0.0);
|
|
|
+}
|
|
|
+
|
|
|
+#pragma anki technique_end comp SkyLut
|
Panagiotis Christopoulos Charitos