path: root/src/transformations/defmodel.cpp

/******************************************************************************
 * Project:  PROJ
 * Purpose:  Functionality related to deformation model
 * 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 "defmodel.hpp"
#include "filemanager.hpp"
#include "grids.hpp"
#include "proj_internal.h"

#include <assert.h>

#include <map>
#include <memory>
#include <utility>

PROJ_HEAD(defmodel, "Deformation model");

using namespace DeformationModel;

namespace {

struct Grid : public GridPrototype {
    PJ_CONTEXT *ctx;
    const NS_PROJ::GenericShiftGrid *realGrid;
    mutable bool checkedHorizontal = false;
    mutable bool checkedVertical = false;
    mutable int sampleX = 0;
    mutable int sampleY = 1;
    mutable int sampleZ = 2;

    Grid(PJ_CONTEXT *ctxIn, const NS_PROJ::GenericShiftGrid *realGridIn)
        : ctx(ctxIn), realGrid(realGridIn) {
        minx = realGridIn->extentAndRes().west;
        miny = realGridIn->extentAndRes().south;
        resx = realGridIn->extentAndRes().resX;
        resy = realGridIn->extentAndRes().resY;
        width = realGridIn->width();
        height = realGridIn->height();
    }

    bool checkHorizontal(const std::string &expectedUnit) const {
        if (!checkedHorizontal) {
            const auto samplesPerPixel = realGrid->samplesPerPixel();
            if (samplesPerPixel < 2) {
                pj_log(ctx, PJ_LOG_ERROR, "grid %s has not enough samples",
                       realGrid->name().c_str());
                return false;
            }
            bool foundDescX = false;
            bool foundDescY = false;
            bool foundDesc = false;
            for (int i = 0; i < samplesPerPixel; i++) {
                const auto desc = realGrid->description(i);
                if (desc == "east_offset") {
                    sampleX = i;
                    foundDescX = true;
                } else if (desc == "north_offset") {
                    sampleY = i;
                    foundDescY = true;
                }
                if (!desc.empty()) {
                    foundDesc = true;
                }
            }
            if (foundDesc && (!foundDescX || !foundDescY)) {
                pj_log(ctx, PJ_LOG_ERROR, "grid %s : Found band description, "
                                          "but not the ones expected",
                       realGrid->name().c_str());
                return false;
            }
            const auto unit = realGrid->unit(sampleX);
            if (!unit.empty() && unit != expectedUnit) {
                pj_log(ctx, PJ_LOG_ERROR, "grid %s : Only unit=%s "
                                          "currently handled for this mode",
                       realGrid->name().c_str(), expectedUnit.c_str());
                return false;
            }
            checkedHorizontal = true;
        }
        return true;
    }

    bool getLonLatOffset(int ix, int iy, double &lonOffsetRadian,
                         double &latOffsetRadian) const {
        if (!checkHorizontal(STR_DEGREE)) {
            return false;
        }
        float lonOffsetDeg;
        float latOffsetDeg;
        if (!realGrid->valueAt(ix, iy, sampleX, lonOffsetDeg) ||
            !realGrid->valueAt(ix, iy, sampleY, latOffsetDeg)) {
            return false;
        }
        lonOffsetRadian = lonOffsetDeg * DEG_TO_RAD;
        latOffsetRadian = latOffsetDeg * DEG_TO_RAD;
        return true;
    }

    bool getZOffset(int ix, int iy, double &zOffset) const {
        if (!checkedVertical) {
            const auto samplesPerPixel = realGrid->samplesPerPixel();
            if (samplesPerPixel == 1) {
                sampleZ = 0;
            } else if (samplesPerPixel < 3) {
                pj_log(ctx, PJ_LOG_ERROR, "grid %s has not enough samples",
                       realGrid->name().c_str());
                return false;
            }
            bool foundDesc = false;
            bool foundDescZ = false;
            for (int i = 0; i < samplesPerPixel; i++) {
                const auto desc = realGrid->description(i);
                if (desc == "vertical_offset") {
                    sampleZ = i;
                    foundDescZ = true;
                }
                if (!desc.empty()) {
                    foundDesc = true;
                }
            }
            if (foundDesc && !foundDescZ) {
                pj_log(ctx, PJ_LOG_ERROR, "grid %s : Found band description, "
                                          "but not the ones expected",
                       realGrid->name().c_str());
                return false;
            }
            const auto unit = realGrid->unit(sampleZ);
            if (!unit.empty() && unit != STR_METRE) {
                pj_log(ctx, PJ_LOG_ERROR, "grid %s : Only unit=metre currently "
                                          "handled for this mode",
                       realGrid->name().c_str());
                return false;
            }
            checkedVertical = true;
        }
        float zOffsetFloat = 0.0f;
        const bool ret = realGrid->valueAt(ix, iy, sampleZ, zOffsetFloat);
        zOffset = zOffsetFloat;
        return ret;
    }

    bool getEastingNorthingOffset(int ix, int iy, double &eastingOffset,
                                  double &northingOffset) const {
        if (!checkHorizontal(STR_METRE)) {
            return false;
        }
        float eastingOffsetFloat = 0.0f;
        float northingOffsetFloat = 0.0f;
        const bool ret =
            realGrid->valueAt(ix, iy, sampleX, eastingOffsetFloat) &&
            realGrid->valueAt(ix, iy, sampleY, northingOffsetFloat);
        eastingOffset = eastingOffsetFloat;
        northingOffset = northingOffsetFloat;
        return ret;
    }

    bool getLonLatZOffset(int ix, int iy, double &lonOffsetRadian,
                          double &latOffsetRadian, double &zOffset) const {
        return getLonLatOffset(ix, iy, lonOffsetRadian, latOffsetRadian) &&
               getZOffset(ix, iy, zOffset);
    }

    bool getEastingNorthingZOffset(int ix, int iy, double &eastingOffset,
                                   double &northingOffset,
                                   double &zOffset) const {
        return getEastingNorthingOffset(ix, iy, eastingOffset,
                                        northingOffset) &&
               getZOffset(ix, iy, zOffset);
    }

#ifdef DEBUG_DEFMODEL
    std::string name() const { return realGrid->name(); }
#endif

  private:
    Grid(const Grid &) = delete;
    Grid &operator=(const Grid &) = delete;
};

struct GridSet : public GridSetPrototype<Grid> {

    PJ_CONTEXT *ctx;
    std::unique_ptr<NS_PROJ::GenericShiftGridSet> realGridSet;
    std::map<const NS_PROJ::GenericShiftGrid *, std::unique_ptr<Grid>>
        mapGrids{};

    GridSet(PJ_CONTEXT *ctxIn,
            std::unique_ptr<NS_PROJ::GenericShiftGridSet> &&realGridSetIn)
        : ctx(ctxIn), realGridSet(std::move(realGridSetIn)) {}

    const Grid *gridAt(double x, double y) {
        const NS_PROJ::GenericShiftGrid *realGrid = realGridSet->gridAt(x, y);
        if (!realGrid) {
            return nullptr;
        }
        auto iter = mapGrids.find(realGrid);
        if (iter == mapGrids.end()) {
            iter = mapGrids
                       .insert(std::pair<const NS_PROJ::GenericShiftGrid *,
                                         std::unique_ptr<Grid>>(
                           realGrid,
                           std::unique_ptr<Grid>(new Grid(ctx, realGrid))))
                       .first;
        }
        return iter->second.get();
    }

  private:
    GridSet(const GridSet &) = delete;
    GridSet &operator=(const GridSet &) = delete;
};

struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> {

    PJ_CONTEXT *ctx;
    PJ *cart;

    EvaluatorIface(PJ_CONTEXT *ctxIn, PJ *cartIn) : ctx(ctxIn), cart(cartIn) {}

    ~EvaluatorIface() {
        if (cart)
            cart->destructor(cart, 0);
    }

    std::unique_ptr<GridSet> open(const std::string &filename) {
        auto realGridSet = NS_PROJ::GenericShiftGridSet::open(ctx, filename);
        if (!realGridSet) {
            pj_log(ctx, PJ_LOG_ERROR, "cannot open %s", filename.c_str());
            return nullptr;
        }
        return std::unique_ptr<GridSet>(
            new GridSet(ctx, std::move(realGridSet)));
    }

    bool isGeographicCRS(const std::string &crsDef) {
        PJ *P = proj_create(ctx, crsDef.c_str());
        if (P == nullptr) {
            return true; // reasonable default value
        }
        const auto type = proj_get_type(P);
        bool ret = (type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
                    type == PJ_TYPE_GEOGRAPHIC_3D_CRS);
        proj_destroy(P);
        return ret;
    }

    void geographicToGeocentric(double lam, double phi, double height, double a,
                                double b, double /*es*/, double &X, double &Y,
                                double &Z) {
        (void)a;
        (void)b;
        assert(cart->a == a);
        assert(cart->b == b);
        PJ_LPZ lpz;
        lpz.lam = lam;
        lpz.phi = phi;
        lpz.z = height;
        PJ_XYZ xyz = cart->fwd3d(lpz, cart);
        X = xyz.x;
        Y = xyz.y;
        Z = xyz.z;
    }

    void geocentricToGeographic(double X, double Y, double Z, double a,
                                double b, double /*es*/, double &lam,
                                double &phi, double &height) {
        (void)a;
        (void)b;
        assert(cart->a == a);
        assert(cart->b == b);
        PJ_XYZ xyz;
        xyz.x = X;
        xyz.y = Y;
        xyz.z = Z;
        PJ_LPZ lpz = cart->inv3d(xyz, cart);
        lam = lpz.lam;
        phi = lpz.phi;
        height = lpz.z;
    }

