#include <psemek/app/app.hpp>
#include <psemek/app/main.hpp>

#include <psemek/gfx/gl.hpp>
#include <psemek/gfx/mesh.hpp>
#include <psemek/gfx/program.hpp>
#include <psemek/gfx/painter.hpp>

#include <psemek/geom/camera.hpp>
#include <psemek/geom/math.hpp>
#include <psemek/geom/homogeneous.hpp>
#include <psemek/geom/gauss.hpp>
#include <psemek/geom/orientation.hpp>
#include <psemek/geom/intersection.hpp>
#include <psemek/geom/distance.hpp>

#include <psemek/util/clock.hpp>
#include <psemek/util/to_string.hpp>
#include <psemek/util/moving_average.hpp>
#include <psemek/util/recursive.hpp>

namespace psemek::cg
{

	template <typename T, std::size_t N>
	struct box;

	template <typename T, std::size_t N>
	box(geom::box<T, N>) -> box<T, N>;

	template <typename T>
	struct box<T, 3>
	{
		box() = default;
		box(geom::box<T, 3> const & b);

		std::array<geom::point<T, 3>, 8> vertices;
	};

	template <typename T>
	box<T, 3>::box(geom::box<T, 3> const & b)
	{
		for (std::size_t z = 0; z < 2; ++z)
		{
			for (std::size_t y = 0; y < 2; ++y)
			{
				for (std::size_t x = 0; x < 2; ++x)
				{
					std::size_t i = z * 4 + y * 2 + x;

					vertices[i][0] = (x == 0) ? b[0].min : b[0].max;
					vertices[i][1] = (y == 0) ? b[1].min : b[1].max;
					vertices[i][2] = (z == 0) ? b[2].min : b[2].max;
				}
			}
		}
	}

	template <typename T>
	auto const & vertices(box<T, 3> const & b)
	{
		return b.vertices;
	}

	namespace detail
	{

		inline auto const & cubiod_edges()
		{
			static const std::array<geom::segment<std::uint8_t>, 12> result =
			{{
				{ 0b000, 0b001 },
				{ 0b010, 0b011 },
				{ 0b100, 0b101 },
				{ 0b110, 0b111 },

				{ 0b000, 0b010 },
				{ 0b001, 0b011 },
				{ 0b100, 0b110 },
				{ 0b101, 0b111 },

				{ 0b000, 0b100 },
				{ 0b001, 0b101 },
				{ 0b010, 0b110 },
				{ 0b011, 0b111 },
			}};

			return result;
		}

		inline auto const & cubiod_faces()
		{
			static const std::array<std::array<std::uint8_t, 4>, 6> result =
			{{
				{{ 0b000, 0b100, 0b110, 0b010 }},
				{{ 0b001, 0b011, 0b111, 0b101 }},
				{{ 0b000, 0b001, 0b101, 0b100 }},
				{{ 0b010, 0b110, 0b111, 0b011 }},
				{{ 0b000, 0b010, 0b011, 0b001 }},
				{{ 0b100, 0b101, 0b111, 0b110 }},
			}};

			return result;
		}
	}

	template <typename T>
	auto const & edges(box<T, 3> const &)
	{
		return detail::cubiod_edges();
	}

	template <typename T>
	auto const & faces(box<T, 3> const &)
	{
		return detail::cubiod_faces();
	}

	template <typename T>
	auto const & face_normals(box<T, 3> const &)
	{
		static const std::array<geom::vector<T, 3>, 6> result =
		{{
			{-1,  0,  0},
			{ 1,  0,  0},
			{ 0, -1,  0},
			{ 0,  1,  0},
			{ 0,  0, -1},
			{ 0,  0,  1},
		}};

		return result;
	}

	template <typename T>
	auto const & edge_directions(box<T, 3> const &)
	{
		static const std::array<geom::vector<T, 3>, 3> result =
		{{
			{ 1,  0,  0},
			{ 0,  1,  0},
			{ 0,  0,  1},
		}};

		return result;
	}

	template <typename T, std::size_t N>
	struct frustum;

	template <typename T, std::size_t N>
	frustum(geom::matrix<T, N, N>) -> frustum<T, N-1>;

	template <typename T>
	struct frustum<T, 3>
	{
		frustum(geom::matrix<T, 4, 4> const & m);

