mvpa.datasets.nifti

1 #emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- 2 #ex: set 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 externals.exists('nifti', raiseException=True) 15 16 import sys 17 import numpy as N 18 from mvpa.support.copy import deepcopy 19 20 if __debug__: 21 from mvpa.base import debug 22 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 name) 28 oldname = __name__ 29 # crazy name with close to zero possibility to cause whatever 30 __name__ = 'iaugf9zrkjsbdv89' 31 from nifti import NiftiImage 32 # restore old settings 33 __name__ = oldname 34 35 from mvpa.datasets.base import Dataset 36 from mvpa.datasets.mapped import MappedDataset 37 from mvpa.datasets.event import EventDataset 38 from mvpa.mappers.base import CombinedMapper 39 from mvpa.mappers.metric import DescreteMetric, cartesianDistance 40 from mvpa.mappers.array import DenseArrayMapper 41 from mvpa.base import warning 42 43

44 -def getNiftiFromAnySource(src, ensure=False, enforce_dim=None):

45 """Load/access NIfTI data from files or instances. 46 47 :Parameter: 48 src: str | NiftiImage 49 Filename of a NIfTI image or a `NiftiImage` instance. 50 ensure : bool 51 If True, through ValueError exception if cannot be loaded. 52 enforce_dim : int or None 53 If not None, it is the dimensionality of the data to be enforced, 54 commonly 4D for the data, and 3D for the mask in case of fMRI. 55 56 :Returns: 57 NiftiImage | None 58 If the source is not supported None is returned. 59 """ 60 nifti = None 61 62 # figure out what type 63 if isinstance(src, str): 64 # open the nifti file 65 try: 66 nifti = NiftiImage(src) 67 except RuntimeError, e: 68 warning("ERROR: NiftiDatasets: Cannot open NIfTI file %s" \ 69 % src) 70 raise e 71 elif isinstance(src, NiftiImage): 72 # nothing special 73 nifti = src 74 elif (isinstance(src, list) or isinstance(src, tuple)) \ 75 and len(src)>0 \ 76 and (isinstance(src[0], str) or isinstance(src[0], NiftiImage)): 77 # load from a list of given entries 78 if enforce_dim is not None: enforce_dim_ = enforce_dim - 1 79 else: enforce_dim_ = None 80 srcs = [getNiftiFromAnySource(s, ensure=ensure, 81 enforce_dim=enforce_dim_) 82 for s in src] 83 if __debug__: 84 # lets check if they all have the same dimensionality 85 shapes = [s.data.shape for s in srcs] 86 if not N.all([s == shapes[0] for s in shapes]): 87 raise ValueError, \ 88 "Input volumes contain variable number of dimensions:" \ 89 " %s" % (shapes,) 90 # Combine them all into a single beast 91 nifti = NiftiImage(N.array([s.asarray() for s in srcs]), 92 srcs[0].header) 93 elif ensure: 94 raise ValueError, "Cannot load NIfTI from %s" % (src,) 95 96 if nifti is not None and enforce_dim is not None: 97 shape, new_shape = nifti.data.shape, None 98 lshape = len(shape) 99 100 # check if we need to tune up shape 101 if lshape < enforce_dim: 102 # if we are missing required dimension(s) 103 new_shape = (1,)*(enforce_dim-lshape) + shape 104 elif lshape > enforce_dim: 105 # if there are bogus dimensions at the beginning 106 bogus_dims = lshape - enforce_dim 107 if shape[:bogus_dims] != (1,)*bogus_dims: 108 raise ValueError, \ 109 "Cannot enforce %dD on data with shape %s" \ 110 % (enforce_dim, shape) 111 new_shape = shape[bogus_dims:] 112 113 # tune up shape if needed 114 if new_shape is not None: 115 if __debug__: 116 debug('DS_NIFTI', 'Enforcing shape %s for %s data from %s' % 117 (new_shape, shape, src)) 118 nifti.data.shape = new_shape 119 120 return nifti

121 122

123 -def getNiftiData(nim):

