paysages : Preetham approximation for sky (WIP) - Reorganize modules.

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

TODO

@@ -6,6 +6,7 @@ Technology Preview 2 :
    => Find a proper model for night sky (maybe Shirley)
 - InputInt doesn't honor small_step.
 - Disable form fields when no layer is selected.
+- Add buttons to restore "auto" default values in tabs and dialogs.
 - Remove color gradations (replace with automatic boolean and simple colors).
 - Replace zone ranges with curves (with curve input and curve dialog).
 - Interface for textures thickness, slope_range and thickness_transparency (and correct slider ranges).
@@ -33,6 +34,7 @@ Technology Preview 2 :
 
 Technology Preview 3 :
 - Restore render progress.
+- Implement High Dynamic Range.
 - Use bicubic interpolation for antialiasing.
 - Improve large render/antialias memory consumption.
 - Add basic vegetation system ?

gui_qt/formatmosphere.cpp

@@ -66,7 +66,7 @@ FormAtmosphere::FormAtmosphere(QWidget *parent):
     addInputDouble(tr("Start distance"), &_definition.distance_near, -500.0, 500.0, 5.0, 50.0);
     addInputDouble(tr("End distance"), &_definition.distance_far, -500.0, 500.0, 5.0, 50.0);
     addInputDouble(tr("Masking power"), &_definition.full_mask, 0.0, 1.0, 0.01, 0.1);
-    addInputBoolean(tr("Lock color on haze"), &_definition.auto_lock_on_haze);
+    addInputBoolean(tr("Lock on horizon color"), &_definition.auto_lock_on_haze);
     addInputColor(tr("Color"), &_definition.color);
 
     revertConfig();

gui_qt/formsky.cpp

@@ -106,7 +106,7 @@ FormSky::FormSky(QWidget *parent):
     previewEast = new PreviewSkyEast(this);
     addPreview(previewEast, QString(tr("East preview")));
 
-    addInputEnum(tr("Color model"), (int*)&_definition.model, QStringList(tr("Custom model")) << tr("Mixed Preetham/Shirley approximation"));
+    addInputEnum(tr("Color model"), (int*)&_definition.model, QStringList(tr("Custom model")) << tr("Rayleigh/Mie scattering") << tr("Preetham/Shirley analytic model"));
     addInputDouble(tr("Day time"), &_definition.daytime, 0.0, 1.0, 0.002, 0.1);
     addInputColor(tr("Sun color"), &_definition.sun_color);
     addInputDouble(tr("Sun radius"), &_definition.sun_radius, 0.0, 0.4, 0.004, 0.04);

gui_qt/formterrain.cpp

@@ -150,7 +150,7 @@ FormTerrain::FormTerrain(QWidget *parent):
 
     addInputNoise(tr("Noise"), _definition.height_noise);
     addInputDouble(tr("Height"), &_definition.height_factor, 0.0, 20.0, 0.1, 1.0);
-    addInputDouble(tr("Scaling"), &_definition.scaling, 1.0, 20.0, 0.1, 1.0);
+    addInputDouble(tr("Scaling"), &_definition.scaling, 1.0, 50.0, 0.1, 5.0);
     addInputDouble(tr("Shadow smoothing"), &_definition.shadow_smoothing, 0.0, 0.3, 0.003, 0.03);
 
     revertConfig();

i18n/paysages_fr.ts

@@ -335,8 +335,12 @@ Maintenir Ctrl : Plus rapide</translation>
     </message>
     <message>
         <location filename="../gui_qt/formatmosphere.cpp" line="69"/>
+        <source>Lock on horizon color</source>
+        <translation type="unfinished"></translation>
+    </message>
+    <message>
         <source>Lock color on haze</source>
-        <translation>Verrouiller sur la couleur de la brume</translation>
+        <translation type="obsolete">Verrouiller sur la couleur de la brume</translation>
     </message>
     <message>
         <location filename="../gui_qt/formatmosphere.cpp" line="70"/>
@@ -608,12 +612,17 @@ Maintenir Ctrl : Plus rapide</translation>
     </message>
     <message>
         <location filename="../gui_qt/formsky.cpp" line="109"/>
-        <source>Mixed Preetham/Shirley approximation</source>
+        <source>Custom model</source>
         <translation type="unfinished"></translation>
     </message>
     <message>
         <location filename="../gui_qt/formsky.cpp" line="109"/>
-        <source>Custom model</source>
+        <source>Rayleigh/Mie scattering</source>
+        <translation type="unfinished"></translation>
+    </message>
+    <message>
+        <location filename="../gui_qt/formsky.cpp" line="109"/>
+        <source>Preetham/Shirley analytic model</source>
         <translation type="unfinished"></translation>
     </message>
     <message>

