python-nibabel 1.2.2-1 / usr / share / pyshared / 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.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
"""Humble attempt to read images in PAR/REC format.

This is yet another MRI image format generated by Phillips
scanner. It is an ASCII header (PAR) plus a binary blob (REC).

This implementation aims to read version 4.2 of this format. Other versions
could probably be supported, but the author is lacking samples of them.
"""

import warnings
import numpy as np
import copy

from .spatialimages import SpatialImage, Header
from .eulerangles import euler2mat
from .volumeutils import Recoder

# PAR header versions we claim to understand
supported_versions = ['V4.2']

# 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': ('san_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),
    'Max. number of diffusion values': ('max_diffusion_values', int),
    'Max. number of gradient orients': ('max_gradient_orient', int),
    'Number of label types   <0=no ASL>': ('nr_label_types', int),
    }

# header items order per image definition line
image_def_dtd = [
    ('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', int,),
    ('window width', int,),
    ('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),
    ('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,)),
    ('label type', int,),                  # (imagekey!)
    ]
image_def_dtype = np.dtype(image_def_dtd)

# slice orientation codes
slice_orientation_codes = Recoder((# code, label
    (1, 'transversal'),
    (2, 'sagital'),
    (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


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
    -------
    (dict, array)
      The dictionary contains all "General Information" from the header file,
      while the (structured) has the properties of all image definitions in the
      header
    """
    # containers for relevant header lines
    general_info = {}
    image_info = []
    version = None

    # single pass through the header
    for line in fobj:
        # no junk
        line = line.strip()
        if line.startswith('#'):
            # try to get the header version
            if line.count('image export tool'):
                version = line.split()[-1]
                if not version in supported_versions:
                    warnings.warn(
                          "PAR/REC version '%s' 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."
                          % version)
            else:
                # just a comment
                continue
        elif line.startswith('.'):
            # read 'general information' and store in a dict
            first_colon = line[1:].find(':') + 1
            key = line[1:first_colon].strip()
            value = line[first_colon + 1:].strip()
            # 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
        elif line:
            # anything else is an image definition: store for later
            # processing
            image_info.append(line)

    # postproc image def props
    # create an array for all image defs
    image_defs = np.zeros(len(image_info), dtype=image_def_dtype)

    # for every image definition
    for i, line in enumerate(image_info):
        items = line.split()
        item_counter = 0
        # for all image properties we know about
        for props in image_def_dtd:
            if np.issubdtype(image_defs[props[0]].dtype, str):
                # simple string
                image_defs[props[0]][i] = items[item_counter]
                item_counter += 1
            elif len(props) == 2:
                # prop with numerical dtype
                image_defs[props[0]][i] = props[1](items[item_counter])
                item_counter += 1
            elif len(props) == 3:
                # array prop with dtype
                nelements = np.prod(props[2])
                # get as many elements as necessary
                itms = items[item_counter:item_counter + nelements]
                # convert to array with dtype
                value = np.fromstring(" ".join(itms), props[1], sep=' ')
                # store
                image_defs[props[0]][i] = value
                item_counter += nelements

    return general_info, image_defs


class PARRECHeader(Header):
    """PAR/REC header"""
    def __init__(self, info, image_defs):
        """
        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()`).
        """
        self.general_info = info
        self.image_defs = image_defs
        self._slice_orientation = None

        # charge with basic properties to be able to use base class
        # functionality
        # dtype
        dtype = np.typeDict[
                    'int'
                    + str(self._get_unique_image_prop('image pixel size')[0])]
        Header.__init__(self,
                        data_dtype=dtype,
                        shape=self.get_data_shape_in_file(),
                        zooms=self._get_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):
        info, image_defs = parse_PAR_header(fileobj)
        return klass(info, image_defs)


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


    def _get_unique_image_prop(self, name):
        """Scan all image definitions and return the unique value of a property.

        If the requested property is an array this method behave _not_ like
        `np.unique`. It will return the unique combination of all array elements
        for any image definition, and _not_ the unique element values.

        Raises
        ------
        If there is more than a single unique value a `PARRECError` is raised.
        """
        prop = self.image_defs[name]
        if len(prop.shape) > 1:
            uprops = [np.unique(prop[i]) for i in range(len(prop.shape))]
        else:
            uprops = [np.unique(prop)]
        if not np.prod([len(uprop) for uprop in uprops]) == 1:
            raise PARRECError('Varying %s in image sequence (%s). This is not '
                              'suppported.' % (name, uprops))
        else:
            return np.array([uprop[0] for uprop in uprops])


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

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


    def get_ndim(self):
        """Return the number of dimensions of the image data."""
        if self.general_info['max_dynamics'] > 1 \
           or self.general_info['max_gradient_orient'] > 1:
            return 4
        else:
            return 3


    def _get_zooms(self):
        """Compute image zooms from header data.

        Spatial axis are first three.
        """
        # slice orientation for the whole image series
        slice_gap = self._get_unique_image_prop('slice gap')[0]
        # scaling per image axis
        zooms = np.ones(self.get_ndim())
        # spatial axes correspond to voxelsize + inter slice gap
        # voxel size (inplaneX, inplaneY, slices)
        zooms[:3] = self.get_voxel_size()
        zooms[2] += slice_gap
        # time axis?
        if len(zooms) > 3  and self.general_info['max_dynamics'] > 1:
            # DTI also has 4D
            zooms[3] = self.general_info['repetition_time']
        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 ignore 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
        -------
        array
          4x4 array, with axis order corresponding to (x,y,z) or (lr, pa, fh).
        """
        # hdr has deg, we need radian
        # order is [ap, fh, rl]
        ang_rad = self.general_info['angulation'] * np.pi / 180.0
        # need to rotate back from what was given in the file
        ang_rad *= -1

        # R2AGUI approach is this, but it comes with remarks ;-)
        # % trying to incorporate AP FH RL rotation angles: determined using some 
        # % common sense, Chris Rordon's help + source code and trial and error, 
        # % this is considered EXPERIMENTAL!
        rot_rl = np.mat(
                [[1.0, 0.0, 0.0],
                 [0.0, np.cos(ang_rad[2]), -np.sin(ang_rad[2])],
                 [0.0, np.sin(ang_rad[2]), np.cos(ang_rad[2])]]
                )
        rot_ap = np.mat(
                [[np.cos(ang_rad[0]), 0.0, np.sin(ang_rad[0])],
                 [0.0, 1.0, 0.0],
                 [-np.sin(ang_rad[0]), 0.0, np.cos(ang_rad[0])]]
                )
        rot_fh = np.mat(
                [[np.cos(ang_rad[1]), -np.sin(ang_rad[1]), 0.0],
                 [np.sin(ang_rad[1]), np.cos(ang_rad[1]), 0.0],
                 [0.0, 0.0, 1.0]]
                )
        rot_r2agui = rot_rl * rot_ap * rot_fh
        # NiBabel way of doing it
        # order is [ap, fh, rl]
        #           x   y   z
        #           0   1   2
        rot_nibabel = euler2mat(ang_rad[1], ang_rad[0], ang_rad[2])

        # XXX for now put some safety net, until we have recorded proper
        # test data with oblique orientations and different readout directions
        # to verify the order of arguments of euler2mat
        assert(np.all(rot_r2agui == rot_nibabel))
        rot = rot_nibabel

        # FOV (always in ap, fh, rl)
        fov = self.general_info['fov']
        # voxel size always (inplaneX, inplaneY, slicethickness (without gap))
        voxsize = self.get_voxel_size()

        slice_orientation = self.get_slice_orientation()
        if slice_orientation == 'sagital':
            # inplane: AP, FH   slices: RL
            recfg_data_axis = np.mat([[  0,  0,  1],
                                      [ -1,  0,  0],
                                      [  0, -1,  0]])
            # fov is already like the data
            fov = fov
        elif slice_orientation == 'transversal':
            # inplane: RL, AP   slices: FH
            recfg_data_axis = np.mat([[ -1,  0,  0],
                                      [  0, -1,  0],
                                      [  0,  0,  1]])
            # fov is already like the data
            fov = fov[[2,0,1]]
        elif slice_orientation == 'coronal':
            # inplane: RL, FH   slices: AP
            recfg_data_axis = np.mat([[ -1,  0,  0],
                                      [  0,  0, -1],
                                      [  0, -1,  0]])
            # fov is already like the data
            fov = fov[[2,1,0]]
        else:
            raise PARRECError("Unknown slice orientation (%s)."
                              % slice_orientation)

        rot = rot * recfg_data_axis

        # ijk origin should be: Anterior, Right, Foot
        # qform should point to the center of the voxel
        fov_center_offset = self.get_voxel_size() / 2 - fov / 2

        # need to rotate this offset into scanner space
        fov_center_offset = np.dot(rot, fov_center_offset)

        # get the scaling by voxelsize and slice thickness (incl. gap)
        scaled = rot * np.mat(np.diag(self.get_zooms()[:3]))

        # compose the affine
        aff = np.eye(4)
        aff[:3,:3] = scaled
        # offset
        aff[:3,3] = fov_center_offset
        if origin == 'fov':
            pass
        elif origin == 'scanner':
            # offset to scanner's iso center (always in ap, fh, rl)
            # -- turn into rl, ap, fh and then lr, pa, fh
            iso_offset = self.general_info['off_center'][[2,0,1]] * [-1,-1,0]
            aff[:3,3] += iso_offset
        return aff


    def get_data_shape_in_file(self):
        """Return the shape of the binary blob in the REC file.

        Returns
        -------
        tuple
          (inplaneX, inplaneY, nslices, ndynamics/ndirections)
        """
        # e.g. number of volumes
        ndynamics = len(np.unique(self.image_defs['dynamic scan number']))
        # DTI volumes (b-values-1 x directions)
        # there is some awkward exception to this rule for b-values > 2
        # XXX need to get test image...
        ndtivolumes = (self.general_info['max_diffusion_values'] - 1) \
                        * self.general_info['max_gradient_orient']
        nslices = len(np.unique(self.image_defs['slice number']))
        if not nslices == self.general_info['max_slices']:
            raise PARRECError("Header inconsistency: Found %i slices, "
                              "but header claims to have %i."
                              % (nslices, self.general_info['max_slices']))

        inplane_shape = tuple(self._get_unique_image_prop('recon resolution'))

        # there should not be both: multiple dynamics and DTI
        if ndynamics > 1:
            return inplane_shape + (nslices, ndynamics)
        elif ndtivolumes > 1:
            return inplane_shape + (nslices, ndtivolumes)
        else:
            return tuple(inplane_shape) + (nslices,)


    def get_data_scaling(self, method="dv"):
        """Returns scaling slope and intercept.

        Parameters
        ----------
        method : {'fp', 'dv'}
          Scaling settings to be reported -- see notes below.

        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)
        """
        # XXX: FP tends to become HUGE, DV seems to be more reasonable -> figure
        #      out which one means what

        # although the is a per-image scaling in the header, it looks like
        # there is just one unique factor and intercept per whole image series
        scale_slope = self._get_unique_image_prop('scale slope')
        rescale_slope = self._get_unique_image_prop('rescale slope')
        rescale_intercept = self._get_unique_image_prop('rescale intercept')

        if method == 'dv':
            slope = rescale_slope
            intercept = rescale_intercept
        elif method == 'fp':
            # actual slopes per definition above
            slope = 1.0 / scale_slope
            # actual intercept per definition above
            intercept = rescale_intercept / (rescale_slope * scale_slope)
        else:
            raise ValueError("Unknown scling method '%s'." % method)
        return (slope, intercept)


    def get_slice_orientation(self):
        """Returns the slice orientation label.

        Returns
        -------
        {'transversal', 'sagital', 'coronal'}
        """
        if self._slice_orientation is None:
            self._slice_orientation = \
                slice_orientation_codes.label[
                    self._get_unique_image_prop('slice orientation')[0]]
        return self._slice_orientation


    def raw_data_from_fileobj(self, fileobj):
        """Returns memmap array of raw unscaled image data.

        Array axes correspond to x,y,z,t.
        """
        # memmap the data -- it is guaranteed to be uncompressed and all
        # properties are known
        # read in Fortran order to have spatial axes first
        data = np.memmap(fileobj,
                         dtype=self.get_data_dtype(),
                         mode='c', # copy-on-write
                         shape=self.get_data_shape_in_file(),
                         order='F')
        return data


    def data_from_fileobj(self, fileobj):
        """Returns scaled image data.

        Behaves identical to `PARRECHeader.raw_data_from_fileobj()`, but
        returns scaled image data. This causes the images data to be loaded into
        memory.
        """
        unscaled = self.raw_data_from_fileobj(fileobj)
        slope, intercept = self.get_data_scaling()
        scaled = unscaled * slope
        scaled += intercept
        return scaled



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

    class ImageArrayProxy(object):
        def __init__(self, rec_fobj, hdr):
            self._rec_fobj = rec_fobj
            self._hdr = hdr
            self._data = None

        def __array__(self):
            ''' Cached read of data from file '''
            if self._data is None:
                self._data = self._hdr.data_from_fileobj(self._rec_fobj)
            return self._data

        @property
        def shape(self):
            # embedded header knows it, without having to touch the data
            return self._hdr.get_data_shape()


    @classmethod
    def from_file_map(klass, file_map):
        hdr_fobj = file_map['header'].get_prepare_fileobj()
        rec_fobj = file_map['image'].get_prepare_fileobj()
        hdr = PARRECHeader.from_fileobj(hdr_fobj)
        data = klass.ImageArrayProxy(rec_fobj, hdr)
        return klass(data,
                     hdr.get_affine(),
                     header=hdr,
                     extra=None,
                     file_map=file_map)


load = PARRECImage.load