Package mvpa :: Package atlases :: Module base
[hide private]
[frames] | no frames]

Source Code for Module mvpa.atlases.base

  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  """Base classes for Anatomy atlases support 
 10   
 11  TODOs: 
 12   * major optimization. Now code is sloppy and slow -- plenty of checks etc 
 13  """ 
 14   
 15  from mvpa.base import externals 
 16   
 17  externals.exists('lxml', raiseException=True) 
 18  from lxml import etree, objectify 
 19   
 20  import os, re 
 21  import numpy as N 
 22  from numpy.linalg import norm 
 23   
 24  from mvpa.atlases.transformation import SpaceTransformation, Linear 
 25  from mvpa.misc.support import reuseAbsolutePath 
 26   
 27  externals.exists('nifti', raiseException=True) 
 28  from nifti import NiftiImage 
 29   
 30  from mvpa.base import warning 
 31  if __debug__: 
 32          from mvpa.base import debug 
33 34 35 -def checkRange(coord, range):
36 """ 37 Check if coordinates are within range (0,0,0) - (range) 38 Return True on success 39 """ 40 # TODO: optimize 41 if len(coord) != len(range): 42 raise ValueError("Provided coordinate %s and given range %s" % \ 43 (`coord`, `range`) + \ 44 " have different dimensionality" 45 ) 46 for c,r in zip(coord, range): 47 if c<0 or c>=r: 48 return False 49 return True
50
51 52 -class BaseAtlas(object):
53 """Base class for the atlases. 54 """ 55
56 - def __init__ (self):
57 """ 58 Create an atlas object based on the... XXX 59 """ 60 self.__name = "blank" # XXX use or remove
61
62 63 -class XMLAtlasException(Exception):
64 """ Exception to be thrown if smth goes wrong dealing with XML based atlas 65 """
66 - def __init__(self, msg=""):
67 self.__msg = msg
68 - def __repr__(self):
69 return self.__msg
70
71 72 -class XMLBasedAtlas(BaseAtlas):
73
74 - def __init__(self, filename=None, resolution=None, query_voxel=False, 75 coordT=None, levels=None):
76 """ 77 :Parameters: 78 filename : string 79 Filename for the xml definition of the atlas 80 resolution : None or float 81 Some atlases link to multiple images at different 82 resolutions. if None -- best resolution is selected 83 using 0th dimension resolution 84 query_voxel : bool 85 By default [x,y,z] assumes coordinates in space, but if 86 query_voxel is True, they are assumed to be voxel coordinates 87 coordT 88 Optional transformation to apply first 89 levels : None or slice or list of int 90 What levels by default to operate on 91 """ 92 BaseAtlas.__init__(self) 93 self.__version = None 94 self.__type = None # XXX use or remove 95 self._imagefile = None 96 self.__atlas = None 97 self._filename = filename 98 self._resolution = resolution 99 self.query_voxel = query_voxel 100 self.levels = levels 101 102 if filename: 103 self.loadAtlas(filename) 104 105 # common sanity checks 106 if not self._checkVersion(self.version): 107 raise IOError("Version %s is not recognized to be native to class %s" % \ 108 (self.version, self.__name__)) 109 110 if not set(['header', 'data']) == set([i.tag for i in self.getchildren()]): 111 raise IOError("No header or data were defined in %s" % filename) 112 113 header = self.header 114 headerChildrenTags = XMLBasedAtlas._children_tags(header) 115 if not ('images' in headerChildrenTags) or \ 116 not ('imagefile' in XMLBasedAtlas._children_tags(header.images)): 117 raise XMLAtlasException("Atlas requires image/imagefile header fields") 118 119 # Load and post-process images 120 self._image = None 121 self._loadImages() 122 if self._image is not None: 123 self._extent = N.abs(N.asanyarray(self._image.extent[0:3])) 124 self._voxdim = N.asanyarray(self._image.voxdim) 125 self.relativeToOrigin = True 126 # Assign transformation to get into voxel coordinates, 127 # spaceT will be set accordingly 128 self.setCoordT(coordT) 129 self._loadData()
130 131
132 - def _checkRange(self, c):
133 """ check and adjust the voxel coordinates""" 134 # check range 135 # list(c) for consistent appearance... some times c might be ndarray 136 if __debug__: debug('ATL__', "Querying for voxel %s" % `list(c)`) 137 if not checkRange(c, self.extent): 138 msg = "Coordinates %s are not within the extent %s." \ 139 "Reset to (0,0,0)" % ( `c`, `self.extent` ) 140 if __debug__: debug('ATL_', msg) 141 # assume that voxel [0,0,0] is blank 142 c = [0]*3; 143 return c
144 145 146 @staticmethod
147 - def _checkVersion(version):
148 """To be overriden in the derived classes. By default anything is good""" 149 return True
150 151
152 - def _loadImages(self):
153 """To be overriden in the derived classes. By default does nothing""" 154 pass
155 156
157 - def _loadData(self):
158 """To be overriden in the derived classes. By default does nothing""" 159 pass
160 161
162 - def loadAtlas(self, filename):
163 if __debug__: debug('ATL_', "Loading atlas definition xml file " + filename) 164 # Create objectify parser first 165 parser = etree.XMLParser(remove_blank_text=True) 166 lookup = objectify.ObjectifyElementClassLookup() 167 parser.setElementClassLookup(lookup) 168 try: 169 self.__atlas = etree.parse(filename, parser).getroot() 170 except IOError: 171 raise XMLAtlasException("Failed to load XML file %s" % filename)
172 173 @property
174 - def version(self):
175 if not self.__atlas is None \ 176 and ("version" in self.__atlas.attrib.keys()): 177 return self.__atlas.get("version") 178 else: 179 return None
180 181 @staticmethod
182 - def _compare_lists(checkitems, neededitems):
183 raise RuntimeError, "DEPRECATED _compare_lists" 184 checkitems.sort() 185 neededitems.sort() 186 return (checkitems == neededitems)
187 188 189 @staticmethod
190 - def _children_tags(root):
191 return [i.tag for i in root.getchildren()]
192 193
194 - def __getattr__(self, attr):
195 """ 196 Lazy way to provide access to the definitions in the atlas 197 """ 198 if not self.__atlas is None: 199 return getattr(self.__atlas, attr) 200 else: 201 raise XMLAtlasException("Atlas in " + self.__name__ + " was not read yet")
202 203
204 - def setCoordT(self, coordT):
205 """Set coordT transformation. 206 207 spaceT needs to be adjusted since we glob those two 208 transformations together 209 """ 210 self._coordT = coordT # lets store for debugging etc 211 if self._image is not None: 212 # Combine with the image's qform 213 coordT = Linear(N.linalg.inv(self._image.qform), 214 previous=coordT) 215 self._spaceT = SpaceTransformation( 216 previous=coordT, toRealSpace=False 217 )
218 219
220 - def labelPoint(self, coord, levels=None):
221 """Return labels for the given spatial point at specified levels 222 223 so we first transform point into the voxel space 224 """ 225 coord_ = N.asarray(coord) # or we would alter what should be constant 226 #if not isinstance(coord, N.numpy): 227 #c = self.getVolumeCoordinate(coord) 228 #c = self.spaceT.toVoxelSpace(coord_) 229 #if self.coordT: 230 # coord_t = self.coordT[coord_] 231 #else: 232 # coord_t = coord_ 233 234 c = self.spaceT(coord_) 235 236 result = self.labelVoxel(c, levels) 237 result['coord_queried'] = coord 238 #result['coord_trans'] = coord_t 239 result['voxel_atlas'] = c 240 return result
241 242
243 - def levelsListing(self):
244 lkeys = range(self.Nlevels) 245 return '\n'.join(['%d: ' % k + str(self._levels_dict[k]) 246 for k in lkeys])
247 248
249 - def _getLevels(self, levels=None):
250 """Helper to provide list of levels to operate on 251 252 Depends on given `levels` as well as self.levels 253 """ 254 if levels is None: 255 levels = [ i for i in xrange(self.Nlevels) ] 256 elif (isinstance(levels, slice)): 257 # levels are given as a range 258 if levels.step: step = levels.step 259 else: step = 1 260 261 if levels.start: start = levels.start 262 else: start = 0 263 264 if levels.stop: stop = levels.stop 265 else: stop = self.Nlevels 266 267 levels = [ i for i in xrange(start, stop, step) ] 268 269 elif isinstance(levels, list) or isinstance(levels, tuple): 270 # levels given as list 271 levels = list(levels) 272 273 elif isinstance(levels, int): 274 levels = [ levels ] 275 276 else: 277 raise TypeError('Given levels "%s" are of unsupported type' % `levels`) 278 279 # test given values 280 levels_dict = self.levels_dict 281 for level in levels: 282 if not level in levels_dict: 283 raise ValueError, \ 284 "Levels %s is not known (out of range?). Known levels are:\n%s" \ 285 % (level, self.levelsListing()) 286 287 return levels
288 289
290 - def __getitem__(self, index):
291 """ 292 Accessing the elements via simple indexing. Examples: 293 print atlas[ 0, -7, 20, [1,2,3] ] 294 print atlas[ (0, -7, 20), 1:2 ] 295 print atlas[ (0, -7, 20) ] 296 print atlas[ (0, -7, 20), : ] 297 """ 298 if len(index) in [2, 4]: 299 levels_slice = index[-1] 300 else: 301 if self.levels is None: 302 levels_slice = slice(None,None,None) 303 else: 304 levels_slice = self.levels 305 306 levels = self._getLevels(levels=levels_slice) 307 308 if len(index) in [3, 4]: 309 # we got coordinates 1 by 1 + may be a level 310 coord = index[0:3] 311 312 elif len(index) in [1, 2]: 313 coord = index[0] 314 if isinstance(coord, list) or isinstance(coord, tuple): 315 if len(coord) != 3: 316 raise TypeError("Given coordinates must be in 3D") 317 else: 318 raise TypeError("Given coordinates must be a list or a tuple") 319 320 else: 321 raise TypeError("Unknown shape of parameters `%s`" % `index`) 322 323 if self.query_voxel: 324 return self.labelVoxel(coord, levels) 325 else: 326 return self.labelPoint(coord, levels)
327 328 329 # REDO in some sane fashion so referenceatlas returns levels for the base
330 - def _getLevelsDict(self):
331 return self._getLevelsDict_virtual()
332
333 - def _getLevelsDict_virtual(self):
334 return self._levels_dict
335 336 levels_dict = property(fget=_getLevelsDict) 337 338 339 origin = property(fget=lambda self:self._origin) 340 extent = property(fget=lambda self:self._extent) 341 voxdim = property(fget=lambda self:self._voxdim) 342 spaceT = property(fget=lambda self:self._spaceT) 343 coordT = property(fget=lambda self:self._spaceT, 344 fset=setCoordT)
345
346 -class Label(object):
347 """Represents a label. Just to bring all relevant information together 348 """
349 - def __init__ (self, text, abbr=None, coord=(None, None,None), 350 count=0, index=0):
351 """ 352 :Parameters: 353 text : basestring 354 fullname of the label 355 abbr : basestring 356 abbreviated name (optional) 357 coord : tuple of float 358 coordinates (optional) 359 count : int 360 count of those labels in the atlas (optional) 361 362 """ 363 self.__text = text.strip() 364 if abbr is not None: 365 abbr = abbr.strip() 366 self.__abbr = abbr 367 self.__coord = coord 368 self.__count = count 369 self.__index = int(index)
370 371 372 @property
373 - def index(self):
374 return self.__index
375
376 - def __repr__(self):
377 return "Label(%s%s, coord=(%s, %s, %s), count=%s, index=%s)" % \ 378 ((self.__text, 379 (', abbr=%s' % repr(self.__abbr), '')[int(self.__abbr is None)]) 380 + tuple(self.__coord) + (self.__count, self.__index))
381
382 - def __str__(self):
383 return self.__text
384 385 @staticmethod
386 - def generateFromXML(Elabel):
387 kwargs = {} 388 if Elabel.attrib.has_key('x'): 389 kwargs['coord'] = ( Elabel.attrib.get('x'), 390 Elabel.attrib.get('y'), 391 Elabel.attrib.get('z') ) 392 for l in ('count', 'abbr', 'index'): 393 if Elabel.attrib.has_key(l): 394 kwargs[l] = Elabel.attrib.get(l) 395 return Label(Elabel.text.strip(), **kwargs)
396 397 @property
398 - def count(self): return self.__count
399 @property
400 - def coord(self): return self.__coord
401 @property
402 - def text(self): return self.__text
403 @property
404 - def abbr(self):
405 """Returns abbreviated version if such is available 406 """ 407 if self.__abbr in [None, ""]: 408 return self.__text 409 else: 410 return self.__abbr
411
412 413 -class Level(object):
414 """Represents a level. Just to bring all relevant information together 415 """
416 - def __init__ (self, description):
417 self.description = description 418 self._type = "Base"
419
420 - def __repr__(self):
421 return "%s Level: %s" % \ 422 (self.levelType, self.description)
423
424 - def __str__(self):
425 return self.description
426 427 @staticmethod
428 - def generateFromXML(Elevel, levelType=None):
429 """ 430 Simple factory of levels 431 """ 432 if levelType is None: 433 if not Elevel.attrib.has_key("type"): 434 raise XMLAtlasException("Level must have type specified. Level: " + `Elevel`) 435 levelType = Elevel.get("type") 436 437 levelTypes = { 'label': LabelsLevel, 438 'reference': ReferencesLevel } 439 440 if levelTypes.has_key(levelType): 441 return levelTypes[levelType].generateFromXML(Elevel) 442 else: 443 raise XMLAtlasException("Unknown level type " + levelType)
444 445 levelType = property(lambda self: self._type)
446
447 448 -class LabelsLevel(Level):
449 """Level of labels. 450 451 XXX extend 452 """
453 - def __init__ (self, description, index=None, labels=[]):
454 Level.__init__(self, description) 455 self.__index = index 456 self.__labels = labels 457 self._type = "Labels"
458
459 - def __repr__(self):
460 return Level.__repr__(self) + " [%d] " % \ 461 (self.__index)
462 463 @staticmethod
464 - def generateFromXML(Elevel, levelIndex=[0]):
465 # XXX this is just for label type of level. For distance we need to ... 466 # we need to assure the right indexing 467 468 index = 0 469 if Elevel.attrib.has_key("index"): 470 index = int(Elevel.get("index")) 471 472 maxindex = max([int(i.get('index')) \ 473 for i in Elevel.label[:]]) 474 labels = [ None for i in xrange(maxindex+1) ] 475 for label in Elevel.label[:]: 476 labels[ int(label.get('index')) ] = Label.generateFromXML(label) 477 478 levelIndex[0] = max(levelIndex[0], index) + 1 # assign next one 479 480 return LabelsLevel(Elevel.get('description'), 481 index, 482 labels)
483 484 @property
485 - def index(self): return self.__index
486 487 @property
488 - def labels(self): return self.__labels
489
490 - def __getitem__(self, index):
491 return self.__labels[index]
492
493 494 -class ReferencesLevel(Level):
495 """Level which carries reference points 496 """
497 - def __init__ (self, description, indexes=[]):
498 Level.__init__(self, description) 499 self.__indexes = indexes 500 self._type = "References"
501 502 @staticmethod
503 - def generateFromXML(Elevel):
504 # XXX should probably do the same for the others? 505 requiredAttrs = ['x', 'y', 'z', 'type', 'description'] 506 if not set(requiredAttrs) == set(Elevel.attrib.keys()): 507 raise XMLAtlasException("ReferencesLevel has to have " + 508 "following attributes defined " + 509 `requiredAttrs`) 510 511 indexes = ( int(Elevel.get("x")), int(Elevel.get("y")), 512 int(Elevel.get("z")) ) 513 514 return ReferencesLevel(Elevel.get('description'), 515 indexes)
516 517 @property
518 - def indexes(self): return self.__indexes
519
520 521 -class PyMVPAAtlas(XMLBasedAtlas):
522 """Base class for PyMVPA atlases, such as LabelsAtlas and ReferenceAtlas 523 """ 524 525 source = 'PyMVPA' 526
527 - def __init__(self, *args, **kwargs):
528 XMLBasedAtlas.__init__(self, *args, **kwargs) 529 530 # sanity checks 531 header = self.header 532 headerChildrenTags = XMLBasedAtlas._children_tags(header) 533 if not ('space' in headerChildrenTags) or \ 534 not ('space-flavor' in headerChildrenTags): 535 raise XMLAtlasException("PyMVPA Atlas requires specification of" + 536 " the space in which atlas resides") 537 538 self.__space = header.space.text 539 self.__spaceFlavor = header['space-flavor'].text
540 541
542 - def _loadImages(self):
543 # shortcut 544 imagefile = self.header.images.imagefile 545 #self.Nlevels = len(self._levels_by_id) 546 547 # Set offset if defined in XML file 548 # XXX: should just take one from the qoffset... now that one is 549 # defined... this origin might be misleading actually 550 self._origin = N.array( (0,0,0) ) 551 if imagefile.attrib.has_key('offset'): 552 self._origin = N.array( map(int, 553 imagefile.get('offset').split(',')) ) 554 555 # Load the image file which has labels 556 imagefilename = reuseAbsolutePath(self._filename, imagefile.text) 557 558 try: 559 self._image = NiftiImage(imagefilename) 560 except RuntimeError, e: 561 raise RuntimeError, " Cannot open file " + imagefilename 562 563 self._data = self._image.data 564 565 # remove bogus dimensions on top of 4th 566 if len(self._data.shape[0:-4]) > 0: 567 bogus_dims = self._data.shape[0:-4] 568 if max(bogus_dims)>1: 569 raise RuntimeError, "Atlas %s has more than 4 of non-singular" \ 570 "dimensions" % imagefilename 571 new_shape = self._data.shape[-4:] 572 self._data.reshape(new_shape)
573 574 #if self._image.extent[3] != self.Nlevels: 575 # raise XMLAtlasException("Atlas %s has %d levels defined whenever %s has %d volumes" % \ 576 # ( filename, self.Nlevels, imagefilename, self._image.extent[3] )) 577 578
579 - def _loadData(self):
580 # Load levels 581 self._levels_dict = {} 582 # preprocess labels for different levels 583 self._Nlevels = 0 584 index_incr = 0 585 for index, child in enumerate(self.data.getchildren()): 586 if child.tag == 'level': 587 level = Level.generateFromXML(child) 588 self._levels_dict[level.description] = level 589 if hasattr(level, 'index'): 590 index = level.index 591 else: 592 # to avoid collision if some levels do 593 # have indexes 594 while index_incr in self._levels_dict: 595 index_incr += 1 596 index, index_incr = index_incr, index_incr+1 597 self._levels_dict[index] = level 598 else: 599 raise XMLAtlasException( 600 "Unknown child '%s' within data" % child.tag) 601 self._Nlevels += 1
602 603
604 - def _getNLevelsVirtual(self):
605 return self._Nlevels
606
607 - def _getNLevels(self):
608 return self._getNLevelsVirtual()
609 610 @staticmethod
611 - def _checkVersion(version):
612 # For compatibility lets support "RUMBA" atlases 613 return version.startswith("pymvpa-") or version.startswith("rumba-")
614 615 616 space = property(fget=lambda self:self.__space) 617 spaceFlavor = property(fget=lambda self:self.__spaceFlavor) 618 Nlevels = property(fget=_getNLevels)
619
620 621 -class LabelsAtlas(PyMVPAAtlas):
622 """ 623 Atlas which provides labels for the given coordinate 624 """ 625
626 - def labelVoxel(self, c, levels=None):
627 """ 628 Return labels for the given voxel at specified levels specified by index 629 """ 630 levels = self._getLevels(levels=levels) 631 632 result = {'voxel_queried' : c} 633 634 # check range 635 c = self._checkRange(c) 636 637 resultLevels = [] 638 for level in levels: 639 if self._levels_dict.has_key(level): 640 level_ = self._levels_dict[ level ] 641 else: 642 raise IndexError( 643 "Unknown index or description for level %d" % level) 644 645 resultIndex = int(self._data[ level_.index, \ 646 c[2], c[1], c[0] ]) 647 648 resultLevels += [ {'index': level_.index, 649 'id': level_.description, 650 'label' : level_[ resultIndex ]} ] 651 652 result['labels'] = resultLevels 653 return result
654
655 656 -class ReferencesAtlas(PyMVPAAtlas):
657 """ 658 Atlas which provides references to the other atlases. 659 660 Example: the atlas which has references to the closest points 661 (closest Gray, etc) in another atlas. 662 """ 663
664 - def __init__(self, distance=0, *args, **kwargs):
665 """Initialize `ReferencesAtlas` 666 """ 667 PyMVPAAtlas.__init__(self, *args, **kwargs) 668 # sanity checks 669 if not ('reference-atlas' in XMLBasedAtlas._children_tags(self.header)): 670 raise XMLAtlasException( 671 "ReferencesAtlas must refer to a some other atlas") 672 673 referenceAtlasName = self.header["reference-atlas"].text 674 675 # uff -- another evil import but we better use the factory method 676 from mvpa.atlases.warehouse import Atlas 677 self.__referenceAtlas = Atlas(filename=reuseAbsolutePath( 678 self._filename, referenceAtlasName)) 679 680 if self.__referenceAtlas.space != self.space or \ 681 self.__referenceAtlas.spaceFlavor != self.spaceFlavor: 682 raise XMLAtlasException( 683 "Reference and original atlases should be in the same space") 684 685 self.__referenceLevel = None 686 self.setDistance(distance)
687 688 689 # number of levels must be of the referenced atlas due to 690 # handling of that in __getitem__ 691 #Nlevels = property(fget=lambda self:self.__referenceAtlas.Nlevels)
692 - def _getNLevelsVirtual(self):
693 return self.__referenceAtlas.Nlevels
694 695
696 - def setReferenceLevel(self, level):
697 """ 698 Set the level which will be queried 699 """ 700 if self._levels_dict.has_key(level): 701 self.__referenceLevel = self._levels_dict[level] 702 else: 703 raise IndexError("Unknown reference level " + `level` + 704 ". Known are " + `self._levels_dict.keys()`)
705 706
707 - def labelVoxel(self, c, levels = None):
708 709 if self.__referenceLevel is None: 710 warning("You did not provide what level to use " 711 "for reference. Assigning 0th level -- '%s'" 712 % (self._levels_dict[0],)) 713 self.setReferenceLevel(0) 714 # return self.__referenceAtlas.labelVoxel(c, levels) 715 716 c = self._checkRange(c) 717 718 # obtain coordinates of the closest voxel 719 cref = self._data[ self.__referenceLevel.indexes, c[2], c[1], c[0] ] 720 dist = norm( (cref - c) * self.voxdim ) 721 if __debug__: 722 debug('ATL__', "Closest referenced point for %s is " 723 "%s at distance %3.2f" % (`c`, `cref`, dist)) 724 if (self.distance - dist) >= 1e-3: # neglect everything smaller 725 result = self.__referenceAtlas.labelVoxel(cref, levels) 726 result['voxel_referenced'] = c 727 result['distance'] = dist 728 else: 729 result = self.__referenceAtlas.labelVoxel(c, levels) 730 if __debug__: 731 debug('ATL__', "Closest referenced point is " 732 "further than desired distance %.2f" % self.distance) 733 result['voxel_referenced'] = None 734 result['distance'] = 0 735 return result
736 737
738 - def levelsListing(self):
739 return self.__referenceAtlas.levelsListing()
740
741 - def _getLevelsDict_virtual(self):
742 return self.__referenceAtlas.levels_dict
743
744 - def setDistance(self, distance):
745 """ 746 Set desired maximal distance for the reference 747 """ 748 if distance < 0: 749 raise ValueError("Distance should not be negative. " 750 " Thus '%f' is not a legal value" % distance) 751 if __debug__: 752 debug('ATL__', 753 "Setting maximal distance for queries to be %d" % distance) 754 self.__distance = distance
755 756 distance = property(fget=lambda self:self.__distance, fset=setDistance)
757