Add a +proj=tinshift for triangulation-based transformations

Implements RFC-6

8 files changed, 1261 insertions, 0 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 86444df3..68aeadf2 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -201,6 +201,10 @@ libproj_la_SOURCES = \
 	transformations/defmodel.hpp \
 	transformations/defmodel_exceptions.hpp \
 	transformations/defmodel_impl.hpp \
+	transformations/tinshift.cpp \
+	transformations/tinshift.hpp \
+	transformations/tinshift_exceptions.hpp \
+	transformations/tinshift_impl.hpp \
 	\
 	aasincos.cpp adjlon.cpp \
 	dmstor.cpp auth.cpp \
@@ -213,6 +217,7 @@ libproj_la_SOURCES = \
 	release.cpp gauss.cpp \
 	fileapi.cpp \
 	generic_inverse.cpp \
+	quadtree.hpp \
 	\
 	datums.cpp datum_set.cpp transform.cpp \
 	utils.cpp \
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 3d41fe9a..1cd7f421 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -186,6 +186,7 @@ set(SRC_LIBPROJ_TRANSFORMATIONS
   transformations/vgridshift.cpp
   transformations/xyzgridshift.cpp
   transformations/defmodel.cpp
+  transformations/tinshift.cpp
 )
 
 set(SRC_LIBPROJ_ISO19111
diff --git a/src/pj_list.h b/src/pj_list.h
index c7e6d209..5b91af9b 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -150,6 +150,7 @@ PROJ_HEAD(gstmerc, "Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reuni
 PROJ_HEAD(tcc, "Transverse Central Cylindrical")
 PROJ_HEAD(tcea, "Transverse Cylindrical Equal Area")
 PROJ_HEAD(times, "Times Projection")
+PROJ_HEAD(tinshift, "Triangulation based transformation")
 PROJ_HEAD(tissot, "Tissot Conic")
 PROJ_HEAD(tmerc, "Transverse Mercator")
 PROJ_HEAD(tobmerc, "Tobler-Mercator")
diff --git a/src/quadtree.hpp b/src/quadtree.hpp
new file mode 100644
index 00000000..b744c83d
--- /dev/null
+++ b/src/quadtree.hpp
@@ -0,0 +1,257 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Implementation of quadtree building and searching functions.
+ *           Derived from shapelib, mapserver and GDAL implementations
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 1999-2008, Frank Warmerdam
+ * Copyright (c) 2008-2020, Even Rouault <even dot rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#ifndef QUADTREE_HPP
+#define QUADTREE_HPP
+
+#include "proj/util.hpp"
+
+#include <functional>
+#include <vector>
+
+//! @cond Doxygen_Suppress
+
+NS_PROJ_START
+
+namespace QuadTree {
+
+/* -------------------------------------------------------------------- */
+/*      If the following is 0.5, psNodes will be split in half.  If it  */
+/*      is 0.6 then each apSubNode will contain 60% of the parent       */
+/*      psNode, with 20% representing overlap.  This can be help to     */
+/*      prevent small objects on a boundary from shifting too high      */
+/*      up the hQuadTree.                                               */
+/* -------------------------------------------------------------------- */
+constexpr double DEFAULT_SPLIT_RATIO = 0.55;
+
+/** Describe a rectangle */
+struct RectObj {
+    double minx = 0; /**< Minimum x */
+    double miny = 0; /**< Minimum y */
+    double maxx = 0; /**< Maximum x */
+    double maxy = 0; /**< Maximum y */
+
+    /* Returns whether this rectangle is contained by other */
+    inline bool isContainedBy(const RectObj &other) const {
+        return minx >= other.minx && maxx <= other.maxx && miny >= other.miny &&
+               maxy <= other.maxy;
+    }
+
+    /* Returns whether this rectangles overlaps other */
+    inline bool overlaps(const RectObj &other) const {
+        return minx <= other.maxx && maxx >= other.minx && miny <= other.maxy &&
+               maxy >= other.miny;
+    }
+
+    /* Returns whether this rectangles contains the specified point */
+    inline bool contains(double x, double y) const {
+        return minx <= x && maxx >= x && miny <= y && maxy >= y;
+    }
+
+    /* Return whether this rectangles is different from other */
+    inline bool operator!=(const RectObj &other) const {
+        return minx != other.minx || miny != other.miny || maxx != other.maxx ||
+               maxy != other.maxy;
+    }
+};
+
+/** Quadtree */
+template <class Feature> class QuadTree {
+
+    struct Node {
+        RectObj rect{}; /* area covered by this psNode */
+
+        /* list of shapes stored at this node. */
+        std::vector<std::pair<Feature, RectObj>> features{};
+
+        std::vector<Node> subnodes{};
+
+        explicit Node(const RectObj &rectIn) : rect(rectIn) {}
+    };
+
+    Node root{};
+    unsigned nBucketCapacity = 8;
+    double dfSplitRatio = DEFAULT_SPLIT_RATIO;
+
+  public:
+    /** Construct a new quadtree with the global bounds of all objects to be
+    * inserted */
+    explicit QuadTree(const RectObj &globalBounds) : root(globalBounds) {}
+
+    /** Add a new feature, with its bounds specified in featureBounds */
+    void insert(const Feature &feature, const RectObj &featureBounds) {
+        insert(root, feature, featureBounds);
+    }
+
+#ifdef UNUSED
+    /** Retrieve all features whose bounds intersects aoiRect */
+    void
+    search(const RectObj &aoiRect,
+           std::vector<std::reference_wrapper<const Feature>> &features) const {
+        search(root, aoiRect, features);
+    }
+#endif
+
+    /** Retrieve all features whose bounds contains (x,y) */
+    void search(double x, double y, std::vector<Feature> &features) const {
+        search(root, x, y, features);
+    }
+
+  private:
+    void splitBounds(const RectObj &in, RectObj &out1, RectObj &out2) {
+        // The output bounds will be very similar to the input bounds,
+        // so just copy over to start.
+        out1 = in;
+        out2 = in;
+
+        // Split in X direction.
+        if ((in.maxx - in.minx) > (in.maxy - in.miny)) {
+            const double range = in.maxx - in.minx;
+
+            out1.maxx = in.minx + range * dfSplitRatio;
+            out2.minx = in.maxx - range * dfSplitRatio;
+        }
+
+        // Otherwise split in Y direction.
+        else {
+            const double range = in.maxy - in.miny;
+
+            out1.maxy = in.miny + range * dfSplitRatio;
+            out2.miny = in.maxy - range * dfSplitRatio;
+        }
+    }
+
+    void insert(Node &node, const Feature &feature,
+                const RectObj &featureBounds) {
+        if (node.subnodes.empty()) {
+            // If we have reached the max bucket capacity, try to insert
+            // in a subnode if possible.
+            if (node.features.size() >= nBucketCapacity) {
+                RectObj half1;
+                RectObj half2;
+                RectObj quad1;
+                RectObj quad2;
+                RectObj quad3;
+                RectObj quad4;
+
+                splitBounds(node.rect, half1, half2);
+                splitBounds(half1, quad1, quad2);
+                splitBounds(half2, quad3, quad4);
+
+                if (node.rect != quad1 && node.rect != quad2 &&
+                    node.rect != quad3 && node.rect != quad4 &&
+                    (featureBounds.isContainedBy(quad1) ||
+                     featureBounds.isContainedBy(quad2) ||
+                     featureBounds.isContainedBy(quad3) ||
+                     featureBounds.isContainedBy(quad4))) {
+                    node.subnodes.reserve(4);
+                    node.subnodes.emplace_back(Node(quad1));
+                    node.subnodes.emplace_back(Node(quad2));
+                    node.subnodes.emplace_back(Node(quad3));
+                    node.subnodes.emplace_back(Node(quad4));
+
+                    auto features = std::move(node.features);
+                    node.features.clear();
+                    for (auto &pair : features) {
+                        insert(node, pair.first, pair.second);
+                    }
+
+                    /* recurse back on this psNode now that it has apSubNodes */
+                    insert(node, feature, featureBounds);
+                    return;
+                }
+            }
+        } else {
+            // If we have sub nodes, then consider whether this object will
+            // fit in them.
+            for (auto &subnode : node.subnodes) {
+                if (featureBounds.isContainedBy(subnode.rect)) {
+                    insert(subnode, feature, featureBounds);
+                    return;
+                }
+            }
+        }
+
+        // If none of that worked, just add it to this nodes list.
+        node.features.push_back(
+            std::pair<Feature, RectObj>(feature, featureBounds));
+    }
+
+#ifdef UNUSED
+    void
+    search(const Node &node, const RectObj &aoiRect,
+           std::vector<std::reference_wrapper<const Feature>> &features) const {
+        // Does this node overlap the area of interest at all?  If not,
+        // return without adding to the list at all.
+        if (!node.rect.overlaps(aoiRect))
+            return;
+
+        // Add the local features to the list.
+        for (const auto &pair : node.features) {
+            if (pair.second.overlaps(aoiRect)) {
+                features.push_back(
+                    std::reference_wrapper<const Feature>(pair.first));
+            }
+        }
+
+        // Recurse to subnodes if they exist.
+        for (const auto &subnode : node.subnodes) {
+            search(subnode, aoiRect, features);
+        }
+    }
+#endif
+
+    void search(const Node &node, double x, double y,
+                std::vector<Feature> &features) const {
+        // Does this node overlap the area of interest at all?  If not,
+        // return without adding to the list at all.
+        if (!node.rect.contains(x, y))
+            return;
+
+        // Add the local features to the list.
+        for (const auto &pair : node.features) {
+            if (pair.second.contains(x, y)) {
+                features.push_back(pair.first);
+            }
+        }
+
+        // Recurse to subnodes if they exist.
+        for (const auto &subnode : node.subnodes) {
+            search(subnode, x, y, features);
+        }
+    }
+};
+
+} // namespace QuadTree
+
+NS_PROJ_END
+
+//! @endcond
+
+#endif //  QUADTREE_HPP
diff --git a/src/transformations/tinshift.cpp b/src/transformations/tinshift.cpp
new file mode 100644
index 00000000..96e0ea4f
--- /dev/null
+++ b/src/transformations/tinshift.cpp
@@ -0,0 +1,133 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Functionality related to triangulation transformation
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#define PJ_LIB__
+#define PROJ_COMPILATION
+
+#include "tinshift.hpp"
+#include "filemanager.hpp"
+#include "proj_internal.h"
+
+PROJ_HEAD(tinshift, "Triangulation based transformation");
+
+using namespace TINSHIFT_NAMESPACE;
+
+namespace {
+
+struct tinshiftData {
+    std::unique_ptr<Evaluator> evaluator{};
+
+    tinshiftData() = default;
+
+    tinshiftData(const tinshiftData &) = delete;
+    tinshiftData &operator=(const tinshiftData &) = delete;
+};
+
+} // namespace
+
+static PJ *destructor(PJ *P, int errlev) {
+    if (nullptr == P)
+        return nullptr;
+
+    auto Q = static_cast<struct tinshiftData *>(P->opaque);
+    delete Q;
+    P->opaque = nullptr;
+
+    return pj_default_destructor(P, errlev);
+}
+
+static PJ_COORD tinshift_forward_4d(PJ_COORD in, PJ *P) {
+    auto *Q = (struct tinshiftData *)P->opaque;
+
+    PJ_COORD out = in;
+    if (!Q->evaluator->forward(in.xyz.x, in.xyz.y, in.xyz.z, out.xyz.x,
+                               out.xyz.y, out.xyz.z)) {
+        return proj_coord_error();
+    }
+    return out;
+}
+
+static PJ_COORD tinshift_reverse_4d(PJ_COORD in, PJ *P) {
+    auto *Q = (struct tinshiftData *)P->opaque;
+
+    PJ_COORD out = in;
+    if (!Q->evaluator->inverse(in.xyz.x, in.xyz.y, in.xyz.z, out.xyz.x,
+                               out.xyz.y, out.xyz.z)) {
+        return proj_coord_error();
+    }
+    return out;
+}
+
+PJ *TRANSFORMATION(tinshift, 1) {
+
+    const char *filename = pj_param(P->ctx, P->params, "sfile").s;
+    if (!filename) {
+        proj_log_error(P, "tinshift: +file= should be specified.");
+        return destructor(P, PJD_ERR_NO_ARGS);
+    }
+
+    auto file = NS_PROJ::FileManager::open_resource_file(P->ctx, filename);
+    if (nullptr == file) {
+        proj_log_error(P, "tinshift: Cannot open %s", filename);
+        return destructor(P, PJD_ERR_INVALID_ARG);
+    }
+    file->seek(0, SEEK_END);
+    unsigned long long size = file->tell();
+    // Arbitrary threshold to avoid ingesting an arbitrarily large JSON file,
+    // that could be a denial of service risk. 10 MB should be sufficiently
+    // large for any valid use !
+    if (size > 10 * 1024 * 1024) {
+        proj_log_error(P, "tinshift: File %s too large", filename);
+        return destructor(P, PJD_ERR_INVALID_ARG);
+    }
+    file->seek(0);
+    std::string jsonStr;
+    jsonStr.resize(static_cast<size_t>(size));
+    if (file->read(&jsonStr[0], jsonStr.size()) != jsonStr.size()) {
+        proj_log_error(P, "tinshift: Cannot read %s", filename);
+        return destructor(P, PJD_ERR_INVALID_ARG);
+    }
+
+    auto Q = new tinshiftData();
+    P->opaque = (void *)Q;
+    P->destructor = destructor;
+
+    try {
+        Q->evaluator.reset(new Evaluator(TINShiftFile::parse(jsonStr)));
+    } catch (const std::exception &e) {
+        proj_log_error(P, "tinshift: invalid model: %s", e.what());
+        return destructor(P, PJD_ERR_INVALID_ARG);
+    }
+
+    P->destructor = destructor;
+    P->fwd4d = tinshift_forward_4d;
+    P->inv4d = tinshift_reverse_4d;
+    P->left = PJ_IO_UNITS_WHATEVER;
+    P->right = PJ_IO_UNITS_WHATEVER;
+
+    return P;
+}
diff --git a/src/transformations/tinshift.hpp b/src/transformations/tinshift.hpp
new file mode 100644
index 00000000..7cebab9f
--- /dev/null
+++ b/src/transformations/tinshift.hpp
@@ -0,0 +1,257 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Functionality related to TIN based transformations
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#ifndef TINSHIFT_HPP
+#define TINSHIFT_HPP
+
+#ifdef PROJ_COMPILATION
+#include "proj/internal/include_nlohmann_json.hpp"
+#else
+#include "nlohmann/json.hpp"
+#endif
+
+#include <algorithm>
+#include <cmath>
+#include <exception>
+#include <functional>
+#include <limits>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "quadtree.hpp"
+
+#ifndef TINSHIFT_NAMESPACE
+#define TINSHIFT_NAMESPACE TINShift
+#endif
+
+#include "tinshift_exceptions.hpp"
+
+namespace TINSHIFT_NAMESPACE {
+
+using json = nlohmann::json;
+
+// ---------------------------------------------------------------------------
+
+/** Content of a TINShift file. */
+class TINShiftFile {
+  public:
+    /** Parse the provided serialized JSON content and return an object.
+     *
+     * @throws ParsingException
+     */
+    static std::unique_ptr<TINShiftFile> parse(const std::string &text);
+
+    /** Get file type. Should always be "triangulation_file" */
+    const std::string &fileType() const { return mFileType; }
+
+    /** Get the version of the format. At time of writing, only "1.0" is known
+     */
+    const std::string &formatVersion() const { return mFormatVersion; }
+
+    /** Get brief descriptive name of the deformation model. */
+    const std::string &name() const { return mName; }
+
+    /** Get a string identifying the version of the deformation model.
+     * The format for specifying version is defined by the agency
+     * responsible for the deformation model. */
+    const std::string &version() const { return mVersion; }
+
+    /** Get a string identifying the license of the file.
+     * e.g "Create Commons Attribution 4.0 International" */
+    const std::string &license() const { return mLicense; }
+
+    /** Get a text description of the model. Intended to be longer than name()
+     */
+    const std::string &description() const { return mDescription; }
+
+    /** Get a text description of the model. Intended to be longer than name()
+     */
+    const std::string &publicationDate() const { return mPublicationDate; }
+
+    /** Basic information on the agency responsible for the model. */
+    struct Authority {
+        std::string name{};
+        std::string url{};
+        std::string address{};
+        std::string email{};
+    };
+
+    /** Get basic information on the agency responsible for the model. */
+    const Authority &authority() const { return mAuthority; }
+
+    /** Hyperlink related to the model. */
+    struct Link {
+        /** URL holding the information */
+        std::string href{};
+
+        /** Relationship to the dataset. e.g. "about", "source", "license",
+         * "metadata" */
+        std::string rel{};
+
+        /** Mime type */
+        std::string type{};
+
+        /** Description of the link */
+        std::string title{};
+    };
+
+    /** Get links to related information. */
+    const std::vector<Link> links() const { return mLinks; }
+
+    /** Get a string identifying the CRS of source coordinates in the
+     * vertices. Typically "EPSG:XXXX". If the transformation is for vertical
+     * component, this should be the code for a compound CRS (can be
+     * EPSG:XXXX+YYYY where XXXX is the code of the horizontal CRS and YYYY
+     * the code of the vertical CRS).
+     * For example, for the KKJ->ETRS89 transformation, this is EPSG:2393
+     * ("KKJ / Finland Uniform Coordinate System").
+     * The input coordinates are assumed to
+     * be passed in the "normalized for visualisation" / "GIS friendly" order,
+     * that is longitude, latitude for geographic coordinates and
+     * easting, northing for projected coordinates.
+     * This may be empty for unspecified CRS.
+     */
+    const std::string &inputCRS() const { return mInputCRS; }
+
+    /** Get a string identifying the CRS of target coordinates in the
+     * vertices. Typically "EPSG:XXXX". If the transformation is for vertical
+     * component, this should be the code for a compound CRS (can be
+     * EPSG:XXXX+YYYY where XXXX is the code of the horizontal CRS and YYYY
+     * the code of the vertical CRS).
+     * For example, for the KKJ->ETRS89 transformation, this is EPSG:3067
+     * ("ETRS89 / TM35FIN(E,N)").
+     * The output coordinates will be
+     * returned in the "normalized for visualisation" / "GIS friendly" order,
+     * that is longitude, latitude for geographic coordinates and
+     * easting, northing for projected coordinates.
+     * This may be empty for unspecified CRS.
+     */
+    const std::string &outputCRS() const { return mOutputCRS; }
+
+    /** Return whether horizontal coordinates are transformed. */
+    bool transformHorizontalComponent() const {
+        return mTransformHorizontalComponent;
+    }
+
+    /** Return whether vertical coordinates are transformed. */
+    bool transformVerticalComponent() const {
+        return mTransformVerticalComponent;
+    }
+
+    /** Indices of vertices of a triangle */
+    struct VertexIndices {
+        /** Index of first vertex */
+        unsigned idx1;
+        /** Index of second vertex */
+        unsigned idx2;
+        /** Index of third vertex */
+        unsigned idx3;
+    };
+
+    /** Return number of elements per vertex of vertices() */
+    unsigned verticesColumnCount() const { return mVerticesColumnCount; }
+
+    /** Return description of triangulation vertices.
+     * Each vertex is described by verticesColumnCount() consecutive values.
+     * They are respectively:
+     * - the source X value
+     * - the source Y value
+     * - (if transformHorizontalComponent() is true) the target X value
+     * - (if transformHorizontalComponent() is true) the target Y value
+     * - (if transformVerticalComponent() is true) the delta Z value (to go from
+     * source to target Z)
+     *
+     * X is assumed to be a longitude (in degrees) or easting value.
+     * Y is assumed to be a latitude (in degrees) or northing value.
+     */
+    const std::vector<double> &vertices() const { return mVertices; }
+
+    /** Return triangles*/
+    const std::vector<VertexIndices> &triangles() const { return mTriangles; }
+
+  private:
+    TINShiftFile() = default;
+
+    std::string mFileType{};
+    std::string mFormatVersion{};
+    std::string mName{};
+    std::string mVersion{};
+    std::string mLicense{};
+    std::string mDescription{};
+    std::string mPublicationDate{};
+    Authority mAuthority{};
+    std::vector<Link> mLinks{};
+    std::string mInputCRS{};
+    std::string mOutputCRS{};
+    bool mTransformHorizontalComponent = false;
+    bool mTransformVerticalComponent = false;
+    unsigned mVerticesColumnCount = 0;
+    std::vector<double> mVertices{};
+    std::vector<VertexIndices> mTriangles{};
+};
+
+// ---------------------------------------------------------------------------
+
+/** Class to evaluate the transformation of a coordinate */
+class Evaluator {
+  public:
+    /** Constructor. */
+    explicit Evaluator(std::unique_ptr<TINShiftFile> &&fileIn);
+
+    /** Get file */
+    const TINShiftFile &file() const { return *(mFile.get()); }
+
+    /** Evaluate displacement of a position given by (x,y,z,t) and
+     * return it in (x_out,y_out_,z_out).
+     */
+    bool forward(double x, double y, double z, double &x_out, double &y_out,
+                 double &z_out);
+
+    /** Apply inverse transformation. */
+    bool inverse(double x, double y, double z, double &x_out, double &y_out,
+                 double &z_out);
+
+  private:
+    std::unique_ptr<TINShiftFile> mFile;
+
+    // Reused between invokations to save memory allocations
+    std::vector<unsigned> mTriangleIndices{};
+
+    std::unique_ptr<NS_PROJ::QuadTree::QuadTree<unsigned>> mQuadTreeForward{};
+    std::unique_ptr<NS_PROJ::QuadTree::QuadTree<unsigned>> mQuadTreeInverse{};
+};
+
+// ---------------------------------------------------------------------------
+
+} // namespace TINSHIFT_NAMESPACE
+
+// ---------------------------------------------------------------------------
+
+#include "tinshift_impl.hpp"
+
+#endif // TINSHIFT_HPP
diff --git a/src/transformations/tinshift_exceptions.hpp b/src/transformations/tinshift_exceptions.hpp
new file mode 100644
index 00000000..d6d71302
--- /dev/null
+++ b/src/transformations/tinshift_exceptions.hpp
@@ -0,0 +1,52 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Functionality related to TIN based transformations
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#ifndef TINSHIFT_NAMESPACE
+#error "Should be included only by tinshift.hpp"
+#endif
+
+#include <exception>
+
+namespace TINSHIFT_NAMESPACE {
+
+// ---------------------------------------------------------------------------
+
+/** Parsing exception. */
+class ParsingException : public std::exception {
+  public:
+    explicit ParsingException(const std::string &msg) : msg_(msg) {}
+    const char *what() const noexcept override;
+
+  private:
+    std::string msg_;
+};
+
+const char *ParsingException::what() const noexcept { return msg_.c_str(); }
+
+// ---------------------------------------------------------------------------
+
+} // namespace DEFORMATON_MODEL_NAMESPACE
diff --git a/src/transformations/tinshift_impl.hpp b/src/transformations/tinshift_impl.hpp
new file mode 100644
index 00000000..90f2b2f8
--- /dev/null
+++ b/src/transformations/tinshift_impl.hpp
@@ -0,0 +1,555 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Functionality related to TIN based transformations
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#ifndef TINSHIFT_NAMESPACE
+#error "Should be included only by tinshift.hpp"
+#endif
+
+#include <limits>
+
+namespace TINSHIFT_NAMESPACE {
+
+// ---------------------------------------------------------------------------
+
+static std::string getString(const json &j, const char *key, bool optional) {
+    if (!j.contains(key)) {
+        if (optional) {
+            return std::string();
+        }
+        throw ParsingException(std::string("Missing \"") + key + "\" key");
+    }
+    const json v = j[key];
+    if (!v.is_string()) {
+        throw ParsingException(std::string("The value of \"") + key +
+                               "\" should be a string");
+    }
+    return v.get<std::string>();
+}
+
+static std::string getReqString(const json &j, const char *key) {
+    return getString(j, key, false);
+}
+
+static std::string getOptString(const json &j, const char *key) {
+    return getString(j, key, true);
+}
+
+// ---------------------------------------------------------------------------
+
+static json getArrayMember(const json &j, const char *key) {
+    if (!j.contains(key)) {
+        throw ParsingException(std::string("Missing \"") + key + "\" key");
+    }
+    const json obj = j[key];
+    if (!obj.is_array()) {
+        throw ParsingException(std::string("The value of \"") + key +
+                               "\" should be a array");
+    }
+    return obj;
+}
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<TINShiftFile> TINShiftFile::parse(const std::string &text) {
+    std::unique_ptr<TINShiftFile> tinshiftFile(new TINShiftFile());
+    json j;
+    try {
+        j = json::parse(text);
+    } catch (const std::exception &e) {
+        throw ParsingException(e.what());
+    }
+    if (!j.is_object()) {
+        throw ParsingException("Not an object");
+    }
+    tinshiftFile->mFileType = getReqString(j, "file_type");
+    tinshiftFile->mFormatVersion = getReqString(j, "format_version");
+    tinshiftFile->mName = getOptString(j, "name");
+    tinshiftFile->mVersion = getOptString(j, "version");
+    tinshiftFile->mLicense = getOptString(j, "license");
+    tinshiftFile->mDescription = getOptString(j, "description");
+    tinshiftFile->mPublicationDate = getOptString(j, "publication_date");
+
+    if (j.contains("authority")) {
+        const json jAuthority = j["authority"];
+        if (!jAuthority.is_object()) {
+            throw ParsingException("authority is not a object");
+        }
+        tinshiftFile->mAuthority.name = getOptString(jAuthority, "name");
+        tinshiftFile->mAuthority.url = getOptString(jAuthority, "url");
+        tinshiftFile->mAuthority.address = getOptString(jAuthority, "address");
+        tinshiftFile->mAuthority.email = getOptString(jAuthority, "email");
+    }
+
+    if (j.contains("links")) {
+        const json jLinks = j["links"];
+        if (!jLinks.is_array()) {
+            throw ParsingException("links is not an array");
+        }
+        for (const json &jLink : jLinks) {
+            if (!jLink.is_object()) {
+                throw ParsingException("links[] item is not an object");
+            }
+            Link link;
+            link.href = getOptString(jLink, "href");
+            link.rel = getOptString(jLink, "rel");
+            link.type = getOptString(jLink, "type");
+            link.title = getOptString(jLink, "title");
+            tinshiftFile->mLinks.emplace_back(std::move(link));
+        }
+    }
+    tinshiftFile->mInputCRS = getOptString(j, "input_crs");
+    tinshiftFile->mOutputCRS = getOptString(j, "output_crs");
+
+    const auto jTransformedComponents =
+        getArrayMember(j, "transformed_components");
+    for (const json &jComp : jTransformedComponents) {
+        if (!jComp.is_string()) {
+            throw ParsingException(
+                "transformed_components[] item is not a string");
+        }
+        const auto jCompStr = jComp.get<std::string>();
+        if (jCompStr == "horizontal") {
+            tinshiftFile->mTransformHorizontalComponent = true;
+        } else if (jCompStr == "vertical") {
+            tinshiftFile->mTransformVerticalComponent = true;
+        } else {
+            throw ParsingException("transformed_components[] = " + jCompStr +
+                                   " is not handled");
+        }
+    }
+
+    const auto jVerticesColumns = getArrayMember(j, "vertices_columns");
+    int sourceXCol = -1;
+    int sourceYCol = -1;
+    int sourceZCol = -1;
+    int targetXCol = -1;
+    int targetYCol = -1;
+    int targetZCol = -1;
+    int offsetZCol = -1;
+    for (size_t i = 0; i < jVerticesColumns.size(); ++i) {
+        const json &jColumn = jVerticesColumns[i];
+        if (!jColumn.is_string()) {
+            throw ParsingException("vertices_columns[] item is not a string");
+        }
+        const auto jColumnStr = jColumn.get<std::string>();
+        if (jColumnStr == "source_x") {
+            sourceXCol = static_cast<int>(i);
+        } else if (jColumnStr == "source_y") {
+            sourceYCol = static_cast<int>(i);
+        } else if (jColumnStr == "source_z") {
+            sourceZCol = static_cast<int>(i);
+        } else if (jColumnStr == "target_x") {
+            targetXCol = static_cast<int>(i);
+        } else if (jColumnStr == "target_y") {
+            targetYCol = static_cast<int>(i);
+        } else if (jColumnStr == "target_z") {
+            targetZCol = static_cast<int>(i);
+        } else if (jColumnStr == "offset_z") {
+            offsetZCol = static_cast<int>(i);
+        }
+    }
+    if (sourceXCol < 0) {
+        throw ParsingException(
+            "source_x must be specified in vertices_columns[]");
+    }
+    if (sourceYCol < 0) {
+        throw ParsingException(
+            "source_y must be specified in vertices_columns[]");
+    }
+    if (tinshiftFile->mTransformHorizontalComponent) {
+        if (targetXCol < 0) {
+            throw ParsingException(
+                "target_x must be specified in vertices_columns[]");
+        }
+        if (targetYCol < 0) {
+            throw ParsingException(
+                "target_y must be specified in vertices_columns[]");
+        }
+    }
+    if (tinshiftFile->mTransformVerticalComponent) {
+        if (offsetZCol >= 0) {
+            // do nothing
+        } else {
+            if (sourceZCol < 0) {
+                throw ParsingException("source_z or delta_z must be specified "
+                                       "in vertices_columns[]");
+            }
+            if (targetZCol < 0) {
+                throw ParsingException(
+                    "target_z must be specified in vertices_columns[]");
+            }
+        }
+    }
+
+    const auto jTrianglesColumns = getArrayMember(j, "triangles_columns");
+    int idxVertex1Col = -1;
+    int idxVertex2Col = -1;
+    int idxVertex3Col = -1;
+    for (size_t i = 0; i < jTrianglesColumns.size(); ++i) {
+        const json &jColumn = jTrianglesColumns[i];
+        if (!jColumn.is_string()) {
+            throw ParsingException("triangles_columns[] item is not a string");
+        }
+        const auto jColumnStr = jColumn.get<std::string>();
+        if (jColumnStr == "idx_vertex1") {
+            idxVertex1Col = static_cast<int>(i);
+        } else if (jColumnStr == "idx_vertex2") {
+            idxVertex2Col = static_cast<int>(i);
+        } else if (jColumnStr == "idx_vertex3") {
+            idxVertex3Col = static_cast<int>(i);
+        }
+    }
+    if (idxVertex1Col < 0) {
+        throw ParsingException(
+            "idx_vertex1 must be specified in triangles_columns[]");
+    }
+    if (idxVertex2Col < 0) {
+        throw ParsingException(
+            "idx_vertex2 must be specified in triangles_columns[]");
+    }
+    if (idxVertex3Col < 0) {
+        throw ParsingException(
+            "idx_vertex3 must be specified in triangles_columns[]");
+    }
+
+    const auto jVertices = getArrayMember(j, "vertices");
+    tinshiftFile->mVerticesColumnCount = 2;
+    if (tinshiftFile->mTransformHorizontalComponent)
+        tinshiftFile->mVerticesColumnCount += 2;
+    if (tinshiftFile->mTransformVerticalComponent)
+        tinshiftFile->mVerticesColumnCount += 1;
+
+    tinshiftFile->mVertices.reserve(tinshiftFile->mVerticesColumnCount *
+                                    jVertices.size());
+    for (const auto &jVertex : jVertices) {
+        if (!jVertex.is_array()) {
+            throw ParsingException("vertices[] item is not an array");
+        }
+        if (jVertex.size() != jVerticesColumns.size()) {
+            throw ParsingException(
+                "vertices[] item has not expected number of elements");
+        }
+        if (!jVertex[sourceXCol].is_number()) {
+            throw ParsingException("vertices[][] item is not a number");
+        }
+        tinshiftFile->mVertices.push_back(jVertex[sourceXCol].get<double>());
+        if (!jVertex[sourceYCol].is_number()) {
+            throw ParsingException("vertices[][] item is not a number");
+        }
+        tinshiftFile->mVertices.push_back(jVertex[sourceYCol].get<double>());
+        if (tinshiftFile->mTransformHorizontalComponent) {
+            if (!jVertex[targetXCol].is_number()) {
+                throw ParsingException("vertices[][] item is not a number");
+            }
+            tinshiftFile->mVertices.push_back(
+                jVertex[targetXCol].get<double>());
+            if (!jVertex[targetYCol].is_number()) {
+                throw ParsingException("vertices[][] item is not a number");
+            }
+            tinshiftFile->mVertices.push_back(
+                jVertex[targetYCol].get<double>());
+        }
+        if (tinshiftFile->mTransformVerticalComponent) {
+            if (offsetZCol >= 0) {
+                if (!jVertex[offsetZCol].is_number()) {
+                    throw ParsingException("vertices[][] item is not a number");
+                }
+                tinshiftFile->mVertices.push_back(
+                    jVertex[offsetZCol].get<double>());
+            } else {
+                if (!jVertex[sourceZCol].is_number()) {
+                    throw ParsingException("vertices[][] item is not a number");
+                }
+                const double sourceZ = jVertex[sourceZCol].get<double>();
+                if (!jVertex[targetZCol].is_number()) {
+                    throw ParsingException("vertices[][] item is not a number");
+                }
+                const double targetZ = jVertex[targetZCol].get<double>();
+                tinshiftFile->mVertices.push_back(targetZ - sourceZ);
+            }
+        }
+    }
+
+    const auto jTriangles = getArrayMember(j, "triangles");
+    tinshiftFile->mTriangles.reserve(jTriangles.size());
+    for (const auto &jTriangle : jTriangles) {
+        if (!jTriangle.is_array()) {
+            throw ParsingException("triangles[] item is not an array");
+        }
+        if (jTriangle.size() != jTrianglesColumns.size()) {
+            throw ParsingException(
+                "triangles[] item has not expected number of elements");
+        }
+
+        if (jTriangle[idxVertex1Col].type() != json::value_t::number_unsigned) {
+            throw ParsingException("triangles[][] item is not an integer");
+        }
+        const unsigned vertex1 = jTriangle[idxVertex1Col].get<unsigned>();
+        if (vertex1 >= jVertices.size()) {
+            throw ParsingException("Invalid value for a vertex index");
+        }
+
+        if (jTriangle[idxVertex2Col].type() != json::value_t::number_unsigned) {
+            throw ParsingException("triangles[][] item is not an integer");
+        }
+        const unsigned vertex2 = jTriangle[idxVertex2Col].get<unsigned>();
+        if (vertex2 >= jVertices.size()) {
+            throw ParsingException("Invalid value for a vertex index");
+        }
+
+        if (jTriangle[idxVertex3Col].type() != json::value_t::number_unsigned) {
+            throw ParsingException("triangles[][] item is not an integer");
+        }
+        const unsigned vertex3 = jTriangle[idxVertex3Col].get<unsigned>();
+        if (vertex3 >= jVertices.size()) {
+            throw ParsingException("Invalid value for a vertex index");
+        }
+
+        VertexIndices vi;
+        vi.idx1 = vertex1;
+        vi.idx2 = vertex2;
+        vi.idx3 = vertex3;
+        tinshiftFile->mTriangles.push_back(vi);
+    }
+
+    return tinshiftFile;
+}
+
+// ---------------------------------------------------------------------------
+
+static NS_PROJ::QuadTree::RectObj GetBounds(const TINShiftFile &file,
+                                            bool forward) {
+    NS_PROJ::QuadTree::RectObj rect;
+    rect.minx = std::numeric_limits<double>::max();
+    rect.miny = std::numeric_limits<double>::max();
+    rect.maxx = -std::numeric_limits<double>::max();
+    rect.maxy = -std::numeric_limits<double>::max();
+    const auto &vertices = file.vertices();
+    const unsigned colCount = file.verticesColumnCount();
+    const int idxX = file.transformHorizontalComponent() && !forward ? 2 : 0;
+    const int idxY = file.transformHorizontalComponent() && !forward ? 3 : 1;
+    for (size_t i = 0; i + colCount - 1 < vertices.size(); i += colCount) {
+        const double x = vertices[i + idxX];
+        const double y = vertices[i + idxY];
+        rect.minx = std::min(rect.minx, x);
+        rect.miny = std::min(rect.miny, y);
+        rect.maxx = std::max(rect.maxx, x);
+        rect.maxy = std::max(rect.maxy, y);
+    }
+    return rect;
+}
+
+// ---------------------------------------------------------------------------
+
+static std::unique_ptr<NS_PROJ::QuadTree::QuadTree<unsigned>>
+BuildQuadTree(const TINShiftFile &file, bool forward) {
+    auto quadtree = std::unique_ptr<NS_PROJ::QuadTree::QuadTree<unsigned>>(
+        new NS_PROJ::QuadTree::QuadTree<unsigned>(GetBounds(file, forward)));
+    const auto &triangles = file.triangles();
+    const auto &vertices = file.vertices();
+    const int idxX = file.transformHorizontalComponent() && !forward ? 2 : 0;
+    const int idxY = file.transformHorizontalComponent() && !forward ? 3 : 1;
+    const unsigned colCount = file.verticesColumnCount();
+    for (size_t i = 0; i < triangles.size(); ++i) {
+        const unsigned i1 = triangles[i].idx1;
+        const unsigned i2 = triangles[i].idx2;
+        const unsigned i3 = triangles[i].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];
+        NS_PROJ::QuadTree::RectObj rect;
+        rect.minx = x1;
+        rect.miny = y1;
+        rect.maxx = x1;
+        rect.maxy = y1;
+        rect.minx = std::min(rect.minx, x2);
+        rect.miny = std::min(rect.miny, y2);
+        rect.maxx = std::max(rect.maxx, x2);
+        rect.maxy = std::max(rect.maxy, y2);
+        rect.minx = std::min(rect.minx, x3);
+        rect.miny = std::min(rect.miny, y3);
+        rect.maxx = std::max(rect.maxx, x3);
+        rect.maxy = std::max(rect.maxy, y3);
+        quadtree->insert(static_cast<unsigned>(i), rect);
+    }
+
+    return quadtree;
+}
+
+// ---------------------------------------------------------------------------
+
+Evaluator::Evaluator(std::unique_ptr<TINShiftFile> &&fileIn)
+    : mFile(std::move(fileIn)) {}
+
+// ---------------------------------------------------------------------------
+
+static const TINShiftFile::VertexIndices *
+FindTriangle(const TINShiftFile &file,
+             const NS_PROJ::QuadTree::QuadTree<unsigned> &quadtree,
+             std::vector<unsigned> &triangleIndices, double x, double y,
+             bool forward, double &lambda1, double &lambda2, double &lambda3) {
+#define USE_QUADTREE
+#ifdef USE_QUADTREE
+    triangleIndices.clear();
+    quadtree.search(x, y, triangleIndices);
+#endif
+    const auto &triangles = file.triangles();
+    const auto &vertices = file.vertices();
+    constexpr double EPS = 1e-10;
+    const int idxX = file.transformHorizontalComponent() && !forward ? 2 : 0;
+    const int idxY = file.transformHorizontalComponent() && !forward ? 3 : 1;
+    const unsigned colCount = file.verticesColumnCount();
+#ifdef USE_QUADTREE
+    for (unsigned i : triangleIndices)
+#else
+    for (size_t i = 0; i < triangles.size(); ++i)
+#endif
+    {
+        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];
+        const double det_T = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3);
+        lambda1 = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) / det_T;
+        lambda2 = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) / det_T;
+        if (lambda1 >= -EPS && lambda1 <= 1 + EPS && lambda2 >= -EPS &&
+            lambda2 <= 1 + EPS) {
+            lambda3 = 1 - lambda1 - lambda2;
+            if (lambda3 >= 0) {
+                return &triangle;
+            }
+        }
+    }
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+bool Evaluator::forward(double x, double y, double z, double &x_out,
+                        double &y_out, double &z_out) {
+    if (!mQuadTreeForward)
+        mQuadTreeForward = BuildQuadTree(*(mFile.get()), true);
+
+    double lambda1 = 0.0;
+    double lambda2 = 0.0;
+    double lambda3 = 0.0;
+    const auto *triangle =
+        FindTriangle(*mFile, *mQuadTreeForward, mTriangleIndices, x, y, true,
+                     lambda1, lambda2, lambda3);
+    if (!triangle)
+        return false;
+    const auto &vertices = mFile->vertices();
+    const unsigned i1 = triangle->idx1;
+    const unsigned i2 = triangle->idx2;
+    const unsigned i3 = triangle->idx3;
+    const unsigned colCount = mFile->verticesColumnCount();
+    if (mFile->transformHorizontalComponent()) {
+        constexpr unsigned idxTargetColX = 2;
+        constexpr unsigned idxTargetColY = 3;
+        x_out = vertices[i1 * colCount + idxTargetColX] * lambda1 +
+                vertices[i2 * colCount + idxTargetColX] * lambda2 +
+                vertices[i3 * colCount + idxTargetColX] * lambda3;
+        y_out = vertices[i1 * colCount + idxTargetColY] * lambda1 +
+                vertices[i2 * colCount + idxTargetColY] * lambda2 +
+                vertices[i3 * colCount + idxTargetColY] * lambda3;
+    } else {
+        x_out = x;
+        y_out = y;
+    }
+    if (mFile->transformVerticalComponent()) {
+        const int idxCol = mFile->transformHorizontalComponent() ? 4 : 2;
+        z_out = z + (vertices[i1 * colCount + idxCol] * lambda1 +
+                     vertices[i2 * colCount + idxCol] * lambda2 +
+                     vertices[i3 * colCount + idxCol] * lambda3);
+    } else {
+        z_out = z;
+    }
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+bool Evaluator::inverse(double x, double y, double z, double &x_out,
+                        double &y_out, double &z_out) {
+    NS_PROJ::QuadTree::QuadTree<unsigned> *quadtree;
+    if (!mFile->transformHorizontalComponent() &&
+        mFile->transformVerticalComponent()) {
+        if (!mQuadTreeForward)
+            mQuadTreeForward = BuildQuadTree(*(mFile.get()), true);
+        quadtree = mQuadTreeForward.get();
+    } else {
+        if (!mQuadTreeInverse)
+            mQuadTreeInverse = BuildQuadTree(*(mFile.get()), false);
+        quadtree = mQuadTreeInverse.get();
+    }
+
+    double lambda1 = 0.0;
+    double lambda2 = 0.0;
+    double lambda3 = 0.0;
+    const auto *triangle = FindTriangle(*mFile, *quadtree, mTriangleIndices, x,
+                                        y, false, lambda1, lambda2, lambda3);
+    if (!triangle)
+        return false;
+    const auto &vertices = mFile->vertices();
+    const unsigned i1 = triangle->idx1;
+    const unsigned i2 = triangle->idx2;
+    const unsigned i3 = triangle->idx3;
+    const unsigned colCount = mFile->verticesColumnCount();
+    if (mFile->transformHorizontalComponent()) {
+        constexpr unsigned idxTargetColX = 0;
+        constexpr unsigned idxTargetColY = 1;
+        x_out = vertices[i1 * colCount + idxTargetColX] * lambda1 +
+                vertices[i2 * colCount + idxTargetColX] * lambda2 +
+                vertices[i3 * colCount + idxTargetColX] * lambda3;
+        y_out = vertices[i1 * colCount + idxTargetColY] * lambda1 +
+                vertices[i2 * colCount + idxTargetColY] * lambda2 +
+                vertices[i3 * colCount + idxTargetColY] * lambda3;
+    } else {
+        x_out = x;
+        y_out = y;
+    }
+    if (mFile->transformVerticalComponent()) {
+        const int idxCol = mFile->transformHorizontalComponent() ? 4 : 2;
+        z_out = z - (vertices[i1 * colCount + idxCol] * lambda1 +
+                     vertices[i2 * colCount + idxCol] * lambda2 +
+                     vertices[i3 * colCount + idxCol] * lambda3);
+    } else {
+        z_out = z;
+    }
+    return true;
+}
+
+} // namespace TINSHIFT_NAMESPACE