Stellarium 0.13.1
Skylight.hpp
1 /*
2  * Copyright (C) 2003 Fabien Chereau
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
17  */
18 
19 // Class which computes the daylight sky color
20 // Fast implementation of the algorithm from the article
21 // "A Practical Analytic Model for Daylight" by A. J. Preetham, Peter Shirley and Brian Smits.
22 
23 #ifndef _SKYLIGHT_HPP_
24 #define _SKYLIGHT_HPP_
25 
26 #include "StelUtils.hpp"
27 
28 #include <cmath>
29 #include <QDebug>
30 
31 typedef struct {
32  float zenithAngle; // zenithAngle : angular distance to the zenith in radian
33  float distSun; // distSun : angular distance to the sun in radian
34  float color[3]; // 3 component color, can be RGB or CIE color system
36 
37 typedef struct {
38  float pos[3]; // Vector to the position (vertical = pos[2])
39  float color[3]; // 3 component color, can be RGB or CIE color system
41 
42 class Skylight
43 {
44 public:
45  Skylight();
46  virtual ~Skylight();
47  // Set the fixed parameters and precompute what can be
48  // This function has to be called once before any call to get_*_value()
49  void setParams(float sunZenithAngle, float turbidity);
50  // Compute the sky color at the given position in the xyY color system and store it in position.color
51  // void getxyYValue(skylightStruct * position);
52  // Return the current zenith color
53  inline void getZenithColor(float * v) const;
54 
55  // Same functions but in vector mode : faster because prevents extra cosine calculations
56  // The position vectors MUST be normalized, and the vertical z component is the third one
57  void setParamsv(const float * sunPos, float turbidity);
58 
59  // Compute the sky color at the given position in the CIE color system and store it in p.color
60  // p.color[0] is CIE x color component
61  // p.color[1] is CIE y color component
62  // p.color[2] is undefined (CIE Y color component (luminance) if uncommented)
63  void getxyYValuev(skylightStruct2& p) const
64  {
65  const float cosDistSun = sunPos[0]*p.pos[0] + sunPos[1]*p.pos[1] + sunPos[2]*p.pos[2];
66  const float distSun = StelUtils::fastAcos(cosDistSun);
67  const float cosDistSun_q = cosDistSun*cosDistSun;
68 
69  Q_ASSERT(p.pos[2] >= 0.f);
70  const float oneOverCosZenithAngle = (p.pos[2]==0.) ? 1e99 : 1.f / p.pos[2];
71  p.color[0] = term_x * (1.f + Ax * std::exp(Bx*oneOverCosZenithAngle))
72  * (1.f + Cx * std::exp(Dx*distSun) + Ex * cosDistSun_q);
73 
74  p.color[1] = term_y * (1.f + Ay * std::exp(By*oneOverCosZenithAngle))
75  * (1.f + Cy * std::exp(Dy*distSun) + Ey * cosDistSun_q);
76 
77 // p.color[2] = term_Y * (1.f + AY * std::exp(BY*oneOverCosZenithAngle))
78 // * (1.f + CY * std::exp(DY*distSun) + EY * cosDistSun_q);
79 
80 
81  if (/*p.color[2] < 0. || */p.color[0] < 0. || p.color[1] < 0.)
82  {
83  p.color[0] = 0.25;
84  p.color[1] = 0.25;
85  p.color[2] = 0.;
86  }
87  }
88 
89  void getShadersParams(Vec3f& asunPos, float& aterm_x, float& aAx, float& aBx, float& aCx, float& aDx, float& aEx,
90  float& aterm_y, float& aAy, float& aBy, float& aCy, float& aDy, float& aEy) const
91  {
92  asunPos=sunPos;
93  aterm_x=term_x;aAx=Ax;aBx=Bx;aCx=Cx;aDx=Dx;aEx=Ex;
94  aterm_y=term_y;aAy=Ay;aBy=By;aCy=Cy;aDy=Dy;aEy=Ey;
95  }
96 
97 
98 
99 private:
100  float thetas; // angular distance between the zenith and the sun in radian
101  float T; // Turbidity : i.e. sky "clarity"
102  // 1 : pure air
103  // 2 : exceptionnally clear
104  // 4 : clear
105  // 8 : light haze
106  // 25 : haze
107  // 64 : thin fog
108 
109  // Computed variables depending on the 2 above
110  float zenithLuminance; // Y color component of the CIE color at zenith (luminance)
111  float zenithColorX; // x color component of the CIE color at zenith
112  float zenithColorY; // y color component of the CIE color at zenith
113 
114  float eyeLumConversion; // luminance conversion for an eye adapted to screen luminance
115  // (around 40 cd/m^2)
116 
117  float AY, BY, CY, DY, EY; // Distribution coefficients for the luminance distribution function
118  float Ax, Bx, Cx, Dx, Ex; // Distribution coefficients for x distribution function
119  float Ay, By, Cy, Dy, Ey; // Distribution coefficients for y distribution function
120 
121  float term_x; // Precomputed term for x calculation
122  float term_y; // Precomputed term for y calculation
123  float term_Y; // Precomputed term for luminance calculation
124 
125  float sunPos[3];
126 
127  // Compute CIE Y (luminance) for zenith in cd/m^2
128  inline void computeZenithLuminance(void);
129  // Compute CIE x and y color components
130  inline void computeZenithColor(void);
131  // Compute the luminance distribution coefficients
132  inline void computeLuminanceDistributionCoefs(void);
133  // Compute the color distribution coefficients
134  inline void computeColorDistributionCoefs(void);
135 };
136 
137 // Return the current zenith color in xyY color system
138 inline void Skylight::getZenithColor(float * v) const
139 {
140  v[0] = zenithColorX;
141  v[1] = zenithColorY;
142  v[2] = zenithLuminance;
143 }
144 
145 // Compute CIE luminance for zenith in cd/m^2
146 inline void Skylight::computeZenithLuminance(void)
147 {
148  zenithLuminance = 1000.f * ((4.0453f*T - 4.9710f) * std::tan( (0.4444f - T/120.f) * (M_PI-2.f*thetas) ) -
149  0.2155f*T + 2.4192f);
150  if (zenithLuminance<=0.f) zenithLuminance=0.00000000001;
151 }
152 
153 // Compute CIE x and y color components
154 // Edit: changed some coefficients to get new sky color
155 inline void Skylight::computeZenithColor(void)
156 {
157  static float thetas2;
158  static float thetas3;
159  static float T2;
160 
161  thetas2 = thetas * thetas;
162  thetas3 = thetas2 * thetas;
163  T2 = T * T;
164 
165  zenithColorX = ( 0.00216f*thetas3 - 0.00375f*thetas2 + 0.00209f*thetas) * T2 +
166  (-0.02903f*thetas3 + 0.06377f*thetas2 - 0.03202f*thetas + 0.00394f) * T +
167  ( 0.10169f*thetas3 - 0.21196f*thetas2 + 0.06052f*thetas + 0.25886f);
168 
169  zenithColorY = ( 0.00275f*thetas3 - 0.00610f*thetas2 + 0.00317f*thetas) * T2 +
170  (-0.04214f*thetas3 + 0.08970f*thetas2 - 0.04153f*thetas + 0.00516f) * T +
171  ( 0.14535f*thetas3 - 0.26756f*thetas2 + 0.06670f*thetas + 0.26688f);
172 
173 }
174 
175 // Compute the luminance distribution coefficients
176 // Edit: changed some coefficients to get new sky color
177 inline void Skylight::computeLuminanceDistributionCoefs(void)
178 {
179  AY = 0.2787f*T - 1.0630f;
180  BY =-0.3554f*T + 0.4275f;
181  CY =-0.0227f*T + 6.3251f;
182  DY = 0.1206f*T - 2.5771f;
183  EY =-0.0670f*T + 0.3703f;
184  // with BY>0 the formulas in getxyYValuev make no sense
185  Q_ASSERT(BY <= 0.0);
186 }
187 
188 // Compute the color distribution coefficients
189 // Edit: changed some coefficients to get new sky color
190 inline void Skylight::computeColorDistributionCoefs(void)
191 {
192  Ax =-0.0148f*T - 0.1703f;
193  Bx =-0.0664f*T + 0.0011f;
194  Cx =-0.0005f*T + 0.2127f;
195  Dx =-0.0641f*T - 0.8992f;
196  Ex =-0.0035f*T + 0.0453f;
197 
198  Ay =-0.0131f*T - 0.2498f;
199  By =-0.0951f*T + 0.0092f;
200  Cy =-0.0082f*T + 0.2404f;
201  Dy =-0.0438f*T - 1.0539f;
202  Ey =-0.0109f*T + 0.0531f;
203  // with Bx,By>0 the formulas in getxyYValuev make no sense
204  Q_ASSERT(Bx <= 0.0);
205  Q_ASSERT(By <= 0.0);
206 }
207 
208 
209 #endif // _SKYLIGHT_H_
210 
float fastAcos(const float x)
Compute acos(x) The taylor serie is not accurate around x=1 and x=-1.
Definition: StelUtils.hpp:291