Add sky light, physical atmosphere, auto exposure and bloom to the GL 4.6 renderer

shaders/atmosphere_ambient.comp

@@ -0,0 +1,127 @@
+// what the cpu needs from the atmosphere, worked out from the sky-view table in one work group: the light arriving at the
+// observer from every direction as spherical harmonics (the ground below the horizon lit by the sun and the sky), the
+// transmittance towards the sun, and the average light of the sky above the horizon and over the whole sky.
+// run after the sky-view table; the result is read back by the cpu a few frames later
+
+#include <common>
+#include <atmosphere_common.glsl>
+
+layout(binding = ATMOSPHERE_TRANSMITTANCE_UNIT) uniform sampler2D atmo_transmittance_lut;
+layout(binding = ATMOSPHERE_SKYVIEW_UNIT) uniform sampler2D atmo_skyview_lut;
+
+layout(local_size_x = ATMOSPHERE_AMBIENT_GROUP_SIZE) in;
+
+layout(std430, binding = ATMOSPHERE_RESULT_SSBO_BINDING) writeonly buffer atmosphere_result_ssbo
+{
+	vec4 result_radiance_sh[9];
+	vec4 result_sun_transmittance;
+	vec4 result_horizon;
+	vec4 result_average;
+	vec4 result_zenith_transmittance;
+};
+
+// directions of the upper hemisphere each invocation takes
+const uint DIRECTIONS_PER_INVOCATION = 4U;
+const uint DIRECTIONS = ATMOSPHERE_AMBIENT_GROUP_SIZE * DIRECTIONS_PER_INVOCATION;
+// the sky up to this sine of the elevation is the horizon (about 6 degrees)
+const float HORIZON_BAND = 0.1;
+
+// partial sums: the 9 coefficients, the light on an upward facing surface, the horizon and its weight, the sky average
+const int SUMS = 12;
+shared vec4 shared_sums[ATMOSPHERE_AMBIENT_GROUP_SIZE];
+
+void sh_basis(vec3 n, out float basis[9])
+{
+	basis[0] = 0.282095;
+	basis[1] = 0.488603 * n.y;
+	basis[2] = 0.488603 * n.z;
+	basis[3] = 0.488603 * n.x;
+	basis[4] = 1.092548 * n.x * n.y;
+	basis[5] = 1.092548 * n.y * n.z;
+	basis[6] = 0.315392 * (3.0 * n.z * n.z - 1.0);
+	basis[7] = 1.092548 * n.x * n.z;
+	basis[8] = 0.546274 * (n.x * n.x - n.y * n.y);
+}
+
+// sum of the value of every invocation, valid in invocation 0
+vec4 reduce(vec4 value)
+{
+	uint index = gl_LocalInvocationID.x;
+	barrier();
+	shared_sums[index] = value;
+	barrier();
+	for (uint stride = ATMOSPHERE_AMBIENT_GROUP_SIZE / 2U; stride > 0U; stride >>= 1)
+	{
+		if (index < stride)
+			shared_sums[index] += shared_sums[index + stride];
+		barrier();
+	}
+	return shared_sums[0];
+}
+
+void main()
+{
+	uint index = gl_LocalInvocationID.x;
+	vec3 sundirection = normalize(atmo_sun_direction.xyz);
+
+	vec3 sums[SUMS];
+	for (int i = 0; i < SUMS; ++i)
+		sums[i] = vec3(0.0);
+
+	// even spread over the upper hemisphere: a fibonacci spiral in the sine of the elevation
+	const float goldenangle = 2.39996323;
+	const float solidangle = 2.0 * ATMO_PI / float(DIRECTIONS);
+	for (uint k = 0U; k < DIRECTIONS_PER_INVOCATION; ++k)
+	{
+		uint i = index * DIRECTIONS_PER_INVOCATION + k;
+		float y = (float(i) + 0.5) / float(DIRECTIONS);
+		float r = sqrt(max(1.0 - y * y, 0.0));
+		float azimuth = goldenangle * float(i);
+		vec3 direction = vec3(r * cos(azimuth), y, r * sin(azimuth));
+
+		vec3 radiance = textureLod(atmo_skyview_lut, atmo_skyview_uv_world(direction, sundirection), 0.0).rgb;
+		float basis[9];
+		sh_basis(direction, basis);
+		for (int c = 0; c < 9; ++c)
+			sums[c] += radiance * (basis[c] * solidangle);
+		sums[9] += radiance * (y * solidangle);
+		if (y < HORIZON_BAND)
+		{
+			sums[10] += radiance;
+			sums[11] += vec3(1.0);
+		}
+	}
+
+	vec3 total[SUMS];
+	for (int i = 0; i < SUMS; ++i)
+		total[i] = reduce(vec4(sums[i], 0.0)).rgb;
+
+	if (index != 0U)
+		return;
+
+	// the transmittance towards the sun, zero once the planet is in the way
+	float radius = atmo_observer_radius;
+	vec3 origin = vec3(0.0, radius, 0.0);
+	vec3 suntransmittance = texture(atmo_transmittance_lut, atmo_transmittance_uv(radius, sundirection.y)).rgb;
+	if (atmo_ray_sphere_nearest(origin, sundirection, atmo_bottom_radius) >= 0.0)
+		suntransmittance = vec3(0.0);
+
+	// the ground, a lambertian surface lit by the sun and the sky, fills the lower hemisphere with uniform light. its
+	// spherical harmonics have only the constant term and the vertical one: for a uniform radiance L over the lower
+	// hemisphere, 0.282095 * 2 pi * L and 0.488603 * (-pi) * L
+	// the band 0 coefficient of the sky alone over the solid angle of the hemisphere
+	vec3 skyaverage = total[0] / (0.282095 * 2.0 * ATMO_PI);
+	vec3 irradiance = suntransmittance * max(sundirection.y, 0.0) + total[9];
+	vec3 groundradiance = atmo_ground.rgb / ATMO_PI * irradiance;
+
+	total[0] += groundradiance * (0.282095 * 2.0 * ATMO_PI);
+	total[1] += groundradiance * (0.488603 * -ATMO_PI);
+	for (int c = 0; c < 9; ++c)
+		result_radiance_sh[c] = vec4(total[c], 0.0);
+
+	result_sun_transmittance = vec4(suntransmittance, 0.0);
+	result_horizon = vec4(total[10] / max(total[11].x, 1.0), 0.0);
+	result_average = vec4(skyaverage, 0.0);
+	// what the sun high in the sky looks like at the observer, for the white balance
+	result_zenith_transmittance = vec4(texture(atmo_transmittance_lut, atmo_transmittance_uv(radius, 1.0)).rgb, 0.0);
+}

shaders/atmosphere_common.glsl

