python-nibabel 2.1.0-1 / usr / lib / python2.7 / dist-packages / nibabel / parrec.py

# 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.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
# Disable line length checking for PAR fragments in module docstring
# flake8: noqa E501
"""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".
"""
from __future__ import print_function, division

import warnings
import numpy as np
from copy import deepcopy
import re
from io import StringIO
from locale import getpreferredencoding
from nibabel.externals import OrderedDict

from .keywordonly import kw_only_meth
from .spatialimages import SpatialHeader, SpatialImage
from .eulerangles import euler2mat
from .volumeutils import Recoder, array_from_file
from .affines import from_matvec, dot_reduce, apply_affine
from .nifti1 import unit_codes
from .fileslice import fileslice, strided_scalar
from .openers import ImageOpener
from .deprecated import deprecate_with_version

# 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.

# 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),
}

# 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
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!)
]

#: 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 = {}
    for key, value in gen_dict.items():
        # get props for this hdr field
        props = _hdr_key_dict[key]
        # 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.shape = 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(
            'Slice numbers outside inclusive range {0} to {1}'.format(
                slice_min, 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(
                "Header inconsistency: Found {0} {1} values, "
                "but expected {2}".format(n_have, idef_name, 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(
            """ PAR/REC version '{0}' 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.
            """.format(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(object):

    @kw_only_meth(2)
    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 dtype(self):
        return self._dtype

    @property
    def is_proxy(self):
        return True

    def get_unscaled(self):
        with ImageOpener(self.file_like) as fileobj:
            return _data_from_rec(fileobj, self._rec_shape, self._dtype,
                                  self._slice_indices, self._shape,
                                  mmap=self._mmap)

    def __array__(self):
        with ImageOpener(self.file_like) as fileobj:
            return _data_from_rec(fileobj,
                                  self._rec_shape,
                                  self._dtype,
                                  self._slice_indices,
                                  self._shape,
                                  scalings=self._slice_scaling,
                                  mmap=self._mmap)

    def __getitem__(self, slicer):
        indices = self._slice_indices
        if indices[0] != 0 or np.any(np.diff(indices) != 1):
            # We can't load direct from REC file, use inefficient slicing
            return np.asanyarray(self)[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:
            raw_data = fileslice(fileobj, slicer, self._shape, self._dtype, 0,
                                 'F')
        # Broadcast scaling to shape of original data
        slopes, inters = self._slice_scaling
        fake_data = strided_scalar(self._shape)
        _, slopes, inters = np.broadcast_arrays(fake_data, slopes, inters)
        # Slice scaling to give output shape
        return raw_data * slopes[slicer] + inters[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('Only 8- and 16-bit data supported (not %s)'
                              'please report this to the nibabel developers'
                              % bitpix)
        # REC data always little endian
        dt = np.dtype('uint' + str(bitpix)).newbyteorder('<')
        super(PARRECHeader, self).__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 = ("%s;%s;%s;%s"
                 % (self.general_info['exam_name'],
                    self.general_info['patient_name'],
                    self.general_info['exam_date'].replace(' ', ''),
                    self.general_info['protocol_name']))[:80]  # max len
        is_fmri = (self.general_info['max_dynamics'] > 1)
        t = 'msec' 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()
        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('Varying {0} in image sequence ({1}). This is '
                              'not suppported.'.format(name, props))
        return props[0]

    @deprecate_with_version('get_voxel_size deprecated. '
                            'Please use "get_zooms" instead.',
                            '2.0', '4.0')
    def get_voxel_size(self):
        """Returns the spatial extent of a voxel.

        Does not include the slice gap in the slice extent.

        If you need the slice thickness not including the slice gap, use
        ``self.image_defs['slice thickness']``.

        Returns
        -------
        vox_size: shape (3,) ndarray
        """
        # slice orientation for the whole image series
        slice_thickness = self._get_unique_image_prop('slice thickness')
        voxsize_inplane = self._get_unique_image_prop('pixel spacing')
        voxsize = np.array((voxsize_inplane[0],
                            voxsize_inplane[1],
                            slice_thickness))
        return voxsize

    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.
        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.)
        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(
                "Unknown slice orientation ({0}).".format(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("Unknown scaling method '%s'." % 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
    valid_exts = ('.rec', '.par')
    files_types = (('image', '.rec'), ('header', '.par'))

    makeable = False
    rw = False

    ImageArrayProxy = PARRECArrayProxy

    @classmethod
    @kw_only_meth(1)
    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
    @kw_only_meth(1)
    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


load = PARRECImage.load