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