test: rework depth algorithm

data/shaders/opengl/geosphere_sky.frag

@@ -44,14 +44,96 @@ 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);
 	vec3 betaR = vec3(3.8e-6f, 13.5e-6f, 33.1e-6f);
 	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;