From 394ba7d4c596cfbe9e823be26e0c1bebc73ea71f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Lemaire?= Date: Sun, 24 Jun 2012 12:33:59 +0000 Subject: [PATCH] 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 | 2 + gui_qt/formatmosphere.cpp | 2 +- gui_qt/formsky.cpp | 2 +- gui_qt/formterrain.cpp | 2 +- i18n/paysages_fr.ts | 15 +++- lib_paysages/atmosphere.c | 2 +- lib_paysages/preetham.c | 137 ++++++++++++++++++++++++++++++++++ lib_paysages/preetham.h | 20 +++++ lib_paysages/rayleigh.c | 30 ++++++++ lib_paysages/rayleigh.h | 20 +++++ lib_paysages/sky.c | 151 ++++---------------------------------- lib_paysages/sky.h | 4 +- 12 files changed, 244 insertions(+), 143 deletions(-) create mode 100644 lib_paysages/preetham.c create mode 100644 lib_paysages/preetham.h create mode 100644 lib_paysages/rayleigh.c create mode 100644 lib_paysages/rayleigh.h diff --git a/TODO b/TODO index ab18220..7d1cf07 100644 --- a/TODO +++ b/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 ? diff --git a/gui_qt/formatmosphere.cpp b/gui_qt/formatmosphere.cpp index 8c8942e..abd1dc4 100644 --- a/gui_qt/formatmosphere.cpp +++ b/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(); diff --git a/gui_qt/formsky.cpp b/gui_qt/formsky.cpp index c24a148..0632337 100644 --- a/gui_qt/formsky.cpp +++ b/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); diff --git a/gui_qt/formterrain.cpp b/gui_qt/formterrain.cpp index d3cb653..c780432 100644 --- a/gui_qt/formterrain.cpp +++ b/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(); diff --git a/i18n/paysages_fr.ts b/i18n/paysages_fr.ts index 9138f49..2b232f6 100644 --- a/i18n/paysages_fr.ts +++ b/i18n/paysages_fr.ts @@ -335,8 +335,12 @@ Maintenir Ctrl : Plus rapide + Lock on horizon color + + + Lock color on haze - Verrouiller sur la couleur de la brume + Verrouiller sur la couleur de la brume @@ -608,12 +612,17 @@ Maintenir Ctrl : Plus rapide - Mixed Preetham/Shirley approximation + Custom model - Custom model + Rayleigh/Mie scattering + + + + + Preetham/Shirley analytic model diff --git a/lib_paysages/atmosphere.c b/lib_paysages/atmosphere.c index ded713f..8b18604 100644 --- a/lib_paysages/atmosphere.c +++ b/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); } } diff --git a/lib_paysages/preetham.c b/lib_paysages/preetham.c new file mode 100644 index 0000000..e74c433 --- /dev/null +++ b/lib_paysages/preetham.c @@ -0,0 +1,137 @@ +#include "preetham.h" + +#include + +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; +} diff --git a/lib_paysages/preetham.h b/lib_paysages/preetham.h new file mode 100644 index 0000000..2e13f4d --- /dev/null +++ b/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 diff --git a/lib_paysages/rayleigh.c b/lib_paysages/rayleigh.c new file mode 100644 index 0000000..969418d --- /dev/null +++ b/lib_paysages/rayleigh.c @@ -0,0 +1,30 @@ +#include "rayleigh.h" + +#include + +/*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; +} diff --git a/lib_paysages/rayleigh.h b/lib_paysages/rayleigh.h new file mode 100644 index 0000000..645d524 --- /dev/null +++ b/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 diff --git a/lib_paysages/sky.c b/lib_paysages/sky.c index d245009..ba987e2 100644 --- a/lib_paysages/sky.c +++ b/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); } diff --git a/lib_paysages/sky.h b/lib_paysages/sky.h index 497274b..4912e1b 100644 --- a/lib_paysages/sky.h +++ b/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 }