From 9e21a37c048f0e7c0bbf6d1021d11731462055c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Lemaire?= Date: Wed, 20 Jun 2012 20:29:58 +0000 Subject: [PATCH] paysages : Preetham approximation for sky (WIP). git-svn-id: https://subversion.assembla.com/svn/thunderk/paysages@354 b1fd45b6-86a6-48da-8261-f70d1f35bdcc --- lib_paysages/sky.c | 176 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 174 insertions(+), 2 deletions(-) diff --git a/lib_paysages/sky.c b/lib_paysages/sky.c index 2b7d004..df7b4d4 100644 --- a/lib_paysages/sky.c +++ b/lib_paysages/sky.c @@ -150,18 +150,190 @@ int skyGetLights(SkyDefinition* sky, LightDefinition* lights, int max_lights) return nblights; } +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 double _perezFunction(double A, double B, double C, double D, double E, double Theta, double Gamma) +{ + double cosGamma = cos(Gamma); + return (1.0 + A * exp(B / cos(Theta))) * (1.0 + C * exp(D * Gamma) + E * cosGamma * cosGamma); +} + +typedef struct +{ + double x; + double y; + double Y; +} ColorxyY; + +/* Distribution coefficients for the luminance(Y) distribution function */ +static double YDC[5][2] = { { 0.1787, - 1.4630}, + {-0.3554, 0.4275}, + {-0.0227, 5.3251}, + { 0.1206, - 2.5771}, + {-0.0670, 0.3703} }; + +/* Distribution coefficients for the x distribution function */ +static double xDC[5][2] = { {-0.0193, -0.2592}, + {-0.0665, 0.0008}, + {-0.0004, 0.2125}, + {-0.0641, -0.8989}, + {-0.0033, 0.0452} }; + +/* Distribution coefficients for the y distribution function */ +static double yDC[5][2] = { {-0.0167, -0.2608}, + {-0.0950, 0.0092}, + {-0.0079, 0.2102}, + {-0.0441, -1.6537}, + {-0.0109, 0.0529} }; + +/* Zenith x value */ +static double xZC[3][4] = { {0.00166, -0.00375, 0.00209, 0}, + {-0.02903, 0.06377, -0.03203, 0.00394}, + {0.11693, -0.21196, 0.06052, 0.25886} }; +/* Zenith y value */ +static double yZC[3][4] = { { 0.00275, -0.00610, 0.00317, 0}, + {-0.04214, 0.08970, -0.04153, 0.00516}, + {0.15346, -0.26756, 0.06670, 0.26688} }; + +static double _distribution(double A, double B, double C, double D, double E, double Theta, double Gamma, double ThetaSun) +{ + double f0 = _perezFunction(A,B,C,D,E,Theta,Gamma); + double f1 = _perezFunction(A,B,C,D,E,0,ThetaSun); + return(f0/f1); +} + +static double _chromaticity( double ZC[3][4], double ThetaSun, double turbidity) +{ + double t1 = ThetaSun; + double t2 = t1*t1; + double t3 = t2 * t1; + + double c = (ZC[0][0]*t3 + ZC[0][1]*t2 + ZC[0][2]*t1 + ZC[0][3])* turbidity * turbidity + + (ZC[1][0]*t3 + ZC[1][1]*t2 + ZC[1][2]*t1 + ZC[1][3])* turbidity + + (ZC[2][0]*t3 + ZC[2][1]*t2 + ZC[2][2]*t1 + ZC[2][3]); + return c; +} + +static Color _toColor(ColorxyY src) +{ + double X, Y, Z; + Color dest; + + /* Convert to XYZ space */ + X = src.x * (src.Y / src.y); + Y = src.Y; + Z = (1.0 - src.x - src.y)* (src.Y/src.y); + + /* Convert to RGB */ + dest.r = 3.240479 * X - 1.537150 * Y - 0.498535 * Z; + dest.g = - 0.969256 * X + 1.875991 * Y + 0.041556 * Z; + dest.b = 0.055648 * X - 0.204043 * Y + 1.057311 * Z; + dest.a = 1.0; + + return dest; +} + +static Color _preethamApproximate(SkyDefinition* definition, double theta, double phi) +{ + double thetaSun; + double phiSun; + double gamma; + double turbidity = 3.0; + ColorxyY skycolor; + ColorxyY Zenith; + double A,B,C,D,E; + double d,chi; + + /* Handle angles */ + if (definition->daytime <= 0.5) + { + if (definition->daytime <= 0.25) + { + thetaSun = M_PI / 2.0 - 0.00001; + } + else + { + thetaSun = M_PI - definition->daytime * 2.0 * M_PI; + } + phiSun = 0.0; + } + else + { + if (definition->daytime >= 0.75) + { + thetaSun = M_PI / 2.0 - 0.00001; + } + else + { + thetaSun = (definition->daytime - 0.5) * 2.0 * M_PI; + } + phiSun = M_PI; + } + gamma = _angleBetween(theta, phi, thetaSun, phiSun); + + /* Zenith luminance */ + chi = (4.0/9.0 - turbidity/120.0)*(M_PI - 2*thetaSun); + Zenith.Y = (4.0453*turbidity - 4.9710)*tan(chi) - 0.2155*turbidity + 2.4192; + if (Zenith.Y < 0.0) Zenith.Y = -Zenith.Y; + + A = YDC[0][0]*turbidity + YDC[0][1]; + B = YDC[1][0]*turbidity + YDC[1][1]; + C = YDC[2][0]*turbidity + YDC[2][1]; + D = YDC[3][0]*turbidity + YDC[3][1]; + E = YDC[4][0]*turbidity + YDC[4][1]; + + /* Sky luminance */ + d = _distribution(A, B, C, D, E, theta, gamma, thetaSun); + skycolor.Y = Zenith.Y * d; + + /* Zenith x */ + Zenith.x = _chromaticity( xZC, thetaSun, turbidity ); + A = xDC[0][0]*turbidity + xDC[0][1]; + B = xDC[1][0]*turbidity + xDC[1][1]; + C = xDC[2][0]*turbidity + xDC[2][1]; + D = xDC[3][0]*turbidity + xDC[3][1]; + E = xDC[4][0]*turbidity + xDC[4][1]; + + /* Sky x */ + d = _distribution(A, B, C, D, E, theta, gamma, thetaSun); + skycolor.x = Zenith.x * d; + + /* Zenith y */ + Zenith.y = _chromaticity( yZC, thetaSun, turbidity ); + A = yDC[0][0]*turbidity + yDC[0][1]; + B = yDC[1][0]*turbidity + yDC[1][1]; + C = yDC[2][0]*turbidity + yDC[2][1]; + D = yDC[3][0]*turbidity + yDC[3][1]; + E = yDC[4][0]*turbidity + yDC[4][1]; + + /* Sky y */ + d = _distribution(A, B, C, D, E, theta, gamma, thetaSun); + skycolor.y = Zenith.y * d; + + /* Small hack on luminance value */ + //skycolor.Y = 1 - exp(-(1.0/25.0) * skycolor.Y); + + return _toColor(skycolor); +} + Color skyGetColor(SkyDefinition* definition, Renderer* renderer, Vector3 eye, Vector3 look) { double dist; Vector3 sun_position; Color sun_color, sky_color; - + sun_position = skyGetSunDirection(definition); look = v3Normalize(look); dist = v3Norm(v3Sub(look, sun_position)); - sky_color = colorGradationGet(definition->_sky_gradation, look.y * 0.5 + 0.5); + sky_color = _preethamApproximate(definition, M_PI/2.0 - asin(look.y), acos(v3Dot(look, sun_position))); if (dist < definition->sun_radius + definition->sun_halo_size) { sun_color = colorGradationGet(definition->sun_color, definition->daytime);