@@ -0,0 +1,194 @@
+// physically based sky after S. Hillaire, "A Scalable and Production Ready Sky and Atmosphere Rendering Technique"
+// (EGSR 2020): transmittance, multiple scattering and sky-view look-up tables. distances in kilometres, coefficients per
+// kilometre. the planet is centred at the origin; the look-up tables are worked out in a frame with z up at the observer,
+// the world (y up) is mapped to it by the callers. radiance is that of a sun of illuminance 1
+
+layout(std140, binding = ATMOSPHERE_UBO_BINDING) uniform atmosphere_ubo
+{
+	// rgb: rayleigh scattering, w: rayleigh scale height
+	vec4 atmo_rayleigh;
+	// x: mie scattering, y: mie extinction, z: mie scale height, w: mie phase asymmetry
+	vec4 atmo_mie;
+	// rgb: ozone absorption, w: altitude of the middle of the ozone layer
+	vec4 atmo_ozone;
+	// rgb: albedo of the ground, w: half width of the ozone layer
+	vec4 atmo_ground;
+	// world space direction towards the sun (y up)
+	vec4 atmo_sun_direction;
+	float atmo_bottom_radius;
+	float atmo_top_radius;
+	// distance of the observer from the centre of the planet
+	float atmo_observer_radius;
+};
+
+const float ATMO_PI = 3.14159265359;
+
+struct atmo_medium_s
+{
+	vec3 scattering;
+	vec3 extinction;
+	vec3 rayleigh_scattering;
+	vec3 mie_scattering;
+};
+
+atmo_medium_s atmo_sample_medium(float radius)
+{
+	float height = max(radius - atmo_bottom_radius, 0.0);
+	float rayleighdensity = exp(-height / atmo_rayleigh.w);
+	float miedensity = exp(-height / atmo_mie.z);
+	float ozonedensity = max(0.0, 1.0 - abs(height - atmo_ozone.w) / atmo_ground.w);
+
+	atmo_medium_s medium;
+	medium.rayleigh_scattering = atmo_rayleigh.rgb * rayleighdensity;
+	medium.mie_scattering = vec3(atmo_mie.x * miedensity);
+	medium.scattering = medium.rayleigh_scattering + medium.mie_scattering;
+	medium.extinction = medium.rayleigh_scattering + vec3(atmo_mie.y * miedensity) + atmo_ozone.rgb * ozonedensity;
+	return medium;
+}
+
+// distance along a ray (unit direction) to the nearest point ahead of the origin where it meets the sphere, -1 if it does not
+float atmo_ray_sphere_nearest(vec3 origin, vec3 direction, float radius)
+{
+	float b = dot(origin, direction);
+	float c = dot(origin, origin) - radius * radius;
+	float discriminant = b * b - c;
+	if (discriminant < 0.0)
+		return -1.0;
+	float root = sqrt(discriminant);
+	float t0 = -b - root;
+	float t1 = -b + root;
+	if (t0 < 0.0 && t1 < 0.0)
+		return -1.0;
+	if (t0 < 0.0)
+		return t1;
+	if (t1 < 0.0)
+		return t0;
+	return min(t0, t1);
+}
+
+float atmo_rayleigh_phase(float cosangle)
+{
+	return 3.0 / (16.0 * ATMO_PI) * (1.0 + cosangle * cosangle);
+}
+
+// Cornette-Shanks
+float atmo_mie_phase(float cosangle)
+{
+	float g = atmo_mie.w;
+	float k = 3.0 / (8.0 * ATMO_PI) * (1.0 - g * g) / (2.0 + g * g);
+	return k * (1.0 + cosangle * cosangle) / pow(max(1.0 + g * g - 2.0 * g * cosangle, 1e-4), 1.5);
+}
+
+// texel centres of a table of the given resolution span [0, 1] of the parameter
+float atmo_unit_to_sub_uv(float u, float resolution)
+{
+	return (u + 0.5 / resolution) * (resolution / (resolution + 1.0));
+}
+
+float atmo_sub_uv_to_unit(float u, float resolution)
+{
+	return (u - 0.5 / resolution) * (resolution / (resolution - 1.0));
+}
+
+// transmittance table: x the distance to the top of the atmosphere along the view, y the height
+vec2 atmo_transmittance_uv(float radius, float cosviewzenith)
+{
+	float H = sqrt(max(0.0, atmo_top_radius * atmo_top_radius - atmo_bottom_radius * atmo_bottom_radius));
+	float rho = sqrt(max(0.0, radius * radius - atmo_bottom_radius * atmo_bottom_radius));
+	float discriminant = radius * radius * (cosviewzenith * cosviewzenith - 1.0) + atmo_top_radius * atmo_top_radius;
+	float d = max(0.0, -radius * cosviewzenith + sqrt(max(discriminant, 0.0)));
+	float dmin = atmo_top_radius - radius;
+	float dmax = rho + H;
+	return vec2((d - dmin) / max(dmax - dmin, 1e-6), rho / H);
+}
+
+void atmo_transmittance_params(vec2 uv, out float radius, out float cosviewzenith)
+{
+	float H = sqrt(max(0.0, atmo_top_radius * atmo_top_radius - atmo_bottom_radius * atmo_bottom_radius));
+	float rho = H * uv.y;
+	radius = sqrt(rho * rho + atmo_bottom_radius * atmo_bottom_radius);
+	float dmin = atmo_top_radius - radius;
+	float dmax = rho + H;
+	float d = dmin + uv.x * (dmax - dmin);
+	cosviewzenith = (d == 0.0) ? 1.0 : (H * H - rho * rho - d * d) / (2.0 * radius * d);
+	cosviewzenith = clamp(cosviewzenith, -1.0, 1.0);
+}
+
+// multiple scattering table: x the cosine of the sun zenith angle, y the height
+vec2 atmo_multiscattering_uv(float radius, float cossunzenith)
+{
+	vec2 uv = vec2(cossunzenith * 0.5 + 0.5, (radius - atmo_bottom_radius) / (atmo_top_radius - atmo_bottom_radius));
+	return clamp(vec2(atmo_unit_to_sub_uv(uv.x, float(ATMOSPHERE_MULTISCATTERING_SIZE)), atmo_unit_to_sub_uv(uv.y, float(ATMOSPHERE_MULTISCATTERING_SIZE))), 0.0, 1.0);
+}
+
+// sky-view table: x the angle between the view and the sun around the zenith (sqrt spaced), y the angle of the view from
+// the zenith, packed towards the horizon, which splits the table in the sky (lower half) and the ground (upper half)
+vec2 atmo_skyview_uv(bool intersectsground, float cosviewzenith, float coslightview, float radius)
+{
+	float vhorizon = sqrt(max(radius * radius - atmo_bottom_radius * atmo_bottom_radius, 0.0));
+	float cosbeta = vhorizon / radius;
+	float beta = acos(clamp(cosbeta, -1.0, 1.0));
+	float zenithhorizonangle = ATMO_PI - beta;
+
+	vec2 uv;
+	if (!intersectsground)
+	{
+		float coord = acos(clamp(cosviewzenith, -1.0, 1.0)) / zenithhorizonangle;
+		coord = 1.0 - coord;
+		coord = sqrt(max(coord, 0.0));
+		coord = 1.0 - coord;
+		uv.y = coord * 0.5;
+	}
+	else
+	{
+		float coord = (acos(clamp(cosviewzenith, -1.0, 1.0)) - zenithhorizonangle) / max(beta, 1e-6);
+		coord = sqrt(max(coord, 0.0));
+		uv.y = coord * 0.5 + 0.5;
+	}
+	uv.x = sqrt(clamp(-coslightview * 0.5 + 0.5, 0.0, 1.0));
+	return vec2(atmo_unit_to_sub_uv(uv.x, float(ATMOSPHERE_SKYVIEW_WIDTH)), atmo_unit_to_sub_uv(uv.y, float(ATMOSPHERE_SKYVIEW_HEIGHT)));
+}
+
+void atmo_skyview_params(vec2 uv, float radius, out float cosviewzenith, out float coslightview)
+{
+	uv = vec2(atmo_sub_uv_to_unit(uv.x, float(ATMOSPHERE_SKYVIEW_WIDTH)), atmo_sub_uv_to_unit(uv.y, float(ATMOSPHERE_SKYVIEW_HEIGHT)));
+	uv = clamp(uv, 0.0, 1.0);
+
+	float vhorizon = sqrt(max(radius * radius - atmo_bottom_radius * atmo_bottom_radius, 0.0));
+	float cosbeta = vhorizon / radius;
+	float beta = acos(clamp(cosbeta, -1.0, 1.0));
+	float zenithhorizonangle = ATMO_PI - beta;
+
+	if (uv.y < 0.5)
+	{
+		float coord = 2.0 * uv.y;
+		coord = 1.0 - coord;
+		coord *= coord;
+		coord = 1.0 - coord;
+		cosviewzenith = cos(zenithhorizonangle * coord);
+	}
+	else
+	{
+		float coord = uv.y * 2.0 - 1.0;
+		coord *= coord;
+		cosviewzenith = cos(zenithhorizonangle + beta * coord);
+	}
+	float coord = uv.x * uv.x;
+	coslightview = -(coord * 2.0 - 1.0);
+}
+
+// the same mapping for a world space view direction and sun direction (y up), as seen by the observer
+vec2 atmo_skyview_uv_world(vec3 viewdirection, vec3 sundirection)
+{
+	float radius = atmo_observer_radius;
+	vec3 origin = vec3(0.0, 0.0, radius);
+	float cosviewzenith = viewdirection.y;
+	// angle between the view and the sun around the zenith
+	vec2 viewflat = viewdirection.xz;
+	vec2 sunflat = sundirection.xz;
+	float lengths = length(viewflat) * length(sunflat);
+	float coslightview = (lengths > 1e-6) ? clamp(dot(viewflat, sunflat) / lengths, -1.0, 1.0) : 1.0;
+	vec3 localview = vec3(sqrt(max(1.0 - cosviewzenith * cosviewzenith, 0.0)), 0.0, cosviewzenith);
+	bool intersectsground = atmo_ray_sphere_nearest(origin, localview, atmo_bottom_radius) >= 0.0;
+	return atmo_skyview_uv(intersectsground, cosviewzenith, coslightview, radius);
+}

