From a2941046b43b95299f298a0366d59b3dae55c325 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Denis=20Rivi=C3=A8re?= <denis.riviere.gm@gmail.com>
Date: Tue, 7 Feb 2023 14:42:15 +0100
Subject: [PATCH] mark if x is the first index when building a volume from an
 array

(#92)
---
 pyaims/src/libpyaims/vector/numpyarrayfunc.cc | 57 ++++++++++++++-----
 pyaims/src/sip/volume.tpl                     | 12 ++--
 2 files changed, 50 insertions(+), 19 deletions(-)

diff --git a/pyaims/src/libpyaims/vector/numpyarrayfunc.cc b/pyaims/src/libpyaims/vector/numpyarrayfunc.cc
index 8cd7d154..d6ae3f26 100644
--- a/pyaims/src/libpyaims/vector/numpyarrayfunc.cc
+++ b/pyaims/src/libpyaims/vector/numpyarrayfunc.cc
@@ -39,28 +39,38 @@ using namespace std;
 namespace
 {
 
-  PyObject* transposeNumpyArray( PyObject* a, bool xyzorder )
+  PyObject* transposeNumpyArray( PyObject* a, bool xyzorder,
+                                 int current_x_first = -1 )
   {
     PyArrayObject *arr = (PyArrayObject *) a;
     int i, ndim = PyArray_NDIM( arr );
     if( ndim < 2 )
       return a;
-    bool firstindexinc = xyzorder;
-    for( i=0; i<ndim - 1; ++i )
+    bool transpose = false;
+    if( current_x_first >= 0 )  // order is given
+      transpose = ( bool( current_x_first ) != xyzorder );
+    else
     {
-      if( PyArray_STRIDES( arr )[i] < PyArray_STRIDES( arr )[i + 1] )
-      {
-        firstindexinc = true;
-        break;
-      }
-      else if( PyArray_STRIDES( arr )[i] > PyArray_STRIDES( arr )[i + 1] )
+      bool firstindexinc = xyzorder;
+      for( i=0; i<ndim - 1; ++i )
       {
-        firstindexinc = false;
-        break;
+        if( PyArray_STRIDES( arr )[i] < PyArray_STRIDES( arr )[i + 1] )
+        {
+          firstindexinc = true;
+          break;
+        }
+        else if( PyArray_STRIDES( arr )[i] > PyArray_STRIDES( arr )[i + 1] )
+        {
+          firstindexinc = false;
+          break;
+        }
+        // otherwise equal, check next dim
       }
-      // otherwise equal, check next dim
+      if( xyzorder != firstindexinc )
+        transpose = true;
     }
-    if( xyzorder == firstindexinc )
+
+    if( !transpose )
       return a;
     PyArray_Dims  adims;
     vector<npy_intp> dims( ndim );
@@ -91,7 +101,16 @@ namespace aims
     {
       PyObject *arr = PyObject_GetAttrString( sipSelf, "_arrayext" );
       if( arr )
-        sipRes = transposeNumpyArray( arr, xyzorder );
+      {
+        int current_x_first = -1;
+        if( PyObject_HasAttrString( sipSelf, "_array_x_is_first" ) )
+        {
+          PyObject *x = PyObject_GetAttrString( sipSelf, "_array_x_is_first" );
+          if( x )
+            current_x_first = int( bool( PyLong_AsLong( x ) ) );
+        }
+        sipRes = transposeNumpyArray( arr, xyzorder, current_x_first );
+      }
     }
     // else look if an array has already been built on the object
     else if( PyObject_HasAttrString( sipSelf, "_arrayref" ) )
@@ -103,7 +122,15 @@ namespace aims
         if( arr )
         {
           Py_INCREF( arr );
-          sipRes = transposeNumpyArray( arr, xyzorder );
+          int current_x_first = -1;
+          if( PyObject_HasAttrString( sipSelf, "_array_x_is_first" ) )
+          {
+            PyObject *x = PyObject_GetAttrString( sipSelf,
+                                                  "_array_x_is_first" );
+            if( x )
+              current_x_first = int( bool( PyLong_AsLong( x ) ) );
+          }
+          sipRes = transposeNumpyArray( arr, xyzorder, current_x_first );
         }
         else
         {
diff --git a/pyaims/src/sip/volume.tpl b/pyaims/src/sip/volume.tpl
index 56e69b82..8fddfaef 100644
--- a/pyaims/src/sip/volume.tpl
+++ b/pyaims/src/sip/volume.tpl
@@ -715,20 +715,20 @@ The header contains all meta-data.
     int nd = PyArray_NDIM( arr );
     if( nd < 4 )
       nd = 4;
+    bool xyzorder = false;
     std::vector<int> dims( nd, 1 );
     int inc = 1, start = 0;
-    /*
     // TODO: retreive exact strides and react accordingly
     if( PyArray_NDIM( arr ) >= 2 )
     {
       if( PyArray_STRIDES( arr )[1] < PyArray_STRIDES( arr )[0] )
       {
+        xyzorder = true;
         // increments higher index first
-        inc = -1;
-        start = PyArray_NDIM( arr )-1;
+        //  inc = -1;
+        // start = PyArray_NDIM( arr )-1;
       }
     }
-    */
 
     std::vector<size_t> strides( nd, 1 );
 
@@ -748,6 +748,10 @@ The header contains all meta-data.
        is still living.
     */
     PyObject_SetAttrString( (PyObject *) sipSelf, "_arrayext", a0 );
+    PyObject *xfirst = PyBool_FromLong( long( xyzorder ) );
+    Py_DECREF( xfirst );
+    PyObject_SetAttrString( (PyObject *) sipSelf, "_array_x_is_first",
+                            xfirst );
   }
   else
   {