From e415d8b17ee734316543cf9a21f3a288f5f97a22 Mon Sep 17 00:00:00 2001 From: jmmuller Date: Tue, 3 Sep 2024 13:49:02 +0200 Subject: [PATCH] SysCo up --- MMVII/include/MMVII_SysCo.h | 15 ++++++++++++ MMVII/src/SysCo/SysCo.cpp | 46 +++++++++++++++++++++++++------------ 2 files changed, 46 insertions(+), 15 deletions(-) diff --git a/MMVII/include/MMVII_SysCo.h b/MMVII/include/MMVII_SysCo.h index e40dd304c9..5d4a941e79 100644 --- a/MMVII/include/MMVII_SysCo.h +++ b/MMVII/include/MMVII_SysCo.h @@ -4,6 +4,7 @@ #include "MMVII_DeclareCste.h" #include "MMVII_Mappings.h" #include "MMVII_AllClassDeclare.h" +#include "MMVII_Geom3D.h" struct pj_ctx; typedef struct pj_ctx PJ_CONTEXT; @@ -39,6 +40,7 @@ void AddData(const cAuxAr2007 & anAux, cSysCoData & aSysCoData); class cSysCo : public cDataInvertibleMapping { public : + typedef cIsometry3D tPoseR; virtual ~cSysCo(); // do not copy and move because of PJ* (not needed via tPtrSysCo) @@ -50,9 +52,14 @@ public : virtual tPt Value(const tPt &) const override = 0; //< to GeoC virtual tPt Inverse(const tPt &) const override = 0; //< from GeoC + + tPt toGeoG(const tPt & aPtIn) const; //< uses mPJ_GeoC2Geog + tPt fromGeoG(const tPt &aPtInGeoG) const; //< uses mPJ_GeoC2Geog + virtual cRotation3D getRot2Vertical(const tPt &) const; //< get rotation from SysCo origin to vertical at point virtual tREAL8 getRadiusApprox(const tPt &in) const; //< approximate earth total curvature radius at a point virtual tREAL8 getDistHzApprox(const tPt & aPtA, const tPt & aPtB) const; //< approximate horizontal distance (along ellipsoid) from one point to an other + virtual const tPoseR *getTranfo2GeoC() const; static tPtrSysCo MakeSysCo(const std::string &aDef, bool aDebug=false); //< factory from a SysCo definition static tPtrSysCo makeRTL(const cPt3dr & anOrigin, const std::string & aSysCoInDef); @@ -64,6 +71,9 @@ public : eSysCo getType() const { return mType; } bool isEuclidian() const; + tREAL8 getEllipsoid_a() const {return semi_axis;} + tREAL8 getEllipsoid_e2() const {return e2;} + tREAL8 getEllipsoid_b() const {return b;} protected : cSysCo(bool aDebug); cSysCo(const std::string & def, bool aDebug); @@ -72,6 +82,11 @@ protected : PJ_CONTEXT* mPJContext; PJ* mPJ_GeoC2Geog; //< for generic use bool mDebug; //< show debug messages + + //GRS80 + const tREAL8 semi_axis = 6378137; + const tREAL8 e2 = 0.00669438; + const tREAL8 b = semi_axis * sqrt(1.-e2); }; inline bool operator==(const cSysCo& lhs, const cSysCo& rhs) { return lhs.Def() == rhs.Def(); } diff --git a/MMVII/src/SysCo/SysCo.cpp b/MMVII/src/SysCo/SysCo.cpp index b46fd6e502..85eca20f5a 100644 --- a/MMVII/src/SysCo/SysCo.cpp +++ b/MMVII/src/SysCo/SysCo.cpp @@ -1,8 +1,6 @@ #include "MMVII_SysCo.h" #include -#include "MMVII_Geom3D.h" - namespace MMVII { @@ -152,7 +150,6 @@ class cSysCoLEuc : public cSysCo { friend class cSysCo; public : - typedef cIsometry3D tPoseR; virtual ~cSysCoLEuc(); // do not copy and move because of PJ* (not needed via tPtrSysCo) @@ -164,7 +161,7 @@ public : tPt Value(const tPt &) const override; //< to GeoC tPt Inverse(const tPt &) const override; //< from GeoC - const tPoseR& getTranfo2GeoC() const { return mTranfo2GeoC; } + virtual const tPoseR* getTranfo2GeoC() const override { return &mTranfo2GeoC; } protected: cSysCoLEuc(const std::string & def, bool aDebug); //< construct from a definition, starting with LEuc cSysCoLEuc(bool aDebug); //< default constructor for derived classes @@ -257,9 +254,6 @@ tREAL8 cSysCo::getRadiusApprox(const tPt &in) const PJ_COORD pj_geog = proj_trans(mPJ_GeoC2Geog, PJ_FWD, pj_geoc); tREAL8 lat = pj_geog.lp.phi/AngleInRad(eTyUnitAngle::eUA_degree); - //GRS80 - const tREAL8 semi_axis = 6378137; - const tREAL8 e2 = 0.00669438; return semi_axis*sqrt(1-e2)/(1-e2*Square(sin(lat)));//total curvature sphere } @@ -281,6 +275,28 @@ tREAL8 cSysCo::getDistHzApprox(const tPt & aPtA, const tPt & aPtB) const return alpha*(radius + aPtAgeog.z()); } +tPt3dr cSysCo::toGeoG(const tPt & aPtIn) const +{ + auto inGeoc = Value(aPtIn); + PJ_COORD pj_geoc = toPjCoord(inGeoc); + PJ_COORD pj_geog = proj_trans(mPJ_GeoC2Geog, PJ_FWD, pj_geoc); + return fromPjCoord(pj_geog); +} + +tPt3dr cSysCo::fromGeoG(const tPt &aPtInGeoG) const +{ + PJ_COORD pj_geog = toPjCoord(aPtInGeoG); + PJ_COORD pj_geoc = proj_trans(mPJ_GeoC2Geog, PJ_INV, pj_geog); + return Inverse(fromPjCoord(pj_geoc)); +} + +const tPoseR* cSysCo::getTranfo2GeoC() const +{ + MMVII_INTERNAL_ASSERT_User(false, eTyUEr::eSysCo, + std::string("Error: getTranfo2GeoC() not defined for SysCo type ") + E2Str(mType)); + return nullptr; +} + tPtrSysCo cSysCo::MakeSysCo(const std::string &aDef, bool aDebug) { if (starts_with(aDef,MMVII_SysCoLocal)) @@ -405,12 +421,12 @@ cSysCoLEuc::~cSysCoLEuc() tPt3dr cSysCoLEuc::Value(const tPt & in) const //< to GeoC { - return getTranfo2GeoC().Rot().Mat() * in + getTranfo2GeoC().Tr(); + return getTranfo2GeoC()->Rot().Mat() * in + getTranfo2GeoC()->Tr(); } tPt3dr cSysCoLEuc::Inverse(const tPt & in) const //< from GeoC { - return getTranfo2GeoC().Rot().Mat().Transpose() * (in - getTranfo2GeoC().Tr()); + return getTranfo2GeoC()->Rot().Mat().Transpose() * (in - getTranfo2GeoC()->Tr()); } @@ -504,7 +520,7 @@ cRotation3D cSysCoRTL::getRot2Vertical(const tPt & aPtIn) const auto anOtherRTL = cSysCoLEuc::makeRTL(ptGeoC, MMVII_SysCoDefGeoC); auto anOtherRTL_asRTL = static_cast(anOtherRTL.get()); // TODO: add vertical deflection - return cRotation3D(anOtherRTL_asRTL->getTranfo2GeoC().Rot().Mat().Transpose(),false) * getTranfo2GeoC().Rot(); + return cRotation3D(anOtherRTL_asRTL->getTranfo2GeoC()->Rot().Mat().Transpose(),false)* getTranfo2GeoC()->Rot(); } //------------------------------------------------------------ @@ -624,16 +640,16 @@ void BenchSysCo(cParamExeBench & aParam) // 0 0.6580976792477065 0.7529325630949846 4779236.016271434 auto aPose = aSysCoRTL_asRTL->getTranfo2GeoC(); - MMVII_INTERNAL_ASSERT_bench(Norm2(aPose.Tr()-tPt3dr(4201661.926785135,177860.1878016033,4779236.016271434))<0.00001,"SysCo RTL"); - MMVII_INTERNAL_ASSERT_bench(Norm2(aPose.Rot().AxeI()-tPt3dr(-0.042293037933441094,0.9991052491816668,0.))<0.00001,"SysCo RTL"); - MMVII_INTERNAL_ASSERT_bench(Norm2(aPose.Rot().AxeJ()-tPt3dr(-0.7522588760680056,-0.031843805452299215,0.6580976792477065))<0.00001,"SysCo RTL"); - MMVII_INTERNAL_ASSERT_bench(Norm2(aPose.Rot().AxeK()-tPt3dr(0.6575088458106564,0.027832950112332798,0.7529325630949846))<0.00001,"SysCo RTL"); + MMVII_INTERNAL_ASSERT_bench(Norm2(aPose->Tr()-tPt3dr(4201661.926785135,177860.1878016033,4779236.016271434))<0.00001,"SysCo RTL"); + MMVII_INTERNAL_ASSERT_bench(Norm2(aPose->Rot().AxeI()-tPt3dr(-0.042293037933441094,0.9991052491816668,0.))<0.00001,"SysCo RTL"); + MMVII_INTERNAL_ASSERT_bench(Norm2(aPose->Rot().AxeJ()-tPt3dr(-0.7522588760680056,-0.031843805452299215,0.6580976792477065))<0.00001,"SysCo RTL"); + MMVII_INTERNAL_ASSERT_bench(Norm2(aPose->Rot().AxeK()-tPt3dr(0.6575088458106564,0.027832950112332798,0.7529325630949846))<0.00001,"SysCo RTL"); // RTL to GeoC tPtrSysCo aSysCoGeoC = cSysCo::MakeSysCo("GeoC"); cChangeSysCo aRTL2GeoC(aSysCoRTL, aSysCoGeoC); - tPt3dr aPtGeoC = aSysCoRTL_asRTL->getTranfo2GeoC().Tr(); + tPt3dr aPtGeoC = aSysCoRTL_asRTL->getTranfo2GeoC()->Tr(); MMVII_INTERNAL_ASSERT_bench(Norm2(aRTL2GeoC.Inverse(aPtGeoC)-tPt3dr(0.,0.,0.))<0.00001,"SysCo RTL2GeoC"); tPt3dr aPtRTL = {100.,10,1.};