lib_paysages/atmosphere.c

@@ -77,7 +77,7 @@ void atmosphereValidateDefinition(AtmosphereDefinition* definition)
     {
         sky = skyCreateDefinition();
         sceneryGetSky(&sky);
-        definition->color = sky.model_custom.haze_color;
+        definition->color = skyGetHorizonColor(&sky);
         skyDeleteDefinition(&sky);
     }
 }

lib_paysages/preetham.c

@@ -0,0 +1,137 @@
+#include "preetham.h"
+
+#include <math.h>
+
+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;
+    if (cospsi < -1.0) return M_PI;
+    return acos(cospsi);
+}
+
+static inline Color _xyYToRGB(double x, double y, double Y)
+{
+	double fX, fY, fZ;
+        Color result;
+
+	fY = Y;
+	fX = x / y * Y;
+	fZ = ((1.0 - x - y) / y) * Y;
+
+	double r, g, b;
+	r =  3.2404 * fX - 1.5371 * fY - 0.4985 * fZ;
+	g = -0.9692 * fX + 1.8759 * fY + 0.0415 * fZ;
+	b =  0.0556 * fX - 0.2040 * fY + 1.0573 * fZ;
+
+	double expo = -(1.0 / 10000.0);
+	r = 1.0 - exp(expo * r);
+	g = 1.0 - exp(expo * g);
+	b = 1.0 - exp(expo * b);
+
+	if (r < 0.0) r = 0.0;
+	if (g < 0.0) g = 0.0;
+	if (b < 0.0) b = 0.0;
+
+	result.r = r;
+	result.g = g;
+	result.b = b; 
+        result.a = 1.0;
+        
+        colorNormalize(&result);
+        
+        return result;
+}
+
+static inline void _directionToThetaPhi(Vector3 direction, double* theta, double* phi)
+{
+    direction = v3Normalize(direction);
+    *phi = atan2(-direction.z, direction.x);
+    *theta = M_PI_2 - asin(direction.y);
+}
+
+Color preethamGetSkyColor(Vector3 viewer, Vector3 direction, Vector3 sun_direction, double turbidity)
+{
+    double theta, phi;
+    double thetaSun, phiSun;
+    
+    _directionToThetaPhi(direction, &theta, &phi);
+    _directionToThetaPhi(sun_direction, &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)
+    {
+        cosTheta = 0.0000001;
+    }
+    else
+    {
+        cosTheta = cos(theta);
+    }
+    
+    double T = turbidity;
+    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 Cx = -0.00041 * T + 0.21247;
+    double Dx = -0.06409 * T - 0.89887;
+    double Ex = -0.00325 * T + 0.04517;
+
+    double Ay = -0.01669 * T - 0.26078;
+    double By = -0.09495 * T + 0.00921;
+    double Cy = -0.00792 * T + 0.21023;
+    double Dy = -0.04405 * T - 1.65369;
+    double Ey = -0.01092 * T + 0.05291;
+
+    double AY =  0.17872 * T - 1.46303;
+    double BY = -0.35540 * T + 0.42749;
+    double CY = -0.02266 * T + 5.32505;
+    double DY =  0.12064 * T - 2.57705;
+    double EY = -0.06696 * T + 0.37027;
+
+    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.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 Yz = ((4.0453 * T - 4.9710) * tan(X) - 0.2155 * T + 2.4192) * 1000.0;
+
+    double val1, val2;
+
+    val1 = (1.0 + Ax * exp(Bx / cosTheta)) * (1.0 + Cx * exp(Dx * gamma) + Ex * sqrt(cosGamma));
+    val2 = (1.0 + Ax * exp(Bx)) * (1.0 + Cx * exp(Dx * suntheta) + Ex * sqrt(cosSTheta));
+    double x = xz * val1 / val2;
+
+    val1 = (1.0 + Ay * exp(By / cosTheta)) * (1.0 + Cy * exp(Dy * gamma) + Ey * sqrt(cosGamma));
+    val2 = (1.0 + Ay * exp(By)) * (1.0 + Cy * exp(Dy * suntheta) + Ey * sqrt(cosSTheta));
+    double y = yz * val1 / val2;
+
+    val1 = (1.0 + AY * exp(BY / cosTheta)) * (1.0 + CY * exp(DY * gamma) + EY * sqrt(cosGamma));
+    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 preethamApplyToObject(Vector3 viewer, Vector3 object_location, Vector3 sun_direction, double turbidity, Color object_color)
+{
+    /* TODO Aerial perspective */
+    return object_color;
+}

lib_paysages/preetham.h

@@ -0,0 +1,20 @@
+#ifndef _PAYSAGES_PREETHAM_H_
+#define _PAYSAGES_PREETHAM_H_
+
+/* Implementation of Preetham/Shirley light scattering */
+
+#include "color.h"
+#include "euclid.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+Color preethamGetSkyColor(Vector3 viewer, Vector3 direction, Vector3 sun_direction, double turbidity);
+Color preethamApplyToObject(Vector3 viewer, Vector3 object_location, Vector3 sun_direction, double turbidity, Color object_color);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

lib_paysages/rayleigh.c

@@ -0,0 +1,30 @@
+#include "rayleigh.h"
+
+#include <math.h>
+
+/*static double _phase(double g, double theta)
+{
+    double g2 = g * g;
+    double costheta = cos(theta);
+    return ((3.0 * (1.0 - g2)) / (2.0 * (2.0 + g2))) * ((1.0 + costheta * costheta) / exp(1.0 + g2 - 2.0 * g * costheta, 1.5));
+}
+
+static double _intOpticalDepthCallback(double h, double* h0)
+{
+    return -h / *h0;
+}
+
+static double _opticalDepth(double h0, double lambda, double k, Vector3 start, Vector3 end)
+{
+    return 4.0 * M_PI * k * toolsIntegrate(_intOpticalDepthCallback);
+}*/
+
+Color rayleighGetSkyColor(Vector3 viewer, Vector3 direction, Vector3 sun_direction)
+{
+    return COLOR_BLACK;
+}
+
+Color rayleighApplyToObject(Vector3 viewer, Vector3 object_location, Vector3 sun_direction, Color object_color)
+{
+    return COLOR_BLACK;
+}

lib_paysages/rayleigh.h

@@ -0,0 +1,20 @@
+#ifndef _PAYSAGES_RAYLEIGH_H_
+#define _PAYSAGES_RAYLEIGH_H_
+
+/* Implementation of Rayleigh/Mie atmospheric light scattering */
+
+#include "color.h"
+#include "euclid.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+Color rayleighGetSkyColor(Vector3 viewer, Vector3 direction, Vector3 sun_direction);
+Color rayleighApplyToObject(Vector3 viewer, Vector3 object_location, Vector3 sun_direction, Color object_color);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

lib_paysages/sky.c

@@ -10,6 +10,8 @@
 #include "lighting.h"
 #include "render.h"
 #include "tools.h"
+#include "preetham.h"
+#include "rayleigh.h"
 
 #define SPHERE_SIZE 1000.0
 
@@ -183,7 +185,7 @@ int skyGetLights(SkyDefinition* sky, LightDefinition* lights, int max_lights)
             lights[1].direction.x = 0.0;
             lights[1].direction.y = -1.0;
             lights[1].direction.z = 0.0;
-            lights[1].color = sky->model_custom.zenith_color;
+            lights[1].color = skyGetZenithColor(sky);
             lights[1].color.r *= 0.6;
             lights[1].color.g *= 0.6;
             lights[1].color.b *= 0.6;
@@ -198,138 +200,6 @@ int skyGetLights(SkyDefinition* sky, LightDefinition* lights, int max_lights)
     return nblights;
 }
 
