#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace psemek; template 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 const & v); struct datum_id_hash { std::size_t operator() (std::pair const & p) const { return ((p.first + 180) << 16) | (p.second + 90); } }; std::mutex datums_mutex; std::unordered_map, std::unique_ptr, datum_id_hash> datums; std::mutex no_datums_mutex; std::unordered_set, datum_id_hash> no_datums; }; float height_provider::height_at(geom::vector const & v) { // static std::string const data_path = "/home/lisyarus/data/srtm/dem/"; static std::string const data_path = "/home/lisyarus/data/srtm/test/"; 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); { std::lock_guard lock{no_datums_mutex}; if (no_datums.count(id) > 0) return 0.f; } std::uint16_t const * values = nullptr; { std::lock_guard lock{datums_mutex}; auto it = datums.find(id); if (it != datums.end()) values = it->second.get(); } if (!values) { { std::lock_guard lock{datums_mutex}; if (datums.size() > 40) datums.clear(); } std::ostringstream os; // os << "ASTGTMV003_"; 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 << ".gray"; // os << ".zip"; std::string const filename = data_path + os.str(); std::ifstream ifs(filename, std::ios::binary); if (!ifs) { std::lock_guard lock{no_datums_mutex}; no_datums.insert(id); return 0.f; } // return 10.f; std::unique_ptr data{new std::uint16_t[3601 * 3601]}; ifs.read(reinterpret_cast(data.get()), 3601 * 3601 * 2); { std::lock_guard lock{datums_mutex}; auto res = datums.insert(std::make_pair(id, std::move(data))); values = res.first->second.get(); } } auto at = [values](int lat, int lon) -> float { return (int)values[(3600 - lat) * 3601 + lon] - 32768; }; int const tlat = geom::clamp(std::floor((lat - ilat) * 3600), {0, 3600 - 1}); int const tlon = geom::clamp(std::floor((lon - ilon) * 3600), {0, 3600 - 1}); float const mlat = (lat - ilat) * 3600.f - tlat; float const mlon = (lon - ilon) * 3600.f - tlon; float h00 = at(tlat, tlon); float h01 = at(tlat + 1, tlon); float h10 = at(tlat, tlon + 1); float h11 = at(tlat + 1, tlon + 1); return geom::lerp(geom::lerp(h00, h01, mlat), geom::lerp(h10, h11, mlat), mlon); } static constexpr int node_size_log2 = 8; static constexpr int node_size = 1 << node_size_log2; static constexpr int node_child_depth = 1; static constexpr int node_child_size = 1 << node_child_depth; static constexpr int node_child_count = node_child_size * node_child_size; static constexpr int max_child_level = 20 - node_size_log2 - node_child_depth; struct node { geom::vector v[3]; geom::interval height_range; virtual bool 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_; } std::size_t loader_queue_size() const { return loader_.queue_size(); } private: struct node_impl : node { node_controller * controller; gfx::array array; gfx::buffer height_buffer; std::unique_ptr children[node_child_count]; std::vector height_data; std::atomic height_data_ready = false; bool height_data_loaded = false; bool height_data_loader_dispatched = false; bool draw(int level) override; node * child(int id) override; void load_heights(); }; cg::icosahedron icosahedron_; gfx::buffer index_buffer_; std::size_t index_counts_[node_size_log2 + 2]; std::unique_ptr roots_[20]; height_provider height_provider_; util::threadpool loader_{"load", 4}; std::size_t node_count_ = 0; std::unique_ptr make_node(); }; node_controller::node_controller() : icosahedron_{geom::point::zero(), 1.f} { std::vector indices; index_counts_[0] = 0; for (std::size_t N = 0; N <= node_size_log2; ++N) { std::size_t step = 256 >> 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); } 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::zero(); n->v[1] = icosahedron_.vertices[face[1]] - geom::point::zero(); n->v[2] = icosahedron_.vertices[face[2]] - geom::point::zero(); roots_[f] = std::move(n); } return roots_[f].get(); } std::unique_ptr node_controller::make_node() { auto n = std::make_unique(); 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; } bool node_controller::node_impl::draw(int level) { if (!height_data_loaded) { if (!height_data_loader_dispatched) { height_data_loader_dispatched = true; load_heights(); } if (!height_data_ready) return false; for (auto h : height_data) height_range |= static_cast(h); 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_SHORT, (std::uint16_t const *)(nullptr) + offset); return true; } node * node_controller::node_impl::child(int id) { if (!children[id]) { auto n = controller->make_node(); int i0, j0, i1, j1, i2, j2; static constexpr int type_1_count = (node_child_size * (node_child_size + 1)) / 2; if (id < type_1_count) { 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 - type_1_count)) - 1.f))); int j = id - type_1_count - (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) / (1.f * node_child_size); float t1 = (1.f * j) / (1.f * node_child_size); 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); children[id] = std::move(n); } return children[id].get(); } void node_controller::node_impl::load_heights() { controller->loader_.dispatch([this]{ height_data.assign(((node_size + 2) * (node_size + 1)) / 2, 0); auto * out = height_data.data(); auto at = [this](int i, int j) { float t0 = 1.f - (1.f * i) / (1.f * node_size); float t1 = (1.f * j) / (1.f * node_size); float t2 = 1.f - t0 - t1; return geom::normalized(v[0] * t0 + v[1] * t1 + v[2] * t2); }; for (int i = 0; i <= node_size; ++i) { for (int j = 0; j <= i; ++j) { *out++ = static_cast(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; uniform sampler1D u_colormap; uniform sampler1D u_colormap_neg; layout (location = 0) in float in_height; out vec3 color; out vec3 pos; 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); pos = p; gl_Position = u_transform * vec4(p, 1.0); color = (in_height > 0.0) ? texture(u_colormap, in_height / 8000.0).rgb : texture(u_colormap_neg, -in_height / 10000.0).rgb; })"; static char const tile_fs[] = R"(#version 330 uniform vec3 u_light; in vec3 color; in vec3 pos; out vec4 out_color; void main() { vec3 dx = dFdx(pos); vec3 dy = dFdy(pos); vec3 n = normalize(cross(dx, dy)); float l = (0.5 + dot(n, u_light) * 0.5); out_color = vec4(color * l, 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; geom::matrix camera_transform; bool camera_forward = false; node_controller nodes; gfx::program tile_program{tile_vs, tile_fs}; gfx::texture_1d color_map; gfx::texture_1d color_map_neg; util::clock> frame_clock; util::moving_average frame_dt_average{32}; gfx::mesh selected_mesh; gfx::simple_renderer simple_renderer; 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.rotateYZ(geom::rad(-90.f)); camera_transform = camera.transform(); selected_mesh.setup>(); { util::array colors({5}); colors(0) = {0, 127, 0}; colors(1) = {127, 127, 0}; colors(2) = {127, 63, 0}; colors(3) = {127, 0, 0}; colors(4) = {191, 191, 191}; color_map.load(colors); color_map.clamp(); color_map.linear_filter(); } { util::array colors({5}); colors(0) = {0, 63, 127}; colors(1) = {0, 0, 127}; colors(2) = {0, 0, 127}; colors(3) = {0, 0, 127}; colors(4) = {0, 127, 127}; color_map_neg.load(colors); color_map_neg.clamp(); color_map_neg.linear_filter(); } } void srtm_app::on_resize(int width, int height) { app::on_resize(width, height); camera.set_fov(camera.fov_y, (1.f * width) / height); camera_transform = camera.transform(); } void srtm_app::on_mouse_move(int x, int y, int dx, int dy) { app::on_mouse_move(x, y, dx, dy); camera.rotateZX(0.01f * dx); camera.rotateYZ(0.01f * dy); } void srtm_app::update() { float dt = frame_clock.restart().count(); frame_dt_average.push(dt); if (is_key_down(SDLK_q)) { camera.rotateXY(-0.04f); } if (is_key_down(SDLK_e)) { camera.rotateXY(0.04f); } float const camera_speed = std::min(5.f, geom::distance(camera.pos, geom::point::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; } camera_transform += (camera.transform() - camera_transform) * std::min(1.f, 10.f * dt); } namespace std { template std::ostream & operator << (std::ostream & os, std::vector 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 icosahedron{geom::point::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 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(0xffffu); { auto d = geom::distance(camera.pos, geom::point::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::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(node_size); tile_program["u_light"] = geom::vector{0.f, 0.f, 1.f}; tile_program["u_colormap"] = 0; tile_program["u_colormap_neg"] = 1; gl::ActiveTexture(gl::TEXTURE0); color_map.bind(); gl::ActiveTexture(gl::TEXTURE1); color_map_neg.bind(); gl::ActiveTexture(gl::TEXTURE0); info.push_back(util::to_string("Camera pos: ", camera_pos)); std::vector> selected_vertices; std::size_t rendered_tiles = 0; std::vector id; auto visit = util::recursive([&](auto & self, node * n, int level = 0) -> bool { auto const & v = n->v; auto const o = geom::point::zero(); { 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 true; (void)culled; } bool const flat = n->height_range.max == n->height_range.min; auto const m3 = (v[0] + v[1] + v[2]) / 3.f; auto const m = geom::normalized(m3); auto const m0 = m * (n->height_range.empty() ? 0.f : n->height_range.min) / 6400000.f; auto const m1 = m * (n->height_range.empty() ? 1.f : flat ? n->height_range.min + 1.f : n->height_range.max) / 6400000.f + (m - m3); geom::triangle> t{o + v[0] + m0, o + v[1] + m0, o + v[2] + m0}; cg::triangular_prism body{t, m1 - m0}; // if (cg::separation(body, frustum).second > 0.f) // return true; 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 { 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::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 * node_child_depth)); int tile_n = std::ceil(std::log2(on_screen_unit * side_length / max_triangle_size)); bool should_draw_children = !flat && level < max_child_level && tile_n > node_size_log2; bool all_children_drawn = true; if (selected && !should_draw_children) { auto const & vs = cg::vertices(body); auto const & es = cg::edges(body); for (auto const & e : es) { selected_vertices.push_back(vs[e[0]]); selected_vertices.push_back(vs[e[1]]); } } if (should_draw_children) { for (int id = 0; id < node_child_count; ++id) all_children_drawn &= self(n->child(id), level + 1); } if (!should_draw_children || !all_children_drawn) { tile_n = geom::clamp(tile_n, {0, node_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]; return n->draw(tile_n); } return true; }); gfx::check_error(); for (int f = 0; f < 20; ++f) { visit(nodes.root(f)); } // selected_mesh.load(selected_vertices, gl::LINES, gl::STREAM_DRAW); // simple_renderer.push(gfx::simple_renderer::render_state{&selected_mesh, gfx::white.as_color_rgba()}); // simple_renderer.render(gfx::simple_renderer::render_options{camera_transform}); 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("Nodes: ", nodes.node_count())); info.push_back(util::to_string("Tasks: ", nodes.loader_queue_size())); info.push_back(util::to_string("Selected: ", selected_mesh.index_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(); }