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