New methods on Atlas image classes for mask-based queries. Not tested. Also

adjusted resample method to force float output (although this might change).

fsl/data/atlases.py

@@ -623,23 +623,60 @@ class LabelAtlas(Atlas):
         Atlas.__init__(self, atlasDesc, resolution, True)
 
 
-    def label(self, worldLoc):
-        """Looks up and returns the label of the region at the given world
-        location, or ``None`` if the location is out of bounds.
+    def label(self, location, *args, **kwargs):
+        """Looks up and returns the label of the region at the given
+        location.
+
+        :arg location: Can be one of the following:
+
+                        - A sequence of three values, interpreted as
+                          atlas coordinates. In this case, :meth:`coordLabel`
+                          is called.
+
+                        - An :class:`.Image` which is interpreted as a
+                          weighted mask. In this case, :meth:`maskLabel` is
+                          called.
+
+        All other arguments are passed through to the :meth:`coordLabel` or
+        :meth:`maskLabel` methods.
+
+
+        :returns: The return value of either :meth:`coordLabel` or
+                  :meth:`maskLabel`.
+        """
+
+        if isinstance(location, fslimage.Image):
+            return self.maskLabel(location, *args, **kwargs)
+        else:
+            return self.coordLabel(location, *args, **kwargs)
+
+
+    def coordLabel(self, loc, voxel=False):
+        """Looks up and returns the label at the given location.
+
+        :arg loc:   A sequence of three values, interpreted as atlas
+                    coordinates. In this case, :meth:`coordLabel` is called.
+
+        :arg voxel: Defaults to ``False``. If ``True``, the ``location``
+                    is interpreted as voxel coordinates.
+
+        :returns:   The label at the given coordinates, or ``None`` if the
+                    coordinates are out of bounds.
         """
 
-        voxelLoc = transform.transform([worldLoc], self.worldToVoxMat)[0]
-        voxelLoc = [int(v) for v in voxelLoc.round()]
+        if not voxel:
+            loc = transform.transform([loc], self.worldToVoxMat)[0]
+            loc = [int(v) for v in loc.round()]
 
-        if voxelLoc[0] <  0             or \
-           voxelLoc[1] <  0             or \
-           voxelLoc[2] <  0             or \
-           voxelLoc[0] >= self.shape[0] or \
-           voxelLoc[1] >= self.shape[1] or \
-           voxelLoc[2] >= self.shape[2]:
+        if loc[0] <  0             or \
+           loc[1] <  0             or \
+           loc[2] <  0             or \
+           loc[0] >= self.shape[0] or \
+           loc[1] >= self.shape[1] or \
+           loc[2] >= self.shape[2]:
             return None
 
-        val = self[voxelLoc[0], voxelLoc[1], voxelLoc[2]]
+        val = self[loc[0], loc[1], loc[2]]
 
         if self.desc.atlasType == 'label':
             return val
@@ -648,6 +685,50 @@ class LabelAtlas(Atlas):
             return val - 1
 
 
+    def maskLabel(self, mask):
+        """Looks up and returns the proportions of all regions that are present
+        in the given ``mask``.
+
+        :arg mask: A 3D :class:`.Image`` which is interpreted as a weighted
+                   mask. If the ``mask`` shape does not match that of this
+                   ``LabelAtlas``, it is resampled using
+                   :meth:`.Image.resample`, with linear interpolation.
+
+        :returns:  A tuple containing:
+
+                     - A sequence of all labels which are present in the mask
+                     - A sequence containing the proportion, within the mask,
+                       of each present label.
+        """
+
+        # Make sure that the mask has the
+        # same number of voxels as the
+        # atlas image
+        mask     = mask.resample(self.shape[:3], order=1)
+        boolmask = mask > 0
+
+        # Extract the labels that are in
+        # the mask, and their corresponding
+        # mask weights
+        vals    = self[boolmask]
+        weights = vals * mask[boolmask]
+        labels  = np.unique(vals)
+        props   = []
+
+        for label in labels:
+
+            # Figure out the number of all voxels
+            # in the mask with this label, weighted
+            # by the mask
+            prop = ((vals == label) * weights).sum()
+
+            # Normalise it to be a proportion
+            # of all voxels in the mask
+            props.append(prop / float(len(vals)))
+
+        return labels, props
+
+
 class ProbabilisticAtlas(Atlas):
     """A 4D atlas which contains one volume for each region.
 
@@ -666,28 +747,98 @@ class ProbabilisticAtlas(Atlas):
         Atlas.__init__(self, atlasDesc, resolution, False)
 
 
-    def proportions(self, worldLoc):
+    def proportions(self, location, *args, **kwargs):
+        """Looks up and returns the proportions of of all regions at the given
+        location.
+
+        :arg location: Can be one of the following:
+
+                        - A sequence of three values, interpreted as atlas
+                          coordinates. In this case, :meth:`coordProportions`
+                          is called.
+
+                        - An :class:`.Image` which is interpreted as a
+                          weighted mask. In this case, :meth:`maskProportions`
+                          is called.
+
+        All other arguments are passed through to the :meth:`coordProportions`
+        or :meth:`maskProportions` methods.
+
+
+        :returns: The return value of either :meth:`coordProportions` or
+                  :meth:`maskProportions`.
+        """
+
+        if isinstance(location, fslimage.Image):
+            return self.maskProportions(location, *args, **kwargs)
+        else:
+            return self.coordProportions(location, *args, **kwargs)
+
+
+    def coordProportions(self, loc, voxel=False):
         """Looks up the region probabilities for the given location.
 
-        :arg worldLoc: Location in the world coordinate system.
+        :arg loc:   A sequence of three values, interpreted as atlas
+                    world or voxel coordinates.
+
+        :arg voxel: Defaults to ``False``. If ``True``, the ``loc``
+                    argument is interpreted as voxel coordinates.
 
         :returns: a list of values, one per region, which represent
                   the probability of each region for the specified
                   location. Returns an empty list if the given
                   location is out of bounds.
         """
-        voxelLoc = transform.transform([worldLoc], self.worldToVoxMat)[0]
-        voxelLoc = [int(v) for v in voxelLoc.round()]
-
-        if voxelLoc[0] <  0             or \
-           voxelLoc[1] <  0             or \
-           voxelLoc[2] <  0             or \
-           voxelLoc[0] >= self.shape[0] or \
-           voxelLoc[1] >= self.shape[1] or \
-           voxelLoc[2] >= self.shape[2]:
+
+        if not voxel:
+            loc = transform.transform([loc], self.worldToVoxMat)[0]
+            loc = [int(v) for v in loc.round()]
+
+        if loc[0] <  0             or \
+           loc[1] <  0             or \
+           loc[2] <  0             or \
+           loc[0] >= self.shape[0] or \
+           loc[1] >= self.shape[1] or \
+           loc[2] >= self.shape[2]:
             return []
 
-        return self[voxelLoc[0], voxelLoc[1], voxelLoc[2], :]
+        return self[loc[0], loc[1], loc[2], :]
+
+
+    def maskProportions(self, mask):
+        """Looks up the probabilities of all regions in the given ``mask``.
+
+        :arg mask: A 3D :class:`.Image`` which is interpreted as a weighted
+                   mask. If the ``mask`` shape does not match that of this
+                   ``ProbabilisticAtlas``, it is resampled using
+                   :meth:`.Image.resample`, with linear interpolation.
+
+        :returns:  A tuple containing:
+
+                     - A sequence of all labels which are present in the mask
+                     - A sequence containing the proportion, within the mask,
+                       of each present label.
+        """
+
+        labels = []
+        props  = []
+
+        # Make sure that the mask has the same
+        # number of voxels as the atlas image
+        mask     = mask.resample(self.shape[:3], order=1)
+        boolmask = mask > 0
+
+        for label in range(self.shape[3]):
+
+            vals = self[..., label]
+            vals = vals[boolmask] * mask[boolmask]
+            prop = vals.mean()
+
+            if prop != 0:
+                labels.append(label)
+                props .append(prop)
+
+        return labels, props
 
 
 registry            = AtlasRegistry()

fsl/data/image.py

@@ -1128,6 +1128,9 @@ class Image(Nifti):
 
         All other arguments are passed through to the ``scipy.ndimage.zoom``
         function.
+
+        :returns: A ``numpy`` array of shape ``shape``, containing an
+                  interpolated copy of the data in this ``Image``.
         """
 
         if sliceobj is None:
@@ -1137,6 +1140,7 @@ class Image(Nifti):
         data  = self[sliceobj]
 
         if tuple(data.shape) != tuple(shape):
+            data  = np.array(data, dtype=np.float, copy=False)
             zooms = [float(shape[i]) / data.shape[i] for i in range(ndims)]
             data  = ndimage.zoom(data, zooms, **kwargs)