Fix the screen space ambient occlusion and add ground truth ambient occlusion to the GL 4.6 renderer

shaders/gtao_denoise.comp

@@ -0,0 +1,50 @@
+// smooths the ambient occlusion of gtao_main.comp: the average of the 4x4 pixels around each one, which takes out the
+// repeating 4x4 pattern of the slices, with the pixels at a clearly different distance left out so the occlusion does
+// not leak across the edges of objects. one invocation per pixel of the half size target
+
+layout(local_size_x = 8, local_size_y = 8) in;
+
+layout(binding = 0) uniform sampler2D depth_tex;
+layout(binding = 1) uniform sampler2D ao_tex;
+layout(binding = 0, r8) writeonly uniform image2D ao_target;
+
+layout(location = 0) uniform mat4 inv_projection;
+
+// relative difference of the distance at which a neighbour stops counting
+const float DEPTH_TOLERANCE = 0.05;
+
+float view_distance(ivec2 texel)
+{
+	ivec2 full_size = textureSize(depth_tex, 0);
+	float depth = texelFetch(depth_tex, min(texel * 2, full_size - 1), 0).r;
+	if (depth <= 0.0)
+		return 1e30;
+	vec2 uv = (vec2(texel * 2) + 0.5) / vec2(full_size);
+	vec4 position = inv_projection * vec4(uv * 2.0 - 1.0, depth, 1.0);
+	return length(position.xyz / position.w);
+}
+
+void main()
+{
+	ivec2 texel = ivec2(gl_GlobalInvocationID.xy);
+	ivec2 size = imageSize(ao_target);
+	if (texel.x >= size.x || texel.y >= size.y)
+		return;
+
+	float distance = view_distance(texel);
+	float sum = 0.0;
+	float weights = 0.0;
+	for (int y = -2; y < 2; ++y)
+	{
+		for (int x = -2; x < 2; ++x)
+		{
+			ivec2 sampletexel = clamp(texel + ivec2(x, y), ivec2(0), size - 1);
+			float sampledistance = view_distance(sampletexel);
+			float weight = clamp(1.0 - abs(sampledistance - distance) / (DEPTH_TOLERANCE * distance), 0.0, 1.0);
+			sum += texelFetch(ao_tex, sampletexel, 0).r * weight;
+			weights += weight;
+		}
+	}
+	// the pixel itself is in the window with the full weight, so weights >= 1
+	imageStore(ao_target, texel, vec4(sum / weights));
+}

shaders/gtao_main.comp