124 """Convenience function to extract the data array from a NiftiImage 125 126 This function will make use of advanced features of PyNIfTI to prevent 127 unnecessary copying if a sufficent version is available. 128 """ 129 if externals.exists('nifti >= 0.20090205.1'): 130 return nim.data 131 else: 132 return nim.asarray()

133 134

135 -class NiftiDataset(MappedDataset):

136 """Dataset loading its samples from a NIfTI image or file. 137 138 Samples can be loaded from a NiftiImage instance or directly from a NIfTI 139 file. This class stores all relevant information from the NIfTI file header 140 and provides information about the metrics and neighborhood information of 141 all voxels. 142 143 Most importantly it allows to map data back into the original data space 144 and format via :meth:`~mvpa.datasets.nifti.NiftiDataset.map2Nifti`. 145 146 This class allows for convenient pre-selection of features by providing a 147 mask to the constructor. Only non-zero elements from this mask will be 148 considered as features. 149 150 NIfTI files are accessed via PyNIfTI. See 151 http://niftilib.sourceforge.net/pynifti/ for more information about 152 pynifti. 153 """ 154 # XXX: Every dataset should really have an example of howto instantiate 155 # it (necessary parameters).

156 - def __init__(self, samples=None, mask=None, dsattr=None, 157 enforce_dim=4, **kwargs):

158 """ 159 :Parameters: 160 samples: str | NiftiImage 161 Filename of a NIfTI image or a `NiftiImage` instance. 162 mask: str | NiftiImage | ndarray 163 Filename of a NIfTI image or a `NiftiImage` instance or an ndarray 164 of appropriate shape. 165 enforce_dim : int or None 166 If not None, it is the dimensionality of the data to be enforced, 167 commonly 4D for the data, and 3D for the mask in case of fMRI. 168 """ 169 # if in copy constructor mode 170 if not dsattr is None and dsattr.has_key('mapper'): 171 MappedDataset.__init__(self, 172 samples=samples, 173 dsattr=dsattr, 174 **kwargs) 175 return 176 177 # 178 # the following code only deals with contructing fresh datasets from 179 # scratch 180 # 181 182 # load the samples 183 niftisamples = getNiftiFromAnySource(samples, ensure=True, 184 enforce_dim=enforce_dim) 185 samples = niftisamples.data 186 187 # do not put the whole NiftiImage in the dict as this will most 188 # likely be deepcopy'ed at some point and ensuring data integrity 189 # of the complex Python-C-Swig hybrid might be a tricky task. 190 # Only storing the header dict should achieve the same and is more 191 # memory efficient and even simpler 192 dsattr = {'niftihdr': niftisamples.header} 193 194 195 # figure out what the mask is, but onyl handle known cases, the rest 196 # goes directly into the mapper which maybe knows more 197 niftimask = getNiftiFromAnySource(mask) 198 if niftimask is None: 199 pass 200 elif isinstance(niftimask, N.ndarray): 201 mask = niftimask 202 else: 203 mask = getNiftiData(niftimask) 204 205 # build an appropriate mapper that knows about the metrics of the NIfTI 206 # data 207 # NiftiDataset uses a DescreteMetric with cartesian 208 # distance and element size from the NIfTI header 209 210 # 'voxdim' is (x,y,z) while 'samples' are (t,z,y,x) 211 elementsize = [i for i in reversed(niftisamples.voxdim)] 212 mapper = DenseArrayMapper(mask=mask, shape=samples.shape[1:], 213 metric=DescreteMetric(elementsize=elementsize, 214 distance_function=cartesianDistance)) 215 216 MappedDataset.__init__(self, 217 samples=samples, 218 mapper=mapper, 219 dsattr=dsattr, 220 **kwargs)

221 222

223 - def map2Nifti(self, data=None):

224 """Maps a data vector into the dataspace and wraps it with a 225 NiftiImage. The header data of this object is used to initialize 226 the new NiftiImage. 227 228 :Parameters: 229 data : ndarray or Dataset 230 The data to be wrapped into NiftiImage. If None (default), it 231 would wrap samples of the current dataset. If it is a Dataset 232 instance -- takes its samples for mapping 233 """ 234 if data is None: 235 data = self.samples 236 elif isinstance(data, Dataset): 237 # ease users life 238 data = data.samples 239 dsarray = self.mapper.reverse(data) 240 return NiftiImage(dsarray, self.niftihdr)