#ifdef DEBUG_DEFMODEL
    void log(const std::string &msg) {
        pj_log(ctx, PJ_LOG_TRACE, "%s", msg.c_str());
    }
#endif

  private:
    EvaluatorIface(const EvaluatorIface &) = delete;
    EvaluatorIface &operator=(const EvaluatorIface &) = delete;
};

struct defmodelData {
    std::unique_ptr<Evaluator<Grid, GridSet, EvaluatorIface>> evaluator{};
    EvaluatorIface evaluatorIface;

    explicit defmodelData(PJ_CONTEXT *ctx, PJ *cart)
        : evaluatorIface(ctx, cart) {}

    defmodelData(const defmodelData &) = delete;
    defmodelData &operator=(const defmodelData &) = delete;
};

} // namespace

static PJ *destructor(PJ *P, int errlev) {
    if (nullptr == P)
        return nullptr;

    auto Q = static_cast<struct defmodelData *>(P->opaque);
    delete Q;
    P->opaque = nullptr;

    return pj_default_destructor(P, errlev);
}

static PJ_COORD forward_4d(PJ_COORD in, PJ *P) {
    auto *Q = (struct defmodelData *)P->opaque;

    PJ_COORD out;
    out.xyzt.t = in.xyzt.t;

    if (!Q->evaluator->forward(Q->evaluatorIface, in.xyzt.x, in.xyzt.y,
                               in.xyzt.z, in.xyzt.t, out.xyzt.x, out.xyzt.y,
                               out.xyzt.z)) {
        return proj_coord_error();
    }

    return out;
}

