diff --git a/nibabel/freesurfer/mghformat.py b/nibabel/freesurfer/mghformat.py index 1bc9071538..c35acc7cf5 100644 --- a/nibabel/freesurfer/mghformat.py +++ b/nibabel/freesurfer/mghformat.py @@ -13,30 +13,39 @@ from os.path import splitext import numpy as np -from ..volumeutils import (array_to_file, array_from_file, Recoder) +from ..affines import voxel_sizes, from_matvec +from ..volumeutils import (array_to_file, array_from_file, endian_codes, + Recoder) from ..spatialimages import HeaderDataError, SpatialImage -from ..fileholders import FileHolder, copy_file_map -from ..arrayproxy import ArrayProxy +from ..fileholders import FileHolder +from ..arrayproxy import ArrayProxy, reshape_dataobj from ..keywordonly import kw_only_meth from ..openers import ImageOpener +from ..batteryrunners import BatteryRunner, Report +from ..wrapstruct import LabeledWrapStruct +from ..deprecated import deprecate_with_version # mgh header # See https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/MghFormat DATA_OFFSET = 284 # Note that mgh data is strictly big endian ( hence the > sign ) header_dtd = [ - ('version', '>i4'), - ('dims', '>i4', (4,)), - ('type', '>i4'), - ('dof', '>i4'), - ('goodRASFlag', '>i2'), - ('delta', '>f4', (3,)), - ('Mdc', '>f4', (3, 3)), - ('Pxyz_c', '>f4', (3,)) + ('version', '>i4'), # 0; must be 1 + ('dims', '>i4', (4,)), # 4; width, height, depth, nframes + ('type', '>i4'), # 20; data type + ('dof', '>i4'), # 24; degrees of freedom + ('goodRASFlag', '>i2'), # 28; Mdc, Pxyz_c fields valid + ('delta', '>f4', (3,)), # 30; zooms (X, Y, Z) + ('Mdc', '>f4', (3, 3)), # 42; TRANSPOSE of direction cosine matrix + ('Pxyz_c', '>f4', (3,)), # 78; mm from (0, 0, 0) RAS to vol center ] # Optional footer. Also has more stuff after this, optionally footer_dtd = [ - ('mrparms', '>f4', (4,)) + ('tr', '>f4'), # 0; repetition time + ('flip_angle', '>f4'), # 4; flip angle + ('te', '>f4'), # 8; echo time + ('ti', '>f4'), # 12; inversion time + ('fov', '>f4'), # 16; field of view (unused) ] header_dtype = np.dtype(header_dtd) @@ -70,7 +79,7 @@ class MGHError(Exception): """ -class MGHHeader(object): +class MGHHeader(LabeledWrapStruct): ''' Class for MGH format header The header also consists of the footer data which MGH places after the data @@ -96,64 +105,35 @@ def __init__(self, Whether to check content of header in initialization. Default is True. ''' - if binaryblock is None: - self._header_data = self._empty_headerdata() - return - # check size - if len(binaryblock) != self.template_dtype.itemsize: - raise HeaderDataError('Binary block is wrong size') - hdr = np.ndarray(shape=(), - dtype=self.template_dtype, - buffer=binaryblock) - # if goodRASFlag, discard delta, Mdc and c_ras stuff - if int(hdr['goodRASFlag']) < 0: - hdr = self._set_affine_default(hdr) - self._header_data = hdr.copy() + min_size = self._hdrdtype.itemsize + full_size = self.template_dtype.itemsize + if binaryblock is not None and len(binaryblock) >= min_size: + # Right zero-pad or truncate binaryblock to appropriate size + # Footer is optional and may contain variable-length text fields, + # so limit to fixed fields + binaryblock = (binaryblock[:full_size] + + b'\x00' * (full_size - len(binaryblock))) + super(MGHHeader, self).__init__(binaryblock=binaryblock, + endianness='big', + check=False) + if not self._structarr['goodRASFlag']: + self._set_affine_default() if check: self.check_fix() - return - def __str__(self): - ''' Print the MGH header object information - ''' - txt = [] - txt.append(str(self.__class__)) - txt.append('Dims: ' + str(self.get_data_shape())) - code = int(self._header_data['type']) - txt.append('MRI Type: ' + self._data_type_codes.mritype[code]) - txt.append('goodRASFlag: ' + str(self._header_data['goodRASFlag'])) - txt.append('delta: ' + str(self._header_data['delta'])) - txt.append('Mdc: ') - txt.append(str(self._header_data['Mdc'])) - txt.append('Pxyz_c: ' + str(self._header_data['Pxyz_c'])) - txt.append('mrparms: ' + str(self._header_data['mrparms'])) - return '\n'.join(txt) - - def __getitem__(self, item): - ''' Return values from header data - ''' - return self._header_data[item] - - def __setitem__(self, item, value): - ''' Set values in header data - ''' - self._header_data[item] = value - - def __iter__(self): - return iter(self.keys()) - - def keys(self): - ''' Return keys from header data''' - return list(self.template_dtype.names) - - def values(self): - ''' Return values from header data''' - data = self._header_data - return [data[key] for key in self.template_dtype.names] + @staticmethod + def chk_version(hdr, fix=False): + rep = Report() + if hdr['version'] != 1: + rep = Report(HeaderDataError, 40) + rep.problem_msg = 'Unknown MGH format version' + if fix: + hdr['version'] = 1 + return hdr, rep - def items(self): - ''' Return items from header data''' - return zip(self.keys(), self.values()) + @classmethod + def _get_checks(klass): + return (klass.chk_version,) @classmethod def from_header(klass, header=None, check=True): @@ -188,58 +168,19 @@ def from_fileobj(klass, fileobj, check=True): int(klass._data_type_codes.bytespervox[tp]) * np.prod(hdr_str_to_np['dims'])) ftr_str = fileobj.read(klass._ftrdtype.itemsize) - return klass(hdr_str + ftr_str, check) - - @property - def binaryblock(self): - ''' binary block of data as string - - Returns - ------- - binaryblock : string - string giving binary data block - - ''' - return self._header_data.tostring() - - def copy(self): - ''' Return copy of header - ''' - return self.__class__(self.binaryblock, check=False) - - def __eq__(self, other): - ''' equality between two MGH format headers - - Examples - -------- - >>> wstr = MGHHeader() - >>> wstr2 = MGHHeader() - >>> wstr == wstr2 - True - ''' - return self.binaryblock == other.binaryblock - - def __ne__(self, other): - return not self == other - - def check_fix(self): - ''' Pass. maybe for now''' + return klass(hdr_str + ftr_str, check=check) def get_affine(self): ''' Get the affine transform from the header information. + MGH format doesn't store the transform directly. Instead it's gleaned from the zooms ( delta ), direction cosines ( Mdc ), RAS centers ( Pxyz_c ) and the dimensions. ''' - hdr = self._header_data - d = np.diag(hdr['delta']) - pcrs_c = hdr['dims'][:3] / 2.0 - Mdc = hdr['Mdc'].T - pxyz_0 = hdr['Pxyz_c'] - np.dot(Mdc, np.dot(d, pcrs_c)) - M = np.eye(4, 4) - M[0:3, 0:3] = np.dot(Mdc, d) - M[0:3, 3] = pxyz_0.T - return M + hdr = self._structarr + MdcD = hdr['Mdc'].T * hdr['delta'] + vol_center = MdcD.dot(hdr['dims'][:3]) / 2 + return from_matvec(MdcD, hdr['Pxyz_c'] - vol_center) # For compatibility with nifti (multiple affines) get_best_affine = get_affine @@ -253,8 +194,8 @@ def get_vox2ras_tkr(self): ''' Get the vox2ras-tkr transform. See "Torig" here: https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems ''' - ds = np.array(self._header_data['delta']) - ns = (np.array(self._header_data['dims'][:3]) * ds) / 2.0 + ds = self._structarr['delta'] + ns = self._structarr['dims'][:3] * ds / 2.0 v2rtkr = np.array([[-ds[0], 0, 0, ns[0]], [0, 0, ds[2], -ns[2]], [0, -ds[1], 0, ns[1]], @@ -271,7 +212,7 @@ def get_data_dtype(self): For examples see ``set_data_dtype`` ''' - code = int(self._header_data['type']) + code = int(self._structarr['type']) dtype = self._data_type_codes.numpy_dtype[code] return dtype @@ -282,44 +223,74 @@ def set_data_dtype(self, datatype): code = self._data_type_codes[datatype] except KeyError: raise MGHError('datatype dtype "%s" not recognized' % datatype) - self._header_data['type'] = code + self._structarr['type'] = code + + def _ndims(self): + ''' Get dimensionality of data + + MGH does not encode dimensionality explicitly, so an image where the + fourth dimension is 1 is treated as three-dimensional. + + Returns + ------- + ndims : 3 or 4 + ''' + return 3 + (self._structarr['dims'][3] > 1) def get_zooms(self): ''' 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 + `_). + + To access only the spatial zooms, use `hdr['delta']`. + Returns ------- z : tuple tuple of header zoom values + + .. _mghformat: https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/MghFormat#line-82 ''' - hdr = self._header_data - zooms = hdr['delta'] - return tuple(zooms[:]) + # Do not return time zoom (TR) if 3D image + tzoom = (self['tr'],)[:self._ndims() > 3] + return tuple(self._structarr['delta']) + tzoom def set_zooms(self, zooms): ''' Set zooms into header fields - See docstring for ``get_zooms`` for examples + 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. + + Parameters + ---------- + zooms : sequence + sequence of floats specifying spatial and (optionally) temporal + zooms ''' - hdr = self._header_data + hdr = self._structarr zooms = np.asarray(zooms) - if len(zooms) != len(hdr['delta']): - raise HeaderDataError('Expecting %d zoom values for ndim' - % hdr['delta']) - if np.any(zooms < 0): + ndims = self._ndims() + if len(zooms) > ndims: + raise HeaderDataError('Expecting %d zoom values' % ndims) + if np.any(zooms <= 0): raise HeaderDataError('zooms must be positive') - delta = hdr['delta'] - delta[:] = zooms[:] + hdr['delta'] = zooms[:3] + if len(zooms) == 4: + hdr['tr'] = zooms[3] def get_data_shape(self): ''' Get shape of data ''' - dims = self._header_data['dims'][:] + shape = tuple(self._structarr['dims']) # If last dimension (nframes) is 1, remove it because # we want to maintain 3D and it's redundant - if int(dims[-1]) == 1: - dims = dims[:-1] - return tuple(int(d) for d in dims) + if shape[3] == 1: + shape = shape[:3] + return shape def set_data_shape(self, shape): ''' Set shape of data @@ -329,26 +300,22 @@ def set_data_shape(self, shape): shape : sequence sequence of integers specifying data array shape ''' - dims = self._header_data['dims'] - # If len(dims) is 3, add a dimension. MGH header always - # needs 4 dimensions. - if len(shape) == 3: - shape = list(shape) - shape.append(1) - shape = tuple(shape) - dims[:] = shape - self._header_data['delta'][:] = 1.0 + shape = tuple(shape) + if len(shape) > 4: + raise ValueError("Shape may be at most 4 dimensional") + self._structarr['dims'] = shape + (1,) * (4 - len(shape)) + self._structarr['delta'] = 1 def get_data_bytespervox(self): ''' Get the number of bytes per voxel of the data ''' return int(self._data_type_codes.bytespervox[ - int(self._header_data['type'])]) + int(self._structarr['type'])]) def get_data_size(self): ''' Get the number of bytes the data chunk occupies. ''' - return self.get_data_bytespervox() * np.prod(self._header_data['dims']) + return self.get_data_bytespervox() * np.prod(self._structarr['dims']) def get_data_offset(self): ''' Return offset into data file to read data @@ -384,32 +351,36 @@ def get_slope_inter(self): """ return None, None - def _empty_headerdata(self): + @classmethod + def guessed_endian(klass, mapping): + """ MGHHeader data must be big-endian """ + return '>' + + @classmethod + def default_structarr(klass, endianness=None): ''' Return header data for empty header + + Ignores byte order; always big endian ''' - dt = self.template_dtype - hdr_data = np.zeros((), dtype=dt) - hdr_data['version'] = 1 - hdr_data['dims'][:] = np.array([1, 1, 1, 1]) - hdr_data['type'] = 3 - hdr_data['goodRASFlag'] = 1 - hdr_data['delta'][:] = np.array([1, 1, 1]) - hdr_data['Mdc'][0][:] = np.array([-1, 0, 0]) # x_ras - hdr_data['Mdc'][1][:] = np.array([0, 0, -1]) # y_ras - hdr_data['Mdc'][2][:] = np.array([0, 1, 0]) # z_ras - hdr_data['Pxyz_c'] = np.array([0, 0, 0]) # c_ras - hdr_data['mrparms'] = np.array([0, 0, 0, 0]) - return hdr_data - - def _set_affine_default(self, hdr): - ''' If goodRASFlag is 0, return the default delta, Mdc and Pxyz_c + if endianness is not None and endian_codes[endianness] != '>': + raise ValueError('MGHHeader must always be big endian') + structarr = super(MGHHeader, + klass).default_structarr(endianness=endianness) + structarr['version'] = 1 + structarr['dims'] = 1 + structarr['type'] = 3 + structarr['goodRASFlag'] = 1 + structarr['delta'] = 1 + structarr['Mdc'] = [[-1, 0, 0], [0, 0, 1], [0, -1, 0]] + return structarr + + def _set_affine_default(self): + ''' If goodRASFlag is 0, set the default affine ''' - hdr['delta'][:] = np.array([1, 1, 1]) - hdr['Mdc'][0][:] = np.array([-1, 0, 0]) # x_ras - hdr['Mdc'][1][:] = np.array([0, 0, -1]) # y_ras - hdr['Mdc'][2][:] = np.array([0, 1, 0]) # z_ras - hdr['Pxyz_c'][:] = np.array([0, 0, 0]) # c_ras - return hdr + self._structarr['goodRASFlag'] = 1 + self._structarr['delta'] = 1 + self._structarr['Mdc'] = [[-1, 0, 0], [0, 0, 1], [0, -1, 0]] + self._structarr['Pxyz_c'] = 0 def writehdr_to(self, fileobj): ''' Write header to fileobj @@ -451,6 +422,81 @@ def writeftr_to(self, fileobj): fileobj.seek(self.get_footer_offset()) fileobj.write(ftr_nd.tostring()) + def copy(self): + ''' Return copy of structure ''' + return self.__class__(self.binaryblock, check=False) + + def as_byteswapped(self, endianness=None): + ''' Return new object with given ``endianness`` + + If big endian, returns a copy of the object. Otherwise raises ValueError. + + Parameters + ---------- + endianness : None or string, optional + endian code to which to swap. None means swap from current + endianness, and is the default + + Returns + ------- + wstr : ``MGHHeader`` + ``MGHHeader`` object + + ''' + if endianness is None or endian_codes[endianness] != '>': + raise ValueError('Cannot byteswap MGHHeader - ' + 'must always be big endian') + return self.copy() + + @classmethod + def diagnose_binaryblock(klass, binaryblock, endianness=None): + if endianness is not None and endian_codes[endianness] != '>': + raise ValueError('MGHHeader must always be big endian') + wstr = klass(binaryblock, check=False) + battrun = BatteryRunner(klass._get_checks()) + reports = battrun.check_only(wstr) + return '\n'.join([report.message + for report in reports if report.message]) + + class _HeaderData: + """ Provide interface to deprecated MGHHeader fields""" + def __init__(self, structarr): + self._structarr = structarr + + def __getitem__(self, item): + sa = self._structarr + if item == 'mrparams': + return np.hstack((sa['tr'], sa['flip_angle'], sa['te'], sa['ti'])) + return sa[item] + + def __setitem__(self, item, val): + sa = self._structarr + if item == 'mrparams': + sa['tr'], sa['flip_angle'], sa['te'], sa['ti'] = val + else: + sa[item] = val + + @property + @deprecate_with_version('_header_data is deprecated.\n' + 'Please use the _structarr interface instead.\n' + 'Note also that some fields have changed name and ' + 'shape.', + '2.3', '4.0') + def _header_data(self): + """ Deprecated field-access interface """ + return self._HeaderData(self._structarr) + + def __getitem__(self, item): + if item == 'mrparams': + return self._header_data[item] + return super(MGHHeader, self).__getitem__(item) + + def __setitem__(self, item, value): + if item == 'mrparams': + self._header_data[item] = value + else: + super(MGHHeader, self).__setitem__(item, value) + class MGHImage(SpatialImage): """ Class for MGH format image @@ -467,6 +513,14 @@ class MGHImage(SpatialImage): ImageArrayProxy = ArrayProxy + def __init__(self, dataobj, affine, header=None, + extra=None, file_map=None): + shape = dataobj.shape + if len(shape) < 3: + dataobj = reshape_dataobj(dataobj, shape + (1,) * (3 - len(shape))) + super(MGHImage, self).__init__(dataobj, affine, header=header, + extra=extra, file_map=file_map) + @classmethod def filespec_to_file_map(klass, filespec): """ Check for compressed .mgz format, then .mgh format """ @@ -516,9 +570,6 @@ def from_file_map(klass, file_map, mmap=True, keep_file_open=None): data = klass.ImageArrayProxy(img_fh.file_like, hdr_copy, mmap=mmap, keep_file_open=keep_file_open) img = klass(data, affine, header, file_map=file_map) - img._load_cache = {'header': hdr_copy, - 'affine': affine.copy(), - 'file_map': copy_file_map(file_map)} return img @classmethod @@ -579,7 +630,7 @@ def to_file_map(self, file_map=None): with file_map['image'].get_prepare_fileobj('wb') as mghf: hdr.writehdr_to(mghf) self._write_data(mghf, data, hdr) - self._write_footer(mghf, hdr) + hdr.writeftr_to(mghf) self._header = hdr self.file_map = file_map @@ -604,33 +655,20 @@ def _write_data(self, mghfile, data, header): out_dtype = header.get_data_dtype() array_to_file(data, mghfile, out_dtype, offset) - def _write_footer(self, mghfile, header): - ''' Utility routine to write header. This write the footer data - which occurs after the data chunk in mgh file - - Parameters - ---------- - mghfile : file-like - file-like object implementing ``write``, open for writing - header : header object - ''' - header.writeftr_to(mghfile) - def _affine2header(self): """ Unconditionally set affine into the header """ hdr = self._header - shape = self._dataobj.shape + shape = np.array(self._dataobj.shape[:3]) + # for more information, go through save_mgh.m in FreeSurfer dist - MdcD = self._affine[:3, :3] - delta = np.sqrt(np.sum(MdcD * MdcD, axis=0)) - Mdc = MdcD / np.tile(delta, (3, 1)) - Pcrs_c = np.array([0, 0, 0, 1], dtype=np.float) - Pcrs_c[:3] = np.array(shape[:3]) / 2.0 - Pxyz_c = np.dot(self._affine, Pcrs_c) - - hdr['delta'][:] = delta - hdr['Mdc'][:, :] = Mdc.T - hdr['Pxyz_c'][:] = Pxyz_c[:3] + voxelsize = voxel_sizes(self._affine) + Mdc = self._affine[:3, :3] / voxelsize + c_ras = self._affine.dot(np.hstack((shape / 2.0, [1])))[:3] + + # Assign after we've had a chance to raise exceptions + hdr['delta'] = voxelsize + hdr['Mdc'] = Mdc.T + hdr['Pxyz_c'] = c_ras load = MGHImage.load diff --git a/nibabel/freesurfer/tests/test_mghformat.py b/nibabel/freesurfer/tests/test_mghformat.py index b6a2e071ac..776c461e18 100644 --- a/nibabel/freesurfer/tests/test_mghformat.py +++ b/nibabel/freesurfer/tests/test_mghformat.py @@ -18,16 +18,22 @@ from ..mghformat import MGHHeader, MGHError, MGHImage from ...tmpdirs import InTemporaryDirectory from ...fileholders import FileHolder +from ...spatialimages import HeaderDataError +from ...volumeutils import sys_is_le +from ...wrapstruct import WrapStructError +from ... import imageglobals from nose.tools import assert_true, assert_false from numpy.testing import (assert_equal, assert_array_equal, assert_array_almost_equal, assert_almost_equal, assert_raises) +from ...testing import assert_not_equal from ...testing import data_path from ...tests import test_spatialimages as tsi +from ...tests.test_wrapstruct import _TestLabeledWrapStruct MGZ_FNAME = os.path.join(data_path, 'test.mgz') @@ -40,6 +46,17 @@ [0.0, -1.0, 0.0, 2.0], [0.0, 0.0, 0.0, 1.0]], dtype=np.float32) +BIG_CODES = ('>', 'big', 'BIG', 'b', 'be', 'B', 'BE') +LITTLE_CODES = ('<', 'little', 'l', 'le', 'L', 'LE') + +if sys_is_le: + BIG_CODES += ('swapped', 's', 'S', '!') + LITTLE_CODES += ('native', 'n', 'N', '=', '|', 'i', 'I') +else: + BIG_CODES += ('native', 'n', 'N', '=', '|', 'i', 'I') + LITTLE_CODES += ('swapped', 's', 'S', '!') + + def test_read_mgh(): # test.mgz was generated by the following command @@ -55,8 +72,11 @@ def test_read_mgh(): assert_equal(h['dof'], 0) assert_equal(h['goodRASFlag'], 1) assert_array_equal(h['dims'], [3, 4, 5, 2]) - assert_array_almost_equal(h['mrparms'], [2.0, 0.0, 0.0, 0.0]) - assert_array_almost_equal(h.get_zooms(), 1) + assert_almost_equal(h['tr'], 2.0) + 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_vox2ras(), v2r) assert_array_almost_equal(h.get_vox2ras_tkr(), v2rtkr) @@ -86,7 +106,11 @@ def test_write_mgh(): assert_equal(h['dof'], 0) assert_equal(h['goodRASFlag'], 1) assert_array_equal(h['dims'], [5, 4, 3, 2]) - assert_array_almost_equal(h['mrparms'], [0.0, 0.0, 0.0, 0.0]) + assert_almost_equal(h['tr'], 0.0) + assert_almost_equal(h['flip_angle'], 0.0) + assert_almost_equal(h['te'], 0.0) + assert_almost_equal(h['ti'], 0.0) + assert_almost_equal(h['fov'], 0.0) assert_array_almost_equal(h.get_vox2ras(), v2r) # data assert_almost_equal(dat, v, 7) @@ -112,15 +136,29 @@ def test_write_noaffine_mgh(): assert_equal(h['dof'], 0) assert_equal(h['goodRASFlag'], 1) assert_array_equal(h['dims'], [7, 13, 3, 22]) - assert_array_almost_equal(h['mrparms'], [0.0, 0.0, 0.0, 0.0]) + assert_almost_equal(h['tr'], 0.0) + assert_almost_equal(h['flip_angle'], 0.0) + assert_almost_equal(h['te'], 0.0) + assert_almost_equal(h['ti'], 0.0) + assert_almost_equal(h['fov'], 0.0) # important part -- whether default affine info is stored - ex_mdc = np.array([[-1, 0, 0], - [0, 0, -1], - [0, 1, 0]], dtype=np.float32) - assert_array_almost_equal(h['Mdc'], ex_mdc) + assert_array_almost_equal(h['Mdc'], [[-1, 0, 0], [0, 0, 1], [0, -1, 0]]) + assert_array_almost_equal(h['Pxyz_c'], [0, 0, 0]) + - ex_pxyzc = np.array([0, 0, 0], dtype=np.float32) - assert_array_almost_equal(h['Pxyz_c'], ex_pxyzc) +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]) + for zooms in ((-1, 1, 1, 1), + (1, -1, 1, 1), + (1, 1, -1, 1), + (1, 1, 1, -1), + (1, 1, 1, 1, 5)): + with assert_raises(HeaderDataError): + h.set_zooms(zooms) def bad_dtype_mgh(): @@ -179,7 +217,7 @@ def test_header_updating(): assert_almost_equal(hdr.get_affine(), exp_aff, 6) # Test that initial wonky header elements have not changed assert_equal(hdr['delta'], 1) - assert_almost_equal(hdr['Mdc'], exp_aff[:3, :3].T) + assert_almost_equal(hdr['Mdc'].T, exp_aff[:3, :3]) # Save, reload, same thing img_fobj = io.BytesIO() mgz2 = _mgh_rt(mgz, img_fobj) @@ -195,7 +233,7 @@ def test_header_updating(): assert_almost_equal(hdr2.get_affine(), exp_aff_d, 6) RZS = exp_aff_d[:3, :3] assert_almost_equal(hdr2['delta'], np.sqrt(np.sum(RZS ** 2, axis=0))) - assert_almost_equal(hdr2['Mdc'], (RZS / hdr2['delta']).T) + assert_almost_equal(hdr2['Mdc'].T, RZS / hdr2['delta']) def test_cosine_order(): @@ -210,7 +248,7 @@ def test_cosine_order(): hdr2 = img2.header RZS = aff[:3, :3] zooms = np.sqrt(np.sum(RZS ** 2, axis=0)) - assert_almost_equal(hdr2['Mdc'], (RZS / zooms).T) + assert_almost_equal(hdr2['Mdc'].T, RZS / zooms) assert_almost_equal(hdr2['delta'], zooms) @@ -251,6 +289,93 @@ def test_mgh_load_fileobj(): assert_array_equal(img.get_data(), img2.get_data()) +def test_mgh_affine_default(): + hdr = MGHHeader() + hdr['goodRASFlag'] = 0 + hdr2 = MGHHeader(hdr.binaryblock) + assert_equal(hdr2['goodRASFlag'], 1) + assert_array_equal(hdr['Mdc'], hdr2['Mdc']) + assert_array_equal(hdr['Pxyz_c'], hdr2['Pxyz_c']) + + +def test_mgh_set_data_shape(): + hdr = MGHHeader() + hdr.set_data_shape((5,)) + assert_array_equal(hdr.get_data_shape(), (5, 1, 1)) + hdr.set_data_shape((5, 4)) + assert_array_equal(hdr.get_data_shape(), (5, 4, 1)) + hdr.set_data_shape((5, 4, 3)) + assert_array_equal(hdr.get_data_shape(), (5, 4, 3)) + hdr.set_data_shape((5, 4, 3, 2)) + assert_array_equal(hdr.get_data_shape(), (5, 4, 3, 2)) + with assert_raises(ValueError): + hdr.set_data_shape((5, 4, 3, 2, 1)) + + +def test_mghheader_default_structarr(): + hdr = MGHHeader.default_structarr() + assert_equal(hdr['version'], 1) + assert_array_equal(hdr['dims'], 1) + assert_equal(hdr['type'], 3) + assert_equal(hdr['dof'], 0) + assert_equal(hdr['goodRASFlag'], 1) + assert_array_equal(hdr['delta'], 1) + assert_array_equal(hdr['Mdc'], [[-1, 0, 0], [0, 0, 1], [0, -1, 0]]) + assert_array_equal(hdr['Pxyz_c'], 0) + assert_equal(hdr['tr'], 0) + assert_equal(hdr['flip_angle'], 0) + assert_equal(hdr['te'], 0) + assert_equal(hdr['ti'], 0) + assert_equal(hdr['fov'], 0) + + for endianness in (None,) + BIG_CODES: + hdr2 = MGHHeader.default_structarr(endianness=endianness) + assert_equal(hdr2, hdr) + assert_equal(hdr2.newbyteorder('>'), hdr) + + for endianness in LITTLE_CODES: + with assert_raises(ValueError): + MGHHeader.default_structarr(endianness=endianness) + + +def test_deprecated_fields(): + hdr = MGHHeader() + hdr_data = MGHHeader._HeaderData(hdr.structarr) + + # mrparams is the only deprecated field at the moment + # Accessing hdr_data is equivalent to accessing hdr, so double all checks + assert_array_equal(hdr['mrparams'], 0) + assert_array_equal(hdr_data['mrparams'], 0) + + hdr['mrparams'] = [1, 2, 3, 4] + assert_array_almost_equal(hdr['mrparams'], [1, 2, 3, 4]) + assert_equal(hdr['tr'], 1) + assert_equal(hdr['flip_angle'], 2) + assert_equal(hdr['te'], 3) + assert_equal(hdr['ti'], 4) + assert_equal(hdr['fov'], 0) + assert_array_almost_equal(hdr_data['mrparams'], [1, 2, 3, 4]) + assert_equal(hdr_data['tr'], 1) + assert_equal(hdr_data['flip_angle'], 2) + assert_equal(hdr_data['te'], 3) + assert_equal(hdr_data['ti'], 4) + assert_equal(hdr_data['fov'], 0) + + hdr['tr'] = 5 + hdr['flip_angle'] = 6 + hdr['te'] = 7 + hdr['ti'] = 8 + assert_array_almost_equal(hdr['mrparams'], [5, 6, 7, 8]) + assert_array_almost_equal(hdr_data['mrparams'], [5, 6, 7, 8]) + + hdr_data['tr'] = 9 + hdr_data['flip_angle'] = 10 + hdr_data['te'] = 11 + hdr_data['ti'] = 12 + assert_array_almost_equal(hdr['mrparams'], [9, 10, 11, 12]) + assert_array_almost_equal(hdr_data['mrparams'], [9, 10, 11, 12]) + + class TestMGHImage(tsi.TestSpatialImage, tsi.MmapImageMixin): """ Apply general image tests to MGHImage """ @@ -262,3 +387,118 @@ def check_dtypes(self, expected, actual): # others may only require the same type # MGH requires the actual to be a big endian version of expected assert_equal(expected.newbyteorder('>'), actual) + + +class TestMGHHeader(_TestLabeledWrapStruct): + header_class = MGHHeader + + def _set_something_into_hdr(self, hdr): + hdr['dims'] = [4, 3, 2, 1] + + def get_bad_bb(self): + return b'\xff' + b'\x00' * self.header_class._hdrdtype.itemsize + + # Update tests to account for big-endian requirement + def test_general_init(self): + hdr = self.header_class() + # binaryblock has length given by header data dtype + binblock = hdr.binaryblock + assert_equal(len(binblock), hdr.structarr.dtype.itemsize) + # Endianness will always be big, and cannot be set + assert_equal(hdr.endianness, '>') + # You can also pass in a check flag, without data this has no + # effect + hdr = self.header_class(check=False) + + def test__eq__(self): + # Test equal and not equal + hdr1 = self.header_class() + hdr2 = self.header_class() + assert_equal(hdr1, hdr2) + self._set_something_into_hdr(hdr1) + assert_not_equal(hdr1, hdr2) + self._set_something_into_hdr(hdr2) + assert_equal(hdr1, hdr2) + # REMOVED as_byteswapped() test + # Check comparing to funny thing says no + assert_not_equal(hdr1, None) + assert_not_equal(hdr1, 1) + + def test_to_from_fileobj(self): + # Successful write using write_to + hdr = self.header_class() + str_io = io.BytesIO() + hdr.write_to(str_io) + str_io.seek(0) + hdr2 = self.header_class.from_fileobj(str_io) + assert_equal(hdr2.endianness, '>') + assert_equal(hdr2.binaryblock, hdr.binaryblock) + + def test_endian_guess(self): + # Check guesses of endian + eh = self.header_class() + assert_equal(eh.endianness, '>') + assert_equal(self.header_class.guessed_endian(eh), '>') + + def test_bytes(self): + # Test get of bytes + hdr1 = self.header_class() + bb = hdr1.binaryblock + hdr2 = self.header_class(hdr1.binaryblock) + assert_equal(hdr1, hdr2) + assert_equal(hdr1.binaryblock, hdr2.binaryblock) + # Do a set into the header, and try again. The specifics of 'setting + # something' will depend on the nature of the bytes object + self._set_something_into_hdr(hdr1) + hdr2 = self.header_class(hdr1.binaryblock) + assert_equal(hdr1, hdr2) + assert_equal(hdr1.binaryblock, hdr2.binaryblock) + # Short binaryblocks give errors (here set through init) + # Long binaryblocks are truncated + assert_raises(WrapStructError, + self.header_class, + bb[:self.header_class._hdrdtype.itemsize - 1]) + # Checking set to true by default, and prevents nonsense being + # set into the header. + bb_bad = self.get_bad_bb() + if bb_bad is None: + return + with imageglobals.LoggingOutputSuppressor(): + assert_raises(HeaderDataError, self.header_class, bb_bad) + # now slips past without check + _ = self.header_class(bb_bad, check=False) + + def test_as_byteswapped(self): + # Check byte swapping + hdr = self.header_class() + assert_equal(hdr.endianness, '>') + # same code just returns a copy + for endianness in BIG_CODES: + hdr2 = hdr.as_byteswapped(endianness) + assert_false(hdr2 is hdr) + assert_equal(hdr2, hdr) + + # Different code raises error + for endianness in (None,) + LITTLE_CODES: + with assert_raises(ValueError): + hdr.as_byteswapped(endianness) + # Note that contents is not rechecked on swap / copy + class DC(self.header_class): + def check_fix(self, *args, **kwargs): + raise Exception + + # Assumes check=True default + assert_raises(Exception, DC, hdr.binaryblock) + hdr = DC(hdr.binaryblock, check=False) + hdr2 = hdr.as_byteswapped('>') + + def test_checks(self): + # Test header checks + hdr_t = self.header_class() + # _dxer just returns the diagnostics as a string + # Default hdr is OK + assert_equal(self._dxer(hdr_t), '') + # Version should be 1 + hdr = hdr_t.copy() + hdr['version'] = 2 + assert_equal(self._dxer(hdr), 'Unknown MGH format version') diff --git a/nibabel/spatialimages.py b/nibabel/spatialimages.py index 2554600649..b88b3e8538 100644 --- a/nibabel/spatialimages.py +++ b/nibabel/spatialimages.py @@ -354,7 +354,7 @@ def __init__(self, dataobj, affine, header=None, ''' super(SpatialImage, self).__init__(dataobj, header=header, extra=extra, file_map=file_map) - if not affine is None: + if affine is not None: # Check that affine is array-like 4,4. Maybe this is too strict at # this abstract level, but so far I think all image formats we know # do need 4,4. diff --git a/nibabel/tests/test_spatialimages.py b/nibabel/tests/test_spatialimages.py index 9cb4759f50..bd8b834b84 100644 --- a/nibabel/tests/test_spatialimages.py +++ b/nibabel/tests/test_spatialimages.py @@ -195,7 +195,7 @@ class DataLike(object): shape = (3,) def __array__(self): - return np.arange(3) + return np.arange(3, dtype=np.int16) class TestSpatialImage(TestCase): @@ -249,8 +249,11 @@ def test_default_header(self): def test_data_api(self): # Test minimal api data object can initialize img = self.image_class(DataLike(), None) - assert_array_equal(img.get_data(), np.arange(3)) - assert_equal(img.shape, (3,)) + # Shape may be promoted to higher dimension, but may not reorder or + # change size + assert_array_equal(img.get_data().flatten(), np.arange(3)) + assert_equal(img.get_shape()[:1], (3,)) + assert_equal(np.prod(img.get_shape()), 3) def check_dtypes(self, expected, actual): # Some images will want dtypes to be equal including endianness, @@ -278,7 +281,10 @@ def test_data_shape(self): # See https://github.com/nipy/nibabel/issues/58 arr = np.arange(4, dtype=np.int16) img = img_klass(arr, np.eye(4)) - assert_equal(img.shape, (4,)) + # Shape may be promoted to higher dimension, but may not reorder or + # change size + assert_equal(img.get_shape()[:1], (4,)) + assert_equal(np.prod(img.get_shape()), 4) img = img_klass(np.zeros((2, 3, 4), dtype=np.float32), np.eye(4)) assert_equal(img.shape, (2, 3, 4)) @@ -290,7 +296,10 @@ def test_str(self): arr = np.arange(5, dtype=np.int16) img = img_klass(arr, np.eye(4)) assert_true(len(str(img)) > 0) - assert_equal(img.shape, (5,)) + # Shape may be promoted to higher dimension, but may not reorder or + # change size + assert_equal(img.shape[:1], (5,)) + assert_equal(np.prod(img.shape), 5) img = img_klass(np.zeros((2, 3, 4), dtype=np.int16), np.eye(4)) assert_true(len(str(img)) > 0) @@ -302,7 +311,10 @@ def test_get_shape(self): # See https://github.com/nipy/nibabel/issues/58 img = img_klass(np.arange(1, dtype=np.int16), np.eye(4)) with suppress_warnings(): - assert_equal(img.get_shape(), (1,)) + # Shape may be promoted to higher dimension, but may not reorder or + # change size + assert_equal(img.get_shape()[:1], (1,)) + assert_equal(np.prod(img.get_shape()), 1) img = img_klass(np.zeros((2, 3, 4), np.int16), np.eye(4)) assert_equal(img.get_shape(), (2, 3, 4))