Repository: openmrslab/suspect

Issue:
ENH: Resample image volumes
It should be possible to derive a new ImageBase from an existing ImageBase object by supplying the coordinate system, size and shape.

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

Coverage increased (+0.4%) to 80.505% when pulling **900415d5cca892040ec52ec50201054a2de4599a on 80_image_resample** into **f5a6a8e7875a9114bab75aba928028cb69d9665e on master**.


Patch:
diff --git a/suspect/_version.py b/suspect/_version.py
index e19434e..334b899 100644
--- a/suspect/_version.py
+++ b/suspect/_version.py
@@ -1,1 +1,1 @@
-__version__ = "0.3.3"
+__version__ = "0.3.4"
diff --git a/suspect/base.py b/suspect/base.py
index 3bdf704..6952aa2 100644
--- a/suspect/base.py
+++ b/suspect/base.py
@@ -1,5 +1,6 @@
-import numpy
+import numpy as np
 import functools
+import scipy.interpolate
 
 from . import _transforms
 
@@ -10,17 +11,19 @@ def requires_transform(func):
         if self.transform is None:
             raise ValueError("No transform set for {} object {}".format(type(self), self))
         return func(self, *args, **kwargs)
+
     return wrapper
 
 
-class ImageBase(numpy.ndarray):
+class ImageBase(np.ndarray):
     """
     numpy.ndarray subclass with an affine transform associated with it
     """
+
     def __new__(cls, input_array, transform=None):
         # input_array is an already formed ndarray
         # we want to make it our class type
-        obj = numpy.asarray(input_array).view(cls)
+        obj = np.asarray(input_array).view(cls)
         if transform is not None:
             obj.transform = transform.copy()
         return obj
@@ -54,9 +57,9 @@ class ImageBase(numpy.ndarray):
         """
         positions = _transforms.normalise_positions_for_transform(*args)
 
-        transformed_point = numpy.einsum("ij,...j", self.transform, positions)
+        transformed_point = np.einsum("ij,...j", self.transform, positions)
 
-        return numpy.squeeze(numpy.asarray(transformed_point))[..., 0:3]
+        return np.squeeze(np.asarray(transformed_point))[..., 0:3]
 
     @requires_transform
     def from_scanner(self, *args):
@@ -84,11 +87,11 @@ class ImageBase(numpy.ndarray):
         """
         positions = _transforms.normalise_positions_for_transform(*args)
 
-        transformed_point = numpy.einsum("ij,...j",
-                                         numpy.linalg.inv(self.transform),
-                                         positions)
+        transformed_point = np.einsum("ij,...j",
+                                      np.linalg.inv(self.transform),
+                                      positions)
 
-        return numpy.squeeze(numpy.asarray(transformed_point))[..., 0:3]
+        return np.squeeze(np.asarray(transformed_point))[..., 0:3]
 
     @property
     @requires_transform
@@ -101,7 +104,7 @@ class ImageBase(numpy.ndarray):
         numpy.ndarray
             The dimensions of a voxel along each axis.
         """
-        return numpy.linalg.norm(self.transform, axis=0)[0:3]
+        return np.linalg.norm(self.transform, axis=0)[0:3]
 
     @property
     @requires_transform
@@ -114,22 +117,22 @@ class ImageBase(numpy.ndarray):
     @property
     @requires_transform
     def slice_vector(self):
-        return self.transform[:3, 2] / numpy.linalg.norm(self.transform[:3, 2])
+        return self.transform[:3, 2] / np.linalg.norm(self.transform[:3, 2])
 
     @property
     @requires_transform
     def row_vector(self):
-        return self.transform[:3, 1] / numpy.linalg.norm(self.transform[:3, 1])
+        return self.transform[:3, 1] / np.linalg.norm(self.transform[:3, 1])
 
     @property
     @requires_transform
     def col_vector(self):
-        return self.transform[:3, 0] / numpy.linalg.norm(self.transform[:3, 0])
+        return self.transform[:3, 0] / np.linalg.norm(self.transform[:3, 0])
 
     @requires_transform
     def _closest_axis(self, target_axis):
-        overlap = numpy.abs(numpy.dot(target_axis, self.transform[:3, :3]))
-        return self.transform[:3, numpy.argmax(overlap)]
+        overlap = np.abs(np.dot(target_axis, self.transform[:3, :3]))
+        return self.transform[:3, np.argmax(overlap)]
 
     @property
     @requires_transform
@@ -147,7 +150,7 @@ class ImageBase(numpy.ndarray):
         """
         # dot the three candidate vectors with (0, 0, 1)
         best_axis = self._closest_axis((0, 0, 1))