@@ -0,0 +1,152 @@
+// ground truth ambient occlusion (Jimenez et al. 2016, "Practical realtime strategies for accurate indirect occlusion";
+// integration as in Intel XeGTAO) at half the resolution of the image. for every pixel, a few slices through the
+// hemisphere around the view direction; in each the horizons on both sides are searched in the depth buffer and the
+// cosine weighted visible arc between them is integrated against the normal rebuilt from the depth.
+// one invocation per pixel of the half size target. the pattern of the slice directions and steps repeats every 4x4
+// pixels; gtao_denoise.comp averages it out
+
+layout(local_size_x = 8, local_size_y = 8) in;
+
+// depth of the scene at full resolution: reversed, 1 at the near plane, 0 at infinity (the sky)
+layout(binding = 0) uniform sampler2D depth_tex;
+layout(binding = 0, r8) writeonly uniform image2D ao_target;
+
+layout(location = 0) uniform mat4 inv_projection;
+// x: radius in metres, y: projection[0][0] times half the width of the target in pixels
+layout(location = 1) uniform vec2 params;
+// the row of the projection which gives the clip w of a view space position
+layout(location = 2) uniform vec4 clipw_row;
+// x: number of slices, y: steps on each side of a slice
+layout(location = 3) uniform ivec2 quality;
+
+const float PI = 3.14159265;
+const float HALF_PI = 1.57079633;
+// share of the radius over which the weight of an occluder falls to nothing
+const float FALLOFF_RANGE = 0.615;
+
+// depth of the pixel of the half size target: the top left one of the four full size pixels it covers, so the
+// application (postfx_ssao_apply.frag) finds the same depth for it
+float half_depth(ivec2 texel)
+{
+	ivec2 full_size = textureSize(depth_tex, 0);
+	return texelFetch(depth_tex, min(texel * 2, full_size - 1), 0).r;
+}
+
+vec3 view_position(ivec2 texel, float depth)
+{
+	ivec2 full_size = textureSize(depth_tex, 0);
+	vec2 uv = (vec2(texel * 2) + 0.5) / vec2(full_size);
+	vec4 position = inv_projection * vec4(uv * 2.0 - 1.0, depth, 1.0);
+	return position.xyz / position.w;
+}
+
+// 4x4 ordered pattern, values 0..15
+const float BAYER[16] = float[](0.0, 8.0, 2.0, 10.0, 12.0, 4.0, 14.0, 6.0, 3.0, 11.0, 1.0, 9.0, 15.0, 7.0, 13.0, 5.0);
+
+void main()
+{
+	ivec2 texel = ivec2(gl_GlobalInvocationID.xy);
+	ivec2 size = imageSize(ao_target);
+	if (texel.x >= size.x || texel.y >= size.y)
+		return;
+
+	float depth = half_depth(texel);
+	if (depth <= 0.0)
+	{
+		imageStore(ao_target, texel, vec4(1.0));
+		return;
+	}
+	vec3 position = view_position(texel, depth);
+	vec3 view = normalize(-position);
+
+	// normal from the neighbours: on each axis the one closer in depth, so edges of objects do not bend it
+	vec3 left = view_position(texel - ivec2(1, 0), half_depth(texel - ivec2(1, 0)));
+	vec3 right = view_position(texel + ivec2(1, 0), half_depth(texel + ivec2(1, 0)));
+	vec3 down = view_position(texel - ivec2(0, 1), half_depth(texel - ivec2(0, 1)));
+	vec3 up = view_position(texel + ivec2(0, 1), half_depth(texel + ivec2(0, 1)));
+	vec3 dx = (abs(right.z - position.z) < abs(position.z - left.z)) ? right - position : position - left;
+	vec3 dy = (abs(up.z - position.z) < abs(position.z - down.z)) ? up - position : position - down;
+	vec3 normal = normalize(cross(dx, dy));
+	if (dot(normal, view) < 0.0)
+		normal = -normal;
+
+	float radius = params.x;
+	float clipw = dot(clipw_row, vec4(position, 1.0));
+	float screenradius = params.y * radius / max(clipw, 1e-6);
+	// smaller than a pixel: nothing to find
+	if (screenradius < 1.0)
+	{
+		imageStore(ao_target, texel, vec4(1.0));
+		return;
+	}
+	screenradius = min(screenradius, 256.0);
+
+	float falloffrange = FALLOFF_RANGE * radius;
+	float falloffmul = -1.0 / falloffrange;
+	float falloffadd = (radius - falloffrange) / falloffrange + 1.0;
+
+	int pattern = (texel.x & 3) + 4 * (texel.y & 3);
+	float slicenoise = (BAYER[pattern] + 0.5) / 16.0;
+	float stepnoise = fract(float(pattern) * 0.618034 + 0.25);
+
+	int slices = max(quality.x, 1);
+	int steps = max(quality.y, 1);
+	float visibility = 0.0;
+	for (int slice = 0; slice < slices; ++slice)
+	{
+		float angle = (float(slice) + slicenoise) * PI / float(slices);
+		vec2 omega = vec2(cos(angle), sin(angle));
+		vec3 direction = vec3(omega, 0.0);
+		vec3 orthodirection = direction - dot(direction, view) * view;
+		vec3 axis = normalize(cross(orthodirection, view));
+		vec3 projectednormal = normal - axis * dot(normal, axis);
+		float projectedlength = length(projectednormal);
+		if (projectedlength < 1e-4)
+			continue;
+		float signnormal = sign(dot(orthodirection, projectednormal));
+		float cosnormal = clamp(dot(projectednormal, view) / projectedlength, 0.0, 1.0);
+		float n = signnormal * acos(cosnormal);
+
+		float lowcos0 = cos(n + HALF_PI);
+		float lowcos1 = cos(n - HALF_PI);
+		float horizoncos0 = lowcos0;
+		float horizoncos1 = lowcos1;
+
+		for (int step = 0; step < steps; ++step)
+		{
+			// denser near the pixel, where the occluders matter most
+			float s = (float(step) + stepnoise) / float(steps);
+			s = s * s;
+			vec2 offset = omega * max(s * screenradius, float(step) + 1.0);
+			for (int side = 0; side < 2; ++side)
+			{
+				ivec2 sampletexel = texel + ivec2(round(side == 0 ? offset : -offset));
+				if (any(lessThan(sampletexel, ivec2(0))) || any(greaterThanEqual(sampletexel, size)))
+					continue;
+				float sampledepth = half_depth(sampletexel);
+				if (sampledepth <= 0.0)
+					continue;
+				vec3 delta = view_position(sampletexel, sampledepth) - position;
+				float distance = length(delta);
+				if (distance < 1e-5)
+					continue;
+				float horizoncos = dot(delta / distance, view);
+				float weight = clamp(distance * falloffmul + falloffadd, 0.0, 1.0);
+				if (side == 0)
+					horizoncos0 = max(horizoncos0, mix(lowcos0, horizoncos, weight));
+				else
+					horizoncos1 = max(horizoncos1, mix(lowcos1, horizoncos, weight));
+			}
+		}
+
+		float h0 = -acos(clamp(horizoncos1, -1.0, 1.0));
+		float h1 = acos(clamp(horizoncos0, -1.0, 1.0));
+		h0 = n + clamp(h0 - n, -HALF_PI, HALF_PI);
+		h1 = n + clamp(h1 - n, -HALF_PI, HALF_PI);
+		float arc0 = (cosnormal + 2.0 * h0 * sin(n) - cos(2.0 * h0 - n)) * 0.25;
+		float arc1 = (cosnormal + 2.0 * h1 * sin(n) - cos(2.0 * h1 - n)) * 0.25;
+		visibility += projectedlength * (arc0 + arc1);
+	}
+	visibility = clamp(visibility / float(slices), 0.0, 1.0);
+	imageStore(ao_target, texel, vec4(visibility));
+}