shaders/atmosphere_integrate.glsl

@@ -0,0 +1,105 @@
+// scattering along a ray through the atmosphere, see atmosphere_common.glsl. the includer declares the samplers
+// atmo_transmittance_lut and (when multiple scattering is used) atmo_multiscattering_lut
+
+// above the ground by this much (km), so the test for the shadow of the planet does not hit the point itself
+const float ATMO_PLANET_RADIUS_OFFSET = 0.01;
+
+struct atmo_scattering_s
+{
+	// light scattered towards the origin
+	vec3 luminance;
+	// transmittance from the origin to the end of the ray
+	vec3 transmittance;
+	// light scattered once more, as if the light arriving at every point was 1 (multiple scattering table)
+	vec3 multiscattering_as_1;
+};
+
+vec3 atmo_transmittance_to_top(float radius, float cossunzenith)
+{
+	return texture(atmo_transmittance_lut, atmo_transmittance_uv(radius, cossunzenith)).rgb;
+}
+
+// Origin and Direction (unit) in the frame of the planet. with Ground the light the ground reflects at the end of the ray
+// is included. Physicalphase: rayleigh and mie phase functions, otherwise isotropic (multiple scattering table).
+// Multiscattering: light scattered more than once comes from the multiple scattering table
+atmo_scattering_s atmo_integrate(vec3 origin, vec3 direction, vec3 sundirection, bool ground, int samples, bool physicalphase, bool multiscattering)
+{
+	atmo_scattering_s result;
+	result.luminance = vec3(0.0);
+	result.transmittance = vec3(1.0);
+	result.multiscattering_as_1 = vec3(0.0);
+
+	float tbottom = atmo_ray_sphere_nearest(origin, direction, atmo_bottom_radius);
+	float ttop = atmo_ray_sphere_nearest(origin, direction, atmo_top_radius);
+	float tmax;
+	if (tbottom < 0.0)
+	{
+		if (ttop < 0.0)
+			return result;
+		tmax = ttop;
+	}
+	else
+	{
+		tmax = (ttop > 0.0) ? min(ttop, tbottom) : tbottom;
+	}
+
+	float cosangle = dot(direction, sundirection);
+	const float uniformphase = 1.0 / (4.0 * ATMO_PI);
+	float miephase = physicalphase ? atmo_mie_phase(cosangle) : uniformphase;
+	float rayleighphase = physicalphase ? atmo_rayleigh_phase(cosangle) : uniformphase;
+
+	vec3 throughput = vec3(1.0);
+	float invsamples = 1.0 / float(samples);
+	for (int s = 0; s < samples; ++s)
+	{
+		// the steps grow with the distance: the medium near the observer matters most
+		float t0 = float(s) * invsamples;
+		float t1 = float(s + 1) * invsamples;
+		t0 = t0 * t0 * tmax;
+		t1 = t1 * t1 * tmax;
+		float dt = t1 - t0;
+		float t = t0 + dt * 0.3;
+
+		vec3 position = origin + t * direction;
+		float radius = length(position);
+		vec3 up = position / radius;
+		atmo_medium_s medium = atmo_sample_medium(radius);
+		vec3 extinction = max(medium.extinction, vec3(1e-7));
+		vec3 sampletransmittance = exp(-extinction * dt);
+
+		float cossunzenith = dot(sundirection, up);
+		vec3 suntransmittance = atmo_transmittance_to_top(radius, cossunzenith);
+		float earthshadow = (atmo_ray_sphere_nearest(position + up * ATMO_PLANET_RADIUS_OFFSET, sundirection, atmo_bottom_radius) >= 0.0) ? 0.0 : 1.0;
+
+		vec3 multiscatteredluminance = vec3(0.0);
+#ifdef ATMO_MULTISCATTERING_LUT
+		if (multiscattering)
+			multiscatteredluminance = texture(atmo_multiscattering_lut, atmo_multiscattering_uv(radius, cossunzenith)).rgb;
+#endif
+
+		vec3 phasescattering = medium.mie_scattering * miephase + medium.rayleigh_scattering * rayleighphase;
+		vec3 S = earthshadow * suntransmittance * phasescattering + multiscatteredluminance * medium.scattering;
+		// analytic integral over the step of the light scattered at the sample, attenuated on the way to its start
+		vec3 Sint = (S - S * sampletransmittance) / extinction;
+		result.luminance += throughput * Sint;
+
+		vec3 MSint = (medium.scattering - medium.scattering * sampletransmittance) / extinction;
+		result.multiscattering_as_1 += throughput * MSint;
+
+		throughput *= sampletransmittance;
+	}
+
+	if (ground && tbottom > 0.0 && tmax == tbottom)
+	{
+		vec3 position = origin + tbottom * direction;
+		float radius = length(position);
+		vec3 up = position / radius;
+		float cossunzenith = dot(sundirection, up);
+		vec3 suntransmittance = atmo_transmittance_to_top(radius, cossunzenith);
+		float NdotL = clamp(dot(up, sundirection), 0.0, 1.0);
+		result.luminance += suntransmittance * throughput * NdotL * atmo_ground.rgb / ATMO_PI;
+	}
+
+	result.transmittance = throughput;
+	return result;
+}

shaders/atmosphere_multiscattering.comp