241 242

243 - def getDt(self):

244 """Return the temporal distance of two samples/volumes. 245 246 This method tries to be clever and always returns `dt` in seconds, by 247 using unit information from the NIfTI header. If such information is 248 not present the assumed unit will also be `seconds`. 249 """ 250 # plain value 251 hdr = self.niftihdr 252 TR = hdr['pixdim'][4] 253 254 # by default assume seconds as unit and do not scale 255 scale = 1.0 256 257 # figure out units, if available 258 if hdr.has_key('time_unit'): 259 unit_code = hdr['time_unit'] 260 elif hdr.has_key('xyzt_unit'): 261 unit_code = int(hdr['xyzt_unit']) / 8 262 else: 263 warning("No information on time units is available. Assuming " 264 "seconds") 265 unit_code = 0 266 267 # handle known units 268 # XXX should be refactored to use actual unit labels from pynifti 269 # when version 0.20090205 or later is assumed to be available on all 270 # machines 271 if unit_code in [0, 1, 2, 3]: 272 if unit_code == 0: 273 warning("Time units were not specified in NiftiImage. " 274 "Assuming seconds.") 275 scale = [ 1.0, 1.0, 1e-3, 1e-6 ][unit_code] 276 else: 277 warning("Time units are incorrectly coded: value %d whenever " 278 "allowed are 8 (sec), 16 (millisec), 24 (microsec). " 279 "Assuming seconds.") 280 return TR * scale

281 282 283 niftihdr = property(fget=lambda self: self._dsattr['niftihdr'], 284 doc='Access to the NIfTI header dictionary.') 285 286 dt = property(fget=getDt, 287 doc='Time difference between two samples (in seconds). ' 288 'AKA TR in fMRI world.') 289 290 samplingrate = property(fget=lambda self: 1.0 / self.dt, 291 doc='Sampling rate (based on .dt).')

292 293

294 -class ERNiftiDataset(EventDataset):

295 """Dataset with event-defined samples from a NIfTI timeseries image. 296 297 This is a convenience dataset to facilitate the analysis of event-related 298 fMRI datasets. Boxcar-shaped samples are automatically extracted from the 299 full timeseries using :class:`~mvpa.misc.support.Event` definition lists. 300 For each event all volumes covering that particular event in time 301 (including partial coverage) are used to form the corresponding sample. 302 303 The class supports the conversion of events defined in 'realtime' into the 304 descrete temporal space defined by the NIfTI image. Moreover, potentially 305 varying offsets between true event onset and timepoint of the first selected 306 volume can be stored as an additional feature in the dataset. 307 308 Additionally, the dataset supports masking. This is done similar to the 309 masking capabilities of :class:`~mvpa.datasets.nifti.NiftiDataset`. However, 310 the mask can either be of the same shape as a single NIfTI volume, or 311 can be of the same shape as the generated boxcar samples, i.e. 312 a samples consisting of three volumes with 24 slices and 64x64 inplane 313 resolution needs a mask with shape (3, 24, 64, 64). In the former case the 314 mask volume is automatically expanded to be identical in a volumes of the 315 boxcar. 316 """

317 - def __init__(self, samples=None, events=None, mask=None, evconv=False, 318 storeoffset=False, tr=None, enforce_dim=4, **kwargs):