		std::array<geom::point<T, 3>, 8> vertices;
		std::array<geom::vector<T, 3>, 6> face_normals;
	};

	template <typename T>
	frustum<T, 3>::frustum(geom::matrix<T, 4, 4> const & m)
	{
		bool flip = (geom::det(m) < 0);

		for (std::size_t z = 0; z < 2; ++z)
		{
			for (std::size_t y = 0; y < 2; ++y)
			{
				for (std::size_t x = 0; x < 2; ++x)
				{
					std::size_t i = z * 4 + y * 2 + (flip ? 1 - x : x);

					geom::vector<T, 4> p;
					p[0] = (x == 0) ? -1 : 1;
					p[1] = (y == 0) ? -1 : 1;
					p[2] = (z == 0) ? -1 : 1;
					p[3] = 1;

					geom::gauss(m, p);

					vertices[i] = geom::as_point(p);
				}
			}
		}
	}

	template <typename T>
	auto const & vertices(frustum<T, 3> const & f)
	{
		return f.vertices;
	}

	template <typename T>
	auto const & edges(frustum<T, 3> const &)
	{
		return detail::cubiod_edges();
	}

	template <typename T>
	auto const & faces(frustum<T, 3> const &)
	{
		return detail::cubiod_faces();
	}

	template <typename T>
	struct triangular_prism
	{
		triangular_prism() = default;
		triangular_prism(geom::triangle<geom::point<T, 3>> const & t, geom::vector<T, 3> const & d);

		std::array<geom::point<T, 3>, 6> vertices;
		std::array<geom::vector<T, 3>, 4> edge_directions;
	};

	template <typename T>
	triangular_prism<T>::triangular_prism(geom::triangle<geom::point<T, 3>> const & t, geom::vector<T, 3> const & d)
	{
		for (std::size_t i = 0; i < 3; ++i)
		{
			vertices[i] = t[i];
			vertices[i + 3] = t[i] + d;
		}

		if (geom::dot(geom::normal(vertices[0], vertices[1], vertices[2]), d) < 0)
		{
			std::swap(vertices[1], vertices[2]);
			std::swap(vertices[4], vertices[5]);
		}

		for (std::size_t i = 0; i < 3; ++i)
		{
			edge_directions[i] = geom::normalized(vertices[(i + 1) % 3] - vertices[i]);
		}
		edge_directions[3] = geom::normalized(d);
	}

	template <typename T>
	auto const & vertices(triangular_prism<T> const & p)
	{
		return p.vertices;
	}

	template <typename T>
	auto const & edges(triangular_prism<T> const &)
	{
		static const std::array<geom::segment<std::uint8_t>, 9> result =
		{{
			{ 0, 1 },
			{ 1, 2 },
			{ 2, 0 },

			{ 3, 4 },
			{ 4, 5 },
			{ 5, 3 },

			{ 0, 3 },
			{ 1, 4 },
			{ 2, 5 },
		}};

		return result;
	}

	template <typename T>
	auto const & faces(triangular_prism<T> const &)
	{
		static std::array<std::vector<std::uint8_t>, 5> result =
		{{
			{ 0, 2, 1 },
			{ 3, 4, 5 },
			{ 0, 1, 4, 3 },
			{ 1, 2, 5, 4 },
			{ 2, 0, 3, 5 },
		}};
		return result;
	}

	template <typename T>
	auto const & edge_directions(triangular_prism<T> const & p)
	{
		return p.edge_directions;
	}

	namespace detail
	{

		template <typename Container>
		struct has_static_size
			: std::false_type
		{};

		template <typename T, std::size_t N>
		struct has_static_size<std::array<T, N>>
			: std::true_type
		{};

		template <typename Container>
		constexpr bool has_static_size_v = has_static_size<Container>::value;

		template <typename Container>
		struct static_size;

		template <typename T, std::size_t N>
		struct static_size<std::array<T, N>>
		{
			static constexpr std::size_t value = N;
		};

		template <typename Container>
		constexpr std::size_t static_size_v = static_size<Container>::value;

	}

