psemek/examples/srtm.cpp

#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/gfx/error.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/cg/body/icosahedron.hpp>
#include <psemek/cg/body/box.hpp>
#include <psemek/cg/body/frustum.hpp>
#include <psemek/cg/body/prism.hpp>
#include <psemek/cg/convex/inside.hpp>
#include <psemek/cg/convex/separation.hpp>

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

#include <fstream>
#include <iomanip>
#include <atomic>
#include <unordered_map>
#include <unordered_set>

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_;
};

struct height_provider
{
	float height_at(geom::vector<float, 3> const & v);

	struct datum_id_hash
	{
		std::size_t operator() (std::pair<int, int> const & p) const
		{
			return ((p.first + 180) << 16) | (p.second + 90);
		}
	};

	std::unordered_map<std::pair<int, int>, std::unique_ptr<std::int16_t[]>, datum_id_hash> datums;
	std::unordered_set<std::pair<int, int>, datum_id_hash> no_datums;
};

static geom::interval<int> height_range;

float height_provider::height_at(geom::vector<float, 3> const & v)
{
	static std::string const data_path = "/home/lisyarus/data/srtm/dem-unpacked/";

	float const lat = geom::deg(std::asin(v[2]));

	if (std::abs(lat) > 60.f) return 0.f;
	return -10.f;

	float const lon = geom::deg(std::atan2(v[0], -v[1]));

	int ilat = std::floor(lat);
	int ilon = std::floor(lon);

	auto id = std::make_pair(ilat, ilon);

	if (no_datums.count(id) > 0)
		return -10.f;

	if (datums.count(id) == 0)
	{
		if (datums.size() > 40)
			datums.clear();

		std::ostringstream os;
		if (ilat >= 0)
			os << 'N' << std::setw(2) << std::setfill('0') << ilat;
		else
			os << 'S' << std::setw(2) << std::setfill('0') << (-ilat);
		if (ilon >= 0)
			os << 'E' << std::setw(3) << std::setfill('0') << ilon;
		else
			os << 'W' << std::setw(3) << std::setfill('0') << (-ilon);

		os << ".hgt";

		std::string const filename = data_path + os.str();
		std::ifstream ifs(filename, std::ios::binary);

		if (!ifs)
		{
			no_datums.insert(id);
			return -10.f;
		}

		return 10.f;

		std::unique_ptr<std::int16_t[]> values{new std::int16_t[3601 * 3601]};
		ifs.read(reinterpret_cast<char *>(values.get()), 3601 * 3601 * 2);
		datums[id] = std::move(values);
	}

	auto const * values = datums[id].get();

	int const tlat = geom::clamp<int>(std::round((lat - ilat) * 3600), {0, 3600});
	int const tlon = geom::clamp<int>(std::round((lon - ilon) * 3600), {0, 3600});

	auto h = values[tlat * 3601 + 3600 - tlon];
	h = ((h & 255) << 8) | (h >> 8);
	height_range |= (int)h;
	return h;
}

struct node
{
	geom::vector<float, 3> v[3];

	virtual void draw(int level) = 0;
	virtual node * child(int id) = 0;
};

struct node_controller
{
	node_controller();

	node * root(int f);

	std::size_t node_count() const { return node_count_; }

private:
	struct node_impl
		: node
	{
		node_controller * controller;
		gfx::array array;
		gfx::buffer height_buffer;
		std::unique_ptr<node> children[256];

		std::vector<std::int16_t> height_data;
		std::atomic<bool> height_data_ready = false;
		bool height_data_loaded = false;

		void draw(int level) override;
		node * child(int id) override;

		void load_heights();
	};

	cg::icosahedron<float> icosahedron_;

	gfx::buffer index_buffer_;
	std::size_t index_counts_[10];
	std::unique_ptr<node_impl> roots_[20];

	height_provider height_provider_;

	util::threadpool loader_{"load", 1};

	std::size_t node_count_ = 0;

	std::unique_ptr<node_impl> make_node();
};

node_controller::node_controller()
	: icosahedron_{geom::point<float, 3>::zero(), 1.f}
{
	std::vector<std::uint32_t> indices;
	index_counts_[0] = 0;
	for (std::size_t N = 0; N <= 8; ++N)
	{
		std::size_t step = 256 >> N;

		auto idx = [step](std::size_t i, std::size_t j) -> std::uint32_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(0xffffffffu);
		}

		index_counts_[N + 1] = indices.size();
	}
	index_buffer_.load(indices, gl::STATIC_DRAW);
}

node * node_controller::root(int f)
{
	if (!roots_[f])
	{
		auto n = make_node();

		auto face = cg::faces(icosahedron_)[f];
		n->v[0] = icosahedron_.vertices[face[0]] - geom::point<float, 3>::zero();
		n->v[1] = icosahedron_.vertices[face[1]] - geom::point<float, 3>::zero();
		n->v[2] = icosahedron_.vertices[face[2]] - geom::point<float, 3>::zero();

		n->load_heights();

		roots_[f] = std::move(n);
	}

	return roots_[f].get();
}

