# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## # # See COPYING file distributed along with the NiBabel package for the # copyright and license terms. # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Read images in PAR/REC format This is yet another MRI image format generated by Philips scanners. It is an ASCII header (PAR) plus a binary blob (REC). This implementation aims to read version 4.0 through 4.2 of this format. Other versions could probably be supported, but we need example images to test against. If you want us to support another version, and have an image we can add to the test suite, let us know. You would make us very happy by submitting a pull request. ############### PAR file format ############### The PAR format appears to have two sections: General information ################### This is a set of lines each giving one key : value pair, examples:: . EPI factor <0,1=no EPI> : 39 . Dynamic scan <0=no 1=yes> ? : 1 . Diffusion <0=no 1=yes> ? : 0 (from ``nibabel/tests/data/phantom_EPI_asc_CLEAR_2_1.PAR``) Image information ################# There is a ``#`` prefixed list of fields under the heading "IMAGE INFORMATION DEFINITION". From the same file, here is the start of this list:: # === IMAGE INFORMATION DEFINITION ============================================= # The rest of this file contains ONE line per image, this line contains the following information: # # slice number (integer) # echo number (integer) # dynamic scan number (integer) There follows a space separated table with values for these fields, each row containing all the named values. Here are the first few lines from the example file above:: # === IMAGE INFORMATION ========================================================== # sl ec dyn ph ty idx pix scan% rec size (re)scale window angulation offcentre thick gap info spacing echo dtime ttime diff avg flip freq RR-int turbo delay b grad cont anis diffusion L.ty 1 1 1 1 0 2 0 16 62 64 64 0.00000 1.29035 4.28404e-003 1070 1860 -13.26 -0.00 -0.00 2.51 -0.81 -8.69 6.000 2.000 0 1 0 2 3.750 3.750 30.00 0.00 0.00 0.00 0 90.00 0 0 0 39 0.0 1 1 8 0 0.000 0.000 0.000 1 2 1 1 1 0 2 1 16 62 64 64 0.00000 1.29035 4.28404e-003 1122 1951 -13.26 -0.00 -0.00 2.51 6.98 -10.53 6.000 2.000 0 1 0 2 3.750 3.750 30.00 0.00 0.00 0.00 0 90.00 0 0 0 39 0.0 1 1 8 0 0.000 0.000 0.000 1 3 1 1 1 0 2 2 16 62 64 64 0.00000 1.29035 4.28404e-003 1137 1977 -13.26 -0.00 -0.00 2.51 14.77 -12.36 6.000 2.000 0 1 0 2 3.750 3.750 30.00 0.00 0.00 0.00 0 90.00 0 0 0 39 0.0 1 1 8 0 0.000 0.000 0.000 1 Orientation ########### PAR files refer to orientations "ap", "fh" and "rl". Nibabel's required affine output axes are RAS (left to Right, posterior to Anterior, inferior to Superior). The correspondence of the PAR file's axes to RAS axes is: * ap = anterior -> posterior = negative A in RAS = P * fh = foot -> head = S in RAS = S * rl = right -> left = negative R in RAS = L We therefore call the PAR file's axis system "PSL" (Posterior, Superior, Left). The orientation of the PAR file axes corresponds to DICOM's LPS coordinate system (right to Left, anterior to Posterior, inferior to Superior), but in a different order. Data type ######### It seems that everyone agrees that Philips stores REC data in little-endian format - see https://github.com/nipy/nibabel/issues/274 Philips XML header files, and some previous experience, suggest that the REC data is always stored as 8 or 16 bit unsigned integers - see https://github.com/nipy/nibabel/issues/275 Data Sorting ############ PAR/REC files have a large number of potential image dimensions. To handle sorting of volumes in PAR/REC files based on these fields and not the order slices first appear in the PAR file, the ``strict_sort`` flag of ``nibabel.load`` (or ``parrec.load``) should be set to ``True``. The fields that are taken into account during sorting are: - slice number - echo number - cardiac phase number - gradient orientation number - diffusion b value number - label type (ASL tag vs. control) - dynamic scan number - image_type_mr (Re, Im, Mag, Phase) Slices are sorted into the third dimension and the order of preference for sorting along the 4th dimension corresponds to the order in the list above. If the image data has more than 4 dimensions these will all be concatenated along the 4th dimension. For example, for a scan with two echos and two dynamics, the 4th dimension will have both echos of dynamic 1 prior to the two echos for dynamic 2. The``get_volume_labels`` method of the header returns a dictionary containing the PAR field labels for this 4th dimension. The volume sorting described above can be enabled in the parrec2nii command utility via the option "--strict-sort". The dimension info can be exported to a CSV file by adding the option "--volume-info". """ import re import warnings from collections import OrderedDict from copy import deepcopy from io import StringIO from locale import getpreferredencoding import numpy as np from .affines import apply_affine, dot_reduce, from_matvec from .eulerangles import euler2mat from .fileslice import fileslice, strided_scalar from .nifti1 import unit_codes from .openers import ImageOpener from .spatialimages import SpatialHeader, SpatialImage from .volumeutils import Recoder, array_from_file # PSL to RAS affine PSL_TO_RAS = np.array( [ [0, 0, -1, 0], # L -> R [-1, 0, 0, 0], # P -> A [0, 1, 0, 0], # S -> S [0, 0, 0, 1], ] ) # Acquisition (tra/sag/cor) to PSL axes # These come from looking at transverse, sagittal, coronal datasets where we # can see the LR, PA, SI orientation of the slice axes from the scanned object ACQ_TO_PSL = dict( transverse=np.array( [ [0, 1, 0, 0], # P [0, 0, 1, 0], # S [1, 0, 0, 0], # L [0, 0, 0, 1], ] ), sagittal=np.diag([1, -1, -1, 1]), coronal=np.array( [ [0, 0, 1, 0], # P [0, -1, 0, 0], # S [1, 0, 0, 0], # L [0, 0, 0, 1], ] ), ) DEG2RAD = np.pi / 180.0 # General information dict definitions # assign props to PAR header entries # values are: (shortname[, dtype[, shape]]) # if shape is None, the number of elements is to be determined on read _hdr_key_dict = { 'Patient name': ('patient_name',), 'Examination name': ('exam_name',), 'Protocol name': ('protocol_name',), 'Examination date/time': ('exam_date',), 'Series Type': ('series_type',), 'Acquisition nr': ('acq_nr', int), 'Reconstruction nr': ('recon_nr', int), 'Scan Duration [sec]': ('scan_duration', float), 'Max. number of cardiac phases': ('max_cardiac_phases', int), 'Max. number of echoes': ('max_echoes', int), 'Max. number of slices/locations': ('max_slices', int), 'Max. number of dynamics': ('max_dynamics', int), 'Max. number of mixes': ('max_mixes', int), 'Patient position': ('patient_position',), 'Preparation direction': ('prep_direction',), 'Technique': ('tech',), 'Scan resolution (x, y)': ('scan_resolution', int, (2,)), 'Scan mode': ('scan_mode',), 'Repetition time [ms]': ('repetition_time', float, None), 'FOV (ap,fh,rl) [mm]': ('fov', float, (3,)), 'Water Fat shift [pixels]': ('water_fat_shift', float), 'Angulation midslice(ap,fh,rl)[degr]': ('angulation', float, (3,)), 'Off Centre midslice(ap,fh,rl) [mm]': ('off_center', float, (3,)), 'Flow compensation <0=no 1=yes> ?': ('flow_compensation', int), 'Presaturation <0=no 1=yes> ?': ('presaturation', int), 'Phase encoding velocity [cm/sec]': ('phase_enc_velocity', float, (3,)), 'MTC <0=no 1=yes> ?': ('mtc', int), 'SPIR <0=no 1=yes> ?': ('spir', int), 'EPI factor <0,1=no EPI>': ('epi_factor', int), 'Dynamic scan <0=no 1=yes> ?': ('dyn_scan', int), 'Diffusion <0=no 1=yes> ?': ('diffusion', int), 'Diffusion echo time [ms]': ('diffusion_echo_time', float), # Lines below added for par / rec versions > 4 'Max. number of diffusion values': ('max_diffusion_values', int), 'Max. number of gradient orients': ('max_gradient_orient', int), # Line below added for par / rec version > 4.1 'Number of label types <0=no ASL>': ('nr_label_types', int), # The following are duplicates of the above fields, but with slightly # different abbreviation, spelling, or capatilization. Both variants have # been observed in the wild in V4.2 PAR files: # https://github.com/nipy/nibabel/issues/505 'Series_data_type': ('series_type',), 'Patient Position': ('patient_position',), 'Repetition time [msec]': ('repetition_time', float, None), 'Diffusion echo time [msec]': ('diffusion_echo_time', float), } # Image information as coded into a numpy structured array # header items order per image definition line image_def_dtds = {} image_def_dtds['V4'] = [ ('slice number', int), ('echo number', int), ('dynamic scan number', int), ('cardiac phase number', int), ('image_type_mr', int), ('scanning sequence', int), ('index in REC file', int), ('image pixel size', int), ('scan percentage', int), ('recon resolution', int, (2)), ('rescale intercept', float), ('rescale slope', float), ('scale slope', float), # Window center, width recorded as integer but can be float ('window center', float), ('window width', float), ('image angulation', float, (3,)), ('image offcentre', float, (3,)), ('slice thickness', float), ('slice gap', float), ('image_display_orientation', int), ('slice orientation', int), ('fmri_status_indication', int), ('image_type_ed_es', int), ('pixel spacing', float, (2,)), ('echo_time', float), ('dyn_scan_begin_time', float), ('trigger_time', float), ('diffusion_b_factor', float), ('number of averages', int), ('image_flip_angle', float), ('cardiac frequency', int), ('minimum RR-interval', int), ('maximum RR-interval', int), ('TURBO factor', int), ('Inversion delay', float), ] # Extra image def fields for 4.1 compared to 4 # fmt: off image_def_dtds['V4.1'] = image_def_dtds['V4'] + [ ('diffusion b value number', int), # (imagekey!) ('gradient orientation number', int), # (imagekey!) ('contrast type', 'S30'), # XXX might be too short? ('diffusion anisotropy type', 'S30'), # XXX might be too short? ('diffusion', float, (3,)), ] # Extra image def fields for 4.2 compared to 4.1 image_def_dtds['V4.2'] = image_def_dtds['V4.1'] + [ ('label type', int), # (imagekey!) ] # fmt: on #: PAR header versions we claim to understand supported_versions = list(image_def_dtds.keys()) #: Deprecated; please don't use image_def_dtype = np.dtype(image_def_dtds['V4.2']) #: slice orientation codes slice_orientation_codes = Recoder( ( # code, label (1, 'transverse'), (2, 'sagittal'), (3, 'coronal'), ), fields=('code', 'label'), ) class PARRECError(Exception): """Exception for PAR/REC format related problems. To be raised whenever PAR/REC is not happy, or we are not happy with PAR/REC. """ # Value after colon may be absent GEN_RE = re.compile(r'.\s+(.*?)\s*:\s*(.*)') def _split_header(fobj): """Split header into `version`, `gen_dict`, `image_lines`""" version = None gen_dict = {} image_lines = [] # Small state-machine state = 'top-header' for line in fobj: line = line.strip() if line == '': continue if state == 'top-header': if not line.startswith('#'): state = 'general-info' elif 'image export tool' in line: version = line.split()[-1] if state == 'general-info': if not line.startswith('.'): state = 'comment-block' else: # Let match raise error for unexpected field format key, value = GEN_RE.match(line).groups() gen_dict[key] = value if state == 'comment-block': if not line.startswith('#'): state = 'image-info' if state == 'image-info': if line.startswith('#'): break image_lines.append(line) return version, gen_dict, image_lines def _process_gen_dict(gen_dict): """Process `gen_dict` key, values into `general_info`""" general_info = {} hdr_key_dict_lower = {key.lower(): val for key, val in _hdr_key_dict.items()} for key, value in gen_dict.items(): # get props for this hdr field props = hdr_key_dict_lower[key.lower()] # turn values into meaningful dtype if len(props) == 2: # only dtype spec and no shape value = props[1](value) elif len(props) == 3: # array with dtype and shape value = np.fromstring(value, props[1], sep=' ') # if shape is None, allow arbitrary length if props[2] is not None: value = value.reshape(props[2]) general_info[props[0]] = value return general_info def _process_image_lines(image_lines, version): """Process image information definition lines according to `version`""" # postproc image def props image_def_dtd = image_def_dtds[version] # create an array for all image defs image_defs = np.zeros(len(image_lines), dtype=image_def_dtd) # for every image definition for i, line in enumerate(image_lines): items = line.split() item_counter = 0 # for all image properties we know about for props in image_def_dtd: if len(props) == 2: name, np_type = props value = items[item_counter] if not np.dtype(np_type).kind == 'S': value = np_type(value) item_counter += 1 elif len(props) == 3: name, np_type, shape = props nelements = np.prod(shape) value = items[item_counter : item_counter + nelements] value = [np_type(v) for v in value] item_counter += nelements image_defs[name][i] = value return image_defs def vol_numbers(slice_nos): """Calculate volume numbers inferred from slice numbers `slice_nos` The volume number for each slice is the number of times this slice number has occurred previously in the `slice_nos` sequence Parameters ---------- slice_nos : sequence Sequence of slice numbers, e.g. ``[1, 2, 3, 4, 1, 2, 3, 4]``. Returns ------- vol_nos : list A list, the same length of `slice_nos` giving the volume number for each corresponding slice number. """ counter = {} vol_nos = [] for s_no in slice_nos: count = counter.setdefault(s_no, 0) vol_nos.append(count) counter[s_no] += 1 return vol_nos def vol_is_full(slice_nos, slice_max, slice_min=1): """Vector with True for slices in complete volume, False otherwise Parameters ---------- slice_nos : sequence Sequence of slice numbers, e.g. ``[1, 2, 3, 4, 1, 2, 3, 4]``. slice_max : int Highest slice number for a full slice set. Slice set will be ``range(slice_min, slice_max+1)``. slice_min : int, optional Lowest slice number for full slice set. Default is 1. Returns ------- is_full : array Bool vector with True for slices in full volumes, False for slices in partial volumes. A full volume is a volume with all slices in the ``slice set`` as defined above. Raises ------ ValueError if any value in `slice_nos` is outside slice set indices. """ slice_set = set(range(slice_min, slice_max + 1)) if not slice_set.issuperset(slice_nos): raise ValueError(f'Slice numbers outside inclusive range {slice_min} to {slice_max}') vol_nos = np.array(vol_numbers(slice_nos)) slice_nos = np.asarray(slice_nos) is_full = np.ones(slice_nos.shape, dtype=bool) for vol_no in set(vol_nos): ours = vol_nos == vol_no if not set(slice_nos[ours]) == slice_set: is_full[ours] = False return is_full def _truncation_checks(general_info, image_defs, permit_truncated): """Check for presence of truncation in PAR file parameters Raise error if truncation present and `permit_truncated` is False. """ def _err_or_warn(msg): if not permit_truncated: raise PARRECError(msg) warnings.warn(msg) def _chk_trunc(idef_name, gdef_max_name): if gdef_max_name not in general_info: return id_values = image_defs[idef_name + ' number'] n_have = len(set(id_values)) n_expected = general_info[gdef_max_name] if n_have != n_expected: _err_or_warn( f'Header inconsistency: Found {n_have} {idef_name} ' f'values, but expected {n_expected}' ) _chk_trunc('slice', 'max_slices') _chk_trunc('echo', 'max_echoes') _chk_trunc('dynamic scan', 'max_dynamics') _chk_trunc('diffusion b value', 'max_diffusion_values') _chk_trunc('gradient orientation', 'max_gradient_orient') # Final check for partial volumes if not np.all(vol_is_full(image_defs['slice number'], general_info['max_slices'])): _err_or_warn('Found one or more partial volume(s)') def one_line(long_str): """Make maybe mutli-line `long_str` into one long line""" return ' '.join(line.strip() for line in long_str.splitlines()) def parse_PAR_header(fobj): """Parse a PAR header and aggregate all information into useful containers. Parameters ---------- fobj : file-object The PAR header file object. Returns ------- general_info : dict Contains all "General Information" from the header file image_info : ndarray Structured array with fields giving all "Image information" in the header """ # single pass through the header version, gen_dict, image_lines = _split_header(fobj) if version not in supported_versions: warnings.warn( one_line( f""" PAR/REC version '{version}' is currently not supported -- making an attempt to read nevertheless. Please email the NiBabel mailing list, if you are interested in adding support for this version. """ ) ) general_info = _process_gen_dict(gen_dict) image_defs = _process_image_lines(image_lines, version) return general_info, image_defs def _data_from_rec( rec_fileobj, in_shape, dtype, slice_indices, out_shape, scalings=None, mmap=True ): """Load and return array data from REC file Parameters ---------- rec_fileobj : file-like The file to process. in_shape : tuple The input shape inferred from the PAR file. dtype : dtype The datatype. slice_indices : array of int The indices used to re-index the resulting array properly. out_shape : tuple The output shape. scalings : {None, sequence}, optional Scalings to use. If not None, a length 2 sequence giving (``slope``, ``intercept``), where ``slope`` and ``intercept`` are arrays that can be broadcast to `out_shape`. mmap : {True, False, 'c', 'r', 'r+'}, optional `mmap` controls the use of numpy memory mapping for reading data. If False, do not try numpy ``memmap`` for data array. If one of {'c', 'r', 'r+'}, try numpy memmap with ``mode=mmap``. A `mmap` value of True gives the same behavior as ``mmap='c'``. If `rec_fileobj` cannot be memory-mapped, ignore `mmap` value and read array from file. Returns ------- data : array The scaled and sorted array. """ rec_data = array_from_file(in_shape, dtype, rec_fileobj, mmap=mmap) rec_data = rec_data[..., slice_indices] rec_data = rec_data.reshape(out_shape, order='F') if scalings is not None: # Don't do in-place b/c this goes int16 -> float64 rec_data = rec_data * scalings[0] rec_data += scalings[1] return rec_data def exts2pars(exts_source): """Parse, return any PAR headers from NIfTI extensions in `exts_source` Parameters ---------- exts_source : sequence or `Nifti1Image`, `Nifti1Header` instance A sequence of extensions, or header containing NIfTI extensions, or an image containing a header with NIfTI extensions. Returns ------- par_headers : list A list of PARRECHeader objects, usually empty or with one element, each element contains a PARRECHeader read from the contained extensions. """ headers = [] exts_source = exts_source.header if hasattr(exts_source, 'header') else exts_source exts_source = exts_source.extensions if hasattr(exts_source, 'extensions') else exts_source for extension in exts_source: content = extension.get_content() content = content.decode(getpreferredencoding(False)) if not content.startswith('# === DATA DESCRIPTION FILE ==='): continue gen_info, image_info = parse_PAR_header(StringIO(content)) headers.append(PARRECHeader(gen_info, image_info)) return headers class PARRECArrayProxy: def __init__(self, file_like, header, *, mmap=True, scaling='dv'): """Initialize PARREC array proxy Parameters ---------- file_like : file-like object Filename or object implementing ``read, seek, tell`` header : PARRECHeader instance Implementing ``get_data_shape, get_data_dtype``, ``get_sorted_slice_indices``, ``get_data_scaling``, ``get_rec_shape``. mmap : {True, False, 'c', 'r'}, optional, keyword only `mmap` controls the use of numpy memory mapping for reading data. If False, do not try numpy ``memmap`` for data array. If one of {'c', 'r'}, try numpy memmap with ``mode=mmap``. A `mmap` value of True gives the same behavior as ``mmap='c'``. If `file_like` cannot be memory-mapped, ignore `mmap` value and read array from file. scaling : {'fp', 'dv'}, optional, keyword only Type of scaling to use - see header ``get_data_scaling`` method. """ if mmap not in (True, False, 'c', 'r'): raise ValueError("mmap should be one of {True, False, 'c', 'r'}") self.file_like = file_like # Copies of values needed to read array self._shape = header.get_data_shape() self._dtype = header.get_data_dtype() self._slice_indices = header.get_sorted_slice_indices() self._mmap = mmap self._slice_scaling = header.get_data_scaling(scaling) self._rec_shape = header.get_rec_shape() @property def shape(self): return self._shape @property def ndim(self): return len(self.shape) @property def dtype(self): return self._dtype @property def is_proxy(self): return True def _get_unscaled(self, slicer): indices = self._slice_indices if slicer == (): with ImageOpener(self.file_like) as fileobj: rec_data = array_from_file(self._rec_shape, self._dtype, fileobj, mmap=self._mmap) rec_data = rec_data[..., indices] return rec_data.reshape(self._shape, order='F') elif indices[0] != 0 or np.any(np.diff(indices) != 1): # We can't load direct from REC file, use inefficient slicing return self._get_unscaled(())[slicer] # Slices all sequential from zero, can use fileslice # This gives more efficient volume by volume loading, for example with ImageOpener(self.file_like) as fileobj: return fileslice(fileobj, slicer, self._shape, self._dtype, 0, 'F') def _get_scaled(self, dtype, slicer): raw_data = self._get_unscaled(slicer) if self._slice_scaling is None: if dtype is None: return raw_data final_type = np.promote_types(raw_data.dtype, dtype) return raw_data.astype(final_type, copy=False) # Broadcast scaling to shape of original data fake_data = strided_scalar(self._shape) _, slopes, inters = np.broadcast_arrays(fake_data, *self._slice_scaling) final_type = np.result_type(raw_data, slopes, inters) if dtype is not None: final_type = np.promote_types(final_type, dtype) # Slice scaling to give output shape return raw_data * slopes[slicer].astype(final_type) + inters[slicer].astype(final_type) def get_unscaled(self): """Read data from file This is an optional part of the proxy API """ return self._get_unscaled(slicer=()) def __array__(self, dtype=None): """Read data from file and apply scaling, casting to ``dtype`` If ``dtype`` is unspecified, the dtype of the returned array is the narrowest dtype that can represent the data without overflow. Generally, it is the wider of the dtypes of the slopes or intercepts. Parameters ---------- dtype : numpy dtype specifier, optional A numpy dtype specifier specifying the type of the returned array. Returns ------- array Scaled image data with type `dtype`. """ arr = self._get_scaled(dtype=dtype, slicer=()) if dtype is not None: arr = arr.astype(dtype, copy=False) return arr def __getitem__(self, slicer): return self._get_scaled(dtype=None, slicer=slicer) class PARRECHeader(SpatialHeader): """PAR/REC header""" def __init__(self, info, image_defs, permit_truncated=False, strict_sort=False): """ Parameters ---------- info : dict "General information" from the PAR file (as returned by `parse_PAR_header()`). image_defs : array Structured array with image definitions from the PAR file (as returned by `parse_PAR_header()`). permit_truncated : bool, optional If True, a warning is emitted instead of an error when a truncated recording is detected. strict_sort : bool, optional, keyword-only If True, a larger number of header fields are used while sorting the REC data array. This may produce a different sort order than `strict_sort=False`, where volumes are sorted by the order in which the slices appear in the .PAR file. """ self.general_info = info.copy() self.image_defs = image_defs.copy() self.permit_truncated = permit_truncated self.strict_sort = strict_sort _truncation_checks(info, image_defs, permit_truncated) # charge with basic properties to be able to use base class # functionality # dtype bitpix = self._get_unique_image_prop('image pixel size') if bitpix not in (8, 16): raise PARRECError( f'Only 8- and 16-bit data supported (not {bitpix}) ' 'please report this to the nibabel developers' ) # REC data always little endian dt = np.dtype('uint' + str(bitpix)).newbyteorder('<') super().__init__(data_dtype=dt, shape=self._calc_data_shape(), zooms=self._calc_zooms()) @classmethod def from_header(klass, header=None): if header is None: raise PARRECError('Cannot create PARRECHeader from air.') if type(header) == klass: return header.copy() raise PARRECError('Cannot create PARREC header from non-PARREC header.') @classmethod def from_fileobj(klass, fileobj, permit_truncated=False, strict_sort=False): info, image_defs = parse_PAR_header(fileobj) return klass(info, image_defs, permit_truncated, strict_sort) def copy(self): return PARRECHeader( deepcopy(self.general_info), self.image_defs.copy(), self.permit_truncated, self.strict_sort, ) def as_analyze_map(self): """Convert PAR parameters to NIFTI1 format""" # Entries in the dict correspond to the parameters found in # the NIfTI1 header, specifically in nifti1.py `header_dtd` defs. # Here we set the parameters we can to simplify PAR/REC # to NIfTI conversion. descr = ( f'{self.general_info["exam_name"]};' f'{self.general_info["patient_name"]};' f'{self.general_info["exam_date"].replace(" ", "")};' f'{self.general_info["protocol_name"]}' )[:80] is_fmri = self.general_info['max_dynamics'] > 1 # PAR/REC uses msec, but in _calc_zooms we convert to sec t = 'sec' if is_fmri else 'unknown' xyzt_units = unit_codes['mm'] + unit_codes[t] return dict(descr=descr, xyzt_units=xyzt_units) # , pixdim=pixdim) def get_water_fat_shift(self): """Water fat shift, in pixels""" return self.general_info['water_fat_shift'] def get_echo_train_length(self): """Echo train length of the recording""" return self.general_info['epi_factor'] def get_q_vectors(self): """Get Q vectors from the data Returns ------- q_vectors : None or array Array of q vectors (bvals * bvecs), or None if not a diffusion acquisition. """ bvals, bvecs = self.get_bvals_bvecs() if bvals is None or bvecs is None: return None return bvecs * bvals[:, np.newaxis] def get_bvals_bvecs(self): """Get bvals and bvecs from data Returns ------- b_vals : None or array Array of b values, shape (n_directions,), or None if not a diffusion acquisition. b_vectors : None or array Array of b vectors, shape (n_directions, 3), or None if not a diffusion acquisition. """ if self.general_info['diffusion'] == 0: return None, None reorder = self.get_sorted_slice_indices() if len(self.get_data_shape()) == 3: # Any original diffusion scans will have >=2 volumes. However, a # single dynamic is possible for a post-processed diffusion volume # such as an ADC map. The b-values are unavailable in this case. return None, None else: n_slices, n_vols = self.get_data_shape()[-2:] bvals = self.image_defs['diffusion_b_factor'][reorder].reshape( (n_slices, n_vols), order='F' ) # All bvals within volume should be the same assert not np.any(np.diff(bvals, axis=0)) bvals = bvals[0] if 'diffusion' not in self.image_defs.dtype.names: return bvals, None bvecs = self.image_defs['diffusion'][reorder].reshape((n_slices, n_vols, 3), order='F') # All 3 values of bvecs should be same within volume assert not np.any(np.diff(bvecs, axis=0)) bvecs = bvecs[0] # rotate bvecs to match stored image orientation permute_to_psl = ACQ_TO_PSL[self.get_slice_orientation()] bvecs = apply_affine(np.linalg.inv(permute_to_psl), bvecs) return bvals, bvecs def get_def(self, name): """Return a single image definition field (or None if missing)""" idef = self.image_defs return idef[name] if name in idef.dtype.names else None def _get_unique_image_prop(self, name): """Scan image definitions and return unique value of a property. * Get array for named field of ``self.image_defs``; * Check that all rows in the array are the same and raise error otherwise; * Return the row. Parameters ---------- name : str Name of the property in ``self.image_defs`` Returns ------- unique_value : scalar or array Raises ------ PARRECError if the rows of ``self.image_defs[name]`` do not all compare equal. """ props = self.image_defs[name] if np.any(np.diff(props, axis=0)): raise PARRECError( f'Varying {name} in image sequence ({props}). This is not supported.' ) return props[0] def get_data_offset(self): """PAR header always has 0 data offset (into REC file)""" return 0 def set_data_offset(self, offset): """PAR header always has 0 data offset (into REC file)""" if offset != 0: raise PARRECError('PAR header assumes offset 0') def _calc_zooms(self): """Compute image zooms from header data. Spatial axis are first three. Returns ------- zooms : array Length 3 array for 3D image, length 4 array for 4D image. Notes ----- This routine gets called in ``__init__``, so may not be able to use some attributes available in the fully initialized object. """ # slice orientation for the whole image series slice_gap = self._get_unique_image_prop('slice gap') # scaling per image axis n_dim = 4 if self._get_n_vols() > 1 else 3 zooms = np.ones(n_dim) # spatial sizes are inplane X mm, inplane Y mm + inter slice gap zooms[:2] = self._get_unique_image_prop('pixel spacing') slice_thickness = self._get_unique_image_prop('slice thickness') zooms[2] = slice_thickness + slice_gap # If 4D dynamic scan, convert time from milliseconds to seconds if len(zooms) > 3 and self.general_info['dyn_scan']: if len(self.general_info['repetition_time']) > 1: warnings.warn('multiple TRs found in .PAR file') zooms[3] = self.general_info['repetition_time'][0] / 1000.0 return zooms def get_affine(self, origin='scanner'): """Compute affine transformation into scanner space. The method only considers global rotation and offset settings in the header and ignores potentially deviating information in the image definitions. Parameters ---------- origin : {'scanner', 'fov'} Transformation origin. By default the transformation is computed relative to the scanner's iso center. If 'fov' is requested the transformation origin will be the center of the field of view instead. Returns ------- aff : (4, 4) array 4x4 array, with output axis order corresponding to RAS or (x,y,z) or (lr, pa, fh). Notes ----- Transformations appear to be specified in (ap, fh, rl) axes. The orientation of data is recorded in the "slice orientation" field of the PAR header "General Information". We need to: * translate to coordinates in terms of the center of the FOV * apply voxel size scaling * reorder / flip the data to Philips' PSL axes * apply the rotations * apply any isocenter scaling offset if `origin` == "scanner" * reorder and flip to RAS axes """ # shape, zooms in original data ordering (ijk ordering) ijk_shape = np.array(self.get_data_shape()[:3]) to_center = from_matvec(np.eye(3), -(ijk_shape - 1) / 2.0) zoomer = np.diag(list(self.get_zooms()[:3]) + [1]) slice_orientation = self.get_slice_orientation() permute_to_psl = ACQ_TO_PSL.get(slice_orientation) if permute_to_psl is None: raise PARRECError(f'Unknown slice orientation ({slice_orientation}).') # hdr has deg, we need radians # Order is [ap, fh, rl] ap_rot, fh_rot, rl_rot = self.general_info['angulation'] * DEG2RAD Mx = euler2mat(x=ap_rot) My = euler2mat(y=fh_rot) Mz = euler2mat(z=rl_rot) # By trial and error, this unexpected order of rotations seem to give # the closest to the observed (converted NIfTI) affine. rot = from_matvec(dot_reduce(Mz, Mx, My)) # compose the PSL affine psl_aff = dot_reduce(rot, permute_to_psl, zoomer, to_center) if origin == 'scanner': # offset to scanner's isocenter (in ap, fh, rl) iso_offset = self.general_info['off_center'] psl_aff[:3, 3] += iso_offset # Currently in PSL; apply PSL -> RAS return np.dot(PSL_TO_RAS, psl_aff) def _get_n_slices(self): """Get number of slices for output data""" return len(set(self.image_defs['slice number'])) def _get_n_vols(self): """Get number of volumes for output data""" slice_nos = self.image_defs['slice number'] vol_nos = vol_numbers(slice_nos) is_full = vol_is_full(slice_nos, self.general_info['max_slices']) return len(set(np.array(vol_nos)[is_full])) def _calc_data_shape(self): """Calculate the output shape of the image data Returns length 3 tuple for 3D image, length 4 tuple for 4D. Returns ------- n_inplaneX : int number of voxels in X direction. n_inplaneY : int number of voxels in Y direction. n_slices : int number of slices. n_vols : int number of volumes or absent for 3D image. Notes ----- This routine gets called in ``__init__``, so may not be able to use some attributes available in the fully initialized object. """ inplane_shape = tuple(self._get_unique_image_prop('recon resolution')) shape = inplane_shape + (self._get_n_slices(),) n_vols = self._get_n_vols() return shape + (n_vols,) if n_vols > 1 else shape def get_data_scaling(self, method='dv'): """Returns scaling slope and intercept. Parameters ---------- method : {'fp', 'dv'} Scaling settings to be reported -- see notes below. Returns ------- slope : array scaling slope intercept : array scaling intercept Notes ----- The PAR header contains two different scaling settings: 'dv' (value on console) and 'fp' (floating point value). Here is how they are defined: DV = PV * RS + RI FP = DV / (RS * SS) where: PV: value in REC RS: rescale slope RI: rescale intercept SS: scale slope """ # These will be 3D or 4D scale_slope = self.image_defs['scale slope'] rescale_slope = self.image_defs['rescale slope'] rescale_intercept = self.image_defs['rescale intercept'] if method == 'dv': slope, intercept = rescale_slope, rescale_intercept elif method == 'fp': slope = 1.0 / scale_slope intercept = rescale_intercept / (rescale_slope * scale_slope) else: raise ValueError(f"Unknown scaling method '{method}'.") reorder = self.get_sorted_slice_indices() slope = slope[reorder] intercept = intercept[reorder] shape = (1, 1) + self.get_data_shape()[2:] slope = slope.reshape(shape, order='F') intercept = intercept.reshape(shape, order='F') return slope, intercept def get_slice_orientation(self): """Returns the slice orientation label. Returns ------- orientation : {'transverse', 'sagittal', 'coronal'} """ lab = self._get_unique_image_prop('slice orientation') return slice_orientation_codes.label[lab] def get_rec_shape(self): inplane_shape = tuple(self._get_unique_image_prop('recon resolution')) return inplane_shape + (len(self.image_defs),) def _strict_sort_order(self): """Determine the sort order based on several image definition fields. The fields taken into consideration, if present, are (in order from slowest to fastest variation after sorting): - image_defs['image_type_mr'] # Re, Im, Mag, Phase - image_defs['dynamic scan number'] # repetition - image_defs['label type'] # ASL tag/control - image_defs['diffusion b value number'] # diffusion b value - image_defs['gradient orientation number'] # diffusion directoin - image_defs['cardiac phase number'] # cardiac phase - image_defs['echo number'] # echo - image_defs['slice number'] # slice Data sorting is done in two stages: 1. an initial sort using the keys described above 2. a resort after generating two additional sort keys: * a key to assign unique volume numbers to any volumes that didn't have a unique sort based on the keys above (see :func:`vol_numbers`). * a sort key based on `vol_is_full` to identify truncated volumes A case where the initial sort may not create a unique label for each volume is diffusion scans acquired in the older V4 .PAR format, where diffusion direction info is not available. """ # sort keys present in all supported .PAR versions idefs = self.image_defs slice_nos = idefs['slice number'] dynamics = idefs['dynamic scan number'] phases = idefs['cardiac phase number'] echos = idefs['echo number'] image_type = idefs['image_type_mr'] # sort keys only present in a subset of .PAR files asl_keys = (idefs['label type'],) if 'label type' in idefs.dtype.names else () if self.general_info['diffusion'] != 0: bvals = self.get_def('diffusion b value number') if bvals is None: bvals = self.get_def('diffusion_b_factor') bvecs = self.get_def('gradient orientation number') if bvecs is None: # no b-vectors available diffusion_keys = (bvals,) else: diffusion_keys = (bvecs, bvals) else: diffusion_keys = () # initial sort (last key is highest precedence) keys = (slice_nos, echos, phases) + diffusion_keys + asl_keys + (dynamics, image_type) initial_sort_order = np.lexsort(keys) # sequentially number the volumes based on the initial sort vol_nos = vol_numbers(slice_nos[initial_sort_order]) # identify truncated volumes is_full = vol_is_full(slice_nos[initial_sort_order], self.general_info['max_slices']) # second stage of sorting return initial_sort_order[np.lexsort((vol_nos, is_full))] def _lax_sort_order(self): """ Sorts by (fast to slow): slice number, volume number. We calculate volume number by looking for repeating slice numbers (see :func:`vol_numbers`). """ slice_nos = self.image_defs['slice number'] is_full = vol_is_full(slice_nos, self.general_info['max_slices']) keys = (slice_nos, vol_numbers(slice_nos), np.logical_not(is_full)) return np.lexsort(keys) def get_sorted_slice_indices(self): """Return indices to sort (and maybe discard) slices in REC file. If the recording is truncated, the returned indices take care of discarding any slice indices from incomplete volumes. If `self.strict_sort` is True, a more complicated sorting based on multiple fields from the .PAR file is used. This may produce a different sort order than `strict_sort=False`, where volumes are sorted by the order in which the slices appear in the .PAR file. Returns ------- slice_indices : list List for indexing into the last (third) dimension of the REC data array, and (equivalently) the only dimension of ``self.image_defs``. """ if not self.strict_sort: sort_order = self._lax_sort_order() else: sort_order = self._strict_sort_order() # Figure out how many we need to remove from the end, and trim them. # Based on our sorting, they should always be last. n_used = np.prod(self.get_data_shape()[2:]) return sort_order[:n_used] def get_volume_labels(self): """Dynamic labels corresponding to the final data dimension(s). This is useful for custom data sorting. A subset of the info in ``self.image_defs`` is returned in an order that matches the final data dimension(s). Only labels that have more than one unique value across the dataset will be returned. Returns ------- sort_info : dict Each key corresponds to volume labels for a dynamically varying sequence dimension. The ordering of the labels matches the volume ordering determined via ``self.get_sorted_slice_indices``. """ sorted_indices = self.get_sorted_slice_indices() image_defs = self.image_defs # define which keys which might vary across image volumes dynamic_keys = [ 'cardiac phase number', 'echo number', 'label type', 'image_type_mr', 'dynamic scan number', 'scanning sequence', 'gradient orientation number', 'diffusion b value number', ] # remove dynamic keys that may not be present in older .PAR versions dynamic_keys = [d for d in dynamic_keys if d in image_defs.dtype.fields] non_unique_keys = [] for key in dynamic_keys: ndim = image_defs[key].ndim if ndim == 1: num_unique = len(np.unique(image_defs[key])) else: raise ValueError('unexpected image_defs shape > 1D') if num_unique > 1: non_unique_keys.append(key) # each key in dynamic keys will be identical across slices, so use # the value at slice 1. sl1_indices = image_defs['slice number'][sorted_indices] == 1 sort_info = OrderedDict() for key in non_unique_keys: sort_info[key] = image_defs[key][sorted_indices][sl1_indices] return sort_info class PARRECImage(SpatialImage): """PAR/REC image""" header_class = PARRECHeader header: PARRECHeader valid_exts = ('.rec', '.par') files_types = (('image', '.rec'), ('header', '.par')) makeable = False rw = False ImageArrayProxy = PARRECArrayProxy @classmethod def from_file_map( klass, file_map, *, mmap=True, permit_truncated=False, scaling='dv', strict_sort=False ): """Create PARREC image from file map `file_map` Parameters ---------- file_map : dict dict with keys ``image, header`` and values being fileholder objects for the respective REC and PAR files. mmap : {True, False, 'c', 'r'}, optional, keyword only `mmap` controls the use of numpy memory mapping for reading image array data. If False, do not try numpy ``memmap`` for data array. If one of {'c', 'r'}, try numpy memmap with ``mode=mmap``. A `mmap` value of True gives the same behavior as ``mmap='c'``. If image data file cannot be memory-mapped, ignore `mmap` value and read array from file. permit_truncated : {False, True}, optional, keyword-only If False, raise an error for an image where the header shows signs that fewer slices / volumes were recorded than were expected. scaling : {'dv', 'fp'}, optional, keyword-only Scaling method to apply to data (see :meth:`PARRECHeader.get_data_scaling`). strict_sort : bool, optional, keyword-only If True, a larger number of header fields are used while sorting the REC data array. This may produce a different sort order than `strict_sort=False`, where volumes are sorted by the order in which the slices appear in the .PAR file. """ with file_map['header'].get_prepare_fileobj('rt') as hdr_fobj: hdr = klass.header_class.from_fileobj( hdr_fobj, permit_truncated=permit_truncated, strict_sort=strict_sort ) rec_fobj = file_map['image'].get_prepare_fileobj() data = klass.ImageArrayProxy(rec_fobj, hdr, mmap=mmap, scaling=scaling) return klass(data, hdr.get_affine(), header=hdr, extra=None, file_map=file_map) @classmethod def from_filename( klass, filename, *, mmap=True, permit_truncated=False, scaling='dv', strict_sort=False ): """Create PARREC image from filename `filename` Parameters ---------- filename : str Filename of "PAR" or "REC" file mmap : {True, False, 'c', 'r'}, optional, keyword only `mmap` controls the use of numpy memory mapping for reading image array data. If False, do not try numpy ``memmap`` for data array. If one of {'c', 'r'}, try numpy memmap with ``mode=mmap``. A `mmap` value of True gives the same behavior as ``mmap='c'``. If image data file cannot be memory-mapped, ignore `mmap` value and read array from file. permit_truncated : {False, True}, optional, keyword-only If False, raise an error for an image where the header shows signs that fewer slices / volumes were recorded than were expected. scaling : {'dv', 'fp'}, optional, keyword-only Scaling method to apply to data (see :meth:`PARRECHeader.get_data_scaling`). strict_sort : bool, optional, keyword-only If True, a larger number of header fields are used while sorting the REC data array. This may produce a different sort order than `strict_sort=False`, where volumes are sorted by the order in which the slices appear in the .PAR file. """ file_map = klass.filespec_to_file_map(filename) return klass.from_file_map( file_map, mmap=mmap, permit_truncated=permit_truncated, scaling=scaling, strict_sort=strict_sort, ) load = from_filename # type: ignore[assignment] load = PARRECImage.from_filename