	template <typename Body>
	auto triangles(Body const & b)
	{
		auto const & fs = faces(b);

		using faces_type = std::remove_cvref_t<decltype(fs)>;
		using face_type = std::remove_cvref_t<decltype(*std::begin(fs))>;
		using index_type = std::remove_cvref_t<decltype((*std::begin(fs))[0])>;

		auto impl = [&fs](auto out)
		{
			for (auto const & f : fs)
			{
				auto it0 = std::begin(f);

				for (auto it = std::next(it0), jt = std::next(it); jt < std::end(f); it = jt++)
				{
					*out++ = {*it0, *it, *jt};
				}
			}
		};

		if constexpr (detail::has_static_size_v<faces_type> && detail::has_static_size_v<face_type>)
		{
			constexpr std::size_t faces_count = detail::static_size_v<faces_type>;
			constexpr std::size_t face_size = detail::static_size_v<face_type>;
			std::array<geom::triangle<index_type>, faces_count * (face_size - 2)> result;
			impl(result.begin());
			return result;
		}
		else
		{
			std::vector<geom::triangle<index_type>> result;
			if constexpr (detail::has_static_size_v<face_type>)
			{
				constexpr std::size_t face_size = detail::static_size_v<face_type>;
				result.reserve(fs.size() * (face_size - 2));
			}
			else
			{
				result.reserve(fs.size());
			}
			impl(std::back_inserter(result));
			return result;
		}
	}

	template <typename Body>
	auto edge_directions(Body const & b)
	{
		auto const & vs = vertices(b);
		auto const & es = edges(b);

		using edges_type = std::remove_cvref_t<decltype(es)>;
		using vector_type = std::remove_cvref_t<decltype(vs[0] - vs[0])>;

		auto impl = [&es, &vs](auto out)
		{
			for (auto const & e : es)
			{
				*out++ = geom::normalized(vs[e[1]] - vs[e[0]]);
			}
		};

		if constexpr (detail::has_static_size_v<edges_type>)
		{
			constexpr std::size_t edges_count = detail::static_size_v<edges_type>;
			std::array<vector_type, edges_count> result;
			impl(result.begin());
			return result;
		}
		else
		{
			std::vector<vector_type> result;
			result.reserve(es.size());
			impl(std::back_inserter(result));
			return result;
		}
	}

	template <typename Body>
	auto faces(Body const & b)
	{
		return triangles(b);
	}

	template <typename Body>
	auto face_normals(Body const & b)
	{
		auto const & vs = vertices(b);
		auto const & fs = faces(b);

		using faces_type = std::remove_cvref_t<decltype(fs)>;
		using vector_type = std::remove_cvref_t<decltype(vs[0] - vs[0])>;

		auto impl = [&fs, &vs](auto out)
		{
			for (auto const & f : fs)
			{
				*out++ = geom::normal(vs[f[0]], vs[f[1]], vs[f[2]]);
			}
		};

		if constexpr (detail::has_static_size_v<faces_type>)
		{
			constexpr std::size_t faces_count = detail::static_size_v<faces_type>;
			std::array<vector_type, faces_count> result;
			impl(result.begin());
			return result;
		}
		else
		{
			std::vector<vector_type> result;
			result.reserve(fs.size());
			impl(std::back_inserter(result));
			return result;
		}
	}

	template <typename P, typename Body>
	bool inside(P const & p, Body const & body)
	{
		auto const & vs = vertices(body);
		auto const & fs = faces(body);

		for (auto const & f : fs)
		{
			if (geom::orientation(vs[f[0]], vs[f[1]], vs[f[2]], p) == geom::sign_t::negative)
				return false;
		}
		return true;
	}

