From f08c7835320f0408a1f18b1ca3c0d64e6caec4bf Mon Sep 17 00:00:00 2001 From: flipflip Date: Thu, 29 Oct 2026 08:25:10 +0100 Subject: [PATCH] Add options, CRS info and thread safety to Transformer - Init() takes TransformerOpts: precise (ALLOW_BALLPARK=NO, ONLY_BEST=YES, instantiable check), min_accuracy (ACCURACY=x), normalize_axes (proj_normalize_for_visualization()), allow_network (default off) and area of interest - SourceCrs() and TargetCrs() return the CRS type, name and compound parts - All methods are thread-safe, Transform() is const - Transform() fails on PROJ errors and non-finite results --- fpsdk_common/include/fpsdk_common/trafo.hpp | 150 +++++++++++-- fpsdk_common/src/trafo.cpp | 221 +++++++++++++++++--- fpsdk_common/test/trafo_test.cpp | 168 +++++++++++++++ 3 files changed, 490 insertions(+), 49 deletions(-) diff --git a/fpsdk_common/include/fpsdk_common/trafo.hpp b/fpsdk_common/include/fpsdk_common/trafo.hpp index 0fb4e720..ae1fd01b 100644 --- a/fpsdk_common/include/fpsdk_common/trafo.hpp +++ b/fpsdk_common/include/fpsdk_common/trafo.hpp @@ -20,7 +20,10 @@ /* LIBC/STL */ #include +#include +#include #include +#include /* EXTERNAL */ #include "fpsdk_common/ext/eigen_core.hpp" @@ -211,9 +214,104 @@ Eigen::Vector3d LlhDegToRad(const Eigen::Vector3d& llh_deg); Eigen::Vector3d LlhRadToDeg(const Eigen::Vector3d& llh_rad); #if FPSDK_USE_PROJ || defined(_DOXYGEN_) +/** + * @brief Coordinate reference system (CRS) types + * + * (Like PROJ's PJ_TYPE_GEOCENTRIC_CRS etc., see their docu for details) + */ +enum class CrsType : int +{ // clang-format off + UNKNOWN = 0, //!< Unknown (or not initialised) + GEOCENTRIC = 1, //!< Geocentric (ECEF x/y/z), e.g. EPSG:4978 + GEOGRAPHIC_2D = 2, //!< Geographic 2D (lat/lon), e.g. EPSG:4326 + GEOGRAPHIC_3D = 3, //!< Geographic 3D (lat/lon/height), e.g. EPSG:4979 + PROJECTED = 4, //!< Projected (east/north), e.g. EPSG:2056 + VERTICAL = 5, //!< Vertical (height), e.g. EPSG:5703 + COMPOUND = 6, //!< Compound (horizontal + vertical), e.g. EPSG:4258+5703 + OTHER = 7, //!< Other type (engineering, temporal, ...) +}; // clang-format on + +/** + * @brief Stringify CRS type + * + * @param[in] type The CRS type + * + * @returns a unique identifier string for the CRS type + */ +const char* CrsTypeStr(const CrsType type); + +/** + * @brief Coordinate reference system (CRS) information + */ +struct CrsInfo +{ + CrsType type_ = CrsType::UNKNOWN; //!< CRS type + std::string name_; //!< CRS name, e.g. "WGS 84" + std::vector sub_types_; //!< Types of the parts of a COMPOUND CRS, e.g. { GEOGRAPHIC_2D, VERTICAL } +}; + +/** + * @brief Transformer options + * + * The defaults select the most suitable operation PROJ can find, which may be an approximate ("ballpark") one. + */ +struct TransformerOpts +{ + /** + * @brief Require a precise transformation + * + * Do not allow approximate ("ballpark") operations (ALLOW_BALLPARK=NO), fail if the best operation is not + * available, e.g. because of missing grids (ONLY_BEST=YES), and check that the operation can be instantiated. + */ + bool precise_ = false; + + /** + * @brief Minimal accuracy [m] of the operation (ACCURACY=x), 0.0 = any accuracy + * + * Operations with unknown accuracy, such as ballpark operations, are not used if this is set. + */ + double min_accuracy_ = 0.0; + + /** + * @brief Use "traditional GIS" axis order + * + * Use longitude/latitude (and east/north) order for all CRS, regardless of the axis order of the CRS (see + * proj_normalize_for_visualization()). For example, EPSG:4979 natively uses latitude/longitude order. + */ + bool normalize_axes_ = false; + + /** + * @brief Allow network access + * + * Allow PROJ to download missing grids from the PROJ content delivery network (see + * https://proj.org/en/9.4/usage/network.html). This overrides the PROJ configuration (proj.ini, PROJ_NETWORK). + */ + bool allow_network_ = false; + + /** + * @brief Area of interest (bounding box) + */ + struct Area + { + double west_ = 0.0; //!< Western-most longitude [deg] + double south_ = 0.0; //!< Southern-most latitude [deg] + double east_ = 0.0; //!< Eastern-most longitude [deg] + double north_ = 0.0; //!< Northern-most latitude [deg] + }; + + /** + * @brief Area of interest, used to select the operation + * + * If unset PROJ may select a different operation for each coordinate, depending on the coordinate's position. + */ + std::optional area_; +}; + /** * @brief "Universal" coordinate transformer, backed by PROJ * + * All methods are thread-safe. Only one thread at a time can transform, other threads wait. + * * @note This is only available if compiled with PROJ, see @ref FPSDK_BUILD_DEPS. */ class Transformer : private types::NoCopyNoMove @@ -231,17 +329,32 @@ class Transformer : private types::NoCopyNoMove /** * @brief Initialise transformer * - * The source and target CRS specification can be anything that proj_create_crs_to_crs() understands (see - * https://proj.org/en/9.4/development/reference/functions.html#c.proj_create_crs_to_crs). You may have to change - * the PROJ configuration (proj.ini) and/or update some PROJ data yourself to support some CRS. See the PROJ - * documentation for details (see https://proj.org/en/9.4/resource_files.html). + * The source and target CRS specification can be anything that proj_create() understands (see + * https://proj.org/en/9.4/development/reference/functions.html#c.proj_create). You may have to change the PROJ + * configuration (proj.ini) and/or update some PROJ data yourself to support some CRS. See the PROJ documentation + * for details (see https://proj.org/en/9.4/resource_files.html). * * @param[in] source_crs Source coordinate reference system, e.g. "EPSG:4326" * @param[in] target_crs Target coordinate reference system, e.g. "EPSG:2056" + * @param[in] opts Options * - * @returns true on success, false otherwise (bad params, missing data, ...) + * @returns true on success, false otherwise (bad params, missing data, already initialised, no PROJ, ...) */ - bool Init(const std::string& source_crs, const std::string& target_crs); + bool Init(const std::string& source_crs, const std::string& target_crs, const TransformerOpts& opts = {}); + + /** + * @brief Get source CRS information + * + * @returns the source CRS information (type UNKNOWN if not initialised) + */ + CrsInfo SourceCrs() const; + + /** + * @brief Get target CRS information + * + * @returns the target CRS information (type UNKNOWN if not initialised) + */ + CrsInfo TargetCrs() const; /** * @brief Transform coordinates @@ -249,9 +362,9 @@ class Transformer : private types::NoCopyNoMove * @param[in,out] inout Coordinates to transform, replaced with result * @param[in] inv Do the inverse transformation (true), default is forward (false) * - * @returns true on success, false otherwise + * @returns true on success, false otherwise (not initialised, transformation failed) */ - bool Transform(Eigen::Vector3d& inout, const bool inv = false); + bool Transform(Eigen::Vector3d& inout, const bool inv = false) const; /** * @brief Transform coordinates @@ -260,9 +373,9 @@ class Transformer : private types::NoCopyNoMove * @param[out] out Transformed coordinates * @param[in] inv Do the inverse transformation (true), default is forward (false) * - * @returns true on success, false otherwise + * @returns true on success, false otherwise (not initialised, transformation failed) */ - bool Transform(const Eigen::Vector3d& in, Eigen::Vector3d& out, const bool inv = false); + bool Transform(const Eigen::Vector3d& in, Eigen::Vector3d& out, const bool inv = false) const; /** * @brief Transform coordinates (with time) @@ -271,15 +384,20 @@ class Transformer : private types::NoCopyNoMove * @param[out] out Transformed coordinates (with time) * @param[in] inv Do the inverse transformation (true), default is forward (false) * - * @returns true on success, false otherwise + * @returns true on success, false otherwise (not initialised, transformation failed) */ - bool Transform(const Eigen::Vector4d& in, Eigen::Vector4d& out, const bool inv = false); + bool Transform(const Eigen::Vector4d& in, Eigen::Vector4d& out, const bool inv = false) const; private: - const std::string name_; //!< Name, for debugging - bool pj_init_; //!< Proj is initialised - void* pj_ctx_; //!< Proj context - void* pj_tf_; //!< Proj transformation object + const std::string name_; //!< Name, for debugging + mutable std::mutex mutex_; //!< One thread at a time uses PROJ + bool pj_init_; //!< Proj is initialised + void* pj_ctx_; //!< Proj context + void* pj_tf_; //!< Proj transformation object + CrsInfo source_crs_; //!< Source CRS information + CrsInfo target_crs_; //!< Target CRS information + + bool _Transform(const Eigen::Vector4d& in, Eigen::Vector4d& out, const bool inv) const; }; #endif // FPSDK_USE_PROJ diff --git a/fpsdk_common/src/trafo.cpp b/fpsdk_common/src/trafo.cpp index aa1b08ac..8221364a 100644 --- a/fpsdk_common/src/trafo.cpp +++ b/fpsdk_common/src/trafo.cpp @@ -12,7 +12,12 @@ */ /* LIBC/STL */ +#include #include +#include +#include +#include +#include /* EXTERNAL */ #if FPSDK_USE_PROJ @@ -236,8 +241,26 @@ Eigen::Vector3d LlhRadToDeg(const Eigen::Vector3d& llh_rad) } /* ****************************************************************************************************************** */ + #if FPSDK_USE_PROJ +const char* CrsTypeStr(const CrsType type) +{ + switch (type) { // clang-format off + case CrsType::UNKNOWN: return "UNKNOWN"; + case CrsType::GEOCENTRIC: return "GEOCENTRIC"; + case CrsType::GEOGRAPHIC_2D: return "GEOGRAPHIC_2D"; + case CrsType::GEOGRAPHIC_3D: return "GEOGRAPHIC_3D"; + case CrsType::PROJECTED: return "PROJECTED"; + case CrsType::VERTICAL: return "VERTICAL"; + case CrsType::COMPOUND: return "COMPOUND"; + case CrsType::OTHER: return "OTHER"; + } // clang-format on + return "?"; +} + +// --------------------------------------------------------------------------------------------------------------------- + Transformer::Transformer(const std::string& name) /* clang-format off */ : name_ { name }, pj_init_ { false }, @@ -248,6 +271,7 @@ Transformer::Transformer(const std::string& name) /* clang-format off */ : Transformer::~Transformer() { + std::unique_lock lock(mutex_); if (pj_init_) { pj_init_ = false; proj_destroy((PJ*)pj_tf_); @@ -270,75 +294,206 @@ static void ProjLogger(void* data, int level, const char* msg) } // clang-format on } -bool Transformer::Init(const std::string& source_crs, const std::string& target_crs) +static const char* ProjErrorStr(PJ_CONTEXT* ctx, const char* fallback = "unknown error") +{ + const int err = proj_context_errno(ctx); + const char* err_str = (err != 0 ? proj_context_errno_string(ctx, err) : NULL); + return ((err_str != NULL) && (err_str[0] != '\0')) ? err_str : fallback; +} + +static CrsType ProjTypeToCrsType(const PJ_TYPE type) +{ + switch (type) { // clang-format off + case PJ_TYPE_UNKNOWN: return CrsType::UNKNOWN; + case PJ_TYPE_GEOCENTRIC_CRS: return CrsType::GEOCENTRIC; + case PJ_TYPE_GEOGRAPHIC_2D_CRS: return CrsType::GEOGRAPHIC_2D; + case PJ_TYPE_GEOGRAPHIC_3D_CRS: return CrsType::GEOGRAPHIC_3D; + case PJ_TYPE_PROJECTED_CRS: return CrsType::PROJECTED; + case PJ_TYPE_VERTICAL_CRS: return CrsType::VERTICAL; + case PJ_TYPE_COMPOUND_CRS: return CrsType::COMPOUND; + default: return CrsType::OTHER; + } // clang-format on +} + +static CrsInfo GetCrsInfo(PJ_CONTEXT* ctx, const PJ* crs) +{ + CrsInfo info; + info.type_ = ProjTypeToCrsType(proj_get_type(crs)); + const char* name = proj_get_name(crs); + info.name_ = (name != NULL ? name : ""); + if (info.type_ == CrsType::COMPOUND) { + for (int ix = 0; true; ix++) { + PJ* sub_crs = proj_crs_get_sub_crs(ctx, crs, ix); + if (sub_crs == NULL) { + break; + } + info.sub_types_.push_back(ProjTypeToCrsType(proj_get_type(sub_crs))); + proj_destroy(sub_crs); + } + } + return info; +} + +bool Transformer::Init(const std::string& source_crs, const std::string& target_crs, const TransformerOpts& opts) { + std::unique_lock lock(mutex_); if (pj_init_) { return false; } - DEBUG("Transformer(%s) [%s] [%s]", name_.c_str(), source_crs.c_str(), target_crs.c_str()); + DEBUG("Transformer(%s) [%s] [%s] precise=%s min_accuracy=%g normalize_axes=%s allow_network=%s area=%s", + name_.c_str(), source_crs.c_str(), target_crs.c_str(), opts.precise_ ? "true" : "false", opts.min_accuracy_, + opts.normalize_axes_ ? "true" : "false", opts.allow_network_ ? "true" : "false", + opts.area_ + ? Sprintf("%g/%g/%g/%g", opts.area_->west_, opts.area_->south_, opts.area_->east_, opts.area_->north_) + .c_str() + : "none"); TRACE("Transformer(%s) create proj context", name_.c_str()); PJ_CONTEXT* ctx = proj_context_create(); + if (ctx == NULL) { + WARNING("Transformer(%s) failed creating PROJ context", name_.c_str()); + return false; + } # ifdef NDEBUG proj_log_level(ctx, PJ_LOG_DEBUG); # else proj_log_level(ctx, PJ_LOG_TRACE); # endif proj_log_func(ctx, (void*)name_.c_str(), ProjLogger); + if ((proj_context_set_enable_network(ctx, opts.allow_network_ ? 1 : 0) == 0) && opts.allow_network_) { + WARNING("Transformer(%s) PROJ has no network support", name_.c_str()); + } + + // Source and target CRS + PJ* source = proj_create(ctx, source_crs.c_str()); + if (source == NULL) { + WARNING("Transformer(%s) bad source CRS %s: %s", name_.c_str(), source_crs.c_str(), ProjErrorStr(ctx)); + proj_context_destroy(ctx); + return false; + } + PJ* target = proj_create(ctx, target_crs.c_str()); + if (target == NULL) { + WARNING("Transformer(%s) bad target CRS %s: %s", name_.c_str(), target_crs.c_str(), ProjErrorStr(ctx)); + proj_destroy(source); + proj_context_destroy(ctx); + return false; + } + CrsInfo source_info = GetCrsInfo(ctx, source); + CrsInfo target_info = GetCrsInfo(ctx, target); + + // Options for proj_create_crs_to_crs_from_pj() + std::vector options; + if (opts.precise_) { + options.push_back("ALLOW_BALLPARK=NO"); + options.push_back("ONLY_BEST=YES"); + } + if (opts.min_accuracy_ > 0.0) { + options.push_back(Sprintf("ACCURACY=%.6f", opts.min_accuracy_)); + } + std::vector options_ptrs; + for (const auto& option : options) { + options_ptrs.push_back(option.c_str()); + } + options_ptrs.push_back(NULL); + + PJ_AREA* area = NULL; + if (opts.area_) { + area = proj_area_create(); + proj_area_set_bbox(area, opts.area_->west_, opts.area_->south_, opts.area_->east_, opts.area_->north_); + } TRACE("Transformer(%s) create proj transformer", name_.c_str()); - PJ* tf = proj_create_crs_to_crs(ctx, source_crs.c_str(), target_crs.c_str(), NULL); + PJ* tf = proj_create_crs_to_crs_from_pj(ctx, source, target, area, options_ptrs.data()); + if (area != NULL) { + proj_area_destroy(area); + } + proj_destroy(source); + proj_destroy(target); if (tf == NULL) { - const int err = proj_context_errno(ctx); - const char* err_str = proj_context_errno_string(ctx, err); - WARNING("Transformer(%s) %s -> %s: %s", name_.c_str(), source_crs.c_str(), target_crs.c_str(), err_str); + WARNING("Transformer(%s) %s -> %s: %s", name_.c_str(), source_crs.c_str(), target_crs.c_str(), + ProjErrorStr(ctx, "no suitable operation")); + proj_context_destroy(ctx); + return false; + } + if (opts.precise_ && !proj_coordoperation_is_instantiable(ctx, tf)) { + WARNING("Transformer(%s) %s -> %s: operation not available (missing grids?)", name_.c_str(), source_crs.c_str(), + target_crs.c_str()); + proj_destroy(tf); proj_context_destroy(ctx); return false; } + if (opts.normalize_axes_) { + PJ* tf_norm = proj_normalize_for_visualization(ctx, tf); + proj_destroy(tf); + if (tf_norm == NULL) { + WARNING("Transformer(%s) %s -> %s: failed normalizing axes: %s", name_.c_str(), source_crs.c_str(), + target_crs.c_str(), ProjErrorStr(ctx)); + proj_context_destroy(ctx); + return false; + } + tf = tf_norm; + } pj_ctx_ = (void*)ctx; pj_tf_ = (void*)tf; + source_crs_ = source_info; + target_crs_ = target_info; pj_init_ = true; TRACE("Transformer(%s) init ok", name_.c_str()); return true; } -bool Transformer::Transform(const Eigen::Vector4d& in, Eigen::Vector4d& out, const bool inv) +// --------------------------------------------------------------------------------------------------------------------- + +CrsInfo Transformer::SourceCrs() const { - if (pj_init_) { - const PJ_COORD src = { { in.x(), in.y(), in.z(), in.w() } }; - const PJ_COORD dst = proj_trans((PJ*)pj_tf_, inv ? PJ_INV : PJ_FWD, src); - out.x() = dst.v[0]; - out.y() = dst.v[1]; - out.z() = dst.v[2]; - out.w() = dst.v[3]; - return true; - } - return false; + std::unique_lock lock(mutex_); + return source_crs_; } -bool Transformer::Transform(Eigen::Vector3d& inout, const bool inv) +CrsInfo Transformer::TargetCrs() const { - if (pj_init_) { - const PJ_COORD src = { { inout.x(), inout.y(), inout.z(), std::numeric_limits::infinity() } }; - const PJ_COORD dst = proj_trans((PJ*)pj_tf_, inv ? PJ_INV : PJ_FWD, src); - inout.x() = dst.v[0]; - inout.y() = dst.v[1]; - inout.z() = dst.v[2]; - return true; + std::unique_lock lock(mutex_); + return target_crs_; +} + +// --------------------------------------------------------------------------------------------------------------------- + +bool Transformer::_Transform(const Eigen::Vector4d& in, Eigen::Vector4d& out, const bool inv) const +{ + std::unique_lock lock(mutex_); + if (!pj_init_) { + return false; } - return false; + PJ* tf = (PJ*)pj_tf_; + proj_errno_reset(tf); + const PJ_COORD src = { { in.x(), in.y(), in.z(), in.w() } }; + const PJ_COORD dst = proj_trans(tf, inv ? PJ_INV : PJ_FWD, src); + const int err = proj_errno(tf); + if ((err != 0) || !std::isfinite(dst.v[0]) || !std::isfinite(dst.v[1]) || !std::isfinite(dst.v[2])) { + TRACE("Transformer(%s) transform failed: %s", name_.c_str(), err != 0 ? proj_errno_string(err) : "bad result"); + return false; + } + out = { dst.v[0], dst.v[1], dst.v[2], dst.v[3] }; + return true; } -bool Transformer::Transform(const Eigen::Vector3d& in, Eigen::Vector3d& out, const bool inv) +bool Transformer::Transform(const Eigen::Vector4d& in, Eigen::Vector4d& out, const bool inv) const { - if (pj_init_) { - const PJ_COORD src = { { in.x(), in.y(), in.z(), std::numeric_limits::infinity() } }; - const PJ_COORD dst = proj_trans((PJ*)pj_tf_, inv ? PJ_INV : PJ_FWD, src); - out.x() = dst.v[0]; - out.y() = dst.v[1]; - out.z() = dst.v[2]; + return _Transform(in, out, inv); +} + +bool Transformer::Transform(Eigen::Vector3d& inout, const bool inv) const +{ + return Transform(inout, inout, inv); +} + +bool Transformer::Transform(const Eigen::Vector3d& in, Eigen::Vector3d& out, const bool inv) const +{ + Eigen::Vector4d out4; + if (_Transform({ in.x(), in.y(), in.z(), std::numeric_limits::infinity() }, out4, inv)) { + out = out4.head<3>(); return true; } return false; diff --git a/fpsdk_common/test/trafo_test.cpp b/fpsdk_common/test/trafo_test.cpp index e1e5d4fe..02cf4cab 100644 --- a/fpsdk_common/test/trafo_test.cpp +++ b/fpsdk_common/test/trafo_test.cpp @@ -13,6 +13,9 @@ /* LIBC/STL */ #include +#include +#include +#include /* EXTERNAL */ #include @@ -324,6 +327,171 @@ TEST(TrafoTest, Trafo) // https://epsg.io/7922 ETRF93 xyz } +TEST(TrafoTest, TrafoCrsInfo) +{ + { + Transformer trafo("uninitialised"); + EXPECT_EQ(trafo.SourceCrs().type_, CrsType::UNKNOWN); + EXPECT_EQ(trafo.TargetCrs().type_, CrsType::UNKNOWN); + } + { + Transformer trafo("ecef_to_llh"); + ASSERT_TRUE(trafo.Init("EPSG:4978", "EPSG:4979")); + EXPECT_EQ(trafo.SourceCrs().type_, CrsType::GEOCENTRIC); + EXPECT_EQ(trafo.SourceCrs().name_, "WGS 84"); + EXPECT_TRUE(trafo.SourceCrs().sub_types_.empty()); + EXPECT_EQ(trafo.TargetCrs().type_, CrsType::GEOGRAPHIC_3D); + } + { + Transformer trafo("llh_to_ch1903"); + ASSERT_TRUE(trafo.Init("EPSG:4326", "EPSG:2056")); + EXPECT_EQ(trafo.SourceCrs().type_, CrsType::GEOGRAPHIC_2D); + EXPECT_EQ(trafo.TargetCrs().type_, CrsType::PROJECTED); + } + { + // https://epsg.io/4936 ETRS89 x/y/z, https://epsg.io/4258 ETRS89 lat/lon, https://epsg.io/5730 EVRF2000 height + Transformer trafo("etrs89_to_compound"); + ASSERT_TRUE(trafo.Init("EPSG:4936", "EPSG:4258+5730")); + EXPECT_EQ(trafo.TargetCrs().type_, CrsType::COMPOUND); + EXPECT_EQ(trafo.TargetCrs().sub_types_, std::vector({ CrsType::GEOGRAPHIC_2D, CrsType::VERTICAL })); + } + EXPECT_STREQ(CrsTypeStr(CrsType::GEOCENTRIC), "GEOCENTRIC"); +} + +TEST(TrafoTest, TrafoInitFail) +{ + { + Transformer trafo("bad_crs"); + EXPECT_FALSE(trafo.Init("EPSG:4978", "not a CRS")); + EXPECT_FALSE(trafo.Init("not a CRS", "EPSG:4979")); + EXPECT_FALSE(trafo.Init("EPSG:4978", "EPSG:4979+3855")); // Geographic 3D + vertical is not allowed + } + { + Transformer trafo("init_twice"); + EXPECT_TRUE(trafo.Init("EPSG:4978", "EPSG:4979")); + EXPECT_FALSE(trafo.Init("EPSG:4978", "EPSG:4979")); + } +} + +TEST(TrafoTest, TrafoOpts) +{ + const Eigen::Vector3d llh_deg = { 47.4002984167, 8.45036571389, 459.410 }; + const Eigen::Vector3d xyz = TfEcefWgs84Llh(LlhDegToRad(llh_deg)); + const Eigen::Vector3d err = { 1e-8, 1e-8, 1e-3 }; + + // Native axis order of EPSG:4979 is lat/lon + { + Transformer trafo("native"); + ASSERT_TRUE(trafo.Init("EPSG:4978", "EPSG:4979")); + Eigen::Vector3d p; + ASSERT_TRUE(trafo.Transform(xyz, p)); + EXPECT_NEAR_EIGEN_VECTOR3D(llh_deg, p, err); + } + // Normalized axis order is lon/lat, for EPSG:4979 as well as for OGC:CRS84h (which natively is lon/lat) + const Eigen::Vector3d lonlath_deg = { llh_deg.y(), llh_deg.x(), llh_deg.z() }; + for (const char* target_crs : { "EPSG:4979", "OGC:CRS84h" }) { + Transformer trafo("normalized"); + TransformerOpts opts; + opts.normalize_axes_ = true; + ASSERT_TRUE(trafo.Init("EPSG:4978", target_crs, opts)); + Eigen::Vector3d p; + ASSERT_TRUE(trafo.Transform(xyz, p)); + EXPECT_NEAR_EIGEN_VECTOR3D(lonlath_deg, p, err); + Eigen::Vector3d q; + ASSERT_TRUE(trafo.Transform(p, q, true)); + EXPECT_NEAR_EIGEN_VECTOR3D(xyz, q, Eigen::Vector3d(1e-3, 1e-3, 1e-3)); + } + // Precise + { + Transformer trafo("precise"); + TransformerOpts opts; + opts.precise_ = true; + EXPECT_TRUE(trafo.Init("EPSG:4978", "EPSG:4979", opts)); + } + // Network (no grids are required, so no network access happens) + { + Transformer trafo("network"); + TransformerOpts opts; + opts.allow_network_ = true; + EXPECT_TRUE(trafo.Init("EPSG:4978", "EPSG:4979", opts)); + } + // ETRS89 to WGS84 has no operation with an accuracy better than 1m + { + Transformer trafo("min_accuracy_ok"); + TransformerOpts opts; + opts.min_accuracy_ = 2.0; + EXPECT_TRUE(trafo.Init("EPSG:4936", "EPSG:4978", opts)); + } + { + Transformer trafo("min_accuracy_fail"); + TransformerOpts opts; + opts.min_accuracy_ = 0.01; + EXPECT_FALSE(trafo.Init("EPSG:4936", "EPSG:4978", opts)); + } + // Area of interest + { + Transformer trafo("area"); + TransformerOpts opts; + opts.area_ = TransformerOpts::Area{ 5.9, 45.8, 10.5, 47.8 }; // Switzerland + ASSERT_TRUE(trafo.Init("EPSG:4326", "EPSG:2056", opts)); + Eigen::Vector3d p; + ASSERT_TRUE(trafo.Transform(llh_deg, p)); + EXPECT_NEAR_EIGEN_VECTOR3D( + Eigen::Vector3d(2676373.211, 1250433.955, 459.410), p, Eigen::Vector3d(1e-3, 1e-3, 1e-3)); + } +} + +TEST(TrafoTest, TrafoTransformFail) +{ + { + Transformer trafo("uninitialised"); + Eigen::Vector3d p = { 1.0, 2.0, 3.0 }; + EXPECT_FALSE(trafo.Transform(p)); + } + { + Transformer trafo("nan"); + ASSERT_TRUE(trafo.Init("EPSG:4978", "EPSG:4979")); + const double nan = std::numeric_limits::quiet_NaN(); + Eigen::Vector3d p = { nan, 1.0, 1.0 }; + EXPECT_FALSE(trafo.Transform(p)); + Eigen::Vector4d q; + EXPECT_FALSE(trafo.Transform(Eigen::Vector4d(nan, 1.0, 1.0, 0.0), q)); + } +} + +TEST(TrafoTest, TrafoThreads) +{ + const Eigen::Vector3d llh_deg = { 47.4002984167, 8.45036571389, 459.410 }; + const Eigen::Vector3d xyz = TfEcefWgs84Llh(LlhDegToRad(llh_deg)); + const Eigen::Vector3d err = { 1e-8, 1e-8, 1e-3 }; + + Transformer trafo_init("threads"); + ASSERT_TRUE(trafo_init.Init("EPSG:4978", "EPSG:4979")); + const Transformer& trafo = trafo_init; + + // Transform in threads different from the one that initialised the transformer, concurrently + constexpr int NUM_THREADS = 4; + constexpr int NUM_TRAFOS = 1000; + std::vector num_ok(NUM_THREADS, 0); + std::vector threads; + for (int t_ix = 0; t_ix < NUM_THREADS; t_ix++) { + threads.emplace_back([&, t_ix]() { + for (int ix = 0; ix < NUM_TRAFOS; ix++) { + Eigen::Vector3d p; + if (trafo.Transform(xyz, p) && ((llh_deg - p).cwiseAbs().array() <= err.array()).all()) { + num_ok[t_ix]++; + } + } + }); + } + for (auto& thread : threads) { + thread.join(); + } + for (int t_ix = 0; t_ix < NUM_THREADS; t_ix++) { + EXPECT_EQ(num_ok[t_ix], NUM_TRAFOS); + } +} + # if 0 # define INFO_VV(n, v1, v2) \ INFO("%-20s v1: %.3f %.3f %.3f v2: %.3f %.3f %.3f d=%.3f", n, (v1).x(), (v1).y(), (v1).z(), (v2).x(), \