#if !defined INCLUDE_SKY_PROJECTION
#define INCLUDE_SKY_PROJECTION

#include "/include/utility/fast_math.glsl"

// Sky map projection from https://sebh.github.io/publications/egsr2020.pdf

const ivec2 sky_map_res = ivec2(191, 108);
const vec2  sky_map_pixel_size = rcp(vec2(sky_map_res));

vec2 project_sky(vec3 direction) {
	vec2 projected_dir = normalize(direction.xz);

	float azimuth_angle = pi + atan(projected_dir.x, -projected_dir.y);
	float altitude_angle = half_pi - fast_acos(direction.y);

	vec2 coord;
	coord.x = azimuth_angle * (1.0 / tau);
	coord.y = 0.5 + 0.5 * sign(altitude_angle) * sqrt(2.0 * rcp_pi * abs(altitude_angle)); // Section 5.3

	// Padding
	const float pad_amount = 2.0;
	const float mul = 1.0 - 2.0 * pad_amount * sky_map_pixel_size.x;
	const float add = pad_amount * sky_map_pixel_size.x;
	coord.x = coord.x * mul + add;

	coord.x *= 191.0 / 192.0;

	return coord;
}

vec3 unproject_sky(vec2 coord) {
	coord.x *= 192.0 / 191.0;

	// Padding
	const float pad_amount = 2.0;
	const float mul = rcp(1.0 - 2.0 * pad_amount * sky_map_pixel_size.x);
	const float add = (pad_amount * sky_map_pixel_size.x) * -mul;
	coord.x = fract(coord.x * mul + add);

	// Non-linear mapping of altitude angle (See section 5.3 of the paper)
	coord.y = (coord.y < 0.5)
		? -sqr(1.0 - 2.0 * coord.y)
		:  sqr(2.0 * coord.y - 1.0);

	float azimuth_angle = coord.x * tau - pi;
	float altitude_angle = coord.y * half_pi;

	float altitude_cos = cos(altitude_angle);
	float altitude_sin = sin(altitude_angle);
	float azimuth_cos = cos(azimuth_angle);
	float azimuth_sin = sin(azimuth_angle);

	return vec3(altitude_cos * azimuth_sin, altitude_sin, -altitude_cos * azimuth_cos);
}

#endif // INCLUDE_SKY_PROJECTION