319 """ 320 :Paramaters: 321 mask: str | NiftiImage | ndarray 322 Filename of a NIfTI image or a `NiftiImage` instance or an ndarray 323 of appropriate shape. 324 evconv: bool 325 Convert event definitions using `onset` and `duration` in some 326 temporal unit into #sample notation. 327 storeoffset: Bool 328 Whether to store temproal offset information when converting 329 Events into descrete time. Only considered when evconv == True. 330 tr: float 331 Temporal distance of two adjacent NIfTI volumes. This can be used 332 to override the corresponding value in the NIfTI header. 333 enforce_dim : int or None 334 If not None, it is the dimensionality of the data to be enforced, 335 commonly 4D for the data, and 3D for the mask in case of fMRI. 336 """ 337 # check if we are in copy constructor mode 338 if events is None: 339 EventDataset.__init__(self, samples=samples, events=events, 340 mask=mask, **kwargs) 341 return 342 343 nifti = getNiftiFromAnySource(samples, ensure=True, 344 enforce_dim=enforce_dim) 345 # no copying 346 samples = nifti.data 347 348 # do not put the whole NiftiImage in the dict as this will most 349 # likely be deepcopy'ed at some point and ensuring data integrity 350 # of the complex Python-C-Swig hybrid might be a tricky task. 351 # Only storing the header dict should achieve the same and is more 352 # memory efficient and even simpler 353 dsattr = {'niftihdr': nifti.header} 354 355 # determine TR, take from NIfTI header by default 356 dt = nifti.rtime 357 # override if necessary 358 if not tr is None: 359 dt = tr 360 361 # NiftiDataset uses a DescreteMetric with cartesian 362 # distance and element size from the NIfTI header 363 # 'voxdim' is (x,y,z) while 'samples' are (t,z,y,x) 364 elementsize = [dt] + [i for i in reversed(nifti.voxdim)] 365 # XXX metric might be inappropriate if boxcar has length 1 366 # might move metric setup after baseclass init and check what has 367 # really happened 368 metric = DescreteMetric(elementsize=elementsize, 369 distance_function=cartesianDistance) 370 371 # convert EVs if necessary -- not altering original 372 if evconv: 373 if dt == 0: 374 raise ValueError, "'dt' cannot be zero when converting Events" 375 376 events = [ev.asDescreteTime(dt, storeoffset) for ev in events] 377 else: 378 # do not touch the original 379 events = deepcopy(events) 380 381 # forcefully convert onset and duration into integers, as expected 382 # by the baseclass 383 for ev in events: 384 oldonset = ev['onset'] 385 oldduration = ev['duration'] 386 ev['onset'] = int(ev['onset']) 387 ev['duration'] = int(ev['duration']) 388 if not oldonset == ev['onset'] \ 389 or not oldduration == ev['duration']: 390 warning("Loosing information during automatic integer " 391 "conversion of EVs. Consider an explicit conversion" 392 " by setting `evconv` in ERNiftiDataset().") 393 394 # pull mask array from NIfTI (if present) 395 if mask is None: 396 pass 397 elif isinstance(mask, N.ndarray): 398 # plain array can be passed on to base class 399 pass 400 else: 401 mask_nim = getNiftiFromAnySource(mask) 402 if not mask_nim is None: 403 mask = getNiftiData(mask_nim) 404 else: 405 raise ValueError, "Cannot load mask from '%s'" % mask 406 407 # finally init baseclass 408 EventDataset.__init__(self, samples=samples, events=events, 409 mask=mask, dametric=metric, dsattr=dsattr, 410 **kwargs)

411 412

413 - def map2Nifti(self, data=None):

414 """Maps a data vector into the dataspace and wraps it with a 415 NiftiImage. The header data of this object is used to initialize 416 the new NiftiImage. 417 418 .. note:: 419 Only the features corresponding to voxels are mapped back -- not 420 any additional features passed via the Event definitions. 421 422 :Parameters: 423 data : ndarray or Dataset 424 The data to be wrapped into NiftiImage. If None (default), it 425 would wrap samples of the current dataset. If it is a Dataset 426 instance -- takes its samples for mapping 427 """ 428 if data is None: 429 data = self.samples 430 elif isinstance(data, Dataset): 431 # ease users life 432 data = data.samples 433 434 mr = self.mapper.reverse(data) 435 436 # trying to determine which part should go into NiftiImage 437 if isinstance(self.mapper, CombinedMapper): 438 # we have additional feature in the dataset -- ignore them 439 mr = mr[0] 440 else: 441 pass 442 443 return NiftiImage(mr, self.niftihdr)

444 445 446 niftihdr = property(fget=lambda self: self._dsattr['niftihdr'], 447 doc='Access to the NIfTI header dictionary.')

448

Source Code for Module mvpa.datasets.nifti