aboutsummaryrefslogtreecommitdiff
path: root/src/transformations/affine.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/transformations/affine.cpp')
-rw-r--r--src/transformations/affine.cpp250
1 files changed, 250 insertions, 0 deletions
diff --git a/src/transformations/affine.cpp b/src/transformations/affine.cpp
new file mode 100644
index 00000000..e2b668d3
--- /dev/null
+++ b/src/transformations/affine.cpp
@@ -0,0 +1,250 @@
+/************************************************************************
+* Copyright (c) 2018, 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__
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj_internal.h"
+#include "proj.h"
+#include "projects.h"
+
+PROJ_HEAD(affine, "Affine transformation");
+PROJ_HEAD(geogoffset, "Geographic Offset");
+
+namespace { // anonymous namespace
+struct pj_affine_coeffs {
+ double s11;
+ double s12;
+ double s13;
+ double s21;
+ double s22;
+ double s23;
+ double s31;
+ double s32;
+ double s33;
+ double tscale;
+};
+} // anonymous namespace
+
+namespace { // anonymous namespace
+struct pj_opaque_affine {
+ double xoff;
+ double yoff;
+ double zoff;
+ double toff;
+ struct pj_affine_coeffs forward;
+ struct pj_affine_coeffs reverse;
+};
+} // anonymous namespace
+
+
+static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
+ PJ_COORD newObs;
+ const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque;
+ const struct pj_affine_coeffs *C = &(Q->forward);
+ newObs.xyzt.x = Q->xoff + C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z;
+ newObs.xyzt.y = Q->yoff + C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z;
+ newObs.xyzt.z = Q->zoff + C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z;
+ newObs.xyzt.t = Q->toff + C->tscale * obs.xyzt.t;
+ return newObs;
+}
+
+static XYZ forward_3d(LPZ lpz, PJ *P) {
+ PJ_COORD point = {{0,0,0,0}};
+ point.lpz = lpz;
+ return forward_4d(point, P).xyz;
+}
+
+
+static XY forward_2d(LP lp, PJ *P) {
+ PJ_COORD point = {{0,0,0,0}};
+ point.lp = lp;
+ return forward_4d(point, P).xy;
+}
+
+
+static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
+ PJ_COORD newObs;
+ const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque;
+ const struct pj_affine_coeffs *C = &(Q->reverse);
+ obs.xyzt.x -= Q->xoff;
+ obs.xyzt.y -= Q->yoff;
+ obs.xyzt.z -= Q->zoff;
+ newObs.xyzt.x = C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z;
+ newObs.xyzt.y = C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z;
+ newObs.xyzt.z = C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z;
+ newObs.xyzt.t = C->tscale * (obs.xyzt.t - Q->toff);
+ return newObs;
+}
+
+static LPZ reverse_3d(XYZ xyz, PJ *P) {
+ PJ_COORD point = {{0,0,0,0}};
+ point.xyz = xyz;
+ return reverse_4d(point, P).lpz;
+}
+
+static LP reverse_2d(XY xy, PJ *P) {
+ PJ_COORD point = {{0,0,0,0}};
+ point.xy = xy;
+ return reverse_4d(point, P).lp;
+}
+
+static struct pj_opaque_affine * initQ() {
+ struct pj_opaque_affine *Q = static_cast<struct pj_opaque_affine *>(pj_calloc(1, sizeof(struct pj_opaque_affine)));
+ if (nullptr==Q)
+ return nullptr;
+
+ /* default values */
+ Q->forward.s11 = 1.0;
+ Q->forward.s22 = 1.0;
+ Q->forward.s33 = 1.0;
+ Q->forward.tscale = 1.0;
+
+ Q->reverse.s11 = 1.0;
+ Q->reverse.s22 = 1.0;
+ Q->reverse.s33 = 1.0;
+ Q->reverse.tscale = 1.0;
+
+ return Q;
+}
+
+static void computeReverseParameters(PJ* P)
+{
+ struct pj_opaque_affine *Q = (struct pj_opaque_affine *) P->opaque;
+
+ /* cf https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3_%C3%97_3_matrices */
+ const double a = Q->forward.s11;
+ const double b = Q->forward.s12;
+ const double c = Q->forward.s13;
+ const double d = Q->forward.s21;
+ const double e = Q->forward.s22;
+ const double f = Q->forward.s23;
+ const double g = Q->forward.s31;
+ const double h = Q->forward.s32;
+ const double i = Q->forward.s33;
+ const double A = e * i - f * h;
+ const double B = -(d * i - f * g);
+ const double C = (d * h - e * g);
+ const double D = -(b * i - c * h);
+ const double E = (a * i - c * g);
+ const double F = -(a * h - b * g);
+ const double G = b * f - c * e;
+ const double H = -(a * f - c * d);
+ const double I = a * e - b * d;
+ const double det = a * A + b * B + c * C;
+ if( det == 0.0 || Q->forward.tscale == 0.0 ) {
+ if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
+ proj_log_debug(P, "Affine: matrix non invertible");
+ }
+ P->inv4d = nullptr;
+ P->inv3d = nullptr;
+ P->inv = nullptr;
+ } else {
+ Q->reverse.s11 = A / det;
+ Q->reverse.s12 = D / det;
+ Q->reverse.s13 = G / det;
+ Q->reverse.s21 = B / det;
+ Q->reverse.s22 = E / det;
+ Q->reverse.s23 = H / det;
+ Q->reverse.s31 = C / det;
+ Q->reverse.s32 = F / det;
+ Q->reverse.s33 = I / det;
+ Q->reverse.tscale = 1.0 / Q->forward.tscale;
+ }
+}
+
+PJ *TRANSFORMATION(affine,0 /* no need for ellipsoid */) {
+ struct pj_opaque_affine *Q = initQ();
+ if (nullptr==Q)
+ return pj_default_destructor(P, ENOMEM);
+ P->opaque = (void *) Q;
+
+ P->fwd4d = forward_4d;
+ P->inv4d = reverse_4d;
+ P->fwd3d = forward_3d;
+ P->inv3d = reverse_3d;
+ P->fwd = forward_2d;
+ P->inv = reverse_2d;
+
+ P->left = PJ_IO_UNITS_WHATEVER;
+ P->right = PJ_IO_UNITS_WHATEVER;
+
+ /* read args */
+ Q->xoff = pj_param(P->ctx, P->params, "dxoff").f;
+ Q->yoff = pj_param(P->ctx, P->params, "dyoff").f;
+ Q->zoff = pj_param(P->ctx, P->params, "dzoff").f;
+ Q->toff = pj_param(P->ctx, P->params, "dtoff").f;
+
+ if(pj_param (P->ctx, P->params, "ts11").i) {
+ Q->forward.s11 = pj_param(P->ctx, P->params, "ds11").f;
+ }
+ Q->forward.s12 = pj_param(P->ctx, P->params, "ds12").f;
+ Q->forward.s13 = pj_param(P->ctx, P->params, "ds13").f;
+ Q->forward.s21 = pj_param(P->ctx, P->params, "ds21").f;
+ if(pj_param (P->ctx, P->params, "ts22").i) {
+ Q->forward.s22 = pj_param(P->ctx, P->params, "ds22").f;
+ }
+ Q->forward.s23 = pj_param(P->ctx, P->params, "ds23").f;
+ Q->forward.s31 = pj_param(P->ctx, P->params, "ds31").f;
+ Q->forward.s32 = pj_param(P->ctx, P->params, "ds32").f;
+ if(pj_param (P->ctx, P->params, "ts33").i) {
+ Q->forward.s33 = pj_param(P->ctx, P->params, "ds33").f;
+ }
+ if(pj_param (P->ctx, P->params, "ttscale").i) {
+ Q->forward.tscale = pj_param(P->ctx, P->params, "dtscale").f;
+ }
+
+ computeReverseParameters(P);
+
+ return P;
+}
+
+
+/* Arcsecond to radians */
+#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0)
+
+
+PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) {
+ struct pj_opaque_affine *Q = initQ();
+ if (nullptr==Q)
+ return pj_default_destructor(P, ENOMEM);
+ P->opaque = (void *) Q;
+
+ P->fwd4d = forward_4d;
+ P->inv4d = reverse_4d;
+ P->fwd3d = forward_3d;
+ P->inv3d = reverse_3d;
+ P->fwd = forward_2d;
+ P->inv = reverse_2d;
+
+ P->left = PJ_IO_UNITS_ANGULAR;
+ P->right = PJ_IO_UNITS_ANGULAR;
+
+ /* read args */
+ Q->xoff = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD;
+ Q->yoff = pj_param(P->ctx, P->params, "ddlat").f * ARCSEC_TO_RAD;
+ Q->zoff = pj_param(P->ctx, P->params, "ddh").f;
+
+ return P;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 13:01:05 +0000