@@ -0,0 +1,69 @@
+// light scattered more than once, for the multiple scattering table (Hillaire 2020, section 5.5): for a point at a height
+// lit by the sun at an angle, the second order scattering from 64 directions around it, summed as a geometric series of
+// all the orders. one work group per texel, one invocation per direction. run after the transmittance table
+
+#include <common>
+#include <atmosphere_common.glsl>
+
+layout(binding = ATMOSPHERE_TRANSMITTANCE_UNIT) uniform sampler2D atmo_transmittance_lut;
+#include <atmosphere_integrate.glsl>
+
+layout(local_size_x = 64) in;
+
+layout(binding = ATMOSPHERE_MULTISCATTERING_UNIT, rgba16f) writeonly uniform image2D multiscattering_target;
+
+const int MULTISCATTERING_SAMPLES = 20;
+const int SQRT_DIRECTIONS = 8;
+
+shared vec3 shared_luminance[64];
+shared vec3 shared_multiscattering[64];
+
+void main()
+{
+	ivec2 texel = ivec2(gl_WorkGroupID.xy);
+	uint index = gl_LocalInvocationID.x;
+
+	vec2 uv = (vec2(texel) + 0.5) / float(ATMOSPHERE_MULTISCATTERING_SIZE);
+	uv = vec2(atmo_sub_uv_to_unit(uv.x, float(ATMOSPHERE_MULTISCATTERING_SIZE)), atmo_sub_uv_to_unit(uv.y, float(ATMOSPHERE_MULTISCATTERING_SIZE)));
+	uv = clamp(uv, 0.0, 1.0);
+
+	float cossunzenith = uv.x * 2.0 - 1.0;
+	vec3 sundirection = vec3(0.0, sqrt(max(1.0 - cossunzenith * cossunzenith, 0.0)), cossunzenith);
+	float radius = atmo_bottom_radius + clamp(uv.y + 1e-4, 0.0, 1.0) * (atmo_top_radius - atmo_bottom_radius - 1e-3);
+	vec3 origin = vec3(0.0, 0.0, radius);
+
+	// directions spread evenly over the sphere
+	float a = (0.5 + float(index / uint(SQRT_DIRECTIONS))) / float(SQRT_DIRECTIONS);
+	float b = (0.5 + float(index % uint(SQRT_DIRECTIONS))) / float(SQRT_DIRECTIONS);
+	float theta = 2.0 * ATMO_PI * a;
+	float phi = acos(clamp(1.0 - 2.0 * b, -1.0, 1.0));
+	vec3 direction = vec3(cos(theta) * sin(phi), sin(theta) * sin(phi), cos(phi));
+
+	atmo_scattering_s result = atmo_integrate(origin, direction, sundirection, true, MULTISCATTERING_SAMPLES, false, false);
+
+	const float spheresolidangle = 4.0 * ATMO_PI;
+	const float directions = float(SQRT_DIRECTIONS * SQRT_DIRECTIONS);
+	shared_luminance[index] = result.luminance * spheresolidangle / directions;
+	shared_multiscattering[index] = result.multiscattering_as_1 * spheresolidangle / directions;
+	barrier();
+
+	for (uint stride = 32U; stride > 0U; stride >>= 1)
+	{
+		if (index < stride)
+		{
+			shared_luminance[index] += shared_luminance[index + stride];
+			shared_multiscattering[index] += shared_multiscattering[index + stride];
+		}
+		barrier();
+	}
+
+	if (index != 0U)
+		return;
+
+	const float isotropicphase = 1.0 / spheresolidangle;
+	vec3 inscatteredluminance = shared_luminance[0] * isotropicphase;
+	vec3 multiscatteringas1 = shared_multiscattering[0] * isotropicphase;
+	// every further order of scattering is the previous one times multiscatteringas1
+	vec3 luminance = inscatteredluminance / max(vec3(1.0) - multiscatteringas1, vec3(1e-3));
+	imageStore(multiscattering_target, texel, vec4(luminance, 1.0));
+}

shaders/atmosphere_skyview.comp

@@ -0,0 +1,39 @@
+// light of the sky seen by the observer in every direction, for the sky-view table (atmosphere_common.glsl), which the
+// sky, the ambient light and the fog colour are read from. one invocation per texel. run every frame the sky is drawn with
+// the atmosphere, after the transmittance and multiple scattering tables
+
+#include <common>
+#include <atmosphere_common.glsl>
+
+layout(binding = ATMOSPHERE_TRANSMITTANCE_UNIT) uniform sampler2D atmo_transmittance_lut;
+layout(binding = ATMOSPHERE_MULTISCATTERING_UNIT) uniform sampler2D atmo_multiscattering_lut;
+#define ATMO_MULTISCATTERING_LUT
+#include <atmosphere_integrate.glsl>
+
+layout(local_size_x = 8, local_size_y = 8) in;
+
+layout(binding = ATMOSPHERE_SKYVIEW_UNIT, rgba16f) writeonly uniform image2D skyview_target;
+
+const int SKYVIEW_SAMPLES = 32;
+
+void main()
+{
+	ivec2 texel = ivec2(gl_GlobalInvocationID.xy);
+	if (texel.x >= int(ATMOSPHERE_SKYVIEW_WIDTH) || texel.y >= int(ATMOSPHERE_SKYVIEW_HEIGHT))
+		return;
+
+	vec2 uv = (vec2(texel) + 0.5) / vec2(ATMOSPHERE_SKYVIEW_WIDTH, ATMOSPHERE_SKYVIEW_HEIGHT);
+	float radius = atmo_observer_radius;
+	float cosviewzenith;
+	float coslightview;
+	atmo_skyview_params(uv, radius, cosviewzenith, coslightview);
+
+	float cossunzenith = clamp(atmo_sun_direction.y, -1.0, 1.0);
+	vec3 sundirection = normalize(vec3(sqrt(max(1.0 - cossunzenith * cossunzenith, 0.0)), 0.0, cossunzenith));
+	float sinviewzenith = sqrt(max(1.0 - cosviewzenith * cosviewzenith, 0.0));
+	vec3 direction = vec3(sinviewzenith * coslightview, sinviewzenith * sqrt(max(1.0 - coslightview * coslightview, 0.0)), cosviewzenith);
+	vec3 origin = vec3(0.0, 0.0, radius);
+
+	atmo_scattering_s result = atmo_integrate(origin, direction, sundirection, false, SKYVIEW_SAMPLES, true, true);
+	imageStore(skyview_target, texel, vec4(result.luminance, 1.0));
+}

shaders/atmosphere_transmittance.comp

@@ -0,0 +1,40 @@
+// transmittance from a point in the atmosphere to its top, for the transmittance table (atmosphere_common.glsl).
+// one invocation per texel. run when the parameters of the atmosphere change
+
+#include <common>
+#include <atmosphere_common.glsl>
+
+layout(local_size_x = 8, local_size_y = 8) in;
+
+layout(binding = ATMOSPHERE_TRANSMITTANCE_UNIT, rgba16f) writeonly uniform image2D transmittance_target;
+
+const int TRANSMITTANCE_SAMPLES = 40;
+
+void main()
+{
+	ivec2 texel = ivec2(gl_GlobalInvocationID.xy);
+	if (texel.x >= int(ATMOSPHERE_TRANSMITTANCE_WIDTH) || texel.y >= int(ATMOSPHERE_TRANSMITTANCE_HEIGHT))
+		return;
+
+	vec2 uv = (vec2(texel) + 0.5) / vec2(ATMOSPHERE_TRANSMITTANCE_WIDTH, ATMOSPHERE_TRANSMITTANCE_HEIGHT);
+	float radius;
+	float cosviewzenith;
+	atmo_transmittance_params(uv, radius, cosviewzenith);
+
+	vec3 origin = vec3(0.0, 0.0, radius);
+	vec3 direction = vec3(sqrt(max(1.0 - cosviewzenith * cosviewzenith, 0.0)), 0.0, cosviewzenith);
+
+	float tbottom = atmo_ray_sphere_nearest(origin, direction, atmo_bottom_radius);
+	float ttop = atmo_ray_sphere_nearest(origin, direction, atmo_top_radius);
+	float tmax = (tbottom > 0.0) ? ((ttop > 0.0) ? min(ttop, tbottom) : tbottom) : max(ttop, 0.0);
+
+	vec3 opticaldepth = vec3(0.0);
+	float dt = tmax / float(TRANSMITTANCE_SAMPLES);
+	for (int s = 0; s < TRANSMITTANCE_SAMPLES; ++s)
+	{
+		vec3 position = origin + (float(s) + 0.5) * dt * direction;
+		opticaldepth += atmo_sample_medium(length(position)).extinction * dt;
+	}
+
+	imageStore(transmittance_target, texel, vec4(exp(-opticaldepth), 1.0));
+}

