diff --git a/nibabel/cifti2/__init__.py b/nibabel/cifti2/__init__.py index 3025a6f991..0c80e4033b 100644 --- a/nibabel/cifti2/__init__.py +++ b/nibabel/cifti2/__init__.py @@ -25,3 +25,4 @@ Cifti2TransformationMatrixVoxelIndicesIJKtoXYZ, Cifti2Vertices, Cifti2Volume, CIFTI_BRAIN_STRUCTURES, Cifti2HeaderError, CIFTI_MODEL_TYPES, load, save) +from .cifti2_axes import (Axis, BrainModel, Parcels, Series, Label, Scalar) \ No newline at end of file diff --git a/nibabel/cifti2/cifti2.py b/nibabel/cifti2/cifti2.py index 0ffe45a169..04c44433b5 100644 --- a/nibabel/cifti2/cifti2.py +++ b/nibabel/cifti2/cifti2.py @@ -1268,6 +1268,40 @@ def get_index_map(self, index): ''' return self.matrix.get_index_map(index) + def get_axis(self, index): + ''' + Generates the Cifti2 axis for a given dimension + + Parameters + ---------- + index : int + Dimension for which we want to obtain the mapping. + + Returns + ------- + axis : cifti2_axes.Axis + ''' + from . import cifti2_axes + return cifti2_axes.from_mapping(self.matrix.get_index_map(index)) + + @classmethod + def from_axes(cls, axes): + ''' + Creates a new Cifti2 header based on the Cifti2 axes + + Parameters + ---------- + axes : Tuple[cifti2_axes.Axis] + sequence of Cifti2 axes describing each row/column of the matrix to be stored + + Returns + ------- + header : Cifti2Header + new header describing the rows/columns in a format consistent with Cifti2 + ''' + from . import cifti2_axes + return cifti2_axes.to_header(axes) + class Cifti2Image(DataobjImage): """ Class for single file CIFTI2 format image @@ -1297,8 +1331,10 @@ def __init__(self, Object containing image data. It should be some object that returns an array from ``np.asanyarray``. It should have a ``shape`` attribute or property. - header : Cifti2Header instance + header : Cifti2Header instance or Sequence[cifti2_axes.Axis] Header with data for / from XML part of CIFTI2 format. + Alternatively a sequence of cifti2_axes.Axis objects can be provided + describing each dimension of the array. nifti_header : None or mapping or NIfTI2 header instance, optional Metadata for NIfTI2 component of this format. extra : None or mapping @@ -1306,6 +1342,8 @@ def __init__(self, file_map : mapping, optional Mapping giving file information for this image format. ''' + if not isinstance(header, Cifti2Header) and header: + header = Cifti2Header.from_axes(header) super(Cifti2Image, self).__init__(dataobj, header=header, extra=extra, file_map=file_map) self._nifti_header = Nifti2Header.from_header(nifti_header) diff --git a/nibabel/cifti2/cifti2_axes.py b/nibabel/cifti2/cifti2_axes.py new file mode 100644 index 0000000000..09faf0e341 --- /dev/null +++ b/nibabel/cifti2/cifti2_axes.py @@ -0,0 +1,1224 @@ +import numpy as np +from nibabel.cifti2 import cifti2 +from six import string_types +from operator import xor + + +def from_mapping(mim): + """ + Parses the MatrixIndicesMap to find the appropriate CIFTI axis describing the rows or columns + + Parameters + ---------- + mim : cifti2.Cifti2MatrixIndicesMap + + Returns + ------- + subtype of Axis + """ + return_type = {'CIFTI_INDEX_TYPE_SCALARS': Scalar, + 'CIFTI_INDEX_TYPE_LABELS': Label, + 'CIFTI_INDEX_TYPE_SERIES': Series, + 'CIFTI_INDEX_TYPE_BRAIN_MODELS': BrainModel, + 'CIFTI_INDEX_TYPE_PARCELS': Parcels} + return return_type[mim.indices_map_to_data_type].from_mapping(mim) + + +def to_header(axes): + """ + Converts the axes describing the rows/columns of a CIFTI vector/matrix to a Cifti2Header + + Parameters + ---------- + axes : iterable[Axis] + one or more axes describing each dimension in turn + + Returns + ------- + cifti2.Cifti2Header + """ + axes = list(axes) + mims_all = [] + matrix = cifti2.Cifti2Matrix() + for dim, ax in enumerate(axes): + if ax in axes[:dim]: + dim_prev = axes.index(ax) + mims_all[dim_prev].applies_to_matrix_dimension.append(dim) + mims_all.append(mims_all[dim_prev]) + else: + mim = ax.to_mapping(dim) + mims_all.append(mim) + matrix.append(mim) + return cifti2.Cifti2Header(matrix) + + +class Axis(object): + """ + Abstract class for any object describing the rows or columns of a CIFTI vector/matrix + + Attributes + ---------- + arr : np.ndarray + (N, ) typed array with the actual information on each row/column + """ + _use_dtype = None + arr = None + + def __init__(self, arr): + self.arr = np.asarray(arr, dtype=self._use_dtype) + + def get_element(self, index): + """ + Extracts a single element from the axis + + Parameters + ---------- + index : int + Indexes the row/column of interest + + Returns + ------- + Description of the row/column + """ + return self.arr[index] + + def __getitem__(self, item): + if isinstance(item, int): + return self.get_element(item) + if isinstance(item, string_types): + raise IndexError("Can not index an Axis with a string (except for Parcels)") + return type(self)(self.arr[item]) + + @property + def size(self, ): + return self.arr.size + + def __len__(self): + return self.size + + def __eq__(self, other): + return (type(self) == type(other) and + len(self) == len(other) and + (self.arr == other.arr).all()) + + def __add__(self, other): + """ + Concatenates two Axes of the same type + + Parameters + ---------- + other : Axis + axis to be appended to the current one + + Returns + ------- + Axis of the same subtype as self and other + """ + if type(self) == type(other): + return type(self)(np.append(self.arr, other.arr)) + return NotImplemented + + +class BrainModel(Axis): + """ + Each row/column in the CIFTI vector/matrix represents a single vertex or voxel + + This Axis describes which vertex/voxel is represented by each row/column. + + Attributes + ---------- + voxel : np.ndarray + (N, 3) array with the voxel indices + vertex : np.ndarray + (N, ) array with the vertex indices + name : np.ndarray + (N, ) array with the brain structure objects + """ + _use_dtype = np.dtype([('vertex', 'i4'), ('voxel', ('i4', 3)), + ('name', 'U%i' % max(len(name) for name in cifti2.CIFTI_BRAIN_STRUCTURES))]) + _affine = None + _volume_shape = None + + def __init__(self, arr, affine=None, volume_shape=None, nvertices=None): + """ + Creates a new BrainModel axis defining the vertices and voxels represented by each row/column + + Parameters + ---------- + arr : np.ndarray + (N, ) structured array with for every element a tuple with 3 elements: + - vertex index (-1 for voxels) + - 3 voxel indices (-1 for vertices) + - string (name of brain structure) + affine : np.ndarray + (4, 4) array mapping voxel indices to mm space (not needed for CIFTI files only covering the surface) + volume_shape : Tuple[int, int, int] + shape of the volume in which the voxels were defined (not needed for CIFTI files only covering the surface) + nvertices : dict[String -> int] + maps names of surface elements to integers + """ + super(BrainModel, self).__init__(arr) + self.name = self.name # correct names to CIFTI brain structures + if nvertices is None: + self.nvertices = {} + else: + self.nvertices = dict(nvertices) + for name in list(self.nvertices.keys()): + if name not in self.name: + del self.nvertices[name] + if self.is_surface.all(): + self.affine = None + self.volume_shape = None + else: + self.affine = affine + self.volume_shape = volume_shape + + @classmethod + def from_mapping(cls, mim): + """ + Creates a new BrainModel axis based on a CIFTI dataset + + Parameters + ---------- + mim : cifti2.Cifti2MatrixIndicesMap + + Returns + ------- + BrainModel + """ + nbm = np.sum([bm.index_count for bm in mim.brain_models]) + arr = np.zeros(nbm, dtype=cls._use_dtype) + arr['voxel'] = -1 + arr['vertex'] = -1 + nvertices = {} + affine, shape = None, None + for bm in mim.brain_models: + index_end = bm.index_offset + bm.index_count + is_surface = bm.model_type == 'CIFTI_MODEL_TYPE_SURFACE' + arr['name'][bm.index_offset: index_end] = bm.brain_structure + if is_surface: + arr['vertex'][bm.index_offset: index_end] = bm.vertex_indices + nvertices[bm.brain_structure] = bm.surface_number_of_vertices + else: + arr['voxel'][bm.index_offset: index_end, :] = bm.voxel_indices_ijk + if affine is None: + shape = mim.volume.volume_dimensions + affine = mim.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix + else: + if shape != mim.volume.volume_dimensions: + raise ValueError("All volume masks should be defined in the same volume") + if (affine != mim.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix).any(): + raise ValueError("All volume masks should have the same affine") + return cls(arr, affine, shape, nvertices) + + @classmethod + def from_mask(cls, mask, name='other', affine=None): + """ + Creates a new BrainModel axis describing the provided mask + + Parameters + ---------- + mask : np.ndarray + all non-zero voxels will be included in the BrainModel axis + should be (Nx, Ny, Nz) array for volume mask or (Nvertex, ) array for surface mask + name : str + Name of the brain structure (e.g. 'CortexRight', 'thalamus_left' or 'brain_stem') + affine : np.ndarray + (4, 4) array with the voxel to mm transformation (defaults to identity matrix) + Argument will be ignored for surface masks + + Returns + ------- + BrainModel which covers the provided mask + """ + if affine is None: + affine = np.eye(4) + if np.asarray(affine).shape != (4, 4): + raise ValueError("Affine transformation should be a 4x4 array or None, not %r" % affine) + if mask.ndim == 1: + return cls.from_surface(np.where(mask != 0)[0], mask.size, name=name) + elif mask.ndim == 3: + voxels = np.array(np.where(mask != 0)).T + arr = np.zeros(len(voxels), dtype=cls._use_dtype) + arr['vertex'] = -1 + arr['voxel'] = voxels + arr['name'] = cls.to_cifti_brain_structure_name(name) + return cls(arr, affine=affine, volume_shape=mask.shape) + else: + raise ValueError("Mask should be either 1-dimensional (for surfaces) or " + "3-dimensional (for volumes), not %i-dimensional" % mask.ndim) + + @classmethod + def from_surface(cls, vertices, nvertex, name='Cortex'): + """ + Creates a new BrainModel axis describing the vertices on a surface + + Parameters + ---------- + vertices : np.ndarray + indices of the vertices on the surface + nvertex : int + total number of vertices on the surface + name : str + Name of the brain structure (e.g. 'CortexLeft' or 'CortexRight') + + Returns + ------- + BrainModel which covers (part of) the surface + """ + arr = np.zeros(len(vertices), dtype=cls._use_dtype) + arr['voxel'] = -1 + arr['vertex'] = vertices + arr['name'] = cls.to_cifti_brain_structure_name(name) + return cls(arr, nvertices={arr['name'][0]: nvertex}) + + def get_element(self, index): + """ + Describes a single element from the axis + + Parameters + ---------- + index : int + Indexes the row/column of interest + + Returns + ------- + tuple with 3 elements + - boolean, which is True if it is a surface element + - vertex index if it is a surface element, otherwise array with 3 voxel indices + - structure.BrainStructure object describing the brain structure the element was taken from + """ + elem = self.arr[index] + is_surface = elem['name'] in self.nvertices.keys() + name = 'vertex' if is_surface else 'voxel' + return is_surface, elem[name], elem['name'] + + def to_mapping(self, dim): + """ + Converts the brain model axis to a MatrixIndicesMap for storage in CIFTI format + + Parameters + ---------- + dim : int + which dimension of the CIFTI vector/matrix is described by this dataset (zero-based) + + Returns + ------- + cifti2.Cifti2MatrixIndicesMap + """ + mim = cifti2.Cifti2MatrixIndicesMap([dim], 'CIFTI_INDEX_TYPE_BRAIN_MODELS') + for name, to_slice, bm in self.iter_structures(): + is_surface = name in self.nvertices.keys() + if is_surface: + voxels = None + vertices = cifti2.Cifti2VertexIndices(bm.vertex) + nvertex = self.nvertices[name] + else: + voxels = cifti2.Cifti2VoxelIndicesIJK(bm.voxel) + vertices = None + nvertex = None + if mim.volume is None: + affine = cifti2.Cifti2TransformationMatrixVoxelIndicesIJKtoXYZ(-3, matrix=self.affine) + mim.volume = cifti2.Cifti2Volume(self.volume_shape, affine) + cifti_bm = cifti2.Cifti2BrainModel(to_slice.start, len(bm), + 'CIFTI_MODEL_TYPE_SURFACE' if is_surface else 'CIFTI_MODEL_TYPE_VOXELS', + name, nvertex, voxels, vertices) + mim.append(cifti_bm) + return mim + + def iter_structures(self, ): + """ + Iterates over all brain structures in the order that they appear along the axis + + Yields + ------ + tuple with + - CIFTI brain structure name + - slice to select the data associated with the brain structure from the tensor + - brain model covering that specific brain structure + """ + idx_start = 0 + start_name = self.name[idx_start] + for idx_current, name in enumerate(self.name): + if start_name != name: + yield start_name, slice(idx_start, idx_current), self[idx_start: idx_current] + idx_start = idx_current + start_name = self.name[idx_start] + yield start_name, slice(idx_start, None), self[idx_start:] + + @property + def affine(self, ): + return self._affine + + @affine.setter + def affine(self, value): + if value is not None: + value = np.asarray(value) + if value.shape != (4, 4): + raise ValueError('Affine transformation should be a 4x4 array') + self._affine = value + + @property + def volume_shape(self, ): + return self._volume_shape + + @volume_shape.setter + def volume_shape(self, value): + if value is not None: + value = tuple(value) + if len(value) != 3: + raise ValueError("Volume shape should be a tuple of length 3") + self._volume_shape = value + + @property + def is_surface(self, ): + """True for any element on the surface + """ + return np.vectorize(lambda name: name in self.nvertices.keys())(self.name) + + @property + def voxel(self, ): + """The voxel represented by each row or column + """ + return self.arr['voxel'] + + @voxel.setter + def voxel(self, values): + self.arr['voxel'] = values + + @property + def vertex(self, ): + """The vertex represented by each row or column + """ + return self.arr['vertex'] + + @vertex.setter + def vertex(self, values): + self.arr['vertex'] = values + + @property + def name(self, ): + """The brain structure to which the voxel/vertices of belong + """ + return self.arr['name'] + + @name.setter + def name(self, values): + self.arr['name'] = [self.to_cifti_brain_structure_name(name) for name in values] + + @staticmethod + def to_cifti_brain_structure_name(name): + """ + Attempts to convert the name of an anatomical region in a format recognized by CIFTI + + This function returns: + * the name if it is in the CIFTI format already + * if the name is a tuple the first element is assumed to be the structure name while + the second is assumed to be the hemisphere (left, right or both). The latter will default + to both. + * names like left_cortex, cortex_left, LeftCortex, or CortexLeft will be converted to + CIFTI_STRUCTURE_CORTEX_LEFT + + see ``nibabel.cifti2.tests.test_name`` for examples of which conversions are possible + + Parameters + ---------- + name: (str, tuple) + input name of an anatomical region + + Returns + ------- + CIFTI2 compatible name + + Raises + ------ + ValueError: raised if the input name does not match a known anatomical structure in CIFTI + """ + if name in cifti2.CIFTI_BRAIN_STRUCTURES: + return name + if not isinstance(name, string_types): + if len(name) == 1: + structure = name[0] + orientation = 'both' + else: + structure, orientation = name + if structure.lower() in ('left', 'right', 'both'): + orientation, structure = name + else: + orient_names = ('left', 'right', 'both') + for poss_orient in orient_names: + idx = len(poss_orient) + if poss_orient == name.lower()[:idx]: + orientation = poss_orient + if name[idx] in '_ ': + structure = name[idx + 1:] + else: + structure = name[idx:] + break + if poss_orient == name.lower()[-idx:]: + orientation = poss_orient + if name[-idx - 1] in '_ ': + structure = name[:-idx - 1] + else: + structure = name[:-idx] + break + else: + orientation = 'both' + structure = name + if orientation.lower() == 'both': + proposed_name = 'CIFTI_STRUCTURE_%s' % structure.upper() + else: + proposed_name = 'CIFTI_STRUCTURE_%s_%s' % (structure.upper(), orientation.upper()) + if proposed_name not in cifti2.CIFTI_BRAIN_STRUCTURES: + raise ValueError('%s was interpreted as %s, which is not a valid CIFTI brain structure' % + (name, proposed_name)) + return proposed_name + + def __getitem__(self, item): + if isinstance(item, int): + return self.get_element(item) + if isinstance(item, string_types): + raise IndexError("Can not index an Axis with a string (except for Parcels)") + return type(self)(self.arr[item], self.affine, self.volume_shape, self.nvertices) + + def __eq__(self, other): + if type(self) != type(other) or len(self) != len(other): + return False + if xor(self.affine is None, other.affine is None): + return False + return (((self.affine is None and other.affine is None) or + (abs(self.affine - other.affine).max() < 1e-8 and + self.volume_shape == other.volume_shape)) and + (self.nvertices == other.nvertices) and + (self.arr == other.arr).all()) + + def __add__(self, other): + """ + Concatenates two BrainModels + + Parameters + ---------- + other : BrainModel + brain model to be appended to the current one + + Returns + ------- + BrainModel + """ + if type(self) == type(other): + if self.affine is None: + affine, shape = other.affine, other.volume_shape + else: + affine, shape = self.affine, self.volume_shape + if other.affine is not None and ((other.affine != affine).all() or + other.volume_shape != shape): + raise ValueError("Trying to concatenate two BrainModels defined in a different brain volume") + nvertices = dict(self.nvertices) + for name, value in other.nvertices.items(): + if name in nvertices.keys() and nvertices[name] != value: + raise ValueError("Trying to concatenate two BrainModels with inconsistent number of vertices for %s" + % name) + nvertices[name] = value + return type(self)(np.append(self.arr, other.arr), affine, shape, nvertices) + return NotImplemented + + +class Parcels(Axis): + """ + Each row/column in the CIFTI vector/matrix represents a parcel of voxels/vertices + + This Axis describes which parcel is represented by each row/column. + + Attributes + ---------- + name : np.ndarray + (N, ) string array with the parcel names + parcel : np.ndarray + (N, ) array with the actual get_parcels (each of which is a BrainModel object) + + Individual get_parcels can also be accessed based on their name, using + ``parcel = parcel_axis[name]`` + """ + _use_dtype = np.dtype([('name', 'U60'), ('voxels', 'object'), ('vertices', 'object')]) + _affine = None + _volume_shape = None + + def __init__(self, arr, affine=None, volume_shape=None, nvertices=None): + """ + Creates a new BrainModel axis defining the vertices and voxels represented by each row/column + + Parameters + ---------- + arr : np.ndarray + (N, ) structured array with for every element a tuple with 3 elements: + - string (name of parcel) + - (M, 3) int array with the M voxel indices in the parcel + - Dict[String -> (K, ) int array] mapping surface brain structure names to vertex indices + affine : np.ndarray + (4, 4) array mapping voxel indices to mm space (not needed for CIFTI files only covering the surface) + volume_shape : Tuple[int, int, int] + shape of the volume in which the voxels were defined (not needed for CIFTI files only covering the surface) + nvertices : dict[String -> int] + maps names of surface elements to integers + """ + super(Parcels, self).__init__(arr) + self.affine = affine + self.volume_shape = volume_shape + if nvertices is None: + self.nvertices = {} + else: + self.nvertices = dict(nvertices) + + @classmethod + def from_brain_models(cls, named_brain_models): + """ + Creates a Parcel axis from a list of BrainModel axes with names + + Parameters + ---------- + named_brain_models : List[Tuple[String, BrainModel]] + list of (parcel name, brain model representation) pairs defining each parcel + + Returns + ------- + Parcels + """ + affine = None + volume_shape = None + arr = np.zeros(len(named_brain_models), dtype=cls._use_dtype) + nvertices = {} + for idx_parcel, (parcel_name, bm) in enumerate(named_brain_models): + voxels = bm.voxel[~bm.is_surface] + if voxels.shape[0] != 0: + if affine is None: + affine = bm.affine + volume_shape = bm.volume_shape + else: + if (affine != bm.affine).any() or (volume_shape != bm.volume_shape): + raise ValueError( + "Can not combine brain models defined in different volumes into a single Parcel axis") + vertices = {} + for name, _, bm_part in bm.iter_structures(): + if name in bm.nvertices.keys(): + if name in nvertices.keys() and nvertices[name] != bm.nvertices[name]: + raise ValueError("Got multiple conflicting number of vertices for surface structure %s" % name) + nvertices[name] = bm.nvertices[name] + vertices[name] = bm_part.vertex + arr[idx_parcel] = (parcel_name, voxels, vertices) + return Parcels(arr, affine, volume_shape, nvertices) + + @classmethod + def from_mapping(cls, mim): + """ + Creates a new Parcels axis based on a CIFTI dataset + + Parameters + ---------- + mim : cifti2.Cifti2MatrixIndicesMap + + Returns + ------- + Parcels + """ + nparcels = len(list(mim.parcels)) + arr = np.zeros(nparcels, dtype=cls._use_dtype) + volume_shape = None if mim.volume is None else mim.volume.volume_dimensions + affine = None if mim.volume is None else mim.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix + nvertices = {} + for surface in mim.surfaces: + nvertices[surface.brain_structure] = surface.surface_number_of_vertices + for idx_parcel, parcel in enumerate(mim.parcels): + nvoxels = 0 if parcel.voxel_indices_ijk is None else len(parcel.voxel_indices_ijk) + voxels = np.zeros((nvoxels, 3), dtype='i4') + if nvoxels != 0: + voxels[()] = parcel.voxel_indices_ijk + vertices = {} + for vertex in parcel.vertices: + name = vertex.brain_structure + vertices[vertex.brain_structure] = np.array(vertex) + if name not in nvertices.keys(): + raise ValueError("Number of vertices for surface structure %s not defined" % name) + arr[idx_parcel]['voxels'] = voxels + arr[idx_parcel]['vertices'] = vertices + arr[idx_parcel]['name'] = parcel.name + return cls(arr, affine, volume_shape, nvertices) + + def to_mapping(self, dim): + """ + Converts the get_parcels to a MatrixIndicesMap for storage in CIFTI format + + Parameters + ---------- + dim : int + which dimension of the CIFTI vector/matrix is described by this dataset (zero-based) + + Returns + ------- + cifti2.Cifti2MatrixIndicesMap + """ + mim = cifti2.Cifti2MatrixIndicesMap([dim], 'CIFTI_INDEX_TYPE_PARCELS') + if self.affine is not None: + affine = cifti2.Cifti2TransformationMatrixVoxelIndicesIJKtoXYZ(-3, matrix=self.affine) + mim.volume = cifti2.Cifti2Volume(self.volume_shape, affine) + for name, nvertex in self.nvertices.items(): + mim.append(cifti2.Cifti2Surface(name, nvertex)) + for name, voxels, vertices in self.arr: + cifti_voxels = cifti2.Cifti2VoxelIndicesIJK(voxels) + element = cifti2.Cifti2Parcel(name, cifti_voxels) + for name, idx_vertices in vertices.items(): + element.vertices.append(cifti2.Cifti2Vertices(name, idx_vertices)) + mim.append(element) + return mim + + def get_element(self, index): + """ + Describes a single element from the axis + + Parameters + ---------- + index : int + Indexes the row/column of interest + + Returns + ------- + tuple with 3 elements + - unicode name of the parcel + - (M, 3) int array with voxel indices + - Dict[String -> (K, ) int array] with vertex indices for a specific surface brain structure + """ + return self.name[index], self.voxels[index], self.vertices[index] + + @property + def affine(self, ): + return self._affine + + @affine.setter + def affine(self, value): + if value is not None: + value = np.asarray(value) + if value.shape != (4, 4): + raise ValueError('Affine transformation should be a 4x4 array') + self._affine = value + + @property + def volume_shape(self, ): + return self._volume_shape + + @volume_shape.setter + def volume_shape(self, value): + if value is not None: + value = tuple(value) + if len(value) != 3: + raise ValueError("Volume shape should be a tuple of length 3") + self._volume_shape = value + + @property + def name(self, ): + return self.arr['name'] + + @name.setter + def name(self, values): + self.arr['name'] = values + + @property + def voxels(self, ): + return self.arr['voxels'] + + @voxels.setter + def voxels(self, values): + self.arr['voxels'] = values + + @property + def vertices(self, ): + return self.arr['vertices'] + + @vertices.setter + def vertices(self, values): + self.arr['vertices'] = values + + def __getitem__(self, item): + if isinstance(item, string_types): + idx = np.where(self.name == item)[0] + if len(idx) == 0: + raise IndexError("Parcel %s not found" % item) + if len(idx) > 1: + raise IndexError("Multiple get_parcels with name %s found" % item) + return self.voxels[idx[0]], self.vertices[idx[0]] + if isinstance(item, int): + return self.get_element(item) + if isinstance(item, string_types): + raise IndexError("Can not index an Axis with a string (except for Parcels)") + return type(self)(self.arr[item], self.affine, self.volume_shape, self.nvertices) + + def __eq__(self, other): + if (type(self) != type(other) or len(self) != len(other) or + (self.name != other.name).all() or self.nvertices != other.nvertices or + any((vox1 != vox2).any() for vox1, vox2 in zip(self.voxels, other.voxels))): + return False + if self.affine is not None: + if ( other.affine is None or + abs(self.affine - other.affine).max() > 1e-8 or + self.volume_shape != other.volume_shape): + return False + elif other.affine is not None: + return False + for vert1, vert2 in zip(self.vertices, other.vertices): + if len(vert1) != len(vert2): + return False + for name in vert1.keys(): + if name not in vert2 or (vert1[name] != vert2[name]).all(): + return False + return True + + def __add__(self, other): + """ + Concatenates two Parcels + + Parameters + ---------- + other : Parcel + parcel to be appended to the current one + + Returns + ------- + Parcel + """ + if type(self) == type(other): + if self.affine is None: + affine, shape = other.affine, other.volume_shape + else: + affine, shape = self.affine, self.volume_shape + if other.affine is not None and ((other.affine != affine).all() or + other.volume_shape != shape): + raise ValueError("Trying to concatenate two Parcels defined in a different brain volume") + nvertices = dict(self.nvertices) + for name, value in other.nvertices.items(): + if name in nvertices.keys() and nvertices[name] != value: + raise ValueError("Trying to concatenate two Parcels with inconsistent number of vertices for %s" + % name) + nvertices[name] = value + return type(self)(np.append(self.arr, other.arr), affine, shape, nvertices) + return NotImplemented + + +class Scalar(Axis): + """ + Along this axis of the CIFTI vector/matrix each row/column has been given a unique name and optionally metadata + + Attributes + ---------- + name : np.ndarray + (N, ) string array with the parcel names + meta : np.ndarray + (N, ) array with a dictionary of metadata for each row/column + """ + _use_dtype = np.dtype([('name', 'U60'), ('meta', 'object')]) + + def __init__(self, arr): + """ + Creates a new Scalar axis from (name, meta-data) pairs + + Parameters + ---------- + arr : Iterable[Tuple[str, dict[str -> str]] + iterates over all rows/columns assigning a name and a dictionary of metadata to each + """ + super(Scalar, self).__init__(arr) + + @classmethod + def from_mapping(cls, mim): + """ + Creates a new get_scalar axis based on a CIFTI dataset + + Parameters + ---------- + mim : cifti2.Cifti2MatrixIndicesMap + + Returns + ------- + Scalar + """ + res = np.zeros(len(list(mim.named_maps)), dtype=cls._use_dtype) + res['name'] = [nm.map_name for nm in mim.named_maps] + res['meta'] = [{} if nm.metadata is None else dict(nm.metadata) for nm in mim.named_maps] + return cls(res) + + @classmethod + def from_names(cls, names): + """ + Creates a new get_scalar axis with the given row/column names + + Parameters + ---------- + names : List[str] + gives a unique name to every row/column in the matrix + + Returns + ------- + Scalar + """ + res = np.zeros(len(names), dtype=cls._use_dtype) + res['name'] = names + res['meta'] = [{} for _ in names] + return cls(res) + + def to_mapping(self, dim): + """ + Converts the hcp_labels to a MatrixIndicesMap for storage in CIFTI format + + Parameters + ---------- + dim : int + which dimension of the CIFTI vector/matrix is described by this dataset (zero-based) + + Returns + ------- + cifti2.Cifti2MatrixIndicesMap + """ + mim = cifti2.Cifti2MatrixIndicesMap([dim], 'CIFTI_INDEX_TYPE_SCALARS') + for elem in self.arr: + meta = None if len(elem['meta']) == 0 else elem['meta'] + named_map = cifti2.Cifti2NamedMap(elem['name'], cifti2.Cifti2MetaData(meta)) + mim.append(named_map) + return mim + + def get_element(self, index): + """ + Describes a single element from the axis + + Parameters + ---------- + index : int + Indexes the row/column of interest + + Returns + ------- + tuple with 2 elements + - unicode name of the get_scalar + - dictionary with the element metadata + """ + return self.arr['name'][index], self.arr['meta'][index] + + def to_label(self, labels): + """ + Creates a new Label axis based on the Scalar axis + + Parameters + ---------- + labels : list[dict] + mapping from integers to (name, (R, G, B, A)), where `name` is a string and R, G, B, and A are floats + between 0 and 1 giving the colour and alpha (transparency) + + Returns + ------- + Label + """ + res = np.zeros(self.size, dtype=Label._use_dtype) + res['name'] = self.arr['name'] + res['meta'] = self.arr['meta'] + res['get_label'] = labels + return Label(res) + + @property + def name(self, ): + return self.arr['name'] + + @name.setter + def name(self, values): + self.arr['name'] = values + + @property + def meta(self, ): + return self.arr['meta'] + + @meta.setter + def meta(self, values): + self.arr['meta'] = values + + +class Label(Axis): + """ + Along this axis of the CIFTI vector/matrix each row/column has been given a unique name, + get_label table, and optionally metadata + + Attributes + ---------- + name : np.ndarray + (N, ) string array with the parcel names + meta : np.ndarray + (N, ) array with a dictionary of metadata for each row/column + get_label : sp.ndarray + (N, ) array with dictionaries mapping integer values to get_label names and RGBA colors + """ + _use_dtype = np.dtype([('name', 'U60'), ('get_label', 'object'), ('meta', 'object')]) + + def __init__(self, arr): + """ + Creates a new Scalar axis from (name, meta-data) pairs + + Parameters + ---------- + arr : Iterable[Tuple[str, dict[int -> (str, (float, float, float, float)), dict(str->str)]] + iterates over all rows/columns assigning a name, dictionary mapping integers to get_label names and rgba colours + and a dictionary of metadata to each + """ + super(Label, self).__init__(arr) + + @classmethod + def from_mapping(cls, mim): + """ + Creates a new get_scalar axis based on a CIFTI dataset + + Parameters + ---------- + mim : cifti2.Cifti2MatrixIndicesMap + + Returns + ------- + Scalar + """ + tables = [{key: (value.label, value.rgba) for key, value in nm.label_table.items()} + for nm in mim.named_maps] + return Scalar.from_mapping(mim).to_label(tables) + + def to_mapping(self, dim): + """ + Converts the hcp_labels to a MatrixIndicesMap for storage in CIFTI format + + Parameters + ---------- + dim : int + which dimension of the CIFTI vector/matrix is described by this dataset (zero-based) + + Returns + ------- + cifti2.Cifti2MatrixIndicesMap + """ + mim = cifti2.Cifti2MatrixIndicesMap([dim], 'CIFTI_INDEX_TYPE_LABELS') + for elem in self.arr: + label_table = cifti2.Cifti2LabelTable() + for key, value in elem['get_label'].items(): + label_table[key] = (value[0],) + tuple(value[1]) + meta = None if len(elem['meta']) == 0 else elem['meta'] + named_map = cifti2.Cifti2NamedMap(elem['name'], cifti2.Cifti2MetaData(meta), + label_table) + mim.append(named_map) + return mim + + def get_element(self, index): + """ + Describes a single element from the axis + + Parameters + ---------- + index : int + Indexes the row/column of interest + + Returns + ------- + tuple with 2 elements + - unicode name of the get_scalar + - dictionary with the get_label table + - dictionary with the element metadata + """ + return self.arr['name'][index], self.arr['get_label'][index], self.arr['meta'][index] + + @property + def name(self, ): + return self.arr['name'] + + @name.setter + def name(self, values): + self.arr['name'] = values + + @property + def meta(self, ): + return self.arr['meta'] + + @meta.setter + def meta(self, values): + self.arr['meta'] = values + + @property + def label(self, ): + return self.arr['get_label'] + + @label.setter + def label(self, values): + self.arr['get_label'] = values + + +class Series(Axis): + """ + Along this axis of the CIFTI vector/matrix the rows/columns increase monotonously in time + + This Axis describes the time point of each row/column. + + Attributes + ---------- + start : float + starting time point + step : float + sampling time (TR) + size : int + number of time points + """ + size = None + _unit = None + + def __init__(self, start, step, size, unit="SECOND"): + """ + Creates a new Series axis + + Parameters + ---------- + start : float + Time of the first datapoint + step : float + Step size between data points + size : int + Number of data points + unit : str + Unit of the step size (one of 'second', 'hertz', 'meter', or 'radian') + """ + self.unit = unit + self.start = start + self.step = step + self.size = size + + @property + def unit(self, ): + return self._unit + + @unit.setter + def unit(self, value): + if value.upper() not in ("SECOND", "HERTZ", "METER", "RADIAN"): + raise ValueError("Series unit should be one of ('second', 'hertz', 'meter', or 'radian'") + self._unit = value.upper() + + @property + def arr(self, ): + return np.arange(self.size) * self.step + self.start + + @classmethod + def from_mapping(cls, mim): + """ + Creates a new Series axis based on a CIFTI dataset + + Parameters + ---------- + mim : cifti2.Cifti2MatrixIndicesMap + + Returns + ------- + Series + """ + start = mim.series_start * 10 ** mim.series_exponent + step = mim.series_step * 10 ** mim.series_exponent + return cls(start, step, mim.number_of_series_points, mim.series_unit) + + def to_mapping(self, dim): + """ + Converts the get_series to a MatrixIndicesMap for storage in CIFTI format + + Parameters + ---------- + dim : int + which dimension of the CIFTI vector/matrix is described by this dataset (zero-based) + + Returns + ------- + cifti2.Cifti2MatrixIndicesMap + """ + mim = cifti2.Cifti2MatrixIndicesMap([dim], 'CIFTI_INDEX_TYPE_SERIES') + mim.series_exponent = 0 + mim.series_start = self.start + mim.series_step = self.step + mim.number_of_series_points = self.size + mim.series_unit = self.unit + return mim + + def extend(self, other_axis): + """ + Concatenates two get_series + + Note: this will ignore the start point of the other axis + + Parameters + ---------- + other_axis : Series + other axis + + Returns + ------- + Series + """ + if other_axis.step != self.step: + raise ValueError('Can only concatenate get_series with the same step size') + if other_axis.unit != self.unit: + raise ValueError('Can only concatenate get_series with the same unit') + return Series(self.start, self.step, self.size + other_axis.size, self.unit) + + def __getitem__(self, item): + if isinstance(item, slice): + step = 1 if item.step is None else item.step + idx_start = ((self.size - 1 if step < 0 else 0) + if item.start is None else + (item.start if item.start >= 0 else self.size + item.start)) + idx_end = ((-1 if step < 0 else self.size) + if item.stop is None else + (item.stop if item.stop >= 0 else self.size + item.stop)) + if idx_start > self.size: + idx_start = self.size - 1 + if idx_end > self.size: + idx_end = self.size + nelements = (idx_end - idx_start) // step + if nelements < 0: + nelements = 0 + return Series(idx_start * self.step + self.start, self.step * step, nelements) + elif isinstance(item, int): + return self.get_element(item) + raise IndexError('Series can only be indexed with integers or slices without breaking the regular structure') + + def get_element(self, index): + """ + Gives the time point of a specific row/column + + Parameters + ---------- + index : int + Indexes the row/column of interest + + Returns + ------- + float + """ + if index < 0: + index = self.size + index + if index >= self.size: + raise IndexError("index %i is out of range for get_series with size %i" % (index, self.size)) + return self.start + self.step * index + + def __add__(self, other): + """ + Concatenates two Series + + Parameters + ---------- + other : Series + Time get_series to append at the end of the current time get_series. + Note that the starting time of the other time get_series is ignored. + + Returns + ------- + Series + New time get_series with the concatenation of the two + + Raises + ------ + ValueError + raised if the repetition time of the two time get_series is different + """ + if isinstance(other, Series): + return self.extend(other) + return NotImplemented diff --git a/nibabel/cifti2/tests/test_axes.py b/nibabel/cifti2/tests/test_axes.py new file mode 100644 index 0000000000..d7b24f03ec --- /dev/null +++ b/nibabel/cifti2/tests/test_axes.py @@ -0,0 +1,245 @@ +import numpy as np +from nose.tools import assert_raises +from .test_cifti2io_axes import check_rewrite +import nibabel.cifti2.cifti2_axes as axes + + +rand_affine = np.random.randn(4, 4) +vol_shape = (5, 10, 3) + + +def get_brain_models(): + """ + Generates a set of practice BrainModel axes + + Yields + ------ + BrainModel axis + """ + mask = np.zeros(vol_shape) + mask[0, 1, 2] = 1 + mask[0, 4, 2] = True + mask[0, 4, 0] = True + yield axes.BrainModel.from_mask(mask, 'ThalamusRight', rand_affine) + mask[0, 0, 0] = True + yield axes.BrainModel.from_mask(mask, affine=rand_affine) + + yield axes.BrainModel.from_surface([0, 5, 10], 15, 'CortexLeft') + yield axes.BrainModel.from_surface([0, 5, 10, 13], 15) + + surface_mask = np.zeros(15, dtype='bool') + surface_mask[[2, 9, 14]] = True + yield axes.BrainModel.from_mask(surface_mask, name='CortexRight') + + +def get_parcels(): + """ + Generates a practice Parcel axis out of all practice brain models + + Returns + ------- + Parcel axis + """ + bml = list(get_brain_models()) + return axes.Parcels.from_brain_models([('mixed', bml[0] + bml[2]), ('volume', bml[1]), ('surface', bml[3])]) + + +def get_scalar(): + """ + Generates a practice Scalar axis with names ('one', 'two', 'three') + + Returns + ------- + Scalar axis + """ + return axes.Scalar.from_names(['one', 'two', 'three']) + + +def get_label(): + """ + Generates a practice Label axis with names ('one', 'two', 'three') and two labels + + Returns + ------- + Label axis + """ + return axes.Scalar.from_names(['one', 'two', 'three']).to_label({0: ('something', (0.2, 0.4, 0.1, 0.5)), + 1: ('even better', (0.3, 0.8, 0.43, 0.9))}) + +def get_series(): + """ + Generates a set of 4 practice Series axes with different starting times/lengths/time steps and units + + Yields + ------ + Series axis + """ + yield axes.Series(3, 10, 4) + yield axes.Series(8, 10, 3) + yield axes.Series(3, 2, 4) + yield axes.Series(5, 10, 5, "HERTZ") + + +def get_axes(): + """ + Iterates through all of the practice axes defined in the functions above + + Yields + ------ + Cifti2 axis + """ + yield get_parcels() + yield get_scalar() + yield get_label() + for elem in get_brain_models(): + yield elem + for elem in get_series(): + yield elem + + +def test_brain_models(): + """ + Tests the introspection and creation of CIFTI2 BrainModel axes + """ + bml = list(get_brain_models()) + assert len(bml[0]) == 3 + assert (bml[0].vertex == -1).all() + assert (bml[0].voxel == [[0, 1, 2], [0, 4, 0], [0, 4, 2]]).all() + assert bml[0][1][0] == False + assert (bml[0][1][1] == [0, 4, 0]).all() + assert bml[0][1][2] == axes.BrainModel.to_cifti_brain_structure_name('thalamus_right') + assert len(bml[1]) == 4 + assert (bml[1].vertex == -1).all() + assert (bml[1].voxel == [[0, 0, 0], [0, 1, 2], [0, 4, 0], [0, 4, 2]]).all() + assert len(bml[2]) == 3 + assert (bml[2].voxel == -1).all() + assert (bml[2].vertex == [0, 5, 10]).all() + assert bml[2][1] == (True, 5, 'CIFTI_STRUCTURE_CORTEX_LEFT') + assert len(bml[3]) == 4 + assert (bml[3].voxel == -1).all() + assert (bml[3].vertex == [0, 5, 10, 13]).all() + assert bml[4][1] == (True, 9, 'CIFTI_STRUCTURE_CORTEX_RIGHT') + assert len(bml[4]) == 3 + assert (bml[4].voxel == -1).all() + assert (bml[4].vertex == [2, 9, 14]).all() + + for bm, label in zip(bml, ['ThalamusRight', 'Other', 'cortex_left', 'cortex']): + structures = list(bm.iter_structures()) + assert len(structures) == 1 + name = structures[0][0] + assert name == axes.BrainModel.to_cifti_brain_structure_name(label) + if 'CORTEX' in name: + assert bm.nvertices[name] == 15 + else: + assert name not in bm.nvertices + assert (bm.affine == rand_affine).all() + assert bm.volume_shape == vol_shape + + bmt = bml[0] + bml[1] + bml[2] + bml[3] + assert len(bmt) == 14 + structures = list(bmt.iter_structures()) + assert len(structures) == 4 + for bm, (name, _, bm_split) in zip(bml, structures): + assert bm == bm_split + assert (bm_split.name == name).all() + assert bm == bmt[bmt.name == bm.name[0]] + assert bm == bmt[np.where(bmt.name == bm.name[0])] + + bmt = bmt + bml[3] + assert len(bmt) == 18 + structures = list(bmt.iter_structures()) + assert len(structures) == 4 + assert len(structures[-1][2]) == 8 + + +def test_parcels(): + """ + Test the introspection and creation of CIFTI2 Parcel axes + """ + prc = get_parcels() + assert isinstance(prc, axes.Parcels) + assert prc['mixed'][0].shape == (3, 3) + assert len(prc['mixed'][1]) == 1 + assert prc['mixed'][1]['CIFTI_STRUCTURE_CORTEX_LEFT'].shape == (3, ) + + assert prc['volume'][0].shape == (4, 3) + assert len(prc['volume'][1]) == 0 + + assert prc['surface'][0].shape == (0, 3) + assert len(prc['surface'][1]) == 1 + assert prc['surface'][1]['CIFTI_STRUCTURE_CORTEX'].shape == (4, ) + + prc2 = prc + prc + assert len(prc2) == 6 + assert (prc2.affine == prc.affine).all() + assert (prc2.nvertices == prc.nvertices) + assert (prc2.volume_shape == prc.volume_shape) + assert prc2[:3] == prc + assert prc2[3:] == prc + + assert prc2[3:]['mixed'][0].shape == (3, 3) + assert len(prc2[3:]['mixed'][1]) == 1 + assert prc2[3:]['mixed'][1]['CIFTI_STRUCTURE_CORTEX_LEFT'].shape == (3, ) + + +def test_scalar(): + """ + Test the introspection and creation of CIFTI2 Scalar axes + """ + sc = get_scalar() + assert len(sc) == 3 + assert isinstance(sc, axes.Scalar) + assert (sc.name == ['one', 'two', 'three']).all() + assert sc[1] == ('two', {}) + sc2 = sc + sc + assert len(sc2) == 6 + assert (sc2.name == ['one', 'two', 'three', 'one', 'two', 'three']).all() + assert sc2[:3] == sc + assert sc2[3:] == sc + + +def test_series(): + """ + Test the introspection and creation of CIFTI2 Series axes + """ + sr = list(get_series()) + assert sr[0].unit == 'SECOND' + assert sr[1].unit == 'SECOND' + assert sr[2].unit == 'SECOND' + assert sr[3].unit == 'HERTZ' + + assert (sr[0].arr == np.arange(4) * 10 + 3).all() + assert (sr[1].arr == np.arange(3) * 10 + 8).all() + assert (sr[2].arr == np.arange(4) * 2 + 3).all() + assert ((sr[0] + sr[1]).arr == np.arange(7) * 10 + 3).all() + assert ((sr[1] + sr[0]).arr == np.arange(7) * 10 + 8).all() + assert ((sr[1] + sr[0] + sr[0]).arr == np.arange(11) * 10 + 8).all() + assert sr[1][2] == 28 + assert sr[1][-2] == sr[1].arr[-2] + assert_raises(ValueError, lambda: sr[0] + sr[2]) + assert_raises(ValueError, lambda: sr[2] + sr[1]) + assert_raises(ValueError, lambda: sr[0] + sr[3]) + assert_raises(ValueError, lambda: sr[3] + sr[1]) + assert_raises(ValueError, lambda: sr[3] + sr[2]) + + # test slicing + assert (sr[0][1:3].arr == sr[0].arr[1:3]).all() + assert (sr[0][1:].arr == sr[0].arr[1:]).all() + assert (sr[0][:-2].arr == sr[0].arr[:-2]).all() + assert (sr[0][1:-1].arr == sr[0].arr[1:-1]).all() + assert (sr[0][1:-1:2].arr == sr[0].arr[1:-1:2]).all() + assert (sr[0][::2].arr == sr[0].arr[::2]).all() + assert (sr[0][:10:2].arr == sr[0].arr[::2]).all() + assert (sr[0][10::-1].arr == sr[0].arr[::-1]).all() + assert (sr[0][3:1:-1].arr == sr[0].arr[3:1:-1]).all() + assert (sr[0][1:3:-1].arr == sr[0].arr[1:3:-1]).all() + + +def test_writing(): + """ + Tests the writing and reading back in of custom created CIFTI2 axes + """ + for ax1 in get_axes(): + for ax2 in get_axes(): + arr = np.random.randn(len(ax1), len(ax2)) + check_rewrite(arr, (ax1, ax2)) diff --git a/nibabel/cifti2/tests/test_cifti2io_axes.py b/nibabel/cifti2/tests/test_cifti2io_axes.py new file mode 100644 index 0000000000..22f4c27253 --- /dev/null +++ b/nibabel/cifti2/tests/test_cifti2io_axes.py @@ -0,0 +1,176 @@ +from nibabel.cifti2 import cifti2_axes, cifti2 +from nibabel.tests.nibabel_data import get_nibabel_data, needs_nibabel_data +import nibabel as nib +import os +import numpy as np +import tempfile + +test_directory = os.path.join(get_nibabel_data(), 'nitest-cifti2') + +hcp_labels = ['CortexLeft', 'CortexRight', 'AccumbensLeft', 'AccumbensRight', 'AmygdalaLeft', 'AmygdalaRight', + 'brain_stem', 'CaudateLeft', 'CaudateRight', 'CerebellumLeft', 'CerebellumRight', + 'Diencephalon_ventral_left', 'Diencephalon_ventral_right', 'HippocampusLeft', 'HippocampusRight', + 'PallidumLeft', 'PallidumRight', 'PutamenLeft', 'PutamenRight', 'ThalamusLeft', 'ThalamusRight'] + +hcp_n_elements = [29696, 29716, 135, 140, 315, 332, 3472, 728, 755, 8709, 9144, 706, + 712, 764, 795, 297, 260, 1060, 1010, 1288, 1248] + +hcp_affine = np.array([[ -2., 0., 0., 90.], + [ 0., 2., 0., -126.], + [ 0., 0., 2., -72.], + [ 0., 0., 0., 1.]]) + + +def check_hcp_grayordinates(brain_model): + """Checks that a BrainModel matches the expected 32k HCP grayordinates + """ + assert isinstance(brain_model, cifti2_axes.BrainModel) + structures = list(brain_model.iter_structures()) + assert len(structures) == len(hcp_labels) + idx_start = 0 + for idx, (name, _, bm), label, nel in zip(range(len(structures)), structures, hcp_labels, hcp_n_elements): + if idx < 2: + assert name in bm.nvertices.keys() + assert (bm.voxel == -1).all() + assert (bm.vertex != -1).any() + assert bm.nvertices[name] == 32492 + else: + assert name not in bm.nvertices.keys() + assert (bm.voxel != -1).any() + assert (bm.vertex == -1).all() + assert (bm.affine == hcp_affine).all() + assert bm.volume_shape == (91, 109, 91) + assert name == cifti2_axes.BrainModel.to_cifti_brain_structure_name(label) + assert len(bm) == nel + assert (bm.arr == brain_model.arr[idx_start:idx_start + nel]).all() + idx_start += nel + assert idx_start == len(brain_model) + + assert (brain_model.arr[:5]['vertex'] == np.arange(5)).all() + assert structures[0][2].vertex[-1] == 32491 + assert structures[1][2].vertex[0] == 0 + assert structures[1][2].vertex[-1] == 32491 + assert (structures[-1][2].arr[-1] == brain_model.arr[-1]).all() + assert (brain_model.arr[-1]['voxel'] == [38, 55, 46]).all() + assert (brain_model.arr[70000]['voxel'] == [56, 22, 19]).all() + + +def check_Conte69(brain_model): + """Checks that the BrainModel matches the expected Conte69 surface coordinates + """ + assert isinstance(brain_model, cifti2_axes.BrainModel) + structures = list(brain_model.iter_structures()) + assert len(structures) == 2 + assert structures[0][0] == 'CIFTI_STRUCTURE_CORTEX_LEFT' + assert structures[0][2].is_surface.all() + assert structures[1][0] == 'CIFTI_STRUCTURE_CORTEX_RIGHT' + assert structures[1][2].is_surface.all() + assert (brain_model.voxel == -1).all() + + assert (brain_model.arr[:5]['vertex'] == np.arange(5)).all() + assert structures[0][2].vertex[-1] == 32491 + assert structures[1][2].vertex[0] == 0 + assert structures[1][2].vertex[-1] == 32491 + + +def check_rewrite(arr, axes, extension='.nii'): + """ + Checks wheter writing the Cifti2 array to disc and reading it back in gives the same object + + Parameters + ---------- + arr : array + N-dimensional array of data + axes : Sequence[cifti2_axes.Axis] + sequence of length N with the meaning of the rows/columns along each dimension + extension : str + custom extension to use + """ + (fd, name) = tempfile.mkstemp(extension) + cifti2.Cifti2Image(arr, header=axes).to_filename(name) + img = nib.load(name) + arr2 = img.get_data() + assert (arr == arr2).all() + for idx in range(len(img.shape)): + assert (axes[idx] == img.header.get_axis(idx)) + return img + + +@needs_nibabel_data('nitest-cifti2') +def test_read_ones(): + img = nib.load(os.path.join(test_directory, 'ones.dscalar.nii')) + arr = img.get_data() + axes = [img.header.get_axis(dim) for dim in range(2)] + assert (arr == 1).all() + assert isinstance(axes[0], cifti2_axes.Scalar) + assert len(axes[0]) == 1 + assert axes[0].name[0] == 'ones' + assert axes[0].meta[0] == {} + check_hcp_grayordinates(axes[1]) + img = check_rewrite(arr, axes) + check_hcp_grayordinates(img.header.get_axis(1)) + + +@needs_nibabel_data('nitest-cifti2') +def test_read_conte69_dscalar(): + img = nib.load(os.path.join(test_directory, 'Conte69.MyelinAndCorrThickness.32k_fs_LR.dscalar.nii')) + arr = img.get_data() + axes = [img.header.get_axis(dim) for dim in range(2)] + assert isinstance(axes[0], cifti2_axes.Scalar) + assert len(axes[0]) == 2 + assert axes[0].name[0] == 'MyelinMap_BC_decurv' + assert axes[0].name[1] == 'corrThickness' + assert axes[0].meta[0] == {'PaletteColorMapping': '\n MODE_AUTO_SCALE_PERCENTAGE\n 98.000000 2.000000 2.000000 98.000000\n -100.000000 0.000000 0.000000 100.000000\n ROY-BIG-BL\n true\n true\n false\n true\n THRESHOLD_TEST_SHOW_OUTSIDE\n THRESHOLD_TYPE_OFF\n false\n -1.000000 1.000000\n -1.000000 1.000000\n -1.000000 1.000000\n \n PALETTE_THRESHOLD_RANGE_MODE_MAP\n'} + check_Conte69(axes[1]) + check_rewrite(arr, axes) + + +@needs_nibabel_data('nitest-cifti2') +def test_read_conte69_dtseries(): + img = nib.load(os.path.join(test_directory, 'Conte69.MyelinAndCorrThickness.32k_fs_LR.dtseries.nii')) + arr = img.get_data() + axes = [img.header.get_axis(dim) for dim in range(2)] + assert isinstance(axes[0], cifti2_axes.Series) + assert len(axes[0]) == 2 + assert axes[0].start == 0 + assert axes[0].step == 1 + assert axes[0].size == arr.shape[0] + assert (axes[0].arr == [0, 1]).all() + check_Conte69(axes[1]) + check_rewrite(arr, axes) + + +@needs_nibabel_data('nitest-cifti2') +def test_read_conte69_dlabel(): + img = nib.load(os.path.join(test_directory, 'Conte69.parcellations_VGD11b.32k_fs_LR.dlabel.nii')) + arr = img.get_data() + axes = [img.header.get_axis(dim) for dim in range(2)] + assert isinstance(axes[0], cifti2_axes.Label) + assert len(axes[0]) == 3 + assert (axes[0].name == ['Composite Parcellation-lh (FRB08_OFP03_retinotopic)', + 'Brodmann lh (from colin.R via pals_R-to-fs_LR)', 'MEDIAL WALL lh (fs_LR)']).all() + assert axes[0].label[1][70] == ('19_B05', (1.0, 0.867, 0.467, 1.0)) + assert (axes[0].meta == [{}] * 3).all() + check_Conte69(axes[1]) + check_rewrite(arr, axes) + + +@needs_nibabel_data('nitest-cifti2') +def test_read_conte69_ptseries(): + img = nib.load(os.path.join(test_directory, 'Conte69.MyelinAndCorrThickness.32k_fs_LR.ptseries.nii')) + arr = img.get_data() + axes = [img.header.get_axis(dim) for dim in range(2)] + assert isinstance(axes[0], cifti2_axes.Series) + assert len(axes[0]) == 2 + assert axes[0].start == 0 + assert axes[0].step == 1 + assert axes[0].size == arr.shape[0] + assert (axes[0].arr == [0, 1]).all() + + assert len(axes[1]) == 54 + voxels, vertices = axes[1]['ER_FRB08'] + assert voxels.shape == (0, 3) + assert len(vertices) == 2 + assert vertices['CIFTI_STRUCTURE_CORTEX_LEFT'].shape == (206 // 2, ) + assert vertices['CIFTI_STRUCTURE_CORTEX_RIGHT'].shape == (206 // 2, ) + check_rewrite(arr, axes) diff --git a/nibabel/cifti2/tests/test_cifti2io.py b/nibabel/cifti2/tests/test_cifti2io_header.py similarity index 100% rename from nibabel/cifti2/tests/test_cifti2io.py rename to nibabel/cifti2/tests/test_cifti2io_header.py diff --git a/nibabel/cifti2/tests/test_name.py b/nibabel/cifti2/tests/test_name.py new file mode 100644 index 0000000000..a73c5e8c46 --- /dev/null +++ b/nibabel/cifti2/tests/test_name.py @@ -0,0 +1,19 @@ +from nibabel.cifti2 import cifti2_axes + +equivalents = [('CIFTI_STRUCTURE_CORTEX_LEFT', ('CortexLeft', 'LeftCortex', 'left_cortex', 'Left Cortex', + 'Cortex_Left', 'cortex left', 'CORTEX_LEFT', 'LEFT CORTEX', + ('cortex', 'left'), ('CORTEX', 'Left'), ('LEFT', 'coRTEX'))), + ('CIFTI_STRUCTURE_CORTEX', ('Cortex', 'CortexBOTH', 'Cortex_both', 'both cortex', + 'BOTH_CORTEX', 'cortex', 'CORTEX', ('cortex', ), + ('COrtex', 'Both'), ('both', 'cortex')))] + + +def test_name_conversion(): + """ + Tests the automatic name conversion to a format recognized by CIFTI2 + """ + func = cifti2_axes.BrainModel.to_cifti_brain_structure_name + for base_name, input_names in equivalents: + assert base_name == func(base_name) + for name in input_names: + assert base_name == func(name) \ No newline at end of file