static PJ_COORD reverse_4d(PJ_COORD in, PJ *P) {
    auto *Q = (struct defmodelData *)P->opaque;

    PJ_COORD out;
    out.xyzt.t = in.xyzt.t;

    if (!Q->evaluator->inverse(Q->evaluatorIface, in.xyzt.x, in.xyzt.y,
                               in.xyzt.z, in.xyzt.t, out.xyzt.x, out.xyzt.y,
                               out.xyzt.z)) {
        return proj_coord_error();
    }

    return out;
}

// Function called by proj_assign_context() when a new context is assigned to
// an existing PJ object. Mostly to deal with objects being passed between
// threads.
static void reassign_context(PJ *P, PJ_CONTEXT *ctx) {
    auto *Q = (struct defmodelData *)P->opaque;
    if (Q->evaluatorIface.ctx != ctx) {
        Q->evaluator->clearGridCache();
        Q->evaluatorIface.ctx = ctx;
    }
}

PJ *TRANSFORMATION(defmodel, 1) {
    // Pass a dummy ellipsoid definition that will be overridden just afterwards
    auto cart = proj_create(P->ctx, "+proj=cart +a=1");
    if (cart == nullptr)
        return destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);

    /* inherit ellipsoid definition from P to Q->cart */
    pj_inherit_ellipsoid_def(P, cart);

    auto Q = new defmodelData(P->ctx, cart);
    P->opaque = (void *)Q;
    P->destructor = destructor;
    P->reassign_context = reassign_context;

    const char *model = pj_param(P->ctx, P->params, "smodel").s;
    if (!model) {
        proj_log_error(P, _("+model= should be specified."));
        return destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }

    auto file = NS_PROJ::FileManager::open_resource_file(P->ctx, model);
    if (nullptr == file) {
        proj_log_error(P, _("Cannot open %s"), model);
        return destructor(P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID);
    }
    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, _("File %s too large"), model);
        return destructor(P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID);
    }
    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, _("Cannot read %s"), model);
        return destructor(P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID);
    }

    try {
        Q->evaluator.reset(new Evaluator<Grid, GridSet, EvaluatorIface>(
            MasterFile::parse(jsonStr), Q->evaluatorIface, P->a, P->b));
    } catch (const std::exception &e) {
        proj_log_error(P, _("invalid model: %s"), e.what());
        return destructor(P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID);
    }

    P->fwd4d = forward_4d;
    P->inv4d = reverse_4d;

    if (Q->evaluator->isGeographicCRS()) {
        P->left = PJ_IO_UNITS_RADIANS;
        P->right = PJ_IO_UNITS_RADIANS;
    } else {
        P->left = PJ_IO_UNITS_PROJECTED;
        P->right = PJ_IO_UNITS_PROJECTED;
    }

    return P;
}