From 3b13f91f5779f7882e8960e1eb07b6145a8d219e Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sat, 13 Jan 2024 22:04:00 +0100
Subject: [PATCH] OGRProjCT::TransformWithErrorCodes(): Improve performance of
 axis swapping (inspired by @SunBlack proposal) (fixes #9073)

---
 ogr/ogrct.cpp | 101 +++++++++++++++++++++++++++++++++-----------------
 1 file changed, 67 insertions(+), 34 deletions(-)

diff --git a/ogr/ogrct.cpp b/ogr/ogrct.cpp
index 782fa6943f25..2fe411a29d30 100644
--- a/ogr/ogrct.cpp
+++ b/ogr/ogrct.cpp
@@ -2308,25 +2308,43 @@ int OGRProjCT::TransformWithErrorCodes(size_t nCount, double *x, double *y,
     /* -------------------------------------------------------------------- */
     /*      Apply data axis to source CRS mapping.                          */
     /* -------------------------------------------------------------------- */
+
+    // Since we may swap the x and y pointers, but cannot tell the caller about this swap,
+    // we save the original pointer. The same axis swap code is executed for poSRSTarget.
+    // If this nullifies, we save the swap of both axes
+    const auto xOriginal = x;
+
     if (poSRSSource)
     {
         const auto &mapping = poSRSSource->GetDataAxisToSRSAxisMapping();
-        if (mapping.size() >= 2 && (mapping[0] != 1 || mapping[1] != 2))
+        if (mapping.size() >= 2)
         {
-            for (size_t i = 0; i < nCount; i++)
+            if (std::abs(mapping[0]) == 2 && std::abs(mapping[1]) == 1)
+            {
+                std::swap(x, y);
+            }
+            const bool bNegateX = mapping[0] < 0;
+            if (bNegateX)
+            {
+                for (size_t i = 0; i < nCount; i++)
+                {
+                    x[i] = -x[i];
+                }
+            }
+            const bool bNegateY = mapping[1] < 0;
+            if (bNegateY)
+            {
+                for (size_t i = 0; i < nCount; i++)
+                {
+                    y[i] = -y[i];
+                }
+            }
+            if (z && mapping.size() >= 3 && mapping[2] == -3)
             {
-                double newX = (mapping[0] == 1)    ? x[i]
-                              : (mapping[0] == -1) ? -x[i]
-                              : (mapping[0] == 2)  ? y[i]
-                                                   : -y[i];
-                double newY = (mapping[1] == 2)    ? y[i]
-                              : (mapping[1] == -2) ? -y[i]
-                              : (mapping[1] == 1)  ? x[i]
-                                                   : -x[i];
-                x[i] = newX;
-                y[i] = newY;
-                if (z && mapping.size() >= 3 && mapping[2] == -3)
+                for (size_t i = 0; i < nCount; i++)
+                {
                     z[i] = -z[i];
+                }
             }
         }
     }
@@ -2374,10 +2392,7 @@ int OGRProjCT::TransformWithErrorCodes(size_t nCount, double *x, double *y,
 
         if (m_eSourceFirstAxisOrient != OAO_East)
         {
-            for (size_t i = 0; i < nCount; i++)
-            {
-                std::swap(x[i], y[i]);
-            }
+            std::swap(x, y);
         }
 
         double y0 = y[0];
@@ -2458,10 +2473,7 @@ int OGRProjCT::TransformWithErrorCodes(size_t nCount, double *x, double *y,
 
         if (m_eTargetFirstAxisOrient != OAO_East)
         {
-            for (size_t i = 0; i < nCount; i++)
-            {
-                std::swap(x[i], y[i]);
-            }
+            std::swap(x, y);
         }
 
         bTransformDone = true;
@@ -2830,26 +2842,47 @@ int OGRProjCT::TransformWithErrorCodes(size_t nCount, double *x, double *y,
     if (poSRSTarget)
     {
         const auto &mapping = poSRSTarget->GetDataAxisToSRSAxisMapping();
-        if (mapping.size() >= 2 && (mapping[0] != 1 || mapping[1] != 2))
+        if (mapping.size() >= 2)
         {
-            for (size_t i = 0; i < nCount; i++)
+            if (std::abs(mapping[0]) == 2 && std::abs(mapping[1]) == 1)
+            {
+                std::swap(x, y);
+            }
+            const bool bNegateX = mapping[0] < 0;
+            if (bNegateX)
+            {
+                for (size_t i = 0; i < nCount; i++)
+                {
+                    x[i] = -x[i];
+                }
+            }
+            const bool bNegateY = mapping[1] < 0;
+            if (bNegateY)
+            {
+                for (size_t i = 0; i < nCount; i++)
+                {
+                    y[i] = -y[i];
+                }
+            }
+            if (z && mapping.size() >= 3 && mapping[2] == -3)
             {
-                double newX = (mapping[0] == 1)    ? x[i]
-                              : (mapping[0] == -1) ? -x[i]
-                              : (mapping[0] == 2)  ? y[i]
-                                                   : -y[i];
-                double newY = (mapping[1] == 2)    ? y[i]
-                              : (mapping[1] == -2) ? -y[i]
-                              : (mapping[1] == 1)  ? x[i]
-                                                   : -x[i];
-                x[i] = newX;
-                y[i] = newY;
-                if (z && mapping.size() >= 3 && mapping[2] == -3)
+                for (size_t i = 0; i < nCount; i++)
+                {
                     z[i] = -z[i];
+                }
             }
         }
     }
 
+    // Check whether final "genuine" axis swap is really necessary
+    if (x != xOriginal)
+    {
+        for (size_t i = 0; i < nCount; i++)
+        {
+            std::swap(x[i], y[i]);
+        }
+    }
+
 #ifdef DEBUG_VERBOSE
     if (bDebugCT)
     {