	// returns pair(normalized vector from 1 to 2, signed distance)
	// distance <= 0 means intersection
	template <typename Body1, typename Body2>
	auto separation(Body1 const & b1, Body2 const & b2)
	{
		auto const & vs1 = vertices(b1);
		auto const & vs2 = vertices(b2);

		auto const & fs1 = faces(b1);
		auto const & fs2 = faces(b2);

		auto const & eds1 = edge_directions(b1);
		auto const & eds2 = edge_directions(b2);

		using vector_type = std::remove_cvref_t<decltype(vs1[0] - vs1[0])>;
		using scalar_type = std::remove_cvref_t<decltype(vs1[0][0])>;

		vector_type res_n = vector_type::zero();
		auto res_d = -std::numeric_limits<scalar_type>::infinity();

		auto process_faces = [](auto const & vs1, auto const & fs1, auto const & vs2)
		{
			vector_type res_n = vector_type::zero();
			scalar_type res_d = -std::numeric_limits<scalar_type>::infinity();

			for (auto const & f : fs1)
			{
				auto const face_n = geom::normal(vs1[f[0]], vs1[f[1]], vs1[f[2]]);
				scalar_type face_d = std::numeric_limits<scalar_type>::infinity();

				auto const face_p = vs1[f[0]];

				for (auto const & v : vs2)
				{
					auto const d = geom::dot(face_n, v - face_p);
					face_d = std::min(d, face_d);
				}

				if (face_d == std::numeric_limits<scalar_type>::infinity())
				{
					throw 42;
				}

				if (face_d > res_d)
				{
					res_d = face_d;
					res_n = face_n;
				}
			}

			return std::make_pair(res_n, res_d);
		};

		auto process_edges = [](auto const & vs1, auto const & eds1, auto const & vs2, auto const & eds2)
		{
			vector_type res_n = vector_type::zero();
			scalar_type res_d = -std::numeric_limits<scalar_type>::infinity();

			for (auto const & ed1 : eds1)
			{
				for (auto const & ed2 : eds2)
				{
					auto edge_n = geom::cross(ed1, ed2);
					auto l = geom::length(edge_n);
					if (l == 0) continue;
					edge_n /= l;

					geom::interval<scalar_type> i1, i2;

					for (auto const & v : vs1)
					{
						i1 |= geom::dot(geom::homogeneous(v), geom::homogeneous(edge_n));
					}

					for (auto const & v : vs2)
					{
						i2 |= geom::dot(geom::homogeneous(v), geom::homogeneous(edge_n));
					}

					auto edge_d12 = i2.min - i1.max;
					auto edge_d21 = i1.min - i2.max;

					scalar_type edge_d;

					if (edge_d12 > edge_d21)
					{
						edge_d = edge_d12;
					}
					else
					{
						edge_d = edge_d21;
						edge_n = -edge_n;
					}

					if (edge_d > res_d)
					{
						res_d = edge_d;
						res_n = edge_n;
					}
				}
			}

			return std::make_pair(res_n, res_d);
		};

		{
			auto res12 = process_faces(vs1, fs1, vs2);
			if (res12.second > res_d)
			{
				res_d = res12.second;
				res_n = res12.first;
			}
		}

		{
			auto res21 = process_faces(vs2, fs2, vs1);
			if (res21.second > res_d)
			{
				res_d = res21.second;
				res_n = -res21.first;
			}
		}

		{
			auto rese = process_edges(vs1, eds1, vs2, eds2);
			if (rese.second > res_d)
			{
				res_d = rese.second;
				res_n = rese.first;
			}
		}

		return std::make_pair(res_n, res_d);
	}

}

using namespace psemek;

template <typename Value, typename T>
struct smooth_updater
{
	smooth_updater(Value & value, T speed)
		: value_{value}
		, target_value_{value}
		, speed_{speed}
	{}

	smooth_updater & operator = (Value const & value)
	{
		target_value_ = value;
		return *this;
	}

	operator Value const & () const
	{
		return target_value_;
	}

	void update(T dt)
	{
		value_ += std::min(T{1}, dt * speed_) * (target_value_ - value_);
	}

private:
	Value & value_;
	Value target_value_;
	T speed_;
};

static const float icosa_a = 0.850650808; // phi / sqrt(1 + phi^2)
static const float icosa_b = 0.525731112; //   1 / sqrt(1 + phi^2)

static const float icosa_side = 1.05146222; // 2 / sqrt(phi * sqrt(5))

static const geom::vector<float, 3> icosa_vertices[12] =
{
	{-icosa_a, -icosa_b, 0.0}, //  0
	{-icosa_a,  icosa_b, 0.0}, //  1
	{ icosa_a, -icosa_b, 0.0}, //  2
	{ icosa_a,  icosa_b, 0.0}, //  3

	{0.0, -icosa_a, -icosa_b}, //  4
	{0.0, -icosa_a,  icosa_b}, //  5
	{0.0,  icosa_a, -icosa_b}, //  6
	{0.0,  icosa_a,  icosa_b}, //  7

	{-icosa_b, 0.0, -icosa_a}, //  8
	{ icosa_b, 0.0, -icosa_a}, //  9
	{-icosa_b, 0.0,  icosa_a}, // 10
	{ icosa_b, 0.0,  icosa_a}, // 11
};