-/******************************** Preetham Model ********************************/
-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;
-    if (cospsi < -1.0) return M_PI;
-    return acos(cospsi);
-}
-
-static inline Color _xyYToRGB(double x, double y, double Y)
-{
-	double fX, fY, fZ;
-        Color result;
-
-	fY = Y;
-	fX = x / y * Y;
-	fZ = ((1.0 - x - y) / y) * Y;
-
-	double r, g, b;
-	r =  3.2404 * fX - 1.5371 * fY - 0.4985 * fZ;
-	g = -0.9692 * fX + 1.8759 * fY + 0.0415 * fZ;
-	b =  0.0556 * fX - 0.2040 * fY + 1.0573 * fZ;
-
-	double expo = -(1.0 / 10000.0);
-	r = 1.0 - exp(expo * r);
-	g = 1.0 - exp(expo * g);
-	b = 1.0 - exp(expo * b);
-
-	if (r < 0.0) r = 0.0;
-	if (g < 0.0) g = 0.0;
-	if (b < 0.0) b = 0.0;
-
-	result.r = r;
-	result.g = g;
-	result.b = b; 
-        result.a = 1.0;
-        
-        colorNormalize(&result);
-        
-        return result;
-}
-
-static Color _preethamApproximate(SkyDefinition* definition, double theta, double phi)
-{
-    double thetaSun;
-    double phiSun;
-    double gamma;
-    double turbidity = definition->model_preetham.turbidity;
-
-    /* Handle angles */
-    if (theta > M_PI / 2.0)
-    {
-        theta = M_PI / 2.0;
-    }
-    if (definition->daytime <= 0.5)
-    {
-        thetaSun = M_PI - definition->daytime * 2.0 * M_PI;
-        phiSun = 0.0;
-    }
-    else
-    {
-        thetaSun = (definition->daytime - 0.5) * 2.0 * M_PI;
-        phiSun = M_PI;
-    }
-    gamma = _angleBetween(theta, phi, thetaSun, phiSun);
-    
-    double cosTheta;
-    if (theta > M_PI / 2.0)
-    {
-        cosTheta = 0.0000001;
-    }
-    else
-    {
-        cosTheta = cos(theta);
-    }
-    
-    double T = turbidity;
-    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 Cx = -0.00041 * T + 0.21247;
-    double Dx = -0.06409 * T - 0.89887;
-    double Ex = -0.00325 * T + 0.04517;
-
-    double Ay = -0.01669 * T - 0.26078;
-    double By = -0.09495 * T + 0.00921;
-    double Cy = -0.00792 * T + 0.21023;
-    double Dy = -0.04405 * T - 1.65369;
-    double Ey = -0.01092 * T + 0.05291;
-
-    double AY =  0.17872 * T - 1.46303;
-    double BY = -0.35540 * T + 0.42749;
-    double CY = -0.02266 * T + 5.32505;
-    double DY =  0.12064 * T - 2.57705;
-    double EY = -0.06696 * T + 0.37027;
-
-    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.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 Yz = ((4.0453 * T - 4.9710) * tan(X) - 0.2155 * T + 2.4192) * 1000.0;
-
-    double val1, val2;
-
-    val1 = (1.0 + Ax * exp(Bx / cosTheta)) * (1.0 + Cx * exp(Dx * gamma) + Ex * sqrt(cosGamma));
-    val2 = (1.0 + Ax * exp(Bx)) * (1.0 + Cx * exp(Dx * suntheta) + Ex * sqrt(cosSTheta));
-    double x = xz * val1 / val2;
-
-    val1 = (1.0 + Ay * exp(By / cosTheta)) * (1.0 + Cy * exp(Dy * gamma) + Ey * sqrt(cosGamma));
-    val2 = (1.0 + Ay * exp(By)) * (1.0 + Cy * exp(Dy * suntheta) + Ey * sqrt(cosSTheta));
-    double y = yz * val1 / val2;
-
-    val1 = (1.0 + AY * exp(BY / cosTheta)) * (1.0 + CY * exp(DY * gamma) + EY * sqrt(cosGamma));
-    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);
-}
-
 /******************************** Rendering ********************************/
 Color skyGetColor(SkyDefinition* definition, Renderer* renderer, Vector3 eye, Vector3 look)
 {
@@ -344,7 +214,11 @@ Color skyGetColor(SkyDefinition* definition, Renderer* renderer, Vector3 eye, Ve
 
     if (definition->model == SKY_MODEL_PREETHAM)
     {
-        sky_color = _preethamApproximate(definition, M_PI/2.0 - asin(look.y), atan2(-look.z, look.x));
+        sky_color = preethamGetSkyColor(eye, look, sun_position, definition->model_preetham.turbidity);
+    }
+    else if (definition->model == SKY_MODEL_RAYLEIGH_MIE)
+    {
+        sky_color = rayleighGetSkyColor(eye, look, sun_position);
     }
     else
     {
@@ -458,5 +332,12 @@ Color skyGetSunColor(SkyDefinition* definition)
 
 Color skyGetZenithColor(SkyDefinition* definition)
 {
-    return definition->model_custom.zenith_color;
+    Vector3 look = {0.0, 1.0, 0.0};
+    return skyGetColor(definition, NULL, VECTOR_ZERO, look);
+}
+
+Color skyGetHorizonColor(SkyDefinition* definition)
+{
+    Vector3 look = {0.0, 0.0, 1.0};
+    return skyGetColor(definition, NULL, VECTOR_ZERO, look);
 }

lib_paysages/sky.h

@@ -13,7 +13,8 @@ extern "C" {
 typedef enum
 {
     SKY_MODEL_CUSTOM = 0,
-    SKY_MODEL_PREETHAM = 1
+    SKY_MODEL_RAYLEIGH_MIE = 1,
+    SKY_MODEL_PREETHAM = 2
 } SkyModel;
 
 typedef struct
@@ -53,6 +54,7 @@ void skyRender(SkyDefinition* definition, Renderer* renderer);
 Vector3 skyGetSunDirection(SkyDefinition* definition);
 Color skyGetSunColor(SkyDefinition* definition);
 Color skyGetZenithColor(SkyDefinition* definition);
+Color skyGetHorizonColor(SkyDefinition* definition);
 
 #ifdef __cplusplus
 }