From eafeb61ce59aeb34dabf38f55f70ba9a3b779c4b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 6 Jan 2020 23:03:29 +0100 Subject: Grid refactoring: address review comments of https://github.com/OSGeo/PROJ/pull/1777 - Move content of legacy apply_gridshift.cpp and apply_vgridshift.cpp in grids.cpp - Rename nad_ functions to pj_hgrid_ - Rename internal proj_hgrid_/proj_vgrid_ functions to pj_ --- src/transformations/deformation.cpp | 10 +++++----- src/transformations/hgridshift.cpp | 10 +++++----- src/transformations/vgridshift.cpp | 10 +++++----- src/transformations/xyzgridshift.cpp | 4 ++-- 4 files changed, 17 insertions(+), 17 deletions(-) (limited to 'src/transformations') diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp index eb109826..993647fc 100644 --- a/src/transformations/deformation.cpp +++ b/src/transformations/deformation.cpp @@ -226,8 +226,8 @@ static PJ_XYZ get_grid_shift(PJ* P, const PJ_XYZ& cartesian) { } else { - shift.lp = proj_hgrid_value(P, Q->hgrids, geodetic.lp); - shift.enu.u = proj_vgrid_value(P, Q->vgrids, geodetic.lp, 1.0); + shift.lp = pj_hgrid_value(P, Q->hgrids, geodetic.lp); + shift.enu.u = pj_vgrid_value(P, Q->vgrids, geodetic.lp, 1.0); if (proj_errno(P) == PJD_ERR_GRID_AREA) proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model", @@ -425,7 +425,7 @@ PJ *TRANSFORMATION(deformation,1) { if( has_grids ) { - Q->grids = proj_generic_grid_init(P, "grids"); + Q->grids = pj_generic_grid_init(P, "grids"); /* Was gridlist compiled properly? */ if ( proj_errno(P) ) { proj_log_error(P, "deformation: could not find required grid(s)."); @@ -434,13 +434,13 @@ PJ *TRANSFORMATION(deformation,1) { } else { - Q->hgrids = proj_hgrid_init(P, "xy_grids"); + Q->hgrids = pj_hgrid_init(P, "xy_grids"); if (proj_errno(P)) { proj_log_error(P, "deformation: could not find requested xy_grid(s)."); return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); } - Q->vgrids = proj_vgrid_init(P, "z_grids"); + Q->vgrids = pj_vgrid_init(P, "z_grids"); if (proj_errno(P)) { proj_log_error(P, "deformation: could not find requested z_grid(s)."); return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp index 24da4dde..122a7ab2 100644 --- a/src/transformations/hgridshift.cpp +++ b/src/transformations/hgridshift.cpp @@ -28,7 +28,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_hgrid_init(P, "grids"); + Q->grids = pj_hgrid_init(P, "grids"); if ( proj_errno(P) ) { return proj_coord_error().xyz; } @@ -37,7 +37,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { if (!Q->grids.empty()) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - point.lp = proj_hgrid_apply(P, Q->grids, point.lp, PJ_FWD); + point.lp = pj_hgrid_apply(P->ctx, Q->grids, point.lp, PJ_FWD); } return point.xyz; @@ -51,7 +51,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_hgrid_init(P, "grids"); + Q->grids = pj_hgrid_init(P, "grids"); if ( proj_errno(P) ) { return proj_coord_error().lpz; } @@ -60,7 +60,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { if (!Q->grids.empty()) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - point.lp = proj_hgrid_apply(P, Q->grids, point.lp, PJ_INV); + point.lp = pj_hgrid_apply(P->ctx, Q->grids, point.lp, PJ_INV); } return point.lpz; @@ -165,7 +165,7 @@ PJ *TRANSFORMATION(hgridshift,0) { Q->defer_grid_opening = true; } else { - Q->grids = proj_hgrid_init(P, "grids"); + Q->grids = pj_hgrid_init(P, "grids"); /* Was gridlist compiled properly? */ if ( proj_errno(P) ) { proj_log_error(P, "hgridshift: could not find required grid(s)."); diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp index 3e7a015e..121b795a 100644 --- a/src/transformations/vgridshift.cpp +++ b/src/transformations/vgridshift.cpp @@ -56,7 +56,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_vgrid_init(P, "grids"); + Q->grids = pj_vgrid_init(P, "grids"); deal_with_vertcon_gtx_hack(P); if ( proj_errno(P) ) { return proj_coord_error().xyz; @@ -66,7 +66,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { if (!Q->grids.empty()) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - point.xyz.z += proj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier); + point.xyz.z += pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier); } return point.xyz; @@ -80,7 +80,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_vgrid_init(P, "grids"); + Q->grids = pj_vgrid_init(P, "grids"); deal_with_vertcon_gtx_hack(P); if ( proj_errno(P) ) { return proj_coord_error().lpz; @@ -90,7 +90,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { if (!Q->grids.empty()) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - point.xyz.z -= proj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier); + point.xyz.z -= pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier); } return point.lpz; @@ -193,7 +193,7 @@ PJ *TRANSFORMATION(vgridshift,0) { } else { /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */ - Q->grids = proj_vgrid_init(P, "grids"); + Q->grids = pj_vgrid_init(P, "grids"); /* Was gridlist compiled properly? */ if ( proj_errno(P) ) { diff --git a/src/transformations/xyzgridshift.cpp b/src/transformations/xyzgridshift.cpp index a76f3255..3ec3863c 100644 --- a/src/transformations/xyzgridshift.cpp +++ b/src/transformations/xyzgridshift.cpp @@ -77,7 +77,7 @@ static bool get_grid_values(PJ* P, { if ( Q->defer_grid_opening ) { Q->defer_grid_opening = false; - Q->grids = proj_generic_grid_init(P, "grids"); + Q->grids = pj_generic_grid_init(P, "grids"); if ( proj_errno(P) ) { return false; } @@ -341,7 +341,7 @@ PJ *TRANSFORMATION(xyzgridshift,0) { Q->defer_grid_opening = true; } else { - Q->grids = proj_generic_grid_init(P, "grids"); + Q->grids = pj_generic_grid_init(P, "grids"); /* Was gridlist compiled properly? */ if ( proj_errno(P) ) { proj_log_error(P, "xyzgridshift: could not find required grid(s)."); -- cgit v1.2.3 From 237296b7e84a8bb270e3be06a690737a601d73e7 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 7 Jan 2020 03:37:23 +0100 Subject: Remote grid: add mechanism to re-open a grid if it has changed while being opened --- src/transformations/deformation.cpp | 48 +++++++++++++++++++----------------- src/transformations/xyzgridshift.cpp | 48 +++++++++++++++++++----------------- 2 files changed, 50 insertions(+), 46 deletions(-) (limited to 'src/transformations') diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp index 993647fc..778c1a2e 100644 --- a/src/transformations/deformation.cpp +++ b/src/transformations/deformation.cpp @@ -81,13 +81,16 @@ struct deformationData { // --------------------------------------------------------------------------- static const GenericShiftGrid* findGrid(const ListOfGenericGrids& grids, - const PJ_LP& input) + const PJ_LP& input, + GenericShiftGridSet *&gridSetOut) { for( const auto& gridset: grids ) { auto grid = gridset->gridAt(input.lam, input.phi); - if( grid ) + if( grid ) { + gridSetOut = gridset.get(); return grid; + } } return nullptr; } @@ -101,7 +104,8 @@ static bool get_grid_values(PJ* P, double& vy, double& vz) { - auto grid = findGrid(Q->grids, lp); + GenericShiftGridSet* gridset = nullptr; + auto grid = findGrid(Q->grids, lp, gridset); if( !grid ) { return false; } @@ -145,30 +149,28 @@ static bool get_grid_values(PJ* P, int iy2 = std::min(iy + 1, grid->height() - 1); float dx1, dy1, dz1; - if( !grid->valueAt(ix, iy, sampleE, dx1) || - !grid->valueAt(ix, iy, sampleN, dy1) || - !grid->valueAt(ix, iy, sampleU, dz1) ) { - return false; - } - float dx2, dy2, dz2; - if( !grid->valueAt(ix2, iy, sampleE, dx2) || - !grid->valueAt(ix2, iy, sampleN, dy2) || - !grid->valueAt(ix2, iy, sampleU, dz2) ) { - return false; - } - float dx3, dy3, dz3; - if( !grid->valueAt(ix, iy2, sampleE, dx3) || - !grid->valueAt(ix, iy2, sampleN, dy3) || - !grid->valueAt(ix, iy2, sampleU, dz3) ) { - return false; - } - float dx4, dy4, dz4; - if( !grid->valueAt(ix2, iy2, sampleE, dx4) || + bool error =( !grid->valueAt(ix, iy, sampleE, dx1) || + !grid->valueAt(ix, iy, sampleN, dy1) || + !grid->valueAt(ix, iy, sampleU, dz1) || + !grid->valueAt(ix2, iy, sampleE, dx2) || + !grid->valueAt(ix2, iy, sampleN, dy2) || + !grid->valueAt(ix2, iy, sampleU, dz2) || + !grid->valueAt(ix, iy2, sampleE, dx3) || + !grid->valueAt(ix, iy2, sampleN, dy3) || + !grid->valueAt(ix, iy2, sampleU, dz3) || + !grid->valueAt(ix2, iy2, sampleE, dx4) || !grid->valueAt(ix2, iy2, sampleN, dy4) || - !grid->valueAt(ix2, iy2, sampleU, dz4) ) { + !grid->valueAt(ix2, iy2, sampleU, dz4) ); + if( grid->hasChanged() ) { + if( gridset->reopen(P->ctx) ) { + return get_grid_values( P, Q, lp, vx, vy, vz); + } + error = true; + } + if( error ) { return false; } diff --git a/src/transformations/xyzgridshift.cpp b/src/transformations/xyzgridshift.cpp index 3ec3863c..4a85fb4c 100644 --- a/src/transformations/xyzgridshift.cpp +++ b/src/transformations/xyzgridshift.cpp @@ -55,13 +55,16 @@ struct xyzgridshiftData { // --------------------------------------------------------------------------- static const GenericShiftGrid* findGrid(const ListOfGenericGrids& grids, - const PJ_LP& input) + const PJ_LP& input, + GenericShiftGridSet *&gridSetOut) { for( const auto& gridset: grids ) { auto grid = gridset->gridAt(input.lam, input.phi); - if( grid ) + if( grid ) { + gridSetOut = gridset.get(); return grid; + } } return nullptr; } @@ -83,7 +86,8 @@ static bool get_grid_values(PJ* P, } } - auto grid = findGrid(Q->grids, lp); + GenericShiftGridSet* gridset = nullptr; + auto grid = findGrid(Q->grids, lp, gridset); if( !grid ) { return false; } @@ -127,30 +131,28 @@ static bool get_grid_values(PJ* P, int iy2 = std::min(iy + 1, grid->height() - 1); float dx1, dy1, dz1; - if( !grid->valueAt(ix, iy, sampleX, dx1) || - !grid->valueAt(ix, iy, sampleY, dy1) || - !grid->valueAt(ix, iy, sampleZ, dz1) ) { - return false; - } - float dx2, dy2, dz2; - if( !grid->valueAt(ix2, iy, sampleX, dx2) || - !grid->valueAt(ix2, iy, sampleY, dy2) || - !grid->valueAt(ix2, iy, sampleZ, dz2) ) { - return false; - } - float dx3, dy3, dz3; - if( !grid->valueAt(ix, iy2, sampleX, dx3) || - !grid->valueAt(ix, iy2, sampleY, dy3) || - !grid->valueAt(ix, iy2, sampleZ, dz3) ) { - return false; - } - float dx4, dy4, dz4; - if( !grid->valueAt(ix2, iy2, sampleX, dx4) || + bool error =( !grid->valueAt(ix, iy, sampleX, dx1) || + !grid->valueAt(ix, iy, sampleY, dy1) || + !grid->valueAt(ix, iy, sampleZ, dz1) || + !grid->valueAt(ix2, iy, sampleX, dx2) || + !grid->valueAt(ix2, iy, sampleY, dy2) || + !grid->valueAt(ix2, iy, sampleZ, dz2) || + !grid->valueAt(ix, iy2, sampleX, dx3) || + !grid->valueAt(ix, iy2, sampleY, dy3) || + !grid->valueAt(ix, iy2, sampleZ, dz3) || + !grid->valueAt(ix2, iy2, sampleX, dx4) || !grid->valueAt(ix2, iy2, sampleY, dy4) || - !grid->valueAt(ix2, iy2, sampleZ, dz4) ) { + !grid->valueAt(ix2, iy2, sampleZ, dz4) ); + if( grid->hasChanged() ) { + if( gridset->reopen(P->ctx) ) { + return get_grid_values( P, Q, lp, dx, dy, dz); + } + error = true; + } + if( error ) { return false; } -- cgit v1.2.3 From 408cb1e4279786e3729e942c2f8b2ee2ad435c94 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 7 Jan 2020 15:01:36 +0100 Subject: deformation.cpp and xyzgridshift.cpp: factor out common code (addresses https://github.com/OSGeo/PROJ/pull/1790#discussion_r363748615) --- src/transformations/deformation.cpp | 74 ++++++------------------------------ src/transformations/xyzgridshift.cpp | 74 ++++++------------------------------ 2 files changed, 22 insertions(+), 126 deletions(-) (limited to 'src/transformations') diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp index 778c1a2e..8aee50c9 100644 --- a/src/transformations/deformation.cpp +++ b/src/transformations/deformation.cpp @@ -80,23 +80,6 @@ struct deformationData { // --------------------------------------------------------------------------- -static const GenericShiftGrid* findGrid(const ListOfGenericGrids& grids, - const PJ_LP& input, - GenericShiftGridSet *&gridSetOut) -{ - for( const auto& gridset: grids ) - { - auto grid = gridset->gridAt(input.lam, input.phi); - if( grid ) { - gridSetOut = gridset.get(); - return grid; - } - } - return nullptr; -} - -// --------------------------------------------------------------------------- - static bool get_grid_values(PJ* P, deformationData* Q, const PJ_LP& lp, @@ -105,7 +88,7 @@ static bool get_grid_values(PJ* P, double& vz) { GenericShiftGridSet* gridset = nullptr; - auto grid = findGrid(Q->grids, lp, gridset); + auto grid = pj_find_generic_grid(Q->grids, lp, gridset); if( !grid ) { return false; } @@ -140,55 +123,20 @@ static bool get_grid_values(PJ* P, return false; } - const auto& extent = grid->extentAndRes(); - double grid_x = (lp.lam - extent.westLon) / extent.resLon; - double grid_y = (lp.phi - extent.southLat) / extent.resLat; - int ix = static_cast(grid_x); - int iy = static_cast(grid_y); - int ix2 = std::min(ix + 1, grid->width() - 1); - int iy2 = std::min(iy + 1, grid->height() - 1); - - float dx1, dy1, dz1; - float dx2, dy2, dz2; - float dx3, dy3, dz3; - float dx4, dy4, dz4; - bool error =( !grid->valueAt(ix, iy, sampleE, dx1) || - !grid->valueAt(ix, iy, sampleN, dy1) || - !grid->valueAt(ix, iy, sampleU, dz1) || - !grid->valueAt(ix2, iy, sampleE, dx2) || - !grid->valueAt(ix2, iy, sampleN, dy2) || - !grid->valueAt(ix2, iy, sampleU, dz2) || - !grid->valueAt(ix, iy2, sampleE, dx3) || - !grid->valueAt(ix, iy2, sampleN, dy3) || - !grid->valueAt(ix, iy2, sampleU, dz3) || - !grid->valueAt(ix2, iy2, sampleE, dx4) || - !grid->valueAt(ix2, iy2, sampleN, dy4) || - !grid->valueAt(ix2, iy2, sampleU, dz4) ); - if( grid->hasChanged() ) { - if( gridset->reopen(P->ctx) ) { + bool must_retry = false; + if( !pj_bilinear_interpolation_three_samples(grid, lp, + sampleE, sampleN, sampleU, + vx, vy, vz, + must_retry) ) + { + if( must_retry ) return get_grid_values( P, Q, lp, vx, vy, vz); - } - error = true; - } - if( error ) { return false; } - - double frct_lam = grid_x - ix; - double frct_phi = grid_y - iy; - double m10 = frct_lam; - double m11 = m10; - double m01 = 1. - frct_lam; - double m00 = m01; - m11 *= frct_phi; - m01 *= frct_phi; - frct_phi = 1. - frct_phi; - m00 *= frct_phi; - m10 *= frct_phi; // divide by 1000 to get m/year - vx = (m00 * dx1 + m10 * dx2 + m01 * dx3 + m11 * dx4) / 1000; - vy = (m00 * dy1 + m10 * dy2 + m01 * dy3 + m11 * dy4) / 1000; - vz = (m00 * dz1 + m10 * dz2 + m01 * dz3 + m11 * dz4) / 1000; + vx /= 1000; + vy /= 1000; + vz /= 1000; return true; } diff --git a/src/transformations/xyzgridshift.cpp b/src/transformations/xyzgridshift.cpp index 4a85fb4c..e1c76518 100644 --- a/src/transformations/xyzgridshift.cpp +++ b/src/transformations/xyzgridshift.cpp @@ -51,24 +51,6 @@ struct xyzgridshiftData { }; } // anonymous namespace - -// --------------------------------------------------------------------------- - -static const GenericShiftGrid* findGrid(const ListOfGenericGrids& grids, - const PJ_LP& input, - GenericShiftGridSet *&gridSetOut) -{ - for( const auto& gridset: grids ) - { - auto grid = gridset->gridAt(input.lam, input.phi); - if( grid ) { - gridSetOut = gridset.get(); - return grid; - } - } - return nullptr; -} - // --------------------------------------------------------------------------- static bool get_grid_values(PJ* P, @@ -87,7 +69,7 @@ static bool get_grid_values(PJ* P, } GenericShiftGridSet* gridset = nullptr; - auto grid = findGrid(Q->grids, lp, gridset); + auto grid = pj_find_generic_grid(Q->grids, lp, gridset); if( !grid ) { return false; } @@ -122,54 +104,20 @@ static bool get_grid_values(PJ* P, return false; } - const auto& extent = grid->extentAndRes(); - double grid_x = (lp.lam - extent.westLon) / extent.resLon; - double grid_y = (lp.phi - extent.southLat) / extent.resLat; - int ix = static_cast(grid_x); - int iy = static_cast(grid_y); - int ix2 = std::min(ix + 1, grid->width() - 1); - int iy2 = std::min(iy + 1, grid->height() - 1); - - float dx1, dy1, dz1; - float dx2, dy2, dz2; - float dx3, dy3, dz3; - float dx4, dy4, dz4; - bool error =( !grid->valueAt(ix, iy, sampleX, dx1) || - !grid->valueAt(ix, iy, sampleY, dy1) || - !grid->valueAt(ix, iy, sampleZ, dz1) || - !grid->valueAt(ix2, iy, sampleX, dx2) || - !grid->valueAt(ix2, iy, sampleY, dy2) || - !grid->valueAt(ix2, iy, sampleZ, dz2) || - !grid->valueAt(ix, iy2, sampleX, dx3) || - !grid->valueAt(ix, iy2, sampleY, dy3) || - !grid->valueAt(ix, iy2, sampleZ, dz3) || - !grid->valueAt(ix2, iy2, sampleX, dx4) || - !grid->valueAt(ix2, iy2, sampleY, dy4) || - !grid->valueAt(ix2, iy2, sampleZ, dz4) ); - if( grid->hasChanged() ) { - if( gridset->reopen(P->ctx) ) { + bool must_retry = false; + if( !pj_bilinear_interpolation_three_samples(grid, lp, + sampleX, sampleY, sampleZ, + dx, dy, dz, + must_retry) ) + { + if( must_retry ) return get_grid_values( P, Q, lp, dx, dy, dz); - } - error = true; - } - if( error ) { return false; } - double frct_lam = grid_x - ix; - double frct_phi = grid_y - iy; - double m10 = frct_lam; - double m11 = m10; - double m01 = 1. - frct_lam; - double m00 = m01; - m11 *= frct_phi; - m01 *= frct_phi; - frct_phi = 1. - frct_phi; - m00 *= frct_phi; - m10 *= frct_phi; - dx = (m00 * dx1 + m10 * dx2 + m01 * dx3 + m11 * dx4) * Q->multiplier; - dy = (m00 * dy1 + m10 * dy2 + m01 * dy3 + m11 * dy4) * Q->multiplier; - dz = (m00 * dz1 + m10 * dz2 + m01 * dz3 + m11 * dz4) * Q->multiplier; + dx *= Q->multiplier; + dy *= Q->multiplier; + dz *= Q->multiplier; return true; } -- cgit v1.2.3