Extends SIRF data algebra possible operands to numpy arrays (#1358)

* added numpy array to possible SIRF image data algebra operands (on the "right")

* added to User Guide a subsection on SIRF/numpy data algebra peculiarities

* updated CHANGES.md

* added tests

Author: evgueni-ovtchinnikov

CHANGES.md

@@ -1,5 +1,10 @@
 # ChangeLog
 
+## v3.9.1
+
+* Python interface
+  - Restored functionality for algebraic operations mixing STIR.ImageData and numpy arrays. (Note that sirf objects need to be on the "left" of the operation.)
+
 ## v3.9.0
 
 * Python interface

doc/UserGuide.md

@@ -108,6 +108,10 @@ The mutators are also responsible for basic error checking.
 
 Some classes are _derived_ from other classes, which means that they have (_inherit_) all the methods of the classes they are derived from. If class B derives from class A, then A is called its _base_ class. <!-- we say that these _derived_ class methods are _inherited_ from the _base_ class.--> For example, class `AcquisitionModelUsingRayTracingMatrix` is derived from `AcquisitionModelUsingMatrix`, which in turn is derived from `AcquisitionModel`, and so it inherits all the methods of the latter two base classes.
 
+### SIRF data algebra <a name="SIRF_data_algebra"></a>
+
+SIRF Python interface supports algebraic operations (`+`, `-`, `*` and `/`): e.g. elements of the data array stored in the object `a*b` are the products of the respective elements in `a` and `b`. Either or both `a` and `b` can be SIRF data objects of the same kind (either both `ImageData` or both `AcquisitionData`) or `numpy` arrays or scalars. One should be aware though that if `a` is a SIRF object then, just as one would expect, the product `a*b` will be a SIRF object of the same kind, but if `a` is a `numpy` object (array or scalar) then Python will try to convert `b` to a `numpy` object before computing `a*b`, and only if this fails it will compute `b*a` instead. To avoid confusion, the users are advised to check the type of `a*b` or, better still, always place a SIRF object on the left side of the product.
+
 ### Error handling <a name="error_handling"></a>
 
 Error handling is via exceptions, i.e. functions do not return an error status, but throw an error if something did not work. The user can catch these exceptions if required as illustrated in the demos. 

src/Registration/pReg/Reg.py

@@ -32,6 +32,7 @@
 
 import sirf.Reg_params as parms
 import numpy
+from numbers import Number
 from sirf.config import SIRF_HAS_SPM
 
 if sys.version_info[0] >= 3 and sys.version_info[1] >= 4:
@@ -264,6 +265,64 @@ def get_voxel_sizes(self):
             self.handle, out.ctypes.data))
         return out
 
+    def dot(self, other):
+        '''
+        Returns the dot product of the container data with another container
+        data or numpy array viewed as vectors.
+        other: NiftiImageData or numpy array.
+        '''
+        if not (issubclass(type(other), type(self))):
+            self_copy = self.clone()
+            self_copy.fill(other)
+            other = self_copy
+        return super().dot(other)
+
+    def add(self, other, out=None):
+        '''
+        Addition for NiftiImageData containers.
+
+        If other is a NiftiData or numpy array, returns the sum of data
+        stored in self and other viewed as vectors.
+        If other is a scalar, returns the same with the second vector filled
+        with the value of other.
+        other: NiftiImageData or numpy array or scalar.
+        out:   NiftiImageData to store the result to.
+        '''
+        if not (issubclass(type(other), type(self)) or isinstance(other, (Number, numpy.number))):
+            self_copy = self.clone()
+            self_copy.fill(other)
+            other = self_copy
+        return super().add(other, out)
+
+    def subtract(self, other, out=None):
+        '''
+        Subtraction for NiftiImageData containers.
+
+        If other is a NiftiData or numpy array, returns the sum of data
+        stored in self and other viewed as vectors.
+        If other is a scalar, returns the same with the second vector filled
+        with the value of other.
+        other: NiftiImageData or numpy array or scalar.
+        out:   NiftiImageData to store the result to.
+        '''
+        if not (issubclass(type(other), type(self)) or isinstance(other, (Number, numpy.number))):
+            self_copy = self.clone()
+            self_copy.fill(other)
+            other = self_copy
+        return super().subtract(other, out)
+
+    def binary(self, other, f, out=None):
+        '''Applies function f(x,y) element-wise to self and other.
+
+        other: NiftiImageData or numpy array or Number
+        f: the name of the function to apply, Python str.
+        '''
+        if not (issubclass(type(other), type(self)) or isinstance(other, (Number, numpy.number))):
+            self_copy = self.clone()
+            self_copy.fill(other)
+            other = self_copy
+        return super().binary(other, f, out)
+
     def fill(self, val):
         """Fill image with single value or numpy array."""
         if self.handle is None:

