diff options
| -rw-r--r-- | docs/source/operations/transformations/affine.rst | 123 | ||||
| -rw-r--r-- | docs/source/operations/transformations/geogoffset.rst | 2 | ||||
| -rw-r--r-- | docs/source/operations/transformations/index.rst | 1 | ||||
| -rw-r--r-- | src/Makefile.am | 2 | ||||
| -rw-r--r-- | src/PJ_affine.c | 246 | ||||
| -rw-r--r-- | src/PJ_geogoffset.c | 111 | ||||
| -rw-r--r-- | src/lib_proj.cmake | 2 | ||||
| -rw-r--r-- | src/makefile.vc | 2 | ||||
| -rw-r--r-- | src/pj_list.h | 1 | ||||
| -rw-r--r-- | test/gie/more_builtins.gie | 46 |
10 files changed, 421 insertions, 115 deletions
diff --git a/docs/source/operations/transformations/affine.rst b/docs/source/operations/transformations/affine.rst new file mode 100644 index 00000000..8b1681b3 --- /dev/null +++ b/docs/source/operations/transformations/affine.rst @@ -0,0 +1,123 @@ +.. _affine: + +================================================================================ +Affine transformation +================================================================================ + +.. versionadded:: 6.0.0 + +The affine transformation applies translation and scaling/rotation terms on the +x,y,z coordinates, and translation and scaling on the temporal cordinate. + ++---------------------+----------------------------------------------------------+ +| **Alias** | affine | ++---------------------+----------------------------------------------------------+ +| **Domain** | 4D | ++---------------------+----------------------------------------------------------+ +| **Input type** | XYZT | ++---------------------+----------------------------------------------------------+ +| **output type** | XYZT | ++---------------------+----------------------------------------------------------+ + +By default, the parameters are set for an identity transforms. The transformation +is reversible unless the determinant of the sji matrix is 0, or `tscale` is 0 + + +Parameters +################################################################################ + +Optional +------------------------------------------------------------------------------- + +.. option:: +xoff=<value> + + Offset in X. Default value: 0 + +.. option:: +yoff=<value> + + Offset in Y. Default value: 0 + +.. option:: +zoff=<value> + + Offset in Z. Default value: 0 + +.. option:: +toff=<value> + + Offset in T. Default value: 0 + +.. option:: +s11=<value> + + Rotation/scaling term. Default value: 1 + +.. option:: +s12=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s13=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s21=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s22=<value> + + Rotation/scaling term. Default value: 1 + +.. option:: +s23=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s31=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s32=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s33=<value> + + Rotation/scaling term. Default value: 1 + +.. option:: +tscale=<value> + + Time scaling term. Default value: 1 + + + +Mathematical description ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. math:: + :label: formula + + \begin{align} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + T \\ + \end{bmatrix}^{dest} = + \begin{bmatrix} + xoff \\ + yoff \\ + zoff \\ + toff \\ + \end{bmatrix} + + \begin{bmatrix} + s11 & s12 & s13 & 0 \\ + s21 & s22 & s23 & 0 \\ + s31 & s32 & s33 & 0 \\ + 0 & 0 & 0 & tscale \\ + \end{bmatrix} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + T \\ + \end{bmatrix}^{source} + \end{align} + + diff --git a/docs/source/operations/transformations/geogoffset.rst b/docs/source/operations/transformations/geogoffset.rst index 91c91c83..f643485e 100644 --- a/docs/source/operations/transformations/geogoffset.rst +++ b/docs/source/operations/transformations/geogoffset.rst @@ -32,6 +32,8 @@ It can also be used to implement the method Vertical offset (code 9616) The reverse transformation simply consists in subtracting the offsets. +This method is a conveniency wrapper for the more general :ref:`affine`. + Examples ############################################################################### diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst index c2fe6baa..121215f7 100644 --- a/docs/source/operations/transformations/index.rst +++ b/docs/source/operations/transformations/index.rst @@ -10,6 +10,7 @@ systems are based on different datums. .. toctree:: :maxdepth: 1 + affine deformation geogoffset helmert diff --git a/src/Makefile.am b/src/Makefile.am index 2320e812..c1fbc12e 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -89,7 +89,7 @@ libproj_la_SOURCES = \ \ proj_4D_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \ PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c PJ_molodensky.c \ - PJ_deformation.c pj_internal.c PJ_axisswap.c PJ_geogoffset.c + PJ_deformation.c pj_internal.c PJ_axisswap.c PJ_affine.c install-exec-local: rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT) diff --git a/src/PJ_affine.c b/src/PJ_affine.c new file mode 100644 index 00000000..0d8b6374 --- /dev/null +++ b/src/PJ_affine.c @@ -0,0 +1,246 @@ +/************************************************************************ +* 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"); + +struct pj_affine_coeffs { + double s11; + double s12; + double s13; + double s21; + double s22; + double s23; + double s31; + double s32; + double s33; + double tscale; +}; + +struct pj_opaque_affine { + double xoff; + double yoff; + double zoff; + double toff; + struct pj_affine_coeffs forward; + struct pj_affine_coeffs reverse; +}; + + +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 = pj_calloc(1, sizeof(struct pj_opaque_affine)); + if (0==Q) + return 0; + + /* 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 = NULL; + P->inv3d = NULL; + P->inv = NULL; + } 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 (0==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 (0==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; +} diff --git a/src/PJ_geogoffset.c b/src/PJ_geogoffset.c deleted file mode 100644 index 7f917fd7..00000000 --- a/src/PJ_geogoffset.c +++ /dev/null @@ -1,111 +0,0 @@ -/************************************************************************ -* 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.h" -#include "projects.h" - -PROJ_HEAD(geogoffset, "Geographic Offset"); - -struct pj_opaque_geogoffset { - double dlon; - double dlat; - double dh; -}; - - -static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { - struct pj_opaque_geogoffset *Q = (struct pj_opaque_geogoffset *) P->opaque; - /* offset coordinate */ - obs.lpz.lam += Q->dlon; - obs.lpz.phi += Q->dlat; - obs.lpz.z += Q->dh; - return obs; -} - -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) { - struct pj_opaque_geogoffset *Q = (struct pj_opaque_geogoffset *) P->opaque; - /* offset coordinate */ - obs.lpz.lam -= Q->dlon; - obs.lpz.phi -= Q->dlat; - obs.lpz.z -= Q->dh; - return obs; -} - -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; -} - - -/* Arcsecond to radians */ -#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0) - - -PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) { - struct pj_opaque_geogoffset *Q = pj_calloc(1, sizeof(struct pj_opaque_geogoffset)); - if (0==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->dlon = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD; - Q->dlat = pj_param(P->ctx, P->params, "ddlat").f * ARCSEC_TO_RAD; - Q->dh = pj_param(P->ctx, P->params, "ddh").f; - - return P; -} diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index f1e76ff3..a359369e 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -40,6 +40,7 @@ SET(SRC_LIBPROJ_PJ nad_init.c PJ_aea.c PJ_aeqd.c + PJ_affine.c PJ_airy.c PJ_aitoff.c PJ_august.c @@ -72,7 +73,6 @@ SET(SRC_LIBPROJ_PJ PJ_fouc_s.c PJ_gall.c PJ_geoc.c - PJ_geogoffset.c PJ_geos.c PJ_gins8.c PJ_gnom.c diff --git a/src/makefile.vc b/src/makefile.vc index 21db1362..73c68c0f 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -61,7 +61,7 @@ support = \ pipeline = \ proj_4D_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \ PJ_vgridshift.obj PJ_hgridshift.obj PJ_unitconvert.obj PJ_molodensky.obj \ - PJ_deformation.obj PJ_axisswap.obj PJ_geogoffset.obj + PJ_deformation.obj PJ_axisswap.obj PJ_affine.obj geodesic = geodesic.obj diff --git a/src/pj_list.h b/src/pj_list.h index 6d94f4fd..72a9af7f 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -7,6 +7,7 @@ static const char PJ_LIST_H_ID[] = "@(#)pj_list.h 4.5 95/08/09 GIE REL"; */ PROJ_HEAD(aea, "Albers Equal Area") PROJ_HEAD(aeqd, "Azimuthal Equidistant") +PROJ_HEAD(affine, "Affine transformation") PROJ_HEAD(airy, "Airy") PROJ_HEAD(aitoff, "Aitoff") PROJ_HEAD(alsk, "Mod. Stererographics of Alaska") diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie index d3fe64f6..10feea70 100644 --- a/test/gie/more_builtins.gie +++ b/test/gie/more_builtins.gie @@ -559,7 +559,7 @@ expect 0 -90 ------------------------------------------------------------------------------- -Test for PJ_geogoffset +Test for PJ_affine ------------------------------------------------------------------------------- ------------------------------------------------------------------------------- operation +proj=geogoffset @@ -590,6 +590,50 @@ expect 11 19 33 40 roundtrip 1 ------------------------------------------------------------------------------- +operation +proj=affine +------------------------------------------------------------------------------- +direction forward +tolerance 1mm + +accept 10 20 30 40 +expect 10 20 30 40 +roundtrip 1 + +------------------------------------------------------------------------------- +operation +proj=affine +xoff=1 +yoff=2 +zoff=3 +toff=4 +s11=11 +s12=12 +s13=13 +s21=21 +s22=22 +s23=23 +s31=-31 +s32=32 +s33=33 +tscale=34 +------------------------------------------------------------------------------- +direction forward +tolerance 1mm + +accept 2 49 10 100 +expect 741.0000 1352.0000 1839.0000 3404.0000 +roundtrip 1 + +accept 2 49 10 +expect 741.0000 1352.0000 1839.0000 +roundtrip 1 + +accept 2 49 +expect 611.0000 1122.0000 +roundtrip 1 + +------------------------------------------------------------------------------- +# Non invertible +operation +proj=affine +s11=0 +s22=0 +s23=0 +------------------------------------------------------------------------------- +direction reverse +accept 0 0 0 0 +expect failure + +------------------------------------------------------------------------------- +# Non invertible +operation +proj=affine +tscale=0 +------------------------------------------------------------------------------- +direction reverse +accept 0 0 0 0 +expect failure + +------------------------------------------------------------------------------- run the few gie-builtin tests, which are currently either awkward or impossible to express in the gie command set ------------------------------------------------------------------------------- |