Rayleigh v2

data/shaders/opengl/basesphere_uniforms.glsl

@@ -26,7 +26,7 @@ layout(std140) uniform BaseSphereData {
 	vec4 atmosColor;
 	vec3 coefficientsR;			// coefficients for approximating the Rayleigh contribution
 	vec3 coefficientsM;			// coefficients for approximating the Mie contribution
-	vec2 scaleHeight;			// height for (R, M) in km, at which density will be reduced by e
+	vec2 scaleHeight;			// height for (R, M) in m, at which density will be reduced by e
 
 	// Eclipse data
 	Eclipse eclipse;

data/shaders/opengl/geosphere_terrain.frag

@@ -4,12 +4,14 @@
 #include "attributes.glsl"
 #include "lib.glsl"
 #include "basesphere_uniforms.glsl"
+#include "rayleigh.glsl"
 
 #define WATER_SHINE 16.0
 
 uniform sampler2D texture0;
 uniform sampler2D texture1;
 in vec2 texCoord0;
+uniform sampler2D scatterLUT;
 
 uniform int NumShadows;
 
@@ -40,6 +42,9 @@ void main(void)
 	vec3 tnorm = normalize(varyingNormal);
 	vec4 diff = vec4(0.0);
 
+	vec2 planetDist = raySphereIntersect(geosphereCenter, eyenorm, geosphereRadius);
+	vec3 planetIntersect = geosphereCenter + (eyenorm * planetDist.x);
+
 	float surfaceDist = dist * geosphereInvRadius;
 
 	// calculate the detail texture contribution from hi and lo textures
@@ -49,11 +54,19 @@ void main(void)
 	vec4 detailMul = mix(vec4(1.0), detailVal, detailMix);
 
 #ifdef TERRAIN_WITH_WATER
-	float specularReflection=0.0;
+	vec4 waterSpecular = vec4(0.0);
 #endif
 
 #if (NUM_LIGHTS > 0)
 	vec3 V = normalize(eyeposScaled - geosphereCenter);
+	vec3 I = normalize(eyeposScaled - planetIntersect);
+
+	float AU = 149598000000.0;
+	frag_color = vec4(0.f);
+
+	// coordinates, in planet radius
+	vec4 planet = vec4(geosphereCenter, geosphereRadius);
+	vec4 atmosphere = vec4(geosphereCenter, geosphereAtmosTopRad);
 
 	for (int i=0; i<NUM_LIGHTS; ++i) {
 		vec3 L = normalize(uLight[i].position.xyz);
@@ -63,7 +76,11 @@ void main(void)
 #ifdef TERRAIN_WITH_WATER
 		//water only for specular
 	    if (vertexColor.b > 0.05 && vertexColor.r < 0.05) {
-			CalcPlanetSpec(specularReflection, uLight[i], L, tnorm, eyenorm, WATER_SHINE);
+			vec3 lightPosAU = uLight[i].position.xyz / AU;
+			float intensity = 1.f / dot(lightPosAU, lightPosAU); // magic to avoid calculating length and then squaring it
+
+			waterSpecular.xyz += calculateAtmosphereColor(planet, atmosphere, toLinear(uLight[i].diffuse), reflect(L, I), (I - geosphereCenter) * geosphereRadius, eyenorm, uneclipsed, scatterLUT) * intensity;
+			//waterSpecular.xyz += computeIncidentLight(reflect(L, I), eyenorm, I * geosphereRadius, toLinear(uLight[i].diffuse), uneclipsed) * intensity;
 		}
 #endif
 	}
@@ -72,16 +89,12 @@ void main(void)
 	vec4 ambient = scene.ambient * vertexColor;
 	vec4 final = vertexColor * detailMul * diff;
 
-#ifdef ATMOSPHERE
-	float ldprod=0.0;
-	float fogFactor=0.0;
-	CalcPlanetFogFactor(ldprod, fogFactor, eyenorm, eyeposScaled - geosphereCenter, surfaceDist);
-
-	//calculate sunset tone red when passing through more atmosphere, clamp everything.
-	float atmpower = (diff.r+diff.g+diff.b)/3.0;
-	vec4 sunset = vec4(0.8,clamp(pow(atmpower,0.8),0.0,1.0),clamp(pow(atmpower,1.2),0.0,1.0),1.0);
+#ifdef TERRAIN_WITH_WATER
+	vec4 waterColor = waterSpecular * 20.f;
+#endif
 
-	final = fogFactor * (ambient + final);
+#ifdef ATMOSPHERE
+	final = ambient + final;
 #else
 	// add extra brightness to atmosphere-less planetoids and dim stars
 	final = ambient + final * 2.0;
@@ -95,13 +108,10 @@ void main(void)
 		final;
 
 #ifdef ATMOSPHERE
-	frag_color +=
-		(1.0-fogFactor) * (diff*atmosColor) +
+	frag_color += diff * atmosColor;
 #ifdef TERRAIN_WITH_WATER
-		  diff * specularReflection * sunset +
+	frag_color += diff * waterColor;
 #endif
-		  (0.02-clamp(fogFactor,0.0,0.01))*diff*ldprod*sunset +	      //increase fog scatter
-		  (pow((1.0-pow(fogFactor,0.75)),256.0)*0.4*diff*atmosColor)*sunset;  //distant fog.
 #endif //ATMOSPHERE
 
 #else // NUM_LIGHTS > 0 -- unlit rendering - stars

data/shaders/opengl/geosphere_terrain.shaderdef

@@ -5,6 +5,7 @@ Shader geosphere_terrain
 
 Texture sampler2D texture0
 Texture sampler2D texture1
+Texture sampler2D scatterLUT
 
 Buffer LightData binding=0
 Buffer DrawData binding=1

data/shaders/opengl/rayleigh.glsl

@@ -1,3 +1,4 @@
+// get height in m
 float height(const in vec3 orig, const in vec3 center)
 {
 	vec3 r = orig - center;
@@ -6,16 +7,44 @@ float height(const in vec3 orig, const in vec3 center)
 	return height;
 }
 
-void scatter(out vec2 density, const in vec3 orig, const in vec3 center)
+// get density at given point
+vec2 getDensityAtPointOld(const in vec3 orig, const in vec3 center, const in sampler2D texture_LUT)
 {
-	float height = height(orig, center);
+	float height = height(orig, center); // in meters
 
-	density = -height / scaleHeight;
+	vec2 density = -height / scaleHeight;
 
 	// earth atmospheric density: 1.225 kg/m^3, divided by 1e5
 	// 1/1.225e-5 = 81632.65306
-	float earthDensities = geosphereAtmosFogDensity * 81632.65306f;
-	density /= earthDensities;
+	float earthDensities = log(geosphereAtmosFogDensity * 81632.65306f);
+	density += earthDensities;
+	return exp(density);
+}
+
+// get density at given point
+vec2 getDensityAtPointNew(const in vec3 orig, const in vec3 center, const in sampler2D texture_LUT)
+{
+	float height = height(orig, center); // in meters
+	float maxHeight = (geosphereAtmosTopRad - 1.f) * geosphereRadius;
+
+	float ratio = height / maxHeight;
+
+	vec2 density = vec2(0.f);
+	// normalized to [0; 1) - remap back to [-128; 128)
+	density.x = texture(texture_LUT, vec2(0.0, ratio)).x;
+	density.y = texture(texture_LUT, vec2(1.0, ratio)).x;
+
+	density = density * 256 - 128;
+	return exp(density);
+}
+
+vec2 getDensityAtPoint(const in vec3 orig, const in vec3 center, const in sampler2D texture_LUT)
+{
+#if 1
+    return getDensityAtPointNew(orig, center, texture_LUT);
+#else
+    return getDensityAtPointOld(orig, center, texture_LUT);
+#endif
 }
 
 // orig: ray origin
@@ -112,7 +141,7 @@ float predictDensityInOut(const in vec3 sample, const in vec3 dir, const in vec3
     return opticalDepth;
 }
 
-void skipRay(inout vec2 opticalDepth, const in vec3 dir, const in vec2 boundaries, const in vec3 center)
+void skipRay(inout vec2 opticalDepth, const in vec3 dir, const in vec2 boundaries, const in vec3 center, const in sampler2D texture_LUT)
 {
 	if (boundaries.y == boundaries.x)
 		return;
@@ -125,15 +154,14 @@ void skipRay(inout vec2 opticalDepth, const in vec3 dir, const in vec2 boundarie
 		vec3 samplePosition = vec3(tCurrent + segmentLength * 0.5f) * dir;
 
 		// primary ray is approximated by (density * isegmentLength)
-		vec2 density;
-		scatter(density, samplePosition, center);
-		opticalDepth += exp(density) * segmentLength;
+		vec2 density = getDensityAtPoint(samplePosition, center, texture_LUT);
+		opticalDepth += density * segmentLength;
 
 		tCurrent += segmentLength;
 	}
 }
 
-void processRay(inout vec3 sumR, inout vec3 sumM, inout vec2 opticalDepth, const in vec3 sunDirection, const in vec3 dir, const in vec2 boundaries, const in vec3 center, const in vec4 diffuse, const in float uneclipsed)
+void processRayFast(inout vec3 sumR, inout vec3 sumM, inout vec2 opticalDepth, const in vec3 sunDirection, const in vec3 dir, const in vec2 boundaries, const in vec3 center, const in vec4 diffuse, const in float uneclipsed, const in sampler2D texture_LUT)
 {
 	if (boundaries.y == boundaries.x)
 		return;
@@ -151,9 +179,9 @@ void processRay(inout vec3 sumR, inout vec3 sumM, inout vec2 opticalDepth, const
 	for (int i = 0; i < numSamples; ++i) {
 		vec3 samplePosition = vec3(tCurrent + segmentLength * 0.5f) * dir;
 
-		vec2 density;
-		scatter(density, samplePosition, center);
-		opticalDepth += exp(density) * segmentLength;
+		vec2 density = getDensityAtPoint(samplePosition, center, texture_LUT);
+		vec2 depthAdd = density * segmentLength;
+		opticalDepth += depthAdd;
 
 		// light optical depth
 		vec2 opticalDepthLight = vec2(0.f);
@@ -167,17 +195,71 @@ void processRay(inout vec3 sumR, inout vec3 sumM, inout vec2 opticalDepth, const
 		vec4 atmosDiffuse = vec4(0.f);
 		CalcPlanetDiffuse(atmosDiffuse, diffuse, sunDirection, surfaceNorm, uneclipsed);
 
-		vec3 tau = -(betaR * (opticalDepth.x + opticalDepthLight.x) + betaM * 1.1f * (opticalDepth.y + opticalDepthLight.y));
-		vec3 tauR = tau + vec3(density.x);
-		vec3 tauM = tau + vec3(density.y);
-		vec3 attenuationR = exp(tauR) * segmentLength;
-		vec3 attenuationM = exp(tauM) * segmentLength;
-		sumR += attenuationR * atmosDiffuse.xyz;
-		sumM += attenuationM * atmosDiffuse.xyz;
+		vec3 tau = exp(-(betaR * (opticalDepth.x + opticalDepthLight.x) + betaM * 1.1f * (opticalDepth.y + opticalDepthLight.y)));
+		sumR += tau * depthAdd.x * atmosDiffuse.xyz;
+		sumM += tau * depthAdd.y * atmosDiffuse.xyz;
 		tCurrent += segmentLength;
 	}
 }
 
+void processRayFull(inout vec3 sumR, inout vec3 sumM, inout vec2 opticalDepth, const in vec3 sunDirection, const in vec3 dir, const in vec2 boundaries, const in vec3 center, const in vec4 diffuse, const in float uneclipsed, const in sampler2D texture_LUT)
+{
+	if (boundaries.y == boundaries.x)
+		return;
+
+	vec3 betaR = vec3(3.8e-6f, 13.5e-6f, 33.1e-6f);
+	vec3 betaM = vec3(21e-6f);
+
+	int numSamples = 16;
+
+	float atmosphereRadius = geosphereRadius * geosphereAtmosTopRad,
+	      atmosphereHeight = atmosphereRadius - geosphereRadius;
+
+	float tCurrent = boundaries.x;
+	float segmentLength = (boundaries.y - boundaries.x) / numSamples;
+	for (int i = 0; i < numSamples; ++i) {
+		vec3 samplePosition = vec3(tCurrent + segmentLength * 0.5f) * dir;
+
+		vec2 density = getDensityAtPoint(samplePosition, center, texture_LUT);
+		vec2 depthAdd = density * segmentLength;
+		opticalDepth += depthAdd;
+
+		// light optical depth
+		vec2 opticalDepthLight = vec2(0.f);
+		vec3 sampleGeoCenter = center - samplePosition;
+
+		int numSamplesLight = 8;
+		vec2 boundariesLight = raySphereIntersect(sampleGeoCenter, sunDirection, geosphereRadius * geosphereAtmosTopRad);
+		float segmentLengthLight = boundariesLight.y / numSamplesLight;
+		float tCurrentLight = 0.f;
+		for (int j = 0; j < numSamplesLight; ++j) {
+			vec3 samplePositionLight = vec3(segmentLengthLight * 0.5f + tCurrentLight) * sunDirection + samplePosition;
+			vec2 densityLDir = getDensityAtPoint(samplePositionLight, center, texture_LUT);
+			opticalDepthLight += densityLDir * segmentLengthLight;
+
+			tCurrentLight += segmentLengthLight;
+		}
+
+		vec3 surfaceNorm = -normalize(sampleGeoCenter);
+		vec4 atmosDiffuse = vec4(0.f);
+		CalcPlanetDiffuse(atmosDiffuse, diffuse, sunDirection, surfaceNorm, uneclipsed);
+
+		vec3 tau = exp(-(betaR * (opticalDepth.x + opticalDepthLight.x) + betaM * 1.1f * (opticalDepth.y + opticalDepthLight.y)));
+		sumR += tau * depthAdd.x * atmosDiffuse.xyz;
+		sumM += tau * depthAdd.y * atmosDiffuse.xyz;
+		tCurrent += segmentLength;
+	}
+}
+
+void processRay(inout vec3 sumR, inout vec3 sumM, inout vec2 opticalDepth, const in vec3 sunDirection, const in vec3 dir, const in vec2 boundaries, const in vec3 center, const in vec4 diffuse, const in float uneclipsed, const in sampler2D texture_LUT)
+{
+#if 1
+    processRayFull(sumR, sumM, opticalDepth, sunDirection, dir, boundaries, center, diffuse, uneclipsed, texture_LUT);
+#else
+    processRayFast(sumR, sumM, opticalDepth, sunDirection, dir, boundaries, center, diffuse, uneclipsed, texture_LUT);
+#endif
+}
+
 // replace (a, b) by (b, a) if a > b
 vec2 sortAscending(const in vec2 segment)
 {
@@ -200,14 +282,19 @@ vec4 segmentSubtraction(const in vec2 a, const in vec2 b)
 	return c;
 }
 
-vec3 computeIncidentLight(const in vec3 sunDirection, const in vec3 dir, const in vec3 center, const in vec2 atmosDist, const in vec4 diffuse, const in float uneclipsed)
+/*
+ * given:
+ * sunDirection - direction from camera to light source, normalized
+ * dir - direction from camera to pixel being rendered, normalized
+ * center - position of planet relative to camera, in absolute scale
+ *
+ * splits given ray into parts by planet shadow:
+ * - xy: ray before shadow
+ * - yz: ray inside shadow
+ * - zw: ray after shadow
+ */
+vec4 getRaySegment(const in vec3 sunDirection, const in vec3 dir, const in vec3 center)
 {
-	vec3 betaR = vec3(3.8e-6f, 13.5e-6f, 33.1e-6f);
-	vec3 betaM = vec3(21e-6f);
-
-	float atmosMin = atmosDist.x * geosphereRadius;
-	float atmosMax = atmosDist.y * geosphereRadius;
-
 	// solve Cylinder entry/exit dist
 	vec2 cylinder_intersect = rayCylinderIntersect(dir, center, sunDirection, geosphereRadius);
 	bool hasIntersect = cylinder_intersect.x != 0 || cylinder_intersect.y != 0;
@@ -243,13 +330,23 @@ vec3 computeIncidentLight(const in vec3 sunDirection, const in vec3 dir, const i
 		atmosphere_minus_shadow.xyz = max(atmosphere_minus_shadow.xyz, ground_intersect.y);
 	}
 
+	return atmosphere_minus_shadow;
+}
+
+vec3 computeIncidentLight(const in vec3 sunDirection, const in vec3 dir, const in vec3 center, const in vec4 diffuse, const in float uneclipsed, const in sampler2D texture_LUT)
+{
+	vec3 betaR = vec3(3.8e-6f, 13.5e-6f, 33.1e-6f);
+	vec3 betaM = vec3(21e-6f);
+
 	vec3 sumR = vec3(0.f);
 	vec3 sumM = vec3(0.f);
 	vec2 opticalDepth = vec2(0.f);
 
-	processRay(sumR, sumM, opticalDepth, sunDirection, dir, atmosphere_minus_shadow.xy, center, diffuse, uneclipsed);
-	skipRay(opticalDepth, dir, atmosphere_minus_shadow.yz, center);
-	processRay(sumR, sumM, opticalDepth, sunDirection, dir, atmosphere_minus_shadow.zw, center, diffuse, uneclipsed);
+	vec4 segment = getRaySegment(sunDirection, dir, center);
+
+	processRay(sumR, sumM, opticalDepth, sunDirection, dir, segment.xy, center, diffuse, uneclipsed, texture_LUT);
+	skipRay(opticalDepth, dir, segment.yz, center, texture_LUT);
+	processRay(sumR, sumM, opticalDepth, sunDirection, dir, segment.zw, center, diffuse, uneclipsed, texture_LUT);
 
 	float mu = dot(dir, sunDirection); // mu in the paper which is the cosine of the angle between the sun direction and the ray direction
 	float phaseR = rayleighPhaseFunction(mu);
@@ -258,3 +355,13 @@ vec3 computeIncidentLight(const in vec3 sunDirection, const in vec3 dir, const i
 	vec3 ret = (sumR * betaR * phaseR + sumM * betaM * phaseM);
 	return ret;
 }
+
+vec3 calculateAtmosphereColor(const in vec4 planet, const in vec4 atmosphere, const in vec4 lightColor, const in vec3 lightDir, const in vec3 rayStart, const in vec3 rayDir, const in float uneclipsed, const in sampler2D texture_LUT)
+{
+    // rayStart is already multiplied by planet radius
+    vec3 planetPosition = planet.xyz * planet.w + rayStart;
+
+    vec3 atmospherePosition = atmosphere.xyz * atmosphere.w;
+
+    return computeIncidentLight(lightDir, rayDir, planetPosition, lightColor, uneclipsed, texture_LUT);
+}