From 7e945f807bc48045582306fa57709ebba2f6a7ad Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 25 Oct 2024 13:12:42 +0200 Subject: [PATCH] proj_factors(): enhance speed when called repeatedly on same compound or projected CRS by caching internal objects within internal state of provided P --- .../development/reference/functions.rst | 4 + src/4D_api.cpp | 167 +++++++++--------- src/factors.cpp | 50 +++--- src/malloc.cpp | 2 + src/proj_internal.h | 6 +- test/unit/gie_self_tests.cpp | 14 +- 6 files changed, 133 insertions(+), 110 deletions(-) diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst index e7e9f6b807..3aa3c303c0 100644 --- a/docs/source/development/reference/functions.rst +++ b/docs/source/development/reference/functions.rst @@ -826,6 +826,10 @@ Various instantiated from a EPSG CRS code. The factors computed will be those of the map projection implied by the transformation from the base geographic CRS of the projected CRS to the projected CRS. + Starting with PROJ 9.6, to improve performance on repeated calls on a + projected CRS object, the above steps will modify the internal state of the + provided P object, and thus calling this function concurrently from multiple + threads on the same P object will no longer be supported. The input geodetic coordinate lp should be such that lp.lam is the longitude in radian, and lp.phi the latitude in radian (thus independently of the diff --git a/src/4D_api.cpp b/src/4D_api.cpp index 1ce26537da..478dab8fa8 100644 --- a/src/4D_api.cpp +++ b/src/4D_api.cpp @@ -2882,90 +2882,99 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) { if (nullptr == P) return factors; - const auto type = proj_get_type(P); - - if (type == PJ_TYPE_COMPOUND_CRS) { - auto ctx = P->ctx; - auto horiz = proj_crs_get_sub_crs(ctx, P, 0); - if (horiz) { - auto ret = proj_factors(horiz, lp); - proj_destroy(horiz); - return ret; + auto pj = P; + auto type = proj_get_type(pj); + + PJ *horiz = nullptr; + if (pj->cached_op_for_proj_factors) { + pj = pj->cached_op_for_proj_factors; + } else { + if (type == PJ_TYPE_COMPOUND_CRS) { + horiz = proj_crs_get_sub_crs(pj->ctx, pj, 0); + pj = horiz; + type = proj_get_type(pj); } - } - if (type == PJ_TYPE_PROJECTED_CRS) { - // If it is a projected CRS, then compute the factors on the conversion - // associated to it. We need to start from a temporary geographic CRS - // using the same datum as the one of the projected CRS, and with - // input coordinates being in longitude, latitude order in radian, - // to be consistent with the expectations of the lp input parameter. - // We also need to create a modified projected CRS with a normalized - // easting/northing axis order in metre, so the resulting operation is - // just a single step pipeline with no axisswap or unitconvert steps. - - auto ctx = P->ctx; - auto geodetic_crs = proj_get_source_crs(ctx, P); - assert(geodetic_crs); - auto pm = proj_get_prime_meridian(ctx, geodetic_crs); - double pm_longitude = 0; - proj_prime_meridian_get_parameters(ctx, pm, &pm_longitude, nullptr, - nullptr); - proj_destroy(pm); - PJ *geogCRSNormalized; - auto cs = proj_create_ellipsoidal_2D_cs( - ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0); - if (pm_longitude != 0) { - auto ellipsoid = proj_get_ellipsoid(ctx, geodetic_crs); - double semi_major_metre = 0; - double inv_flattening = 0; - proj_ellipsoid_get_parameters(ctx, ellipsoid, &semi_major_metre, - nullptr, nullptr, &inv_flattening); - geogCRSNormalized = proj_create_geographic_crs( - ctx, "unname crs", "unnamed datum", proj_get_name(ellipsoid), - semi_major_metre, inv_flattening, "reference prime meridian", 0, - nullptr, 0, cs); - proj_destroy(ellipsoid); - } else { - auto datum = proj_crs_get_datum(ctx, geodetic_crs); - auto datum_ensemble = - proj_crs_get_datum_ensemble(ctx, geodetic_crs); - geogCRSNormalized = proj_create_geographic_crs_from_datum( - ctx, "unnamed crs", datum ? datum : datum_ensemble, cs); - proj_destroy(datum); - proj_destroy(datum_ensemble); + if (type == PJ_TYPE_PROJECTED_CRS) { + // If it is a projected CRS, then compute the factors on the + // conversion associated to it. We need to start from a temporary + // geographic CRS using the same datum as the one of the projected + // CRS, and with input coordinates being in longitude, latitude + // order in radian, to be consistent with the expectations of the lp + // input parameter. We also need to create a modified projected CRS + // with a normalized easting/northing axis order in metre, so the + // resulting operation is just a single step pipeline with no + // axisswap or unitconvert steps. + + auto ctx = pj->ctx; + auto geodetic_crs = proj_get_source_crs(ctx, pj); + assert(geodetic_crs); + auto pm = proj_get_prime_meridian(ctx, geodetic_crs); + double pm_longitude = 0; + proj_prime_meridian_get_parameters(ctx, pm, &pm_longitude, nullptr, + nullptr); + proj_destroy(pm); + PJ *geogCRSNormalized; + auto cs = proj_create_ellipsoidal_2D_cs( + ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0); + if (pm_longitude != 0) { + auto ellipsoid = proj_get_ellipsoid(ctx, geodetic_crs); + double semi_major_metre = 0; + double inv_flattening = 0; + proj_ellipsoid_get_parameters(ctx, ellipsoid, &semi_major_metre, + nullptr, nullptr, + &inv_flattening); + geogCRSNormalized = proj_create_geographic_crs( + ctx, "unname crs", "unnamed datum", + proj_get_name(ellipsoid), semi_major_metre, inv_flattening, + "reference prime meridian", 0, nullptr, 0, cs); + proj_destroy(ellipsoid); + } else { + auto datum = proj_crs_get_datum(ctx, geodetic_crs); + auto datum_ensemble = + proj_crs_get_datum_ensemble(ctx, geodetic_crs); + geogCRSNormalized = proj_create_geographic_crs_from_datum( + ctx, "unnamed crs", datum ? datum : datum_ensemble, cs); + proj_destroy(datum); + proj_destroy(datum_ensemble); + } + proj_destroy(cs); + auto conversion = proj_crs_get_coordoperation(ctx, pj); + auto projCS = proj_create_cartesian_2D_cs( + ctx, PJ_CART2D_EASTING_NORTHING, "metre", 1.0); + auto projCRSNormalized = proj_create_projected_crs( + ctx, nullptr, geodetic_crs, conversion, projCS); + assert(projCRSNormalized); + proj_destroy(geodetic_crs); + proj_destroy(conversion); + proj_destroy(projCS); + auto newOp = proj_create_crs_to_crs_from_pj( + ctx, geogCRSNormalized, projCRSNormalized, nullptr, nullptr); + proj_destroy(geogCRSNormalized); + proj_destroy(projCRSNormalized); + assert(newOp); + // For debugging: + // printf("%s\n", proj_as_proj_string(ctx, newOp, PJ_PROJ_5, + // nullptr)); + + P->cached_op_for_proj_factors = newOp; + pj = newOp; + } else if (type != PJ_TYPE_CONVERSION && + type != PJ_TYPE_TRANSFORMATION && + type != PJ_TYPE_CONCATENATED_OPERATION && + type != PJ_TYPE_OTHER_COORDINATE_OPERATION) { + proj_log_error(P, _("Invalid type for P object")); + proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); + if (horiz) + proj_destroy(horiz); + return factors; } - proj_destroy(cs); - auto conversion = proj_crs_get_coordoperation(ctx, P); - auto projCS = proj_create_cartesian_2D_cs( - ctx, PJ_CART2D_EASTING_NORTHING, "metre", 1.0); - auto projCRSNormalized = proj_create_projected_crs( - ctx, nullptr, geodetic_crs, conversion, projCS); - assert(projCRSNormalized); - proj_destroy(geodetic_crs); - proj_destroy(conversion); - proj_destroy(projCS); - auto newOp = proj_create_crs_to_crs_from_pj( - ctx, geogCRSNormalized, projCRSNormalized, nullptr, nullptr); - proj_destroy(geogCRSNormalized); - proj_destroy(projCRSNormalized); - assert(newOp); - // For debugging: - // printf("%s\n", proj_as_proj_string(ctx, newOp, PJ_PROJ_5, nullptr)); - auto ret = proj_factors(newOp, lp); - proj_destroy(newOp); - return ret; - } - - if (type != PJ_TYPE_CONVERSION && type != PJ_TYPE_TRANSFORMATION && - type != PJ_TYPE_CONCATENATED_OPERATION && - type != PJ_TYPE_OTHER_COORDINATE_OPERATION) { - proj_log_error(P, _("Invalid type for P object")); - proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); - return factors; } - if (pj_factors(lp.lp, P, 0.0, &f)) + const int ret = pj_factors(lp.lp, P, pj, 0.0, &f); + if (horiz) + proj_destroy(horiz); + if (ret) return factors; factors.meridional_scale = f.h; diff --git a/src/factors.cpp b/src/factors.cpp index 1d975055ab..58ed4a3d52 100644 --- a/src/factors.cpp +++ b/src/factors.cpp @@ -12,38 +12,31 @@ #define EPS 1.0e-12 -int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) { +int pj_factors(PJ_LP lp, PJ *toplevel, const PJ *internal, double h, + struct FACTORS *fac) { double cosphi, t, n, r; int err; PJ_COORD coo = {{0, 0, 0, 0}}; coo.lp = lp; - /* Failing the 3 initial checks will most likely be due to */ - /* earlier errors, so we leave errno alone */ - if (nullptr == fac) - return 1; - - if (nullptr == P) - return 1; - if (HUGE_VAL == lp.lam) return 1; /* But from here, we're ready to make our own mistakes */ - err = proj_errno_reset(P); + err = proj_errno_reset(toplevel); /* Indicate that all factors are numerical approximations */ fac->code = 0; /* Check for latitude or longitude overange */ if ((fabs(lp.phi) - M_HALFPI) > EPS) { - proj_log_error(P, _("Invalid latitude")); - proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD); + proj_log_error(toplevel, _("Invalid latitude")); + proj_errno_set(toplevel, PROJ_ERR_COORD_TRANSFM_INVALID_COORD); return 1; } if (fabs(lp.lam) > 10.) { - proj_log_error(P, _("Invalid longitude")); - proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD); + proj_log_error(toplevel, _("Invalid longitude")); + proj_errno_set(toplevel, PROJ_ERR_COORD_TRANSFM_INVALID_COORD); return 1; } @@ -53,8 +46,8 @@ int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) { h = DEFAULT_H; /* If input latitudes are geocentric, convert to geographic */ - if (P->geoc) - lp = pj_geocentric_latitude(P, PJ_INV, coo).lp; + if (internal->geoc) + lp = pj_geocentric_latitude(internal, PJ_INV, coo).lp; /* If latitude + one step overshoots the pole, move it slightly inside, */ /* so the numerical derivative still exists */ @@ -62,14 +55,14 @@ int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) { lp.phi = lp.phi < 0. ? -(M_HALFPI - h) : (M_HALFPI - h); /* Longitudinal distance from central meridian */ - lp.lam -= P->lam0; - if (!P->over) + lp.lam -= internal->lam0; + if (!internal->over) lp.lam = adjlon(lp.lam); /* Derivatives */ - if (pj_deriv(lp, h, P, &(fac->der))) { - proj_log_error(P, _("Invalid latitude or longitude")); - proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD); + if (pj_deriv(lp, h, internal, &(fac->der))) { + proj_log_error(toplevel, _("Invalid latitude or longitude")); + proj_errno_set(toplevel, PROJ_ERR_COORD_TRANSFM_INVALID_COORD); return 1; } @@ -78,13 +71,13 @@ int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) { fac->h = hypot(fac->der.x_p, fac->der.y_p); fac->k = hypot(fac->der.x_l, fac->der.y_l) / cosphi; - if (P->es != 0.0) { + if (internal->es != 0.0) { t = sin(lp.phi); - t = 1. - P->es * t * t; + t = 1. - internal->es * t * t; n = sqrt(t); - fac->h *= t * n / P->one_es; + fac->h *= t * n / internal->one_es; fac->k *= n; - r = t * t / P->one_es; + r = t * t / internal->one_es; } else r = 1.; @@ -96,7 +89,7 @@ int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) { cosphi; /* Meridian-parallel angle (theta prime) */ - fac->thetap = aasin(P->ctx, fac->s / (fac->h * fac->k)); + fac->thetap = aasin(internal->ctx, fac->s / (fac->h * fac->k)); /* Tissot ellipse axis */ t = fac->k * fac->k + fac->h * fac->h; @@ -107,8 +100,9 @@ int pj_factors(PJ_LP lp, const PJ *P, double h, struct FACTORS *fac) { fac->a = 0.5 * (fac->a + t); /* Angular distortion */ - fac->omega = 2. * aasin(P->ctx, (fac->a - fac->b) / (fac->a + fac->b)); + fac->omega = + 2. * aasin(internal->ctx, (fac->a - fac->b) / (fac->a + fac->b)); - proj_errno_restore(P, err); + proj_errno_restore(toplevel, err); return 0; } diff --git a/src/malloc.cpp b/src/malloc.cpp index 5b79d60a27..401a0a11d2 100644 --- a/src/malloc.cpp +++ b/src/malloc.cpp @@ -166,6 +166,8 @@ PJ *pj_default_destructor(PJ *P, int errlev) { /* Destructor */ proj_destroy(P->hgridshift); proj_destroy(P->vgridshift); + proj_destroy(P->cached_op_for_proj_factors); + free(static_cast(P->opaque)); delete P; return nullptr; diff --git a/src/proj_internal.h b/src/proj_internal.h index 0384d0b61e..f51337296b 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -675,6 +675,9 @@ struct PJconsts { true; /* to remove in PROJ 10? */ bool skipNonInstantiable = true; + // Used internally by proj_factors() + PJ *cached_op_for_proj_factors = nullptr; + /************************************************************************************* E N D O F G E N E R A L P A R A M E T E R S T R U C T @@ -916,7 +919,8 @@ COMPLEX pj_zpoly1(COMPLEX, const COMPLEX *, int); COMPLEX pj_zpolyd1(COMPLEX, const COMPLEX *, int, COMPLEX *); int pj_deriv(PJ_LP, double, const PJ *, struct DERIVS *); -int pj_factors(PJ_LP, const PJ *, double, struct FACTORS *); +int pj_factors(PJ_LP, PJ *toplevel, const PJ *internal, double, + struct FACTORS *); void *proj_mdist_ini(double); double proj_mdist(double, double, double, const void *); diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp index 3807b452a4..259eca263f 100644 --- a/test/unit/gie_self_tests.cpp +++ b/test/unit/gie_self_tests.cpp @@ -611,9 +611,19 @@ TEST(gie, info_functions) { PJ_COORD c; c.lp.lam = proj_torad(10.729030600); c.lp.phi = proj_torad(59.916494500); - const auto factors2 = proj_factors(P, c); - EXPECT_NEAR(factors2.meridional_scale, 1 - 28.54587730 * 1e-5, 1e-8); + { + const auto factors2 = proj_factors(P, c); + EXPECT_NEAR(factors2.meridional_scale, 1 - 28.54587730 * 1e-5, + 1e-8); + } + + // Try again to test caching of internal objects + { + const auto factors2 = proj_factors(P, c); + EXPECT_NEAR(factors2.meridional_scale, 1 - 28.54587730 * 1e-5, + 1e-8); + } proj_destroy(P); }