Skip to content

Commit

Permalink
rayleigh: move common subroutines to AtmosphereParameters
Browse files Browse the repository at this point in the history
  • Loading branch information
Mc-Pain committed Nov 13, 2023
1 parent d5eb9bd commit ce42421
Show file tree
Hide file tree
Showing 6 changed files with 76 additions and 38 deletions.
2 changes: 2 additions & 0 deletions data/shaders/opengl/basesphere_uniforms.glsl
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
45 changes: 8 additions & 37 deletions data/shaders/opengl/rayleigh_sky.frag
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down
6 changes: 5 additions & 1 deletion src/BaseSphere.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@ struct BaseSphereDataBlock {
float geosphereAtmosFogDensity;
float geosphereAtmosInvScaleHeight;
Color4f atmosColor;
vector3f coefficientsR;
vector3f coefficientsM;

// Eclipse struct data
alignas(16) vector3f shadowCentreX;
Expand All @@ -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),
Expand Down Expand Up @@ -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();
Expand Down
2 changes: 2 additions & 0 deletions src/galaxy/AtmosphereParameters.h
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@ struct AtmosphereParameters {
float planetRadius;
Color atmosCol;
vector3d center;
vector3f rayleighCoefficients;
vector3f mieCoefficients;
};

#endif // ATMOSPHEREPARAMETERS_H_INCLUDED
55 changes: 55 additions & 0 deletions src/galaxy/SystemBody.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down Expand Up @@ -320,6 +370,11 @@ AtmosphereParameters SystemBody::CalcAtmosphereParams() const

params.planetRadius = static_cast<float>(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;
}

Expand Down
4 changes: 4 additions & 0 deletions src/galaxy/SystemBody.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit ce42421

Please sign in to comment.