paysages : Restored pseudo aerial perspective for preetham atmosphere model.

git-svn-id: https://subversion.assembla.com/svn/thunderk/paysages@470 b1fd45b6-86a6-48da-8261-f70d1f35bdcc

TODO

@@ -10,7 +10,6 @@ Technology Preview 2 :
    => Apply model to atmosphere (aerial perspective)
    => Find a proper model for night sky (maybe Shirley)
 - Clouds should keep distance to ground.
-- Restore aerial perspective.
 - Implement Bruneton's scattering model.
 - Add clouds to explorer with 3d textures.
 - Start using OpenCL to optimize rendering.

gui_qt/formatmosphere.cpp

@@ -6,7 +6,7 @@
 #include <QSlider>
 #include <math.h>
 
-#include "../lib_paysages/atmosphere/atmosphere.h"
+#include "../lib_paysages/atmosphere/public.h"
 #include "../lib_paysages/scenery.h"
 #include "../lib_paysages/renderer.h"
 
@@ -20,7 +20,7 @@ public:
         BasePreview(parent)
     {
         _renderer = rendererCreate();
-        
+
         configScaling(0.5, 5.0, 0.5, 2.5);
     }
 protected:
@@ -57,7 +57,7 @@ public:
         BasePreview(parent)
     {
         _renderer = rendererCreate();
-        
+
         configScaling(0.5, 5.0, 0.5, 2.5);
     }
 protected:
@@ -92,7 +92,7 @@ FormAtmosphere::FormAtmosphere(QWidget *parent):
     BaseForm(parent)
 {
     BaseInput* input;
-    
+
     _definition = (AtmosphereDefinition*)AtmosphereDefinitionClass.create();
 
     previewWest = new PreviewSkyWest(this);
@@ -107,7 +107,7 @@ FormAtmosphere::FormAtmosphere(QWidget *parent):
     addInputDouble(tr("Sun halo radius"), &_definition->sun_halo_size, 0.0, 0.4, 0.004, 0.04);
     addInputCurve(tr("Sun halo profile"), _definition->sun_halo_profile, 0.0, 1.0, 0.0, 1.0, tr("Distance to center of the sun"), tr("Light influence (halo opacity)"));
     addInputDouble(tr("Influence of skydome on lighting"), &_definition->dome_lighting, 0.0, 2.0, 0.01, 0.1);
-    input = addInputDouble(tr("Humidity"), &_definition->humidity, 1.8, 6.0, 0.05, 0.5);
+    addInputDouble(tr("Humidity"), &_definition->humidity, 0.0, 1.0, 0.01, 0.1);
 
     revertConfig();
 }

gui_qt/formrender.cpp

@@ -33,9 +33,9 @@ public:
         _renderer.customData[1] = &_textures;
         _renderer.customData[2] = &_lighting;
         _renderer.customData[3] = &_water;
-        
+
         addOsd(QString("geolocation"));
-        
+
         configScaling(0.5, 200.0, 3.0, 50.0);
         configScrolling(-1000.0, 1000.0, 0.0, -1000.0, 1000.0, 0.0);
     }
@@ -45,7 +45,7 @@ protected:
         Vector3 down = {0.0, -1.0, 0.0};
         Vector3 location;
         double height = terrainGetHeight(&_terrain, x, y);
-        
+
         if (height < _water.height)
         {
             location.x = x;
@@ -64,9 +64,10 @@ protected:
         sceneryGetLighting(&_lighting);
         sceneryGetTextures(&_textures);
         sceneryGetWater(&_water);
-        
+
         sceneryGetAtmosphere(_atmosphere);
         AtmosphereRendererClass.bind(_renderer.atmosphere, _atmosphere);
+        _renderer.atmosphere->applyAerialPerspective = _applyAerialPerspective;
     }
 private:
     Renderer _renderer;
@@ -95,6 +96,11 @@ private:
     {
         lightingGetStatus((LightingDefinition*)renderer->customData[2], renderer, location, status);
     }
+
+    static Color _applyAerialPerspective(Renderer* renderer, Vector3 location, Color base)
+    {
+        return base;
+    }
 };
 
 /**************** Form ****************/
@@ -108,14 +114,14 @@ FormRender::FormRender(QWidget *parent) :
     _params.height = 600;
     _params.antialias = 1;
     _camera = cameraCreateDefinition();
-    
+
     _renderer_inited = false;
-    
+
     disablePreviewsUpdate();
 
     _preview_landscape = new PreviewRenderLandscape(this);
     addPreview(_preview_landscape, QString(tr("Top-down preview")));
-    
+
     addInput(new InputCamera(this, tr("Camera"), &_camera));
     addInputInt(tr("Quality"), &_params.quality, 1, 10, 1, 1);
     addInputInt(tr("Image width"), &_params.width, 100, 2000, 10, 100);
