From 3cebba327034ed22919652d9e1ed187ac39a7f50 Mon Sep 17 00:00:00 2001 From: Johannes Schauer Marin Rodrigues Date: Sat, 16 Oct 2021 15:38:17 +0200 Subject: [PATCH] Add fallback_strategy to tinshift transform - this bumps format_version of tinshift JSON to 1.1 for the new field fallback_strategy - the default behaviour without that field is retained but if fallback_strategy is set to "nearest", then points that do not fall into any of the triangles will be transformed according to the nearest triangle --- data/Makefile.am | 1 + data/tests/tinshift_fallback_strategy.json | 15 +++ data/triangulation.schema.json | 9 +- .../operations/transformations/tinshift.rst | 11 ++ src/transformations/tinshift.hpp | 10 ++ src/transformations/tinshift_impl.hpp | 123 +++++++++++++++++- test/gie/tinshift.gie | 6 + test/unit/test_tinshift.cpp | 67 ++++++++++ 8 files changed, 240 insertions(+), 2 deletions(-) create mode 100644 data/tests/tinshift_fallback_strategy.json diff --git a/data/Makefile.am b/data/Makefile.am index 4990c45cee..f5f3e75f51 100644 --- a/data/Makefile.am +++ b/data/Makefile.am @@ -129,6 +129,7 @@ EXTRA_DIST = proj.ini GL27 nad.lst nad27 nad83 \ tests/tinshift_crs_implicit.json \ tests/tinshift_simplified_kkj_etrs.json \ tests/tinshift_simplified_n60_n2000.json \ + tests/tinshift_fallback_strategy.json \ generate_proj_db.cmake sql_filelist.cmake \ $(SQL_ORDERED_LIST) diff --git a/data/tests/tinshift_fallback_strategy.json b/data/tests/tinshift_fallback_strategy.json new file mode 100644 index 0000000000..91b57af9dd --- /dev/null +++ b/data/tests/tinshift_fallback_strategy.json @@ -0,0 +1,15 @@ +{ + "file_type": "triangulation_file", + "format_version": "1.1", + "fallback_strategy": "nearest", + "transformed_components": [ "horizontal" ], + "vertices_columns": [ "source_x", "source_y", "target_x", "target_y" ], + "triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ], + "vertices": [ + [0, 0, 0, 0], + [1, 0, 2, 0], + [1, 1, 2, 2], + [0, 1, 0, 2] + ], + "triangles": [ [0, 1, 2], [0, 2, 3] ] +} diff --git a/data/triangulation.schema.json b/data/triangulation.schema.json index d3f8004d5b..28ced040a0 100644 --- a/data/triangulation.schema.json +++ b/data/triangulation.schema.json @@ -13,7 +13,7 @@ "format_version": { "type": "string", "enum": [ - "1.0" + "1.0", "1.1" ] }, "name": { @@ -28,6 +28,13 @@ "$ref": "#/definitions/datetime", "description": "The date on which this version of the triangulation was published (or possibly the date on which it takes effect?)" }, + "fallback_strategy": { + "type": "string", + "enum": [ + "none", + "nearest" + ] + }, "license": { "type": "string", "description": "License under which the file is published" diff --git a/docs/source/operations/transformations/tinshift.rst b/docs/source/operations/transformations/tinshift.rst index ee7ae4a6c7..efc85a0a7a 100644 --- a/docs/source/operations/transformations/tinshift.rst +++ b/docs/source/operations/transformations/tinshift.rst @@ -166,6 +166,17 @@ transformed_components vertical component is transformed. +fallback_strategy + String identifying how to treat points that do not fall into any of the + specified triangles. This item is available for ``format_version`` >= 1.1. + Possible values are ``none`` and ``nearest``. The default is ``none`` and + signifies, that points that fall outside the specified triangles are not + transformed. This is also the behaviour for ``format_version`` before 1.1. + If ``fallback_strategy`` is set to ``nearest``, then points that do not fall + into any triangle are transformed according to the triangle closest to them + by euclidean distance. + + vertices_columns Specify the name of the columns of the rows in the ``vertices`` array. There must be exactly as many elements in ``vertices_columns`` as in a diff --git a/src/transformations/tinshift.hpp b/src/transformations/tinshift.hpp index cc8771b307..31fce139fc 100644 --- a/src/transformations/tinshift.hpp +++ b/src/transformations/tinshift.hpp @@ -53,6 +53,11 @@ namespace TINSHIFT_NAMESPACE { +enum FallbackStrategy { + FALLBACK_NONE = 0, + FALLBACK_NEAREST = 1, +}; + using json = nlohmann::json; // --------------------------------------------------------------------------- @@ -93,6 +98,10 @@ class TINShiftFile { */ const std::string &publicationDate() const { return mPublicationDate; } + const enum FallbackStrategy &fallbackStrategy() const { + return mFallbackStrategy; + } + /** Basic information on the agency responsible for the model. */ struct Authority { std::string name{}; @@ -204,6 +213,7 @@ class TINShiftFile { std::string mLicense{}; std::string mDescription{}; std::string mPublicationDate{}; + enum FallbackStrategy mFallbackStrategy {}; Authority mAuthority{}; std::vector mLinks{}; std::string mInputCRS{}; diff --git a/src/transformations/tinshift_impl.hpp b/src/transformations/tinshift_impl.hpp index 90f2b2f8c0..453fc813f2 100644 --- a/src/transformations/tinshift_impl.hpp +++ b/src/transformations/tinshift_impl.hpp @@ -93,6 +93,22 @@ std::unique_ptr TINShiftFile::parse(const std::string &text) { tinshiftFile->mDescription = getOptString(j, "description"); tinshiftFile->mPublicationDate = getOptString(j, "publication_date"); + tinshiftFile->mFallbackStrategy = FALLBACK_NONE; + if (j.contains("fallback_strategy")) { + if (tinshiftFile->mFormatVersion != "1.1") { + throw ParsingException( + "fallback_strategy needs format_version 1.1"); + } + const auto fallback_strategy = getOptString(j, "fallback_strategy"); + if (fallback_strategy == "nearest") { + tinshiftFile->mFallbackStrategy = FALLBACK_NEAREST; + } else if (fallback_strategy == "none") { + tinshiftFile->mFallbackStrategy = FALLBACK_NONE; + } else { + throw ParsingException("invalid fallback_strategy"); + } + } + if (j.contains("authority")) { const json jAuthority = j["authority"]; if (!jAuthority.is_object()) { @@ -410,6 +426,28 @@ Evaluator::Evaluator(std::unique_ptr &&fileIn) // --------------------------------------------------------------------------- +static inline double sqr(double x) { return x * x; } +static inline double squared_distance(double x1, double y1, double x2, + double y2) { + return sqr(x1 - x2) + sqr(y1 - y2); +} +static double distance_point_segment(double x, double y, double x1, double y1, + double x2, double y2, double dist12) { + // squared distance of point x/y to line segment x1/y1 -- x2/y2 + double t = ((x - x1) * (x2 - x1) + (y - y1) * (y2 - y1)) / dist12; + if (t <= 0) { + // closest to x1/y1 + return squared_distance(x, y, x1, y1); + } + if (t >= 1) { + // closest to y2/y2 + return squared_distance(x, y, x2, y2); + } + + // closest to line segment x1/y1 -- x2/y2 + return squared_distance(x, y, x1 + t * (x2 - x1), y1 + t * (y2 - y1)); +} + static const TINShiftFile::VertexIndices * FindTriangle(const TINShiftFile &file, const NS_PROJ::QuadTree::QuadTree &quadtree, @@ -453,7 +491,90 @@ FindTriangle(const TINShiftFile &file, } } } - return nullptr; + if (file.fallbackStrategy() != FALLBACK_NEAREST) { + return nullptr; + } + // find triangle with the shortest squared distance + // + // TODO: extend quadtree to support nearest neighbor search + double closest_dist = std::numeric_limits::infinity(); + double closest_dist2 = std::numeric_limits::infinity(); + size_t closest_i = 0; + for (size_t i = 0; i < triangles.size(); ++i) { + const auto &triangle = triangles[i]; + const unsigned i1 = triangle.idx1; + const unsigned i2 = triangle.idx2; + const unsigned i3 = triangle.idx3; + const double x1 = vertices[i1 * colCount + idxX]; + const double y1 = vertices[i1 * colCount + idxY]; + const double x2 = vertices[i2 * colCount + idxX]; + const double y2 = vertices[i2 * colCount + idxY]; + const double x3 = vertices[i3 * colCount + idxX]; + const double y3 = vertices[i3 * colCount + idxY]; + + // don't check this triangle if the query point plusminus the + // currently closest found distance is outside the triangle's AABB + if (x + closest_dist < std::min(x1, std::min(x2, x3)) || + x - closest_dist > std::max(x1, std::max(x2, x3)) || + y + closest_dist < std::min(y1, std::min(y2, y3)) || + y - closest_dist > std::max(y1, std::max(y2, y3))) { + continue; + } + + double dist12 = squared_distance(x1, y1, x2, y2); + double dist23 = squared_distance(x2, y2, x3, y3); + double dist13 = squared_distance(x1, y1, x3, y3); + if (dist12 < EPS || dist23 < EPS || dist13 < EPS) { + // do not use degenerate triangles + continue; + } + double dist2; + // we don't know whether the points of the triangle are given + // clockwise or counter-clockwise, so we have to check the distance of + // the point to all three sides of the triangle + dist2 = distance_point_segment(x, y, x1, y1, x2, y2, dist12); + if (dist2 < closest_dist2) { + closest_dist2 = dist2; + closest_dist = sqrt(dist2); + closest_i = i; + } + dist2 = distance_point_segment(x, y, x2, y2, x3, y3, dist23); + if (dist2 < closest_dist2) { + closest_dist2 = dist2; + closest_dist = sqrt(dist2); + closest_i = i; + } + dist2 = distance_point_segment(x, y, x1, y1, x3, y3, dist13); + if (dist2 < closest_dist2) { + closest_dist2 = dist2; + closest_dist = sqrt(dist2); + closest_i = i; + } + } + if (std::isinf(closest_dist)) { + // nothing was found due to empty triangle list or only degenerate + // triangles + return nullptr; + } + const auto &triangle = triangles[closest_i]; + const unsigned i1 = triangle.idx1; + const unsigned i2 = triangle.idx2; + const unsigned i3 = triangle.idx3; + const double x1 = vertices[i1 * colCount + idxX]; + const double y1 = vertices[i1 * colCount + idxY]; + const double x2 = vertices[i2 * colCount + idxX]; + const double y2 = vertices[i2 * colCount + idxY]; + const double x3 = vertices[i3 * colCount + idxX]; + const double y3 = vertices[i3 * colCount + idxY]; + const double det_T = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3); + if (det_T < EPS) { + // the nearest triangle is degenerate + return nullptr; + } + lambda1 = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) / det_T; + lambda2 = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) / det_T; + lambda3 = 1 - lambda1 - lambda2; + return ▵ } // --------------------------------------------------------------------------- diff --git a/test/gie/tinshift.gie b/test/gie/tinshift.gie index f16c16bbfa..3e135c2626 100644 --- a/test/gie/tinshift.gie +++ b/test/gie/tinshift.gie @@ -47,4 +47,10 @@ accept 3210000.0000 6700000.0000 10.0 expect 3210000.0000 6700000.0000 10.2886 roundtrip 1 +# Test fallback strategy nearest +operation +proj=tinshift +file=tests/tinshift_fallback_strategy.json +accept 2 3 +expect 4 6 +roundtrip 1 + diff --git a/test/unit/test_tinshift.cpp b/test/unit/test_tinshift.cpp index 96da6e2572..876d41aae9 100644 --- a/test/unit/test_tinshift.cpp +++ b/test/unit/test_tinshift.cpp @@ -65,6 +65,7 @@ TEST(tinshift, basic) { EXPECT_EQ(f->formatVersion(), "1.0"); EXPECT_EQ(f->inputCRS(), "EPSG:2393"); EXPECT_EQ(f->outputCRS(), "EPSG:3067"); + EXPECT_EQ(f->fallbackStrategy(), FALLBACK_NONE); auto eval = Evaluator(std::move(f)); double x_out = 0; @@ -206,6 +207,72 @@ TEST(tinshift, basic) { EXPECT_EQ(y_out, 0.75); EXPECT_EQ(z_out, 1000.0); } + + // invalid fallback_strategy field with 1.0 version + { + auto j(jMinValid); + j["fallback_strategy"] = "none"; + EXPECT_THROW(TINShiftFile::parse(j.dump()), ParsingException); + } + + // invalid fallback_strategy field with 1.1 version + { + auto j(jMinValid); + j["format_version"] = "1.1"; + j["fallback_strategy"] = "invalid"; + EXPECT_THROW(TINShiftFile::parse(j.dump()), ParsingException); + } + + // fail with no triangles and fallback nearest + { + auto j(jMinValid); + j["format_version"] = "1.1"; + j["fallback_strategy"] = "nearest"; + j["triangles"] = json::array(); // empty + + auto f = TINShiftFile::parse(j.dump()); + auto eval = Evaluator(std::move(f)); + double x_out = 0; + double y_out = 0; + double z_out = 0; + + EXPECT_FALSE(eval.forward(1.0, 1.0, 1.0, x_out, y_out, z_out)); + } + + // fail with only degenerate triangles (one size is zero length) and + // fallback nearest + { + auto j(jMinValid); + j["format_version"] = "1.1"; + j["fallback_strategy"] = "nearest"; + j["vertices"] = {{0, 0, 101, 101}, {0, 1, 100, 101}, {0, 1, 100, 100}}; + + auto f = TINShiftFile::parse(j.dump()); + auto eval = Evaluator(std::move(f)); + double x_out = 0; + double y_out = 0; + double z_out = 0; + + EXPECT_FALSE(eval.forward(1.0, 1.0, 1.0, x_out, y_out, z_out)); + } + + // fail with only degenerate triangles (two angles are 0°) and fallback + // nearest + { + auto j(jMinValid); + j["format_version"] = "1.1"; + j["fallback_strategy"] = "nearest"; + j["vertices"] = { + {0, 0, 101, 101}, {0, 0.5, 100, 101}, {0, 1, 100, 100}}; + + auto f = TINShiftFile::parse(j.dump()); + auto eval = Evaluator(std::move(f)); + double x_out = 0; + double y_out = 0; + double z_out = 0; + + EXPECT_FALSE(eval.forward(1.0, 1.0, 1.0, x_out, y_out, z_out)); + } } } // namespace