shaders/color.frag

@@ -10,22 +10,7 @@ layout(location = 0) out vec4 out_color;
 layout(location = 1) out vec4 out_motion;
 #endif
 
-// the skydome vertex colours carry Perez luminance scaled by a fixed exposure
-// (see CSkyDome::SetExposure), so the circumsolar region and the horizon band
-// routinely land well above 1.0. hard-clipping them, as this shader used to do,
-// flattened everything within ~20 degrees of the sun - and most of the lower sky -
-// into one solid white patch. compress the top end instead: below the knee the
-// response is unchanged, above it it rolls off towards SKY_SHOULDER_WHITE, so the
-// brightness gradient around the sun survives and the patch keeps a soft edge.
-const float SKY_SHOULDER_KNEE = 0.25;
-const float SKY_SHOULDER_WHITE = 0.95;
-
-vec3 sky_shoulder(vec3 color)
-{
-	const float range = SKY_SHOULDER_WHITE - SKY_SHOULDER_KNEE;
-	vec3 excess = max(color - vec3(SKY_SHOULDER_KNEE), vec3(0.0));
-	return min(color, vec3(SKY_SHOULDER_KNEE)) + range * (vec3(1.0) - exp(-excess / range));
-}
+#include <sky_shoulder.glsl>
 
 void main()
 {

shaders/envmap_prefilter.comp

@@ -0,0 +1,100 @@
+// prefilters the environment map with the GGX lobe for rough reflections (split sum approximation, Karis 2013), one level
+// of the prefiltered map per dispatch. filtered importance sampling (Krivanek, Colbert 2008): each sample reads the mipmap
+// of the source whose texels are about as large as the solid angle the sample stands for, so a few dozen samples give a
+// smooth result. run when an update of the environment map is complete
+
+#include <common>
+
+layout(local_size_x = ENVMAP_GGX_GROUP_SIZE, local_size_y = ENVMAP_GGX_GROUP_SIZE) in;
+
+layout(binding = 0, rgba16f) writeonly uniform imageCube prefilter_target;
+// the plain map with its box filtered mipmaps, on the unit gl::program binds envmap to
+uniform samplerCube envmap;
+
+// perceptual roughness of the level being written
+layout(location = 0) uniform float roughness;
+// size of the level being written, and of the first level of the source
+layout(location = 1) uniform int target_size;
+layout(location = 2) uniform float source_size;
+layout(location = 3) uniform int sample_count;
+
+const float PI = 3.14159265359;
+
+// direction through the centre of a texel of a cube map face, in the order and orientation of GL_TEXTURE_CUBE_MAP_POSITIVE_X..
+vec3 cube_direction(ivec3 texel, int size)
+{
+	vec2 uv = (vec2(texel.xy) + 0.5) / float(size) * 2.0 - 1.0;
+	switch (texel.z)
+	{
+	case 0: return vec3(1.0, -uv.y, -uv.x);
+	case 1: return vec3(-1.0, -uv.y, uv.x);
+	case 2: return vec3(uv.x, 1.0, uv.y);
+	case 3: return vec3(uv.x, -1.0, -uv.y);
+	case 4: return vec3(uv.x, -uv.y, 1.0);
+	default: return vec3(-uv.x, -uv.y, -1.0);
+	}
+}
+
+vec2 hammersley(uint i, uint n)
+{
+	uint bits = bitfieldReverse(i);
+	return vec2(float(i) / float(n), float(bits) * 2.3283064365386963e-10);
+}
+
+// half vector around N for the GGX distribution of linear roughness a
+vec3 importance_sample_ggx(vec2 xi, float a, vec3 N)
+{
+	float phi = 2.0 * PI * xi.x;
+	float costheta = sqrt((1.0 - xi.y) / (1.0 + (a * a - 1.0) * xi.y));
+	float sintheta = sqrt(max(1.0 - costheta * costheta, 0.0));
+	vec3 H = vec3(sintheta * cos(phi), sintheta * sin(phi), costheta);
+	vec3 up = abs(N.z) < 0.999 ? vec3(0.0, 0.0, 1.0) : vec3(1.0, 0.0, 0.0);
+	vec3 tangentx = normalize(cross(up, N));
+	vec3 tangenty = cross(N, tangentx);
+	return tangentx * H.x + tangenty * H.y + N * H.z;
+}
+
+float distribution_ggx(float NdotH, float a)
+{
+	float a2 = a * a;
+	float d = NdotH * NdotH * (a2 - 1.0) + 1.0;
+	return a2 / (PI * d * d);
+}
+
+void main()
+{
+	ivec3 texel = ivec3(gl_GlobalInvocationID);
+	if (texel.x >= target_size || texel.y >= target_size)
+		return;
+
+	vec3 N = normalize(cube_direction(texel, target_size));
+	// the view is taken along the normal, which is what the split sum approximation assumes
+	vec3 V = N;
+	float a = max(roughness * roughness, 1e-3);
+	float texelsolidangle = 4.0 * PI / (6.0 * source_size * source_size);
+
+	vec3 sum = vec3(0.0);
+	float weight = 0.0;
+	uint count = uint(sample_count);
+	for (uint i = 0U; i < count; ++i)
+	{
+		vec3 H = importance_sample_ggx(hammersley(i, count), a, N);
+		float VdotH = dot(V, H);
+		vec3 L = 2.0 * VdotH * H - V;
+		float NdotL = dot(N, L);
+		if (NdotL <= 0.0)
+			continue;
+		float NdotH = max(dot(N, H), 0.0);
+		// pdf of L: D * NdotH / (4 * VdotH)
+		float pdf = distribution_ggx(NdotH, a) * NdotH / max(4.0 * VdotH, 1e-4);
+		float samplesolidangle = 1.0 / (float(count) * max(pdf, 1e-6));
+		float level = max(0.5 * log2(samplesolidangle / texelsolidangle) + 1.0, 0.0);
+		vec3 color = textureLod(envmap, L, level).rgb;
+		if (any(isnan(color)) || any(isinf(color)))
+			color = vec3(0.0);
+		sum += max(color, vec3(0.0)) * NdotL;
+		weight += NdotL;
+	}
+
+	imageStore(prefilter_target, texel, vec4(sum / max(weight, 1e-4), 1.0));
+}

shaders/envmapping.glsl

@@ -1,5 +1,7 @@
 #if ENVMAP_ENABLED
 uniform samplerCube envmap;
+// the same map prefiltered with the GGX lobe (envmap_prefilter.comp), used when env_ggx is 1
+uniform samplerCube envmap_ggx;
 #endif
 
 vec3 envmap_color( vec3 normal )
@@ -39,3 +41,25 @@ vec3 envmap_color_lod(vec3 fragnormal, float lod)
     return vec3(0.5); // was vec3(0.0), match the non-LOD fallback
 #endif
 }
