From ba9c7a4418fa71d734a07c143a6a64db97e36559 Mon Sep 17 00:00:00 2001 From: lisyarus Date: Sun, 8 Nov 2020 16:19:25 +0300 Subject: [PATCH] SRTM example (wip) --- examples/srtm.cpp | 305 +++++++++++++++++++++++++++++++++------------- 1 file changed, 221 insertions(+), 84 deletions(-) diff --git a/examples/srtm.cpp b/examples/srtm.cpp index d5cec440..6351d05b 100644 --- a/examples/srtm.cpp +++ b/examples/srtm.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -27,6 +28,7 @@ #include #include #include +#include #include #include @@ -79,20 +81,22 @@ struct height_provider } }; - std::unordered_map, std::unique_ptr, datum_id_hash> datums; + 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; }; -static geom::interval height_range; - float height_provider::height_at(geom::vector const & v) { - static std::string const data_path = "/home/lisyarus/data/srtm/dem-unpacked/"; +// 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; +// if (std::abs(lat) > 60.f) return 0.f; +// return -10.f; float const lon = geom::deg(std::atan2(v[0], -v[1])); @@ -101,15 +105,30 @@ float height_provider::height_at(geom::vector const & v) 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::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 @@ -119,40 +138,63 @@ float height_provider::height_at(geom::vector const & v) else os << 'W' << std::setw(3) << std::setfill('0') << (-ilon); - os << ".hgt"; + 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 -10.f; + return 0.f; } - return 10.f; +// return 10.f; - std::unique_ptr values{new std::int16_t[3601 * 3601]}; - ifs.read(reinterpret_cast(values.get()), 3601 * 3601 * 2); - datums[id] = std::move(values); + 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 const * values = datums[id].get(); + auto at = [values](int lat, int lon) -> float + { + return (int)values[(3600 - lat) * 3601 + lon] - 32768; + }; - int const tlat = geom::clamp(std::round((lat - ilat) * 3600), {0, 3600}); - int const tlon = geom::clamp(std::round((lon - ilon) * 3600), {0, 3600}); + 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}); - auto h = values[tlat * 3601 + 3600 - tlon]; - h = ((h & 255) << 8) | (h >> 8); - height_range |= (int)h; - return h; + 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 void draw(int level) = 0; + virtual bool draw(int level) = 0; virtual node * child(int id) = 0; }; @@ -164,6 +206,8 @@ struct node_controller 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 @@ -171,13 +215,14 @@ private: node_controller * controller; gfx::array array; gfx::buffer height_buffer; - std::unique_ptr children[256]; + 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; - void draw(int level) override; + bool draw(int level) override; node * child(int id) override; void load_heights(); @@ -186,12 +231,12 @@ private: cg::icosahedron icosahedron_; gfx::buffer index_buffer_; - std::size_t index_counts_[10]; + std::size_t index_counts_[node_size_log2 + 2]; std::unique_ptr roots_[20]; height_provider height_provider_; - util::threadpool loader_{"load", 1}; + util::threadpool loader_{"load", 4}; std::size_t node_count_ = 0; @@ -201,13 +246,13 @@ private: node_controller::node_controller() : icosahedron_{geom::point::zero(), 1.f} { - std::vector indices; + std::vector indices; index_counts_[0] = 0; - for (std::size_t N = 0; N <= 8; ++N) + 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::uint32_t + auto idx = [step](std::size_t i, std::size_t j) -> std::uint16_t { i *= step; j *= step; @@ -222,7 +267,7 @@ node_controller::node_controller() indices.push_back(idx(i, j)); } indices.push_back(idx(i + 1, i + 1)); - indices.push_back(0xffffffffu); + indices.push_back(0xffffu); } index_counts_[N + 1] = indices.size(); @@ -241,8 +286,6 @@ node * node_controller::root(int f) n->v[1] = icosahedron_.vertices[face[1]] - geom::point::zero(); n->v[2] = icosahedron_.vertices[face[2]] - geom::point::zero(); - n->load_heights(); - roots_[f] = std::move(n); } @@ -265,11 +308,20 @@ std::unique_ptr node_controller::make_node() return n; } -void node_controller::node_impl::draw(int level) +bool node_controller::node_impl::draw(int level) { if (!height_data_loaded) { - if (!height_data_ready) return; + 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(); @@ -279,7 +331,9 @@ void node_controller::node_impl::draw(int level) 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); + 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) @@ -290,7 +344,9 @@ node * node_controller::node_impl::child(int id) int i0, j0, i1, j1, i2, j2; - if (id < 136) + 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; @@ -304,8 +360,8 @@ node * node_controller::node_impl::child(int id) } 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; + 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; @@ -317,8 +373,8 @@ node * node_controller::node_impl::child(int id) auto at = [this](int i, int j) { - float t0 = 1.f - (1.f * i) / 16.f; - float t1 = (1.f * j) / 16.f; + 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); @@ -328,8 +384,6 @@ node * node_controller::node_impl::child(int id) n->v[1] = at(i1, j1); n->v[2] = at(i2, j2); - n->load_heights(); - children[id] = std::move(n); } @@ -339,20 +393,20 @@ node * node_controller::node_impl::child(int id) void node_controller::node_impl::load_heights() { controller->loader_.dispatch([this]{ - height_data.assign((257 * 258) / 2, 0); + 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) / 256.f; - float t1 = (1.f * j) / 256.f; + 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 <= 256; ++i) + for (int i = 0; i <= node_size; ++i) { for (int j = 0; j <= i; ++j) { @@ -374,9 +428,13 @@ 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() { @@ -388,19 +446,27 @@ void main() 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) ? vec3(0.0, 0.0, 0.5) : vec3(0.0, 0.5, 0.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() { - out_color = vec4(color, 1.0); + 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 @@ -416,17 +482,21 @@ struct srtm_app void present() override; geom::free_camera camera; - smooth_updater camera_azimuthal_angle_updater{camera.azimuthal_angle, 20.f}; - smooth_updater camera_elevation_angle_updater{camera.elevation_angle, 20.f}; + 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; }; @@ -440,25 +510,52 @@ srtm_app::srtm_app() 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.rotateYZ(geom::rad(-90.f)); - camera_azimuthal_angle_updater = camera.azimuthal_angle; - camera_elevation_angle_updater = camera.elevation_angle; + 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_azimuthal_angle_updater = camera_azimuthal_angle_updater - 0.01f * dx; - camera_elevation_angle_updater = camera_elevation_angle_updater + 0.01f * dy; + camera.rotateZX(0.01f * dx); + camera.rotateYZ(0.01f * dy); } void srtm_app::update() @@ -468,15 +565,14 @@ void srtm_app::update() if (is_key_down(SDLK_q)) { + camera.rotateXY(-0.04f); } if (is_key_down(SDLK_e)) { + camera.rotateXY(0.04f); } - 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::zero()) - 1.f); auto const camera_forward = camera.direction(); @@ -512,6 +608,8 @@ void srtm_app::update() { camera.pos -= camera_speed * dt * camera_up; } + + camera_transform += (camera.transform() - camera_transform) * std::min(1.f, 10.f * dt); } namespace std @@ -556,7 +654,7 @@ void srtm_app::present() gl::DepthFunc(gl::LEQUAL); gl::Enable(gl::PRIMITIVE_RESTART); - gl::PrimitiveRestartIndex(0xffffffffu); + gl::PrimitiveRestartIndex(0xffffu); { auto d = geom::distance(camera.pos, geom::point::zero()); @@ -564,9 +662,9 @@ void srtm_app::present() camera.near_clip = (d > 2.f) ? d - 2.f : 0.0001f; } - auto const camera_transform = camera.transform(); +// auto const camera_transform = camera.transform(); auto const camera_pos = camera.position(); -// auto const camera_direction = camera.direction(); + 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")); @@ -575,18 +673,26 @@ void srtm_app::present() tile_program.bind(); tile_program["u_transform"] = camera_transform; - tile_program["u_N"] = static_cast(256); + 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) -> void + auto visit = util::recursive([&](auto & self, node * n, int level = 0) -> bool { auto const & v = n->v; auto const o = geom::point::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; @@ -599,17 +705,24 @@ void srtm_app::present() } } if (culled) - return; + return true; (void)culled; } - geom::triangle> t{o + v[0], o + v[1], o + v[2]}; - cg::triangular_prism body{t, m}; + bool const flat = n->height_range.max == n->height_range.min; - if (cg::separation(body, frustum).second > 0.f) - return; + 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); -// bool const selected = false && geom::intersect(geom::ray{camera_pos, camera_direction}, t); + 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; @@ -644,18 +757,34 @@ void srtm_app::present() float const max_triangle_size = 5.f; // pixels - float const side_length = icosa_side / (1 << (level * 4)); + 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)); - if (level < 3 && tile_n > 8) + bool should_draw_children = !flat && level < max_child_level && tile_n > node_size_log2; + bool all_children_drawn = true; + + if (selected && !should_draw_children) { - for (int id = 0; id < 256; ++id) - self(n->child(id), level + 1); + 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]]); + } } - else + + if (should_draw_children) { - tile_n = geom::clamp(tile_n, {0, 8}); + 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; @@ -672,8 +801,10 @@ void srtm_app::present() tile_program["u_p2"] = v[2]; tile_program["u_color"] = colors[tile_n % 4]; - n->draw(tile_n); + return n->draw(tile_n); } + + return true; }); gfx::check_error(); @@ -683,14 +814,20 @@ void srtm_app::present() 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("Heights: ", height_range)); 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()));