experimental: pre-calculate integral density

Author: Mc-Pain

data/shaders/opengl/rayleigh_accurate.frag

@@ -8,6 +8,8 @@
 
 uniform int NumShadows;
 uniform sampler2D scatterLUT;
+uniform sampler2D rayleighLUT;
+uniform sampler2D mieLUT;
 
 in vec4 varyingEyepos;
 in vec4 vertexColor;
@@ -54,4 +56,21 @@ void main(void)
 
 	frag_color = toSRGB(1 - exp(-color));
 
+#if 1
+    frag_color = vec4(0.f, 0.f, 0.f, 1.f);
+
+    // gl_FragCoord: [0-1280], [0, 720]
+
+    frag_color.x = texelFetch(rayleighLUT, ivec2(gl_FragCoord.x - 512, gl_FragCoord.y - 232), 0).x;
+    //frag_color.y = texture(mieLUT, vec2(gl_FragCoord.x / (64 * 2), gl_FragCoord.y / (36 * 256))).x;
+    // draw border
+    bool outer_range = (136 < gl_FragCoord.x) || (620 < gl_FragCoord.y);
+    bool inner_range = (gl_FragCoord.x < 128) && (gl_FragCoord.y < 612);
+
+    if (!inner_range && !outer_range) {
+        frag_color.xyz = vec3(0.f, 0.f, 0.25f);
+    }
+
+    //frag_color.xyz = texelFetch(scatterLUT, ivec2(gl_FragCoord.x / 64, gl_FragCoord.y / 36), 0).x == 0 ? vec3(0.1f, 0.0f, 0.1f) : vec3(0.0f, 0.0f, 0.1f);
+#endif
 }

data/shaders/opengl/rayleigh_accurate.shaderdef

@@ -4,6 +4,8 @@
 Shader rayleigh_accurate
 
 Texture sampler2D scatterLUT
+Texture sampler2D rayleighLUT
+Texture sampler2D mieLUT
 
 Buffer LightData binding=0
 Buffer DrawData binding=1

src/BaseSphere.cpp

@@ -143,15 +143,46 @@ void BaseSphere::SetMaterialParameters(const matrix4x4d &trans, const float radi
 		static_cast<void *>(buf.get()), texSize,
 		Graphics::TEXTURE_R32);
 
+	std::unique_ptr<float, FreeDeleter> bufR(
+		static_cast<float *>(malloc((DENSITY_STEPS + 1) * (DENSITY_STEPS + 1) * sizeof(float))));
+	std::unique_ptr<float, FreeDeleter> bufM(
+		static_cast<float *>(malloc((DENSITY_STEPS + 1) * (DENSITY_STEPS + 1) * sizeof(float))));
+
+	for (int i = 0; i <= DENSITY_STEPS; ++i) {
+		for (int j = 0; j <= DENSITY_STEPS; ++j) {
+			int index = i * (DENSITY_STEPS + 1) + j;
+
+			bufR.get()[index] = 0.5 + (ap.logDensityMapR[index] / 256);
+			bufM.get()[index] = 0.5 + (ap.logDensityMapM[index] / 256);
+		}
+	}
+
+	const vector3f texSizeR(DENSITY_STEPS + 1, DENSITY_STEPS + 1, 0.f);
+	const vector3f texSizeM(DENSITY_STEPS + 1, DENSITY_STEPS + 1, 0.f);
+	const Graphics::TextureDescriptor texDescR(
+		Graphics::TEXTURE_R32, texSizeR, Graphics::LINEAR_CLAMP, false, true, true, 0, Graphics::TEXTURE_2D);
+	const Graphics::TextureDescriptor texDescM(
+		Graphics::TEXTURE_R32, texSizeM, Graphics::LINEAR_CLAMP, false, true, true, 0, Graphics::TEXTURE_2D);
+
+	m_scatteringTextureR.Reset(Pi::renderer->CreateTexture(texDescR));
+	m_scatteringTextureR->Update(static_cast<void *>(bufR.get()), texSizeR, Graphics::TEXTURE_R32);
+
+	m_scatteringTextureM.Reset(Pi::renderer->CreateTexture(texDescM));
+	m_scatteringTextureM->Update(static_cast<void *>(bufM.get()), texSizeM, Graphics::TEXTURE_R32);
+
 	// FIXME: these two should share the same buffer data instead of making two separate allocs
 	m_surfaceMaterial->SetBufferDynamic(s_baseSphereData, &matData);
 	m_surfaceMaterial->SetPushConstant(s_numShadows, int(shadows.size()));
 	m_surfaceMaterial->SetTexture("scatterLUT"_hash, m_scatteringTexture.Get());
