From ce424215c708ba57d186be8982384c90b6c5620a Mon Sep 17 00:00:00 2001 From: Oleg Solovyov Date: Mon, 13 Nov 2023 12:32:12 +0300 Subject: [PATCH] rayleigh: move common subroutines to AtmosphereParameters --- data/shaders/opengl/basesphere_uniforms.glsl | 2 + data/shaders/opengl/rayleigh_sky.frag | 45 +++------------- src/BaseSphere.cpp | 6 ++- src/galaxy/AtmosphereParameters.h | 2 + src/galaxy/SystemBody.cpp | 55 ++++++++++++++++++++ src/galaxy/SystemBody.h | 4 ++ 6 files changed, 76 insertions(+), 38 deletions(-) diff --git a/data/shaders/opengl/basesphere_uniforms.glsl b/data/shaders/opengl/basesphere_uniforms.glsl index f10ea051e07..88abb32d033 100644 --- a/data/shaders/opengl/basesphere_uniforms.glsl +++ b/data/shaders/opengl/basesphere_uniforms.glsl @@ -13,6 +13,8 @@ layout(std140) uniform BaseSphereData { float geosphereAtmosFogDensity; // TODO documentation float geosphereAtmosInvScaleHeight; // TODO documentation vec4 atmosColor; + vec3 coefficientsR; + vec3 coefficientsM; // Eclipse data Eclipse eclipse; diff --git a/data/shaders/opengl/rayleigh_sky.frag b/data/shaders/opengl/rayleigh_sky.frag index ffc2011425e..6b377025789 100644 --- a/data/shaders/opengl/rayleigh_sky.frag +++ b/data/shaders/opengl/rayleigh_sky.frag @@ -53,27 +53,6 @@ void findClosestHeight(out float h, out float t, const in vec3 orig, const in ve t = dot(tangent, dir); } -// 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) - radius; - density += step * exp(-h / baricStep); - } - - return density; -} - // error function float erf(const in float x) { @@ -116,18 +95,6 @@ float predictDensityIn(const in float radius, const in float atmosphereHeight, c } } -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 = 1.f; - 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, const in vec3 dir, const in vec3 center, in float tmin, in float tmax) { vec3 betaR = vec3(3.8e-6f, 13.5e-6f, 33.1e-6f); @@ -156,10 +123,14 @@ vec3 computeIncidentLight(const in vec3 sunDirection, const in vec3 dir, const i 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); + + kR = coefficientsR.x; + bR = coefficientsR.y; + cR = coefficientsR.z; + + kM = coefficientsM.x; + bM = coefficientsM.y; + cM = coefficientsM.z; for (int i = 0; i < numSamples; ++i) { vec3 samplePosition = vec3(tCurrent + segmentLength * 0.5f) * dir; diff --git a/src/BaseSphere.cpp b/src/BaseSphere.cpp index f9f7e8c6a6f..7a41e8d35dd 100644 --- a/src/BaseSphere.cpp +++ b/src/BaseSphere.cpp @@ -19,6 +19,8 @@ struct BaseSphereDataBlock { float geosphereAtmosFogDensity; float geosphereAtmosInvScaleHeight; Color4f atmosColor; + vector3f coefficientsR; + vector3f coefficientsM; // Eclipse struct data alignas(16) vector3f shadowCentreX; @@ -28,7 +30,7 @@ struct BaseSphereDataBlock { alignas(16) vector3f lrad; alignas(16) vector3f sdivlrad; }; -static_assert(sizeof(BaseSphereDataBlock) == 144, ""); +static_assert(sizeof(BaseSphereDataBlock) == 176, ""); BaseSphere::BaseSphere(const SystemBody *body) : m_sbody(body), @@ -97,6 +99,8 @@ void BaseSphere::SetMaterialParameters(const matrix4x4d &trans, const float radi matData.geosphereAtmosFogDensity = ap.atmosDensity; matData.geosphereAtmosInvScaleHeight = ap.atmosInvScaleHeight; matData.atmosColor = ap.atmosCol.ToColor4f(); + matData.coefficientsR = ap.rayleighCoefficients; + matData.coefficientsM = ap.mieCoefficients; // we handle up to three shadows at a time auto it = shadows.cbegin(), itEnd = shadows.cend(); diff --git a/src/galaxy/AtmosphereParameters.h b/src/galaxy/AtmosphereParameters.h index 16afd0c10fe..2fa37c001e4 100644 --- a/src/galaxy/AtmosphereParameters.h +++ b/src/galaxy/AtmosphereParameters.h @@ -11,6 +11,8 @@ struct AtmosphereParameters { float planetRadius; Color atmosCol; vector3d center; + vector3f rayleighCoefficients; + vector3f mieCoefficients; }; #endif // ATMOSPHEREPARAMETERS_H_INCLUDED diff --git a/src/galaxy/SystemBody.cpp b/src/galaxy/SystemBody.cpp index 9dab74cc4b1..356668dba04 100644 --- a/src/galaxy/SystemBody.cpp +++ b/src/galaxy/SystemBody.cpp @@ -262,6 +262,56 @@ bool SystemBody::IsScoopable() const return false; } +// error function: used in scattering prediction +static double erf(const double x) +{ + double a = 0.14001228904; + double four_over_pi = 1.27323954; + + double x2 = x*x; + double 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; +} + +// all input units are in km +double SystemBody::computeDensity(const double radius, const double atmosphereHeight, const double h, const double baricStep) const +{ + int numSamples = 16; + double totalHeight = radius + atmosphereHeight; + double minHeight = radius + h; + double tmax = sqrt(totalHeight*totalHeight - minHeight*minHeight); // maximum t + double tmin = -tmax; + double length = tmax - tmin; + double step = length / numSamples; + + double density = 0.0; + for (int i = 0; i < numSamples; ++i) { + double t = tmin + step * (i + 0.5); + double h = sqrt(minHeight*minHeight + t*t) - radius; + density += step * exp(-h / baricStep); + } + + return density; +} + +// all input units are in km +vector3f SystemBody::GetCoefficients(const double radius, const double atmHeight, const double baricStep) const +{ + float k, b, c; + k = computeDensity(radius, atmHeight, 0.f, baricStep); + b = log(k / computeDensity(radius, atmHeight, 1.f, baricStep)); + + float s = 1.f; + float sHeight = sqrt(radius*radius + s*s) - radius; + float c1 = exp(-sHeight / baricStep) * s; + float c2 = k * (erf(s) - erf(0)); + c = c1 / c2; + return vector3f(k, b, c); +} + // Calculate parameters used in the atmospheric model for shaders AtmosphereParameters SystemBody::CalcAtmosphereParams() const { @@ -320,6 +370,11 @@ AtmosphereParameters SystemBody::CalcAtmosphereParams() const params.planetRadius = static_cast(radiusPlanet_in_m); + const float radiusPlanet_in_km = radiusPlanet_in_m / 1000; + const float atmosHeight_in_km = radiusPlanet_in_km * (params.atmosRadius - 1); + params.rayleighCoefficients = GetCoefficients(radiusPlanet_in_km, atmosHeight_in_km, 7994); + params.mieCoefficients = GetCoefficients(radiusPlanet_in_km, atmosHeight_in_km, 1200); + return params; } diff --git a/src/galaxy/SystemBody.h b/src/galaxy/SystemBody.h index 70f5b3b372c..ae2e275b658 100644 --- a/src/galaxy/SystemBody.h +++ b/src/galaxy/SystemBody.h @@ -292,6 +292,10 @@ class SystemBody : public RefCounted, public SystemBodyType, protected SystemBod // Calculate atmosphere density at given altitude and pressure (kg/m^3) double GetAtmDensity(double altitude, double pressure) const; + // for rayleigh scattering + double computeDensity(const double radius, const double atmosphereHeight, const double h, const double baricStep) const; + vector3f GetCoefficients(const double radius, const double atmHeight, const double baricStep) const; + AtmosphereParameters CalcAtmosphereParams() const; bool IsScoopable() const;