src/common/SIRF.py

@@ -139,9 +139,7 @@ def __sub__(self, other):
         data viewed as vectors.
         other: DataContainer
         '''
-        if isinstance(other, (DataContainer, Number) ):
-            return self.subtract(other)
-        return NotImplemented
+        return self.subtract(other)
 
     def __mul__(self, other):
         '''
@@ -331,6 +329,8 @@ def dot(self, other):
         data viewed as vectors.
         other: DataContainer
         '''
+        if not (issubclass(type(other), type(self))):
+            other = self.clone().fill(other)
         assert_validities(self, other)
         # Check if input are the same size
         if self.size != other.size:
@@ -383,6 +383,8 @@ def add(self, other, out=None):
         other: DataContainer or scalar.
         out:   DataContainer to store the result to.
         '''
+        if not (issubclass(type(other), type(self)) or isinstance(other, (Number, numpy.number))):
+            other = self.clone().fill(other)
         if out is None:
             z = self.same_object()
         else:
@@ -416,10 +418,14 @@ def subtract(self, other, out=None):
         other: DataContainer or scalar.
         other: DataContainer
         '''
+        if not (issubclass(type(other), type(self)) or isinstance(other, (Number, numpy.number))):
+            other = self.clone().fill(other)
         if not isinstance (other, (DataContainer, Number)):
             return NotImplemented
         if isinstance(other, Number):
             return self.add(-other, out=out)
+        if not (issubclass(type(other), type(self)) or isinstance(other, (Number, numpy.number))):
+            other = self.clone().fill(other)
         assert_validities(self, other)
         pl_one = numpy.asarray([1.0, 0.0], dtype = numpy.float32)
         mn_one = numpy.asarray([-1.0, 0.0], dtype = numpy.float32)
@@ -625,6 +631,8 @@ def binary(self, other, f, out=None):
                 try_calling(pysirf.cSIRF_compute_semibinary(self.handle, y.ctypes.data, \
                                                         f, out.handle))
         else:
+            if not (issubclass(type(other), type(self)) or isinstance(other, (Number, numpy.number))):
+                other = self.clone().fill(other)
             assert_validities(self, other)
             if out.handle is None:
                 out.handle = pysirf.cSIRF_binary(self.handle, other.handle, f)

src/common/Utilities.py

@@ -742,6 +742,11 @@ def data_container_algebra_tests(test, x, eps=1e-4):
     # needs increased tolerance for large data size
     test.check_if_equal_within_tolerance(t, s, 0, eps * 10);
 
+    s = x.dot(ay)
+    t = numpy.vdot(ay, ax)
+    # needs increased tolerance for large data size
+    test.check_if_equal_within_tolerance(t, s, 0, eps * 10);
+
     x2 = x.multiply(2)
     ax2 = x2.as_array()
     s = numpy.linalg.norm(ax2 - 2*ax)
@@ -779,12 +784,30 @@ def data_container_algebra_tests(test, x, eps=1e-4):
     t = numpy.linalg.norm(az)
     test.check_if_zero_within_tolerance(s, eps * t)
 
+    z = x*ay
+    az = z.as_array()
+    s = numpy.linalg.norm(az - ax * ay)
+    t = numpy.linalg.norm(az)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
     y = x + 1
     ay = y.as_array()
     s = numpy.linalg.norm(ay - (ax + 1))
     t = numpy.linalg.norm(ay)
     test.check_if_zero_within_tolerance(s, eps * t)
 
+    y = x + ax
+    ay = y.as_array()
+    s = numpy.linalg.norm(ay - (ax + ax))
+    t = numpy.linalg.norm(ay)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
+    t = numpy.linalg.norm(ay)
+    y = x - ax
+    ay = y.as_array()
+    s = numpy.linalg.norm(ay)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
     y *= 0
     x.add(1, out=y)
     ay = y.as_array()
@@ -798,6 +821,12 @@ def data_container_algebra_tests(test, x, eps=1e-4):
     t = numpy.linalg.norm(az)
     test.check_if_zero_within_tolerance(s, eps * t)
 
+    z = x/ay
+    az = z.as_array()
+    s = numpy.linalg.norm(az - ax/ay)
+    t = numpy.linalg.norm(az)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
     z = x/2
     az = z.as_array()
     s = numpy.linalg.norm(az - ax/2)