static const std::size_t icosa_faces[20][3] =
{
	{0, 1, 8,},
	{0, 10, 1,},
	{2, 9, 3,},
	{2, 3, 11,},
	{4, 5, 0,},
	{4, 2, 5,},
	{6, 1, 7,},
	{6, 7, 3,},
	{8, 9, 4,},
	{8, 6, 9,},
	{10, 5, 11,},
	{10, 11, 7,},
	{0, 8, 4,},
	{0, 5, 10,},
	{1, 6, 8,},
	{1, 10, 7,},
	{2, 4, 9,},
	{2, 11, 5,},
	{3, 9, 6,},
	{3, 7, 11},
};

static char const tile_vs[] =
R"(#version 330

uniform mat4 u_transform;

uniform int u_icosa_face;
uniform int u_N;
uniform vec3 u_p0;
uniform vec3 u_p1;
uniform vec3 u_p2;
uniform vec4 u_color;

out vec3 color;

void main()
{
	int i = int(floor(0.5 * (sqrt(1.0 + 8.0 * gl_VertexID) - 1.0)));
	int j = gl_VertexID - (i * (i + 1)) / 2;

	float t0 = 1.0 - float(i) / float(u_N);
	float t1 = float(j) / float(u_N);
	float t2 = 1.0 - t0 - t1;

	vec3 p = normalize(u_p0 * t0 + u_p1 * t1 + u_p2 * t2);
	gl_Position = u_transform * vec4(p, 1.0);
	color = u_color.rgb;
})";

static char const tile_fs[] =
R"(#version 330

in vec3 color;
out vec4 out_color;

void main()
{
	out_color = vec4(color, 1.0);
})";

static char const test_vs[] =
R"(#version 330

uniform mat4 u_transform;

layout (location = 0) in vec4 in_position;

out vec4 color;

void main()
{
	gl_Position = u_transform * in_position;
	color = vec4(in_position.xyz * 0.5 + vec3(0.5), 1.0);
})";

static char const test_fs[] =
R"(#version 330

in vec4 color;
out vec4 out_color;

void main()
{
	out_color = color;
})";

struct srtm_app
	: app::app
{
	srtm_app();

	void on_resize(int width, int height) override;

	void on_mouse_move(int x, int y, int dx, int dy) override;

	void update() override;
	void present() override;

//	geom::spherical_camera camera;
	geom::free_camera camera;
	smooth_updater<float, float> camera_azimuthal_angle_updater{camera.azimuthal_angle, 20.f};
	smooth_updater<float, float> camera_elevation_angle_updater{camera.elevation_angle, 20.f};
	bool camera_forward = false;

	static constexpr int tile_size_log2 = 8;
	static constexpr int tile_size = (1 << tile_size_log2);

	gfx::mesh tile_mesh[tile_size_log2 + 1];

	gfx::program tile_program{tile_vs, tile_fs};
	gfx::program test_program{test_vs, test_fs};

	gfx::mesh test_mesh;

	struct test_object
	{
		cg::box<float, 3> body;
		gfx::mesh mesh;
		int x, y;
	};

	int const test_object_count = 21;
	std::vector<test_object> test_objects;

	util::clock<std::chrono::duration<float>> frame_clock;
	util::moving_average<float> frame_dt_average{32};

	gfx::painter painter;
};