-        norm_axis = best_axis / numpy.linalg.norm(best_axis)
+        norm_axis = best_axis / np.linalg.norm(best_axis)
         # work out if we need to reverse the direction
         return norm_axis if norm_axis[2] > 0 else -1 * norm_axis
 
@@ -167,7 +170,7 @@ class ImageBase(numpy.ndarray):
         """
         # dot the three candidate vectors with (0, 1, 0)
         best_axis = self._closest_axis((0, 1, 0))
-        norm_axis = best_axis / numpy.linalg.norm(best_axis)
+        norm_axis = best_axis / np.linalg.norm(best_axis)
         return norm_axis if norm_axis[1] > 0 else -1 * norm_axis
 
     @property
@@ -186,5 +189,92 @@ class ImageBase(numpy.ndarray):
         """
         # dot the three candidate vectors with (1, 0, 0)
         best_axis = self._closest_axis((1, 0, 0))
-        norm_axis = best_axis / numpy.linalg.norm(best_axis)
+        norm_axis = best_axis / np.linalg.norm(best_axis)
         return norm_axis if norm_axis[0] > 0 else -1 * norm_axis
+
+    @property
+    @requires_transform
+    def centre(self):
+        """
+        Returns the centre of the image volume in scanner coordinates.
+
+        Returns
+        -------
+        numpy.ndarray
+            The centre of the image volume
+        :return:
+        """
+        return self.to_scanner((np.array(self.shape[::-1]) - 1) / 2)
+
+    @requires_transform
+    def resample(self,
+                 row_vector,
+                 col_vector,
+                 shape,
+                 centre=(0, 0, 0),
+                 voxel_size=(1, 1, 1),
+                 method='linear'):
+        """
+        Create a new volume by resampling this one using a different coordinate
+        system.
+
+        Parameters
+        ----------
+        row_vector: array
+            Row direction vector for new volume
+        col_vector: array
+            Column direction vector for new volume
+        shape: array
+            The shape of the new volume, as slices, rows, columns
+        centre: array
+            The position of the centre of the new volume in scanner
+            coordinates, in mm
+        voxel_size: array
+            The size of each voxel in the new volume, in mm
+        method: str
+            The interpolation method to use - either "linear" or "nearest"
+
+        Returns
+        -------
+        suspect.base.ImageBase
+            The resampled volume
+        """
+        # make sure row_vector and col_vector are normalised
+        row_vector = np.asanyarray(row_vector) / np.linalg.norm(row_vector)
+        col_vector = np.asanyarray(col_vector) / np.linalg.norm(col_vector)
+
+        # mgrid produces 3D index grids for the x, y and z coords separately
+        II, JJ, KK = np.mgrid[0:shape[2],
+                              0:shape[1],
+                              0:shape[0]].astype(np.float)
+        # shift the indices from the corner to the centre
+        II -= (shape[2] - 1) / 2
+        JJ -= (shape[1] - 1) / 2
+        KK -= (shape[0] - 1) / 2
+        # scale the indices by the size of the voxel
+        II *= voxel_size[0]
+        JJ *= voxel_size[1]
+        KK *= voxel_size[2]
+
+        slice_vector = np.cross(row_vector, col_vector)
+
+        # combine the x, y and z indices with the row, col and slice vectors
+        # to get the spatial coordinates at each point in the new volume
+        space_coords = II[..., np.newaxis] * row_vector \
+                       + JJ[..., np.newaxis] * col_vector \
+                       + KK[..., np.newaxis] * slice_vector + centre
+
+        resampled = scipy.interpolate.interpn([np.arange(dim) for dim in self.shape],
+                                              self,
+                                              self.from_scanner(space_coords)[..., ::-1],
+                                              method=method,
+                                              bounds_error=False,
+                                              fill_value=0)
+
+        transform = _transforms.transformation_matrix(row_vector,
+                                                      col_vector,
+                                                      space_coords[0, 0, 0],
+                                                      voxel_size)
+
+        # we have to transpose the result to go from x, y, z to row, col, slice
+        return ImageBase(resampled.T, transform=transform)
diff --git a/suspect/mrsobjects.py b/suspect/mrsobjects.py
index 19dcf41..1c976fc 100644
--- a/suspect/mrsobjects.py
+++ b/suspect/mrsobjects.py
@@ -230,6 +230,13 @@ class MRSBase(suspect.base.ImageBase):
         return self.slice_hz(self.ppm_to_hertz(lower_bound),
                              self.ppm_to_hertz(upper_bound))
 