+
+// reflection of the surroundings in a surface of the given perceptual roughness (0 mirror, 1 matte)
+vec3 envmap_color_rough(vec3 fragnormal, float roughness)
+{
+	// the plain map, its mipmaps (box filtered) standing for the roughness
+	const float MAX_REFLECTION_LOD = 8.0;
+	vec3 boxcolor = envmap_color_lod(fragnormal, roughness * MAX_REFLECTION_LOD);
+#if ENVMAP_ENABLED
+	if (env_ggx == 0U || roughness <= ENVMAP_GGX_ROUGHNESS_MIN)
+		return boxcolor;
+	vec3 refvec = reflect(f_pos.xyz, fragnormal);
+	refvec = vec3(inv_view * vec4(refvec, 0.0));
+	float level = (roughness - ENVMAP_GGX_ROUGHNESS_MIN) / (1.0 - ENVMAP_GGX_ROUGHNESS_MIN) * float(ENVMAP_GGX_LEVELS - 1U);
+	vec3 ggxcolor = textureLod(envmap_ggx, refvec, level).rgb;
+	if (any(isnan(ggxcolor)) || any(isinf(ggxcolor))) ggxcolor = vec3(0.0);
+	ggxcolor = max(ggxcolor, vec3(0.0));
+	// the prefiltered map takes over over the first tenth of the roughness above its lowest level
+	return mix(boxcolor, ggxcolor, clamp((roughness - ENVMAP_GGX_ROUGHNESS_MIN) * 10.0, 0.0, 1.0));
+#else
+	return boxcolor;
+#endif
+}

shaders/exposure_average.comp

@@ -0,0 +1,79 @@
+// exposure from the luminance histogram (exposure_histogram.comp): the average log2 luminance of the pixels between the
+// low and high percentile (the darkest and brightest parts of the image do not move the exposure), adapted over time to,
+// faster towards the brighter. the exposure brings the adapted luminance to the key, within the limits. one work group,
+// one invocation per bin; the histogram is cleared for the next frame
+
+#include <common>
+
+layout(local_size_x = EXPOSURE_HISTOGRAM_BINS) in;
+
+layout(std430, binding = EXPOSURE_SSBO_BINDING) buffer exposure_ssbo
+{
+	uint histogram[EXPOSURE_HISTOGRAM_BINS];
+	float adapted_luminance;
+	float exposure;
+	float average_luminance;
+	uint valid;
+};
+
+// x: log2 of the lowest luminance binned, y: log2 luminance per bin
+layout(location = 0) uniform vec2 log_range;
+// x: time since the last frame (s), y: rate of adaptation to a darker image, z: to a brighter one (1/s), w: key
+layout(location = 1) uniform vec4 adaptation;
+// lowest and highest exposure (multipliers)
+layout(location = 2) uniform vec2 exposure_range;
+
+const float LOW_PERCENTILE = 0.4;
+const float HIGH_PERCENTILE = 0.95;
+
+shared float shared_counts[EXPOSURE_HISTOGRAM_BINS];
+
+void main()
+{
+	uint index = gl_LocalInvocationIndex;
+	shared_counts[index] = float(histogram[index]);
+	histogram[index] = 0U;
+	barrier();
+
+	if (index != 0U)
+		return;
+
+	float total = 0.0;
+	for (uint bin = 1U; bin < EXPOSURE_HISTOGRAM_BINS; ++bin)
+		total += shared_counts[bin];
+	// an all black image (or none) leaves the exposure as it was
+	if (total < 1.0)
+		return;
+
+	float low = total * LOW_PERCENTILE;
+	float high = total * HIGH_PERCENTILE;
+	float counted = 0.0;
+	float logsum = 0.0;
+	float weight = 0.0;
+	for (uint bin = 1U; bin < EXPOSURE_HISTOGRAM_BINS; ++bin)
+	{
+		float count = shared_counts[bin];
+		// the part of the bin between the percentiles
+		float inside = max(min(counted + count, high) - max(counted, low), 0.0);
+		counted += count;
+		float binluminance = log_range.x + (float(bin - 1U) + 0.5) * log_range.y;
+		logsum += inside * binluminance;
+		weight += inside;
+	}
+	float averagelog = (weight > 0.0) ? logsum / weight : log_range.x;
+	average_luminance = exp2(averagelog);
+
+	if (valid == 0U || !(adapted_luminance > 0.0))
+	{
+		adapted_luminance = average_luminance;
+		valid = 1U;
+	}
+	else
+	{
+		float adaptedlog = log2(adapted_luminance);
+		float rate = (averagelog < adaptedlog) ? adaptation.y : adaptation.z;
+		adaptedlog = mix(adaptedlog, averagelog, 1.0 - exp(-adaptation.x * rate));
+		adapted_luminance = exp2(adaptedlog);
+	}
+	exposure = clamp(adaptation.w / adapted_luminance, exposure_range.x, exposure_range.y);
+}

shaders/exposure_histogram.comp

@@ -0,0 +1,49 @@
+// luminance histogram of the image for the automatic exposure: every second pixel in both directions, binned by the log2
+// of its luminance. bin 0 takes the black pixels, which the average leaves out. exposure_average.comp reads and clears it
+
+#include <common>
+
+layout(local_size_x = EXPOSURE_HISTOGRAM_GROUP_SIZE, local_size_y = EXPOSURE_HISTOGRAM_GROUP_SIZE) in;
+
+layout(binding = 0) uniform sampler2D hdr_image;
+
+layout(std430, binding = EXPOSURE_SSBO_BINDING) buffer exposure_ssbo
+{
+	uint histogram[EXPOSURE_HISTOGRAM_BINS];
+	float adapted_luminance;
+	float exposure;
+	float average_luminance;
+	uint valid;
+};
+
+// x: log2 of the lowest luminance binned, y: number of bins per unit of log2 luminance
+layout(location = 0) uniform vec2 log_range;
+
+shared uint shared_histogram[EXPOSURE_HISTOGRAM_BINS];
+
+void main()
+{
+	uint index = gl_LocalInvocationIndex;
+	shared_histogram[index] = 0U;
+	barrier();
+
+	ivec2 pixel = ivec2(gl_GlobalInvocationID.xy) * 2;
+	ivec2 size = textureSize(hdr_image, 0);
+	if (pixel.x < size.x && pixel.y < size.y)
+	{
+		vec3 color = texelFetch(hdr_image, pixel, 0).rgb;
+		float luminance = dot(color, vec3(0.2126, 0.7152, 0.0722));
+		uint bin = 0U;
+		// also false for a NaN
+		if (luminance > 1e-6)
+		{
+			float position = clamp((log2(luminance) - log_range.x) * log_range.y, 0.0, float(EXPOSURE_HISTOGRAM_BINS - 2U));
+			bin = 1U + uint(position);
+		}
+		atomicAdd(shared_histogram[bin], 1U);
+	}
+	barrier();
+
+	if (shared_histogram[index] != 0U)
+		atomicAdd(histogram[index], shared_histogram[index]);
+}

shaders/light_common.glsl

@@ -308,12 +308,39 @@ vec3 material_diffuse()
     return (value > 0.0 ? diffusecolor * (param[1].x / value) : vec3(0.0));
 }
 
