diff --git a/nibabel/analyze.py b/nibabel/analyze.py index a1c1cf1d2f..fb07443ebc 100644 --- a/nibabel/analyze.py +++ b/nibabel/analyze.py @@ -399,7 +399,7 @@ def from_header(klass, header=None, check=True): f"but output header {klass} does not support it") obj.set_data_dtype(header.get_data_dtype()) obj.set_data_shape(header.get_data_shape()) - obj.set_zooms(header.get_zooms()) + obj.set_zooms(header.get_zooms(units='raw'), units='raw') if check: obj.check_fix() return obj @@ -661,13 +661,22 @@ def get_base_affine(self): get_best_affine = get_base_affine - def get_zooms(self): - """ Get zooms from header + def get_zooms(self, units='norm', raise_unknown=False): + """ Get zooms (spacing between voxels along each axis) from header + + Parameters + ---------- + units : {'norm', 'raw'}, optional + Return zooms in normalized units of mm/sec for spatial/temporal or + as raw values stored in header. + raise_unkown : bool, optional + If normalized units are requested and the units are ambiguous, raise + a ``ValueError`` Returns ------- - z : tuple - tuple of header zoom values + zooms : tuple + tuple of header zoom values Examples -------- @@ -681,6 +690,8 @@ def get_zooms(self): >>> hdr.get_zooms() (3.0, 4.0) """ + if units not in ('norm', 'raw'): + raise ValueError("`units` parameter must be 'norm' or 'raw'") hdr = self._structarr dims = hdr['dim'] ndim = dims[0] @@ -689,11 +700,22 @@ def get_zooms(self): pixdims = hdr['pixdim'] return tuple(pixdims[1:ndim + 1]) - def set_zooms(self, zooms): + def set_zooms(self, zooms, units='norm'): """ Set zooms into header fields See docstring for ``get_zooms`` for examples + + Parameters + ---------- + zooms : sequence of floats + Zoom values to set in header + units : {'norm', 'raw'}, optional + Zooms are specified in normalized units of mm/sec for + spatial/temporal or as raw values to be interpreted according to + format specification. """ + if units not in ('norm', 'raw'): + raise ValueError("`units` parameter must be 'norm' or 'raw'") hdr = self._structarr dims = hdr['dim'] ndim = dims[0] diff --git a/nibabel/ecat.py b/nibabel/ecat.py index 54f600f147..afb85e78bd 100644 --- a/nibabel/ecat.py +++ b/nibabel/ecat.py @@ -578,8 +578,10 @@ def get_frame_affine(self, frame=0): z_off]) return aff - def get_zooms(self, frame=0): + def get_zooms(self, frame=0, units='norm', raise_unknown=False): """returns zooms ...pixdims""" + if units not in ('norm', 'raw'): + raise ValueError("`units` parameter must be 'norm' or 'raw'") subhdr = self.subheaders[frame] x_zoom = subhdr['x_pixel_size'] * 10 y_zoom = subhdr['y_pixel_size'] * 10 diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index 9d2cdb905b..6bc91db7a3 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -237,40 +237,75 @@ def _ndims(self): """ return 3 + (self._structarr['dims'][3] > 1) - def get_zooms(self): + def get_zooms(self, units='norm', raise_unknown=False): """ Get zooms from header Returns the spacing of voxels in the x, y, and z dimensions. For four-dimensional files, a fourth zoom is included, equal to the - repetition time (TR) in ms (see `The MGH/MGZ Volume Format - `_). + repetition time (TR). + TR is stored in milliseconds (see `The MGH/MGZ Volume Format + `_), + so if ``units == 'raw'``, the fourth zoom will be in ms. + If ``units == 'norm'`` (default), the fourth zoom will be in + seconds. To access only the spatial zooms, use `hdr['delta']`. + Parameters + ---------- + units : {'norm', 'raw'}, optional + Return zooms in normalized units of mm/sec for spatial/temporal or + as raw values stored in header. + raise_unkown : bool, optional + If normalized units are requested and the units are ambiguous, raise + a ``ValueError`` + Returns ------- - z : tuple - tuple of header zoom values + zooms : tuple + tuple of header zoom values .. _mghformat: https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/MghFormat#line-82 """ + if units == 'norm': + tfactor = 0.001 + elif units == 'raw': + tfactor = 1 + else: + raise ValueError("`units` parameter must be 'norm' or 'raw'") + # Do not return time zoom (TR) if 3D image - tzoom = (self['tr'],) if self._ndims() > 3 else () + tzoom = (self['tr'] * tfactor,) if self._ndims() > 3 else () return tuple(self._structarr['delta']) + tzoom - def set_zooms(self, zooms): + def set_zooms(self, zooms, units='norm'): """ Set zooms into header fields Sets the spacing of voxels in the x, y, and z dimensions. - For four-dimensional files, a temporal zoom (repetition time, or TR, in - ms) may be provided as a fourth sequence element. + For four-dimensional files, a temporal zoom (repetition time, or TR) + may be provided as a fourth sequence element. + + TR is stored in milliseconds (see `The MGH/MGZ Volume Format + `_), + so if ``units == 'raw'``, the fourth zoom will be interpreted as ms + and stored unmodified. + If ``units == 'norm'`` (default), the fourth zoom will be interpreted + as seconds, and converted to ms before storing in the header. Parameters ---------- zooms : sequence sequence of floats specifying spatial and (optionally) temporal zooms + units : {'norm', 'raw'}, optional + Zooms are specified in normalized units of mm/sec for + spatial/temporal dimensions or as raw values to be stored in + header. + + .. _mghformat: https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/MghFormat#line-82 """ + if units not in ('norm', 'raw'): + raise ValueError("`units` parameter must be 'norm' or 'raw'") hdr = self._structarr zooms = np.asarray(zooms) ndims = self._ndims() @@ -283,7 +318,13 @@ def set_zooms(self, zooms): if len(zooms) == 4: if zooms[3] < 0: raise HeaderDataError(f'TR must be non-negative; got {zooms[3]}') - hdr['tr'] = zooms[3] + if units == 'norm': + tfactor = 1000 + elif units == 'raw': + tfactor = 1 + else: + raise ValueError("`units` parameter must be 'norm' or 'raw'") + hdr['tr'] = zooms[3] * tfactor def get_data_shape(self): """ Get shape of data diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index 4c812087c2..5e7c61a887 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -74,7 +74,8 @@ def test_read_mgh(): assert_almost_equal(h['flip_angle'], 0.0) assert_almost_equal(h['te'], 0.0) assert_almost_equal(h['ti'], 0.0) - assert_array_almost_equal(h.get_zooms(), [1, 1, 1, 2]) + assert_array_almost_equal(h.get_zooms(units='raw'), [1, 1, 1, 2]) + assert_array_almost_equal(h.get_zooms(units='norm'), [1, 1, 1, 0.002]) assert_array_almost_equal(h.get_vox2ras(), v2r) assert_array_almost_equal(h.get_vox2ras_tkr(), v2rtkr) @@ -147,9 +148,14 @@ def test_write_noaffine_mgh(): def test_set_zooms(): mgz = load(MGZ_FNAME) h = mgz.header - assert_array_almost_equal(h.get_zooms(), [1, 1, 1, 2]) - h.set_zooms([1, 1, 1, 3]) - assert_array_almost_equal(h.get_zooms(), [1, 1, 1, 3]) + assert_array_almost_equal(h.get_zooms(units='raw'), [1, 1, 1, 2]) + assert_array_almost_equal(h.get_zooms(units='norm'), [1, 1, 1, 0.002]) + h.set_zooms([1, 1, 1, 3], units='raw') + assert_array_almost_equal(h.get_zooms(units='raw'), [1, 1, 1, 3]) + assert_array_almost_equal(h.get_zooms(units='norm'), [1, 1, 1, 0.003]) + h.set_zooms([1, 1, 1, 3], units='norm') + assert_array_almost_equal(h.get_zooms(units='raw'), [1, 1, 1, 3000]) + assert_array_almost_equal(h.get_zooms(units='norm'), [1, 1, 1, 3]) for zooms in ((-1, 1, 1, 1), (1, -1, 1, 1), (1, 1, -1, 1), @@ -351,6 +357,63 @@ def check_dtypes(self, expected, actual): # MGH requires the actual to be a big endian version of expected assert expected.newbyteorder('>') == actual + def test_zooms_edge_cases(self): + img_klass = self.image_class + aff = np.eye(4) + arr = np.arange(120, dtype=np.int16).reshape((2, 3, 4, 5)) + img = img_klass(arr, aff) + + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (1, 1, 1, 0)) + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (1, 1, 1, 0)) + + img.header.set_zooms((1, 1, 1, 2000), units='raw') + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (1, 1, 1, 2000)) + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (1, 1, 1, 2)) + assert_array_almost_equal(img.header.get_zooms(), (1, 1, 1, 2)) + + img.header.set_zooms((2, 2, 2, 3), units='norm') + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (2, 2, 2, 3000)) + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (2, 2, 2, 3)) + assert_array_almost_equal(img.header.get_zooms(), (2, 2, 2, 3)) + + # It's legal to set zooms for spatial dimensions only + img.header.set_zooms((3, 3, 3), units='norm') + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (3, 3, 3, 3000)) + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (3, 3, 3, 3)) + assert_array_almost_equal(img.header.get_zooms(), (3, 3, 3, 3)) + + arr = np.arange(24, dtype=np.int16).reshape((2, 3, 4)) + img = img_klass(arr, aff) + + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (1, 1, 1)) + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (1, 1, 1)) + + img.header.set_zooms((2, 2, 2), units='raw') + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (2, 2, 2)) + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (2, 2, 2)) + + img.header.set_zooms((3, 3, 3), units='norm') + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (3, 3, 3)) + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (3, 3, 3)) + + # Cannot set TR as zoom for 3D image + assert_raises(HeaderDataError, img.header.set_zooms, (4, 4, 4, 5), 'raw') + assert_raises(HeaderDataError, img.header.set_zooms, (4, 4, 4, 5), 'norm') + class TestMGHHeader(tws._TestLabeledWrapStruct): header_class = MGHHeader diff --git a/nibabel/minc1.py b/nibabel/minc1.py index c0ae95bd7b..1525ddd004 100644 --- a/nibabel/minc1.py +++ b/nibabel/minc1.py @@ -90,8 +90,10 @@ def get_data_dtype(self): def get_data_shape(self): return self._image.data.shape - def get_zooms(self): + def get_zooms(self, units='norm', raise_unknown=False): """ Get real-world sizes of voxels """ + if units not in ('norm', 'raw'): + raise ValueError("`units` parameter must be 'norm' or 'raw'") # zooms must be positive; but steps in MINC can be negative return tuple([abs(float(dim.step)) if hasattr(dim, 'step') else 1.0 for dim in self._dims]) diff --git a/nibabel/nifti1.py b/nibabel/nifti1.py index 1bffac10ce..4efa3bbc28 100644 --- a/nibabel/nifti1.py +++ b/nibabel/nifti1.py @@ -784,7 +784,7 @@ def get_data_shape(self): Expanding number of dimensions gets default zooms - >>> hdr.get_zooms() + >>> hdr.get_zooms(units='norm') (1.0, 1.0, 1.0) Notes @@ -1671,12 +1671,159 @@ def set_xyzt_units(self, xyz=None, t=None): xyz = 0 if t is None: t = 0 - xyz_code = self.structarr['xyzt_units'] % 8 - t_code = self.structarr['xyzt_units'] - xyz_code xyz_code = unit_codes[xyz] t_code = unit_codes[t] self.structarr['xyzt_units'] = xyz_code + t_code + def get_zooms(self, units=None, raise_unknown=False): + ''' Get zooms (spacing between voxels along each axis) from header + + NIfTI-1 headers may specify that zooms are encoded in units other than + mm and sec (see ``get_xyzt_units``). + Default behavior has been to return the raw zooms, and leave it to the + programmer to handle non-standard units. + However, most files indicate mm/sec units or have unspecified units, + and it is common practice to neglect specified units and assume all + files will be in mm/sec. + + The default behavior for ``get_zooms`` will remain to return the raw + zooms until version 4.0, when it will change to return zooms in + normalized mm/sec units. + Because the default behavior will change, a warning will be given to + prompt programmers to specify whether they intend to retrieve raw + values, or values coerced into normalized units. + + Parameters + ---------- + units : {'norm', 'raw'} + Return zooms in normalized units of mm/sec for spatial/temporal or + as raw values stored in header. + raise_unkown : bool, optional + If normalized units are requested and the units are ambiguous, raise + a ``ValueError`` + + Returns + ------- + zooms : tuple + tuple of header zoom values + + ''' + if units is None: + units = 'raw' + warnings.warn('Units not specified in `{}.get_zooms`. Returning ' + 'raw zooms, but default will change to normalized.\n' + 'Please explicitly specify units parameter.' + ''.format(self.__class__.__name__), + FutureWarning, stacklevel=2) + + raw_zooms = super(Nifti1Header, self).get_zooms(units='raw', + raise_unknown=False) + + if units == 'raw': + return raw_zooms + + elif units != 'norm': + raise ValueError("`units` parameter must be 'norm' or 'raw'") + + xyz_zooms = raw_zooms[:3] + t_zoom = raw_zooms[3:4] # Tuple of length 0 or 1 + + xyz_code, t_code = self.get_xyzt_units() + xyz_msg = t_msg = '' + if xyz_code == 'unknown': + xyz_msg = 'Unknown spatial units' + xyz_code = 'mm' + if t_zoom: + if t_code == 'unknown': + t_msg = 'Unknown time units' + t_code = 'sec' + elif t_code in ('hz', 'ppm', 'rads'): + t_msg = 'Unconvertible temporal units: {}'.format(t_code) + + if raise_unknown and (xyz_msg or t_msg.startswith('Unknown')): + if xyz_msg and t_msg.startswith('Unknown'): + msg = 'Unknown spatial and time units' + else: + msg = xyz_msg or t_msg + raise ValueError("Error: {}".format(msg)) + if xyz_msg: + warnings.warn('{} - assuming mm'.format(xyz_msg)) + if t_msg: + if t_code == 'sec': + warnings.warn('{} - assuming sec'.format(t_msg)) + else: + warnings.warn(t_msg) + + xyz_factor = {'meter': 1000, 'mm': 1, 'micron': 0.001}[xyz_code] + xyz_zooms = tuple(np.array(xyz_zooms) * xyz_factor) + + if t_zoom: + t_factor = {'sec': 1, 'msec': 0.001, 'usec': 0.000001, + 'hz': 1, 'ppm': 1, 'rads': 1}[t_code] + t_zoom = (t_zoom[0] * t_factor,) + + return xyz_zooms + t_zoom + + def set_zooms(self, zooms, units=None): + ''' Set zooms into header fields + + NIfTI-1 headers may specify that zooms are encoded in units other than + mm and sec. + Default behavior has been to set the raw zooms, and leave unit handling + to a separate method (``set_xyzt_units``). + However, most files indicate mm/sec units or have unspecified units, + and it is common practice to neglect specified units and assume all + files will be in mm/sec. + + The default behavior for ``set_zooms`` will remain to update only the + zooms until version 4.0, when it will change to setting units in + normalized mm/sec units. + Because the default behavior will change, a warning will be given to + prompt programmers to specify whether they intend to set only the raw + values, or set normalized units, as well. + + For setting normalized units, the following rules apply: + Spatial units will be set to mm. + The temporal code for a <4D image is considered "unknown". + If the current temporal units are Hz, PPM or Rads, these are left + unchanged. + Otherwise, the fourth zoom is considered to be in seconds. + + Parameters + ---------- + zooms : sequence of floats + Zoom values to set in header + units : {'norm', 'raw', tuple of unit codes}, optional + Zooms are specified in normalized units of mm/sec for + spatial/temporal, as raw values to be interpreted according to + format specification, or according to a tuple of NIFTI unit + codes. + + ''' + if units is None: + units = 'raw' + warnings.warn('Units not specified in `{}.set_zooms`. Setting ' + 'raw zooms, but default will change to normalized.\n' + 'Please explicitly specify units parameter.' + ''.format(self.__class__.__name__), + FutureWarning, stacklevel=2) + + if units not in ('norm', 'raw') and not isinstance(units, tuple): + raise ValueError("`units` parameter must be 'norm', 'raw'," + " or a tuple of unit codes (see set_xyzt_units)") + super(Nifti1Header, self).set_zooms(zooms, units='raw') + + if isinstance(units, tuple): + self.set_xyzt_units(*units) + elif units == 'norm': + _, t_code = self.get_xyzt_units() + xyz_code = 'mm' + if len(zooms) == 3: + t_code = 'unknown' + elif len(zooms) > 3 and t_code in ('unknown', 'sec', 'msec', 'usec'): + t_code = 'sec' + self.set_xyzt_units(xyz_code, t_code) + def _clean_after_mapping(self): """ Set format-specific stuff after converting header from mapping diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 09744d0149..93877e947c 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -224,17 +224,48 @@ def set_data_shape(self, shape): nzs = min(len(self._zooms), ndim) self._zooms = self._zooms[:nzs] + (1.0,) * (ndim - nzs) - def get_zooms(self): + def get_zooms(self, units='norm', raise_unknown=False): + ''' Get zooms (spacing between voxels along each axis) from header + + Parameters + ---------- + units : {'norm', 'raw'}, optional + Return zooms in normalized units of mm/sec for spatial/temporal or + as raw values stored in header. + raise_unknown : bool, optional + If canonical units are requested and the units are ambiguous, raise + a ``ValueError`` + + Returns + ------- + zooms : tuple + Spacing between voxels along each axis + ''' + if units not in ('norm', 'raw'): + raise ValueError("`units` parameter must be 'norm' or 'raw'") return self._zooms - def set_zooms(self, zooms): - zooms = tuple([float(z) for z in zooms]) + def set_zooms(self, zooms, units='norm'): + ''' Set zooms (spacing between voxels along each axis) in header + + Parameters + ---------- + zooms : sequence of floats + Spacing between voxels along each axis + units : {'norm', 'raw'}, optional + Zooms are specified in normalized units of mm/sec for + spatial/temporal or as raw values to be interpreted according to + format specification. + ''' + if units not in ('norm', 'raw'): + raise ValueError("`units` parameter must be 'norm' or 'raw'") + zooms = tuple(float(z) for z in zooms) shape = self.get_data_shape() ndim = len(shape) if len(zooms) != ndim: raise HeaderDataError('Expecting %d zoom values for ndim %d' % (ndim, ndim)) - if len([z for z in zooms if z < 0]): + if any(z < 0 for z in zooms): raise HeaderDataError('zooms must be positive') self._zooms = zooms diff --git a/nibabel/tests/test_analyze.py b/nibabel/tests/test_analyze.py index d91769bc73..de02f3f3aa 100644 --- a/nibabel/tests/test_analyze.py +++ b/nibabel/tests/test_analyze.py @@ -91,7 +91,8 @@ def test_general_init(self): assert_array_equal(np.diag(hdr.get_base_affine()), [-1, 1, 1, 1]) # But zooms only go with number of dimensions - assert hdr.get_zooms() == (1.0,) + assert hdr.get_zooms(units='raw') == (1.0,) + assert hdr.get_zooms(units='norm') == (1.0,) def test_header_size(self): assert self.header_class.template_dtype.itemsize == self.sizeof_hdr @@ -437,7 +438,8 @@ def test_data_shape_zooms_affine(self): else: assert hdr.get_data_shape() == (0,) # Default zoom - for 3D - is 1(()) - assert hdr.get_zooms() == (1,) * L + assert hdr.get_zooms(units='raw') == (1,) * L + assert hdr.get_zooms(units='norm') == (1,) * L # errors if zooms do not match shape if len(shape): with pytest.raises(HeaderDataError): @@ -454,16 +456,19 @@ def test_data_shape_zooms_affine(self): # it again reverts the previously set zoom values to 1.0 hdr = self.header_class() hdr.set_data_shape((1, 2, 3)) - hdr.set_zooms((4, 5, 6)) - assert_array_equal(hdr.get_zooms(), (4, 5, 6)) + hdr.set_zooms((4, 5, 6), units='norm') + assert_array_equal(hdr.get_zooms(units='raw'), (4, 5, 6)) + assert_array_equal(hdr.get_zooms(units='norm'), (4, 5, 6)) hdr.set_data_shape((1, 2)) - assert_array_equal(hdr.get_zooms(), (4, 5)) + assert_array_equal(hdr.get_zooms(units='raw'), (4, 5)) + assert_array_equal(hdr.get_zooms(units='norm'), (4, 5)) hdr.set_data_shape((1, 2, 3)) - assert_array_equal(hdr.get_zooms(), (4, 5, 1)) + assert_array_equal(hdr.get_zooms(units='raw'), (4, 5, 1)) + assert_array_equal(hdr.get_zooms(units='norm'), (4, 5, 1)) # Setting zooms changes affine assert_array_equal(np.diag(hdr.get_base_affine()), [-4, 5, 1, 1]) - hdr.set_zooms((1, 1, 1)) + hdr.set_zooms((1, 1, 1), units='norm') assert_array_equal(np.diag(hdr.get_base_affine()), [-1, 1, 1, 1]) @@ -471,7 +476,7 @@ def test_default_x_flip(self): hdr = self.header_class() hdr.default_x_flip = True hdr.set_data_shape((1, 2, 3)) - hdr.set_zooms((1, 1, 1)) + hdr.set_zooms((1, 1, 1), units='norm') assert_array_equal(np.diag(hdr.get_base_affine()), [-1, 1, 1, 1]) hdr.default_x_flip = False @@ -490,7 +495,7 @@ def test_orientation(self): hdr = self.header_class() assert hdr.default_x_flip hdr.set_data_shape((3, 5, 7)) - hdr.set_zooms((4, 5, 6)) + hdr.set_zooms((4, 5, 6), units='norm') aff = np.diag((-4, 5, 6, 1)) aff[:3, 3] = np.array([1, 2, 3]) * np.array([-4, 5, 6]) * -1 assert_array_equal(hdr.get_base_affine(), aff) @@ -517,7 +522,7 @@ def test_from_header(self): hdr = klass() hdr.set_data_dtype(np.float64) hdr.set_data_shape((1, 2, 3)) - hdr.set_zooms((3.0, 2.0, 1.0)) + hdr.set_zooms((3.0, 2.0, 1.0), units='norm') for check in (True, False): copy = klass.from_header(hdr, check=check) assert hdr == copy @@ -529,18 +534,19 @@ def get_data_dtype(self): return np.dtype('i2') def get_data_shape(self): return (5, 4, 3) - def get_zooms(self): return (10.0, 9.0, 8.0) + def get_zooms(self, units=None, raise_unknown=None): return (10.0, 9.0, 8.0) converted = klass.from_header(C()) assert isinstance(converted, klass) assert converted.get_data_dtype() == np.dtype('i2') assert converted.get_data_shape() == (5, 4, 3) - assert converted.get_zooms() == (10.0, 9.0, 8.0) + assert converted.get_zooms(units='raw') == (10.0, 9.0, 8.0) + assert converted.get_zooms(units='norm') == (10.0, 9.0, 8.0) def test_base_affine(self): klass = self.header_class hdr = klass() hdr.set_data_shape((3, 5, 7)) - hdr.set_zooms((3, 2, 1)) + hdr.set_zooms((3, 2, 1), units='norm') assert hdr.default_x_flip assert_array_almost_equal( hdr.get_base_affine(), @@ -640,12 +646,12 @@ def get_data_shape(self): class H4(H3): - def get_zooms(self): + def get_zooms(self, units=None, raise_unknown=None): return 4., 5., 6. exp_hdr = klass() exp_hdr.set_data_dtype(np.dtype('u1')) exp_hdr.set_data_shape((2, 3, 4)) - exp_hdr.set_zooms((4, 5, 6)) + exp_hdr.set_zooms((4, 5, 6), units='raw') assert klass.from_header(H4()) == exp_hdr # cal_max, cal_min get properly set from ``as_analyze_map`` @@ -682,7 +688,7 @@ def as_analyze_map(self): def test_best_affine(): hdr = AnalyzeHeader() hdr.set_data_shape((3, 5, 7)) - hdr.set_zooms((4, 5, 6)) + hdr.set_zooms((4, 5, 6), units='norm') assert_array_equal(hdr.get_base_affine(), hdr.get_best_affine()) @@ -830,24 +836,24 @@ def test_header_updating(self): # With a None affine - don't overwrite zooms img = img_klass(np.zeros((2, 3, 4)), None) hdr = img.header - hdr.set_zooms((4, 5, 6)) + hdr.set_zooms((4, 5, 6), units='norm') # Save / reload using bytes IO objects for key, value in img.file_map.items(): value.fileobj = BytesIO() img.to_file_map() hdr_back = img.from_file_map(img.file_map).header - assert_array_equal(hdr_back.get_zooms(), (4, 5, 6)) + assert_array_equal(hdr_back.get_zooms(units='norm'), (4, 5, 6)) # With a real affine, update zooms img = img_klass(np.zeros((2, 3, 4)), np.diag([2, 3, 4, 1]), hdr) hdr = img.header - assert_array_equal(hdr.get_zooms(), (2, 3, 4)) + assert_array_equal(hdr.get_zooms(units='norm'), (2, 3, 4)) # Modify affine in-place? Update on save. img.affine[0, 0] = 9 for key, value in img.file_map.items(): value.fileobj = BytesIO() img.to_file_map() hdr_back = img.from_file_map(img.file_map).header - assert_array_equal(hdr.get_zooms(), (9, 3, 4)) + assert_array_equal(hdr.get_zooms(units='norm'), (9, 3, 4)) # Modify data in-place? Update on save data = img.get_fdata() data.shape = (3, 2, 4) diff --git a/nibabel/tests/test_ecat.py b/nibabel/tests/test_ecat.py index 9e56fd73c7..d15d098ecd 100644 --- a/nibabel/tests/test_ecat.py +++ b/nibabel/tests/test_ecat.py @@ -169,6 +169,10 @@ def test_subheader(self): np.array([2.20241979, 2.20241979, 3.125, 1.])) assert self.subhdr.get_zooms()[0] == 2.20241978764534 assert self.subhdr.get_zooms()[2] == 3.125 + assert_array_equal(self.subhdr.get_zooms(units='raw'), + self.subhdr.get_zooms(units='norm')) + with pytest.raises(ValueError): + self.subhdr.get_zooms(units='badarg') assert self.subhdr._get_data_dtype(0) == np.int16 #assert_equal(self.subhdr._get_frame_offset(), 1024) assert self.subhdr._get_frame_offset() == 1536 diff --git a/nibabel/tests/test_minc1.py b/nibabel/tests/test_minc1.py index 4fecf5782e..1bbe7c8746 100644 --- a/nibabel/tests/test_minc1.py +++ b/nibabel/tests/test_minc1.py @@ -117,6 +117,10 @@ def test_mincfile(self): assert mnc.get_data_dtype().type == tp['dtype'] assert mnc.get_data_shape() == tp['shape'] assert mnc.get_zooms() == tp['zooms'] + assert mnc.get_zooms(units='raw') == tp['zooms'] + assert mnc.get_zooms(units='norm') == tp['zooms'] + with pytest.raises(ValueError): + mnc.get_zooms(units='badarg') assert_array_equal(mnc.get_affine(), tp['affine']) data = mnc.get_scaled_data() assert data.shape == tp['shape'] diff --git a/nibabel/tests/test_nifti1.py b/nibabel/tests/test_nifti1.py index 63cf13c103..9993c52b75 100644 --- a/nibabel/tests/test_nifti1.py +++ b/nibabel/tests/test_nifti1.py @@ -718,6 +718,30 @@ def test_recoded_fields(self): hdr['slice_code'] = 4 # alternating decreasing assert hdr.get_value_label('slice_code') == 'alternating decreasing' + def test_general_init(self): + with clear_and_catch_warnings(modules=(nifti1,)) as warns: + warnings.filterwarnings('ignore', 'get_zooms', FutureWarning) + warnings.filterwarnings('ignore', 'set_zooms', FutureWarning) + warnings.filterwarnings('ignore', 'Unknown (spatial|time) units', + UserWarning) + super(TestNifti1PairHeader, self).test_general_init() + + def test_from_header(self): + with clear_and_catch_warnings(modules=(nifti1,)) as warns: + warnings.filterwarnings('ignore', 'get_zooms', FutureWarning) + warnings.filterwarnings('ignore', 'set_zooms', FutureWarning) + warnings.filterwarnings('ignore', 'Unknown (spatial|time) units', + UserWarning) + super(TestNifti1PairHeader, self).test_from_header() + + def test_data_shape_zooms_affine(self): + with clear_and_catch_warnings(modules=(nifti1,)) as warns: + warnings.filterwarnings('ignore', 'get_zooms', FutureWarning) + warnings.filterwarnings('ignore', 'set_zooms', FutureWarning) + warnings.filterwarnings('ignore', 'Unknown (spatial|time) units', + UserWarning) + super(TestNifti1PairHeader, self).test_data_shape_zooms_affine() + def unshear_44(affine): RZS = affine[:3, :3] @@ -881,11 +905,11 @@ def test_set_qform(self): img.set_qform(new_affine, 1, update_affine=False) assert_array_almost_equal(img.affine, aff_affine) # Clear qform using None, zooms unchanged - assert_array_almost_equal(hdr.get_zooms(), [1.1, 1.1, 1.1]) + assert_array_almost_equal(hdr.get_zooms(units='raw'), [1.1, 1.1, 1.1]) img.set_qform(None) qaff, code = img.get_qform(coded=True) assert (qaff, code) == (None, 0) - assert_array_almost_equal(hdr.get_zooms(), [1.1, 1.1, 1.1]) + assert_array_almost_equal(hdr.get_zooms(units='raw'), [1.1, 1.1, 1.1]) # Best affine similarly assert_array_almost_equal(img.affine, hdr.get_best_affine()) # If sform is not set, qform should update affine @@ -942,9 +966,9 @@ def test_set_sform(self): assert_array_almost_equal(img.affine, aff_affine) # zooms do not get updated when qform is 0 assert_array_almost_equal(img.get_qform(), orig_aff) - assert_array_almost_equal(hdr.get_zooms(), [2.2, 3.3, 4.3]) + assert_array_almost_equal(hdr.get_zooms(units='raw'), [2.2, 3.3, 4.3]) img.set_qform(None) - assert_array_almost_equal(hdr.get_zooms(), [2.2, 3.3, 4.3]) + assert_array_almost_equal(hdr.get_zooms(units='raw'), [2.2, 3.3, 4.3]) # Set sform using new_affine when qform is set img.set_qform(qform_affine, 1) img.set_sform(new_affine, 1) @@ -953,7 +977,7 @@ def test_set_sform(self): assert_array_almost_equal(saff, new_affine) assert_array_almost_equal(img.affine, new_affine) # zooms follow qform - assert_array_almost_equal(hdr.get_zooms(), [1.2, 1.2, 1.2]) + assert_array_almost_equal(hdr.get_zooms(units='raw'), [1.2, 1.2, 1.2]) # Clear sform using None, best_affine should fall back on qform img.set_sform(None) assert hdr['sform_code'] == 0 @@ -1042,7 +1066,7 @@ def test_load_pixdims(self): # Check qform, sform, pixdims are the same assert_array_equal(img_hdr.get_qform(), qaff) assert_array_equal(img_hdr.get_sform(), saff) - assert_array_equal(img_hdr.get_zooms(), [2, 3, 4]) + assert_array_equal(img_hdr.get_zooms(units='raw'), [2, 3, 4]) # Save to stringio re_simg = bytesio_round_trip(simg) assert_array_equal(re_simg.get_fdata(), arr) @@ -1050,7 +1074,7 @@ def test_load_pixdims(self): rimg_hdr = re_simg.header assert_array_equal(rimg_hdr.get_qform(), qaff) assert_array_equal(rimg_hdr.get_sform(), saff) - assert_array_equal(rimg_hdr.get_zooms(), [2, 3, 4]) + assert_array_equal(rimg_hdr.get_zooms(units='raw'), [2, 3, 4]) def test_affines_init(self): # Test we are doing vaguely spec-related qform things. The 'spec' here @@ -1064,20 +1088,20 @@ def test_affines_init(self): hdr = img.header assert hdr['qform_code'] == 0 assert hdr['sform_code'] == 2 - assert_array_equal(hdr.get_zooms(), [2, 3, 4]) + assert_array_equal(hdr.get_zooms(units='raw'), [2, 3, 4]) # This is also true for affines with header passed qaff = np.diag([3, 4, 5, 1]) saff = np.diag([6, 7, 8, 1]) hdr.set_qform(qaff, code='scanner') hdr.set_sform(saff, code='talairach') - assert_array_equal(hdr.get_zooms(), [3, 4, 5]) + assert_array_equal(hdr.get_zooms(units='raw'), [3, 4, 5]) img = IC(arr, aff, hdr) new_hdr = img.header # Again affine is sort of anonymous space assert new_hdr['qform_code'] == 0 assert new_hdr['sform_code'] == 2 assert_array_equal(new_hdr.get_sform(), aff) - assert_array_equal(new_hdr.get_zooms(), [2, 3, 4]) + assert_array_equal(new_hdr.get_zooms(units='raw'), [2, 3, 4]) # But if no affine passed, codes and matrices stay the same img = IC(arr, None, hdr) new_hdr = img.header @@ -1086,7 +1110,7 @@ def test_affines_init(self): assert new_hdr['sform_code'] == 3 # Still talairach assert_array_equal(new_hdr.get_sform(), saff) # Pixdims as in the original header - assert_array_equal(new_hdr.get_zooms(), [3, 4, 5]) + assert_array_equal(new_hdr.get_zooms(units='raw'), [3, 4, 5]) def test_read_no_extensions(self): IC = self.image_class @@ -1176,6 +1200,123 @@ def test_static_dtype_aliases(self): img_rt = bytesio_round_trip(img) assert img_rt.get_data_dtype() == effective_dt + def test_zooms_edge_cases(self): + img_klass = self.image_class + arr = np.arange(120, dtype=np.int16).reshape((2, 3, 4, 5)) + aff = np.eye(4) + img = img_klass(arr, aff) + + # Unknown units = 2 warnings + with clear_and_catch_warnings(modules=(nifti1,)) as warns: + warnings.simplefilter('always') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (1, 1, 1, 1)) + assert_equal(len(warns), 2) + assert_raises(ValueError, img.header.get_zooms, + units='norm', raise_unknown=True) + + img.header.set_xyzt_units(xyz='meter') + with clear_and_catch_warnings(modules=(nifti1,)) as warns: + warnings.simplefilter('always') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (1000, 1000, 1000, 1)) + assert_equal(len(warns), 1) + assert_raises(ValueError, img.header.get_zooms, + units='norm', raise_unknown=True) + + img.header.set_xyzt_units(xyz='mm', t='sec') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (1, 1, 1, 1)) + img.header.set_xyzt_units(xyz='micron', t='sec') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (0.001, 0.001, 0.001, 1)) + + img.header.set_xyzt_units(t='sec') + with clear_and_catch_warnings(modules=(nifti1,)) as warns: + warnings.simplefilter('always') + assert_array_equal(img.header.get_zooms(units='norm'), + (1, 1, 1, 1)) + assert_equal(len(warns), 1) + assert_raises(ValueError, img.header.get_zooms, + units='norm', raise_unknown=True) + + img.header.set_xyzt_units(xyz='mm', t='msec') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (1, 1, 1, 0.001)) + + img.header.set_xyzt_units(xyz='mm', t='usec') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (1, 1, 1, 0.000001)) + + img.header.set_xyzt_units(xyz='meter', t='usec') + assert_equal(img.header.get_xyzt_units(), ('meter', 'usec')) + + # Verify `set_zooms(units='raw')` leaves units unchanged + img.header.set_zooms((2, 2, 2, 2.5), units='raw') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (2000, 2000, 2000, 0.0000025)) + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (2, 2, 2, 2.5)) + assert_equal(img.header.get_xyzt_units(), ('meter', 'usec')) + + # Verify `set_zooms(units=)` sets units explicitly + img.header.set_zooms((2, 2, 2, 2.5), units=('micron', 'msec')) + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (0.002, 0.002, 0.002, 0.0025)) + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (2, 2, 2, 2.5)) + assert_equal(img.header.get_xyzt_units(), ('micron', 'msec')) + + # Verify `set_zooms(units='norm')` resets units + img.header.set_zooms((2, 2, 2, 2.5), units='norm') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (2, 2, 2, 2.5)) + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (2, 2, 2, 2.5)) + assert_equal(img.header.get_xyzt_units(), ('mm', 'sec')) + + # Non-temporal t units are not transformed + img.header.set_zooms((1, 1, 1, 1.5), units=('mm', 'ppm')) + with clear_and_catch_warnings(modules=(nifti1,)) as warns: + warnings.simplefilter('always') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (1, 1, 1, 1.5)) + assert_equal(len(warns), 1) + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (1, 1, 1, 1.5)) + + # Non-temporal t units are not normalized + img.header.set_zooms((2, 2, 2, 3.5), units='norm') + with clear_and_catch_warnings(modules=(nifti1,)) as warns: + warnings.simplefilter('always') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (2, 2, 2, 3.5)) + assert_equal(len(warns), 1) + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (2, 2, 2, 3.5)) + assert_equal(img.header.get_xyzt_units(), ('mm', 'ppm')) + + # Unknown t units are normalized to seconds + img.header.set_xyzt_units(xyz='mm', t='unknown') + img.header.set_zooms((2, 2, 2, 3.5), units='norm') + assert_array_almost_equal(img.header.get_zooms(units='norm'), + (2, 2, 2, 3.5)) + assert_array_almost_equal(img.header.get_zooms(units='raw'), + (2, 2, 2, 3.5)) + assert_equal(img.header.get_xyzt_units(), ('mm', 'sec')) + + assert_raises(ValueError, img.header.get_zooms, units='badparam') + assert_raises(ValueError, img.header.set_zooms, (3, 3, 3, 3.5), + units='badparam') + + def test_no_finite_values(self): + with clear_and_catch_warnings(modules=(nifti1,)) as warns: + warnings.filterwarnings('ignore', 'get_zooms', FutureWarning) + warnings.filterwarnings('ignore', 'set_zooms', FutureWarning) + warnings.filterwarnings('ignore', 'Unknown (spatial|time) units', + UserWarning) + super(TestNifti1Pair, self).test_no_finite_values() + class TestNifti1Image(TestNifti1Pair): # Run analyze-flavor spatialimage tests diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index fc11452151..246f92b712 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -214,7 +214,7 @@ def test_isolation(self): # Pass it back in img = img_klass(arr, aff, ihdr) # Check modifying header outside does not modify image - ihdr.set_zooms((4, 5, 6)) + ihdr.set_zooms((4, 5, 6), units='norm') assert img.header != ihdr def test_float_affine(self): @@ -528,6 +528,34 @@ def test_slicer(self): assert (sliced_data == img.get_data()[sliceobj]).all() assert (sliced_data == img.get_fdata()[sliceobj]).all() + def test_zooms(self): + ''' Should be true for all images ''' + img_klass = self.image_class + arr = np.arange(120, dtype=np.int16).reshape((2, 3, 4, 5)) + aff = np.eye(4) + img = img_klass(arr, aff) + img.header.set_zooms((2, 2, 2, 2.5), units='norm') + assert_array_equal(img.header.get_zooms(units='norm'), (2, 2, 2, 2.5)) + with assert_raises(ValueError): + img.header.set_zooms((1, 1, 1, 1), units='badarg') + with assert_raises(ValueError): + img.header.get_zooms(units='badarg') + with assert_raises(HeaderDataError): + img.header.set_zooms((-1, 1, 1, 1), units='norm') + + def test_zooms_edge_cases(self): + ''' Override for classes where *_norm_zooms != *_zooms ''' + img_klass = self.image_class + arr = np.arange(120, dtype=np.int16).reshape((2, 3, 4, 5)) + aff = np.eye(4) + img = img_klass(arr, aff) + img.header.set_zooms((2, 2, 2, 2.5), units='raw') + assert_array_equal(img.header.get_zooms(units='raw'), (2, 2, 2, 2.5)) + assert_array_equal(img.header.get_zooms(units='norm'), (2, 2, 2, 2.5)) + img.header.set_zooms((2, 2, 2, 2.5), units='norm') + assert_array_equal(img.header.get_zooms(units='raw'), (2, 2, 2, 2.5)) + assert_array_equal(img.header.get_zooms(units='norm'), (2, 2, 2, 2.5)) + class MmapImageMixin: """ Mixin for testing images that may return memory maps """