@@ -143,7 +149,7 @@ FormRender::~FormRender()
 void FormRender::savePack(PackStream* stream)
 {
     BaseForm::savePack(stream);
-    
+
     packWriteInt(stream, &_params.width);
     packWriteInt(stream, &_params.height);
     packWriteInt(stream, &_params.antialias);
@@ -153,12 +159,12 @@ void FormRender::savePack(PackStream* stream)
 void FormRender::loadPack(PackStream* stream)
 {
     BaseForm::loadPack(stream);
-    
+
     packReadInt(stream, &_params.width);
     packReadInt(stream, &_params.height);
     packReadInt(stream, &_params.antialias);
     packReadInt(stream, &_params.quality);
-    
+
     revertConfig();
 }
 
@@ -188,7 +194,7 @@ void FormRender::startQuickRender()
     }
     _renderer = sceneryCreateStandardRenderer();
     _renderer_inited = true;
-    
+
     DialogRender* dialog = new DialogRender(this, &_renderer);
     RenderParams params = {400, 300, 1, 3};
     dialog->startRender(params);
@@ -204,7 +210,7 @@ void FormRender::startRender()
     }
     _renderer = sceneryCreateStandardRenderer();
     _renderer_inited = true;
-    
+
     DialogRender* dialog = new DialogRender(this, &_renderer);
     dialog->startRender(_params);
 

lib_paysages/atmosphere/atmosphere.c

@@ -1,4 +1,5 @@
-#include "atmosphere.h"
+#include "public.h"
+#include "private.h"
 
 #include <math.h>
 #include <stdlib.h>
@@ -43,8 +44,8 @@ static AtmosphereDefinition* _createDefinition()
     curveQuickAddPoint(result->sun_halo_profile, 0.1, 0.2);
     curveQuickAddPoint(result->sun_halo_profile, 1.0, 0.0);
     result->dome_lighting = 0.6;
-    result->humidity = 2.0;
+    result->humidity = 0.1;
-    
+
     _validateDefinition(result);
 
     return result;
@@ -65,9 +66,9 @@ static void _copyDefinition(AtmosphereDefinition* source, AtmosphereDefinition*
     destination->sun_halo_size = source->sun_halo_size;
     destination->dome_lighting = source->dome_lighting;
     destination->humidity = source->humidity;
-    
+
     curveCopy(source->sun_halo_profile, destination->sun_halo_profile);
-    
+
     _validateDefinition(destination);
 }
 
@@ -107,13 +108,10 @@ StandardDefinition AtmosphereDefinitionClass = {
 };
 
 /******************** Binding ********************/
-extern Color brunetonGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3 direction, Vector3 sun_position);
+static Color _fakeApplyAerialPerspective(Renderer* renderer, Vector3 location, Color base)
-extern Color preethamGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3 direction, Vector3 sun_position);
-
-static Color _fakeApplyAerialPerspective(Renderer* renderer, Vector3 direction, Color base)
 {
     UNUSED(renderer);
-    UNUSED(direction);
+    UNUSED(location);
     return base;
 }
 static Color _fakeGetSkyColor(Renderer* renderer, Vector3 direction)
@@ -136,9 +134,9 @@ static Color _getSkyColor(Renderer* renderer, Vector3 direction)
     double dist;
     Vector3 sun_position;
     Color sky_color, sun_color;
-    
+
     definition = renderer->atmosphere->definition;
-    
+
     sun_position = renderer->atmosphere->getSunDirection(renderer);
     direction = v3Normalize(direction);
     dist = v3Norm(v3Sub(direction, sun_position));
@@ -155,7 +153,7 @@ static Color _getSkyColor(Renderer* renderer, Vector3 direction)
     default:
         sky_color = COLOR_BLUE;
     }
