From c41195277cf2ba3d26c49f04b776821c6426fc54 Mon Sep 17 00:00:00 2001 From: program-- Date: Mon, 10 Jul 2023 08:52:44 -0700 Subject: [PATCH] rev: fix off-by-one in WKB reader; reorg WKB impl Review changes per https://github.com/NOAA-OWP/ngen/pull/522#discussion_r1256246141. This commit also reorganizes the WKB implementation by moving inline definitions to a source file, since there is no need to inline those functions. --- include/geopackage/WKB.hpp | 283 ++----------------------------- src/geopackage/CMakeLists.txt | 1 + src/geopackage/wkb.cpp | 306 ++++++++++++++++++++++++++++++++++ test/geopackage/WKB_Test.cpp | 64 +++++++ 4 files changed, 383 insertions(+), 271 deletions(-) create mode 100644 src/geopackage/wkb.cpp diff --git a/include/geopackage/WKB.hpp b/include/geopackage/WKB.hpp index 7e0bdc4b8a..d7386324ba 100644 --- a/include/geopackage/WKB.hpp +++ b/include/geopackage/WKB.hpp @@ -88,148 +88,6 @@ struct wkb { static multipolygon_t read_multipolygon(const byte_vector&, int&, uint8_t); }; -inline typename wkb::point_t wkb::read_point(const byte_vector& buffer, int& index, uint8_t order) -{ - double x, y; - utils::copy_from(buffer, index, x, order); - utils::copy_from(buffer, index, y, order); - return point_t{x, y}; -} - -inline typename wkb::linestring_t wkb::read_linestring(const byte_vector& buffer, int& index, uint8_t order) -{ - uint32_t count; - utils::copy_from(buffer, index, count, order); - - linestring_t linestring; - linestring.resize(count); - for (auto& child : linestring) { - child = read_point(buffer, index, order); - } - - return linestring; -} - -inline typename wkb::polygon_t wkb::read_polygon(const byte_vector& buffer, int& index, uint8_t order) -{ - uint32_t count; - utils::copy_from(buffer, index, count, order); - - polygon_t polygon; - - if (count > 1) { - // polygons only have 1 outer ring, - // so any extra vectors are considered to be - // inner rings. - polygon.inners().resize(count); - } - - auto outer = read_linestring(buffer, index, order); - polygon.outer().reserve(outer.size()); - for (auto& p : outer) { - polygon.outer().push_back(p); - } - - for (uint32_t i = 1; i < count; i++) { - auto inner = read_linestring(buffer, index, order); - polygon.inners().at(i).reserve(inner.size()); - for (auto& p : inner) { - polygon.inners().at(i).push_back(p); - } - } - - return polygon; -} - -inline typename wkb::multipoint_t wkb::read_multipoint(const byte_vector& buffer, int& index, uint8_t order) -{ - uint32_t count; - utils::copy_from(buffer, index, count, order); - - multipoint_t mp; - mp.resize(count); - - for (auto& point : mp) { - const byte_t new_order = buffer[index]; - index++; - - uint32_t type; - utils::copy_from(buffer, index, type, new_order); - - point = read_point(buffer, index, new_order); - } - - return mp; -} - -inline typename wkb::multilinestring_t wkb::read_multilinestring(const byte_vector& buffer, int& index, uint8_t order) -{ - uint32_t count; - utils::copy_from(buffer, index, count, order); - - multilinestring_t ml; - ml.resize(count); - for (auto& line : ml) { - const byte_t new_order = buffer[index]; - index++; - - uint32_t type; - utils::copy_from(buffer, index, type, new_order); - - line = read_linestring(buffer, index, new_order); - } - - return ml; -} - -inline typename wkb::multipolygon_t wkb::read_multipolygon(const byte_vector& buffer, int& index, uint8_t order) -{ - uint32_t count; - utils::copy_from(buffer, index, count, order); - - multipolygon_t mpl; - mpl.resize(count); - for (auto& polygon : mpl) { - const byte_t new_order = buffer[index]; - index++; - - uint32_t type; - utils::copy_from(buffer, index, type, new_order); - - polygon = read_polygon(buffer, index, new_order); - } - - return mpl; -} - - -inline typename wkb::geometry wkb::read(const byte_vector& buffer) -{ - if (buffer.size() < 5) { - throw std::runtime_error("buffer reached end before encountering WKB"); - } - - int index = 0; - const byte_t order = buffer[index]; - index++; - - uint32_t type; - utils::copy_from(buffer, index, type, order); - - geometry g; - switch(type) { - case 1: g = read_point(buffer, index, order); break; - case 2: g = read_linestring(buffer, index, order); break; - case 3: g = read_polygon(buffer, index, order); break; - case 4: g = read_multipoint(buffer, index, order); break; - case 5: g = read_multilinestring(buffer, index, order); break; - case 6: g = read_multipolygon(buffer, index, order); break; - default: g = point_t{std::nan("0"), std::nan("0")}; break; - } - - return g; -} - /** * EPSG 5070 projection definition for use with boost::geometry. * @@ -246,6 +104,12 @@ const auto epsg5070 = bg::srs::dpar::parameters<>(bg::srs::dpar::proj_aea) (bg::srs::dpar::x_0, 0) (bg::srs::dpar::y_0, 0); +/** + * EPSG 3857 projection definition for use with boost::geometry. + * + * @note this is required because boost 1.72.0 does not + * have an EPSG definition for 3857 in boost::srs::epsg. + */ const auto epsg3857 = bg::srs::dpar::parameters<>(bg::srs::dpar::proj_merc) (bg::srs::dpar::units_m) (bg::srs::dpar::no_defs) @@ -257,141 +121,18 @@ const auto epsg3857 = bg::srs::dpar::parameters<>(bg::srs::dpar::proj_merc) (bg::srs::dpar::y_0, 0) (bg::srs::dpar::k, 1); -inline bg::srs::dpar::parameters<> wkb::get_prj(uint32_t srid) { - switch(srid) { - case 5070: - return epsg5070; - case 3857: - return epsg3857; - default: - return bg::projections::detail::epsg_to_parameters(srid); - } -} - struct wkb::wgs84 : public boost::static_visitor { wgs84(uint32_t srs, const bg::srs::transformation<>& tr) : srs(srs) , tr(tr) {}; - geojson::geometry operator()(point_t& g) - { - if (this->srs == 4326) { - return geojson::coordinate_t(g.get<0>(), g.get<1>()); - } - - geojson::coordinate_t h; - this->tr.forward(g, h); - return h; - } - - geojson::geometry operator()(linestring_t& g) - { - geojson::linestring_t h; - - if (this->srs == 4326) { - h.reserve(g.size()); - for (auto&& gg : g) { - h.emplace_back( - std::move(gg.get<0>()), - std::move(gg.get<1>()) - ); - } - } else { - this->tr.forward(g, h); - } - return h; - } - - geojson::geometry operator()(polygon_t& g) - { - geojson::polygon_t h; - - if(this->srs == 4326) { - h.outer().reserve(g.outer().size()); - for (auto&& gg : g.outer()) { - h.outer().emplace_back( - std::move(gg.get<0>()), - std::move(gg.get<1>()) - ); - } - - h.inners().resize(g.inners().size()); - auto&& inner_g = g.inners().begin(); - auto&& inner_h = h.inners().begin(); - for (; inner_g != g.inners().end(); inner_g++, inner_h++) { - inner_h->reserve(inner_g->size()); - for (auto&& gg : *inner_g) { - inner_h->emplace_back( - std::move(gg.get<0>()), - std::move(gg.get<1>()) - ); - } - } - } else { - this->tr.forward(g, h); - } - return h; - } - - geojson::geometry operator()(multipoint_t& g) - { - geojson::multipoint_t h; - - if (this->srs == 4326) { - h.reserve(g.size()); - for (auto&& gg : g) { - h.emplace_back( - std::move(gg.get<0>()), - std::move(gg.get<1>()) - ); - } - } else { - this->tr.forward(g, h); - } - - return h; - } - - geojson::geometry operator()(multilinestring_t& g) - { - geojson::multilinestring_t h; - - if (this->srs == 4326) { - h.resize(g.size()); - auto&& line_g = g.begin(); - auto&& line_h = h.begin(); - for(; line_g != g.end(); line_g++, line_h++) { - *line_h = std::move( - boost::get(this->operator()(*line_g)) - ); - } - } else { - this->tr.forward(g, h); - } - - return h; - } - - geojson::geometry operator()(multipolygon_t& g) - { - geojson::multipolygon_t h; - - if (this->srs == 4326) { - h.resize(g.size()); - auto&& polygon_g = g.begin(); - auto&& polygon_h = h.begin(); - for (; polygon_g != g.end(); polygon_g++, polygon_h++) { - *polygon_h = std::move( - boost::get(this->operator()(*polygon_g)) - ); - } - } else { - this->tr.forward(g, h); - } - - return h; - } + geojson::geometry operator()(point_t& g); + geojson::geometry operator()(linestring_t& g); + geojson::geometry operator()(polygon_t& g); + geojson::geometry operator()(multipoint_t& g); + geojson::geometry operator()(multilinestring_t& g); + geojson::geometry operator()(multipolygon_t& g); private: uint32_t srs; diff --git a/src/geopackage/CMakeLists.txt b/src/geopackage/CMakeLists.txt index 2aead82963..5a61c54658 100644 --- a/src/geopackage/CMakeLists.txt +++ b/src/geopackage/CMakeLists.txt @@ -15,6 +15,7 @@ add_library(geopackage STATIC geometry.cpp properties.cpp feature.cpp read.cpp + wkb.cpp sqlite/iterator.cpp sqlite/database.cpp ) diff --git a/src/geopackage/wkb.cpp b/src/geopackage/wkb.cpp new file mode 100644 index 0000000000..c2cabc1d47 --- /dev/null +++ b/src/geopackage/wkb.cpp @@ -0,0 +1,306 @@ +#include "WKB.hpp" + +using namespace geopackage; + +// ---------------------------------------------------------------------------- +// WKB Readers +// ---------------------------------------------------------------------------- + +typename wkb::point_t wkb::read_point(const byte_vector& buffer, int& index, uint8_t order) +{ + double x, y; + utils::copy_from(buffer, index, x, order); + utils::copy_from(buffer, index, y, order); + return point_t{x, y}; +} + +// ---------------------------------------------------------------------------- + +typename wkb::linestring_t wkb::read_linestring(const byte_vector& buffer, int& index, uint8_t order) +{ + uint32_t count; + utils::copy_from(buffer, index, count, order); + + linestring_t linestring; + linestring.resize(count); + for (auto& child : linestring) { + child = read_point(buffer, index, order); + } + + return linestring; +} + +// ---------------------------------------------------------------------------- + +typename wkb::polygon_t wkb::read_polygon(const byte_vector& buffer, int& index, uint8_t order) +{ + uint32_t count; + utils::copy_from(buffer, index, count, order); + + polygon_t polygon; + + if (count > 1) { + // polygons only have 1 outer ring, + // so any extra vectors are considered to be + // inner rings. + polygon.inners().resize(count - 1); + } + + auto outer = read_linestring(buffer, index, order); + polygon.outer().reserve(outer.size()); + for (auto& p : outer) { + polygon.outer().push_back(p); + } + + for (uint32_t i = 0; i < count - 1; i++) { + auto inner = read_linestring(buffer, index, order); + polygon.inners().at(i).reserve(inner.size()); + for (auto& p : inner) { + polygon.inners().at(i).push_back(p); + } + } + + return polygon; +} + +// ---------------------------------------------------------------------------- + +typename wkb::multipoint_t wkb::read_multipoint(const byte_vector& buffer, int& index, uint8_t order) +{ + uint32_t count; + utils::copy_from(buffer, index, count, order); + + multipoint_t mp; + mp.resize(count); + + for (auto& point : mp) { + const byte_t new_order = buffer[index]; + index++; + + uint32_t type; + utils::copy_from(buffer, index, type, new_order); + + point = read_point(buffer, index, new_order); + } + + return mp; +} + +// ---------------------------------------------------------------------------- + +typename wkb::multilinestring_t wkb::read_multilinestring(const byte_vector& buffer, int& index, uint8_t order) +{ + uint32_t count; + utils::copy_from(buffer, index, count, order); + + multilinestring_t ml; + ml.resize(count); + for (auto& line : ml) { + const byte_t new_order = buffer[index]; + index++; + + uint32_t type; + utils::copy_from(buffer, index, type, new_order); + + line = read_linestring(buffer, index, new_order); + } + + return ml; +} + +// ---------------------------------------------------------------------------- + +typename wkb::multipolygon_t wkb::read_multipolygon(const byte_vector& buffer, int& index, uint8_t order) +{ + uint32_t count; + utils::copy_from(buffer, index, count, order); + + multipolygon_t mpl; + mpl.resize(count); + for (auto& polygon : mpl) { + const byte_t new_order = buffer[index]; + index++; + + uint32_t type; + utils::copy_from(buffer, index, type, new_order); + + polygon = read_polygon(buffer, index, new_order); + } + + return mpl; +} + +// ---------------------------------------------------------------------------- + +typename wkb::geometry wkb::read(const byte_vector& buffer) +{ + if (buffer.size() < 5) { + throw std::runtime_error("buffer reached end before encountering WKB"); + } + + int index = 0; + const byte_t order = buffer[index]; + index++; + + uint32_t type; + utils::copy_from(buffer, index, type, order); + + geometry g; + switch(type) { + case 1: g = read_point(buffer, index, order); break; + case 2: g = read_linestring(buffer, index, order); break; + case 3: g = read_polygon(buffer, index, order); break; + case 4: g = read_multipoint(buffer, index, order); break; + case 5: g = read_multilinestring(buffer, index, order); break; + case 6: g = read_multipolygon(buffer, index, order); break; + default: g = point_t{std::nan("0"), std::nan("0")}; break; + } + + return g; +} + +// ---------------------------------------------------------------------------- +// WKB Projection Visitor +// ---------------------------------------------------------------------------- + +bg::srs::dpar::parameters<> wkb::get_prj(uint32_t srid) { + switch(srid) { + case 5070: + return epsg5070; + case 3857: + return epsg3857; + default: + return bg::projections::detail::epsg_to_parameters(srid); + } +} + +// ---------------------------------------------------------------------------- + +geojson::geometry wkb::wgs84::operator()(point_t& g) +{ + if (this->srs == 4326) { + return geojson::coordinate_t(g.get<0>(), g.get<1>()); + } + + geojson::coordinate_t h; + this->tr.forward(g, h); + return h; +} + +// ---------------------------------------------------------------------------- + +geojson::geometry wkb::wgs84::operator()(linestring_t& g) +{ + geojson::linestring_t h; + + if (this->srs == 4326) { + h.reserve(g.size()); + for (auto&& gg : g) { + h.emplace_back( + std::move(gg.get<0>()), + std::move(gg.get<1>()) + ); + } + } else { + this->tr.forward(g, h); + } + return h; +} + +// ---------------------------------------------------------------------------- + +geojson::geometry wkb::wgs84::operator()(polygon_t& g) +{ + geojson::polygon_t h; + + if(this->srs == 4326) { + h.outer().reserve(g.outer().size()); + for (auto&& gg : g.outer()) { + h.outer().emplace_back( + std::move(gg.get<0>()), + std::move(gg.get<1>()) + ); + } + + h.inners().resize(g.inners().size()); + auto&& inner_g = g.inners().begin(); + auto&& inner_h = h.inners().begin(); + for (; inner_g != g.inners().end(); inner_g++, inner_h++) { + inner_h->reserve(inner_g->size()); + for (auto&& gg : *inner_g) { + inner_h->emplace_back( + std::move(gg.get<0>()), + std::move(gg.get<1>()) + ); + } + } + } else { + this->tr.forward(g, h); + } + return h; +} + +// ---------------------------------------------------------------------------- + +geojson::geometry wkb::wgs84::operator()(multipoint_t& g) +{ + geojson::multipoint_t h; + + if (this->srs == 4326) { + h.reserve(g.size()); + for (auto&& gg : g) { + h.emplace_back( + std::move(gg.get<0>()), + std::move(gg.get<1>()) + ); + } + } else { + this->tr.forward(g, h); + } + + return h; +} + +// ---------------------------------------------------------------------------- + +geojson::geometry wkb::wgs84::operator()(multilinestring_t& g) +{ + geojson::multilinestring_t h; + + if (this->srs == 4326) { + h.resize(g.size()); + auto&& line_g = g.begin(); + auto&& line_h = h.begin(); + for(; line_g != g.end(); line_g++, line_h++) { + *line_h = std::move( + boost::get(this->operator()(*line_g)) + ); + } + } else { + this->tr.forward(g, h); + } + + return h; +} + +// ---------------------------------------------------------------------------- + +geojson::geometry wkb::wgs84::operator()(multipolygon_t& g) +{ + geojson::multipolygon_t h; + + if (this->srs == 4326) { + h.resize(g.size()); + auto&& polygon_g = g.begin(); + auto&& polygon_h = h.begin(); + for (; polygon_g != g.end(); polygon_g++, polygon_h++) { + *polygon_h = std::move( + boost::get(this->operator()(*polygon_g)) + ); + } + } else { + this->tr.forward(g, h); + } + + return h; +} diff --git a/test/geopackage/WKB_Test.cpp b/test/geopackage/WKB_Test.cpp index b03b5cf370..f0732a979c 100644 --- a/test/geopackage/WKB_Test.cpp +++ b/test/geopackage/WKB_Test.cpp @@ -107,6 +107,35 @@ class WKB_Test : public ::testing::Test 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x3e, 0x40, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40 }; + + // POLYGON ((0 0, 5 0, 5 5, 0 5, 0 0), (2 2, 3 2, 3 3, 2 3, 2 2)) + this->wkb_bytes["polygon_with_holes"] = { + 0x01, + 0x03, 0x00, 0x00, 0x00, + 0x02, 0x00, 0x00, 0x00, + 0x05, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x14, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x14, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x14, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x14, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x05, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40 + }; } void TearDown() override {} @@ -162,6 +191,41 @@ TEST_F(WKB_Test, wkb_polygon_test) } } +TEST_F(WKB_Test, wkb_polygon_with_holes_test) +{ + const wkb::geometry geom = wkb::read(this->wkb_bytes["polygon_with_holes"]); + EXPECT_EQ(geom.which() + 1, 3); + + const wkb::polygon_t& poly = boost::get(geom); + ASSERT_EQ(poly.inners().size(), 1); + + const std::vector> expected_outer = { + {0, 0}, + {5, 0}, + {5, 5}, + {0, 5}, + {0, 0} + }; + + for (int i = 0; i < expected_outer.size(); i++) { + EXPECT_NEAR(poly.outer()[i].get<0>(), expected_outer[i].first, 0.0001); + EXPECT_NEAR(poly.outer()[i].get<1>(), expected_outer[i].second, 0.0001); + } + + const std::vector> expected_inner = { + {2, 2}, + {3, 2}, + {3, 3}, + {2, 3}, + {2, 2} + }; + + for (int i = 0; i < expected_inner.size(); i++) { + EXPECT_NEAR(poly.inners().at(0)[i].get<0>(), expected_inner[i].first, 0.0001); + EXPECT_NEAR(poly.inners().at(0)[i].get<1>(), expected_inner[i].second, 0.0001); + } +} + TEST_F(WKB_Test, wkb_multipoint_test) { const wkb::geometry geom = wkb::read(this->wkb_bytes["multipoint"]);