$FG_GLSL_VERSION layout(location = 0) out vec4 fragColor; in vec2 texcoord; in vec2 raw_texcoord; // QUAD_TEXCOORD_RAW in vec3 w_pos; uniform sampler3D detailed_tex; uniform sampler3D rough_tex; uniform sampler3D shade_tex; uniform sampler2D depth_tex; uniform sampler3D cloud_noise_tex; uniform sampler1D wind_offset_tex; uniform vec3 fg_SunDirectionWorld; uniform vec3 fg_CameraPositionCart; uniform float fg_EarthRadius; uniform mat4 fg_CameraZUpMatrix; uniform mat4 fg_ViewMatrix[FG_NUM_VIEWS]; FG_VIEW_GLOBAL uniform mat4 fg_ViewMatrixInverse[FG_NUM_VIEWS]; uniform vec4 fg_Viewport[FG_NUM_VIEWS]; uniform vec3 cloud_field_center; uniform bool cloud_field_repeating; uniform bool cloud_field_mirror_u; uniform bool cloud_field_mirror_v; uniform float cloud_base_z_norm; uniform float voxel_resolution_m; uniform int voxel_field_width; uniform int voxel_field_height; uniform float active_voxel_field_height_norm; uniform int rough_field_factor; uniform int rough_size_factor; uniform float voxel_optical_depth; uniform vec4 ground_albedo; // math.glsl float M_1_4PI(); float M_1_PI(); float safe_sqrt(float x); mat3 getRotateRollPitchYaw(vec3 rpw); // exposure.glsl vec3 apply_exposure(vec3 color); // logarithmic_depth.glsl float logdepth_prepare_vs_depth(float z); float logdepth_decode(float z); float logdepth_encode(float z); // sun.glsl vec3 get_sun_radiance_sea_level(); // pos_from_depth.glsl vec3 get_view_space_from_vs_depth(vec2 uv, float vs_depth); // aerial_perspective.glsl vec4 get_aerial_perspective(vec2 raw_coord, vec3 P); vec4 diff_aerial_perspective(vec4 apFar, vec4 apNear); const float MIN_DIST = 0.0001; const float MAX_DIST = 4.0; const float EPSILON = 0.000001; const float NOISE_SCALE = 48.0; const float NOISE_SKIP_THRESHOLD = 0.95; const int MAX_LIGHT_STEPS = 5; int MAX_MARCHING_STEPS = 4 * voxel_field_width; float IN_CLOUD_STEP_SIZE = 0.1f / float(voxel_field_width); float IN_CLOUD_SUN_RAY_STEP_SIZE = 1.0 / float(voxel_field_width); float DENSITY_SCALE = IN_CLOUD_STEP_SIZE * float(voxel_field_height) / 0.1; float VOXEL_FIELD_WIDTH_M = float(voxel_field_width * voxel_resolution_m); float VOXEL_FIELD_HEIGHT_M = float(voxel_field_height * voxel_resolution_m); float CURVATURE_DENOM = 2.0 * fg_EarthRadius * VOXEL_FIELD_HEIGHT_M; float MAX_CLOUD_DENSITY = 5.0; // Scaling factor to balance between direct scattering and the ambient/multiscattering/ground-bounce vec3 PIPELINE_RADIANCE_SCALE = vec3(0.1); // Direct light desaturates very slightly, simulating the eye partially // adapting to the overall warmth of the direct sunlight. float DIRECT_DESATURATION = 0.15; // Ground-bounced light only meaningfully reaches a few hundred metres // above the reflecting surface. However, our origin is at 0 elevation, // so we need to consider the mean elevation of the Earth surface (244m) float GROUND_BOUNCE_MAX_HEIGHT_M = 1244.0; // Scaling factor to account for the voxel space not being a cube. vec3 VOXEL_SCALE = vec3(1.0, 1.0, float(voxel_field_width) / float(voxel_field_height)); vec3 EYE_SCALE = vec3(VOXEL_FIELD_WIDTH_M, VOXEL_FIELD_WIDTH_M, VOXEL_FIELD_HEIGHT_M); // Combined ratio between the rough field's physical size and the detailed field's. Both the // voxel-count ratio (rough_field_factor) and the per-voxel physical-size ratio // (rough_size_factor) contribute; using either alone under-scales the rough field's UV frame. float ROUGH_FIELD_SCALE = float(rough_field_factor) * float(rough_size_factor); float hash12(vec2 p) { vec3 p3 = fract(vec3(p.xyx) * 0.1031); p3 += dot(p3, p3.yzx + 33.33); return fract((p3.x + p3.y) * p3.z); } // // Function to erode a value given an erosion amount. A simplified version of SetRange. // float ValueErosion(float inValue, float inOldMin) { float old_min_max_range = (1.0 - inOldMin); return clamp((inValue - inOldMin) / old_min_max_range, 0.0, 1.0); } // Get a position in a curved reference frame to make the voxel space curve // around the earth. We achieve this by adjusting the position upwards by the sagitta // appromixation of the earth's curvature. // This should be used for all lookups into the voxel spaces. // eye is the eyepoint in voxel coordinates, and p is the position to be adjusted // in voxel space. vec3 getCurvedP(vec3 eye, vec3 p) { vec2 l = p.xy - eye.xy; float horiz_dist_m_2 = dot(l, l) * VOXEL_FIELD_WIDTH_M * VOXEL_FIELD_WIDTH_M; float sagitta_norm = horiz_dist_m_2 / CURVATURE_DENOM; return vec3(p.xy, p.z + sagitta_norm); } // HenyeyGreenstein forward phase scattering float HenyeyGreenstein(float inCosAngle, float inG) { float num = 1.0 - inG * inG; float denom = 1.0 + inG * inG - 2.0 * inG * inCosAngle; float rsqrt_denom = 1.0 / safe_sqrt(denom); return num * rsqrt_denom * rsqrt_denom * rsqrt_denom * M_1_4PI(); } // Calculate any silver lining. This is direct scattering right at the edge of the clouds looking // directly at the sun. // Peaks for partially transparent clouds, falls to zero for fully opaque float getSilverLining(const float light_absorption, const float cosTheta) { const float MIN_DENSITY = 0.1; // The minimum density of cloud we will generate a lining for. Anything less than this is assumed transparent. const float SUN_DIRECTIONALITY = 40.0; // How directional the silver lining is to the sun. Higher values limit it closer to the sun. float cloudPresence = clamp((light_absorption - MIN_DENSITY) * 4.0, 0.0, 1.0); // ramps up quickly float cloudEdge = cloudPresence * (1.0 - light_absorption); // falls to opaque float rim = pow(clamp((1.0 + cosTheta) * 0.5, 0.0, 1.0), SUN_DIRECTIONALITY); return rim * 0.8 * cloudEdge; } /** * Get cloud information for a given sample point. This could be * outside the voxel space entirely, or inside a repeating or non-repeating * space. */ vec4 getCloud(vec3 cameraEye, vec3 samplePoint, vec3 dir) { // If outside the voxel space then calculate an SDF directly. This is basically the // z coordinate plus a little bit to ensure it ends up within the voxel space, vec3 curvedPoint = getCurvedP(cameraEye, samplePoint); if (curvedPoint.z < 0.0) return vec4(0.0, 0.0, 0.0, -curvedPoint.z / length(dir)); if (curvedPoint.z > active_voxel_field_height_norm) return vec4(0.0, 0.0, 0.0, (curvedPoint.z - active_voxel_field_height_norm) / length(dir)); if (cloud_field_repeating) { // Within a repeating cloud field we need to adjust the UV coordinates manually // to mirror through the voxel space. if (cloud_field_mirror_u) curvedPoint.x = (1.0 - curvedPoint.x); if (cloud_field_mirror_v) curvedPoint.y = (1.0 - curvedPoint.y); return texture(detailed_tex, curvedPoint); } // Non-repeating field: a detailed voxel space plus a much larger, coarser rough voxel space // beyond it, sharing the same center. samplePoint/curvedPoint are already normalized to the // detailed field's own UV space (see EYE_SCALE above), so "in the detailed field" is simply // "inside its own [0,1] range". if (abs(samplePoint.x - 0.5) < 0.5 && abs(samplePoint.y - 0.5) < 0.5) { // Within the detailed field's own footprint - sample it directly, as the repeating // branch above does (mirroring aside, which doesn't apply to a finite field). return texture(detailed_tex, curvedPoint); } else { // Outside the detailed field - reproject into the rough field's own, much larger UV // frame before sampling. Both fields share the same center and Z extent, so this is a // linear rescale around 0.5 on X/Y only. vec3 roughPoint = vec3(0.5 + (curvedPoint.x - 0.5) / ROUGH_FIELD_SCALE, 0.5 + (curvedPoint.y - 0.5) / ROUGH_FIELD_SCALE, curvedPoint.z); vec4 roughSample = texture(rough_tex, roughPoint); // The rough texture's SDF (alpha) is baked as a fraction of its own (much larger) UV // width, but the raymarch's `distance` accumulator is in detailed-field-normalized // units - rescale before returning (sign is unaffected either way). roughSample.a *= ROUGH_FIELD_SCALE; return roughSample; } } /** * Calculate the density of a given samplePoint */ float calculateDensity(vec3 samplePoint, vec3 dir, vec4 cloud) { float cloudDimension = cloud.x; float cloudType = cloud.y; float cloudDensity = cloud.z; if (cloudDimension > NOISE_SKIP_THRESHOLD) { // Skip expensive noise if we're deep inside cloud, but apply the sharpening. float uprezzed_density = cloudDensity * cloudDimension * voxel_optical_depth; float powered_density_scale = pow(clamp(cloudDensity, 0.0, 1.0), 4.0); return pow(uprezzed_density, mix(0.3, 0.6, max(EPSILON, powered_density_scale))); } vec4 noise = texture(cloud_noise_tex, samplePoint * NOISE_SCALE / VOXEL_SCALE); float wispy_noise = mix(noise.r, noise.g, cloudDimension); // Define billowy noise float billowy_type_gradient = sqrt(sqrt(cloudDimension)); float billowy_noise = mix(noise.b * 0.3, noise.a * 0.3, billowy_type_gradient); // Define Noise composite: blend to wispy depending on the cloud type float noise_composite = mix(wispy_noise, billowy_noise, cloudType); // Ultra-HF noise. //float hhf_wisps = 1.0 - pow(abs(abs(noise.g * 2.0 - 1.0) * 2.0 - 1.0), 4.0); //float hhf_billows = pow(abs(abs(noise.a * 2.0 - 1.0) * 2.0 - 1.0), 2.0); //float hhf_noise_composite = mix(hhf_wisps, hhf_billows, cloudType); //noise_composite = mix(hhf_noise_composite, noise_composite, depth * 10.0); float uprezzed_density = noise_composite; // Composite Noises and use as a Value Erosion float dim = cloudDimension; dim = smoothstep(0.0, 1.0, dim); uprezzed_density = ValueErosion(dim, noise_composite); // Apply Density Scale Data to Result uprezzed_density *= cloudDensity * voxel_optical_depth; // Sharpen result float powered_density_scale = pow(clamp(cloudDensity, 0.0, 1.0), 4.0); uprezzed_density = pow(uprezzed_density, mix(0.3, 0.6, max(EPSILON, powered_density_scale))); return uprezzed_density; } float getRayDensity(vec3 cameraEye, vec3 eye, vec3 marchingDirection, float start, float end) { float density = 0.0; float distance = start; vec3 p = eye + distance * marchingDirection; // Sample shade_tex at the ray origin (eye), not at start offset vec3 windOffset = texture(wind_offset_tex, eye.z).xyz; float storedDensity = texture(shade_tex, getCurvedP(cameraEye, eye + windOffset)).r; // If already heavily shadowed, skip local march if (storedDensity > MAX_CLOUD_DENSITY) { return storedDensity; } for (int i = 0; i < MAX_LIGHT_STEPS; i++) { p = eye + distance * marchingDirection; windOffset = texture(wind_offset_tex, p.z).xyz; if (p.z > active_voxel_field_height_norm) return density; // Reached the top of the actual cloud space vec4 t = getCloud(cameraEye, p + windOffset, marchingDirection); if (t.a < EPSILON) { // Inside a cloud, so add density density += calculateDensity(p + windOffset, marchingDirection, t); distance += IN_CLOUD_SUN_RAY_STEP_SIZE; } else { distance += t.a; } if (density > MAX_CLOUD_DENSITY) { // Reached maximum density, so no point in marching further. return density; } if (distance > end) { // Reached the end of the raymarch. return density; } } // At this point just look up the pre-calculated shade texture p = eye + (distance + IN_CLOUD_SUN_RAY_STEP_SIZE) * marchingDirection; // The shade texture contains the density towards the sun. return density + texture(shade_tex, getCurvedP(cameraEye, p)).r; } /** * Absolute value of the SDF value indicates the distance to the surface. * Sign indicates whether the point is inside or outside the surface, * negative indicating inside. */ struct sample_information { float sdf; float density; float direct_scattering; float ambient_scattering; }; struct ray_data { float light_absorption; vec3 intensity; float first_hit; }; ray_data cloudRayMarch(vec3 eye, vec3 marchingDirection, float start, float end, float zscaleFactor) { ray_data lreturn; lreturn.light_absorption = 0.0; lreturn.intensity = vec3(0.0); lreturn.first_hit = -1.0; float dirZ = marchingDirection.z; // ------------------------------------------------------- // Intersect ray with Z slab [0, active_voxel_field_height_norm] // ------------------------------------------------------- float tEnter = start; float tExit = end; // For rays above the slab pointing upward, nothing to render. if (dirZ > EPSILON && eye.z > active_voxel_field_height_norm) return lreturn; // For rays with meaningful dirZ, clamp tExit to avoid marching forever above the slab. if (dirZ > EPSILON && eye.z < active_voxel_field_height_norm) { float tTop = (active_voxel_field_height_norm - eye.z) / dirZ; tExit = min(tExit, tTop); } // Fast-forward tEnter to approximate curved cloud base entry point // to avoid wasting march steps on empty space below the slab. float h_norm_2 = dot(marchingDirection.xy, marchingDirection.xy); if (h_norm_2 > EPSILON && eye.z < 0.0) { float curve_coeff = h_norm_2 * VOXEL_FIELD_WIDTH_M * VOXEL_FIELD_WIDTH_M / CURVATURE_DENOM; float tCurvedBase = safe_sqrt(-eye.z / curve_coeff); // Also compute the flat-geometry entry (via dirZ climbing to z=0) float tFlatBase = (dirZ > EPSILON) ? (-eye.z / dirZ) : tExit; // Use whichever gets us to the cloud base sooner float tFastForward = min(tCurvedBase, tFlatBase) * 0.9; tEnter = max(tEnter, tFastForward); if (tEnter >= tExit) return lreturn; } // ------------------------------------------------------- // Standard raymarch, now guaranteed inside slab. // Calculate some dynamic limits and step sizes. // ------------------------------------------------------- // Stable per-pixel jitter (world-space stable) vec3 entryPoint = eye + tEnter * marchingDirection; float jitter = hash12(entryPoint.xy); float stepSize = IN_CLOUD_STEP_SIZE; // Offset initial march position slightly float distance = tEnter + jitter * stepSize; float maxTravel = tExit - tEnter; int dynamicMaxSteps = int(maxTravel / IN_CLOUD_STEP_SIZE) + 1; dynamicMaxSteps = min(dynamicMaxSteps, MAX_MARCHING_STEPS); float cachedSunDensity = -1.0; float cachedDensity = -1.0; float multiScatterTerm = 0.0; float lastAir = 0.0; float lastCloud = 0.0; vec4 lastAP = vec4(0.0, 0.0, 0.0, 1.0); vec3 sunRadiance = get_sun_radiance_sea_level(); vec3 sunLuminance = vec3(dot(sunRadiance, vec3(0.2126, 0.7152, 0.0722))); // For direct scatter: keep almost the full warm colour of sunRadiance, // desaturating very slightly to simulate the eye adjusting. vec3 directColor = mix(sunRadiance, sunLuminance, DIRECT_DESATURATION); // For ambient scale: use magnitude only, colour comes from skyAmbientColor. // Clamped to 1.0: at any normal sunRadiance magnitude (~100-300) this is // equivalent to dropping the multiplier entirely; it only lets the real // value through below 1.0, preserving ambient dimming near twilight. float sunMagnitude = min(length(sunRadiance),1.0); // Multi-scatter desaturates toward white so we don't use raw sunRadiance colour vec3 multiScatterColor = mix(sunRadiance, sunLuminance, 0.90); // nearly grey vec3 groundBounce = vec3(0.0); // fg_SunDirectionWorld points FROM world TOWARD sun mat3 zup = mat3(fg_CameraZUpMatrix); vec3 toSunDir = zup * fg_SunDirectionWorld; float dayFactor = smoothstep(-0.2, 0.1, toSunDir.z); // cosTheta > 0 = looking toward sun (forward scatter) // cosTheta < 0 = looking away from sun (back scatter) // sample -> camera vec3 V = normalize(-marchingDirection); float cosTheta = clamp(dot(V, -toSunDir), -1.0, 1.0); // Prevent over-warming at low sun angles by blending toward neutral blue vec3 skyAmbientColor = mix(vec3(0.4, 0.5, 0.7), vec3(0.7, 0.8, 1.0), dayFactor); for (int i = 0; i < dynamicMaxSteps; i++) { if (distance >= tExit) break; vec3 p = eye + distance * marchingDirection; vec3 windOffset = texture(wind_offset_tex, getCurvedP(eye, p).z ).xyz; sample_information s; vec4 cloud = getCloud(eye, p + windOffset, marchingDirection); s.sdf = cloud.a; s.density = 0.0; s.direct_scattering = 0.0; s.ambient_scattering = 0.0; if (s.sdf < EPSILON) { // We're inside a cloud, so calculate the density etc. s.density = calculateDensity(p + windOffset, marchingDirection, cloud) * min(1.0, (tExit - distance) / IN_CLOUD_STEP_SIZE); s.sdf = IN_CLOUD_STEP_SIZE; // March from sample toward sun vec3 sunDirMarch = normalize(toSunDir / EYE_SCALE * VOXEL_SCALE); float densityToSun; // Cache density information and save on expensive ray cast when the cloud density hasn't // changed much. This is a major performance bottleneck. if (cachedSunDensity < 0.0 || abs(s.density - cachedDensity) > 0.02) { // World-space stable jitter - use the sample position in metres // hash12 needs a 2D input so combine two world-space axes vec3 worldP = p * EYE_SCALE; // convert from voxel-norm to metres float jitter = texture(cloud_noise_tex, worldP * 0.0001).r * IN_CLOUD_SUN_RAY_STEP_SIZE; densityToSun = getRayDensity(eye, p + sunDirMarch * jitter, sunDirMarch, IN_CLOUD_SUN_RAY_STEP_SIZE, 1.0); cachedSunDensity = densityToSun; cachedDensity = s.density; } else { densityToSun = cachedSunDensity; } float transmittance = exp(-densityToSun); float phaseForward = HenyeyGreenstein(cosTheta, 0.5); float phaseBackward = HenyeyGreenstein(cosTheta, -0.3); float phase = mix(phaseBackward, phaseForward, 0.5); s.direct_scattering = transmittance * phase * s.density; // sun colored - only direct scatter // multiScatter is indirect/diffuse - should be ambient colored, not sun colored float multiScatter = (1.0 - transmittance) * 0.25 * (0.3 + 0.7 * s.density); float skyTransmittance = exp(-texture(shade_tex, getCurvedP(eye, p + windOffset)).g); float dimensionalProfile = cloud.r; // ambient gets both sky terms AND the indirect multiple scatter float skyAmbient = pow(1.0 - dimensionalProfile, 1.2) * skyTransmittance; float multiScatterAmbient = 0.2 * dimensionalProfile; // p.z is offset so 0 = cloud base; add cloud_base_z_norm back to recover // height above cloud_field_center, i.e. height above the reflecting ground/sea surface. float heightAboveGround_m = (p.z + cloud_base_z_norm) * VOXEL_FIELD_HEIGHT_M; float groundBounceFalloff = 1.0 - smoothstep(0.0, GROUND_BOUNCE_MAX_HEIGHT_M, heightAboveGround_m); groundBounce = ground_albedo.rgb * (1.0 - skyTransmittance) * groundBounceFalloff * 0.1 * dayFactor; s.ambient_scattering = (skyAmbient + multiScatterAmbient) * dayFactor; multiScatterTerm = multiScatter * dayFactor; } if (s.density > 0.0) { if (lastAir >= lastCloud) { // Enter cloud, so apply any adjustment for the aerial pespective of the air // between this entry point and the last cloud float z = distance * VOXEL_FIELD_WIDTH_M / zscaleFactor; vec3 P = get_view_space_from_vs_depth(texcoord, z); vec4 ap = get_aerial_perspective(raw_texcoord, P); vec4 partialAP = diff_aerial_perspective(ap, lastAP); float occlusion = (1.0 - clamp(lreturn.light_absorption, 0.0, 1.0)); lreturn.light_absorption += (1.0 - partialAP.a) * occlusion; lreturn.intensity += partialAP.rgb * occlusion; lastAP = ap; } float occlusion = (1.0 - clamp(lreturn.light_absorption, 0.0, 1.0)); lreturn.light_absorption += s.density * occlusion * DENSITY_SCALE; // Combine the direct and ambient calculations to determine the overall intensity // Direct term - uses (almost) sun colour lreturn.intensity += directColor * s.direct_scattering * occlusion * DENSITY_SCALE; // Ambient term - uses the neutral sky colour multiplied by the sun magnitude lreturn.intensity += PIPELINE_RADIANCE_SCALE * sunMagnitude * skyAmbientColor * s.ambient_scattering * s.density * occlusion * DENSITY_SCALE; // Multiscatter term - completely neutral, scaled by sun magnitude lreturn.intensity += PIPELINE_RADIANCE_SCALE * multiScatterColor * multiScatterTerm * occlusion * DENSITY_SCALE; // Ground bounce lreturn.intensity += PIPELINE_RADIANCE_SCALE * sunRadiance * groundBounce * s.density * occlusion * DENSITY_SCALE; if (lreturn.first_hit < 0.0) lreturn.first_hit = distance; lastCloud = distance; } else { lastAir = distance; } if (lreturn.light_absorption > 0.98) { // Reached maximum density, so don't march further. An opaque cloud will have no silver lining, so no need to add it here. lreturn.light_absorption = 1.0; return lreturn; } distance += s.sdf; } // We've reached the end of the ray march, so add any silver lining effect. float sl = getSilverLining(lreturn.light_absorption, cosTheta); lreturn.intensity += sunRadiance * sl; lreturn.light_absorption += sl * 0.3; // Cancel any aerial perspective showing through from the background that // we've already applied between clouds. if (lastAP.a > 0.0) { vec4 partialAP = diff_aerial_perspective(vec4(0.0, 0.0, 0.0, 1.0), lastAP); float occlusion = (1.0 - clamp(lreturn.light_absorption, 0.0, 1.0)); lreturn.light_absorption += (1.0 - partialAP.a) * occlusion; lreturn.intensity += partialAP.rgb * occlusion; } return lreturn; } void main() { mat3 zup = mat3(fg_CameraZUpMatrix); vec3 wdir = w_pos - fg_CameraPositionCart; vec3 dir = normalize(zup * wdir) * VOXEL_SCALE; // Take into account that the voxel space is not a cube by increasing the Z-factor float zscaleFactor = length(dir); // Get a Z-up eyepoint relative to the center of the cloud field in the X-Y plane, and offset to place the bottom of the field at the cloudbase. vec3 cameraEye = (zup * (fg_CameraPositionCart - cloud_field_center)) / EYE_SCALE + vec3(0.5, 0.5, - cloud_base_z_norm); vec4 color = vec4(0.0, 0.0, 0.0, 0.0); // Convert the logarithmic depth value into metres, and then scale to the voxel resolution // so we know how far to search before we reach something solid. This needs to take into account that the voxel space is not a cube by adjusting for the // actual length of the "normalized" direction. float max_depth_m = logdepth_decode(texture(depth_tex, texcoord).r); float max_depth_vx = min(max_depth_m / VOXEL_FIELD_WIDTH_M, MAX_DIST) * zscaleFactor; ray_data ray = cloudRayMarch(cameraEye, dir, MIN_DIST, max_depth_vx, zscaleFactor); if (ray.light_absorption > 0.001) { color.rgb = ray.intensity; color.a = ray.light_absorption; float z = logdepth_prepare_vs_depth(ray.first_hit * VOXEL_FIELD_WIDTH_M / zscaleFactor); gl_FragDepth = logdepth_encode(z); } else { gl_FragDepth = 1.0; } // Only pre-expose when not rendering to the environment map. // We want the non-exposed radiance values for IBL. color.rgb = apply_exposure(color.rgb); fragColor = color; }