diff options
Diffstat (limited to 'src/transformations/defmodel.cpp')
| -rw-r--r-- | src/transformations/defmodel.cpp | 451 |
1 files changed, 451 insertions, 0 deletions
diff --git a/src/transformations/defmodel.cpp b/src/transformations/defmodel.cpp new file mode 100644 index 00000000..3d0f2a58 --- /dev/null +++ b/src/transformations/defmodel.cpp @@ -0,0 +1,451 @@ +/****************************************************************************** + * 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, + "defmodel: 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, + "defmodel: 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, "defmodel: 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, + "defmodel: 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, + "defmodel: 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, + "defmodel: 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, "defmodel: 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, 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, "defmodel: +model= should be specified."); + return destructor(P, PJD_ERR_NO_ARGS); + } + + auto file = NS_PROJ::FileManager::open_resource_file(P->ctx, model); + if (nullptr == file) { + proj_log_error(P, "defmodel: Cannot open %s", model); + 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, "defmodel: File %s too large", model); + 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, "defmodel: Cannot read %s", model); + return destructor(P, PJD_ERR_INVALID_ARG); + } + + 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, "defmodel: invalid model: %s", e.what()); + return destructor(P, PJD_ERR_INVALID_ARG); + } + + 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; +} |