From dfada9c5efe5954e747996cfcbb157f819523c77 Mon Sep 17 00:00:00 2001 From: Oleg Solovyov Date: Fri, 1 Sep 2023 19:21:21 +0300 Subject: [PATCH] test: rework depth algorithm --- data/shaders/opengl/geosphere_sky.frag | 135 ++++++++++++++++++++++++- 1 file changed, 134 insertions(+), 1 deletion(-) diff --git a/data/shaders/opengl/geosphere_sky.frag b/data/shaders/opengl/geosphere_sky.frag index ea75ec61f2..5cb5fb350f 100644 --- a/data/shaders/opengl/geosphere_sky.frag +++ b/data/shaders/opengl/geosphere_sky.frag @@ -44,6 +44,87 @@ void scatter(out float hr, out float hm, const in vec3 orig, const in vec3 cente hm = -height / Hm; } +void findClosestHeight(out float h, out float t, const in vec3 orig, const in vec3 dir, const in vec3 center) +{ + vec3 radiusVector = orig - center; + float tangent = dot(dir, radiusVector); + vec3 normal = radiusVector - (tangent * dir); + h = dot(normal, normal); + t = tangent; +} + +// compute exact density, used for obtaining coefficients +float computeDensity(const in float radius, const in float atmosphereHeight, const in float h, const in float baricStep) +{ + int numSamples = 16; + float totalHeight = radius + atmosphereHeight; + float minHeight = radius + h; + float tmax = sqrt(totalHeight*totalHeight - minHeight*minHeight); // maximum t + float tmin = -tmax; + float length = tmax - tmin; + float step = length / numSamples; + + float density = 0.0f; + for (int i = 0; i < numSamples; ++i) { + float t = tmin + step * (i + 0.5f); + float h = sqrt(minHeight*minHeight + t*t); + density += step * exp(-h / baricStep); + } + + return density; +} + +// error function +float erf(const in float x) +{ + float a = 0.14001228904; + float four_over_pi = 1.27323954; + + float x2 = x*x; + float r = -x2 * (four_over_pi + a * x2) / (1 + a * x2); + if (x > 0) + return sqrt(1-exp(r)) * 0.5 + 0.5; + else + return -sqrt(1-exp(r)) * 0.5 + 0.5; +} + +// predict out-scattering density based on coefficients +float predictDensityOut(const in float atmosphereHeight, const in float height, const in float k, const in float b) +{ + if (height > atmosphereHeight) + return 0.f; + + return k * exp(-height * b); +} + +// predict in-scattering rate: from 0 to 1 +float predictDensityIn(const in float radius, const in float atmosphereHeight, const in float height, const in float c, const in float t) +{ + float minHeight = radius + height; + float h = sqrt(minHeight*minHeight + t*t); // height for our position + + if (h > radius + atmosphereHeight) { + if (t > 0) + return 1.f; // no in-scattering, looking towards atmosphere + else + return 0.f; // looking from atmosphere + } else { + return erf(c * t); // erf is monotonic ascending + } +} + +void getCoefficients(out float k, out float b, out float c, const in float radius, const in float atmHeight, const in float baricStep) +{ + k = computeDensity(radius, atmHeight, 0.f, baricStep); + b = log(k / computeDensity(radius, atmHeight, 1.f, baricStep)); + + float s = 0.001f; + float sHeight = sqrt(radius*radius + s*s) - radius; + float c1 = exp(-sHeight / baricStep) * s; + float c2 = k * (erf(s) - erf(0)); + c = c1 / c2; +} + vec3 computeIncidentLight(const in vec3 sunDirection, in vec3 orig, const in vec3 dir, const in vec3 center, in float tmin, in float tmax) { vec3 nSunDirection = normalize(sunDirection); @@ -51,7 +132,8 @@ vec3 computeIncidentLight(const in vec3 sunDirection, in vec3 orig, const in vec vec3 betaM = vec3(21e-6f); float earthRadius = geosphereRadius, - atmosphereRadius = geosphereRadius * geosphereAtmosTopRad; + atmosphereRadius = earthRadius * geosphereAtmosTopRad, + atmosphereHeight = atmosphereRadius - earthRadius; float t0 = 0.f; float t1 = 0.f; @@ -71,6 +153,56 @@ vec3 computeIncidentLight(const in vec3 sunDirection, in vec3 orig, const in vec float g = 0.76f; float phaseM = 3.f / (8.f * 3.141592) * ((1.f - g * g) * (1.f + mu * mu)) / ((2.f + g * g) * pow(1.f + g * g - 2.f * g * mu, 1.5f)); + float kR, bR, cR; + float kM, bM, cM; + float Hr = 7994; + float Hm = 1200; + getCoefficients(kR, bR, cR, earthRadius, atmosphereHeight, Hr); + getCoefficients(kM, bM, cM, earthRadius, atmosphereHeight, Hm); + + float h = 0.f; + float t = 0.f; + findClosestHeight(h, t, orig, dir, center); + h -= earthRadius; + + // find out-scattering density + opticalDepthR = predictDensityOut(atmosphereHeight, h, kR, bR); + opticalDepthM = predictDensityOut(atmosphereHeight, h, kM, bM); + + // apply in-scattering filter + opticalDepthR *= predictDensityIn(earthRadius, atmosphereHeight, h, cR, t); + opticalDepthM *= predictDensityIn(earthRadius, atmosphereHeight, h, cM, t); + + for (int i = 0; i < numSamples; ++i) { + vec3 samplePosition = orig + vec3(tCurrent + segmentLength * 0.5f) * dir; + float opticalDepthLightR = 0, opticalDepthLightM = 0; + + float hr, hm; + scatter(hr, hm, samplePosition, center); + + float hl = 0.f; + float tl = 0.f; + findClosestHeight(hl, tl, samplePosition, sunDirection, center); + hl -= earthRadius; + + // find out-scattering density + opticalDepthLightR = predictDensityOut(atmosphereHeight, hl, kR, bR); + opticalDepthLightM = predictDensityOut(atmosphereHeight, hl, kM, bM); + + // apply in-scattering filter + opticalDepthLightR *= predictDensityIn(earthRadius, atmosphereHeight, hl, cR, tl); + opticalDepthLightM *= predictDensityIn(earthRadius, atmosphereHeight, hl, cM, tl); + + vec3 tau = -(betaR * (opticalDepthR + opticalDepthLightR) + betaM * 1.1f * (opticalDepthM + opticalDepthLightM)); + vec3 tauR = tau + vec3(hr); + vec3 tauM = tau + vec3(hm); + vec3 attenuationR = vec3(exp(tauR.x), exp(tauR.y), exp(tauR.z)) * segmentLength; + vec3 attenuationM = vec3(exp(tauM.x), exp(tauM.y), exp(tauM.z)) * segmentLength; + sumR += attenuationR; + sumM += attenuationM; + } + + /* for (int i = 0; i < numSamples; ++i) { vec3 samplePosition = orig + vec3(tCurrent + segmentLength * 0.5f) * dir; float hr, hm; @@ -102,6 +234,7 @@ vec3 computeIncidentLight(const in vec3 sunDirection, in vec3 orig, const in vec } tCurrent += segmentLength; } + */ vec3 ret = (sumR * betaR * phaseR + sumM * betaM * phaseM); return ret;