Shaders/water_lightfield_lr.frag

   1 // This shader is mostly an adaptation of the shader found at
   2 //  http://www.bonzaisoftware.com/water_tut.html and its glsl conversion
   3 //  available at http://forum.bonzaisoftware.com/viewthread.php?tid=10
   4 //  © Michael Horsch - 2005
   5 //  Major update and revisions - 2011-10-07
   6 //  © Emilian Huminiuc and Vivian Meazza
   7 // ported to lightfield shading Thorsten Renk 2012
   8
   9 #version 120
  10
  11 uniform sampler2D water_normalmap;
  12 uniform sampler2D water_dudvmap;
  13 uniform sampler2D sea_foam;
  14 uniform sampler2D perlin_normalmap;
  15
  16 uniform sampler3D Noise;
  17
  18 uniform float saturation, Overcast, WindE, WindN;
  19 uniform float osg_SimulationTime;
  20
  21 varying vec4 waterTex1; //moving texcoords
  22 varying vec4 waterTex2; //moving texcoords
  23 varying vec4 waterTex4; //viewts
  24 varying vec3 viewerdir;
  25 varying vec3 lightdir;
  26 //varying vec3 specular_light;
  27 varying vec3 relPos;
  28
  29 varying float earthShade;
  30 varying float yprime_alt;
  31 varying float mie_angle;
  32
  33 uniform    float WaveFreq ;
  34 uniform    float WaveAmp ;
  35 uniform    float WaveSharp ;
  36 uniform    float WaveAngle ;
  37 uniform    float WaveFactor ;
  38 uniform    float WaveDAngle ;
  39 uniform    float normalmap_dds;
  40
  41
  42 uniform float hazeLayerAltitude;
  43 uniform float terminator;
  44 uniform float terrain_alt;
  45 uniform float avisibility;
  46 uniform float visibility;
  47 uniform float overcast;
  48 uniform float scattering;
  49 uniform float ground_scattering;
  50 uniform float cloud_self_shading;
  51 uniform float eye_alt;
  52 uniform float sea_r;
  53 uniform float sea_g;
  54 uniform float sea_b;
  55
  56
  57 vec3 specular_light;
  58
  59 //uniform int wquality_level;
  60
  61 const float terminator_width = 200000.0;
  62 const float EarthRadius = 5800000.0;
  63 ////fog "include" /////
  64 //uniform int fogType;
  65
  66 vec3 fog_Func(vec3 color, int type);
  67 //////////////////////
  68
  69 /////// functions /////////
  70
  71 void rotationmatrix(in float angle, out mat4 rotmat)
  72         {
  73         rotmat = mat4( cos( angle ), -sin( angle ), 0.0, 0.0,
  74                 sin( angle ),  cos( angle ), 0.0, 0.0,
  75                 0.0         ,  0.0         , 1.0, 0.0,
  76                 0.0         ,  0.0         , 0.0, 1.0 );
  77         }
  78
  79 // wave functions ///////////////////////
  80
  81 struct Wave {
  82         float freq;  // 2*PI / wavelength
  83         float amp;   // amplitude
  84         float phase; // speed * 2*PI / wavelength
  85         vec2 dir;
  86         };
  87
  88 Wave wave0 = Wave(1.0, 1.0, 0.5, vec2(0.97, 0.25));
  89 Wave wave1 = Wave(2.0, 0.5, 1.3, vec2(0.97, -0.25));
  90 Wave wave2 = Wave(1.0, 1.0, 0.6, vec2(0.95, -0.3));
  91 Wave wave3 = Wave(2.0, 0.5, 1.4, vec2(0.99, 0.1));
  92
  93
  94
  95
  96 float evaluateWave(in Wave w, vec2 pos, float t)
  97         {
  98         return w.amp * sin( dot(w.dir, pos) * w.freq + t * w.phase);
  99         }
 100
 101 // derivative of wave function
 102 float evaluateWaveDeriv(Wave w, vec2 pos, float t)
 103         {
 104         return w.freq * w.amp * cos( dot(w.dir, pos)*w.freq + t*w.phase);
 105         }
 106
 107 // sharp wave functions
 108 float evaluateWaveSharp(Wave w, vec2 pos, float t, float k)
 109         {
 110         return w.amp * pow(sin( dot(w.dir, pos)*w.freq + t*w.phase)* 0.5 + 0.5 , k);
 111         }
 112
 113 float evaluateWaveDerivSharp(Wave w, vec2 pos, float t, float k)
 114         {
 115         return k*w.freq*w.amp * pow(sin( dot(w.dir, pos)*w.freq + t*w.phase)* 0.5 + 0.5 , k - 1) * cos( dot(w.dir, pos)*w.freq + t*w.phase);
 116         }
 117
 118 void sumWaves(float angle, float dangle, float windScale, float factor, out float ddx, float ddy)
 119         {
 120         mat4 RotationMatrix;
 121         float deriv;
 122         vec4 P = waterTex1 * 1024;
 123
 124         rotationmatrix(radians(angle + dangle * windScale + 0.6 * sin(P.x * factor)), RotationMatrix);
 125         P *= RotationMatrix;
 126
 127         P.y += evaluateWave(wave0, P.xz, osg_SimulationTime);
 128         deriv = evaluateWaveDeriv(wave0, P.xz, osg_SimulationTime );
 129         ddx = deriv * wave0.dir.x;
 130         ddy = deriv * wave0.dir.y;
 131
 132         //P.y += evaluateWave(wave1, P.xz, osg_SimulationTime);
 133         //deriv = evaluateWaveDeriv(wave1, P.xz, osg_SimulationTime);
 134         //ddx += deriv * wave1.dir.x;
 135         //ddy += deriv * wave1.dir.y;
 136
 137         P.y += evaluateWaveSharp(wave2, P.xz, osg_SimulationTime, WaveSharp);
 138         deriv = evaluateWaveDerivSharp(wave2, P.xz, osg_SimulationTime, WaveSharp);
 139         ddx += deriv * wave2.dir.x;
 140         ddy += deriv * wave2.dir.y;
 141
 142         //P.y += evaluateWaveSharp(wave3, P.xz, osg_SimulationTime, WaveSharp);
 143         //deriv = evaluateWaveDerivSharp(wave3, P.xz, osg_SimulationTime, WaveSharp);
 144         //ddx += deriv * wave3.dir.x;
 145         //ddy += deriv * wave3.dir.y;
 146         }
 147
 148
 149 float light_func (in float x, in float a, in float b, in float c, in float d, in float e)
 150 {
 151 x = x - 0.5;
 152
 153 // use the asymptotics to shorten computations
 154 if (x > 30.0) {return e;}
 155 if (x < -15.0) {return 0.0;}
 156
 157 return e / pow((1.0 + a * exp(-b * (x-c)) ),(1.0/d));
 158 }
 159
 160 // this determines how light is attenuated in the distance
 161 // physically this should be exp(-arg) but for technical reasons we use a sharper cutoff
 162 // for distance > visibility
 163
 164 float fog_func (in float targ)
 165 {
 166
 167
 168 float fade_mix;
 169
 170 // for large altitude > 30 km, we switch to some component of quadratic distance fading to
 171 // create the illusion of improved visibility range
 172
 173 targ = 1.25 * targ; // need to sync with the distance to which terrain is drawn
 174
 175
 176 if (eye_alt < 30000.0)
 177         {return exp(-targ - targ * targ * targ * targ);}
 178 else if (eye_alt < 50000.0)
 179         {
 180         fade_mix = (eye_alt - 30000.0)/20000.0;
 181         return fade_mix * exp(-targ*targ - pow(targ,4.0)) + (1.0 - fade_mix) * exp(-targ - pow(targ,4.0));
 182         }
 183 else
 184         {
 185         return exp(- targ * targ - pow(targ,4.0));
 186         }
 187
 188 }
 189
 190 void main(void)
 191         {
 192
 193
 194         vec3 shadedFogColor = vec3(0.65, 0.67, 0.78);
 195         float effective_scattering = min(scattering, cloud_self_shading);
 196
 197         float dist = length(relPos);
 198         const vec4 sca = vec4(0.005, 0.005, 0.005, 0.005);
 199         const vec4 sca2 = vec4(0.02, 0.02, 0.02, 0.02);
 200         const vec4 tscale = vec4(0.25, 0.25, 0.25, 0.25);
 201
 202         mat4 RotationMatrix;
 203
 204         // compute direction to viewer
 205         vec3 E = normalize(viewerdir);
 206
 207         // compute direction to light source
 208         vec3 L = lightdir; // normalize(lightdir);
 209
 210         // half vector
 211         vec3 Hv = normalize(L + E);
 212
 213         //vec3 Normal = normalize(normal);
 214         vec3 Normal = vec3 (0.0, 0.0, 1.0);
 215
 216         const float water_shininess = 240.0;
 217
 218         // approximate cloud cover
 219         //float cover = 0.0;
 220         //bool Status = true;
 221
 222         float windEffect = sqrt( WindE*WindE + WindN*WindN ) * 0.6;                             //wind speed in kt
 223         float windScale =  15.0/(3.0 + windEffect);                                                                                                     //wave scale
 224         float windEffect_low = 0.3 + 0.7 * smoothstep(0.0, 5.0, windEffect);                                    //low windspeed wave filter
 225         float waveRoughness = 0.01 + smoothstep(0.0, 40.0, windEffect);                                         //wave roughness filter
 226
 227         float mixFactor = 0.2 + 0.02 * smoothstep(0.0, 50.0, windEffect);
 228         //mixFactor = 0.2;
 229         mixFactor = clamp(mixFactor, 0.3, 0.8);
 230
 231         // there's no need to do wave patterns or foam for pixels which are so far away that we can't actually see them
 232         // we only need detail in the near zone or where the sun reflection is
 233
 234         int detail_flag;
 235         if ((dist > 15000.0) && (dot(normalize(vec3 (lightdir.x, lightdir.y, 0.0) ), normalize(relPos)) < 0.7 ))  {detail_flag = 0;}
 236         else {detail_flag = 1;}
 237
 238         //detail_flag = 1;
 239
 240         // sine waves
 241         float ddx, ddx1, ddx2, ddx3, ddy, ddy1, ddy2, ddy3;
 242         float angle;
 243
 244
 245         ddx = 0.0, ddy = 0.0;
 246         ddx1 = 0.0, ddy1 = 0.0;
 247         ddx2 = 0.0, ddy2 = 0.0;
 248         ddx3 = 0.0, ddy3 = 0.0;
 249         if (detail_flag == 1)
 250         {
 251         angle = 0.0;
 252
 253         wave0.freq = WaveFreq ;
 254         wave0.amp = WaveAmp;
 255         wave0.dir =  vec2 (0.0, 1.0); //vec2(cos(radians(angle)), sin(radians(angle)));
 256
 257         angle -= 45;
 258         wave1.freq = WaveFreq * 2.0 ;
 259         wave1.amp = WaveAmp * 1.25;
 260         wave1.dir =  vec2(0.70710, -0.7071); //vec2(cos(radians(angle)), sin(radians(angle)));
 261
 262         angle += 30;
 263         wave2.freq = WaveFreq * 3.5;
 264         wave2.amp = WaveAmp * 0.75;
 265         wave2.dir =  vec2(0.96592, -0.2588);// vec2(cos(radians(angle)), sin(radians(angle)));
 266
 267         angle -= 50;
 268         wave3.freq = WaveFreq * 3.0 ;
 269         wave3.amp = WaveAmp * 0.75;
 270         wave3.dir =  vec2(0.42261, -0.9063); //vec2(cos(radians(angle)), sin(radians(angle)));
 271
 272         // sum waves
 273
 274         sumWaves(WaveAngle, -1.5, windScale, WaveFactor, ddx, ddy);
 275         sumWaves(WaveAngle, 1.5, windScale, WaveFactor, ddx1, ddy1);
 276
 277         //reset the waves
 278         angle = 0.0;
 279         float waveamp = WaveAmp * 0.75;
 280
 281         wave0.freq = WaveFreq ;
 282         wave0.amp = waveamp;
 283         wave0.dir =  vec2 (0.0, 1.0); //vec2(cos(radians(angle)), sin(radians(angle)));
 284
 285         angle -= 20;
 286         wave1.freq = WaveFreq * 2.0 ;
 287         wave1.amp = waveamp * 1.25;
 288         wave1.dir =  vec2(0.93969, -0.34202);// vec2(cos(radians(angle)), sin(radians(angle)));
 289
 290         angle += 35;
 291         wave2.freq = WaveFreq * 3.5;
 292         wave2.amp = waveamp * 0.75;
 293         wave2.dir =  vec2(0.965925, 0.25881);  //vec2(cos(radians(angle)), sin(radians(angle)));
 294
 295         angle -= 45;
 296         wave3.freq = WaveFreq * 3.0 ;
 297         wave3.amp = waveamp * 0.75;
 298         wave3.dir =  vec2(0.866025, -0.5); //vec2(cos(radians(angle)), sin(radians(angle)));
 299
 300
 301         //sumWaves(WaveAngle + WaveDAngle, -1.5, windScale, WaveFactor, ddx2, ddy2);
 302         //sumWaves(WaveAngle + WaveDAngle, 1.5, windScale, WaveFactor, ddx3, ddy3);
 303
 304         }
 305         // end sine stuff
 306
 307         //cover = 5.0 * smoothstep(0.6, 1.0, scattering);
 308         //cover = 5.0 * ground_scattering;
 309
 310         vec4 viewt = normalize(waterTex4);
 311
 312         vec4 disdis = texture2D(water_dudvmap, vec2(waterTex2 * tscale)* windScale) * 2.0 - 1.0;
 313
 314         vec4 vNorm;
 315
 316
 317         //normalmaps
 318         vec4 nmap   = texture2D(water_normalmap, vec2(waterTex1 + disdis * sca2) * windScale) * 2.0 - 1.0;
 319         vec4 nmap1  = texture2D(perlin_normalmap, vec2(waterTex1 + disdis * sca2) * windScale) * 2.0 - 1.0;
 320
 321         rotationmatrix(radians(3.0 * sin(osg_SimulationTime * 0.0075)), RotationMatrix);
 322         nmap  += texture2D(water_normalmap, vec2(waterTex2 * RotationMatrix * tscale) * windScale) * 2.0 - 1.0;
 323         //nmap1 += texture2D(perlin_normalmap, vec2(waterTex2 * RotationMatrix * tscale) * windScale) * 2.0 - 1.0;
 324
 325         nmap  *= windEffect_low;
 326         nmap1 *= windEffect_low;
 327
 328         // mix water and noise, modulated by factor
 329         vNorm = normalize(mix(nmap, nmap1, mixFactor) * waveRoughness);
 330         vNorm.r += ddx + ddx1 + ddx2 + ddx3;
 331
 332
 333         if (normalmap_dds > 0)
 334                 {vNorm = -vNorm;}               //dds fix
 335
 336
 337         //load reflection
 338
 339         vec4 refl ;
 340
 341         refl.r = sea_r;
 342         refl.g = sea_g;
 343         refl.b = sea_b;
 344         refl.a = 1.0;
 345
 346
 347         float intensity;
 348         // de-saturate for reduced light
 349         refl.rgb = mix(refl.rgb,  vec3 (0.248, 0.248, 0.248), 1.0 - smoothstep(0.1, 0.8, ground_scattering));
 350
 351         // de-saturate light for overcast haze
 352         intensity = length(refl.rgb);
 353         refl.rgb = mix(refl.rgb,  intensity * vec3 (1.0, 1.0, 1.0), 0.5 * smoothstep(0.1, 0.9, overcast));
 354
 355         vec3 N;
 356
 357
 358
 359
 360         vec3 N0 = vec3(texture2D(water_normalmap, vec2(waterTex1 + disdis * sca2) * windScale) * 2.0 - 1.0);
 361         vec3 N1 = vec3(texture2D(perlin_normalmap, vec2(waterTex1 + disdis * sca) * windScale) * 2.0 - 1.0);
 362
 363         N0 += vec3(texture2D(water_normalmap, vec2(waterTex1 * tscale) * windScale) * 2.0 - 1.0);
 364         N1 += vec3(texture2D(perlin_normalmap, vec2(waterTex2 * tscale) * windScale) * 2.0 - 1.0);
 365
 366
 367
 368         rotationmatrix(radians(2.0 * sin(osg_SimulationTime * 0.005)), RotationMatrix);
 369         N0 += vec3(texture2D(water_normalmap, vec2(waterTex2 * RotationMatrix * (tscale + sca2)) * windScale) * 2.0 - 1.0);
 370         N1 += vec3(texture2D(perlin_normalmap, vec2(waterTex2 * RotationMatrix * (tscale + sca2)) * windScale) * 2.0 - 1.0);
 371
 372         rotationmatrix(radians(-4.0 * sin(osg_SimulationTime * 0.003)), RotationMatrix);
 373         N0 += vec3(texture2D(water_normalmap, vec2(waterTex1 * RotationMatrix + disdis * sca2) * windScale) * 2.0 - 1.0);
 374         N1 += vec3(texture2D(perlin_normalmap, vec2(waterTex1 * RotationMatrix + disdis * sca) * windScale) * 2.0 - 1.0);
 375
 376
 377         N0 *= windEffect_low;
 378         N1 *= windEffect_low;
 379
 380         N0.r += (ddx + ddx1 + ddx2 + ddx3);
 381         N0.g += (ddy + ddy1 + ddy2 + ddy3);
 382
 383         N = normalize(mix(Normal + N0, Normal + N1, mixFactor) * waveRoughness);
 384
 385          if (normalmap_dds > 0)
 386                 {N = -N;} //dds fix
 387
 388
 389
 390
 391
 392        specular_light = gl_Color.rgb;
 393
 394
 395         vec3 specular_color = vec3(specular_light)
 396                 * pow(max(0.0, dot(N, Hv)), water_shininess) * 6.0;
 397         vec4 specular = vec4(specular_color, 0.5);
 398
 399         specular = specular * saturation * 0.3  * earthShade  ;
 400
 401         //calculate fresnel
 402         vec4 invfres = vec4( dot(vNorm, viewt) );
 403         vec4 fres = vec4(1.0) + invfres;
 404         refl *= fres;
 405
 406
 407
 408         vec4 ambient_light;
 409         //intensity = length(specular_light.rgb);
 410         ambient_light.rgb = max(specular_light.rgb, vec3(0.05, 0.05, 0.05));
 411         //ambient_light.rgb = max(intensity * normalize(vec3 (0.33, 0.4, 0.5)), vec3 (0.1,0.1,0.1));
 412         ambient_light.a = 1.0;
 413
 414
 415         vec4 finalColor;
 416
 417
 418
 419         finalColor = refl + specular * smoothstep(0.3, 0.6, ground_scattering);
 420
 421         //add foam
 422         vec4 foam_texel = texture2D(sea_foam, vec2(waterTex2 * tscale) * 25.0);
 423
 424         if (dist < 10000.0)
 425         {
 426         float foamSlope = 0.10 + 0.1 * windScale;
 427         float waveSlope = N.g;
 428
 429         if (windEffect >= 8.0)
 430                 if (waveSlope >= foamSlope){
 431                         finalColor = mix(finalColor, max(finalColor, finalColor + foam_texel), smoothstep(0.01, 0.50, N.g));
 432                         }
 433         }
 434
 435
 436
 437                 finalColor *= ambient_light;
 438
 439
 440
 441 // here comes the terrain haze model
 442
 443
 444 float delta_z = hazeLayerAltitude - eye_alt;
 445
 446
 447
 448 if (dist > 40.0)
 449 {
 450
 451
 452 float transmission;
 453 float vAltitude;
 454 float delta_zv;
 455 float H;
 456 float distance_in_layer;
 457 float transmission_arg;
 458
 459
 460 // angle with horizon
 461 float ct = dot(vec3(0.0, 0.0, 1.0), relPos)/dist;
 462
 463
 464 // we solve the geometry what part of the light path is attenuated normally and what is through the haze layer
 465
 466 if (delta_z > 0.0) // we're inside the layer
 467         {
 468         if (ct < 0.0) // we look down
 469                 {
 470                 distance_in_layer = dist;
 471                 vAltitude = min(distance_in_layer,min(visibility,avisibility)) * ct;
 472                 delta_zv = delta_z - vAltitude;
 473                 }
 474         else    // we may look through upper layer edge
 475                 {
 476                 H = dist * ct;
 477                 if (H > delta_z) {distance_in_layer = dist/H * delta_z;}
 478                 else {distance_in_layer = dist;}
 479                 vAltitude = min(distance_in_layer,visibility) * ct;
 480                 delta_zv = delta_z - vAltitude;
 481                 }
 482         }
 483   else // we see the layer from above, delta_z < 0.0
 484         {
 485         H = dist * -ct;
 486         if (H  < (-delta_z)) // we don't see into the layer at all, aloft visibility is the only fading
 487                 {
 488                 distance_in_layer = 0.0;
 489                 delta_zv = 0.0;
 490                 }
 491         else
 492                 {
 493                 vAltitude = H + delta_z;
 494                 distance_in_layer = vAltitude/H * dist;
 495                 vAltitude = min(distance_in_layer,visibility) * (-ct);
 496                 delta_zv = vAltitude;
 497                 }
 498         }
 499
 500
 501 // ground haze cannot be thinner than aloft visibility in the model,
 502 // so we need to use aloft visibility otherwise
 503
 504
 505 transmission_arg = (dist-distance_in_layer)/avisibility;
 506
 507
 508 float eqColorFactor;
 509
 510
 511 if (visibility < avisibility)
 512         {
 513         transmission_arg = transmission_arg + (distance_in_layer/visibility);
 514         // this combines the Weber-Fechner intensity
 515         eqColorFactor = 1.0 - 0.1 * delta_zv/visibility - (1.0 -effective_scattering);
 516
 517         }
 518 else
 519         {
 520         transmission_arg = transmission_arg + (distance_in_layer/avisibility);
 521         // this combines the Weber-Fechner intensity
 522         eqColorFactor = 1.0 - 0.1 * delta_zv/avisibility - (1.0 -effective_scattering);
 523         }
 524
 525
 526 transmission =  fog_func(transmission_arg);
 527
 528 // there's always residual intensity, we should never be driven to zero
 529 if (eqColorFactor < 0.2) eqColorFactor = 0.2;
 530
 531
 532 float lightArg = (terminator-yprime_alt)/100000.0;
 533
 534 vec3 hazeColor;
 535 hazeColor.b = light_func(lightArg, 1.330e-05, 0.264, 2.527, 1.08e-05, 1.0);
 536 hazeColor.g = light_func(lightArg, 3.931e-06, 0.264, 3.827, 7.93e-06, 1.0);
 537 hazeColor.r = light_func(lightArg, 8.305e-06, 0.161, 3.827, 3.04e-05, 1.0);
 538
 539 // now dim the light for haze
 540 float eShade = 0.9 * smoothstep(terminator_width+ terminator, -terminator_width + terminator, yprime_alt) + 0.1;
 541
 542 // Mie-like factor
 543
 544 if (lightArg < 10.0)
 545         {intensity = length(hazeColor);
 546         float mie_magnitude = 0.5 * smoothstep(350000.0, 150000.0, terminator-sqrt(2.0 * EarthRadius * terrain_alt));
 547         hazeColor = intensity * ((1.0 - mie_magnitude) + mie_magnitude * mie_angle) * normalize(mix(hazeColor,  vec3 (0.5, 0.58, 0.65), mie_magnitude * (0.5 - 0.5 * mie_angle)) );
 548         }
 549
 550 // high altitude desaturation of the haze color
 551
 552 intensity = length(hazeColor);
 553
 554
 555 if (intensity > 0.0) // this needs to be a condition, because otherwise hazeColor doesn't come out correctly
 556         {
 557         hazeColor = intensity * normalize (mix(hazeColor, intensity * vec3 (1.0,1.0,1.0), 0.7* smoothstep(5000.0, 50000.0, eye_alt)));
 558
 559         // blue hue of haze
 560
 561         hazeColor.x = hazeColor.x * 0.83;
 562         hazeColor.y = hazeColor.y * 0.9;
 563
 564
 565         // additional blue in indirect light
 566         float fade_out = max(0.65 - 0.3 *overcast, 0.45);
 567         intensity = length(hazeColor);
 568         hazeColor = intensity * normalize(mix(hazeColor,  1.5* shadedFogColor, 1.0 -smoothstep(0.25, fade_out,eShade) ));
 569
 570         // change haze color to blue hue for strong fogging
 571         hazeColor = intensity * normalize(mix(hazeColor,  shadedFogColor, (1.0-smoothstep(0.5,0.9,eqColorFactor))));
 572         }
 573
 574
 575
 576         finalColor.rgb = mix(eqColorFactor * hazeColor * eShade, finalColor.rgb,transmission);
 577
 578
 579         }
 580         gl_FragColor = finalColor;
 581
 582 }