Package mvpa :: Package datasets :: Module nifti

[hide private]
[frames] | no frames]

# Source Code for Module mvpa.datasets.nifti

      1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
      2  # vi: set ft=python sts=4 ts=4 sw=4 et:
      3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
      4  #
      5  #   See COPYING file distributed along with the PyMVPA package for the
      6  #   copyright and license terms.
      7  #
      8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
      9  """Dataset that gets its samples from a NIfTI file"""
     10
     11  __docformat__ = 'restructuredtext'
     12
     13  from mvpa.base import externals
     14
     15  import sys
     16  import numpy as N
     17  from mvpa.support.copy import deepcopy
     18
     19  if __debug__:
     20      from mvpa.base import debug
     21
     22  if externals.exists('nifti', raiseException=True):
     23      if sys.version_info[:2] >= (2, 5):
     24          # enforce absolute import
     25          NiftiImage = __import__('nifti', globals(), locals(), [], 0).NiftiImage
     26      else:
     27          # little trick to be able to import 'nifti' package (which has same
     28          # name)
     29          oldname = __name__
     30          # crazy name with close to zero possibility to cause whatever
     31          __name__ = 'iaugf9zrkjsbdv89'
     32          from nifti import NiftiImage
     33          # restore old settings
     34          __name__ = oldname
     35
     36  from mvpa.datasets.base import Dataset
     37  from mvpa.datasets.mapped import MappedDataset
     38  from mvpa.datasets.event import EventDataset
     39  from mvpa.mappers.base import CombinedMapper
     40  from mvpa.mappers.metric import DescreteMetric, cartesianDistance
     41  from mvpa.mappers.array import DenseArrayMapper
     42  from mvpa.base import warning
     43
     44
     45 -def getNiftiFromAnySource(src, ensure=False, enforce_dim=None, scale_data=True):
     46      """Load/access NIfTI data from files or instances.
     47
     48      :Parameters:
     49        src: str | NiftiImage
     50          Filename of a NIfTI image or a `NiftiImage` instance.
     51        ensure : bool
     52          If True, through ValueError exception if cannot be loaded.
     53        enforce_dim : int or None
     54          If not None, it is the dimensionality of the data to be enforced,
     55          commonly 4D for the data, and 3D for the mask in case of fMRI.
     56        scale_data : bool
     57          NIfTI header specifies scl_slope and scl_inter for scaling and
     58          offsetting the data.  By default those will get applied to the data
     59          (change in behavior post 0.4.6).
     60
     61      :Returns:
     62        NiftiImage | None
     63          If the source is not supported None is returned.
     64      """
     65      nifti = None
     66
     67      # figure out what type
     68      if isinstance(src, str):
     69          # open the nifti file
     70          try:
     71              nifti = NiftiImage(src)
     72          except RuntimeError, e:
     73              warning("ERROR: NiftiDatasets: Cannot open NIfTI file %s" \
     74                      % src)
     75              raise e
     76      elif isinstance(src, NiftiImage):
     77          # nothing special
     78          nifti = src
     79      elif (isinstance(src, list) or isinstance(src, tuple)) \
     80          and len(src)>0 \
     81          and (isinstance(src[0], str) or isinstance(src[0], NiftiImage)):
     82          # load from a list of given entries
     83          if enforce_dim is not None: enforce_dim_ = enforce_dim - 1
     84          else:                       enforce_dim_ = None
     85          if __debug__:
     86              debug('DS_NIFTI', 'Loading from a sequence of sources: %s' % (src,))
     87          srcs = [getNiftiFromAnySource(s, ensure=ensure,
     88                                        enforce_dim=enforce_dim_)
     89                  for s in src]
     90          if __debug__:
     91              # lets check if they all have the same dimensionality
     92              shapes = [s.data.shape for s in srcs]
     93              if not N.all([s == shapes[0] for s in shapes]):
     94                  raise ValueError, \
     95                        "Input volumes contain variable number of dimensions:" \
     96                        " %s" % (shapes,)
     97          # Combine them all into a single beast
     98          # And since they all could have varying scl_* - reset those
     99          hdr = srcs[0].header
    100          hdr['scl_slope'] = 1.
    101          hdr['scl_inter'] = 0.
    102          nifti = NiftiImage(N.array([s.asarray() for s in srcs]), hdr)
    103      elif ensure:
    104          raise ValueError, "Cannot load NIfTI from %s" % (src,)
    105
    106      if nifti is not None and enforce_dim is not None:
    107          shape, new_shape = nifti.data.shape, None
    108          lshape = len(shape)
    109
    110          # check if we need to tune up shape
    111          if lshape < enforce_dim:
    112              # if we are missing required dimension(s)
    113              new_shape = (1,)*(enforce_dim-lshape) + shape
    114          elif lshape > enforce_dim:
    115              # if there are bogus dimensions at the beginning
    116              bogus_dims = lshape - enforce_dim
    117              if shape[:bogus_dims] != (1,)*bogus_dims:
    118                  raise ValueError, \
    119                        "Cannot enforce %dD on data with shape %s" \
    120                        % (enforce_dim, shape)
    121              new_shape = shape[bogus_dims:]
    122
    123          # tune up shape if needed
    124          if new_shape is not None:
    125              if __debug__:
    126                  debug('DS_NIFTI', 'Enforcing shape %s for %s data from %s' %
    127                        (new_shape, shape, src))
    128              nifti.data.shape = new_shape
    129
    130      if nifti is not None and scale_data:
    131          if nifti.slope and not (nifti.slope == 1.0 and nifti.intercept == 0.0):
    132              if __debug__:
    133                  debug('DS_NIFTI', 'Scaling the data from %s' % (src,))
    134              nifti.data = nifti.getScaledData()
    135          else:
    136              # Do nothing -- just debug message
    137              if __debug__:
    138                  debug('DS_NIFTI', 'Although scaling was requested, data from %s'
    139                        ' has no scaling parameters set -- thus no scaling' % (src,))
    140      return nifti
    141
    142
    143 -def getNiftiData(nim):
    144      """Convenience function to extract the data array from a NiftiImage
    145
    146      This function will make use of advanced features of PyNIfTI to prevent
    147      unnecessary copying if a sufficent version is available.
    148      """
    149      if externals.exists('nifti ge 0.20090205.1'):
    150          return nim.data
    151      else:
    152          return nim.asarray()
    153
    154 -def _get_safe_header(nids):
    155      """Given *NiftiDataset, returns a copy of NIfTI header with reset scl_ fields
    156      """
    157      hdr = nids.niftihdr.copy()
    158      hdr['scl_slope'] = 1.
    159      hdr['scl_inter'] = 0.
    160      return hdr
    161
    162 -class NiftiDataset(MappedDataset):
    163      """Dataset loading its samples from a NIfTI image or file.
    164
    165      Samples can be loaded from a NiftiImage instance or directly from a NIfTI
    166      file. This class stores all relevant information from the NIfTI file header
    167      and provides information about the metrics and neighborhood information of
    168      all voxels.
    169
    170      Most importantly it allows to map data back into the original data space
    171      and format via :meth:`~mvpa.datasets.nifti.NiftiDataset.map2Nifti`.
    172
    173      This class allows for convenient pre-selection of features by providing a
    174      mask to the constructor. Only non-zero elements from this mask will be
    175      considered as features.
    176
    177      NIfTI files are accessed via PyNIfTI. See
    178      http://niftilib.sourceforge.net/pynifti/ for more information about
    179      pynifti.
    180      """
    181      # XXX: Every dataset should really have an example of howto instantiate
    182      #      it (necessary parameters).
    183 -    def __init__(self, samples=None, mask=None, dsattr=None,
    184                   enforce_dim=4, scale_data=True, **kwargs):
    185          """
    186          :Parameters:
    187            samples: str | NiftiImage
    188              Filename of a NIfTI image or a `NiftiImage` instance.
    189            mask: str | NiftiImage | ndarray
    190              Filename of a NIfTI image or a `NiftiImage` instance or an ndarray
    191              of appropriate shape.
    192            enforce_dim : int or None
    193              If not None, it is the dimensionality of the data to be enforced,
    194              commonly 4D for the data, and 3D for the mask in case of fMRI.
    195            scale_data : bool
    196              NIfTI header specifies scl_slope and scl_inter for scaling and
    197              offsetting the data.  By default those will get applied to the data
    198              (change in behavior post 0.4.6).
    199          """
    200          # if in copy constructor mode
    201          if not dsattr is None and dsattr.has_key('mapper'):
    202              MappedDataset.__init__(self,
    203                                     samples=samples,
    204                                     dsattr=dsattr,
    205                                     **kwargs)
    206              return
    207
    208          #
    209          # the following code only deals with contructing fresh datasets from
    210          # scratch
    211          #
    212
    213          # load the samples
    214          niftisamples = getNiftiFromAnySource(samples, ensure=True,
    215                                               enforce_dim=enforce_dim,
    216                                               scale_data=scale_data)
    217          samples = niftisamples.data
    218
    219          # do not put the whole NiftiImage in the dict as this will most
    220          # likely be deepcopy'ed at some point and ensuring data integrity
    221          # of the complex Python-C-Swig hybrid might be a tricky task.
    222          # Only storing the header dict should achieve the same and is more
    223          # memory efficient and even simpler
    224          dsattr = {'niftihdr': niftisamples.header}
    225
    226
    227          # figure out what the mask is, but onyl handle known cases, the rest
    228          # goes directly into the mapper which maybe knows more
    229          niftimask = getNiftiFromAnySource(mask, scale_data=scale_data)
    230          if niftimask is None:
    231              pass
    232          elif isinstance(niftimask, N.ndarray):
    233              mask = niftimask
    234          else:
    235              mask = getNiftiData(niftimask)
    236
    237          # build an appropriate mapper that knows about the metrics of the NIfTI
    238          # data
    239          # NiftiDataset uses a DescreteMetric with cartesian
    240          # distance and element size from the NIfTI header
    241
    242          # 'voxdim' is (x,y,z) while 'samples' are (t,z,y,x)
    243          elementsize = [i for i in reversed(niftisamples.voxdim)]
    244          mapper = DenseArrayMapper(mask=mask, shape=samples.shape[1:],
    245                      metric=DescreteMetric(elementsize=elementsize,
    246                                            distance_function=cartesianDistance))
    247
    248          MappedDataset.__init__(self,
    249                                 samples=samples,
    250                                 mapper=mapper,
    251                                 dsattr=dsattr,
    252                                 **kwargs)
    253
    254
    255 -    def map2Nifti(self, data=None):
    256          """Maps a data vector into the dataspace and wraps it with a
    257          NiftiImage. The header data of this object is used to initialize
    258          the new NiftiImage (scl_slope and scl_inter are reset to 1.0 and
    259          0.0 accordingly).
    260
    261          :Parameters:
    262            data : ndarray or Dataset
    263              The data to be wrapped into NiftiImage. If None (default), it
    264              would wrap samples of the current dataset. If it is a Dataset
    265              instance -- takes its samples for mapping
    266          """
    267          if data is None:
    268              data = self.samples
    269          elif isinstance(data, Dataset):
    270              # ease users life
    271              data = data.samples
    272          dsarray = self.mapper.reverse(data)
    273
    274          return NiftiImage(dsarray, _get_safe_header(self))
    275
    276
    277 -    def getDt(self):
    278          """Return the temporal distance of two samples/volumes.
    279
    280          This method tries to be clever and always returns `dt` in seconds, by
    281          using unit information from the NIfTI header. If such information is
    282          not present the assumed unit will also be `seconds`.
    283          """
    284          # plain value
    285          hdr = self.niftihdr
    286          TR = hdr['pixdim'][4]
    287
    288          # by default assume seconds as unit and do not scale
    289          scale = 1.0
    290
    291          # figure out units, if available
    292          if hdr.has_key('time_unit'):
    293              unit_code = hdr['time_unit'] / 8
    294          elif hdr.has_key('xyzt_unit'):
    295              unit_code = int(hdr['xyzt_unit']) / 8
    296          else:
    297              warning("No information on time units is available. Assuming "
    298                      "seconds")
    299              unit_code = 0
    300
    301          # handle known units
    302          # XXX should be refactored to use actual unit labels from pynifti
    303          # when version 0.20090205 or later is assumed to be available on all
    304          # machines
    305          if unit_code in [0, 1, 2, 3]:
    306              if unit_code == 0:
    307                  warning("Time units were not specified in NiftiImage. "
    308                          "Assuming seconds.")
    309              scale = [ 1.0, 1.0, 1e-3, 1e-6 ][unit_code]
    310          else:
    311              warning("Time units are incorrectly coded: value %d whenever "
    312                      "allowed are 8 (sec), 16 (millisec), 24 (microsec). "
    313                      "Assuming seconds." % (unit_code * 8,)
    314                      )
    315          return TR * scale
    316
    317
    318      niftihdr = property(fget=lambda self: self._dsattr['niftihdr'],
    319                          doc='Access to the NIfTI header dictionary.')
    320
    321      dt = property(fget=getDt,
    322                    doc='Time difference between two samples (in seconds). '
    323                    'AKA TR in fMRI world.')
    324
    325      samplingrate = property(fget=lambda self: 1.0 / self.dt,
    326                            doc='Sampling rate (based on .dt).')
    327
    328
    329 -class ERNiftiDataset(EventDataset):
    330      """Dataset with event-defined samples from a NIfTI timeseries image.
    331
    332      This is a convenience dataset to facilitate the analysis of event-related
    333      fMRI datasets. Boxcar-shaped samples are automatically extracted from the
    334      full timeseries using :class:`~mvpa.misc.support.Event` definition lists.
    335      For each event all volumes covering that particular event in time
    336      (including partial coverage) are used to form the corresponding sample.
    337
    338      The class supports the conversion of events defined in 'realtime' into the
    339      descrete temporal space defined by the NIfTI image. Moreover, potentially
    340      varying offsets between true event onset and timepoint of the first selected
    341      volume can be stored as an additional feature in the dataset.
    342
    343      Additionally, the dataset supports masking. This is done similar to the
    344      masking capabilities of :class:`~mvpa.datasets.nifti.NiftiDataset`. However,
    345      the mask can either be of the same shape as a single NIfTI volume, or
    346      can be of the same shape as the generated boxcar samples, i.e.
    347      a samples consisting of three volumes with 24 slices and 64x64 inplane
    348      resolution needs a mask with shape (3, 24, 64, 64). In the former case the
    349      mask volume is automatically expanded to be identical in a volumes of the
    350      boxcar.
    351      """
    352 -    def __init__(self, samples=None, events=None, mask=None, evconv=False,
    353                   storeoffset=False, tr=None, enforce_dim=4,
    354                   scale_data=True, **kwargs):
    355          """
    356          :Parameters:
    357            mask: str | NiftiImage | ndarray
    358              Filename of a NIfTI image or a `NiftiImage` instance or an ndarray
    359              of appropriate shape.
    360            evconv: bool
    361              Convert event definitions using `onset` and `duration` in some
    362              temporal unit into #sample notation.
    363            storeoffset: bool
    364              Whether to store temproal offset information when converting
    365              Events into descrete time. Only considered when evconv == True.
    366            tr: float
    367              Temporal distance of two adjacent NIfTI volumes. This can be used
    368              to override the corresponding value in the NIfTI header.
    369            enforce_dim : int or None
    370              If not None, it is the dimensionality of the data to be enforced,
    371              commonly 4D for the data, and 3D for the mask in case of fMRI.
    372            scale_data : bool
    373              NIfTI header specifies scl_slope and scl_inter for scaling and
    374              offsetting the data.  By default those will get applied to the data
    375              (change in behavior post 0.4.6).
    376          """
    377          # check if we are in copy constructor mode
    378          if events is None:
    379              EventDataset.__init__(self, samples=samples, events=events,
    380                                    mask=mask, **kwargs)
    381              return
    382
    383          nifti = getNiftiFromAnySource(samples, ensure=True,
    384                                        enforce_dim=enforce_dim,
    385                                        scale_data=scale_data)
    386          # no copying
    387          samples = nifti.data
    388
    389          # do not put the whole NiftiImage in the dict as this will most
    390          # likely be deepcopy'ed at some point and ensuring data integrity
    391          # of the complex Python-C-Swig hybrid might be a tricky task.
    392          # Only storing the header dict should achieve the same and is more
    393          # memory efficient and even simpler
    394          dsattr = {'niftihdr': nifti.header}
    395
    396          # determine TR, take from NIfTI header by default
    397          dt = nifti.rtime
    398          # override if necessary
    399          if not tr is None:
    400              dt = tr
    401
    402          # NiftiDataset uses a DescreteMetric with cartesian
    403          # distance and element size from the NIfTI header
    404          # 'voxdim' is (x,y,z) while 'samples' are (t,z,y,x)
    405          elementsize = [dt] + [i for i in reversed(nifti.voxdim)]
    406          # XXX metric might be inappropriate if boxcar has length 1
    407          # might move metric setup after baseclass init and check what has
    408          # really happened
    409          metric = DescreteMetric(elementsize=elementsize,
    410                                  distance_function=cartesianDistance)
    411
    412          # convert EVs if necessary -- not altering original
    413          if evconv:
    414              if dt == 0:
    415                  raise ValueError, "'dt' cannot be zero when converting Events"
    416
    417              events = [ev.asDescreteTime(dt, storeoffset) for ev in events]
    418          else:
    419              # do not touch the original
    420              events = deepcopy(events)
    421
    422              # forcefully convert onset and duration into integers, as expected
    423              # by the baseclass
    424              for ev in events:
    425                  oldonset = ev['onset']
    426                  oldduration = ev['duration']
    427                  ev['onset'] = int(ev['onset'])
    428                  ev['duration'] = int(ev['duration'])
    429                  if not oldonset == ev['onset'] \
    430                     or not oldduration == ev['duration']:
    431                      warning("Loosing information during automatic integer "
    432                              "conversion of EVs. Consider an explicit conversion"
    433                              " by setting `evconv` in ERNiftiDataset().")
    434
    435          # pull mask array from NIfTI (if present)
    436          if mask is None:
    437              pass
    438          elif isinstance(mask, N.ndarray):
    439              # plain array can be passed on to base class
    440              pass
    441          else:
    442              mask_nim = getNiftiFromAnySource(mask, scale_data=scale_data)
    443              if not mask_nim is None:
    444                  mask = getNiftiData(mask_nim)
    445              else:
    446                  raise ValueError, "Cannot load mask from '%s'" % mask
    447
    448          # finally init baseclass
    449          EventDataset.__init__(self, samples=samples, events=events,
    450                                mask=mask, dametric=metric, dsattr=dsattr,
    451                                **kwargs)
    452
    453
    454 -    def map2Nifti(self, data=None):
    455          """Maps a data vector into the dataspace and wraps it with a
    456          NiftiImage. The header data of this object is used to initialize
    457          the new NiftiImage (scl_slope and scl_inter are reset to 1.0 and
    458          0.0 accordingly).
    459
    460          .. note::
    461            Only the features corresponding to voxels are mapped back -- not
    462            any additional features passed via the Event definitions.
    463
    464          :Parameters:
    465            data : ndarray or Dataset
    466              The data to be wrapped into NiftiImage. If None (default), it
    467              would wrap samples of the current dataset. If it is a Dataset
    468              instance -- takes its samples for mapping
    469          """
    470          if data is None:
    471              data = self.samples
    472          elif isinstance(data, Dataset):
    473              # ease users life
    474              data = data.samples
    475
    476          mr = self.mapper.reverse(data)
    477
    478          # trying to determine which part should go into NiftiImage
    479          if isinstance(self.mapper, CombinedMapper):
    480              # we have additional feature in the dataset -- ignore them
    481              mr = mr[0]
    482          else:
    483              pass
    484
    485          return NiftiImage(mr, _get_safe_header(self))
    486
    487
    488      niftihdr = property(fget=lambda self: self._dsattr['niftihdr'],
    489                          doc='Access to the NIfTI header dictionary.')
    490