shaders/postfx_ssao.frag

@@ -9,8 +9,9 @@ uniform sampler2D noise_tex;
 
 #include <common>
 
+// the depth is reversed (glClipControl GL_ZERO_TO_ONE, 1 at the near plane, 0 at infinity), so it is the clip z itself
 vec3 view_pos_from_depth(vec2 uv, float depth) {
-    vec4 clip = vec4(uv * 2.0 - 1.0, depth * 2.0 - 1.0, 1.0);
+    vec4 clip = vec4(uv * 2.0 - 1.0, depth, 1.0);
     vec4 view = inverse(projection) * clip;
     return view.xyz / view.w;
 }
@@ -31,7 +32,7 @@ vec3 kernel_sample(int i) {
 
 void main() {
     float d = texture(depth_tex, f_coords).r;
-    if (d >= 1.0) { out_color = vec4(1.0); return; }   // skybox -> no occlusion
+    if (d <= 0.0) { out_color = vec4(1.0); return; }   // skybox -> no occlusion
 
     vec3 pos = view_pos_from_depth(f_coords, d);
 

shaders/postfx_ssao_apply.frag

@@ -1,22 +1,80 @@
+// applies the ambient occlusion of the half size target to the image at full size, as a multiplier: drawn with the
+// blending dst * src over the image (opengl46_ao::apply). the four nearest pixels of the occlusion are weighed by how
+// close they are and by how close their distance is to that of the pixel, so the occlusion does not spill over the
+// edges of objects. each pixel of the occlusion stands for the top left one of the four full size pixels it covers
+// params.x: 0 ground truth ambient occlusion, 1 the older ssao (its own fixed curve)
+// params.y: strength, a power of the occlusion (1 = as computed); params.z: albedo assumed by the multiple bounce term
+
 in vec2 f_coords;
 layout(location = 0) out vec4 out_color;
 
-#texture (color_tex, 0, RGB)
-uniform sampler2D color_tex;
+layout(binding = 0) uniform sampler2D depth_tex;
+layout(binding = 1) uniform sampler2D ao_tex;
 
-#texture (ssao_tex, 1, R)
-uniform sampler2D ssao_tex;
+layout(location = 0) uniform mat4 inv_projection;
+layout(location = 1) uniform vec4 params;
 
-const float SSAO_STRENGTH = 0.8;
-const float SSAO_MIN      = 0.3;   // floor on darkening so shadows don't go pitch black
+float view_distance(ivec2 pixel, float depth)
+{
+	ivec2 full_size = textureSize(depth_tex, 0);
+	vec2 uv = (vec2(pixel) + 0.5) / vec2(full_size);
+	vec4 position = inv_projection * vec4(uv * 2.0 - 1.0, depth, 1.0);
+	return length(position.xyz / position.w);
+}
 
-void main() {
-    vec3  c   = texture(color_tex, f_coords).rgb;
-    float occ = texture(ssao_tex, f_coords).r;
+// light which comes back after more than one bounce off surfaces of the given albedo, fit of Jimenez et al. 2016
+float multi_bounce(float visibility, float albedo)
+{
+	float a = 2.0404 * albedo - 0.3324;
+	float b = -4.7951 * albedo + 0.6417;
+	float c = 2.7552 * albedo + 0.6903;
+	return max(visibility, ((visibility * a + b) * visibility + c) * visibility);
+}
 
-    // remap occlusion with strength and a minimum floor
-    occ = pow(clamp(occ, 0.0, 1.0), SSAO_STRENGTH);
-    occ = mix(SSAO_MIN, 1.0, occ);
+void main()
+{
+	ivec2 full_size = textureSize(depth_tex, 0);
+	ivec2 half_size = textureSize(ao_tex, 0);
+	ivec2 pixel = ivec2(gl_FragCoord.xy);
+	float depth = texelFetch(depth_tex, pixel, 0).r;
+	if (depth <= 0.0)
+	{
+		// the sky
+		out_color = vec4(1.0);
+		return;
+	}
+	float distance = view_distance(pixel, depth);
 
-    out_color = vec4(c * occ, 1.0);
-}
+	vec2 position = vec2(pixel) * 0.5;
+	ivec2 base = ivec2(floor(position));
+	vec2 fraction = position - vec2(base);
+	float sum = 0.0;
+	float weights = 0.0;
+	for (int y = 0; y < 2; ++y)
+	{
+		for (int x = 0; x < 2; ++x)
+		{
+			ivec2 texel = clamp(base + ivec2(x, y), ivec2(0), half_size - 1);
+			ivec2 texelpixel = min(texel * 2, full_size - 1);
+			float texeldepth = texelFetch(depth_tex, texelpixel, 0).r;
+			float texeldistance = (texeldepth > 0.0) ? view_distance(texelpixel, texeldepth) : 1e30;
+			float bilinear = (x == 0 ? 1.0 - fraction.x : fraction.x) * (y == 0 ? 1.0 - fraction.y : fraction.y);
+			float weight = (bilinear + 1e-3) / (1e-3 + abs(texeldistance - distance) / distance);
+			sum += texelFetch(ao_tex, texel, 0).r * weight;
+			weights += weight;
+		}
+	}
+	float ao = clamp(sum / max(weights, 1e-6), 0.0, 1.0);
+
+	if (params.x > 0.5)
+	{
+		// the curve the older ssao was applied with, a floor so the shadows do not go black
+		ao = mix(0.3, 1.0, pow(ao, 0.8));
+	}
+	else
+	{
+		ao = pow(ao, params.y);
+		ao = multi_bounce(ao, params.z);
+	}
+	out_color = vec4(vec3(ao), 1.0);
+}