4f57a046d68965a61842343f07ab3d5ccc257891
[blender.git] / source / blender / render / intern / source / sunsky.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
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 Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * ***** END GPL LICENSE BLOCK *****
19  */
20
21 /** \file blender/render/intern/source/sunsky.c
22  *  \ingroup render
23  *
24  * This feature comes from Preetham paper on "A Practical Analytic Model for Daylight"
25  * and example code from Brian Smits, another author of that paper in
26  * http://www.cs.utah.edu/vissim/papers/sunsky/code/
27  */
28
29 #include "sunsky.h"
30 #include "BLI_math.h"
31 #include "BKE_global.h"
32
33 /**
34  * These macros are defined for vector operations
35  * */
36
37 /**
38  * compute v1 = v2 op v3
39  * v1, v2 and v3 are vectors contains 3 float
40  * */
41 #define VEC3OPV(v1, v2, op, v3)                                               \
42         {                                                                         \
43                 v1[0] = (v2[0] op v3[0]);                                             \
44                 v1[1] = (v2[1] op v3[1]);                                             \
45                 v1[2] = (v2[2] op v3[2]);                                             \
46         } (void)0
47
48 /**
49  * compute v1 = v2 op f1
50  * v1, v2 are vectors contains 3 float
51  * and f1 is a float
52  * */
53 #define VEC3OPF(v1, v2, op, f1)                                               \
54         {                                                                         \
55                 v1[0] = (v2[0] op(f1));                                              \
56                 v1[1] = (v2[1] op(f1));                                              \
57                 v1[2] = (v2[2] op(f1));                                              \
58         } (void)0
59
60 /**
61  * compute v1 = f1 op v2
62  * v1, v2 are vectors contains 3 float
63  * and f1 is a float
64  * */
65 #define FOPVEC3(v1, f1, op, v2)                                               \
66         {                                                                         \
67                 v1[0] = ((f1) op v2[0]);                                              \
68                 v1[1] = ((f1) op v2[1]);                                              \
69                 v1[2] = ((f1) op v2[2]);                                              \
70         } (void)0
71
72 /**
73  * ClipColor:
74  * clip a color to range [0, 1];
75  * */
76 void ClipColor(float c[3])
77 {
78         if (c[0] > 1.0f) c[0] = 1.0f;
79         if (c[0] < 0.0f) c[0] = 0.0f;
80         if (c[1] > 1.0f) c[1] = 1.0f;
81         if (c[1] < 0.0f) c[1] = 0.0f;
82         if (c[2] > 1.0f) c[2] = 1.0f;
83         if (c[2] < 0.0f) c[2] = 0.0f;
84 }
85
86 /**
87  * AngleBetween:
88  * compute angle between to direction 
89  * all angles are in radians
90  * */
91 static float AngleBetween(float thetav, float phiv, float theta, float phi)
92 {
93         float cospsi = sin(thetav) * sin(theta) * cos(phi - phiv) + cos(thetav) * cos(theta);
94
95         if (cospsi > 1.0f)
96                 return 0;
97         if (cospsi < -1.0f)
98                 return M_PI;
99
100         return acos(cospsi);
101 }
102
103 /**
104  * DirectionToThetaPhi:
105  * this function convert a direction to it's theta and phi value
106  * parameters:
107  * toSun: contains direction information
108  * theta, phi, are return values from this conversion
109  * */
110 static void DirectionToThetaPhi(float *toSun, float *theta, float *phi)
111 {
112         *theta = acos(toSun[2]);
113         if (fabs(*theta) < 1e-5)
114                 *phi = 0;
115         else
116                 *phi = atan2(toSun[1], toSun[0]);
117 }
118
119 /**
120  * PerezFunction:
121  * compute perez function value based on input paramters
122  * */
123 static float PerezFunction(struct SunSky *sunsky, const float *lam, float theta, float gamma, float lvz)
124 {
125         float den, num;
126         
127         den = ((1 + lam[0] * expf(lam[1])) *
128                (1 + lam[2] * expf(lam[3] * sunsky->theta) + lam[4] * cosf(sunsky->theta) * cosf(sunsky->theta)));
129         
130         num = ((1 + lam[0] * expf(lam[1] / cosf(theta))) *
131                (1 + lam[2] * expf(lam[3] * gamma) + lam[4] * cosf(gamma) * cosf(gamma)));
132         
133         return(lvz * num / den);
134 }
135
136 /**
137  * InitSunSky:
138  * this function compute some sun,sky parameters according to input parameters and also initiate some other sun, sky parameters
139  * parameters:
140  * sunSky, is a structure that contains information about sun, sky and atmosphere, in this function, most of its values initiated
141  * turb, is atmosphere turbidity
142  * toSun, contains sun direction
143  * horizon_brighness, controls the brightness of the horizon colors
144  * spread, controls colors spreed at horizon
145  * sun_brightness, controls sun's brightness
146  * sun_size, controls sun's size
147  * back_scatter, controls back scatter light
148  * */
149 void InitSunSky(struct SunSky *sunsky, float turb, const float toSun[3], float horizon_brightness,
150                 float spread, float sun_brightness, float sun_size, float back_scatter,
151                 float skyblendfac, short skyblendtype, float sky_exposure, float sky_colorspace)
152 {
153         float theta2;
154         float theta3;
155         float T;
156         float T2;
157         float chi;
158
159         sunsky->turbidity = turb;
160
161         sunsky->horizon_brightness = horizon_brightness;
162         sunsky->spread = spread;
163         sunsky->sun_brightness = sun_brightness;
164         sunsky->sun_size = sun_size;
165         sunsky->backscattered_light = back_scatter;
166         sunsky->skyblendfac = skyblendfac;
167         sunsky->skyblendtype = skyblendtype;
168         sunsky->sky_exposure = -sky_exposure;
169         sunsky->sky_colorspace = sky_colorspace;
170         
171         sunsky->toSun[0] = toSun[0];
172         sunsky->toSun[1] = toSun[1];
173         sunsky->toSun[2] = toSun[2];
174
175         DirectionToThetaPhi(sunsky->toSun, &sunsky->theta, &sunsky->phi);
176
177         sunsky->sunSolidAngle = 0.25 * M_PI * 1.39 * 1.39 / (150 * 150);   /* = 6.7443e-05 */
178
179         theta2 = sunsky->theta * sunsky->theta;
180         theta3 = theta2 * sunsky->theta;
181         T = turb;
182         T2 = turb * turb;
183
184         chi = (4.0f / 9.0f - T / 120.0f) * ((float)M_PI - 2.0f * sunsky->theta);
185         sunsky->zenith_Y = (4.0453f * T - 4.9710f) * tanf(chi) - 0.2155f * T + 2.4192f;
186         sunsky->zenith_Y *= 1000;   /* conversion from kcd/m^2 to cd/m^2 */
187
188         if (sunsky->zenith_Y <= 0)
189                 sunsky->zenith_Y = 1e-6;
190         
191         sunsky->zenith_x =
192             (+0.00165f * theta3 - 0.00374f * theta2 + 0.00208f * sunsky->theta + 0.0f) * T2 +
193             (-0.02902f * theta3 + 0.06377f * theta2 - 0.03202f * sunsky->theta + 0.00394f) * T +
194             (+0.11693f * theta3 - 0.21196f * theta2 + 0.06052f * sunsky->theta + 0.25885f);
195
196         sunsky->zenith_y =
197             (+0.00275f * theta3 - 0.00610f * theta2 + 0.00316f * sunsky->theta + 0.0f) * T2 +
198             (-0.04214f * theta3 + 0.08970f * theta2 - 0.04153f * sunsky->theta + 0.00515f) * T +
199             (+0.15346f * theta3 - 0.26756f * theta2 + 0.06669f * sunsky->theta + 0.26688f);
200
201         
202         sunsky->perez_Y[0] = 0.17872f * T - 1.46303f;
203         sunsky->perez_Y[1] = -0.35540f * T + 0.42749f;
204         sunsky->perez_Y[2] = -0.02266f * T + 5.32505f;
205         sunsky->perez_Y[3] = 0.12064f * T - 2.57705f;
206         sunsky->perez_Y[4] = -0.06696f * T + 0.37027f;
207
208         sunsky->perez_x[0] = -0.01925f * T - 0.25922f;
209         sunsky->perez_x[1] = -0.06651f * T + 0.00081f;
210         sunsky->perez_x[2] = -0.00041f * T + 0.21247f;
211         sunsky->perez_x[3] = -0.06409f * T - 0.89887f;
212         sunsky->perez_x[4] = -0.00325f * T + 0.04517f;
213
214         sunsky->perez_y[0] = -0.01669f * T - 0.26078f;
215         sunsky->perez_y[1] = -0.09495f * T + 0.00921f;
216         sunsky->perez_y[2] = -0.00792f * T + 0.21023f;
217         sunsky->perez_y[3] = -0.04405f * T - 1.65369f;
218         sunsky->perez_y[4] = -0.01092f * T + 0.05291f;
219         
220         /* suggested by glome in 
221          * http://projects.blender.org/tracker/?func=detail&atid=127&aid=8063&group_id=9*/
222         sunsky->perez_Y[0] *= sunsky->horizon_brightness;
223         sunsky->perez_x[0] *= sunsky->horizon_brightness;
224         sunsky->perez_y[0] *= sunsky->horizon_brightness;
225         
226         sunsky->perez_Y[1] *= sunsky->spread;
227         sunsky->perez_x[1] *= sunsky->spread;
228         sunsky->perez_y[1] *= sunsky->spread;
229
230         sunsky->perez_Y[2] *= sunsky->sun_brightness;
231         sunsky->perez_x[2] *= sunsky->sun_brightness;
232         sunsky->perez_y[2] *= sunsky->sun_brightness;
233         
234         sunsky->perez_Y[3] *= sunsky->sun_size;
235         sunsky->perez_x[3] *= sunsky->sun_size;
236         sunsky->perez_y[3] *= sunsky->sun_size;
237         
238         sunsky->perez_Y[4] *= sunsky->backscattered_light;
239         sunsky->perez_x[4] *= sunsky->backscattered_light;
240         sunsky->perez_y[4] *= sunsky->backscattered_light;
241 }
242
243 /**
244  * GetSkyXYZRadiance:
245  * this function compute sky radiance according to a view parameters `theta' and `phi'and sunSky values
246  * parameters:
247  * sunSky, sontains sun and sky parameters
248  * theta, is sun's theta
249  * phi, is sun's phi
250  * color_out, is computed color that shows sky radiance in XYZ color format
251  * */
252 void GetSkyXYZRadiance(struct SunSky *sunsky, float theta, float phi, float color_out[3])
253 {
254         float gamma;
255         float x, y, Y, X, Z;
256         float hfade = 1, nfade = 1;
257
258
259         if (theta > (0.5f * (float)M_PI)) {
260                 hfade = 1.0f - (theta * (float)M_1_PI - 0.5f) * 2.0f;
261                 hfade = hfade * hfade * (3.0f - 2.0f * hfade);
262                 theta = 0.5 * M_PI;
263         }
264
265         if (sunsky->theta > (0.5f * (float)M_PI)) {
266                 if (theta <= 0.5f * (float)M_PI) {
267                         nfade = 1.0f - (0.5f - theta * (float)M_1_PI) * 2.0f;
268                         nfade *= 1.0f - (sunsky->theta * (float)M_1_PI - 0.5f) * 2.0f;
269                         nfade = nfade * nfade * (3.0f - 2.0f * nfade);
270                 }
271         }
272
273         gamma = AngleBetween(theta, phi, sunsky->theta, sunsky->phi);
274         
275         /* Compute xyY values */
276         x = PerezFunction(sunsky, sunsky->perez_x, theta, gamma, sunsky->zenith_x);
277         y = PerezFunction(sunsky, sunsky->perez_y, theta, gamma, sunsky->zenith_y);
278         Y = 6.666666667e-5f * nfade * hfade * PerezFunction(sunsky, sunsky->perez_Y, theta, gamma, sunsky->zenith_Y);
279
280         if (sunsky->sky_exposure != 0.0f)
281                 Y = 1.0 - exp(Y * sunsky->sky_exposure);
282         
283         X = (x / y) * Y;
284         Z = ((1 - x - y) / y) * Y;
285
286         color_out[0] = X;
287         color_out[1] = Y;
288         color_out[2] = Z;
289 }
290
291 /**
292  * GetSkyXYZRadiancef:
293  * this function compute sky radiance according to a view direction `varg' and sunSky values
294  * parameters:
295  * sunSky, sontains sun and sky parameters
296  * varg, shows direction
297  * color_out, is computed color that shows sky radiance in XYZ color format
298  * */
299 void GetSkyXYZRadiancef(struct SunSky *sunsky, const float varg[3], float color_out[3])
300 {
301         float theta, phi;
302         float v[3];
303
304         copy_v3_v3(v, (float *)varg);
305         normalize_v3(v);
306
307         if (v[2] < 0.001f) {
308                 v[2] = 0.001f;
309                 normalize_v3(v);
310         }
311
312         DirectionToThetaPhi(v, &theta, &phi);
313         GetSkyXYZRadiance(sunsky, theta, phi, color_out);
314 }
315
316 /**
317  * ComputeAttenuatedSunlight:
318  * this function compute attenuated sun light based on sun's theta and atmosphere turbidity
319  * parameters:
320  * theta, is sun's theta
321  * turbidity: is atmosphere turbidity
322  * fTau: contains computed attenuated sun light
323  * */
324 static void ComputeAttenuatedSunlight(float theta, int turbidity, float fTau[3])
325 {
326         float fBeta;
327         float fTauR, fTauA;
328         float m;
329         float fAlpha;
330
331         int i;
332         float fLambda[3]; 
333         fLambda[0] = 0.65f;
334         fLambda[1] = 0.57f;
335         fLambda[2] = 0.475f;
336
337         fAlpha = 1.3f;
338         fBeta = 0.04608365822050f * turbidity - 0.04586025928522f;
339         
340         m =  1.0f / (cosf(theta) + 0.15f * powf(93.885f - theta / (float)M_PI * 180.0f, -1.253f));
341
342         for (i = 0; i < 3; i++) {
343                 /* Rayleigh Scattering */
344                 fTauR = expf(-m * 0.008735f * powf(fLambda[i], (float)(-4.08f)));
345
346                 /* Aerosal (water + dust) attenuation */
347                 fTauA = exp(-m * fBeta * powf(fLambda[i], -fAlpha));
348
349                 fTau[i] = fTauR * fTauA;
350         }
351 }
352
353 /**
354  * InitAtmosphere:
355  * this function initiate sunSky structure with user input parameters.
356  * parameters:
357  * sunSky, contains information about sun, and in this function some atmosphere parameters will initiated
358  * sun_intens, shows sun intensity value
359  * mief, Mie scattering factor this factor currently call with 1.0 
360  * rayf, Rayleigh scattering factor, this factor currently call with 1.0
361  * inscattf, inscatter light factor that range from 0.0 to 1.0, 0.0 means no inscatter light and 1.0 means full inscatter light
362  * extincf, extinction light factor that range from 0.0 to 1.0, 0.0 means no extinction and 1.0 means full extinction
363  * disf, is distance factor, multiplied to pixle's z value to compute each pixle's distance to camera,
364  * */
365 void InitAtmosphere(struct SunSky *sunSky, float sun_intens, float mief, float rayf,
366                     float inscattf, float extincf, float disf)
367 {
368         const float pi = M_PI;
369         const float n = 1.003f;  /* refractive index */
370         const float N = 2.545e25;
371         const float pn = 0.035f;
372         const float T = 2.0f;
373         float fTemp, fTemp2, fTemp3, fBeta, fBetaDash;
374         float c = (6.544f * T - 6.51f) * 1e-17f;
375         float K[3] = {0.685f, 0.679f, 0.670f};
376         float vBetaMieTemp[3];
377         
378         float fLambda[3], fLambda2[3], fLambda4[3];
379         float vLambda2[3];
380         float vLambda4[3];
381         
382         int i;
383
384         sunSky->atm_SunIntensity = sun_intens;
385         sunSky->atm_BetaMieMultiplier  = mief;
386         sunSky->atm_BetaRayMultiplier = rayf;
387         sunSky->atm_InscatteringMultiplier = inscattf;
388         sunSky->atm_ExtinctionMultiplier = extincf;
389         sunSky->atm_DistanceMultiplier = disf;
390                 
391         sunSky->atm_HGg = 0.8;
392
393         fLambda[0]  = 1 / 650e-9f;
394         fLambda[1]  = 1 / 570e-9f;
395         fLambda[2]  = 1 / 475e-9f;
396         for (i = 0; i < 3; i++) {
397                 fLambda2[i] = fLambda[i] * fLambda[i];
398                 fLambda4[i] = fLambda2[i] * fLambda2[i];
399         }
400
401         vLambda2[0] = fLambda2[0];
402         vLambda2[1] = fLambda2[1];
403         vLambda2[2] = fLambda2[2];
404  
405         vLambda4[0] = fLambda4[0];
406         vLambda4[1] = fLambda4[1];
407         vLambda4[2] = fLambda4[2];
408
409         /* Rayleigh scattering constants. */
410         fTemp = pi * pi * (n * n - 1) * (n * n - 1) * (6 + 3 * pn) / (6 - 7 * pn) / N;
411         fBeta = 8 * fTemp * pi / 3;
412
413         VEC3OPF(sunSky->atm_BetaRay, vLambda4, *, fBeta);
414         fBetaDash = fTemp / 2;
415         VEC3OPF(sunSky->atm_BetaDashRay, vLambda4, *, fBetaDash);
416
417
418         /* Mie scattering constants. */
419         fTemp2 = 0.434f * c * (2 * pi) * (2 * pi) * 0.5f;
420         VEC3OPF(sunSky->atm_BetaDashMie, vLambda2, *, fTemp2);
421         
422         fTemp3 = 0.434f * c * pi * (2 * pi) * (2 * pi);
423         
424         VEC3OPV(vBetaMieTemp, K, *, fLambda);
425         VEC3OPF(sunSky->atm_BetaMie, vBetaMieTemp, *, fTemp3);
426         
427 }
428
429 /**
430  * AtmospherePixleShader:
431  * this function apply atmosphere effect on a pixle color `rgb' at distance `s'
432  * parameters:
433  * sunSky, contains information about sun parameters and user values
434  * view, is camera view vector
435  * s, is distance 
436  * rgb, contains rendered color value for a pixle
437  * */
438 void AtmospherePixleShader(struct SunSky *sunSky, float view[3], float s, float rgb[3])
439 {
440         float costheta;
441         float Phase_1;
442         float Phase_2;
443         float sunColor[3];
444         
445         float E[3];
446         float E1[3];
447         
448         
449         float I[3];
450         float fTemp;
451         float vTemp1[3], vTemp2[3];
452
453         float sunDirection[3];
454         
455         s *= sunSky->atm_DistanceMultiplier;
456         
457         sunDirection[0] = sunSky->toSun[0];
458         sunDirection[1] = sunSky->toSun[1];
459         sunDirection[2] = sunSky->toSun[2];
460
461         costheta = dot_v3v3(view, sunDirection); /* cos(theta) */
462         Phase_1 = 1 + (costheta * costheta); /* Phase_1 */
463
464         VEC3OPF(sunSky->atm_BetaRay, sunSky->atm_BetaRay, *, sunSky->atm_BetaRayMultiplier);
465         VEC3OPF(sunSky->atm_BetaMie, sunSky->atm_BetaMie, *, sunSky->atm_BetaMieMultiplier);
466         VEC3OPV(sunSky->atm_BetaRM, sunSky->atm_BetaRay, +, sunSky->atm_BetaMie);
467
468         /* e^(-(beta_1 + beta_2) * s) = E1 */
469         VEC3OPF(E1, sunSky->atm_BetaRM, *, -s / (float)M_LN2);
470         E1[0] = exp(E1[0]);
471         E1[1] = exp(E1[1]);
472         E1[2] = exp(E1[2]);
473
474         copy_v3_v3(E, E1);
475
476         /* Phase2(theta) = (1-g^2)/(1+g-2g*cos(theta))^(3/2) */
477         fTemp = 1 + sunSky->atm_HGg - 2 * sunSky->atm_HGg * costheta;
478         fTemp = fTemp * sqrtf(fTemp);
479         Phase_2 = (1 - sunSky->atm_HGg * sunSky->atm_HGg) / fTemp;
480         
481         VEC3OPF(vTemp1, sunSky->atm_BetaDashRay, *, Phase_1);
482         VEC3OPF(vTemp2, sunSky->atm_BetaDashMie, *, Phase_2);
483
484         VEC3OPV(vTemp1, vTemp1, +, vTemp2);
485         FOPVEC3(vTemp2, 1.0f, -, E1);
486         VEC3OPV(vTemp1, vTemp1, *, vTemp2);
487
488         FOPVEC3(vTemp2, 1.0f, /, sunSky->atm_BetaRM);
489
490         VEC3OPV(I, vTemp1, *, vTemp2);
491                 
492         VEC3OPF(I, I, *, sunSky->atm_InscatteringMultiplier);
493         VEC3OPF(E, E, *, sunSky->atm_ExtinctionMultiplier);
494
495         /* scale to color sun */
496         ComputeAttenuatedSunlight(sunSky->theta, sunSky->turbidity, sunColor);
497         VEC3OPV(E, E, *, sunColor);
498
499         VEC3OPF(I, I, *, sunSky->atm_SunIntensity);
500
501         VEC3OPV(rgb, rgb, *, E);
502         VEC3OPV(rgb, rgb, +, I);
503 }
504
505 #undef VEC3OPV
506 #undef VEC3OPF
507 #undef FOPVEC3
508
509 /* EOF */