Repository: openmrslab/suspect

Issue:
ENH: to/from_scanner() should accept arrays of coords
Currently the to_scanner() and from_scanner() functions take three arguments: x, y, z and return an numpy array. It would be convenient if the functions could handle a numpy array as an input for consistency, this would also allow passing in an array of coords to be transformed all at once, for example from meshgrid(). On the other hand, in the case where the user wants to pass in simple known x, y and z values, this requires them to be wrapped in a tuple, being called like to_scanner((x, y, z)), similarly to the numpy.ones() function. This could be addressed by permitting either 1 or 3 *args and interpreting either as an ndarray of shape [...3] or 3 floats.

Discussion:
coveralls: 
[![Coverage Status](https://coveralls.io/builds/12233612/badge)](https://coveralls.io/builds/12233612)

Coverage increased (+0.8%) to 78.486% when pulling **c1e331e468074797918589bae34635e5166f596b on 63_array_transforms** into **e63218387349081ee3d790ee96d1cb8e0ec74af2 on master**.


Patch:
diff --git a/suspect/_transforms.py b/suspect/_transforms.py
index caddbfe..ee2f9f2 100644
--- a/suspect/_transforms.py
+++ b/suspect/_transforms.py
@@ -68,3 +68,37 @@ def rotation_matrix(angle, axis):
     matrix[2, 1] = axis[2] * axis[1] * (1 - c) + axis[0] * s
     matrix[2, 2] = c + axis[2] ** 2 * (1 - c)
     return matrix
+
+
+def normalise_positions_for_transform(*args):
+    """
+    Takes an input set of arguments which should represent some (x, y, z)
+    coords to be transformed and makes sure they are in a numpy.ndarray,
+    and adds a w dimension of magnitude 1.
+    
+    The two acceptable forms of input arguments are a single array_like
+    with final dimension of 3, or three floating point arguments representing
+    x, y and z. In the first case the returned array will have the same shape
+    for all axes except the last, which will go from size 3 to 4, while in the
+    second case the returned array will be of shape (4,).
+    
+    Parameters
+    ----------
+    args : array_like or 3 separate floats
+        The arguments to be processed
+    
+    Returns
+    -------
+    numpy.ndarray
+        Points ready for transformation by a matrix
+    """
+    if len(args) == 3:
+        positions = [*args, 1]
+    elif len(args) == 1:
+        positions = numpy.atleast_2d(args[0])
+        w_array = numpy.expand_dims(numpy.ones(positions.shape[:-1]), axis=-1)
+        positions = numpy.append(positions, w_array, axis=-1)
+    else:
+        raise ValueError("Unrecognised form for input args")
+
+    return positions
diff --git a/suspect/base.py b/suspect/base.py
index 02025a5..a7eeb60 100644
--- a/suspect/base.py
+++ b/suspect/base.py
@@ -1,5 +1,7 @@
 import numpy
 
+from . import _transforms
+
 
 class ImageBase(numpy.ndarray):
     """
@@ -15,10 +17,14 @@ class ImageBase(numpy.ndarray):
     def __array_finalize__(self, obj):
         self.transform = getattr(obj, 'transform', None)
 
-    def to_scanner(self, x, y, z):
+    def to_scanner(self, *args):
         """
         Converts a 3d position in ImageBase space to the scanner
-        reference space. Raises a ValueError if no transform is set.
+        reference space. Argument can either be 3 individual floats
+        for x, y and z or a numpy array_like with final dimension of
+        size 3.
+        
+        Raises a ValueError if no transform is set.
 
         Parameters
         ----------
@@ -37,14 +43,20 @@ class ImageBase(numpy.ndarray):
         if self.transform is None:
             raise ValueError("No transform set for {} object {}".format(type(self), self))
 
-        transformed_point = self.transform * numpy.matrix([x, y, z, 1]).T
+        positions = _transforms.normalise_positions_for_transform(*args)
+
+        transformed_point = numpy.einsum("ij,...j", self.transform, positions)
 
-        return numpy.squeeze(numpy.asarray(transformed_point))[0:3]
+        return numpy.squeeze(numpy.asarray(transformed_point))[..., 0:3]
 
-    def from_scanner(self, x, y, z):
+    def from_scanner(self, *args):
         """
         Converts a 3d position in scanner space to the ImageBase
-        reference space. Raises a ValueError if no transform is set.
+        reference space. Argument can either be 3 individual floats
+        for x, y and z or a numpy array_like with final dimension of
+        size 3.
+        
+        Raises a ValueError if no transform is set.
 
         Parameters
         ----------
@@ -63,9 +75,13 @@ class ImageBase(numpy.ndarray):
         if self.transform is None:
             raise ValueError("No transform set for {} object {}".format(type(self), self))
 
-        transformed_point = numpy.linalg.inv(self.transform) * numpy.matrix([x, y, z, 1]).T
+        positions = _transforms.normalise_positions_for_transform(*args)
+
+        transformed_point = numpy.einsum("ij,...j",
+                                         numpy.linalg.inv(self.transform),
+                                         positions)
 
-        return numpy.squeeze(numpy.asarray(transformed_point))[0:3]
+        return numpy.squeeze(numpy.asarray(transformed_point))[..., 0:3]
 
     @property
     def voxel_size(self):


Test patch:
diff --git a/tests/test_mrs/test_base.py b/tests/test_mrs/test_base.py
index 58cee68..9d35987 100644
--- a/tests/test_mrs/test_base.py
+++ b/tests/test_mrs/test_base.py
@@ -1,6 +1,8 @@
 import suspect
 import numpy
+import pytest
 import suspect.base
+from suspect import _transforms
 
 
 def test_create_base():
@@ -12,9 +14,29 @@ def test_create_base():
 
 
 def test_base_transform():
-    position = [10, 20, 30]
-    voxel_size = [20, 20, 20]
+    position = numpy.array([10, 20, 30])
+    voxel_size = numpy.array([20, 20, 20])
     transform = suspect.transformation_matrix([1, 0, 0], [0, 1, 0], position, voxel_size)
     base = suspect.base.ImageBase(numpy.zeros(1), transform)
     numpy.testing.assert_equal(base.position, position)
     numpy.testing.assert_equal(base.voxel_size, voxel_size)
+    transformed = base.to_scanner(0, 0, 0)
+    numpy.testing.assert_equal(transformed, position)
+    transformed = base.to_scanner(numpy.array([[0, 0, 0], [1, 1, 1]]))
+    numpy.testing.assert_equal(transformed, [position, position + voxel_size])
+    transformed = base.from_scanner(position)
+    numpy.testing.assert_equal((0, 0, 0), transformed)
+
+
+def test_transforms_fail():
+    base = suspect.base.ImageBase(numpy.zeros(1))
+    with pytest.raises(ValueError):
+        base.to_scanner(0, 0, 0)
+    with pytest.raises(ValueError):
+        base.from_scanner(0, 0, 0)
+    with pytest.raises(ValueError):
+        pos = base.position
+    with pytest.raises(ValueError):
+        vox = base.voxel_size
+    with pytest.raises(ValueError):
+        _transforms.normalise_positions_for_transform(0, 0)