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