From 2516b4f3fe80a2e77895f6bc508f7bcb781e4e24 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Fri, 18 Aug 2017 15:31:33 -0400 Subject: [PATCH 01/36] ENH: Add funcs.crop_image --- nibabel/funcs.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/nibabel/funcs.py b/nibabel/funcs.py index 240b20f802..ed5dd282c5 100644 --- a/nibabel/funcs.py +++ b/nibabel/funcs.py @@ -219,3 +219,39 @@ def _aff_is_diag(aff): ''' Utility function returning True if affine is nearly diagonal ''' rzs_aff = aff[:3, :3] return np.allclose(rzs_aff, np.diag(np.diag(rzs_aff))) + + +def crop_image(img, mask=None): + ''' Crop ``img`` to smallest region that contains all non-zero data + + The image is cropped in the current orientation; no rotations or resampling + are performed. + The affine matrix is updated with the new intercept, so that all values are + found at the same RAS locations. + + Parameters + ---------- + img : ``spatialimage`` + mask : ``spatialimage``, optional + If supplied, use the bounding box of ``mask``, rather than ``img`` + + Returns + ------- + cropped_img : ``spatialimage`` + Version of `img` with cropped data array and updated affine matrix + ''' + if mask is None: + mask = img + elif not np.allclose(img.affine == mask.affine): + raise ValueError('Affine for image does not match affine for mask') + + bounds = np.sort(np.vstack(np.nonzero(mask.get_data())))[:, [0, -1]] + x, y, z = bounds + new_origin = np.vstack((bounds[:, [0]], [1])) + + new_data = img.get_data()[x[0]:x[1] + 1, y[0]:y[1] + 1, z[0]:z[1] + 1] + + new_aff = img.affine.copy() + new_aff[:, [3]] = img.affine.dot(new_origin) + + return img.__class__(new_data, new_aff, img.header) From f031494034687e38bca8df16ab223e0f436d8c06 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 6 Sep 2017 10:38:37 -0400 Subject: [PATCH 02/36] ENH: crop_image to bounds, adding margin --- nibabel/funcs.py | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/nibabel/funcs.py b/nibabel/funcs.py index ed5dd282c5..f953afe56a 100644 --- a/nibabel/funcs.py +++ b/nibabel/funcs.py @@ -221,37 +221,42 @@ def _aff_is_diag(aff): return np.allclose(rzs_aff, np.diag(np.diag(rzs_aff))) -def crop_image(img, mask=None): - ''' Crop ``img`` to smallest region that contains all non-zero data +def crop_image(img, bounds, margin=0): + ''' Crop ``img`` to specified bounds - The image is cropped in the current orientation; no rotations or resampling - are performed. + The image is cropped in the current orientation; no rotation or resampling + is performed. The affine matrix is updated with the new intercept, so that all values are found at the same RAS locations. Parameters ---------- img : ``spatialimage`` - mask : ``spatialimage``, optional - If supplied, use the bounding box of ``mask``, rather than ``img`` + Image to be cropped + bounds : (3, 2) array-like + Minimum and maximum (inclusive) voxel indices, in each dimension, to be + included in cropped image + margin : int, optional + Margin, in voxels, to include on each side of cropped image Returns ------- cropped_img : ``spatialimage`` Version of `img` with cropped data array and updated affine matrix ''' - if mask is None: - mask = img - elif not np.allclose(img.affine == mask.affine): - raise ValueError('Affine for image does not match affine for mask') - bounds = np.sort(np.vstack(np.nonzero(mask.get_data())))[:, [0, -1]] + try: + bounds = np.asanyarray(bounds) + np.array([-margin, margin]) + assert bounds.shape == (3, 2) + except (ValueError, AssertionError) as err: + raise ValueError("bounds must be interpretable as a 3x2 array") + x, y, z = bounds - new_origin = np.vstack((bounds[:, [0]], [1])) + new_origin = np.vstack((bounds[:, [0]], 1)) - new_data = img.get_data()[x[0]:x[1] + 1, y[0]:y[1] + 1, z[0]:z[1] + 1] + bounded_data = img.get_data()[x[0]:x[1] + 1, y[0]:y[1] + 1, z[0]:z[1] + 1] new_aff = img.affine.copy() new_aff[:, [3]] = img.affine.dot(new_origin) - return img.__class__(new_data, new_aff, img.header) + return img.__class__(bounded_data, new_aff, img.header) From 008255ad8974a7ff0ef655a9cc5f7609c8ef044d Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 7 Sep 2017 11:17:45 -0400 Subject: [PATCH 03/36] STY: Do not capture exception in variable --- nibabel/funcs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/funcs.py b/nibabel/funcs.py index f953afe56a..db242e8de6 100644 --- a/nibabel/funcs.py +++ b/nibabel/funcs.py @@ -248,7 +248,7 @@ def crop_image(img, bounds, margin=0): try: bounds = np.asanyarray(bounds) + np.array([-margin, margin]) assert bounds.shape == (3, 2) - except (ValueError, AssertionError) as err: + except (ValueError, AssertionError): raise ValueError("bounds must be interpretable as a 3x2 array") x, y, z = bounds From a6af773d204a8438f84f409a4f0a0f5d66c7ad87 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 21 Sep 2017 11:02:35 -0400 Subject: [PATCH 04/36] ENH: Permit cropping with negative indices, give more specific errors --- nibabel/funcs.py | 35 ++++++++++++++++++++++++++++++----- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/nibabel/funcs.py b/nibabel/funcs.py index db242e8de6..2258d9fe8f 100644 --- a/nibabel/funcs.py +++ b/nibabel/funcs.py @@ -245,16 +245,41 @@ def crop_image(img, bounds, margin=0): Version of `img` with cropped data array and updated affine matrix ''' - try: - bounds = np.asanyarray(bounds) + np.array([-margin, margin]) - assert bounds.shape == (3, 2) - except (ValueError, AssertionError): + shape = np.reshape(img.shape[:3], (3, 1)) + bounds = np.asanyarray(bounds) + if bounds.shape != (3, 2): raise ValueError("bounds must be interpretable as a 3x2 array") + elif np.any(bounds > shape): + raise ValueError("bounds must not exceed image dimensions") + + # Permit negative bounds + if np.any(bounds < 0): + bounds = (bounds + shape) % shape + + if np.any(bounds < 0): + raise ValueError("negative bounds may not exceed image dimensions") + elif np.any(bounds[:, 0] > bounds[:, 1]): + raise ValueError("degenerate (0 width) crops are not permitted") + + # Add margin in all directions + bounds += np.array([-margin, margin]) + + # Set min/max + bounds[bounds < 0] = 0 + over = bounds[:, 1] > shape.reshape(-1) - 1 + bounds[over, 1] = shape[over, 0] - 1 + + # Include upper bounds + bounds[:, 1] += 1 + + # Return original image if no cropping required + if np.array_equal(bounds, np.hstack(([[0], [0], [0]], shape))): + return img x, y, z = bounds new_origin = np.vstack((bounds[:, [0]], 1)) - bounded_data = img.get_data()[x[0]:x[1] + 1, y[0]:y[1] + 1, z[0]:z[1] + 1] + bounded_data = img.get_data()[x[0]:x[1], y[0]:y[1], z[0]:z[1]] new_aff = img.affine.copy() new_aff[:, [3]] = img.affine.dot(new_origin) From 5c1117430065190cbaac1a2a9e3f725f99ce9734 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 26 Sep 2017 20:09:53 -0400 Subject: [PATCH 05/36] TEST: Initial crop_image test --- nibabel/tests/test_funcs.py | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/nibabel/tests/test_funcs.py b/nibabel/tests/test_funcs.py index 6032c08672..91fd83aff0 100644 --- a/nibabel/tests/test_funcs.py +++ b/nibabel/tests/test_funcs.py @@ -11,7 +11,7 @@ import numpy as np -from ..funcs import concat_images, as_closest_canonical, OrientationError +from ..funcs import concat_images, as_closest_canonical, OrientationError, crop_image from ..analyze import AnalyzeImage from ..nifti1 import Nifti1Image from ..loadsave import save @@ -193,3 +193,22 @@ def test_closest_canonical(): img.header.set_dim_info(None, None, 2) xyz_img = as_closest_canonical(img) assert_true(xyz_img.header.get_dim_info() == (None, None, 1)) + + +def test_crop_image(): + # Use 32-bit data so that the AnalyzeImage class doesn't complain + arr = np.arange(60).reshape((5, 3, 4, 1)).astype(np.int32) + + img = AnalyzeImage(arr, np.eye(4)) + + cropped_img = crop_image(img, [[1, 3], [1, 1], [1, 2]]) + assert_equal(cropped_img.shape, (3, 1, 2, 1)) + assert_array_equal(cropped_img.affine, [[1, 0, 0, 1], + [0, 1, 0, 1], + [0, 0, 1, 1], + [0, 0, 0, 1]]) + + cropped_img = crop_image(img, [[1, 3], [1, 1], [1, 2]], margin=1) + assert_equal(cropped_img.shape, (5, 3, 4, 1)) + assert_array_equal(cropped_img.affine, img.affine) + assert_true(cropped_img is img) From 9b7a35cdbd07db2e501c5a6f379ce6b17480bc53 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 27 Sep 2017 12:18:02 -0400 Subject: [PATCH 06/36] ENH: Enable SpatialImage.{slice,__getitem__} --- nibabel/spatialimages.py | 116 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 111 insertions(+), 5 deletions(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index b88b3e8538..285df28f5c 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -461,12 +461,118 @@ def from_image(klass, img): klass.header_class.from_header(img.header), extra=img.extra.copy()) + @staticmethod + def requires_downsampling(idx): + for slicer in idx: + if isinstance(slicer, slice) and slicer.step not in (1, None): + return True + elif np.asanyarray(slicer).shape != (1,): + return True + + return False + + def slice_afffine(self, idx): + """ Retrieve affine for current image, if sliced by a given index + + Applies scaling if down-sampling is applied, and adjusts the intercept + to account for any cropping. + + Parameters + ---------- + idx : numpy-compatible slice index + + Returns + ------- + affine : (4,4) ndarray + Affine with updated scale and intercept + """ + # Expand ellipsis to retrieve bounds + if Ellipsis in idx: + i = idx.index(Ellipsis) + subslices = tuple(slice(None) + for _ in range(len(self.shape) - len(idx.shape))) + idx = idx[:i] + subslices + idx[i + 1:] + + origin = np.array([[0], [0], [0], [1]]) + scales = np.eye(3, dtype=int) + + for i, slicer in enumerate(idx[:3]): + new_origin[i] = (slicer.start or 0 if isinstance(slicer, slice) + else np.asanyarray(slicer)[0]) + if slicer.step is not None: + scales[i, i] = slicer.step + + affine = self.affine.copy() + affine[:3, :3] = scales.dot(self.affine[:3, :3]) + affine[:, [3]] = self.affine.dot(origin) + return affine + + def slice(self, idx, allow_downsampling=False): + """ Slice image to specification + + Nibabel performs no spatial filtering, so, in order to avoid aliasing + data, only slices that specify cropping data are permitted. + If the image has been filtered, or aliasing is not a concern, + ``allow_downsampling=True`` will disable this check and permit more + complicated slices. + + The image is resliced in the current orientation; no rotation or + resampling is performed. + The affine matrix is updated with the new intercept (and scales, if + down-sampling is used), so that all values are found at the same RAS + locations. + + Slicing may include non-spatial dimensions. + However, this method does not currently adjust the repetition time in + the image header. + + Parameters + ---------- + idx : numpy-compatible slice index + allow_downsampling : bool (default: False) + Permit indices that specify down-sampling + + Returns + ------- + resliced_img : ``spatialimage`` + Version of `img` with resliced data array and updated affine matrix + """ + if not isinstance(idx, tuple): + idx = (idx,) + + if not allow_downsampling and self.requires_downsampling(idx): + raise IndexError('slicing requires downsampling, please use ' + 'img.slice(..., allow_downsampling=True) instead') + + # Get bounded data, allow numpy to complain about any indexing errors + dataobj = self.dataobj[slices] + affine = self.slice_affine(idx) + return self.__class__(dataobj, affine, self.header) + def __getitem__(self, idx): - ''' No slicing or dictionary interface for images - ''' - raise TypeError("Cannot slice image objects; consider slicing image " - "array data with `img.dataobj[slice]` or " - "`img.get_data()[slice]`") + """ Crop image to specified slices + + The image is cropped in the current orientation; no rotation or + resampling is performed. + The affine matrix is updated with the new intercept, so that all + values are found at the same RAS locations. + + Cropping may include non-spatial dimensions. + + Down-sampling and other slicing operations may be performed with + ``img.slice(..., allow_downsampling=True)``. + + Parameters + ---------- + idx : tuple containing (slice(step=None), Ellipsis, array-like with size (1,)) + numpy-compatible slice index describing cropping only + + Returns + ------- + cropped_img : ``spatialimage`` + Version of `img` with cropped data array and updated affine matrix + """ + return self.slice(idx) def orthoview(self): """Plot the image using OrthoSlicer3D From 0c376375ca7fcdff3ae65737c9f1326afafdf324 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 28 Sep 2017 15:59:49 -0400 Subject: [PATCH 07/36] ENH: Enable minimally-fancy indexing --- nibabel/fileslice.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/nibabel/fileslice.py b/nibabel/fileslice.py index e55f48c127..40d7a130b6 100644 --- a/nibabel/fileslice.py +++ b/nibabel/fileslice.py @@ -54,6 +54,10 @@ def is_fancy(sliceobj): # slice or Ellipsis or None OK for basic if isinstance(slicer, slice) or slicer in (None, Ellipsis): continue + # Allow single items in a dimension to be sliced + ary = np.asanyarray(slicer) + if ary.shape == (1,): + slicer = ary[0] try: int(slicer) except TypeError: @@ -113,9 +117,16 @@ def canonical_slicers(sliceobj, shape, check_inds=True): can_slicers.extend((slice(None),) * n_ellided) n_real += n_ellided continue - # int / slice indexing cases + # int / [int] / slice indexing cases dim_len = shape[n_real] n_real += 1 + + # Treat [x] as "x:x+1" + ary = np.asanyarray(slicer) + if ary.shape == (1,): + can_slicers.append(slice(int(ary[0]), int(ary[0]) + 1)) + continue + try: # test for integer indexing slicer = int(slicer) except TypeError: # should be slice object From 00f6ed1c52484bdc3e4a53c116607c844a14db11 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 28 Sep 2017 16:00:14 -0400 Subject: [PATCH 08/36] RF: Move to img.slicer approach --- nibabel/spatialimages.py | 114 ++++++++++++++------------------------- 1 file changed, 41 insertions(+), 73 deletions(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 285df28f5c..950fe2794a 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -140,6 +140,7 @@ from .filebasedimages import ImageFileError # flake8: noqa; for back-compat from .viewers import OrthoSlicer3D from .volumeutils import shape_zoom_affine +from .fileslice import canonical_slicers from .deprecated import deprecate_with_version from .orientations import apply_orientation, inv_ornt_aff @@ -325,6 +326,23 @@ class SpatialImage(DataobjImage): ''' Template class for volumetric (3D/4D) images ''' header_class = SpatialHeader + class Slicer(object): + def __init__(self, img): + self.img = img + + def __getitem__(self, idx): + idx = self._validate(idx) + dataobj = self.img.dataobj[idx] + affine = self.img.slice_affine(idx) + return self.img.__class__(dataobj, affine, self.img.header) + + def _validate(self, idx): + idx = canonical_slicers(idx, self.img.shape) + for slicer in idx: + if isinstance(slicer, int): + raise IndexError("Cannot slice image with integer") + return idx + def __init__(self, dataobj, affine, header=None, extra=None, file_map=None): ''' Initialize image @@ -461,17 +479,7 @@ def from_image(klass, img): klass.header_class.from_header(img.header), extra=img.extra.copy()) - @staticmethod - def requires_downsampling(idx): - for slicer in idx: - if isinstance(slicer, slice) and slicer.step not in (1, None): - return True - elif np.asanyarray(slicer).shape != (1,): - return True - - return False - - def slice_afffine(self, idx): + def slice_affine(self, idx): """ Retrieve affine for current image, if sliced by a given index Applies scaling if down-sampling is applied, and adjusts the intercept @@ -486,38 +494,33 @@ def slice_afffine(self, idx): affine : (4,4) ndarray Affine with updated scale and intercept """ - # Expand ellipsis to retrieve bounds - if Ellipsis in idx: - i = idx.index(Ellipsis) - subslices = tuple(slice(None) - for _ in range(len(self.shape) - len(idx.shape))) - idx = idx[:i] + subslices + idx[i + 1:] + idx = canonical_slicers(idx, self.shape, check_inds=False)[:3] origin = np.array([[0], [0], [0], [1]]) scales = np.eye(3, dtype=int) for i, slicer in enumerate(idx[:3]): - new_origin[i] = (slicer.start or 0 if isinstance(slicer, slice) - else np.asanyarray(slicer)[0]) - if slicer.step is not None: - scales[i, i] = slicer.step + if isinstance(slicer, slice): + origin[i] = slicer.start or 0 + if slicer.step is not None: + scales[i, i] = slicer.step + elif isinstance(slicer, int): + origin[i] = slicer + # If slicer is None, nothing to do affine = self.affine.copy() affine[:3, :3] = scales.dot(self.affine[:3, :3]) affine[:, [3]] = self.affine.dot(origin) return affine - def slice(self, idx, allow_downsampling=False): - """ Slice image to specification - - Nibabel performs no spatial filtering, so, in order to avoid aliasing - data, only slices that specify cropping data are permitted. - If the image has been filtered, or aliasing is not a concern, - ``allow_downsampling=True`` will disable this check and permit more - complicated slices. + @property + def slicer(self): + """ Slicer object that returns cropped and subsampled images The image is resliced in the current orientation; no rotation or - resampling is performed. + resampling is performed, and no attempt is made to filter the image + to avoid aliasing. + The affine matrix is updated with the new intercept (and scales, if down-sampling is used), so that all values are found at the same RAS locations. @@ -525,54 +528,19 @@ def slice(self, idx, allow_downsampling=False): Slicing may include non-spatial dimensions. However, this method does not currently adjust the repetition time in the image header. - - Parameters - ---------- - idx : numpy-compatible slice index - allow_downsampling : bool (default: False) - Permit indices that specify down-sampling - - Returns - ------- - resliced_img : ``spatialimage`` - Version of `img` with resliced data array and updated affine matrix """ - if not isinstance(idx, tuple): - idx = (idx,) + return self.Slicer(self) - if not allow_downsampling and self.requires_downsampling(idx): - raise IndexError('slicing requires downsampling, please use ' - 'img.slice(..., allow_downsampling=True) instead') - - # Get bounded data, allow numpy to complain about any indexing errors - dataobj = self.dataobj[slices] - affine = self.slice_affine(idx) - return self.__class__(dataobj, affine, self.header) def __getitem__(self, idx): - """ Crop image to specified slices - - The image is cropped in the current orientation; no rotation or - resampling is performed. - The affine matrix is updated with the new intercept, so that all - values are found at the same RAS locations. + ''' No slicing or dictionary interface for images - Cropping may include non-spatial dimensions. - - Down-sampling and other slicing operations may be performed with - ``img.slice(..., allow_downsampling=True)``. - - Parameters - ---------- - idx : tuple containing (slice(step=None), Ellipsis, array-like with size (1,)) - numpy-compatible slice index describing cropping only - - Returns - ------- - cropped_img : ``spatialimage`` - Version of `img` with cropped data array and updated affine matrix - """ - return self.slice(idx) + Use the slicer attribute to perform cropping and subsampling at your + own risk. + ''' + raise TypeError("Cannot slice image objects; consider slicing image " + "array data with `img.dataobj[slice]` or " + "`img.get_data()[slice]`") def orthoview(self): """Plot the image using OrthoSlicer3D From ee0c2b4e512bd7ab64c2acac4f9d4fd3e3b59019 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 28 Sep 2017 16:34:40 -0400 Subject: [PATCH 09/36] RF: Simplify affine transform --- nibabel/spatialimages.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 950fe2794a..f327ac3c60 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -496,22 +496,24 @@ def slice_affine(self, idx): """ idx = canonical_slicers(idx, self.shape, check_inds=False)[:3] - origin = np.array([[0], [0], [0], [1]]) - scales = np.eye(3, dtype=int) - - for i, slicer in enumerate(idx[:3]): + # Transform: + # sx 0 0 tx + # 0 sy 0 ty + # 0 0 sz tz + # 0 0 0 1 + transform = np.eye(4, dtype=int) + + for i, slicer in enumerate(idx): if isinstance(slicer, slice): - origin[i] = slicer.start or 0 - if slicer.step is not None: - scales[i, i] = slicer.step + if slicer.step == 0: + raise ValueError("slice step cannot be 0") + transform[i, i] = slicer.step if slicer.step is not None else 1 + transform[i, 3] = slicer.start or 0 elif isinstance(slicer, int): - origin[i] = slicer + transform[i, 3] = slicer # If slicer is None, nothing to do - affine = self.affine.copy() - affine[:3, :3] = scales.dot(self.affine[:3, :3]) - affine[:, [3]] = self.affine.dot(origin) - return affine + return self.affine.dot(transform) @property def slicer(self): From 8c113d4b4779b5f556ea79618a17cac655652319 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 3 Oct 2017 08:42:18 -0400 Subject: [PATCH 10/36] REVERT fileslice modifications --- nibabel/fileslice.py | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/nibabel/fileslice.py b/nibabel/fileslice.py index 40d7a130b6..e55f48c127 100644 --- a/nibabel/fileslice.py +++ b/nibabel/fileslice.py @@ -54,10 +54,6 @@ def is_fancy(sliceobj): # slice or Ellipsis or None OK for basic if isinstance(slicer, slice) or slicer in (None, Ellipsis): continue - # Allow single items in a dimension to be sliced - ary = np.asanyarray(slicer) - if ary.shape == (1,): - slicer = ary[0] try: int(slicer) except TypeError: @@ -117,16 +113,9 @@ def canonical_slicers(sliceobj, shape, check_inds=True): can_slicers.extend((slice(None),) * n_ellided) n_real += n_ellided continue - # int / [int] / slice indexing cases + # int / slice indexing cases dim_len = shape[n_real] n_real += 1 - - # Treat [x] as "x:x+1" - ary = np.asanyarray(slicer) - if ary.shape == (1,): - can_slicers.append(slice(int(ary[0]), int(ary[0]) + 1)) - continue - try: # test for integer indexing slicer = int(slicer) except TypeError: # should be slice object From 469376be33fdcb7b7820f74934f411099ad4fce3 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 3 Oct 2017 14:39:07 -0400 Subject: [PATCH 11/36] ENH: Add SpatialImage._spatial_dims slicer --- nibabel/minc1.py | 6 ++++++ nibabel/spatialimages.py | 3 +++ 2 files changed, 9 insertions(+) diff --git a/nibabel/minc1.py b/nibabel/minc1.py index 5eb077ada0..e3a1632199 100644 --- a/nibabel/minc1.py +++ b/nibabel/minc1.py @@ -319,6 +319,12 @@ def from_file_map(klass, file_map): data = klass.ImageArrayProxy(minc_file) return klass(data, affine, header, extra=None, file_map=file_map) + @property + def _spatial_dims(self): + if len(self.shape) > 3: + return slice(1, 4) + return slice(0, 3) + load = Minc1Image.load diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index f327ac3c60..c9f7992f94 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -325,6 +325,9 @@ class ImageDataError(Exception): class SpatialImage(DataobjImage): ''' Template class for volumetric (3D/4D) images ''' header_class = SpatialHeader + # Make a strong assumption that the first three dimensions are spatial + # Subclasses/images for which this is not true should change this attribute + _spatial_dims = slice(0, 3) class Slicer(object): def __init__(self, img): From cfa9fdde61085af33e6b7bfed158d95504c95ce7 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 3 Oct 2017 14:40:07 -0400 Subject: [PATCH 12/36] RF: Move slicer validation into image --- nibabel/spatialimages.py | 82 ++++++++++++++++++++++++++++++---------- 1 file changed, 61 insertions(+), 21 deletions(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index c9f7992f94..59986f6760 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -330,21 +330,36 @@ class SpatialImage(DataobjImage): _spatial_dims = slice(0, 3) class Slicer(object): + ''' Slicing interface that returns a new image with an updated affine + ''' def __init__(self, img): self.img = img - def __getitem__(self, idx): - idx = self._validate(idx) - dataobj = self.img.dataobj[idx] - affine = self.img.slice_affine(idx) - return self.img.__class__(dataobj, affine, self.img.header) + def __getitem__(self, slicer): + klass = self.img.__class__ + if not klass.makeable: + raise NotImplementedError( + "Cannot slice un-makeable image types") + + slicer = self.img._check_slicing(self._arr_to_slice(slicer), + self.img.shape) + dataobj = self.img.dataobj[slicer] + affine = self.img._slice_affine(slicer) + return klass(dataobj.copy(), affine, self.img.header) + + def _arr_to_slice(self, slicer): + ''' Convert single item sequence indices to slices ''' + if not isinstance(slicer, tuple): + slicer = (slicer,) - def _validate(self, idx): - idx = canonical_slicers(idx, self.img.shape) - for slicer in idx: - if isinstance(slicer, int): - raise IndexError("Cannot slice image with integer") - return idx + out = [] + for subslicer in slicer: + arr = np.asarray(subslicer) + if arr.shape == (1,): + subslicer = slice(arr[0], arr[0] + 1) + out.append(subslicer) + + return tuple(out) def __init__(self, dataobj, affine, header=None, extra=None, file_map=None): @@ -482,7 +497,32 @@ def from_image(klass, img): klass.header_class.from_header(img.header), extra=img.extra.copy()) - def slice_affine(self, idx): + def _check_slicing(self, slicer, return_spatial=False): + ''' Canonicalize slicers and check for scalar indices in spatial dims + + Parameters + ---------- + slicer : object + something that can be used to slice an array as in + ``arr[sliceobj]`` + return_spatial : bool + return only slices along spatial dimensions (x, y, z) + + Returns + ------- + slicer : object + Validated slicer object that will slice image's `dataobj` + without collapsing spatial dimensions + ''' + slicer = canonical_slicers(slicer, self.shape) + spatial_slices = slicer[self._spatial_dims] + if any(not isinstance(subslicer, (slice, None)) + for subslicer in spatial_slices): + raise IndexError("Scalar indices disallowed in spatial dimensions; " + "Use `[x]` or `x:x+1`.") + return spatial_slices if return_spatial else slicer + + def _slice_affine(self, slicer): """ Retrieve affine for current image, if sliced by a given index Applies scaling if down-sampling is applied, and adjusts the intercept @@ -490,14 +530,16 @@ def slice_affine(self, idx): Parameters ---------- - idx : numpy-compatible slice index + slicer : object + something that can be used to slice an array as in + ``arr[sliceobj]`` Returns ------- affine : (4,4) ndarray Affine with updated scale and intercept """ - idx = canonical_slicers(idx, self.shape, check_inds=False)[:3] + slicer = self._check_slicing(slicer, return_spatial=True) # Transform: # sx 0 0 tx @@ -506,14 +548,12 @@ def slice_affine(self, idx): # 0 0 0 1 transform = np.eye(4, dtype=int) - for i, slicer in enumerate(idx): - if isinstance(slicer, slice): - if slicer.step == 0: + for i, subslicer in enumerate(slicer): + if isinstance(subslicer, slice): + if subslicer.step == 0: raise ValueError("slice step cannot be 0") - transform[i, i] = slicer.step if slicer.step is not None else 1 - transform[i, 3] = slicer.start or 0 - elif isinstance(slicer, int): - transform[i, 3] = slicer + transform[i, i] = subslicer.step if subslicer.step is not None else 1 + transform[i, 3] = subslicer.start or 0 # If slicer is None, nothing to do return self.affine.dot(transform) From 9bd4d0e5b78c719831f3e43b1d0ad561a24ffd9d Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 3 Oct 2017 14:42:38 -0400 Subject: [PATCH 13/36] RF: Remove crop_image (for now) --- nibabel/funcs.py | 66 ------------------------------------- nibabel/tests/test_funcs.py | 21 +----------- 2 files changed, 1 insertion(+), 86 deletions(-) diff --git a/nibabel/funcs.py b/nibabel/funcs.py index 2258d9fe8f..240b20f802 100644 --- a/nibabel/funcs.py +++ b/nibabel/funcs.py @@ -219,69 +219,3 @@ def _aff_is_diag(aff): ''' Utility function returning True if affine is nearly diagonal ''' rzs_aff = aff[:3, :3] return np.allclose(rzs_aff, np.diag(np.diag(rzs_aff))) - - -def crop_image(img, bounds, margin=0): - ''' Crop ``img`` to specified bounds - - The image is cropped in the current orientation; no rotation or resampling - is performed. - The affine matrix is updated with the new intercept, so that all values are - found at the same RAS locations. - - Parameters - ---------- - img : ``spatialimage`` - Image to be cropped - bounds : (3, 2) array-like - Minimum and maximum (inclusive) voxel indices, in each dimension, to be - included in cropped image - margin : int, optional - Margin, in voxels, to include on each side of cropped image - - Returns - ------- - cropped_img : ``spatialimage`` - Version of `img` with cropped data array and updated affine matrix - ''' - - shape = np.reshape(img.shape[:3], (3, 1)) - bounds = np.asanyarray(bounds) - if bounds.shape != (3, 2): - raise ValueError("bounds must be interpretable as a 3x2 array") - elif np.any(bounds > shape): - raise ValueError("bounds must not exceed image dimensions") - - # Permit negative bounds - if np.any(bounds < 0): - bounds = (bounds + shape) % shape - - if np.any(bounds < 0): - raise ValueError("negative bounds may not exceed image dimensions") - elif np.any(bounds[:, 0] > bounds[:, 1]): - raise ValueError("degenerate (0 width) crops are not permitted") - - # Add margin in all directions - bounds += np.array([-margin, margin]) - - # Set min/max - bounds[bounds < 0] = 0 - over = bounds[:, 1] > shape.reshape(-1) - 1 - bounds[over, 1] = shape[over, 0] - 1 - - # Include upper bounds - bounds[:, 1] += 1 - - # Return original image if no cropping required - if np.array_equal(bounds, np.hstack(([[0], [0], [0]], shape))): - return img - - x, y, z = bounds - new_origin = np.vstack((bounds[:, [0]], 1)) - - bounded_data = img.get_data()[x[0]:x[1], y[0]:y[1], z[0]:z[1]] - - new_aff = img.affine.copy() - new_aff[:, [3]] = img.affine.dot(new_origin) - - return img.__class__(bounded_data, new_aff, img.header) diff --git a/nibabel/tests/test_funcs.py b/nibabel/tests/test_funcs.py index 91fd83aff0..6032c08672 100644 --- a/nibabel/tests/test_funcs.py +++ b/nibabel/tests/test_funcs.py @@ -11,7 +11,7 @@ import numpy as np -from ..funcs import concat_images, as_closest_canonical, OrientationError, crop_image +from ..funcs import concat_images, as_closest_canonical, OrientationError from ..analyze import AnalyzeImage from ..nifti1 import Nifti1Image from ..loadsave import save @@ -193,22 +193,3 @@ def test_closest_canonical(): img.header.set_dim_info(None, None, 2) xyz_img = as_closest_canonical(img) assert_true(xyz_img.header.get_dim_info() == (None, None, 1)) - - -def test_crop_image(): - # Use 32-bit data so that the AnalyzeImage class doesn't complain - arr = np.arange(60).reshape((5, 3, 4, 1)).astype(np.int32) - - img = AnalyzeImage(arr, np.eye(4)) - - cropped_img = crop_image(img, [[1, 3], [1, 1], [1, 2]]) - assert_equal(cropped_img.shape, (3, 1, 2, 1)) - assert_array_equal(cropped_img.affine, [[1, 0, 0, 1], - [0, 1, 0, 1], - [0, 0, 1, 1], - [0, 0, 0, 1]]) - - cropped_img = crop_image(img, [[1, 3], [1, 1], [1, 2]], margin=1) - assert_equal(cropped_img.shape, (5, 3, 4, 1)) - assert_array_equal(cropped_img.affine, img.affine) - assert_true(cropped_img is img) From c686a5a74ee6753d6f692bf1aa7acddf79bc6c5b Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 4 Oct 2017 10:02:02 -0400 Subject: [PATCH 14/36] FIX: Type check, _check_slicing arguments --- nibabel/spatialimages.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 59986f6760..5be3a4afd5 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -341,8 +341,7 @@ def __getitem__(self, slicer): raise NotImplementedError( "Cannot slice un-makeable image types") - slicer = self.img._check_slicing(self._arr_to_slice(slicer), - self.img.shape) + slicer = self.img._check_slicing(self._arr_to_slice(slicer)) dataobj = self.img.dataobj[slicer] affine = self.img._slice_affine(slicer) return klass(dataobj.copy(), affine, self.img.header) @@ -516,10 +515,12 @@ def _check_slicing(self, slicer, return_spatial=False): ''' slicer = canonical_slicers(slicer, self.shape) spatial_slices = slicer[self._spatial_dims] - if any(not isinstance(subslicer, (slice, None)) - for subslicer in spatial_slices): - raise IndexError("Scalar indices disallowed in spatial dimensions; " - "Use `[x]` or `x:x+1`.") + for subslicer in spatial_slices: + if subslicer is None: + raise IndexError("New axis not permitted in spatial dimensions") + elif isinstance(subslicer, int): + raise IndexError("Scalar indices disallowed in spatial dimensions; " + "Use `[x]` or `x:x+1`.") return spatial_slices if return_spatial else slicer def _slice_affine(self, slicer): From dfdbbb7ecd555722ff0e24e2de6df4457f1bae66 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 4 Oct 2017 10:02:39 -0400 Subject: [PATCH 15/36] ENH: Handle new axes in Minc --- nibabel/minc1.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/nibabel/minc1.py b/nibabel/minc1.py index e3a1632199..06afad2dea 100644 --- a/nibabel/minc1.py +++ b/nibabel/minc1.py @@ -325,6 +325,42 @@ def _spatial_dims(self): return slice(1, 4) return slice(0, 3) + def _check_slicing(self, slicer, return_spatial=False): + ''' Canonicalize slicers and check for scalar indices in spatial dims + + Parameters + ---------- + slicer : object + something that can be used to slice an array as in + ``arr[sliceobj]`` + return_spatial : bool + return only slices along spatial dimensions (x, y, z) + + Returns + ------- + slicer : object + Validated slicer object that will slice image's `dataobj` + without collapsing spatial dimensions + ''' + try: + all_slices = super(Minc1Image, self)._check_slicing(slicer, False) + sp_dims = self._spatial_dims + except IndexError: + # Prepending a new axis for 3D images is valid in Minc + if slicer[0] is None and self._spatial_dims == slice(0, 3): + all_slices = (None,) + super(Minc1Image, self)._check_slicing(slicer[1:], False) + sp_dims = slice(1, 4) + else: + raise + # Added complications of first axis being time + if self._spatial_dims == slice(1, 4) and all_slices[0] is None: + raise IndexError("New temporal axis is not permitted in 4D Minc images") + elif (self._spatial_dims == slice(0, 3) and len(all_slices) > 3 and + all_slices[0] is not None): + raise IndexError("New axes cannot be added to 3D Minc image " + "without new temporal axis (first dimension)") + return all_slices[sp_dims] if return_spatial else all_slices + load = Minc1Image.load From 2d8c7b58fae069d12cccaaeb7c76ec7994465bb4 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 4 Oct 2017 16:31:05 -0400 Subject: [PATCH 16/36] TEST: Image slicing --- nibabel/tests/test_spatialimages.py | 107 ++++++++++++++++++++++++++++ 1 file changed, 107 insertions(+) diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index bd8b834b84..4164ad84a4 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -411,6 +411,113 @@ def test_get_data(self): assert_false(rt_img.get_data() is out_data) assert_array_equal(rt_img.get_data(), in_data) + def test_slicer(self): + img_klass = self.image_class + in_data_template = np.arange(240, dtype=np.int16) + base_affine = np.eye(4) + t_axis = None + for dshape in ((4, 5, 6, 2), # Time series + (8, 5, 6)): # Volume + in_data = in_data_template.copy().reshape(dshape) + img = img_klass(in_data, base_affine.copy()) + + # Detect time axis on first loop (4D image) + if t_axis is None: + t_axis = 3 if img._spatial_dims.start == 0 else 0 + + assert_true(hasattr(img.slicer, '__getitem__')) + if not img_klass.makeable: + assert_raises(NotImplementedError, img.slicer[:]) + continue + # Note spatial zooms are always first 3, even when + spatial_zooms = img.header.get_zooms()[:3] + + # Down-sample with [::2, ::2, ::2] along spatial dimensions + sliceobj = [slice(None)] * len(dshape) + sliceobj[img._spatial_dims] = [slice(None, None, 2)] * 3 + downsampled_img = img.slicer[tuple(sliceobj)] + assert_array_equal(downsampled_img.header.get_zooms()[:3], + np.array(spatial_zooms) * 2) + + # Check newaxis errors + if t_axis == 3: + with assert_raises(IndexError): + img.slicer[None] + elif len(img.shape) == 4: + with assert_raises(IndexError): + img.slicer[None] + # 4D Minc to 3D + assert_equal(img.slicer[0].shape, img.shape[1:]) + else: + # 3D Minc to 4D + assert_equal(img.slicer[None].shape, (1,) + img.shape) + # Axes 1 and 2 are always spatial + with assert_raises(IndexError): + img.slicer[:, None] + with assert_raises(IndexError): + img.slicer[:, :, None] + if t_axis == 0: + with assert_raises(IndexError): + img.slicer[:, :, :, None] + elif len(img.shape) == 4: + # Reorder non-spatial axes + assert_equal(img.slicer[:, :, :, None].shape, img.shape[:3] + (1,) + img.shape[3:]) + else: + # 3D Analyze/NIfTI/MGH to 4D + assert_equal(img.slicer[:, :, :, None].shape, img.shape + (1,)) + if len(img.shape) == 3: + # Slices exceed dimensions + with assert_raises(IndexError): + img.slicer[:, :, :, :, None] + else: + assert_equal(img.slicer[:, :, :, :, None].shape, img.shape + (1,)) + + # Crop by one voxel in each dimension + if len(img.shape) == 3 or t_axis == 3: + sliced_i = img.slicer[1:] + sliced_j = img.slicer[:, 1:] + sliced_k = img.slicer[:, :, 1:] + sliced_ijk = img.slicer[1:, 1:, 1:] + else: + # 4D Minc + sliced_i = img.slicer[:, 1:] + sliced_j = img.slicer[:, :, 1:] + sliced_k = img.slicer[:, :, :, 1:] + sliced_ijk = img.slicer[:, 1:, 1:, 1:] + + # No scaling change + assert_array_equal(sliced_i.affine[:3, :3], img.affine[:3, :3]) + assert_array_equal(sliced_j.affine[:3, :3], img.affine[:3, :3]) + assert_array_equal(sliced_k.affine[:3, :3], img.affine[:3, :3]) + assert_array_equal(sliced_ijk.affine[:3, :3], img.affine[:3, :3]) + # Translation + assert_array_equal(sliced_i.affine[:, 3], [1, 0, 0, 1]) + assert_array_equal(sliced_j.affine[:, 3], [0, 1, 0, 1]) + assert_array_equal(sliced_k.affine[:, 3], [0, 0, 1, 1]) + assert_array_equal(sliced_ijk.affine[:, 3], [1, 1, 1, 1]) + + # No change to affines with upper-bound slices + assert_array_equal(img.slicer[:1, :1, :1].affine, img.affine) + + # Check data is consistent with slicing numpy arrays + slice_elems = (None, Ellipsis, 0, 1, -1, slice(None), slice(1), + slice(-1), slice(1, -1)) + for n_elems in range(6): + for _ in range(10): + sliceobj = tuple( + np.random.choice(slice_elems, n_elems).tolist()) + try: + sliced_img = img.slicer[sliceobj] + except IndexError: + # Only checking valid slices + pass + else: + sliced_data = in_data[sliceobj] + assert_array_equal(sliced_data, sliced_img.get_data()) + assert_array_equal(sliced_data, sliced_img.dataobj) + assert_array_equal(sliced_data, img.dataobj[sliceobj]) + assert_array_equal(sliced_data, img.get_data()[sliceobj]) + def test_api_deprecations(self): class FakeImage(self.image_class): From 07187016970bf2aade0500c3e5c006d333566fda Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 4 Oct 2017 17:39:27 -0400 Subject: [PATCH 17/36] Translate ValueErrors to IndexErrors when used --- nibabel/spatialimages.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 5be3a4afd5..02f9a58b91 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -341,7 +341,10 @@ def __getitem__(self, slicer): raise NotImplementedError( "Cannot slice un-makeable image types") - slicer = self.img._check_slicing(self._arr_to_slice(slicer)) + try: + slicer = self.img._check_slicing(self._arr_to_slice(slicer)) + except ValueError as err: + raise IndexError(*err.args) dataobj = self.img.dataobj[slicer] affine = self.img._slice_affine(slicer) return klass(dataobj.copy(), affine, self.img.header) From b3b4a22e81a6400dd4a1cb06714d301cfb8dacf6 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 4 Oct 2017 19:53:39 -0400 Subject: [PATCH 18/36] TEST: Handle MGHFormat (max: 4D) correctly --- nibabel/tests/test_spatialimages.py | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index 4164ad84a4..5aab7737d1 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -439,6 +439,9 @@ def test_slicer(self): assert_array_equal(downsampled_img.header.get_zooms()[:3], np.array(spatial_zooms) * 2) + max4d = (hasattr(img.header, '_header_data') and + 'dims' in img.header._header_data.dtype.fields and + img.header._header_data['dims'].shape == (4,)) # Check newaxis errors if t_axis == 3: with assert_raises(IndexError): @@ -460,8 +463,13 @@ def test_slicer(self): with assert_raises(IndexError): img.slicer[:, :, :, None] elif len(img.shape) == 4: - # Reorder non-spatial axes - assert_equal(img.slicer[:, :, :, None].shape, img.shape[:3] + (1,) + img.shape[3:]) + if max4d: + with assert_raises(ValueError): + img.slicer[:, :, :, None] + else: + # Reorder non-spatial axes + assert_equal(img.slicer[:, :, :, None].shape, + img.shape[:3] + (1,) + img.shape[3:]) else: # 3D Analyze/NIfTI/MGH to 4D assert_equal(img.slicer[:, :, :, None].shape, img.shape + (1,)) @@ -469,8 +477,12 @@ def test_slicer(self): # Slices exceed dimensions with assert_raises(IndexError): img.slicer[:, :, :, :, None] + elif max4d: + with assert_raises(ValueError): + img.slicer[:, :, :, :, None] else: - assert_equal(img.slicer[:, :, :, :, None].shape, img.shape + (1,)) + assert_equal(img.slicer[:, :, :, :, None].shape, + img.shape + (1,)) # Crop by one voxel in each dimension if len(img.shape) == 3 or t_axis == 3: @@ -508,7 +520,7 @@ def test_slicer(self): np.random.choice(slice_elems, n_elems).tolist()) try: sliced_img = img.slicer[sliceobj] - except IndexError: + except (IndexError, ValueError): # Only checking valid slices pass else: From 67efd1c4b03f8fffa143cfd37cacb59872666874 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 5 Oct 2017 23:16:12 -0400 Subject: [PATCH 19/36] TEST: Test single-element list indices --- nibabel/tests/test_spatialimages.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index 5aab7737d1..93d09901e0 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -512,8 +512,8 @@ def test_slicer(self): assert_array_equal(img.slicer[:1, :1, :1].affine, img.affine) # Check data is consistent with slicing numpy arrays - slice_elems = (None, Ellipsis, 0, 1, -1, slice(None), slice(1), - slice(-1), slice(1, -1)) + slice_elems = (None, Ellipsis, 0, 1, -1, [0], [1], [-1], + slice(None), slice(1), slice(-1), slice(1, -1)) for n_elems in range(6): for _ in range(10): sliceobj = tuple( From 526d81c881530bfbe556dff91efffbd629c5302b Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 5 Oct 2017 23:23:47 -0400 Subject: [PATCH 20/36] RF: All SpatialImages are makeable --- nibabel/spatialimages.py | 7 +------ nibabel/tests/test_spatialimages.py | 4 +--- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 02f9a58b91..09655e1cd2 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -336,18 +336,13 @@ def __init__(self, img): self.img = img def __getitem__(self, slicer): - klass = self.img.__class__ - if not klass.makeable: - raise NotImplementedError( - "Cannot slice un-makeable image types") - try: slicer = self.img._check_slicing(self._arr_to_slice(slicer)) except ValueError as err: raise IndexError(*err.args) dataobj = self.img.dataobj[slicer] affine = self.img._slice_affine(slicer) - return klass(dataobj.copy(), affine, self.img.header) + return self.img.__class__(dataobj.copy(), affine, self.img.header) def _arr_to_slice(self, slicer): ''' Convert single item sequence indices to slices ''' diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index 93d09901e0..d4638beb25 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -426,9 +426,7 @@ def test_slicer(self): t_axis = 3 if img._spatial_dims.start == 0 else 0 assert_true(hasattr(img.slicer, '__getitem__')) - if not img_klass.makeable: - assert_raises(NotImplementedError, img.slicer[:]) - continue + # Note spatial zooms are always first 3, even when spatial_zooms = img.header.get_zooms()[:3] From fcf55b58c2f2458375bcd8b23e2a847025077c78 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 5 Oct 2017 23:23:59 -0400 Subject: [PATCH 21/36] TEST: Check step = 0 case --- nibabel/tests/test_spatialimages.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index d4638beb25..530e3ae8a1 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -509,6 +509,12 @@ def test_slicer(self): # No change to affines with upper-bound slices assert_array_equal(img.slicer[:1, :1, :1].affine, img.affine) + # Yell about step = 0 + with assert_raises(ValueError): + img.slicer[:, ::0] + with assert_raises(ValueError): + img._slice_affine((slice(None), slice(None, None, 0))) + # Check data is consistent with slicing numpy arrays slice_elems = (None, Ellipsis, 0, 1, -1, [0], [1], [-1], slice(None), slice(1), slice(-1), slice(1, -1)) From 370d0e1a1c0c7f834ff24332663cd72c5867a933 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Sun, 8 Oct 2017 12:36:31 -0400 Subject: [PATCH 22/36] ENH: Drop fancy indexing --- nibabel/minc1.py | 1 + nibabel/spatialimages.py | 16 +--------------- 2 files changed, 2 insertions(+), 15 deletions(-) diff --git a/nibabel/minc1.py b/nibabel/minc1.py index 06afad2dea..a947949053 100644 --- a/nibabel/minc1.py +++ b/nibabel/minc1.py @@ -342,6 +342,7 @@ def _check_slicing(self, slicer, return_spatial=False): Validated slicer object that will slice image's `dataobj` without collapsing spatial dimensions ''' + slicer = canonical_slicers(slicer, self.shape) try: all_slices = super(Minc1Image, self)._check_slicing(slicer, False) sp_dims = self._spatial_dims diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 09655e1cd2..8120a95276 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -337,27 +337,13 @@ def __init__(self, img): def __getitem__(self, slicer): try: - slicer = self.img._check_slicing(self._arr_to_slice(slicer)) + slicer = self.img._check_slicing(slicer) except ValueError as err: raise IndexError(*err.args) dataobj = self.img.dataobj[slicer] affine = self.img._slice_affine(slicer) return self.img.__class__(dataobj.copy(), affine, self.img.header) - def _arr_to_slice(self, slicer): - ''' Convert single item sequence indices to slices ''' - if not isinstance(slicer, tuple): - slicer = (slicer,) - - out = [] - for subslicer in slicer: - arr = np.asarray(subslicer) - if arr.shape == (1,): - subslicer = slice(arr[0], arr[0] + 1) - out.append(subslicer) - - return tuple(out) - def __init__(self, dataobj, affine, header=None, extra=None, file_map=None): ''' Initialize image From 2d6547c4fbd7115cf11df39cede48cbdb754ea7a Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Sun, 8 Oct 2017 12:37:21 -0400 Subject: [PATCH 23/36] TEST: Test single-slice and fancy indexing raise errors --- nibabel/tests/test_spatialimages.py | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index 530e3ae8a1..211b5f70e4 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -440,10 +440,12 @@ def test_slicer(self): max4d = (hasattr(img.header, '_header_data') and 'dims' in img.header._header_data.dtype.fields and img.header._header_data['dims'].shape == (4,)) - # Check newaxis errors + # Check newaxis and single-slice errors if t_axis == 3: with assert_raises(IndexError): img.slicer[None] + with assert_raises(IndexError): + img.slicer[0] elif len(img.shape) == 4: with assert_raises(IndexError): img.slicer[None] @@ -455,11 +457,17 @@ def test_slicer(self): # Axes 1 and 2 are always spatial with assert_raises(IndexError): img.slicer[:, None] + with assert_raises(IndexError): + img.slicer[:, 0] with assert_raises(IndexError): img.slicer[:, :, None] + with assert_raises(IndexError): + img.slicer[:, :, 0] if t_axis == 0: with assert_raises(IndexError): img.slicer[:, :, :, None] + with assert_raises(IndexError): + img.slicer[:, :, :, 0] elif len(img.shape) == 4: if max4d: with assert_raises(ValueError): @@ -468,6 +476,9 @@ def test_slicer(self): # Reorder non-spatial axes assert_equal(img.slicer[:, :, :, None].shape, img.shape[:3] + (1,) + img.shape[3:]) + # 4D to 3D using ellipsis or slices + assert_equal(img.slicer[..., 0].shape, img.shape[:-1]) + assert_equal(img.slicer[:, :, :, 0].shape, img.shape[:-1]) else: # 3D Analyze/NIfTI/MGH to 4D assert_equal(img.slicer[:, :, :, None].shape, img.shape + (1,)) @@ -515,11 +526,19 @@ def test_slicer(self): with assert_raises(ValueError): img._slice_affine((slice(None), slice(None, None, 0))) + # No fancy indexing + with assert_raises(IndexError): + img.slicer[[0]] + with assert_raises(IndexError): + img.slicer[[-1]] + with assert_raises(IndexError): + img.slicer[[0], [-1]] + # Check data is consistent with slicing numpy arrays slice_elems = (None, Ellipsis, 0, 1, -1, [0], [1], [-1], slice(None), slice(1), slice(-1), slice(1, -1)) for n_elems in range(6): - for _ in range(10): + for _ in range(1 if n_elems == 0 else 10): sliceobj = tuple( np.random.choice(slice_elems, n_elems).tolist()) try: From e8096364c7bb66c6879cdb6c429de0a9e447abe9 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 18 Dec 2017 12:33:36 -0500 Subject: [PATCH 24/36] DOC: Changelog entry --- Changelog | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/Changelog b/Changelog index 4c65cf0bfe..b96ddd40a7 100644 --- a/Changelog +++ b/Changelog @@ -24,6 +24,27 @@ Gerhard (SG) and Eric Larson (EL). References like "pr/298" refer to github pull request numbers. +Upcoming Release +================ + +New features +------------ +* Image slicing for SpatialImages (pr/550) (CM) + +Enhancements +------------ +* Simplfiy MGHImage and add footer fields (pr/569) (CM, reviewed by MB) + +Bug fixes +--------- + +Maintenance +----------- + +API changes and deprecations +---------------------------- + + 2.2.1 (Wednesday 22 November 2017) ================================== From 07bbb250928bb521eace305745dcecd0e1d5e272 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 19 Feb 2018 10:43:00 -0500 Subject: [PATCH 25/36] TEST: Do not use deprecated _header_data --- nibabel/tests/test_spatialimages.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index 211b5f70e4..75ca523940 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -437,9 +437,9 @@ def test_slicer(self): assert_array_equal(downsampled_img.header.get_zooms()[:3], np.array(spatial_zooms) * 2) - max4d = (hasattr(img.header, '_header_data') and - 'dims' in img.header._header_data.dtype.fields and - img.header._header_data['dims'].shape == (4,)) + max4d = (hasattr(img.header, '_structarr') and + 'dims' in img.header._structarr.dtype.fields and + img.header._structarr['dims'].shape == (4,)) # Check newaxis and single-slice errors if t_axis == 3: with assert_raises(IndexError): From 1e15c9852636399b1f8e087883069478e449ccdf Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 19 Feb 2018 13:14:43 -0500 Subject: [PATCH 26/36] DOC: Link to aliasing wiki article [skip ci] --- nibabel/spatialimages.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 8120a95276..d9108a63b6 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -549,7 +549,7 @@ def slicer(self): The image is resliced in the current orientation; no rotation or resampling is performed, and no attempt is made to filter the image - to avoid aliasing. + to avoid `aliasing`_. The affine matrix is updated with the new intercept (and scales, if down-sampling is used), so that all values are found at the same RAS @@ -558,6 +558,8 @@ def slicer(self): Slicing may include non-spatial dimensions. However, this method does not currently adjust the repetition time in the image header. + + .. _aliasing: https://en.wikipedia.org/wiki/Aliasing """ return self.Slicer(self) From fbd9e6e701f395ee8e9bc34221111c8d259dc2df Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 19 Feb 2018 13:40:13 -0500 Subject: [PATCH 27/36] ENH: Improve exception message --- nibabel/spatialimages.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index d9108a63b6..d15f270c42 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -570,9 +570,11 @@ def __getitem__(self, idx): Use the slicer attribute to perform cropping and subsampling at your own risk. ''' - raise TypeError("Cannot slice image objects; consider slicing image " - "array data with `img.dataobj[slice]` or " - "`img.get_data()[slice]`") + raise TypeError( + "Cannot slice image objects; consider using `img.slicer[slice]` " + "to generate a sliced image (see documentation for caveats) or " + "slicing image array data with `img.dataobj[slice]` or " + "`img.get_data()[slice]`") def orthoview(self): """Plot the image using OrthoSlicer3D From f57bfd54d76f2bfb736dac9c99eb23894b975e42 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 19 Feb 2018 14:51:02 -0500 Subject: [PATCH 28/36] RF: Specify spatial dims in each class --- nibabel/analyze.py | 1 + nibabel/ecat.py | 1 + nibabel/freesurfer/mghformat.py | 1 + nibabel/parrec.py | 1 + nibabel/spatialimages.py | 4 +--- 5 files changed, 5 insertions(+), 3 deletions(-) diff --git a/nibabel/analyze.py b/nibabel/analyze.py index 5a43d16773..887523ed95 100644 --- a/nibabel/analyze.py +++ b/nibabel/analyze.py @@ -911,6 +911,7 @@ class AnalyzeImage(SpatialImage): files_types = (('image', '.img'), ('header', '.hdr')) valid_exts = ('.img', '.hdr') _compressed_suffixes = ('.gz', '.bz2') + _spatial_dims = slice(0, 3) makeable = True rw = True diff --git a/nibabel/ecat.py b/nibabel/ecat.py index 3c0957e11d..84cdc6c7c0 100644 --- a/nibabel/ecat.py +++ b/nibabel/ecat.py @@ -740,6 +740,7 @@ class EcatImage(SpatialImage): valid_exts = ('.v',) _subheader = EcatSubHeader files_types = (('image', '.v'), ('header', '.v')) + _spatial_dims = slice(0, 3) ImageArrayProxy = EcatImageArrayProxy diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index 927d6126c0..e444291e55 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -507,6 +507,7 @@ class MGHImage(SpatialImage): ImageOpener.compress_ext_map['.mgz'] = ImageOpener.gz_def files_types = (('image', '.mgh'),) _compressed_suffixes = () + _spatial_dims = slice(0, 3) makeable = True rw = True diff --git a/nibabel/parrec.py b/nibabel/parrec.py index 5fd460b4e1..57ec1b3cba 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -1222,6 +1222,7 @@ class PARRECImage(SpatialImage): header_class = PARRECHeader valid_exts = ('.rec', '.par') files_types = (('image', '.rec'), ('header', '.par')) + _spatial_dims = slice(0, 3) makeable = False rw = False diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index d15f270c42..1396dcf753 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -325,9 +325,7 @@ class ImageDataError(Exception): class SpatialImage(DataobjImage): ''' Template class for volumetric (3D/4D) images ''' header_class = SpatialHeader - # Make a strong assumption that the first three dimensions are spatial - # Subclasses/images for which this is not true should change this attribute - _spatial_dims = slice(0, 3) + _spatial_dims = None class Slicer(object): ''' Slicing interface that returns a new image with an updated affine From 8a6150cf7f273de31a87f06ea6e6be07c4345a56 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Sun, 4 Mar 2018 13:29:59 -0500 Subject: [PATCH 29/36] TEST: Check for _spatial_dims --- nibabel/tests/test_spatialimages.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index 75ca523940..dfc0dde08e 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -422,7 +422,7 @@ def test_slicer(self): img = img_klass(in_data, base_affine.copy()) # Detect time axis on first loop (4D image) - if t_axis is None: + if t_axis is None and img._spatial_dims is not None: t_axis = 3 if img._spatial_dims.start == 0 else 0 assert_true(hasattr(img.slicer, '__getitem__')) @@ -451,7 +451,7 @@ def test_slicer(self): img.slicer[None] # 4D Minc to 3D assert_equal(img.slicer[0].shape, img.shape[1:]) - else: + elif t_axis is not None: # 3D Minc to 4D assert_equal(img.slicer[None].shape, (1,) + img.shape) # Axes 1 and 2 are always spatial @@ -479,7 +479,7 @@ def test_slicer(self): # 4D to 3D using ellipsis or slices assert_equal(img.slicer[..., 0].shape, img.shape[:-1]) assert_equal(img.slicer[:, :, :, 0].shape, img.shape[:-1]) - else: + elif t_axis is not None: # 3D Analyze/NIfTI/MGH to 4D assert_equal(img.slicer[:, :, :, None].shape, img.shape + (1,)) if len(img.shape) == 3: @@ -499,7 +499,7 @@ def test_slicer(self): sliced_j = img.slicer[:, 1:] sliced_k = img.slicer[:, :, 1:] sliced_ijk = img.slicer[1:, 1:, 1:] - else: + elif t_axis is not None: # 4D Minc sliced_i = img.slicer[:, 1:] sliced_j = img.slicer[:, :, 1:] From ea053c88625cf400549e7971c3598f543361089e Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Fri, 6 Apr 2018 23:07:05 -0400 Subject: [PATCH 30/36] TEST: Minor fixes --- nibabel/tests/test_spatialimages.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index dfc0dde08e..067a671b99 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -385,9 +385,10 @@ def test_get_data(self): img[0, 0, 0] # Make sure the right message gets raised: assert_equal(str(exception_manager.exception), - ("Cannot slice image objects; consider slicing image " - "array data with `img.dataobj[slice]` or " - "`img.get_data()[slice]`")) + "Cannot slice image objects; consider using " + "`img.slicer[slice]` to generate a sliced image (see " + "documentation for caveats) or slicing image array data " + "with `img.dataobj[slice]` or `img.get_data()[slice]`") assert_true(in_data is img.dataobj) out_data = img.get_data() assert_true(in_data is out_data) @@ -432,7 +433,8 @@ def test_slicer(self): # Down-sample with [::2, ::2, ::2] along spatial dimensions sliceobj = [slice(None)] * len(dshape) - sliceobj[img._spatial_dims] = [slice(None, None, 2)] * 3 + if img._spatial_dims is not None: + sliceobj[img._spatial_dims] = [slice(None, None, 2)] * 3 downsampled_img = img.slicer[tuple(sliceobj)] assert_array_equal(downsampled_img.header.get_zooms()[:3], np.array(spatial_zooms) * 2) From 6f6e38a80cedb2c5d44062671f412fc1f311e354 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Sat, 5 May 2018 23:31:52 -0400 Subject: [PATCH 31/36] RF: Switch to spatial_axes_first --- nibabel/minc1.py | 43 ---------------------------------------- nibabel/parrec.py | 1 - nibabel/spatialimages.py | 7 +++++-- 3 files changed, 5 insertions(+), 46 deletions(-) diff --git a/nibabel/minc1.py b/nibabel/minc1.py index a947949053..5eb077ada0 100644 --- a/nibabel/minc1.py +++ b/nibabel/minc1.py @@ -319,49 +319,6 @@ def from_file_map(klass, file_map): data = klass.ImageArrayProxy(minc_file) return klass(data, affine, header, extra=None, file_map=file_map) - @property - def _spatial_dims(self): - if len(self.shape) > 3: - return slice(1, 4) - return slice(0, 3) - - def _check_slicing(self, slicer, return_spatial=False): - ''' Canonicalize slicers and check for scalar indices in spatial dims - - Parameters - ---------- - slicer : object - something that can be used to slice an array as in - ``arr[sliceobj]`` - return_spatial : bool - return only slices along spatial dimensions (x, y, z) - - Returns - ------- - slicer : object - Validated slicer object that will slice image's `dataobj` - without collapsing spatial dimensions - ''' - slicer = canonical_slicers(slicer, self.shape) - try: - all_slices = super(Minc1Image, self)._check_slicing(slicer, False) - sp_dims = self._spatial_dims - except IndexError: - # Prepending a new axis for 3D images is valid in Minc - if slicer[0] is None and self._spatial_dims == slice(0, 3): - all_slices = (None,) + super(Minc1Image, self)._check_slicing(slicer[1:], False) - sp_dims = slice(1, 4) - else: - raise - # Added complications of first axis being time - if self._spatial_dims == slice(1, 4) and all_slices[0] is None: - raise IndexError("New temporal axis is not permitted in 4D Minc images") - elif (self._spatial_dims == slice(0, 3) and len(all_slices) > 3 and - all_slices[0] is not None): - raise IndexError("New axes cannot be added to 3D Minc image " - "without new temporal axis (first dimension)") - return all_slices[sp_dims] if return_spatial else all_slices - load = Minc1Image.load diff --git a/nibabel/parrec.py b/nibabel/parrec.py index 57ec1b3cba..5fd460b4e1 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -1222,7 +1222,6 @@ class PARRECImage(SpatialImage): header_class = PARRECHeader valid_exts = ('.rec', '.par') files_types = (('image', '.rec'), ('header', '.par')) - _spatial_dims = slice(0, 3) makeable = False rw = False diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 1396dcf753..8fb99c64d0 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -325,7 +325,6 @@ class ImageDataError(Exception): class SpatialImage(DataobjImage): ''' Template class for volumetric (3D/4D) images ''' header_class = SpatialHeader - _spatial_dims = None class Slicer(object): ''' Slicing interface that returns a new image with an updated affine @@ -496,7 +495,7 @@ def _check_slicing(self, slicer, return_spatial=False): without collapsing spatial dimensions ''' slicer = canonical_slicers(slicer, self.shape) - spatial_slices = slicer[self._spatial_dims] + spatial_slices = slicer[:3] for subslicer in spatial_slices: if subslicer is None: raise IndexError("New axis not permitted in spatial dimensions") @@ -559,6 +558,10 @@ def slicer(self): .. _aliasing: https://en.wikipedia.org/wiki/Aliasing """ + from .imageclasses import spatial_axes_first + if not spatial_axes_first(self): + raise ValueError("Cannot predict position of spatial axes for " + "Image type " + self.__class__.__name__) return self.Slicer(self) From b9201cfe15cbc9fd832b4fee118be58375a874e7 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 9 May 2018 21:57:03 -0400 Subject: [PATCH 32/36] RF/TEST: Purge _spatial_dims and related tests --- nibabel/analyze.py | 1 - nibabel/ecat.py | 1 - nibabel/freesurfer/mghformat.py | 1 - nibabel/tests/test_spatialimages.py | 54 +++++++++-------------------- 4 files changed, 17 insertions(+), 40 deletions(-) diff --git a/nibabel/analyze.py b/nibabel/analyze.py index 887523ed95..5a43d16773 100644 --- a/nibabel/analyze.py +++ b/nibabel/analyze.py @@ -911,7 +911,6 @@ class AnalyzeImage(SpatialImage): files_types = (('image', '.img'), ('header', '.hdr')) valid_exts = ('.img', '.hdr') _compressed_suffixes = ('.gz', '.bz2') - _spatial_dims = slice(0, 3) makeable = True rw = True diff --git a/nibabel/ecat.py b/nibabel/ecat.py index 84cdc6c7c0..3c0957e11d 100644 --- a/nibabel/ecat.py +++ b/nibabel/ecat.py @@ -740,7 +740,6 @@ class EcatImage(SpatialImage): valid_exts = ('.v',) _subheader = EcatSubHeader files_types = (('image', '.v'), ('header', '.v')) - _spatial_dims = slice(0, 3) ImageArrayProxy = EcatImageArrayProxy diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index e444291e55..927d6126c0 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -507,7 +507,6 @@ class MGHImage(SpatialImage): ImageOpener.compress_ext_map['.mgz'] = ImageOpener.gz_def files_types = (('image', '.mgh'),) _compressed_suffixes = () - _spatial_dims = slice(0, 3) makeable = True rw = True diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index 067a671b99..8d8ad2fb81 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -17,6 +17,7 @@ from io import BytesIO from ..spatialimages import (SpatialHeader, SpatialImage, HeaderDataError, Header, ImageDataError) +from ..imageclasses import spatial_axes_first from unittest import TestCase from nose.tools import (assert_true, assert_false, assert_equal, @@ -422,9 +423,10 @@ def test_slicer(self): in_data = in_data_template.copy().reshape(dshape) img = img_klass(in_data, base_affine.copy()) - # Detect time axis on first loop (4D image) - if t_axis is None and img._spatial_dims is not None: - t_axis = 3 if img._spatial_dims.start == 0 else 0 + if not spatial_axes_first(img): + with assert_raises(ValueError): + img.slicer + continue assert_true(hasattr(img.slicer, '__getitem__')) @@ -432,9 +434,8 @@ def test_slicer(self): spatial_zooms = img.header.get_zooms()[:3] # Down-sample with [::2, ::2, ::2] along spatial dimensions - sliceobj = [slice(None)] * len(dshape) - if img._spatial_dims is not None: - sliceobj[img._spatial_dims] = [slice(None, None, 2)] * 3 + sliceobj = [slice(None, None, 2)] * 3 + \ + [slice(None)] * (len(dshape) - 3) downsampled_img = img.slicer[tuple(sliceobj)] assert_array_equal(downsampled_img.header.get_zooms()[:3], np.array(spatial_zooms) * 2) @@ -443,19 +444,10 @@ def test_slicer(self): 'dims' in img.header._structarr.dtype.fields and img.header._structarr['dims'].shape == (4,)) # Check newaxis and single-slice errors - if t_axis == 3: - with assert_raises(IndexError): - img.slicer[None] - with assert_raises(IndexError): - img.slicer[0] - elif len(img.shape) == 4: - with assert_raises(IndexError): - img.slicer[None] - # 4D Minc to 3D - assert_equal(img.slicer[0].shape, img.shape[1:]) - elif t_axis is not None: - # 3D Minc to 4D - assert_equal(img.slicer[None].shape, (1,) + img.shape) + with assert_raises(IndexError): + img.slicer[None] + with assert_raises(IndexError): + img.slicer[0] # Axes 1 and 2 are always spatial with assert_raises(IndexError): img.slicer[:, None] @@ -465,12 +457,7 @@ def test_slicer(self): img.slicer[:, :, None] with assert_raises(IndexError): img.slicer[:, :, 0] - if t_axis == 0: - with assert_raises(IndexError): - img.slicer[:, :, :, None] - with assert_raises(IndexError): - img.slicer[:, :, :, 0] - elif len(img.shape) == 4: + if len(img.shape) == 4: if max4d: with assert_raises(ValueError): img.slicer[:, :, :, None] @@ -481,7 +468,7 @@ def test_slicer(self): # 4D to 3D using ellipsis or slices assert_equal(img.slicer[..., 0].shape, img.shape[:-1]) assert_equal(img.slicer[:, :, :, 0].shape, img.shape[:-1]) - elif t_axis is not None: + else: # 3D Analyze/NIfTI/MGH to 4D assert_equal(img.slicer[:, :, :, None].shape, img.shape + (1,)) if len(img.shape) == 3: @@ -496,17 +483,10 @@ def test_slicer(self): img.shape + (1,)) # Crop by one voxel in each dimension - if len(img.shape) == 3 or t_axis == 3: - sliced_i = img.slicer[1:] - sliced_j = img.slicer[:, 1:] - sliced_k = img.slicer[:, :, 1:] - sliced_ijk = img.slicer[1:, 1:, 1:] - elif t_axis is not None: - # 4D Minc - sliced_i = img.slicer[:, 1:] - sliced_j = img.slicer[:, :, 1:] - sliced_k = img.slicer[:, :, :, 1:] - sliced_ijk = img.slicer[:, 1:, 1:, 1:] + sliced_i = img.slicer[1:] + sliced_j = img.slicer[:, 1:] + sliced_k = img.slicer[:, :, 1:] + sliced_ijk = img.slicer[1:, 1:, 1:] # No scaling change assert_array_equal(sliced_i.affine[:3, :3], img.affine[:3, :3]) From a1394a6b5e72c33ad943e6416f1279d881ba66fa Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Sat, 2 Jun 2018 13:48:24 -0400 Subject: [PATCH 33/36] RF: Move slicing machinery to SpatialFirstSlicer class --- nibabel/spatialimages.py | 170 ++++++++++++++++++++------------------- 1 file changed, 87 insertions(+), 83 deletions(-) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 8fb99c64d0..d14ec16f80 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -322,24 +322,95 @@ class ImageDataError(Exception): pass -class SpatialImage(DataobjImage): - ''' Template class for volumetric (3D/4D) images ''' - header_class = SpatialHeader +class SpatialFirstSlicer(object): + ''' Slicing interface that returns a new image with an updated affine + + Checks that an image's first three axes are spatial + ''' + def __init__(self, img): + from .imageclasses import spatial_axes_first + if not spatial_axes_first(img): + raise ValueError("Cannot predict position of spatial axes for " + "Image type " + img.__class__.__name__) + self.img = img + + def __getitem__(self, slicer): + try: + slicer = self.check_slicing(slicer) + except ValueError as err: + raise IndexError(*err.args) + dataobj = self.img.dataobj[slicer] + affine = self.slice_affine(slicer) + return self.img.__class__(dataobj.copy(), affine, self.img.header) + + def check_slicing(self, slicer, return_spatial=False): + ''' Canonicalize slicers and check for scalar indices in spatial dims + + Parameters + ---------- + slicer : object + something that can be used to slice an array as in + ``arr[sliceobj]`` + return_spatial : bool + return only slices along spatial dimensions (x, y, z) - class Slicer(object): - ''' Slicing interface that returns a new image with an updated affine + Returns + ------- + slicer : object + Validated slicer object that will slice image's `dataobj` + without collapsing spatial dimensions ''' - def __init__(self, img): - self.img = img + slicer = canonical_slicers(slicer, self.img.shape) + spatial_slices = slicer[:3] + for subslicer in spatial_slices: + if subslicer is None: + raise IndexError("New axis not permitted in spatial dimensions") + elif isinstance(subslicer, int): + raise IndexError("Scalar indices disallowed in spatial dimensions; " + "Use `[x]` or `x:x+1`.") + return spatial_slices if return_spatial else slicer - def __getitem__(self, slicer): - try: - slicer = self.img._check_slicing(slicer) - except ValueError as err: - raise IndexError(*err.args) - dataobj = self.img.dataobj[slicer] - affine = self.img._slice_affine(slicer) - return self.img.__class__(dataobj.copy(), affine, self.img.header) + def slice_affine(self, slicer): + """ Retrieve affine for current image, if sliced by a given index + + Applies scaling if down-sampling is applied, and adjusts the intercept + to account for any cropping. + + Parameters + ---------- + slicer : object + something that can be used to slice an array as in + ``arr[sliceobj]`` + + Returns + ------- + affine : (4,4) ndarray + Affine with updated scale and intercept + """ + slicer = self.check_slicing(slicer, return_spatial=True) + + # Transform: + # sx 0 0 tx + # 0 sy 0 ty + # 0 0 sz tz + # 0 0 0 1 + transform = np.eye(4, dtype=int) + + for i, subslicer in enumerate(slicer): + if isinstance(subslicer, slice): + if subslicer.step == 0: + raise ValueError("slice step cannot be 0") + transform[i, i] = subslicer.step if subslicer.step is not None else 1 + transform[i, 3] = subslicer.start or 0 + # If slicer is None, nothing to do + + return self.img.affine.dot(transform) + + +class SpatialImage(DataobjImage): + ''' Template class for volumetric (3D/4D) images ''' + header_class = SpatialHeader + ImageSlicer = SpatialFirstSlicer def __init__(self, dataobj, affine, header=None, extra=None, file_map=None): @@ -477,69 +548,6 @@ def from_image(klass, img): klass.header_class.from_header(img.header), extra=img.extra.copy()) - def _check_slicing(self, slicer, return_spatial=False): - ''' Canonicalize slicers and check for scalar indices in spatial dims - - Parameters - ---------- - slicer : object - something that can be used to slice an array as in - ``arr[sliceobj]`` - return_spatial : bool - return only slices along spatial dimensions (x, y, z) - - Returns - ------- - slicer : object - Validated slicer object that will slice image's `dataobj` - without collapsing spatial dimensions - ''' - slicer = canonical_slicers(slicer, self.shape) - spatial_slices = slicer[:3] - for subslicer in spatial_slices: - if subslicer is None: - raise IndexError("New axis not permitted in spatial dimensions") - elif isinstance(subslicer, int): - raise IndexError("Scalar indices disallowed in spatial dimensions; " - "Use `[x]` or `x:x+1`.") - return spatial_slices if return_spatial else slicer - - def _slice_affine(self, slicer): - """ Retrieve affine for current image, if sliced by a given index - - Applies scaling if down-sampling is applied, and adjusts the intercept - to account for any cropping. - - Parameters - ---------- - slicer : object - something that can be used to slice an array as in - ``arr[sliceobj]`` - - Returns - ------- - affine : (4,4) ndarray - Affine with updated scale and intercept - """ - slicer = self._check_slicing(slicer, return_spatial=True) - - # Transform: - # sx 0 0 tx - # 0 sy 0 ty - # 0 0 sz tz - # 0 0 0 1 - transform = np.eye(4, dtype=int) - - for i, subslicer in enumerate(slicer): - if isinstance(subslicer, slice): - if subslicer.step == 0: - raise ValueError("slice step cannot be 0") - transform[i, i] = subslicer.step if subslicer.step is not None else 1 - transform[i, 3] = subslicer.start or 0 - # If slicer is None, nothing to do - - return self.affine.dot(transform) - @property def slicer(self): """ Slicer object that returns cropped and subsampled images @@ -558,11 +566,7 @@ def slicer(self): .. _aliasing: https://en.wikipedia.org/wiki/Aliasing """ - from .imageclasses import spatial_axes_first - if not spatial_axes_first(self): - raise ValueError("Cannot predict position of spatial axes for " - "Image type " + self.__class__.__name__) - return self.Slicer(self) + return self.ImageSlicer(self) def __getitem__(self, idx): From d4341af010527ea8eb0e8c02eecaafc53528a2ed Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Sat, 2 Jun 2018 13:53:34 -0400 Subject: [PATCH 34/36] ENH: Raise IndexError on empty slices --- nibabel/spatialimages.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index d14ec16f80..27fd7bf476 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -339,7 +339,11 @@ def __getitem__(self, slicer): slicer = self.check_slicing(slicer) except ValueError as err: raise IndexError(*err.args) + dataobj = self.img.dataobj[slicer] + if any(dim == 0 for dim in dataobj.shape): + raise IndexError("Empty slice requested") + affine = self.slice_affine(slicer) return self.img.__class__(dataobj.copy(), affine, self.img.header) From 44e86975eda69c470ae0ba47a6ae115b4a0d8355 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Sat, 2 Jun 2018 14:09:42 -0400 Subject: [PATCH 35/36] TEST: Update slice_affine test, check zero slices --- nibabel/tests/test_spatialimages.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index 8d8ad2fb81..b0f571023d 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -506,7 +506,11 @@ def test_slicer(self): with assert_raises(ValueError): img.slicer[:, ::0] with assert_raises(ValueError): - img._slice_affine((slice(None), slice(None, None, 0))) + img.slicer.slice_affine((slice(None), slice(None, None, 0))) + + # Don't permit zero-length slices + with assert_raises(IndexError): + img.slicer[:0] # No fancy indexing with assert_raises(IndexError): From 919b71a008234ecbd1a524d9b519f58e8cc1b607 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Sat, 2 Jun 2018 18:59:33 -0400 Subject: [PATCH 36/36] DOC: Some clarifying comments in SpatialFirstSlicer --- nibabel/spatialimages.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 27fd7bf476..9c3fef4fbb 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -328,6 +328,7 @@ class SpatialFirstSlicer(object): Checks that an image's first three axes are spatial ''' def __init__(self, img): + # Local import to avoid circular import on module load from .imageclasses import spatial_axes_first if not spatial_axes_first(img): raise ValueError("Cannot predict position of spatial axes for " @@ -365,6 +366,9 @@ def check_slicing(self, slicer, return_spatial=False): without collapsing spatial dimensions ''' slicer = canonical_slicers(slicer, self.img.shape) + # We can get away with this because we've checked the image's + # first three axes are spatial. + # More general slicers will need to be smarter, here. spatial_slices = slicer[:3] for subslicer in spatial_slices: if subslicer is None: