another improvement in isEquivalentTo() when comparing Helmert transformations

2 files changed, 57 insertions, 46 deletions

diff --git a/src/coordinateoperation.cpp b/src/coordinateoperation.cpp
index 90854d2e..a1f305bf 100644
--- a/src/coordinateoperation.cpp
+++ b/src/coordinateoperation.cpp
@@ -1510,59 +1510,44 @@ bool SingleOperation::_isEquivalentTo(const util::IComparable *other,
 
     if (!equivalentMethods) {
         if (criterion == util::IComparable::Criterion::EQUIVALENT) {
-            // _1SP methods can sometimes be equivalent to _2SP ones
-            // Check it by using convertToOtherMethod()
 
-            if ((methodEPSGCode ==
-                     EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOCENTRIC &&
-                 (otherMethodEPSGCode ==
-                      EPSG_CODE_METHOD_POSITION_VECTOR_GEOCENTRIC ||
-                  otherMethodEPSGCode ==
-                      EPSG_CODE_METHOD_COORDINATE_FRAME_GEOCENTRIC)) ||
-                (methodEPSGCode ==
-                     EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D &&
-                 (otherMethodEPSGCode ==
-                      EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
-                  otherMethodEPSGCode ==
-                      EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D)) ||
-                (methodEPSGCode ==
-                     EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_3D &&
-                 (otherMethodEPSGCode ==
-                      EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_3D ||
-                  otherMethodEPSGCode ==
-                      EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_3D))) {
+            const auto isTOWGS84Transf = [](int code) {
+                return code ==
+                           EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOCENTRIC ||
+                       code == EPSG_CODE_METHOD_POSITION_VECTOR_GEOCENTRIC ||
+                       code == EPSG_CODE_METHOD_COORDINATE_FRAME_GEOCENTRIC ||
+                       code ==
+                           EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D ||
+                       code == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
+                       code ==
+                           EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D ||
+                       code ==
+                           EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_3D ||
+                       code == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_3D ||
+                       code == EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_3D;
+            };
+
+            // Translation vs (PV or CF)
+            // or different PV vs CF convention
+            if (isTOWGS84Transf(methodEPSGCode) &&
+                isTOWGS84Transf(otherMethodEPSGCode)) {
                 auto transf = static_cast<const Transformation *>(this);
                 auto otherTransf = static_cast<const Transformation *>(otherSO);
                 auto params = transf->getTOWGS84Parameters();
                 auto otherParams = otherTransf->getTOWGS84Parameters();
-                return params == otherParams;
-            }
-
-            if ((otherMethodEPSGCode ==
-                     EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOCENTRIC &&
-                 (methodEPSGCode ==
-                      EPSG_CODE_METHOD_POSITION_VECTOR_GEOCENTRIC ||
-                  methodEPSGCode ==
-                      EPSG_CODE_METHOD_COORDINATE_FRAME_GEOCENTRIC)) ||
-                (otherMethodEPSGCode ==
-                     EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D &&
-                 (methodEPSGCode ==
-                      EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
-                  methodEPSGCode ==
-                      EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D)) ||
-                (otherMethodEPSGCode ==
-                     EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_3D &&
-                 (methodEPSGCode ==
-                      EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_3D ||
-                  methodEPSGCode ==
-                      EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_3D))) {
-                auto transf = static_cast<const Transformation *>(this);
-                auto otherTransf = static_cast<const Transformation *>(otherSO);
-                auto params = transf->getTOWGS84Parameters();
-                auto otherParams = otherTransf->getTOWGS84Parameters();
-                return params == otherParams;
+                assert(params.size() == 7);
+                assert(otherParams.size() == 7);
+                for (size_t i = 0; i < 7; i++) {
+                    if (std::fabs(params[i] - otherParams[i]) >
+                        1e-10 * std::fabs(params[i])) {
+                        return false;
+                    }
+                }
+                return true;
             }
 
+            // _1SP methods can sometimes be equivalent to _2SP ones
+            // Check it by using convertToOtherMethod()
             if (methodEPSGCode ==
                     EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_1SP &&
                 otherMethodEPSGCode ==
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 66284b59..2695153a 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -6693,6 +6693,32 @@ TEST(operation, three_param_equivalent_to_seven_param) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, position_vector_equivalent_coordinate_frame) {
+
+    auto pv = Transformation::createPositionVector(
+        PropertyMap(), GeographicCRS::EPSG_4269, GeographicCRS::EPSG_4326, 1.0,
+        2.0, 3.0, 4.0, 5.0, 6.0, 7.0, {});
+
+    auto cf = Transformation::createCoordinateFrameRotation(
+        PropertyMap(), GeographicCRS::EPSG_4269, GeographicCRS::EPSG_4326, 1.0,
+        2.0, 3.0, -4 + 1e-11, -5.0, -6.0, 7.0, {});
+
+    auto cf_non_eq = Transformation::createCoordinateFrameRotation(
+        PropertyMap(), GeographicCRS::EPSG_4269, GeographicCRS::EPSG_4326, 1.0,
+        2.0, 3.0, 4.0, 5.0, 6.0, 7.0, {});
+
+    EXPECT_TRUE(
+        pv->isEquivalentTo(cf.get(), IComparable::Criterion::EQUIVALENT));
+
+    EXPECT_TRUE(
+        cf->isEquivalentTo(pv.get(), IComparable::Criterion::EQUIVALENT));
+
+    EXPECT_FALSE(pv->isEquivalentTo(cf_non_eq.get(),
+                                    IComparable::Criterion::EQUIVALENT));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, conversion_missing_parameter) {
 
     auto wkt1 = "PROJCS[\"NAD83(CSRS98) / UTM zone 20N (deprecated)\","