srtm_app::srtm_app()
	: app("SRTM", 4)
{
	vsync(false);
	show_cursor(false);

	camera.fov_y = geom::rad(45.f);
	camera.near_clip = 0.0001f;
	camera.far_clip = 10.f;
//	camera.target = {0.f, 0.f, 0.f};
	camera.pos = {0.f, -10.f, 0.f};
	camera.azimuthal_angle = 0.f;
	camera.elevation_angle = 0.f;
//	camera.distance = 5.f;

	camera_azimuthal_angle_updater = camera.azimuthal_angle;
	camera_elevation_angle_updater = camera.elevation_angle;
//	camera_distance_updater = camera.distance;

	for (std::size_t N = 0; N <= tile_size_log2; ++N)
	{
		std::vector<std::uint16_t> indices;

		std::size_t step = tile_size / (1 << N);

		auto idx = [step](std::size_t i, std::size_t j) -> std::uint16_t
		{
			i *= step;
			j *= step;
			return (i * (i + 1)) / 2 + j;
		};

		for (std::size_t i = 0; i < (1 << N); ++i)
		{
			for (std::size_t j = 0; j <= i; ++j)
			{
				indices.push_back(idx(i + 1, j));
				indices.push_back(idx(i, j));
			}
			indices.push_back(idx(i + 1, i + 1));
			indices.push_back(0xffffu);
		}

		tile_mesh[N].load_index(indices, gl::TRIANGLE_STRIP, gl::STATIC_DRAW);
	}

	for (int x = - test_object_count + 1; x <= test_object_count - 1; x += 2)
	{
		for (int y = - test_object_count + 1; y <= test_object_count - 1; y += 2)
		{
			auto & o = test_objects.emplace_back();

			o.x = (x + test_object_count - 1) / 2;
			o.y = (y + test_object_count - 1) / 2;

			geom::box<float, 3> b {{
				{x - 1.f, x + 1.f},
				{y - 1.f, y + 1.f},
				{0 - 1.f, 0 + 1.f}
			}};

			o.body = cg::box(b);

			auto const & vertices = cg::vertices(o.body);
			auto const & edges = cg::edges(o.body);

			o.mesh.setup<geom::point<float, 3>>();
			o.mesh.load(vertices.data(), vertices.size(), edges.data(), edges.size());
		}
	}
}

void srtm_app::on_resize(int width, int height)
{
	app::on_resize(width, height);
	camera.set_fov(camera.fov_y, (1.f * width) / height);
}

void srtm_app::on_mouse_move(int x, int y, int dx, int dy)
{
	app::on_mouse_move(x, y, dx, dy);

	camera_azimuthal_angle_updater = camera_azimuthal_angle_updater - 0.01f * dx;
	camera_elevation_angle_updater = camera_elevation_angle_updater + 0.01f * dy;
}

void srtm_app::update()
{
	float dt = frame_clock.restart().count();
	frame_dt_average.push(dt);

	if (is_key_down(SDLK_q))
	{
	}

	if (is_key_down(SDLK_e))
	{
	}

	camera_azimuthal_angle_updater.update(dt);
	camera_elevation_angle_updater.update(dt);

	float const camera_speed = std::min(5.f, geom::distance(camera.pos, geom::point<float, 3>::zero()) - 1.f);
	auto const camera_forward = camera.direction();

	auto const camera_up = camera.axis_y();
	auto const camera_right = camera.axis_x();

	if (is_key_down(SDLK_w))
	{
		camera.pos += camera_speed * dt * camera_forward;
	}

	if (is_key_down(SDLK_s))
	{
		camera.pos -= camera_speed * dt * camera_forward;
	}

	if (is_key_down(SDLK_d))
	{
		camera.pos += camera_speed * dt * camera_right;
	}

	if (is_key_down(SDLK_a))
	{
		camera.pos -= camera_speed * dt * camera_right;
	}

	if (is_key_down(SDLK_LSHIFT))
	{
		camera.pos += camera_speed * dt * camera_up;
	}

	if (is_key_down(SDLK_LCTRL))
	{
		camera.pos -= camera_speed * dt * camera_up;
	}
}

namespace std
{

template <typename T>
std::ostream & operator << (std::ostream & os, std::vector<T> const & v)
{
	os << "[";
	bool first = true;
	for (auto const & x : v)
	{
		if (first)
			first = false;
		else
			os << ", ";
		os << x;
	}
	return os << "]";
}

}

