Repository: openmrslab/suspect

Issue:
ENH: Mask function for voxel on image
It should be possible to call a `create_mask()` function which produces a binary mask showing which voxels of a reference image are inside a spectroscopy voxel. In the default case this would be for all the voxels in the spectroscopy, but it should also be possible to specify a specific voxel or list of voxels to be used instead, so that masks for individual CSI voxels are also possible.

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

Coverage increased (+0.2%) to 78.728% when pulling **2e5c6eaa68957507ad19205cfbc872dbb3011ace on 62_voxel_mask** into **a70490c09196d7996ecfd99401eb68765ae25a67 on master**.


Patch:
diff --git a/suspect/image/__init__.py b/suspect/image/__init__.py
index 14e78f7..cffdc25 100644
--- a/suspect/image/__init__.py
+++ b/suspect/image/__init__.py
@@ -1,1 +1,2 @@
-from suspect.image._image import *
\ No newline at end of file
+from suspect.image._image import *
+from ._mask import create_mask
\ No newline at end of file
diff --git a/suspect/image/_mask.py b/suspect/image/_mask.py
new file mode 100644
index 0000000..e31cff1
--- /dev/null
+++ b/suspect/image/_mask.py
@@ -0,0 +1,46 @@
+import numpy as np
+
+
+def create_mask(source_image, ref_image, voxels=None):
+    """
+    Creates a volumetric mask for the source_image voxel in the coordinate
+    system of the ref_image volume.
+    
+    Parameters
+    ----------
+    source_image : MRSBase
+        The spectroscopy volume from which to create the mask.
+    ref_image : ImageBase
+        The reference image volume which defines the coordinate system for
+        the mask.
+    
+    Returns
+    -------
+    numpy.ndarray
+        Boolean array with the same shape as ref_image, True for all voxels
+        inside source_image, false for all others.
+    """
+
+    # create a grid of coordinates for all points in the ref_image
+    # the ref_image has coord index order [z, y, x] so we reverse the shape
+    # to get the indices in (x, y, z) format for the coordinate conversion
+    ref_coords = np.mgrid[[range(0, size) for size in ref_image.shape[::-1]]]
+
+    # mgrid puts the (x, y, z) tuple at the front, we want it at the back
+    ref_coords = np.moveaxis(ref_coords, 0, -1)
+
+    # now we can apply to_scanner and from_scanner to convert from ref coords
+    # into source coords
+    scanner_coords = ref_image.to_scanner(ref_coords)
+    source_coords = source_image.from_scanner(scanner_coords)
+
+    # now check whether the source_coords are in the selected voxel
+    # TODO for now, we assume single voxel data until issue 50 is resolved
+
+    # have to transpose the result to get it to match the shape of ref_image
+    return np.all((source_coords[..., 0] < 0.5,
+                   source_coords[..., 0] >= -0.5,
+                   source_coords[..., 1] >= -0.5,
+                   source_coords[..., 2] >= -0.5,
+                   source_coords[..., 1] < 0.5,
+                   source_coords[..., 2] < 0.5), axis=0).T


Test patch:
diff --git a/tests/test_mrs/test_image.py b/tests/test_mrs/test_image.py
new file mode 100644
index 0000000..70fca64
--- /dev/null
+++ b/tests/test_mrs/test_image.py
@@ -0,0 +1,23 @@
+import suspect
+
+import suspect._transforms
+
+import numpy as np
+
+
+def test_simple_mask():
+    source_transform = suspect._transforms.transformation_matrix([1, 0, 0],
+                                                                 [0, 1, 0],
+                                                                 [5, 0, 0],
+                                                                 [10, 10, 10])
+    ref_transform = suspect._transforms.transformation_matrix([1, 0, 0],
+                                                              [0, 1, 0],
+                                                              [-10, -5, -5],
+                                                              [1, 1, 1])
+    source_volume = suspect.MRSBase(np.ones(1024), 1e-3, 123, transform=source_transform)
+    ref_volume = suspect.base.ImageBase(np.zeros((20, 20, 20)), transform=ref_transform)
+    mask = suspect.image.create_mask(source_volume, ref_volume)
+    assert ref_volume.shape == mask.shape
+    mask_target = np.zeros_like(ref_volume)
+    mask_target[0:10, 0:10, 10:20] = 1
+    np.testing.assert_equal(mask_target.astype('bool'), mask)