Advanced eclipses

config/base/scene/simple_desktop.json

@@ -282,9 +282,13 @@
   "graphics": {
     "enableHDR": true,
     "eclipseShadowMaps": {
-      "Earth": {},
+      "Earth": {
+        "texture":"../share/resources/textures/earthShadow.tif"
+      },
       "Moon": {},
-      "Mars": {},
+      "Mars": {
+        "texture":"../share/resources/textures/marsShadow.tif"
+      },
       "Phobos": {},
       "Deimos": {},
       "Jupiter": {},
@@ -739,6 +743,7 @@
       "atmospheres": {
         "Earth": {
           "cloudTexture": "../share/resources/textures/earth-clouds.jpg",
+          "limbLuminanceTexture": "../share/resources/textures/earthLimbLuminance.tif",
           "topAltitude": 80000,
           "bottomAltitude": -100,
           "model": "Bruneton",

plugins/csp-atmospheres/README.md

@@ -145,15 +145,21 @@ They are not physically based but provide some plausible results.
 The Bruneton model is significantly more advanced.
 It precomputes multiple scattering and is based on [this open-source implementation](https://github.com/ebruneton/precomputed_atmospheric_scattering) (see also the corresponding [Paper](https://inria.hal.science/inria-00288758/en)).
 
-Similar to the `CosmoScoutVR` model, the original implementation by Eric Bruneton uses Rayleigh scattering for molecules and the Cornette-Shanks phase function for aerosols.
-We generalized this implementation by loading phase functions, extinction coefficients, and particle density distributions from CSV files.
+We have significantly extended the original implementation to allow for more flexibility.
+See our paper [Physically Based Real-Time Rendering of Atmospheres using Mie Theory](https://onlinelibrary.wiley.com/doi/full/10.1111/cgf.15010) for more details.
+
+As a first change, we now load phase functions, extinction coefficients, and particle density distributions from CSV files.
 This allows us to simulate arbitrary particle types.
 In particular, we can now use Mie Theory to precompute the scattering behaviour of a wide variety of particle types, including for instance Martian dust.
 
+Next, our implementation can compute refraction.
+This allows for displaced horizons and the simulation of astronomical refraction.
+This is also used for computing light entering the eclipse shadows of celestial bodies.
+
 Another change to the original implementation is that we put the precomputation of the atmospheric scattering into a separate executable.
 This allows us to perform the preprocessing offline with a much higher fidelity than what would be possible during application startup.
 
-As a consequence to the changes mentioned above, **two preprocessing steps are required to use this model**.
+There are **two preprocessing steps are required to use this model**.
 
 #### Preprocessing Step 1: Precompute the Particle-Scattering CSV Tables
 
@@ -197,7 +203,7 @@ Once the multiple scattering textures are precomputed, the configuration for the
 </details>
 
 <details>
-<summary>Example Configuration for Mars (Realistic)</summary>
+<summary>Example Configuration for Mars</summary>
 
 ```javascript
 "Mars": {

plugins/csp-atmospheres/REUSE.toml

@@ -21,6 +21,12 @@ precedence = "aggregate"
 SPDX-FileCopyrightText = "German Aerospace Center (DLR) <[email protected]>"
 SPDX-License-Identifier = "CC0-1.0"
 
+[[annotations]]
+path = "scattering-table-generator/ior-settings/**"
+precedence = "aggregate"
+SPDX-FileCopyrightText = "German Aerospace Center (DLR) <[email protected]>"
+SPDX-License-Identifier = "CC0-1.0"
+
 [[annotations]]
 path = "scattering-table-generator/output/**"
 precedence = "aggregate"

plugins/csp-atmospheres/bruneton-preprocessor/Metadata.cpp

@@ -14,6 +14,7 @@ void to_json(nlohmann::json& j, Metadata const& o) {
   cs::core::Settings::serialize(j, "sunIlluminance", o.mSunIlluminance);
   cs::core::Settings::serialize(j, "scatteringTextureNuSize", o.mScatteringTextureNuSize);
   cs::core::Settings::serialize(j, "maxSunZenithAngle", o.mMaxSunZenithAngle);
+  cs::core::Settings::serialize(j, "refraction", o.mRefraction);
 }
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////

plugins/csp-atmospheres/bruneton-preprocessor/Metadata.hpp

@@ -26,6 +26,9 @@ struct Metadata {
 
   /// The maximum Sun zenith angle for which atmospheric scattering was be precomputed.
   float mMaxSunZenithAngle{};
+
+  /// Whether refraction was used during preprocessing.
+  bool mRefraction{};
 };
 
 void to_json(nlohmann::json& j, Metadata const& o);

plugins/csp-atmospheres/bruneton-preprocessor/Params.cpp

@@ -31,13 +31,18 @@ void from_json(nlohmann::json const& j, Params& o) {
   cs::core::Settings::deserialize(j, "molecules", o.mMolecules);
   cs::core::Settings::deserialize(j, "aerosols", o.mAerosols);
   cs::core::Settings::deserialize(j, "ozone", o.mOzone);
+  cs::core::Settings::deserialize(j, "ior", o.mRefractiveIndex);
   cs::core::Settings::deserialize(j, "minAltitude", o.mMinAltitude);
   cs::core::Settings::deserialize(j, "maxAltitude", o.mMaxAltitude);
+  cs::core::Settings::deserialize(j, "refraction", o.mRefraction);
   cs::core::Settings::deserialize(j, "groundAlbedo", o.mGroundAlbedo);
   cs::core::Settings::deserialize(j, "multiScatteringOrder", o.mMultiScatteringOrder);
   cs::core::Settings::deserialize(j, "sampleCountOpticalDepth", o.mSampleCountOpticalDepth);
+  cs::core::Settings::deserialize(j, "stepSizeOpticalDepth", o.mStepSizeOpticalDepth);
   cs::core::Settings::deserialize(j, "sampleCountSingleScattering", o.mSampleCountSingleScattering);
+  cs::core::Settings::deserialize(j, "stepSizeSingleScattering", o.mStepSizeSingleScattering);
   cs::core::Settings::deserialize(j, "sampleCountMultiScattering", o.mSampleCountMultiScattering);
+  cs::core::Settings::deserialize(j, "stepSizeMultiScattering", o.mStepSizeMultiScattering);
   cs::core::Settings::deserialize(
       j, "sampleCountScatteringDensity", o.mSampleCountScatteringDensity);
   cs::core::Settings::deserialize(

plugins/csp-atmospheres/bruneton-preprocessor/Params.hpp

@@ -61,6 +61,10 @@ struct Params {
   ScatteringComponent               mAerosols;
   std::optional<AbsorbingComponent> mOzone;
 
+  /// To compute the refraction of light in the atmosphere, the refractive index of the atmosphere
+  /// is needed. For increased precision, the refractive index is stored as n-1.
+  float mRefractiveIndex = 0.0002777F;
+
   /// The wavelength values, in nanometers, and sorted in increasing order, for which the
   /// phase functions and extinction coefficients in the atmosphere components are given.
   std::vector<float> mWavelengths;
@@ -69,22 +73,36 @@ struct Params {
   float mMinAltitude = 6371000.F;
   float mMaxAltitude = 6471000.F;
 
+  /// Refract the light when it travels through the atmosphere. This will produce an additional
+  /// look-up texture in the same parameter space as the transmittance texture. For each sample,
+  /// it contains the wavelength-dependent angular deviation of the light ray due to refraction.
+  cs::utils::DefaultProperty<bool> mRefraction{true};
+
   /// The average reflectance of the ground used during multiple scattering.
   cs::utils::DefaultProperty<float> mGroundAlbedo{0.1F};
 
   /// The number of multiple scattering events to precompute. Use zero for single-scattering only.
   cs::utils::DefaultProperty<int32_t> mMultiScatteringOrder{4};
 
-  /// The number of samples to evaluate when precomputing the optical depth.
+  /// The number of samples to evaluate when precomputing the optical depth. If refraction is used,
+  /// the algorithm uses a fixed step size instead of a fixed number of samples. So if mRefraction
+  /// is true, mStepSizeOpticalDepth (in meters) will be used instead of mSampleCountOpticalDepth.
   cs::utils::DefaultProperty<int32_t> mSampleCountOpticalDepth{500};
+  cs::utils::DefaultProperty<int32_t> mStepSizeOpticalDepth{10000};
 
   /// The number of samples to evaluate when precomputing the single scattering. Larger values
-  /// improve the sampling of thin atmospheric layers.
+  /// improve the sampling of thin atmospheric layers. If refraction is used, the algorithm uses a
+  /// fixed step size instead of a fixed number of samples. So if mRefraction is true,
+  /// mStepSizeSingleScattering (in meters) will be used instead of mSampleCountSingleScattering.
   cs::utils::DefaultProperty<int32_t> mSampleCountSingleScattering{50};
+  cs::utils::DefaultProperty<int32_t> mStepSizeSingleScattering{10000};
 
   /// The number of samples to evaluate when precomputing the multiple scattering. Larger values
-  /// tend to darken the horizon for thick atmospheres.
+  /// tend to darken the horizon for thick atmospheres. If refraction is used, the algorithm uses a
+  /// fixed step size instead of a fixed number of samples. So if mRefraction is true,
+  /// mStepSizeMultiScattering (in meters) will be used instead of mSampleCountMultiScattering.
   cs::utils::DefaultProperty<int32_t> mSampleCountMultiScattering{50};
+  cs::utils::DefaultProperty<int32_t> mStepSizeMultiScattering{10000};
 
   /// The number of samples to evaluate when precomputing the scattering density. Larger values
   /// spread out colors in the sky.

plugins/csp-atmospheres/bruneton-preprocessor/Preprocessor.cpp

@@ -31,6 +31,13 @@
 // density distributions are now loaded from CSV files and then later sampled from textures. We also
 // store photometric values instead of radiometric values in the final textures.
 
+// Also, we added the possibility to compute the refraction of light rays through the atmosphere.
+// If enabled, an additional texture is generated which contains the angular deviation of the light
+// rays as well as their closest approach to the planet's surface (this can be negative if the light
+// ray intersects the planet's surface).
+// The texture uses the same parametrization as the transmittance texture. The values are computed
+// by the kComputeTransmittanceShader.
+
 // Below, we will indicate for each group of function whether something has been changed and a link
 // to the original explanations of the methods by Eric Bruneton.
 
@@ -185,9 +192,11 @@ constexpr float XYZ_TO_SRGB[9] = {
 // An explanation of the following shaders is available online:
 // https://ebruneton.github.io/precomputed_atmospheric_scattering/atmosphere/model.cc.html#shaders
 
-// The only functional difference is that the kAtmosphereShader does not provide the radiance API
+// The only functional differences are that the kAtmosphereShader does not provide the radiance API
 // anymore as it is not required by CosmoScout VR. Also, the shadow_length parameters have been
 // removed and the GetSunAndSkyIlluminance() does not require the surface normal anymore.
+// Lastly, the kComputeTransmittanceShader has been extended to compute the refraction of light rays
+// through the atmosphere.
 
 const char kVertexShader[] = R"(
   #version 330
@@ -227,8 +236,20 @@ const char kGeometryShader[] = R"(
 const char kComputeTransmittanceShader[] = R"(
   layout(location = 0) out vec3 oTransmittance;
 
+#if COMPUTE_REFRACTION
+  layout(location = 1) out vec3 oThetaDeviationContactRadius;
+#endif
+
   void main() {
-    oTransmittance = computeTransmittanceToTopAtmosphereBoundaryTexture(ATMOSPHERE, gl_FragCoord.xy);
+    #if COMPUTE_REFRACTION
+      float contactRadius;
+      float thetaDeviation;
+      oTransmittance = computeTransmittanceToTopAtmosphereBoundaryTexture(ATMOSPHERE, gl_FragCoord.xy, thetaDeviation, contactRadius);
+      oThetaDeviationContactRadius = vec3(thetaDeviation, contactRadius, 0.0);
+
+    #else
+      oTransmittance = computeTransmittanceToTopAtmosphereBoundaryTexture(ATMOSPHERE, gl_FragCoord.xy);
+    #endif
   }
 )";
 
@@ -313,7 +334,7 @@ const char kComputeMultipleScatteringShader[] = R"(
 
   void main() {
     float nu;
-    oDeltaMultipleScattering = computeMultipleScatteringTexture(uTransmittanceTexture,
+    oDeltaMultipleScattering = computeMultipleScatteringTexture(ATMOSPHERE, uTransmittanceTexture,
                                                                 uScatteringDensityTexture,
                                                                 vec3(gl_FragCoord.xy, uLayer + 0.5),
                                                                 nu);
@@ -674,6 +695,7 @@ Preprocessor::Preprocessor(Params params)
   mMetadata.mSunIlluminance          = glm::vec3(sunKR, sunKG, sunKB);
   mMetadata.mScatteringTextureNuSize = mParams.mScatteringTextureNuSize.get();
   mMetadata.mMaxSunZenithAngle       = mParams.mMaxSunZenithAngle.get();
+  mMetadata.mRefraction              = mParams.mRefraction.get();
 
   // A lambda that creates a GLSL header containing our atmosphere computation functions,
   // specialized for the given atmosphere parameters and for the 3 wavelengths in 'lambdas'.
@@ -689,8 +711,9 @@ Preprocessor::Preprocessor(Params params)
   // clang-format off
   mGlslHeaderFactory = [=](glm::vec3 const& lambdas) {
     return
-      "#version 330\n" +
+      "#version 400\n" +
       definitions +
+      "#define COMPUTE_REFRACTION "                   + cs::utils::toString(mParams.mRefraction) + "\n" +
       "const int TRANSMITTANCE_TEXTURE_WIDTH = "      + cs::utils::toString(mParams.mTransmittanceTextureWidth) + ";\n" +
       "const int TRANSMITTANCE_TEXTURE_HEIGHT = "     + cs::utils::toString(mParams.mTransmittanceTextureHeight) + ";\n" +
       "const int SCATTERING_TEXTURE_R_SIZE = "        + cs::utils::toString(mParams.mScatteringTextureRSize) + ";\n" +
@@ -700,12 +723,16 @@ Preprocessor::Preprocessor(Params params)
       "const int IRRADIANCE_TEXTURE_WIDTH = "         + cs::utils::toString(mParams.mIrradianceTextureWidth) + ";\n" +
       "const int IRRADIANCE_TEXTURE_HEIGHT = "        + cs::utils::toString(mParams.mIrradianceTextureHeight) + ";\n" +
       "const int SAMPLE_COUNT_OPTICAL_DEPTH = "       + cs::utils::toString(mParams.mSampleCountOpticalDepth) + ";\n" +
+      "const int STEP_SIZE_OPTICAL_DEPTH = "          + cs::utils::toString(mParams.mStepSizeOpticalDepth) + ";\n" +
       "const int SAMPLE_COUNT_SINGLE_SCATTERING = "   + cs::utils::toString(mParams.mSampleCountSingleScattering) + ";\n" +
+      "const int STEP_SIZE_SINGLE_SCATTERING = "      + cs::utils::toString(mParams.mStepSizeSingleScattering) + ";\n" +
       "const int SAMPLE_COUNT_SCATTERING_DENSITY = "  + cs::utils::toString(mParams.mSampleCountScatteringDensity) + ";\n" +
       "const int SAMPLE_COUNT_MULTI_SCATTERING = "    + cs::utils::toString(mParams.mSampleCountMultiScattering) + ";\n" +
+      "const int STEP_SIZE_MULTI_SCATTERING = "       + cs::utils::toString(mParams.mStepSizeMultiScattering) + ";\n" +
       "const int SAMPLE_COUNT_INDIRECT_IRRADIANCE = " + cs::utils::toString(mParams.mSampleCountIndirectIrradiance) + ";\n" +
       "const vec3 SOLAR_IRRADIANCE = "                + extractVec3(WAVELENGTHS, SOLAR_IRRADIANCE, lambdas) + ";\n" +
       "const vec3 GROUND_ALBEDO = vec3("              + cs::utils::toString(mParams.mGroundAlbedo) + ");\n" +
+      "const float INDEX_OF_REFRACTION = "            + cs::utils::toString(mParams.mRefractiveIndex) + ";\n" +
       "const float SUN_ANGULAR_RADIUS = "             + cs::utils::toString(mMetadata.mSunAngularRadius) + ";\n" +
       "const float BOTTOM_RADIUS = "                  + cs::utils::toString(mParams.mMinAltitude) + ";\n" +
       "const float TOP_RADIUS = "                     + cs::utils::toString(mParams.mMaxAltitude) + ";\n" +
@@ -723,6 +750,9 @@ Preprocessor::Preprocessor(Params params)
   mTransmittanceTexture = NewTexture2d(mParams.mTransmittanceTextureWidth.get(),
       mParams.mTransmittanceTextureHeight.get(), GL_RGB32F, GL_RGB, GL_FLOAT);
 
+  mThetaDeviationTexture = NewTexture2d(mParams.mTransmittanceTextureWidth.get(),
+      mParams.mTransmittanceTextureHeight.get(), GL_RGB32F, GL_RGB, GL_FLOAT);
+
   mMultipleScatteringTexture = NewTexture3d(mScatteringTextureWidth, mScatteringTextureHeight,
       mScatteringTextureDepth, GL_RGB32F, GL_RGB, GL_FLOAT);
 
@@ -782,6 +812,7 @@ Preprocessor::~Preprocessor() {
   glDeleteVertexArrays(1, &mFullScreenQuadVAO);
   glDeleteTextures(1, &mPhaseTexture);
   glDeleteTextures(1, &mTransmittanceTexture);
+  glDeleteTextures(1, &mThetaDeviationTexture);
   glDeleteTextures(1, &mMultipleScatteringTexture);
   glDeleteTextures(1, &mSingleAerosolsScatteringTexture);
   glDeleteTextures(1, &mIrradianceTexture);
@@ -872,13 +903,22 @@ void Preprocessor::run(unsigned int numScatteringOrders) {
           numScatteringOrders);
     }
 
-    // After the above iterations, the transmittance texture contains the transmittance for the 3
-    // wavelengths used at the last iteration. But we want the transmittance at kLambdaR, kLambdaG,
-    // kLambdaB instead, so we must recompute it here for these 3 wavelengths:
+    // After the above iterations, the transmittance texture and the theta-deviation texture contain
+    // data for the 3 wavelengths used at the last iteration. But we want data at kLambdaR,
+    // kLambdaG, kLambdaB instead, so we must recompute it here for these 3 wavelengths:
     std::string header = mGlslHeaderFactory({kLambdaR, kLambdaG, kLambdaB});
     Program     computeTransmittance(kVertexShader, header + kComputeTransmittanceShader);
     glFramebufferTexture(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, mTransmittanceTexture, 0);
-    glDrawBuffer(GL_COLOR_ATTACHMENT0);
+
+    if (mParams.mRefraction.get()) {
+      glFramebufferTexture(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, mThetaDeviationTexture, 0);
+
+      const GLuint kDrawBuffers[4] = {GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1};
+      glDrawBuffers(2, kDrawBuffers);
+    } else {
+      glDrawBuffer(GL_COLOR_ATTACHMENT0);
+    }
+
     glViewport(
         0, 0, mParams.mTransmittanceTextureWidth.get(), mParams.mTransmittanceTextureHeight.get());
     glScissor(
@@ -915,8 +955,29 @@ void Preprocessor::run(unsigned int numScatteringOrders) {
 void Preprocessor::save(std::string const& directory) {
   std::cout << "Saving precomputed atmosphere to disk..." << std::endl;
 
-  // Save the precomputed textures to disk. We need to store mMultipleScatteringTexture,
-  // mSingleAerosolsScatteringTexture, mPhaseTexture, mTransmittanceTexture, and mIrradianceTexture.
+  // For debugging purposes, we print the maximum ray deviation in degrees.
+  std::vector<float> pixels(
+      mParams.mTransmittanceTextureWidth.get() * mParams.mTransmittanceTextureHeight.get() * 3);
+  glBindTexture(GL_TEXTURE_2D, mThetaDeviationTexture);
+  glGetTexImage(GL_TEXTURE_2D, 0, GL_RGB, GL_FLOAT, pixels.data());
+  glBindTexture(GL_TEXTURE_2D, 0);
+
+  float maxThetaDeviation = 0.F;
+  for (int x = 0; x < mParams.mTransmittanceTextureWidth.get(); ++x) {
+    for (int y = 0; y < mParams.mTransmittanceTextureHeight.get(); ++y) {
+      int i = 3 * (y * mParams.mTransmittanceTextureWidth.get() + x);
+
+      float thetaDeviation = pixels[i];
+      float contactRadius  = pixels[i + 1];
+
+      if (contactRadius > 0.F) {
+        maxThetaDeviation = std::max(maxThetaDeviation, thetaDeviation);
+      }
+    }
+  }
+
+  std::cout << "Maximum ray deviation: " << maxThetaDeviation * 180.F / glm::pi<float>()
+            << " degrees." << std::endl;
 
   auto write2D = [](std::string const& path, GLuint texture, int width, int height) {
     std::vector<float> data(width * height * 3);
@@ -977,6 +1038,11 @@ void Preprocessor::save(std::string const& directory) {
   write3D(directory + "/single_aerosols_scattering.tif", mSingleAerosolsScatteringTexture,
       mScatteringTextureWidth, mScatteringTextureHeight, mScatteringTextureDepth);
 
+  if (mParams.mRefraction.get()) {
+    write2D(directory + "/theta_deviation.tif", mThetaDeviationTexture,
+        mParams.mTransmittanceTextureWidth.get(), mParams.mTransmittanceTextureHeight.get());
+  }
+
   std::ofstream  out(directory + "/metadata.json");
   nlohmann::json data = mMetadata;
   out << std::setw(2) << data;
@@ -1069,7 +1135,16 @@ void Preprocessor::precompute(GLuint fbo, GLuint deltaIrradianceTexture,
 
   // 1. Compute the transmittance, and store it in mTransmittanceTexture.
   glFramebufferTexture(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, mTransmittanceTexture, 0);
-  glDrawBuffer(GL_COLOR_ATTACHMENT0);
+
+  if (mParams.mRefraction.get()) {
+    glFramebufferTexture(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, mThetaDeviationTexture, 0);
+
+    const GLuint kDrawBuffers[4] = {GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1};
+    glDrawBuffers(2, kDrawBuffers);
+  } else {
+    glDrawBuffer(GL_COLOR_ATTACHMENT0);
+  }
+
   glViewport(
       0, 0, mParams.mTransmittanceTextureWidth.get(), mParams.mTransmittanceTextureHeight.get());
   glScissor(