+// directional part of the sky ambient light (sky_ssbo, bands 1 and 2 of the spherical harmonics) for a world space direction
+vec3 sky_ambient_directional(vec3 n)
+{
+    return sky_ambient_sh[1].rgb * (0.488603 * n.y)
+         + sky_ambient_sh[2].rgb * (0.488603 * n.z)
+         + sky_ambient_sh[3].rgb * (0.488603 * n.x)
+         + sky_ambient_sh[4].rgb * (1.092548 * n.x * n.y)
+         + sky_ambient_sh[5].rgb * (1.092548 * n.y * n.z)
+         + sky_ambient_sh[6].rgb * (0.315392 * (3.0 * n.z * n.z - 1.0))
+         + sky_ambient_sh[7].rgb * (1.092548 * n.x * n.z)
+         + sky_ambient_sh[8].rgb * (0.546274 * (n.x * n.x - n.y * n.y));
+}
+
+// adds the directional part of the sky light to the ambient light the material starts from (the uniform ambient, or a share
+// of it, like the water takes). the same share of the directional part is added, so with the directional part at zero
+// (gfx.ambient.sh off) the result is what it always was
+vec3 add_sky_ambient(vec3 basecolor, vec3 fragnormal)
+{
+    if (sky_ambient_factor == 0.0)
+        return basecolor;
+    float share = dot(basecolor, vec3(1.0)) / max(dot(ambient, vec3(1.0)), 1e-4);
+    vec3 worldnormal = normalize(mat3(inv_view) * fragnormal);
+    return max(basecolor + sky_ambient_directional(worldnormal) * (sky_ambient_factor * share), vec3(0.0));
+}
+
 // [0] - diffuse, [1] - specular
 // do magic here
 vec3 apply_lights(vec3 fragcolor, vec3 fragnormal, vec3 texturecolor, float reflectivity, float specularity, float shadowtone)
 {
     // Diffuse colour is applied once, to the albedo; the light terms below are not scaled by it again
     texturecolor *= material_diffuse();
+    // the sky lights the faces turned up and towards the bright part of the sky more than the ones turned to the ground
+    fragcolor = add_sky_ambient(fragcolor, fragnormal);
     // Scale ambient before it gets tinted by the texture.
     // Sun, headlights and emission are added afterwards so they are NOT
     // attenuated by AMBIENT_SCALE - this only dims the indirect term.
@@ -326,9 +353,8 @@ vec3 apply_lights(vec3 fragcolor, vec3 fragnormal, vec3 texturecolor, float refl
     vec3 F0 = mix(vec3(0.04), texturecolor, metalic);
     vec3 fresnel = F0 + (1.0 - F0) * pow(1.0 - NdotV, 5.0);
 
-    const float MAX_REFLECTION_LOD = 8.0;
     float env_roughness = 1.0 - clamp(glossiness / max(abs(param[1].w), 1.0), 0.0, 1.0);
-    vec3 envcolor = envmap_color_lod(fragnormal, env_roughness * MAX_REFLECTION_LOD);
+    vec3 envcolor = envmap_color_rough(fragnormal, env_roughness);
 
     // Pre-integrated env BRDF: roughness/F0/view-dependent specular scale.
     // Replaces the old raw `fresnel` weighting so matte surfaces stop

shaders/mat_sunlessnormalmap_specgloss.frag

@@ -49,9 +49,8 @@ vec3 apply_lights_sunless(vec3 fragcolor, vec3 fragnormal, vec3 texturecolor, fl
     vec3 F0 = mix(vec3(0.04), texturecolor, metalic);
     vec3 fresnel = F0 + (1.0 - F0) * pow(1.0 - NdotV, 5.0);
 
-    const float MAX_REFLECTION_LOD = 8.0;
     float env_roughness = 1.0 - clamp(glossiness / max(abs(param[1].w), 1.0), 0.0, 1.0);
-    vec3 envcolor = envmap_color_lod(fragnormal, env_roughness * MAX_REFLECTION_LOD);
+    vec3 envcolor = envmap_color_rough(fragnormal, env_roughness);
 
     // Pre-integrated env BRDF (matches apply_lights): roughness-attenuated
     // reflection so matte cab surfaces no longer mirror the environment.

shaders/postfx_bloom_down.frag

@@ -0,0 +1,77 @@
+// one step down the bloom chain: a 13 tap filter (Jimenez, "Next generation post processing in Call of Duty: Advanced
+// Warfare", 2014) from a level twice as large. the first step, from the image itself, weighs the groups of taps by their
+// brightness (Karis average) so single bright pixels do not flicker, and applies the threshold.
+// param[0].x: 1 for the first step, param[0].y: threshold (0: none), param[0].z: soft knee
+
+in vec2 f_coords;
+layout(location = 0) out vec4 out_color;
+
+#texture (source_tex, 0, RGB)
+uniform sampler2D source_tex;
+
+#include <common>
+
+vec3 tap(vec2 uv)
+{
+	vec3 color = texture(source_tex, uv).rgb;
+	if (any(isnan(color)) || any(isinf(color)))
+		return vec3(0.0);
+	return clamp(color, vec3(0.0), vec3(65000.0));
+}
+
+float karis_weight(vec3 color)
+{
+	return 1.0 / (1.0 + dot(color, vec3(0.2126, 0.7152, 0.0722)));
+}
+
+void main()
+{
+	vec2 texel = 1.0 / vec2(textureSize(source_tex, 0));
+	vec2 uv = f_coords;
+
+	vec3 a = tap(uv + texel * vec2(-2.0, 2.0));
+	vec3 b = tap(uv + texel * vec2(0.0, 2.0));
+	vec3 c = tap(uv + texel * vec2(2.0, 2.0));
+	vec3 d = tap(uv + texel * vec2(-2.0, 0.0));
+	vec3 e = tap(uv);
+	vec3 f = tap(uv + texel * vec2(2.0, 0.0));
+	vec3 g = tap(uv + texel * vec2(-2.0, -2.0));
+	vec3 h = tap(uv + texel * vec2(0.0, -2.0));
+	vec3 i = tap(uv + texel * vec2(2.0, -2.0));
+	vec3 j = tap(uv + texel * vec2(-1.0, 1.0));
+	vec3 k = tap(uv + texel * vec2(1.0, 1.0));
+	vec3 l = tap(uv + texel * vec2(-1.0, -1.0));
+	vec3 m = tap(uv + texel * vec2(1.0, -1.0));
+
+	vec3 color;
+	if (param[0].x > 0.5)
+	{
+		vec3 g0 = (j + k + l + m) * 0.25;
+		vec3 g1 = (a + b + d + e) * 0.25;
+		vec3 g2 = (b + c + e + f) * 0.25;
+		vec3 g3 = (d + e + g + h) * 0.25;
+		vec3 g4 = (e + f + h + i) * 0.25;
+		float w0 = 0.5 * karis_weight(g0);
+		float w1 = 0.125 * karis_weight(g1);
+		float w2 = 0.125 * karis_weight(g2);
+		float w3 = 0.125 * karis_weight(g3);
+		float w4 = 0.125 * karis_weight(g4);
+		color = (g0 * w0 + g1 * w1 + g2 * w2 + g3 * w3 + g4 * w4) / (w0 + w1 + w2 + w3 + w4);
+
+		float threshold = param[0].y;
+		if (threshold > 0.0)
+		{
+			// quadratic soft knee below the threshold
+			float knee = max(threshold * param[0].z, 1e-4);
+			float brightness = max(color.r, max(color.g, color.b));
+			float soft = clamp(brightness - threshold + knee, 0.0, 2.0 * knee);
+			soft = soft * soft / (4.0 * knee);
+			color *= max(soft, brightness - threshold) / max(brightness, 1e-4);
+		}
+	}
+	else
+	{
+		color = e * 0.125 + (a + c + g + i) * 0.03125 + (b + d + f + h) * 0.0625 + (j + k + l + m) * 0.125;
+	}
+	out_color = vec4(color, 1.0);
+}

shaders/postfx_bloom_up.frag

@@ -0,0 +1,29 @@
+// one step up the bloom chain: the level below (half as large, with all the levels under it already added) blurred up with
+// a 3x3 tent filter and added to this level of the way down
+
+in vec2 f_coords;
+layout(location = 0) out vec4 out_color;
+
+#texture (lower_tex, 0, RGB)
+uniform sampler2D lower_tex;
+
+#texture (current_tex, 1, RGB)
+uniform sampler2D current_tex;
+
+void main()
+{
+	vec2 texel = 1.0 / vec2(textureSize(lower_tex, 0));
+	vec2 uv = f_coords;
+
+	vec3 sum = texture(lower_tex, uv).rgb * 4.0;
+	sum += (texture(lower_tex, uv + texel * vec2(0.0, 1.0)).rgb
+	      + texture(lower_tex, uv + texel * vec2(-1.0, 0.0)).rgb
+	      + texture(lower_tex, uv + texel * vec2(1.0, 0.0)).rgb
+	      + texture(lower_tex, uv + texel * vec2(0.0, -1.0)).rgb) * 2.0;
+	sum += texture(lower_tex, uv + texel * vec2(-1.0, 1.0)).rgb
+	     + texture(lower_tex, uv + texel * vec2(1.0, 1.0)).rgb
+	     + texture(lower_tex, uv + texel * vec2(-1.0, -1.0)).rgb
+	     + texture(lower_tex, uv + texel * vec2(1.0, -1.0)).rgb;
+
+	out_color = vec4(texture(current_tex, uv).rgb + sum / 16.0, 1.0);
+}

shaders/postfx_tonemapping.frag

@@ -5,11 +5,35 @@ layout(location = 0) out vec4 out_color;
 #texture (tex1, 0, RGB)
 uniform sampler2D tex1;
 
+// glow of the bright parts of the image (postfx_bloom_*.frag)
+#texture (bloom_tex, 1, RGB)
+uniform sampler2D bloom_tex;
+
+#include <common>
+
+// exposure worked out by exposure_average.comp
+layout(std430, binding = EXPOSURE_SSBO_BINDING) readonly buffer exposure_ssbo
+{
+	uint histogram[EXPOSURE_HISTOGRAM_BINS];
+	float adapted_luminance;
+	float exposure;
+	float average_luminance;
+	uint valid;
+};
+
 #include <tonemapping.glsl>
 
+// param[0].x: share of the bloom in the image (0: no bloom, the texture is not read), param[0].y: 1 if the automatic
+// exposure is applied, param[0].z: fixed exposure multiplier, param[0].w: scale of the bloom texture (it holds the sum of
+// the levels of the chain)
 void main()
 {
 	vec2 texcoord = f_coords;
 	vec3 hdr_color = texture(tex1, texcoord).xyz;
-	out_color = tonemap(vec4(hdr_color, 1.0));
-}  
+	if (param[0].x > 0.0)
+		hdr_color = mix(hdr_color, texture(bloom_tex, texcoord).rgb * param[0].w, param[0].x);
+	float exposurescale = param[0].z;
+	if (param[0].y > 0.5 && exposure > 0.0)
+		exposurescale *= exposure;
+	out_color = tonemap(vec4(hdr_color * exposurescale, 1.0));
+}

shaders/sky_atmosphere.frag

@@ -0,0 +1,54 @@
+// the sky dome drawn with the light of the physical atmosphere (atmosphere_common.glsl): the sky-view table looked up in
+// the direction of each fragment. the colour of the old dome (vertex colours, see color.frag) is kept for the night,
+// which the atmosphere does not cover (no moon, no stars), and the two are blended by the height of the sun.
+// param[0].x: share of the atmosphere, param[0].y: brightness of the atmosphere, param[0].z: how much the overcast flattens
+// the sky, param[0].w: luminance the flattened sky tends to. param[1].rgb: white balance of the atmosphere
+
+in vec3 f_color;
+in vec4 f_pos;
+
+#include <common>
+#include <atmosphere_common.glsl>
+#include <apply_fog.glsl>
+#include <tonemapping.glsl>
+#include <sky_shoulder.glsl>
+
+layout(binding = SKY_VIEW_TEX) uniform sampler2D sky_view_lut;
+
+layout(location = 0) out vec4 out_color;
+#if MOTIONBLUR_ENABLED
+layout(location = 1) out vec4 out_motion;
+#endif
+
+// the brightest the sky gets, a guard against the circumsolar peak of the mie scattering blowing up the bloom
+const float SKY_RADIANCE_MAX = 64.0;
+
+void main()
+{
+	vec3 domecolor = sky_shoulder(pow(max(f_color.rgb, vec3(0.0)), vec3(2.2)));
+
+	// the dome is drawn around the camera with the rotation of the view only, so the view space position is the direction
+	vec3 direction = normalize(mat3(inv_view) * f_pos.xyz);
+	// the bottom of the dome reaches a little under the horizon, where the table holds the dark air in front of the ground;
+	// it gets the light of the horizon instead, as the old dome did
+	direction = normalize(vec3(direction.x, max(direction.y, 0.002), direction.z));
+	vec3 sundirection = normalize(atmo_sun_direction.xyz);
+	vec3 atmosphere = textureLod(sky_view_lut, atmo_skyview_uv_world(direction, sundirection), 0.0).rgb;
+	atmosphere *= param[1].rgb * param[0].y;
+	// an overcast sky is greyer and more even than the clear one the atmosphere describes
+	vec3 overcast = vec3(param[0].w) * (0.75 + 0.5 * clamp(direction.y, 0.0, 1.0));
+	atmosphere = mix(atmosphere, overcast, param[0].z);
+	atmosphere = min(atmosphere, vec3(SKY_RADIANCE_MAX));
+
+	vec3 col = mix(domecolor, atmosphere, param[0].x);
+	if (any(isnan(col)) || any(isinf(col)))
+		col = domecolor;
+#if POSTFX_ENABLED
+	out_color = vec4(apply_fog(col), 1.0f);
+#else
+	out_color = tonemap(vec4(apply_fog(col), 1.0f));
+#endif
+#if MOTIONBLUR_ENABLED
+	out_motion = vec4(0.0f);
+#endif
+}

shaders/sky_shoulder.glsl

@@ -0,0 +1,16 @@
+// the skydome vertex colours carry Perez luminance scaled by a fixed exposure
+// (see CSkyDome::SetExposure), so the circumsolar region and the horizon band
+// routinely land well above 1.0. hard-clipping them, as this shader used to do,
+// flattened everything within ~20 degrees of the sun - and most of the lower sky -
+// into one solid white patch. compress the top end instead: below the knee the
+// response is unchanged, above it it rolls off towards SKY_SHOULDER_WHITE, so the
+// brightness gradient around the sun survives and the patch keeps a soft edge.
+const float SKY_SHOULDER_KNEE = 0.25;
+const float SKY_SHOULDER_WHITE = 0.95;
+
+vec3 sky_shoulder(vec3 color)
+{
+	const float range = SKY_SHOULDER_WHITE - SKY_SHOULDER_KNEE;
+	vec3 excess = max(color - vec3(SKY_SHOULDER_KNEE), vec3(0.0));
+	return min(color, vec3(SKY_SHOULDER_KNEE)) + range * (vec3(1.0) - exp(-excess / range));
+}