-    
+
     /* Get sun halo */
     if (dist < definition->sun_radius + definition->sun_halo_size)
     {
@@ -207,7 +205,7 @@ static int _getSkydomeLights(Renderer* renderer, LightDefinition* lights, int ma
     double sun_angle;
     Vector3 sun_direction;
     int nblights = 0;
-    
+
     mutexAcquire(cache->lock);
     if (cache->nblights < 0)
     {
@@ -219,7 +217,7 @@ static int _getSkydomeLights(Renderer* renderer, LightDefinition* lights, int ma
         sun_direction.z = 0.0;
 
         /* TODO Moon light */
-        
+
         if (max_lights > MAX_SKYDOME_LIGHTS)
         {
             max_lights = MAX_SKYDOME_LIGHTS;
@@ -279,7 +277,7 @@ static int _getSkydomeLights(Renderer* renderer, LightDefinition* lights, int ma
         cache->nblights = nblights;
     }
     mutexRelease(cache->lock);
-    
+
     memcpy(lights, cache->lights, sizeof(LightDefinition) * cache->nblights);
     return cache->nblights;
 }
@@ -289,20 +287,20 @@ static AtmosphereRenderer* _createRenderer()
 {
     AtmosphereRenderer* result;
     AtmosphereRendererCache* cache;
-    
+
     result = malloc(sizeof(AtmosphereRenderer));
     result->definition = AtmosphereDefinitionClass.create();
-    
+
     result->getSunDirection = _getSunDirection;
     result->applyAerialPerspective = _fakeApplyAerialPerspective;
     result->getSkydomeLights = _fakeGetSkydomeLights;
     result->getSkyColor = _fakeGetSkyColor;
-    
+
     cache = malloc(sizeof(AtmosphereRendererCache));
     cache->lock = mutexCreate();
     cache->nblights = -1;
     result->_internal_data = cache;
-    
+
     return result;
 }
 
@@ -311,7 +309,7 @@ static void _deleteRenderer(AtmosphereRenderer* renderer)
     AtmosphereRendererCache* cache = (AtmosphereRendererCache*)renderer->_internal_data;
     mutexDestroy(cache->lock);
     free(cache);
-    
+
     AtmosphereDefinitionClass.destroy(renderer->definition);
     free(renderer);
 }
@@ -319,12 +317,21 @@ static void _deleteRenderer(AtmosphereRenderer* renderer)
 static void _bindRenderer(AtmosphereRenderer* renderer, AtmosphereDefinition* definition)
 {
     AtmosphereRendererCache* cache = (AtmosphereRendererCache*)renderer->_internal_data;
-    
+
     AtmosphereDefinitionClass.copy(definition, renderer->definition);
-    
+
     renderer->getSkyColor = _getSkyColor;
     renderer->getSkydomeLights = _getSkydomeLights;
-    
+
+    switch (definition->model)
+    {
+        case ATMOSPHERE_MODEL_PREETHAM:
+            renderer->applyAerialPerspective = preethamApplyAerialPerspective;
+            break;
+        default:
+            renderer->applyAerialPerspective = _fakeApplyAerialPerspective;
+    }
+
     mutexAcquire(cache->lock);
     cache->nblights = -1;
     mutexRelease(cache->lock);
@@ -341,7 +348,7 @@ static Color _postProcessFragment(Renderer* renderer, Vector3 location, void* da
 {
     Vector3 direction;
     Color result;
-    
+
     UNUSED(data);
 
     direction = v3Sub(location, renderer->camera_location);

lib_paysages/atmosphere/bruneton.c

@@ -1,6 +1,356 @@
-#include "atmosphere.h"
+#include "public.h"
 
+/*
+ * Atmospheric scattering, based on E. Bruneton and F.Neyret work.
+ * http://evasion.inrialpes.fr/~Eric.Bruneton/
+ */
+
+#if 0
+
+#include <math.h>
+
+#define TRANSMITTANCE_NON_LINEAR
+#define INSCATTER_NON_LINEAR
+#define FIX
+
+static const double Rg = 6360.0;
+static const double Rt = 6420.0;
+static const double RL = 6421.0;
+static const double exposure = 0.4;
+static const float ISun = 100.0;
+
+// Rayleigh
+static const double HR = 8.0;
+static const Color betaR = {5.8e-3, 1.35e-2, 3.31e-2, 1.0};
+
+// Mie
+// DEFAULT
+static const double HM = 1.2;
+/*static const vec3 betaMSca = vec3(4e-3);
+static const vec3 betaMEx = betaMSca / 0.9;*/
+static const double mieG = 0.8;
+// CLEAR SKY
+/*static const float HM = 1.2;
+static const vec3 betaMSca = vec3(20e-3);
+static const vec3 betaMEx = betaMSca / 0.9;
+static const float mieG = 0.76;*/
+// PARTLY CLOUDY
+/*static const float HM = 3.0;
+static const vec3 betaMSca = vec3(3e-3);
+static const vec3 betaMEx = betaMSca / 0.9;
+static const float mieG = 0.65;*/
+
+#define step(_a_,_b_) ((_a_) < (_b_) ? 0 : 1)
+#define max(_a_,_b_) ((_a_) < (_b_) ? (_a_) : (_b_))
+#define sign(_a_) ((_a_) < 0.0 ? -1.0 : ((_a_) > 0.0 ? 1.0 : 0.0))
+
+/* Rayleigh phase function */
+static double _phaseFunctionR(double mu)
+{
+    return (3.0 / (16.0 * M_PI)) * (1.0 + mu * mu);
+}
+
+/* Mie phase function */
+static double _phaseFunctionM(double mu)
+{
+    return 1.5 * 1.0 / (4.0 * M_PI) * (1.0 - mieG * mieG) * pow(1.0 + (mieG * mieG) - 2.0 * mieG * mu, -3.0 / 2.0) * (1.0 + mu * mu) / (2.0 + mieG * mieG);
+}
+
+/* approximated single Mie scattering (cf. approximate Cm in paragraph "Angular precision") */
+static Color _getMie(Color rayMie)
+{
+    double value = rayMie.a / max(rayMie.r, 1e-4) * (betaR.r / betaR); /* TODO divide by a vector ? */
+    rayMie.r *= value;
+    rayMie.g *= value;
+    rayMie.b *= value;
+    rayMie.a = 1.0;
+    return rayMie;
+}
+
+/* optical depth for ray (r,mu) of length d, using analytic formula
+   (mu=cos(view zenith angle)), intersections with ground ignored
+   H=height scale of exponential density function */
+static double _opticalDepth(double H, double r, double mu, double d)
+{
+    double a = sqrt((0.5 / H) * r);
+    double ax = a * (mu);
+    double ay = a * (mu + d / r);
+    double axs = sign(ax);
+    double ays = sign(ay);
+    double axq = ax * ax;
+    double ayq = ay * ay;
+    double x = ays > axs ? exp(axq) : 0.0;
+    double yx = axs / (2.3193 * fabs(ax) + sqrt(1.52 * axq + 4.0));
+    double yy = ays / (2.3193 * fabs(ay) + sqrt(1.52 * ayq + 4.0)) *  exp(-d / H * (d / (2.0 * r) + mu));
+    return sqrt((6.2831 * H) * r) * exp((Rg - r) / H) * (x + yx - yy);
+}
+
+static Color _texture4D(sampler3D table, float r, float mu, float muS, float nu)
+{
+    float H = sqrt(Rt * Rt - Rg * Rg);
+    float rho = sqrt(r * r - Rg * Rg);
+#ifdef INSCATTER_NON_LINEAR
+    float rmu = r * mu;
+    float delta = rmu * rmu - r * r + Rg * Rg;
+    vec4 cst = rmu < 0.0 && delta > 0.0 ? vec4(1.0, 0.0, 0.0, 0.5 - 0.5 / float(RES_MU)) : vec4(-1.0, H * H, H, 0.5 + 0.5 / float(RES_MU));
+	float uR = 0.5 / float(RES_R) + rho / H * (1.0 - 1.0 / float(RES_R));
+    float uMu = cst.w + (rmu * cst.x + sqrt(delta + cst.y)) / (rho + cst.z) * (0.5 - 1.0 / float(RES_MU));
+    // paper formula
+    //float uMuS = 0.5 / float(RES_MU_S) + max((1.0 - exp(-3.0 * muS - 0.6)) / (1.0 - exp(-3.6)), 0.0) * (1.0 - 1.0 / float(RES_MU_S));
+    // better formula
+    float uMuS = 0.5 / float(RES_MU_S) + (atan(max(muS, -0.1975) * tan(1.26 * 1.1)) / 1.1 + (1.0 - 0.26)) * 0.5 * (1.0 - 1.0 / float(RES_MU_S));
+#else
+    float uR = 0.5 / float(RES_R) + rho / H * (1.0 - 1.0 / float(RES_R));
+    float uMu = 0.5 / float(RES_MU) + (mu + 1.0) / 2.0 * (1.0 - 1.0 / float(RES_MU));
+    float uMuS = 0.5 / float(RES_MU_S) + max(muS + 0.2, 0.0) / 1.2 * (1.0 - 1.0 / float(RES_MU_S));
+#endif
+    float lerp = (nu + 1.0) / 2.0 * (float(RES_NU) - 1.0);
+    float uNu = floor(lerp);
+    lerp = lerp - uNu;
+    return texture3D(table, vec3((uNu + uMuS) / float(RES_NU), uMu, uR)) * (1.0 - lerp) +
+           texture3D(table, vec3((uNu + uMuS + 1.0) / float(RES_NU), uMu, uR)) * lerp;
+}
+
+/* transmittance(=transparency) of atmosphere for ray (r,mu) of length d
+   (mu=cos(view zenith angle)), intersections with ground ignored
+   uses analytic formula instead of transmittance texture */
+static Vector3 _analyticTransmittance(double r, double mu, double d)
+{
+    return exp(-betaR * _opticalDepth(HR, r, mu, d) - betaMEx * _opticalDepth(HM, r, mu, d));
+}
+
+/* inscattered light along ray x+tv, when sun in direction s (=S[L]-T(x,x0)S[L]|x0) */
+static Color _inscatter(Vector3* x, double* t, Vector3 v, Vector3 s, double* _r, double* _mu, Color* attenuation)
+{
+    Color result;
+    double r = v3Norm(x);
+    double mu = v3Dot(x, v) / r;
+    double d = -r * mu - sqrt(r * r * (mu * mu - 1.0) + Rt * Rt);
+    if (d > 0.0)
+    { // if x in space and ray intersects atmosphere
+        // move x to nearest intersection of ray with top atmosphere boundary
+        x += d * v;
+        t -= d;
+        mu = (r * mu + d) / Rt;
+        r = Rt;
+    }
+    if (r <= Rt)
+    { // if ray intersects atmosphere
+        double nu = v3Dot(v, s);
+        double muS = v3Dot(x, s) / r;
+        double phaseR = _phaseFunctionR(nu);
+        double phaseM = _phaseFunctionM(nu);
+        vec4 inscatter = max(texture4D(inscatterSampler, r, mu, muS, nu), 0.0);
+        if (t > 0.0)
+        {
+            Vector3 x0 = v3Add(x, v3Scale(v, t));
+            double r0 = v3Norm(x0);
+            double rMu0 = v3Dot(x0, v);
+            double mu0 = rMu0 / r0;
+            double muS0 = v3Dot(x0, s) / r0;
+#ifdef FIX
+            // avoids imprecision problems in transmittance computations based on textures
+            *attenuation = _analyticTransmittance(r, mu, t);
+#else
+            *attenuation = _transmittance(r, mu, v, x0);
+#endif
+            if (r0 > Rg + 0.01)
+            {
+                // computes S[L]-T(x,x0)S[L]|x0
+                inscatter = max(inscatter - attenuation.rgbr * texture4D(inscatterSampler, r0, mu0, muS0, nu), 0.0);
+#ifdef FIX
+                // avoids imprecision problems near horizon by interpolating between two points above and below horizon
+                const float EPS = 0.004;
+                float muHoriz = -sqrt(1.0 - (Rg / r) * (Rg / r));
+                if (abs(mu - muHoriz) < EPS)
+                {
+                    float a = ((mu - muHoriz) + EPS) / (2.0 * EPS);
+
+                    mu = muHoriz - EPS;
+                    r0 = sqrt(r * r + t * t + 2.0 * r * t * mu);
+                    mu0 = (r * mu + t) / r0;
+                    vec4 inScatter0 = texture4D(inscatterSampler, r, mu, muS, nu);
+                    vec4 inScatter1 = texture4D(inscatterSampler, r0, mu0, muS0, nu);
+                    vec4 inScatterA = max(inScatter0 - attenuation.rgbr * inScatter1, 0.0);
+
+                    mu = muHoriz + EPS;
+                    r0 = sqrt(r * r + t * t + 2.0 * r * t * mu);
+                    mu0 = (r * mu + t) / r0;
+                    inScatter0 = texture4D(inscatterSampler, r, mu, muS, nu);
+                    inScatter1 = texture4D(inscatterSampler, r0, mu0, muS0, nu);
+                    vec4 inScatterB = max(inScatter0 - attenuation.rgbr * inScatter1, 0.0);
+
+                    inscatter = mix(inScatterA, inScatterB, a);
+                }
+#endif
+            }
+        }
+#ifdef FIX
+        // avoids imprecision problems in Mie scattering when sun is below horizon
+        inscatter.w *= smoothstep(0.00, 0.02, muS);
+#endif
+        result = max(inscatter.rgb * phaseR + getMie(inscatter) * phaseM, 0.0);
+    }
+    else
+    { // x in space and ray looking in space
+        result = COLOR_BLACK;
+    }
+
+    *_r = r;
+    *_mu = mu;
+    result.r *= ISun;
+    result.g *= ISun;
+    result.b *= ISun;
+    result.a = 1.0;
+    return result;
+}
+
+/*ground radiance at end of ray x+tv, when sun in direction s
+ *attenuated bewteen ground and viewer (=R[L0]+R[L*]) */
+
+/*static Color _groundColor(Vector3 x, double t, Vector3 v, Vector3 s, double r, double mu, Color attenuation)
+{
+    Color result;
+    if (t > 0.0)
+    { // if ray hits ground surface
+
+        // ground reflectance at end of ray, x0
+        Vector3 x0 = v3Add(x, v3Scale(v, t));
+        float r0 = v3Norm(x0);
+        Vector3 n = v3Scale(x0, 1.0 / r0);
+        vec2 coords = vec2(atan(n.y, n.x), acos(n.z)) * vec2(0.5, 1.0) / M_PI + vec2(0.5, 0.0);
+        Color reflectance;
+        if (r0 > Rg + 0.01)
+        {
+            reflectance = vec4(0.4, 0.4, 0.4, 0.0);
+        }
+        else
+        {
+            reflectance = texture2D(reflectanceSampler, coords) * vec4(0.2, 0.2, 0.2, 1.0);
+        }
+
+        // direct sun light (radiance) reaching x0
+        float muS = v3Dot(n, s);
+        Color sunLight = _transmittanceWithShadow(r0, muS);
+
+        // precomputed sky light (irradiance) (=E[L*]) at x0
+        Color groundSkyLight = irradiance(irradianceSampler, r0, muS);
+
+        // light reflected at x0 (=(R[L0]+R[L*])/T(x,x0))
+        Color groundColor = reflectance.rgb * (max(muS, 0.0) * sunLight + groundSkyLight) * ISun / M_PI;
+
+        // water specular color due to sunLight
+        if (reflectance.w > 0.0)
+        {
+            vec3 h = normalize(s - v);
+            float fresnel = 0.02 + 0.98 * pow(1.0 - dot(-v, h), 5.0);
+            float waterBrdf = fresnel * pow(max(dot(h, n), 0.0), 150.0);
+            groundColor += reflectance.w * max(waterBrdf, 0.0) * sunLight * ISun;
+        }
+
+        result = attenuation * groundColor; //=R[L0]+R[L*]
+    }
+    else
+    { // ray looking at the sky
+        return COLOR_BLACK;
+    }
+    return result;
+}*/
+
+static inline void _getTransmittanceUV(double r, double mu, double* u, double* v)
+{
+#ifdef TRANSMITTANCE_NON_LINEAR
+    *v = sqrt((r - Rg) / (Rt - Rg));
+    *u = atan((mu + 0.15) / (1.0 + 0.15) * tan(1.5)) / 1.5;
+#else
+    *v = (r - Rg) / (Rt - Rg);
+    *u = (mu + 0.15) / (1.0 + 0.15);
+#endif
+}
+
+/* transmittance(=transparency) of atmosphere for infinite ray (r,mu)
+   (mu=cos(view zenith angle)), intersections with ground ignored */
+static Color _transmittance(double r, double mu)
+{
+    double u, v;
+    _getTransmittanceUV(r, mu, u, v);
+    return texture2D(transmittanceSampler, uv).rgb;
+}
+
+/* transmittance(=transparency) of atmosphere for infinite ray (r,mu)
+   (mu=cos(view zenith angle)), or zero if ray intersects ground */
+static Color _transmittanceWithShadow(double r, double mu)
+{
+    return mu < -sqrt(1.0 - (Rg / r) * (Rg / r)) ? COLOR_BLACK : _transmittance(r, mu);
+}
+
+/* direct sun light for ray x+tv, when sun in direction s (=L0) */
+static Color _sunColor(Vector3 x, double t, Vector3 v, Vector3 s, double r, double mu)
+{
+    if (t > 0.0)
+    {
+        return COLOR_BLACK;
+    }
+    else
+    {
+        Color transmittance = r <= Rt ? _transmittanceWithShadow(r, mu) : COLOR_WHITE; // T(x,xo)
+        double isun = step(cos(M_PI / 180.0), v3Dot(v, s)) * ISun; // Lsun
+        return transmittance * isun; // Eq (9)
+    }
+}
+
+static Color _hdr(Color c1, Color c2, Color c3)
+{
+    Color L = {c1.r + c2.r + c3.r, c1.g + c2.g + c3.g, c1.b + c2.b + c3.b, 1.0};
+
+    L = L * exposure;
+    L.r = L.r < 1.413 ? pow(L.r * 0.38317, 1.0 / 2.2) : 1.0 - exp(-L.r);
+    L.g = L.g < 1.413 ? pow(L.g * 0.38317, 1.0 / 2.2) : 1.0 - exp(-L.g);
+    L.b = L.b < 1.413 ? pow(L.b * 0.38317, 1.0 / 2.2) : 1.0 - exp(-L.b);
+
+    return L;
+}
+
+Color brunetonGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3 direction, Vector3 sun_position)
+{
+    Vector3 x = {0.0, Rg + eye.y, 0.0};
+    Vector3 v = v3Normalize(direction);
+    Vector3 s = v3Normalize(v3Sub(sun_position, eye));
+
+    double r = v3Norm(x);
+    double mu = v3Dot(x, v) / r;
+    double t = -r * mu - sqrt(r * r * (mu * mu - 1.0) + Rg * Rg);
+
+    Vector3 g = {0.0, 0.0, Rg + 10.0};
+    g = v3Sub(x, g);
+    double a = v.x * v.x + v.y * v.y - v.z * v.z;
+    double b = 2.0 * (g.x * v.x + g.y * v.y - g.z * v.z);
+    double c = g.x * g.x + g.y * g.y - g.z * g.z;
+    double d = -(b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
+    int cone = d > 0.0 && fabs(x.z + d * v.z - Rg) <= 10.0;
+
+    if (t > 0.0)
+    {
+        if (cone && d < t)
+        {
+            t = d;
+        }
+    }
+    else if (cone)
+    {
+        t = d;
+    }
+
+    Vector3 attenuation;
+    Color inscatterColor = _inscatter(x, t, v, s, &r, &mu, &attenuation); //S[L]-T(x,xs)S[l]|xs
+    /*Color groundColor = _groundColor(x, t, v, s, r, mu, attenuation); //R[L0]+R[L*]*/
+    Color groundColor = COLOR_BLACK;
+    Color sunColor = _sunColor(x, t, v, s, r, mu); //L0
+    return _hdr(sunColor, groundColor, inscatterColor); // Eq (16)
+}
+#else
 Color brunetonGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3 direction, Vector3 sun_position)
 {
     return COLOR_BLACK;
 }
+#endif

lib_paysages/atmosphere/preetham.c

@@ -1,9 +1,11 @@
-#include "atmosphere.h"
+#include "public.h"
+#include "private.h"
 
 #include <math.h>
+#include <assert.h>
 #include "../renderer.h"
 
-static inline double _angleBetween(double thetav, double phiv, double theta, double phi) 
+static inline double _angleBetween(double thetav, double phiv, double theta, double phi)
 {
     double cospsi = sin(thetav) * sin(theta) * cos(phi-phiv) + cos(thetav) * cos(theta);
     if (cospsi > 1.0) return 0.0;
@@ -36,11 +38,11 @@ static inline Color _xyYToRGB(double x, double y, double Y)
 
 	result.r = r;
 	result.g = g;
-	result.b = b; 
+	result.b = b;
         result.a = 1.0;
-        
+
         colorNormalize(&result);
-        
+
         return result;
 }
 
@@ -55,16 +57,16 @@ Color preethamGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3
 {
     double theta, phi;
     double thetaSun, phiSun;
-    
+
     _directionToThetaPhi(direction, &theta, &phi);
     _directionToThetaPhi(sun_position, &thetaSun, &phiSun);
-    
+
     if (theta > M_PI / 2.0)
     {
         theta = M_PI / 2.0;
     }
     double gamma = _angleBetween(theta, phi, thetaSun, phiSun);
-    
+
     double cosTheta;
     if (theta > M_PI / 2.0)
     {
@@ -74,15 +76,15 @@ Color preethamGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3
     {
         cosTheta = cos(theta);
     }
-    
+
-    double T = definition->humidity;
+    double T = definition->humidity * 3.0 + 1.8;
     double T2 = T * T;
     double suntheta = thetaSun;
     double suntheta2 = thetaSun * thetaSun;
     double suntheta3 = thetaSun * suntheta2;
 
     double Ax = -0.01925 * T - 0.25922;
-    double Bx = -0.06651 * T + 0.00081; 
+    double Bx = -0.06651 * T + 0.00081;
     double Cx = -0.00041 * T + 0.21247;
     double Dx = -0.06409 * T - 0.89887;
     double Ex = -0.00325 * T + 0.04517;
@@ -102,16 +104,16 @@ Color preethamGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3
     double cosGamma = cos(gamma);
     cosGamma = cosGamma < 0.0 ? 0.0 : cosGamma;
     double cosSTheta = fabs(cos(thetaSun));
-    
+
     double xz = ( 0.00165 * suntheta3 - 0.00375 * suntheta2 + 0.00209 * suntheta + 0.00000) * T2 +
 		 (-0.02903 * suntheta3 + 0.06377 * suntheta2 - 0.03202 * suntheta + 0.00394) * T +
 		 ( 0.11693 * suntheta3 - 0.21196 * suntheta2 + 0.06052 * suntheta + 0.25886);
 
     double yz = ( 0.00275 * suntheta3 - 0.00610 * suntheta2 + 0.00317 * suntheta + 0.00000) * T2 +
-                (-0.04214 * suntheta3 + 0.08970 * suntheta2 - 0.04153 * suntheta + 0.00516) * T +	
+                (-0.04214 * suntheta3 + 0.08970 * suntheta2 - 0.04153 * suntheta + 0.00516) * T +
                 ( 0.15346 * suntheta3 - 0.26756 * suntheta2 + 0.06670 * suntheta + 0.26688);
-    
+
-    double X = (4.0 / 9.0 - T / 120.0) * (M_PI - 2.0 * suntheta);		
+    double X = (4.0 / 9.0 - T / 120.0) * (M_PI - 2.0 * suntheta);
     double Yz = ((4.0453 * T - 4.9710) * tan(X) - 0.2155 * T + 2.4192) * 1000.0;
 
     double val1, val2;
@@ -128,5 +130,38 @@ Color preethamGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3
     val2 = (1.0 + AY * exp(BY)) * (1.0 + CY * exp(DY * suntheta) + EY * sqrt(cosSTheta));
     double Y = Yz * val1 / val2;
 
-    return _xyYToRGB(x, y, Y);
+    Color result = _xyYToRGB(x, y, Y);
+    Color humidity_color = {0.8, 0.8, 0.8, definition->humidity * 0.8 * colorGetValue(&result)};
+    colorMask(&result, &humidity_color);
+    return result;
+}
+
+Color preethamApplyAerialPerspective(Renderer* renderer, Vector3 location, Color base)
+{
+    Vector3 ray = v3Sub(location, renderer->camera_location);
+    Vector3 direction = v3Normalize(ray);
+    double distance = v3Norm(ray);
+    AtmosphereDefinition* definition = renderer->atmosphere->definition;
+    double near = 10.0 - definition->humidity * 10.0;
+    double far = 100.0 - definition->humidity * 90.0;
+    double max = 0.85 + definition->humidity * 0.12;
+
+    assert(near < far);
+
+    if (distance < near)
+    {
+        return base;
+    }
+    else
+    {
+        Vector3 sun_position = renderer->atmosphere->getSunDirection(renderer);
+        Color sky = preethamGetSkyColor(definition, renderer->camera_location, direction, sun_position);
+        if (distance > far)
+        {
+            distance = far;
+        }
+        sky.a = max * ((distance - near) / (far - near));
+        colorMask(&base, &sky);
+        return base;
+    }
 }

lib_paysages/atmosphere/private.h

@@ -0,0 +1,9 @@
+#ifndef _PAYSAGES_ATMOSPHERE_PRIVATE_H_
+#define _PAYSAGES_ATMOSPHERE_PRIVATE_H_
+
+Color brunetonGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3 direction, Vector3 sun_position);
+
+Color preethamGetSkyColor(AtmosphereDefinition* definition, Vector3 eye, Vector3 direction, Vector3 sun_position);
+Color preethamApplyAerialPerspective(Renderer* renderer, Vector3 location, Color base);
+
+#endif

lib_paysages/atmosphere/public.h

@@ -1,5 +1,5 @@
-#ifndef _PAYSAGES_ATMOSPHERE_INTERFACE_H_
+#ifndef _PAYSAGES_ATMOSPHERE_PUBLIC_H_
-#define _PAYSAGES_ATMOSPHERE_INTERFACE_H_
+#define _PAYSAGES_ATMOSPHERE_PUBLIC_H_
 
 #include "../shared/types.h"
 #include "../euclid.h"
@@ -19,7 +19,7 @@ typedef void (*FuncObjectValidate)(void* object);
 typedef void (*FuncObjectSave)(PackStream* stream, void* object);
 typedef void (*FuncObjectLoad)(PackStream* stream, void* object);
 typedef void (*FuncObjectBind)(void* base, void* sub);
-    
+
 typedef struct {
     FuncObjectCreate create;
     FuncObjectDelete destroy;
@@ -62,12 +62,12 @@ typedef Vector3 (*FuncAtmosphereGetSunDirection)(Renderer* renderer);
 typedef struct
 {
     AtmosphereDefinition* definition;
-    
+
     FuncAtmosphereGetSkydomeLights getSkydomeLights;
     FuncAtmosphereApplyAerialPerspective applyAerialPerspective;
     FuncAtmosphereGetSkyColor getSkyColor;
     FuncAtmosphereGetSunDirection getSunDirection;
-    
+
     void* _internal_data;
 } AtmosphereRenderer;
 

lib_paysages/renderer.h

@@ -2,7 +2,7 @@
 #define _PAYSAGES_RENDERER_H_
 
 #include "shared/types.h"
-#include "atmosphere/atmosphere.h"
+#include "atmosphere/public.h"
 
 #ifdef __cplusplus
 extern "C" {
@@ -40,7 +40,7 @@ struct Renderer
     void (*alterLight)(Renderer* renderer, LightDefinition* light, Vector3 location);
     void (*getLightStatus)(Renderer* renderer, LightStatus* status, Vector3 location);
     Color (*applyLightStatus)(Renderer* renderer, LightStatus* status, Vector3 location, Vector3 normal, SurfaceMaterial material);
-    
+
     /* Autonomous sub-renderers */
     AtmosphereRenderer* atmosphere;
 

lib_paysages/scenery.h

@@ -8,7 +8,7 @@
  * a standard renderer.
  */
 
-#include "atmosphere/atmosphere.h"
+#include "atmosphere/public.h"
 #include "camera.h"
 #include "clouds.h"
 #include "lighting.h"
@@ -23,7 +23,7 @@ extern "C" {
 #endif
 
 typedef void (*SceneryCustomDataCallback)(PackStream* stream, void* data);
-    
+
 void sceneryInit();
 void sceneryQuit();