void srtm_app::present()
{
	std::vector<std::string> info;

	gl::ClearColor(0.9f, 0.9f, 0.9f, 0.f);
	gl::Clear(gl::COLOR_BUFFER_BIT | gl::DEPTH_BUFFER_BIT);

	gl::LineWidth(2.f);
	gl::PolygonMode(gl::FRONT_AND_BACK, gl::LINE);
	gl::PointSize(5.f);
	gl::Enable(gl::CULL_FACE);
	gl::Enable(gl::DEPTH_TEST);
	gl::DepthFunc(gl::LEQUAL);

	gl::Enable(gl::PRIMITIVE_RESTART);
	gl::PrimitiveRestartIndex(0xffffu);

	{
		auto d = geom::distance(camera.pos, geom::point<float, 3>::zero());
		camera.far_clip = std::sqrt(d * d + 1.f);
		camera.near_clip = (d > 2.f) ? d - 2.f : 0.0001f;
	}

	auto const camera_transform = camera.transform();
	auto const camera_pos = camera.position();
	auto const camera_direction = camera.direction();

	info.push_back(util::to_string("Camera height: ", (geom::distance(camera_pos, geom::point<float, 3>::zero()) - 1.f) * 6400000.f, " m"));

	auto const frustum = cg::frustum(camera_transform);

	tile_program.bind();
	tile_program["u_transform"] = camera_transform;
	tile_program["u_N"] = static_cast<int>(tile_size);

	info.push_back(util::to_string("Camera pos: ", camera_pos));

	std::size_t rendered_tiles = 0;
	std::vector<std::size_t> id;
	auto visit = util::recursive([&](auto & self, geom::vector<float, 3> const (&v)[3], int level = 0) -> void
	{
		auto const o = geom::point<float, 3>::zero();
		auto m = (v[0] + v[1] + v[2]) / 3.f;
		m = geom::normalized(m) * (1.f + 1.f / 6400.f) - m;

		{
			bool culled = true;
			for (std::size_t i = 0; i < 3; ++i)
			{
				if (geom::dot(v[i], geom::normalized(camera_pos - o)) >= 0.f)
				{
					culled = false;
					break;
				}
			}
			if (culled)
				return;
			(void)culled;
		}

		geom::triangle<geom::point<float, 3>> t{o + v[0], o + v[1], o + v[2]};
		cg::triangular_prism<float> body{t, m};

		if (cg::separation(body, frustum).second > 0.f)
			return;

		bool const selected = geom::intersect(geom::ray{camera_pos, camera_direction}, t);

		float on_screen_unit;

		{
			auto edge = [](auto const & v0, auto const & v1, auto const & u) -> std::optional<float>
			{
				auto const n = geom::normalized(geom::cross(v0, v1));
				auto v = geom::normalized(u - n * dot(u, n));
				if (geom::dot(geom::cross(v0, v), geom::cross(v, v1)) >= 0.f)
					return geom::length(v - u);
				return std::nullopt;
			};

			float distance = std::numeric_limits<float>::infinity();
			auto c = camera_pos - o;
			if (geom::det(v[0], v[1], c) >= 0.f && geom::det(v[1], v[2], c) >= 0.f && geom::det(v[2], v[0], c) >= 0.f)
				distance = std::min(distance, geom::length(c) - 1.f);
			if (auto d = edge(v[0], v[1], c); d)
				distance = std::min(distance, *d);
			if (auto d = edge(v[1], v[2], c); d)
				distance = std::min(distance, *d);
			if (auto d = edge(v[2], v[0], c); d)
				distance = std::min(distance, *d);
			distance = std::min(distance, geom::length(c - v[0]));
			distance = std::min(distance, geom::length(c - v[1]));
			distance = std::min(distance, geom::length(c - v[2]));

			on_screen_unit = width() / distance / std::tan(camera.fov_x / 2.f);
		}

		assert(on_screen_unit > 0.f);

		float const max_triangle_size = 10.f; // pixels

		float const side_length = icosa_side / (1 << (level * 4));

		int tile_n = std::ceil(std::log2(on_screen_unit * side_length / max_triangle_size));

		if (level < 3 && tile_n > tile_size_log2)
		{
			std::size_t const C = 16;

			auto at = [&](std::size_t i, std::size_t j)
			{
				float t0 = 1.f - (1.f * i) / C;
				float t1 = (1.f * j) / C;
				float t2 = 1.f - t0 - t1;

				return geom::normalized(v[0] * t0 + v[1] * t1 + v[2] * t2);
			};

			std::size_t child_id = 0;

			auto child = [&](std::size_t i0, std::size_t j0, std::size_t i1, std::size_t j1, std::size_t i2, std::size_t j2)
			{
				geom::vector<float, 3> v[3];
				v[0] = at(i0, j0);
				v[1] = at(i1, j1);
				v[2] = at(i2, j2);
				id.push_back(child_id++);
				self(v, level + 1);
				id.pop_back();
			};

			for (std::size_t i = 0; i < C; ++i)
			{
				for (std::size_t j = 0; j < i; ++j)
				{
					child(i + 1, j, i, j, i + 1, j + 1);
					child(i + 1, j + 1, i, j, i, j + 1);
				}
				child(i + 1, i, i, i, i + 1, i + 1);
			}
		}
		else
		{
			tile_n = geom::clamp(tile_n, {0, tile_size_log2});

			++rendered_tiles;

			static gfx::color_4f colors[4]
			{
				gfx::black,
				gfx::dark(gfx::red).as_color_4f(),
				gfx::dark(gfx::green).as_color_4f(),
				gfx::blue
			};

			tile_program["u_p0"] = v[0];
			tile_program["u_p1"] = v[1];
			tile_program["u_p2"] = v[2];
			tile_program["u_color"] = colors[tile_n % 4];

			tile_mesh[tile_n].draw();

			if (selected)
			{
				gl::Disable(gl::DEPTH_TEST);
				tile_program["u_color"] = gfx::white.as_color_4f();
				tile_mesh[0].draw();
				gl::Enable(gl::DEPTH_TEST);

				info.push_back(util::to_string("Selected id: ", id));
				info.push_back(util::to_string("Selected separation: ", cg::separation(body, frustum).second));
			}
		}
	});

//	if(false)
//	for (int f = 4; f < 5; ++f)
	for (int f = 0; f < 20; ++f)
	{
		geom::vector<float, 3> v[3];
		v[0] = icosa_vertices[icosa_faces[f][0]];
		v[1] = icosa_vertices[icosa_faces[f][1]];
		v[2] = icosa_vertices[icosa_faces[f][2]];
		id.push_back(f);
		visit(v);
		id.pop_back();
	}

	info.push_back(util::to_string("Tiles: ", rendered_tiles));

	if (false)
	{
		gl::PolygonMode(gl::FRONT_AND_BACK, gl::FILL);
		gl::Disable(gl::CULL_FACE);

		test_program.bind();
		test_program["u_transform"] = camera_transform;
		int visible_count = 0;
		for (auto const & o : test_objects)
		{
			gfx::color_rgba color;

			if (cg::separation(o.body, frustum).second < 0.f)
			{
				++visible_count;
				o.mesh.draw();

				color = gfx::red;
			}
			else
			{
				color = gfx::black;
			}

			float sx = width() - 10 * test_object_count + o.x * 10;
			float sy = height() - 10 * test_object_count + (test_object_count - 1 - o.y) * 10;

			painter.rect({{{sx, sx + 10.f}, {sy, sy + 10.f}}}, color);
		}

		info.push_back(util::to_string("Visible count: ", visible_count));
	}

	{
		float s = 10.f;
		painter.line({width() / 2.f - s, height() / 2.f}, {width() / 2.f + s, height() / 2.f}, 3.f, gfx::cyan, false);
		painter.line({width() / 2.f, height() / 2.f - s}, {width() / 2.f, height() / 2.f + s}, 3.f, gfx::cyan, false);
	}

	{
		info.insert(info.begin(), util::to_string("FPS: ", 1.f / frame_dt_average.average()));

		gfx::painter::text_options opts;
		opts.x = gfx::painter::x_align::left;
		opts.y = gfx::painter::y_align::top;
		opts.f = gfx::painter::font::font_9x12;
		opts.c = gfx::cyan;
		opts.scale = 2.f;

		for (int l = 0; l < info.size(); ++l)
		{
			painter.text({10.f + 1.f, 10.f + 24.f * l + 1.f}, info[l], opts);
		}

		opts.c = gfx::black;

		for (int l = 0; l < info.size(); ++l)
		{
			painter.text({10.f, 10.f + 24.f * l}, info[l], opts);
		}
	}

	gl::PolygonMode(gl::FRONT_AND_BACK, gl::FILL);
	gl::Enable(gl::BLEND);
	gl::BlendFunc(gl::SRC_ALPHA, gl::ONE_MINUS_SRC_ALPHA);
	gl::Disable(gl::DEPTH_TEST);
	painter.render(geom::window_camera{width(), height()}.transform());
}

int main()
{
	return app::main<srtm_app>();
}