Better lighting balance for ALS water shader in very low light /night
[fg:toms-fgdata.git] / 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 }