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