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