std::unique_ptr<node_controller::node_impl> node_controller::make_node()
{
	auto n = std::make_unique<node_controller::node_impl>();

	n->controller = this;
	n->array.bind();
	n->height_buffer.bind();
	gl::EnableVertexAttribArray(0);
	gl::VertexAttribPointer(0, 1, gl::SHORT, gl::FALSE, 0, nullptr);
	gl::BindBuffer(gl::ELEMENT_ARRAY_BUFFER, index_buffer_.id());

	++node_count_;

	return n;
}

void node_controller::node_impl::draw(int level)
{
	if (!height_data_loaded)
	{
		if (!height_data_ready) return;

		height_buffer.load(height_data, gl::STATIC_DRAW);
		height_data.clear();
		height_data_loaded = true;
	}

	array.bind();
	std::size_t offset = controller->index_counts_[level];
	std::size_t count = controller->index_counts_[level + 1] - offset;
	gl::DrawElements(gl::TRIANGLE_STRIP, count, gl::UNSIGNED_INT, (std::uint32_t const *)(nullptr) + offset);
}

node * node_controller::node_impl::child(int id)
{
	if (!children[id])
	{
		auto n = controller->make_node();

		int i0, j0, i1, j1, i2, j2;

		if (id < 136)
		{
			int i = int(std::floor(0.5f * (sqrt(1.f + 8.f * id) - 1.f)));
			int j = id - (i * (i + 1)) / 2;

			i0 = i;
			j0 = j;
			i1 = i + 1;
			j1 = j + 1;
			i2 = i + 1;
			j2 = j;
		}
		else
		{
			int i = int(std::floor(0.5f * (sqrt(1.f + 8.f * (id - 136)) - 1.f)));
			int j = id - 136 - (i * (i + 1)) / 2;

			i0 = i + 1;
			j0 = j;
			i1 = i + 1;
			j1 = j + 1;
			i2 = i + 2;
			j2 = j + 1;
		}

		auto at = [this](int i, int j)
		{
			float t0 = 1.f - (1.f * i) / 16.f;
			float t1 = (1.f * j) / 16.f;
			float t2 = 1.f - t0 - t1;

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

		n->v[0] = at(i0, j0);
		n->v[1] = at(i1, j1);
		n->v[2] = at(i2, j2);

		n->load_heights();

		children[id] = std::move(n);
	}

	return children[id].get();
}

void node_controller::node_impl::load_heights()
{
	controller->loader_.dispatch([this]{
		height_data.assign((257 * 258) / 2, 0);

		auto * out = height_data.data();

		auto at = [this](int i, int j)
		{
			float t0 = 1.f - (1.f * i) / 256.f;
			float t1 = (1.f * j) / 256.f;
			float t2 = 1.f - t0 - t1;

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

		for (int i = 0; i <= 256; ++i)
		{
			for (int j = 0; j <= i; ++j)
			{
				*out++ = static_cast<std::int16_t>(controller->height_provider_.height_at(at(i, j)));
			}
		}

		height_data_ready = true;
	});
}

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

uniform mat4 u_transform;

uniform int u_N;
uniform vec3 u_p0;
uniform vec3 u_p1;
uniform vec3 u_p2;

layout (location = 0) in float in_height;

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) * (1.0 + in_height / 6400000.0);
	gl_Position = u_transform * vec4(p, 1.0);
	color = (in_height < 0.0) ? vec3(0.0, 0.0, 0.5) : vec3(0.0, 0.5, 0.0);
})";

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

in vec3 color;
out vec4 out_color;

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

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::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;

	node_controller nodes;

	gfx::program tile_program{tile_vs, tile_fs};

	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(true);
	show_cursor(false);

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

	camera_azimuthal_angle_updater = camera.azimuthal_angle;
	camera_elevation_angle_updater = camera.elevation_angle;
}

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()
{
	cg::icosahedron<float> icosahedron{geom::point<float, 3>::zero(), 1.f};

	auto const & icosa_vertices = cg::vertices(icosahedron);
	auto const & icosa_faces = cg::faces(icosahedron);
	auto const icosa_side = geom::distance(icosa_vertices[icosa_faces[0][0]], icosa_vertices[icosa_faces[0][1]]);

	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::FILL);
	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(0xffffffffu);

	{
		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);
	(void)frustum;

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

	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, node * n, int level = 0) -> void
	{
		auto const & v = n->v;
		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 = false && 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 = 5.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 > 8)
		{
			for (int id = 0; id < 256; ++id)
				self(n->child(id), level + 1);
		}
		else
		{
			tile_n = geom::clamp(tile_n, {0, 8});

			++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];

			n->draw(tile_n);
		}
	});

	gfx::check_error();

	for (int f = 0; f < 20; ++f)
	{
		visit(nodes.root(f));
	}

	info.push_back(util::to_string("Tiles: ", rendered_tiles));
	{
		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.push_back(util::to_string("Heights: ", height_range));
	info.push_back(util::to_string("Nodes: ", nodes.node_count()));

	{
		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>();
}