builtin_data/_sky.glsl

   1 struct AtmosphericEvents
   2 {
   3         vec3 rayleigh_scatter;
   4         vec3 mie_scatter;
   5         vec3 mie_absorb;
   6         vec3 ozone_absorb;
   7 };
   8
   9 layout(set=2) uniform Atmosphere
  10 {
  11         AtmosphericEvents events;
  12         float rayleigh_density_decay;
  13         float mie_density_decay;
  14         float ozone_band_center;
  15         float ozone_band_extent;
  16         float planet_radius;
  17         float atmosphere_thickness;
  18         vec3 ground_albedo;
  19         int n_steps;
  20 };
  21
  22 layout(set=2) uniform View
  23 {
  24         float view_height;
  25         vec4 light_color;
  26         vec3 light_dir;
  27 };
  28
  29 struct OpticalPathInfo
  30 {
  31         vec3 optical_depth;
  32         vec3 luminance;
  33 };
  34
  35 const float mie_asymmetry = 0.8;
  36
  37 layout(set=2) uniform sampler2D transmittance_lookup;
  38
  39 vec3 rayleigh_density(vec3 base, float height)
  40 {
  41         return base*exp(height/rayleigh_density_decay);
  42 }
  43
  44 float rayleigh_phase(float cos_theta)
  45 {
  46         return 3.0*(1.0+cos_theta*cos_theta)/(16.0*PI);
  47 }
  48
  49 vec3 mie_density(vec3 base, float height)
  50 {
  51         return base*exp(height/mie_density_decay);
  52 }
  53
  54 float mie_phase(float cos_theta)
  55 {
  56         float g = mie_asymmetry;
  57         float num = (1.0-g*g)*(1.0+cos_theta*cos_theta);
  58         float denom = (2.0+g*g)*pow(1.0+g*g-2.0*g*cos_theta, 1.5);
  59         return 3.0/(8.0*PI)*num/denom;
  60 }
  61
  62 vec3 ozone_density(vec3 base, float height)
  63 {
  64         return base*max(1.0-abs(height-ozone_band_center)/ozone_band_extent, 0.0);
  65 }
  66
  67 AtmosphericEvents calculate_events(float height)
  68 {
  69         AtmosphericEvents ev;
  70         ev.rayleigh_scatter = rayleigh_density(events.rayleigh_scatter, height);
  71         ev.mie_scatter = mie_density(events.mie_scatter, height);
  72         ev.mie_absorb = mie_density(events.mie_absorb, height);
  73         ev.ozone_absorb = ozone_density(events.ozone_absorb, height);
  74         return ev;
  75 }
  76
  77 vec3 total_extinction(AtmosphericEvents ev)
  78 {
  79         return ev.rayleigh_scatter+ev.mie_scatter+ev.mie_absorb+ev.ozone_absorb;
  80 }
  81
  82 float ray_sphere_intersect(vec3 ray_start, vec3 ray_dir, vec3 sphere_center, float sphere_radius)
  83 {
  84         float t = dot(sphere_center-ray_start, ray_dir);
  85         vec3 nearest = ray_start+t*ray_dir-sphere_center;
  86         float d_sq = dot(nearest, nearest);
  87         float r_sq = sphere_radius*sphere_radius;
  88         if(d_sq>r_sq)
  89                 return -1.0;
  90
  91         float offset = sqrt(r_sq-d_sq)*sign(sphere_radius);
  92         return t-offset;
  93 }
  94
  95 #pragma MSP stage(fragment)
  96 OpticalPathInfo raymarch_path(float start_height, vec3 look_dir)
  97 {
  98         float cos_theta = dot(look_dir, light_dir);
  99         float p_rayleigh = rayleigh_phase(cos_theta);
 100         float p_mie = mie_phase(cos_theta);
 101
 102         vec3 planet_center = vec3(0.0, 0.0, -planet_radius);
 103         vec3 pos = vec3(0.0, 0.0, start_height);
 104         vec3 path_luminance = vec3(0.0);
 105         vec3 path_extinction = vec3(0.0);
 106
 107         float ground_t = ray_sphere_intersect(pos, look_dir, planet_center, planet_radius);
 108         float space_t = ray_sphere_intersect(pos, look_dir, planet_center, -(planet_radius+atmosphere_thickness));
 109         float ray_length = (ground_t>0.0 ? ground_t : space_t);
 110         float step_size = ray_length/n_steps;
 111
 112         for(int i=0; i<=n_steps; ++i)
 113         {
 114                 vec3 from_center = pos-planet_center;
 115                 float height = length(from_center);
 116                 float light_z = dot(from_center/height, light_dir);
 117                 height -= planet_radius;
 118
 119                 AtmosphericEvents ev = calculate_events(height);
 120                 vec3 transmittance = exp(-path_extinction);
 121                 float in_ground_t = ray_sphere_intersect(pos, light_dir, planet_center, planet_radius);
 122                 if(in_ground_t<0.0)
 123                 {
 124                         vec3 in_transmittance = texture(transmittance_lookup, vec2(sqrt(height/atmosphere_thickness), light_z)).rgb;
 125                         vec3 in_luminance = (ev.rayleigh_scatter*p_rayleigh+ev.mie_scatter*p_mie)*step_size;
 126                         if(i==n_steps && ground_t>0.0 && light_z>0.0)
 127                                 in_luminance += ground_albedo*light_z/PI;
 128                         path_luminance += transmittance*in_transmittance*in_luminance;
 129                 }
 130
 131                 path_extinction += total_extinction(ev)*step_size;
 132                 pos += look_dir*step_size;
 133         }
 134
 135         return OpticalPathInfo(path_extinction, path_luminance);
 136 }