+	m_surfaceMaterial->SetTexture("rayleighLUT"_hash, m_scatteringTextureR.Get());
+	m_surfaceMaterial->SetTexture("mieLUT"_hash, m_scatteringTextureM.Get());
 
 	if (m_atmosphereMaterial.Valid() && ap.atmosDensity > 0.0) {
 		m_atmosphereMaterial->SetBufferDynamic(s_baseSphereData, &matData);
 		m_atmosphereMaterial->SetPushConstant(s_numShadows, int(shadows.size()));
 		m_atmosphereMaterial->SetTexture("scatterLUT"_hash, m_scatteringTexture.Get());
+		m_atmosphereMaterial->SetTexture("rayleighLUT"_hash, m_scatteringTextureR.Get());
+		m_atmosphereMaterial->SetTexture("mieLUT"_hash, m_scatteringTextureM.Get());
 	}
 }
 

src/BaseSphere.h

@@ -66,6 +66,8 @@ class BaseSphere {
 	AtmosphereParameters m_atmosphereParameters;
 
 	RefCountedPtr<Graphics::Texture> m_scatteringTexture;
+	RefCountedPtr<Graphics::Texture> m_scatteringTextureR;
+	RefCountedPtr<Graphics::Texture> m_scatteringTextureM;
 
 private:
 	static void ResetAtmosphereGeometry(Graphics::Renderer *renderer);

src/galaxy/AtmosphereParameters.h

@@ -17,6 +17,8 @@ struct AtmosphereParameters {
 	vector3f mieCoefficients;
 	vector2f scaleHeight;
 	vector2f logDensityMap[DENSITY_STEPS + 1];
+	float logDensityMapR[(DENSITY_STEPS + 1)*(DENSITY_STEPS + 1)];
+	float logDensityMapM[(DENSITY_STEPS + 1)*(DENSITY_STEPS + 1)];
 };
 
 #endif // ATMOSPHEREPARAMETERS_H_INCLUDED

src/galaxy/StarSystemGenerator.cpp

@@ -429,7 +429,9 @@ void StarSystemCustomGenerator::CustomGetKidsOf(RefCountedPtr<StarSystem::Genera
 
 		if (kid->GetType() != SystemBody::TYPE_STARPORT_SURFACE) {
 			if (kid->GetSuperType() == SystemBody::SUPERTYPE_STARPORT) {
-				fixed lowestOrbit = fixed().FromDouble(parent->CalcAtmosphereParams().atmosRadius + 500000.0 / EARTH_RADIUS);
+				const double parentAtmosRadius = 1.0 + parent->GetAtmRadius() / parent->GetRadius();
+
+				fixed lowestOrbit = fixed().FromDouble(parentAtmosRadius + 500000.0 / EARTH_RADIUS);
 				if (kid->GetOrbit().GetSemiMajorAxis() < lowestOrbit.ToDouble()) {
 					Error("%s's orbit is too close to its parent (%.2f/%.2f)", kid->m_name.c_str(), kid->GetOrbit().GetSemiMajorAxis(), lowestOrbit.ToFloat());
 				}
@@ -1661,7 +1663,9 @@ void PopulateStarSystemGenerator::PopulateAddStations(SystemBody *sbody, StarSys
 
 	if (sbody->GetPopulationAsFixed() < fixed(1, 1000)) return;
 	fixed orbMaxS = fixed(1, 4) * fixed(CalcHillRadius(sbody));
-	fixed orbMinS = fixed().FromDouble((sbody->CalcAtmosphereParams().atmosRadius + +500000.0 / EARTH_RADIUS)) * AU_EARTH_RADIUS;
+
+	const double sbodyAtmosRadius = 1.0 + sbody->GetAtmRadius() / sbody->GetRadius();
+	fixed orbMinS = fixed().FromDouble((sbodyAtmosRadius + +500000.0 / EARTH_RADIUS)) * AU_EARTH_RADIUS;
 	if (sbody->GetNumChildren() > 0)
 		orbMaxS = std::min(orbMaxS, fixed(1, 2) * sbody->GetChildren()[0]->GetOrbMinAsFixed());
 

src/galaxy/SystemBody.cpp

@@ -413,6 +413,33 @@ vector3f SystemBody::GetCoefficients(const double radius, const double atmHeight
 	return vector3f(k, b, c);
 }
 
+/*
+ * given ray and a circle, calculate points where ray intersects a circle
+ * r - radius of a circle
+ * v - ray start, given circle center is at origin
+ * dir - ray direction, pre-normalized
+ * x1, x2 - solutions
+ *
+ * returns false if circle is not intersected
+ */
+bool rayCircleIntersection(double *x1, double *x2, const vector2d v, const double radius, const vector2d dir)
+{
+	double b = -(v.x * dir.x + v.y * dir.y);
+	double det = (b * b) - v.LengthSqr() + (radius * radius);
+	//Output("DEBUG: v = [%lf, %lf], radius = %lf, dir = [%lf, %lf], b = %lf, det = %lf\n", v.x, v.y, radius, dir.x, dir.y, b, det);
+
+	if (det >= 0) {
+		double sdet = sqrt(det);
+
+		*x1 = b - sdet;
+		*x2 = b + sdet;
+
+		return true;
+	} else {
+		return false;
+	}
+}
+
 // Calculate parameters used in the atmospheric model for shaders
 AtmosphereParameters SystemBody::CalcAtmosphereParams() const
 {
@@ -493,6 +520,72 @@ AtmosphereParameters SystemBody::CalcAtmosphereParams() const
 		params.logDensityMap[i] = vector2f(rLogDensity, mLogDensity);
 	}
 
+	// tests
+	double tx1, tx2;
+	if (GetName() == "Mars") {
+		rayCircleIntersection(&tx1, &tx2, vector2d(0.0, 4.0), 5.0, vector2d(0.0, 1.0));
+		//Output("%s: test intersections %lf, %lf\n", GetName().c_str(), tx1, tx2);
+		rayCircleIntersection(&tx1, &tx2, vector2d(0.0, 4.0), 5.0, vector2d(1.0, 0.0));
+		//Output("%s: test intersections %lf, %lf\n", GetName().c_str(), tx1, tx2);
+	}
+
+	Output("%s: generating density model\n", GetName());
+
+	// pre-calculate density LUT based on current height and angle to horizon
+	for (int i = 0; i <= DENSITY_STEPS; ++i) {
+		double startHeight = i * atmosHeight / DENSITY_STEPS;
+
+		for (int j = 0; j <= DENSITY_STEPS; ++j) {
+			// pitch in radians
+			double pitch = M_PI * ((1.0 * j / DENSITY_STEPS) - 0.5);
+			vector2d direction = vector2d(1.0, 0.0).Rotate(pitch); // calculate in 2d space
+			vector2d tCurrent = vector2d(0.0, startHeight + radiusPlanet_in_m);
+
+			double atm1 = 0.0, atm2 = 0.0;
+			bool atmIntersection = rayCircleIntersection(&atm1, &atm2, tCurrent, radiusPlanet_in_m * params.atmosRadius, direction);
+
+			if (!atmIntersection)
+				continue;
+
+			double planet1 = 0.0, planet2 = 0.0;
+			bool planetIntersection = rayCircleIntersection(&planet1, &planet2, tCurrent, radiusPlanet_in_m - 0.01, direction);
+
+			double startRay = 0.0, finishRay = atm2;
+			if (planetIntersection) {
+				if (planet1 >= 0) finishRay = planet1; // we hit the planet
+
+				// otherwise we would hit the planet while tracing *backwards*
+			}
+
+			int numSamples = 16;
+			double segment = (finishRay - startRay) / numSamples;
+			double rDensity = 0.0, mDensity = 0.0;
+			for (int k = 0; k < numSamples; ++k) {
+				vector2d samplePosition = tCurrent + (segment * 0.5) * direction;
+				double sampleHeight = samplePosition.Length() - radiusPlanet_in_m;
+
+				double sampleDensityR = GetAtmDensity(sampleHeight, GetAtmPressure(sampleHeight));
+				double sampleDensityM = GetAtmDensity(6.66 * sampleHeight, GetAtmPressure(6.66 * sampleHeight));
+				if (std::isnan(sampleDensityM))
+					sampleDensityM = 0.0;
+
+				rDensity += segment * sampleDensityR;
+				mDensity += segment * sampleDensityM;
+
+				tCurrent += direction * segment;
+			}
+
+			double rLogDensity, mLogDensity;
+
+			rLogDensity = rDensity > 0.0 ? log(rDensity) : -128.0;
+			mLogDensity = mDensity > 0.0 ? log(mDensity) : -128.0;
+
+			if (GetName() == "Mars") {
+				//Output("%s: calculated density for ray (height = %lf, pitch = %lf): log(r) = %lf, log(m) = %lf\n\n", GetName().c_str(), startHeight, pitch, rLogDensity, mLogDensity);
+			}
+		}
+	}
+
 	return params;
 }