+    @property
+    @suspect.base.requires_transform
+    def centre(self):
+        # unlike the ImageBase class, the centre of spectroscopy volumes is
+        # already encoded in the position of the volume
+        return self.position
+
 
 class MRSData(MRSBase):
     """


Test patch:
diff --git a/tests/test_mrs/test_base.py b/tests/test_mrs/test_base.py
index 59a36a8..72a8d6d 100644
--- a/tests/test_mrs/test_base.py
+++ b/tests/test_mrs/test_base.py
@@ -67,4 +67,36 @@ def test_find_axes_reversed():
     base = suspect.base.ImageBase(np.zeros(1), transform=transform)
     np.testing.assert_equal(base.axial_vector, base.slice_vector)
     np.testing.assert_equal(base.coronal_vector, -base.row_vector)
-    np.testing.assert_equal(base.sagittal_vector, -base.col_vector)
\ No newline at end of file
+    np.testing.assert_equal(base.sagittal_vector, -base.col_vector)
+
+
+def test_centre():
+    transform = _transforms.transformation_matrix([1, 0, 0],
+                                                  [0, 1, 0],
+                                                  [-127.5, -99.5, -99],
+                                                  [1, 1, 2])
+    base = suspect.base.ImageBase(np.zeros((128, 256, 256)), transform=transform)
+    np.testing.assert_equal(base.centre, [0, 28, 28])
+
+
+def test_resample():
+    transform = _transforms.transformation_matrix([1, 0, 0],
+                                                  [0, 1, 0],
+                                                  [-1.5, -0.5, -1.5],
+                                                  [1, 1, 1])
+    data = np.random.rand(4, 2, 4)
+    base = suspect.base.ImageBase(data, transform=transform)
+    resampled_1 = base.resample([1, 0, 0], [0, 1, 0], [1, 2, 4],
+                                centre=[0, 0, 0.5])
+    np.testing.assert_equal(base[2], resampled_1)
+    resampled_1_transform = np.array([[1, 0, 0, -1.5],
+                                      [0, 1, 0, -0.5],
+                                      [0, 0, 1, 0.5],
+                                      [0, 0, 0, 1]])
+    np.testing.assert_equal(resampled_1_transform, resampled_1.transform)
+    resampled_2 = base.resample([0, 1, 0], [0, 0, 2], [1, 4, 2],
+                                centre=[0.5, 0, 0])
+    np.testing.assert_equal(base[:, :, 2], resampled_2)
+    resampled_3 = base.resample([1, 0, 0], [0, 0, 1], [1, 4, 4],
+                                centre=[0, 0.5, 0])
+    np.testing.assert_equal(base[:, 1], resampled_3)
diff --git a/tests/test_mrs/test_fitting.py b/tests/test_mrs/test_fitting.py
index 76287ad..609feb5 100644
--- a/tests/test_mrs/test_fitting.py
+++ b/tests/test_mrs/test_fitting.py
@@ -47,7 +47,7 @@ def test_gaussian(fixed_fid):
 
     numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["fwhm"], 50.0, rtol=1e-2)
     numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["amplitude"], 1.0, rtol=2e-2)
-    numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["frequency"], 0.0, atol=5e-2)
+    numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["frequency"], 0.0, atol=1e-1)
 
     numpy.testing.assert_allclose(fitting_result["fit"], fixed_fid, atol=0.001)
 
@@ -123,7 +123,7 @@ def test_missing_peak_phase(fixed_fid):
 
     numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["fwhm"], 50.0, rtol=5e-2)
     numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["amplitude"], 1.0, rtol=5e-2)
-    numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["frequency"], 0.0, atol=5e-2)
+    numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["frequency"], 0.0, atol=1e-1)
 
     numpy.testing.assert_allclose(fitting_result["fit"], fixed_fid, atol=0.001)
 
@@ -241,7 +241,7 @@ def test_dependencies(fixed_fid_sum):
 
     numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["fwhm"], 50.0, rtol=1e-2)
     numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["amplitude"], 1.0, rtol=2e-2)
-    numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["frequency"], 0.0, atol=5e-2)
+    numpy.testing.assert_allclose(fitting_result["model"]["pcr"]["frequency"], 0.0, atol=1e-1)
 
     numpy.testing.assert_allclose(fitting_result["fit"], fixed_fid_sum, atol=0.001)
 
diff --git a/tests/test_mrs/test_mrsdata.py b/tests/test_mrs/test_mrsdata.py
index cb34c84..e98f1b0 100644
--- a/tests/test_mrs/test_mrsdata.py
+++ b/tests/test_mrs/test_mrsdata.py
@@ -2,6 +2,7 @@ import suspect
 
 import numpy
 import pytest
+from suspect import _transforms
 
 
 def test_create_mrs():
@@ -40,3 +41,14 @@ def test_spectrum_2_fid():
     numpy.testing.assert_equal(spectrum, numpy.fft.fftshift(numpy.fft.fft(data)))
     fid = spectrum.fid()
     numpy.testing.assert_equal(data, fid)
+
+
+def test_centre():
+    position = (12.4, -9.8, 11.0)
+    transform = _transforms.transformation_matrix((1, 0, 0),
+                                                  (0, 1, 0),
+                                                  (12.4, -9.8, 11.0),
+                                                  (10, 10, 10))
+    data = suspect.MRSData(numpy.ones(1024, 'complex'), 5e-4, 123, transform=transform)
+    numpy.testing.assert_equal(data.centre, position)
+    numpy.testing.assert_equal(data.position, position)
diff --git a/tests/test_mrs/test_processing/test_phasing.py b/tests/test_mrs/test_processing/test_phasing.py
index 7889c45..6d23666 100644
--- a/tests/test_mrs/test_processing/test_phasing.py
+++ b/tests/test_mrs/test_processing/test_phasing.py
@@ -40,8 +40,8 @@ def test_acme_first():
 
     out_0, out_1 = suspect.processing.phase.acme(sample_data)
 
-    np.testing.assert_allclose(in_0, -out_0, rtol=0.15)
-    np.testing.assert_allclose(in_1, -out_1, rtol=0.15)
+    np.testing.assert_allclose(in_0, -out_0, rtol=0.25)
+    np.testing.assert_allclose(in_1, -out_1, rtol=0.25)
 
 
 def test_acme_range_hz():