Take the Sun angular radius into account in single scattering and dir…

…ect irradiance computations. Also fix a bug in GetSunAndSkyIrradiance (altitude was not taken into account to compute the visible Sun disc fraction).

atmosphere/definitions.glsl

@@ -205,7 +205,8 @@ The atmosphere parameters are then defined by the following struct:
 struct AtmosphereParameters {
   // The solar irradiance at the top of the atmosphere.
   IrradianceSpectrum solar_irradiance;
-  // The sun's angular radius.
+  // The sun's angular radius. Warning: the implementation uses approximations
+  // that are valid only if this angle is smaller than 0.1 radians.
   Angle sun_angular_radius;
   // The distance between the planet center and the bottom of the atmosphere.
   Length bottom_radius;

atmosphere/functions.glsl

@@ -526,6 +526,43 @@ very close to the horizon, due to the finite precision and rounding errors of
 floating point operations. And also because the caller generally has more robust
 ways to know whether a ray intersects the ground or not (see below).
 
+<p>Finally, we will also need the transmittance between a point in the
+atmosphere and the Sun. The Sun is not a point light source, so this is an
+integral of the transmittance over the Sun disc. Here we consider that the
+transmittance is constant over this disc, except below the horizon, where the
+transmittance is 0. As a consequence, the transmittance to the Sun can be
+computed with <code>GetTransmittanceToTopAtmosphereBoundary</code>, times the
+fraction of the Sun disc which is above the horizon.
+
+<p>This fraction varies from 0 when the Sun zenith angle $\theta_s$ is larger
+than the horizon zenith angle $\theta_h$ plus the Sun angular radius $\alpha_s$,
+to 1 when $\theta_s$ is smaller than $\theta_h-\alpha_s$. Equivalently, it
+varies from 0 when $\mu_s=\cos\theta_s$ is smaller than
+$\cos(\theta_h+\alpha_s)\approx\cos\theta_h-\alpha_s\sin\theta_h$ to 1 when
+$\mu_s$ is larger than
+$\cos(\theta_h-\alpha_s)\approx\cos\theta_h+\alpha_s\sin\theta_h$. In between,
+the visible Sun disc fraction varies approximately like a smoothstep (this can
+be verified by plotting the area of <a
+href="https://en.wikipedia.org/wiki/Circular_segment">circular segment</a> as a
+function of its <a href="https://en.wikipedia.org/wiki/Sagitta_(geometry)"
+>sagitta</a>). Therefore, since $\sin\theta_h=r_{\mathrm{bottom}}/r$, we can
+approximate the transmittance to the Sun with the following function:
+*/
+
+DimensionlessSpectrum GetTransmittanceToSun(
+    IN(AtmosphereParameters) atmosphere,
+    IN(TransmittanceTexture) transmittance_texture,
+    Length r, Number mu_s) {
+  Number sin_theta_h = atmosphere.bottom_radius / r;
+  Number cos_theta_h = -sqrt(max(1.0 - sin_theta_h * sin_theta_h, 0.0));
+  return GetTransmittanceToTopAtmosphereBoundary(
+          atmosphere, transmittance_texture, r, mu_s) *
+      smoothstep(-sin_theta_h * atmosphere.sun_angular_radius / rad,
+                 sin_theta_h * atmosphere.sun_angular_radius / rad,
+                 mu_s - cos_theta_h);
+}
+
+/*
 <h3 id="single_scattering">Single scattering</h3>
 
 <p>The single scattered radiance is the light arriving from the Sun at some
@@ -579,8 +616,8 @@ below):
 <p>The radiance arriving at $\bp$ is the product of:
 <ul>
 <li>the solar irradiance at the top of the atmosphere,</li>
-<li>the transmittance between the top of the atmosphere and $\bq$ (i.e. the
-fraction of the light at the top of the atmosphere that reaches $\bq$),</li>
+<li>the transmittance between the Sun and $\bq$ (i.e. the fraction of the Sun
+light at the top of the atmosphere that reaches $\bq$),</li>
 <li>the Rayleigh scattering coefficient at $\bq$ (i.e. the fraction of the
 light arriving at $\bq$ which is scattered, in any direction),</li>
 <li>the Rayleigh phase function (i.e. the fraction of the scattered light at
@@ -618,22 +655,16 @@ void ComputeSingleScatteringIntegrand(
     OUT(DimensionlessSpectrum) rayleigh, OUT(DimensionlessSpectrum) mie) {
   Length r_d = ClampRadius(atmosphere, sqrt(d * d + 2.0 * r * mu * d + r * r));
   Number mu_s_d = ClampCosine((r * mu_s + d * nu) / r_d);
-
-  if (RayIntersectsGround(atmosphere, r_d, mu_s_d)) {
-    rayleigh = DimensionlessSpectrum(0.0);
-    mie = DimensionlessSpectrum(0.0);
-  } else {
-    DimensionlessSpectrum transmittance =
-        GetTransmittance(
-            atmosphere, transmittance_texture, r, mu, d,
-            ray_r_mu_intersects_ground) *
-        GetTransmittanceToTopAtmosphereBoundary(
-            atmosphere, transmittance_texture, r_d, mu_s_d);
-    rayleigh = transmittance * GetProfileDensity(
-        atmosphere.rayleigh_density, r_d - atmosphere.bottom_radius);
-    mie = transmittance * GetProfileDensity(
-        atmosphere.mie_density, r_d - atmosphere.bottom_radius);
-  }
+  DimensionlessSpectrum transmittance =
+      GetTransmittance(
+          atmosphere, transmittance_texture, r, mu, d,
+          ray_r_mu_intersects_ground) *
+      GetTransmittanceToSun(
+          atmosphere, transmittance_texture, r_d, mu_s_d);
+  rayleigh = transmittance * GetProfileDensity(
+      atmosphere.rayleigh_density, r_d - atmosphere.bottom_radius);
+  mie = transmittance * GetProfileDensity(
+      atmosphere.mie_density, r_d - atmosphere.bottom_radius);
 }
 
 /*
@@ -1393,11 +1424,14 @@ direct irradiance.
 <p>The irradiance is the integral over an hemisphere of the incident radiance,
 times a cosine factor. For the direct ground irradiance, the incident radiance
 is the Sun radiance at the top of the atmosphere, times the transmittance
-through the atmosphere. And, since this radiance is zero outside the small solid
-angle of the Sun, we can approximate the irradiance integral with the Sun
-radiance, times the Sun solid angle (yielding the solar irradiance), times the
-transmittance and the cosine factor for the Sun direction, i.e. $\mu_s$. This
-yields the following implementation:
+through the atmosphere. And, since the Sun solid angle is small, we can
+approximate the transmittance with a constant, i.e. we can move it outside the
+irradiance integral, which can be performed over (the visible fraction of) the
+Sun disc rather than the hemisphere. Then the integral becomes equivalent to the
+ambient occlusion due to a sphere, also called a view factor, which is given in
+<a href="http://webserver.dmt.upm.es/~isidoro/tc3/Radiation%20View%20factors.pdf
+">Radiative view factors</a> (page 10). For a small solid angle, these complex
+equations can be simplified as follows:
 */
 
 IrradianceSpectrum ComputeDirectIrradiance(
@@ -1407,9 +1441,17 @@ IrradianceSpectrum ComputeDirectIrradiance(
   assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
   assert(mu_s >= -1.0 && mu_s <= 1.0);
 
+  Number alpha_s = atmosphere.sun_angular_radius / rad;
+  // Approximate average of the cosine factor mu_s over the visible fraction of
+  // the Sun disc.
+  Number average_cosine_factor =
+    mu_s < -alpha_s ? 0.0 : (mu_s > alpha_s ? mu_s :
+        (mu_s + alpha_s) * (mu_s + alpha_s) / (4.0 * alpha_s));
+
   return atmosphere.solar_irradiance *
       GetTransmittanceToTopAtmosphereBoundary(
-          atmosphere, transmittance_texture, r, mu_s) * max(mu_s, 0.0);
+          atmosphere, transmittance_texture, r, mu_s) * average_cosine_factor;
+
 }
 
 /*
@@ -1817,39 +1859,12 @@ RadianceSpectrum GetSkyRadianceToPoint(
 
 <p>To render the ground we need the irradiance received on the ground after 0 or
 more bounce(s) in the atmosphere or on the ground. The direct irradiance can be
-computed with a lookup in the transmittance texture, while the indirect
-irradiance is given by a lookup in the precomputed irradiance texture (this
-texture only contains the irradiance for horizontal surfaces; we use the
-approximation defined in our
+computed with a lookup in the transmittance texture,
+via <code>GetTransmittanceToSun</code>, while the indirect irradiance is given
+by a lookup in the precomputed irradiance texture (this texture only contains
+the irradiance for horizontal surfaces; we use the approximation defined in our
 <a href="https://hal.inria.fr/inria-00288758/en">paper</a> for the other cases).
-
-<p>Note that it is useful here to take the angular size of the sun into account.
-With a punctual light source (as we assumed in all the above functions), the
-direct irradiance on a slanted surface would be discontinuous when the sun
-moves across the horizon. With an area light source this discontinuity issue
-disappears because the visible sun area decreases continously as the sun moves
-across the horizon.
-
-<p>Taking the angular size of the sun into account, without approximations, is
-quite complex because the visible sun area is restricted both by the distant
-horizon and by the local surface. Here we ignore the masking by the local
-surface, and we approximate the masking by the horizon with a
-<code>smoothstep</code>.
-
-<p>The smoothstep approximation is justified as follows. When the sun, of
-angular radius $\alpha_s$, is at an angle $\alpha$ above the horizon (we assume
-that $\alpha$ is between $-\alpha_s$ and $\alpha_s$), the fraction $f$ of its
-surface which is visible can be computed from the area of a <a
-href="https://en.wikipedia.org/wiki/Circular_segment">circular segment</a>:
-$f(\alpha)=(\theta-\sin\theta)/2\pi$, with
-$\theta=2\arccos(-\alpha/\alpha_s)$. The smoothstep approximation is justified
-by the fact that $f$, expressed as a function of the cosine of the sun zenith
-angle, $\mu_s=\sin\alpha\approx\alpha$, is quite similar to
-<code>smoothstep(-alpha_s, alpha_s, mu_s)</code>.
-
-<p>The function below returns the direct and indirect irradiances separately,
-and takes the angular size of the sun into account by using the above
-approximation:
+The function below returns the direct and indirect irradiances separately:
 */
 
 IrradianceSpectrum GetSunAndSkyIrradiance(
@@ -1867,10 +1882,7 @@ IrradianceSpectrum GetSunAndSkyIrradiance(
 
   // Direct irradiance.
   return atmosphere.solar_irradiance *
-      GetTransmittanceToTopAtmosphereBoundary(
+      GetTransmittanceToSun(
           atmosphere, transmittance_texture, r, mu_s) *
-      smoothstep(-atmosphere.sun_angular_radius / rad,
-                 atmosphere.sun_angular_radius / rad,
-                 mu_s) *
       max(dot(normal, sun_direction), 0.0);
 }

atmosphere/model.h

@@ -190,7 +190,8 @@ class Model {
     // The solar irradiance at the top of the atmosphere, in W/m^2/nm. This
     // vector must have the same size as the wavelengths parameter.
     const std::vector<double>& solar_irradiance,
-    // The sun's angular radius, in radians.
+    // The sun's angular radius, in radians. Warning: the implementation uses
+    // approximations that are valid only if this value is smaller than 0.1.
     double sun_angular_radius,
     // The distance between the planet center and the bottom of the atmosphere,
     // in m.

atmosphere/reference/definitions.h

@@ -240,7 +240,8 @@ struct DensityProfile {
 struct AtmosphereParameters {
   // The solar irradiance at the top of the atmosphere.
   IrradianceSpectrum solar_irradiance;
-  // The sun's angular radius.
+  // The sun's angular radius. Warning: the implementation uses approximations
+  // that are valid only if this angle is smaller than 0.1 radians.
   Angle sun_angular_radius;
   // The distance between the planet center and the bottom of the atmosphere.
   Length bottom_radius;