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

@@ -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)
 

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] ]
+}

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"

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

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<Link> mLinks{};
     std::string mInputCRS{};

src/transformations/tinshift_impl.hpp

@@ -93,6 +93,22 @@ std::unique_ptr<TINShiftFile> 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<TINShiftFile> &&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<unsigned> &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<double>::infinity();
+    double closest_dist2 = std::numeric_limits<double>::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 &triangle;
 }
 
 // ---------------------------------------------------------------------------

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
+
 </gie-strict>

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