Fixed sirf.Reg.ImageData.asarray().

CHANGES.md

@@ -5,8 +5,11 @@
 * SIRF/STIR
   - The implementation of the creation of `sirf.STIR.ImageData` from `sirf.STIR.AcquisitionData` has been revised to ensure compatibility of `ImageData` dimensions and voxel sizes with `AcquisitionData`.
   - Missing `__del__` added to sirf.STIR.AcquisitionModel.
+* SIRF/Registration
+  - sirf.Reg.ImageData stores voxels values in a 3D Fortran-style array, hence strides had to be set in `ImageData.__array_interface__` accordingly (other SIRF ImageData objects use default C-style).
 * Python interface
   - Restored functionality for algebraic operations mixing SIRF data containers and numpy arrays and corrected the description of the result type in User Guide.
+  - Adjoint operator added in `common/SIRF.py`.
   - Error raised if `AcquisitionModel.adjoint` ran when the model is not linear.
 * SIRF
   - `common/SIRF.py` adding adjoint operator

src/Registration/pReg/Reg.py

@@ -488,9 +488,17 @@ def __array_interface__(self):
         """As per https://numpy.org/doc/stable/reference/arrays.interface.html"""
         if not self.supports_array_view:
             raise ContiguousError("please make an array-copy first with `asarray(copy=True)` or `as_array()`")
-        return {
-            'shape': self.shape, 'typestr': '<f4', 'version': 3,
-            'data': (parms.size_t_par(self.handle, 'NiftiImageData', 'address'), False)}
+        # NiftiImageData stores voxels values in a 3D Fortran-style array,
+        # hence strides need to be arranged accordingly (default is C-style).
+        shape = self.shape
+        dims = len(shape)
+        strides = ()
+        stride = 4
+        for i in range(dims):
+            strides += (stride,)
+            stride *= shape[i]
+        return {'shape': self.shape, 'typestr': '<f4', 'strides': strides, 'version': 3,
+                'data': (parms.size_t_par(self.handle, 'NiftiImageData', 'address'), False)}
 
 
 class NiftiImageData3D(NiftiImageData):

src/common/Utilities.py

@@ -731,6 +731,7 @@ def is_operator_adjoint(operator, num_tests=5, max_err=10e-5, verbose=True):
 def data_container_algebra_tests(test, x, eps=1e-4):
 
     ax = x.as_array()
+
     ay = numpy.ones_like(ax)
     y = x.clone()
     y.fill(ay)