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

Source Code for Module mvpa.datasets.nifti