python3-nibabel 2.0.2-2 / usr / lib / python3 / 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.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
"""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 and 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 nibabe/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's 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
* fh = foot -> head = S in RAS
* rl = right -> left = negative R in RAS

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.

We call the PAR file's axis system "PSL" (Posterior, Superior, Left)

#########
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
"""
from __future__ import print_function, division

import warnings
import numpy as np
from copy import deepcopy
import re

from .keywordonly import kw_only_meth
from .spatialimages import SpatialImage, Header
from .eulerangles import euler2mat
from .volumeutils import Recoder, array_from_file, BinOpener
from .affines import from_matvec, dot_reduce, apply_affine
from .nifti1 import unit_codes
from .fileslice import fileslice, strided_scalar

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

# General information dict definitions
# assign props to PAR header entries
# values are: (shortname[, dtype[, shape]])
_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),
    '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.
    """
    pass


# 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=' ')
            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 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
        Lowest slice number for full slice set.

    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 `slice_nos` value is outside slice set.
    """
    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 not gdef_max_name 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):
    """Get 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] + scalings[1]
    return rec_data


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 BinOpener(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 BinOpener(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 BinOpener(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(Header):
    """PAR/REC header"""
    def __init__(self, info, image_defs, permit_truncated=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.
        """
        self.general_info = info.copy()
        self.image_defs = image_defs.copy()
        self.permit_truncated = permit_truncated
        _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('<')
        Header.__init__(self,
                        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):
        info, image_defs = parse_PAR_header(fileobj)
        return klass(info, image_defs, permit_truncated)

    def copy(self):
        return PARRECHeader(deepcopy(self.general_info),
                            self.image_defs.copy(),
                            self.permit_truncated)

    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 and 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]
        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_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]

    def get_voxel_size(self):
        """Returns the spatial extent of a voxel.

        Does not include the slice gap in the slice extent.

        This function is deprecated and we will remove it in future versions of
        nibabel.  Please use ``get_zooms`` instead.  If you need the slice
        thickness not including the slice gap, use ``self.image_defs['slice
        thickness']``.

        Returns
        -------
        vox_size: shape (3,) ndarray
        """
        warnings.warn('Please use "get_zooms" instead of "get_voxel_size"',
                      DeprecationWarning,
                      stacklevel=2)
        # 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 called in ``__init__``, so may not be able to use
        some attributes available in the fully initalized 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']:
            zooms[3] = self.general_info['repetition_time'] / 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]
        ang_rad = self.general_info['angulation'] * np.pi / 180.0
        # euler2mat accepts z, y, x angles and does rotation around z, y, x
        # axes in that order. It's possible that PAR assumes rotation in a
        # different order, we still need some relevant data to test this
        rot = from_matvec(euler2mat(*ang_rad[::-1]), [0, 0, 0])
        # 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 called in ``__init__``, so may not be able to use
        some attributes available in the fully initalized 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:

        PV: value in REC
        RS: rescale slope
        RI: rescale intercept
        SS: scale slope

        DV = PV * RS + RI
        FP = DV / (RS * SS)
        """
        # 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 get_sorted_slice_indices(self):
        """Indices to sort (and maybe discard) slices in REC file

        Returns list for indexing into the last (third) dimension of the REC
        data array, and (equivalently) the only dimension of
        ``self.image_defs``.

        If the recording is truncated, the returned indices take care of
        discarding any indices that are not meant to be used.
        """
        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))
        # 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 np.lexsort(keys)[:n_used]


class PARRECImage(SpatialImage):
    """PAR/REC image"""
    header_class = PARRECHeader
    files_types = (('image', '.rec'), ('header', '.par'))

    ImageArrayProxy = PARRECArrayProxy

    @classmethod
    @kw_only_meth(1)
    def from_file_map(klass, file_map, mmap=True, permit_truncated=False,
                      scaling='dv'):
        """ 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`).
        """
        with file_map['header'].get_prepare_fileobj('rt') as hdr_fobj:
            hdr = klass.header_class.from_fileobj(
                hdr_fobj,
                permit_truncated=permit_truncated)
        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'):
        """ 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`).
        """
        file_map = klass.filespec_to_file_map(filename)
        return klass.from_file_map(file_map,
                                   mmap=mmap,
                                   permit_truncated=permit_truncated,
                                   scaling=scaling)

    load = from_filename


load = PARRECImage.load