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 ====== 13 14 * major optimization. Now code is sloppy and slow -- plenty of checks etc 15 16 Module Organization 17 =================== 18 mvpa.atlases.base module contains support for various atlases 19 20 .. packagetree:: 21 :style: UML 22 23 :group Base: BaseAtlas XMLBasedAtlas Label Level LabelsLevel 24 :group Talairach: PyMVPAAtlas LabelsAtlas ReferencesAtlas 25 :group Exceptions: XMLAtlasException 26 27 """ 28 29 from mvpa.base import externals 30 31 externals.exists('lxml', raiseException=True) 32 from lxml import etree, objectify 33 34 import os, re 35 import numpy as N 36 from numpy.linalg import norm 37 38 from mvpa.atlases.transformation import SpaceTransformation, Linear 39 from mvpa.misc.support import reuseAbsolutePath 40 41 externals.exists('nifti', raiseException=True) 42 from nifti import NiftiImage 43 44 from mvpa.base import warning 45 if __debug__: 46 from mvpa.base import debug 47 48

49 -def checkRange(coord, range):

50 """ 51 Check if coordinates are within range (0,0,0) - (range) 52 Return True on success 53 """ 54 # TODO: optimize 55 if len(coord) != len(range): 56 raise ValueError("Provided coordinate %s and given range %s" % \ 57 (`coord`, `range`) + \ 58 " have different dimensionality" 59 ) 60 for c,r in zip(coord, range): 61 if c<0 or c>=r: 62 return False 63 return True

64 65

66 -class BaseAtlas(object):

67 """Base class for the atlases. 68 """ 69

70 - def __init__ (self):

71 """ 72 Create an atlas object based on the... XXX 73 """ 74 self.__name = "blank" # XXX use or remove

75 76

77 -class XMLAtlasException(Exception):

78 """ Exception to be thrown if smth goes wrong dealing with XML based atlas 79 """

80 - def __init__(self, msg=""):

81 self.__msg = msg

82 - def __repr__(self):

83 return self.__msg

84 85

86 -class XMLBasedAtlas(BaseAtlas):

87

88 - def __init__(self, filename=None, resolution=None, query_voxel=False, 89 coordT=None, levels=None):

90 """ 91 :Parameters: 92 filename : string 93 Filename for the xml definition of the atlas 94 resolution : None or float 95 Some atlases link to multiple images at different 96 resolutions. if None -- best resolution is selected 97 using 0th dimension resolution 98 query_voxel : bool 99 By default [x,y,z] assumes coordinates in space, but if 100 query_voxel is True, they are assumed to be voxel coordinates 101 coordT 102 Optional transformation to apply first 103 levels : None or slice or list of int 104 What levels by default to operate on 105 """ 106 BaseAtlas.__init__(self) 107 self.__version = None 108 self.__type = None # XXX use or remove 109 self._imagefile = None 110 self.__atlas = None 111 self._filename = filename 112 self._resolution = resolution 113 self.query_voxel = query_voxel 114 self.levels = levels 115 116 if filename: 117 self.loadAtlas(filename) 118 119 # common sanity checks 120 if not self._checkVersion(self.version): 121 raise IOError("Version %s is not recognized to be native to class %s" % \ 122 (self.version, self.__name__)) 123 124 if not set(['header', 'data']) == set([i.tag for i in self.getchildren()]): 125 raise IOError("No header or data were defined in %s" % filename) 126 127 header = self.header 128 headerChildrenTags = XMLBasedAtlas._children_tags(header) 129 if not ('images' in headerChildrenTags) or \ 130 not ('imagefile' in XMLBasedAtlas._children_tags(header.images)): 131 raise XMLAtlasException("Atlas requires image/imagefile header fields") 132 133 # Load and post-process images 134 self._image = None 135 self._loadImages() 136 if self._image is not None: 137 self._extent = N.abs(N.asanyarray(self._image.extent[0:3])) 138 self._voxdim = N.asanyarray(self._image.voxdim) 139 self.relativeToOrigin = True 140 # Assign transformation to get into voxel coordinates, 141 # spaceT will be set accordingly 142 self.setCoordT(coordT) 143 self._loadData()

144 145

146 - def _checkRange(self, c):

147 """ check and adjust the voxel coordinates""" 148 # check range 149 # list(c) for consistent appearance... some times c might be ndarray 150 if __debug__: debug('ATL__', "Querying for voxel %s" % `list(c)`) 151 if not checkRange(c, self.extent): 152 msg = "Coordinates %s are not within the extent %s." \ 153 "Reset to (0,0,0)" % ( `c`, `self.extent` ) 154 if __debug__: debug('ATL_', msg) 155 # assume that voxel [0,0,0] is blank 156 c = [0]*3; 157 return c

158 159 160 @staticmethod

161 - def _checkVersion(version):

162 """To be overriden in the derived classes. By default anything is good""" 163 return True

164 165

166 - def _loadImages(self):

167 """To be overriden in the derived classes. By default does nothing""" 168 pass

169 170

171 - def _loadData(self):

172 """To be overriden in the derived classes. By default does nothing""" 173 pass

174 175

176 - def loadAtlas(self, filename):

177 if __debug__: debug('ATL_', "Loading atlas definition xml file " + filename) 178 # Create objectify parser first 179 parser = etree.XMLParser(remove_blank_text=True) 180 lookup = objectify.ObjectifyElementClassLookup() 181 parser.setElementClassLookup(lookup) 182 try: 183 self.__atlas = etree.parse(filename, parser).getroot() 184 except IOError: 185 raise XMLAtlasException("Failed to load XML file %s" % filename)

186 187 @property

188 - def version(self):

189 if not self.__atlas is None \ 190 and ("version" in self.__atlas.attrib.keys()): 191 return self.__atlas.get("version") 192 else: 193 return None

194 195 @staticmethod

196 - def _compare_lists(checkitems, neededitems):

197 raise RuntimeError, "DEPRECATED _compare_lists" 198 checkitems.sort() 199 neededitems.sort() 200 return (checkitems == neededitems)

201 202 203 @staticmethod

204 - def _children_tags(root):

205 return [i.tag for i in root.getchildren()]

206 207

208 - def __getattr__(self, attr):

209 """ 210 Lazy way to provide access to the definitions in the atlas 211 """ 212 if not self.__atlas is None: 213 return getattr(self.__atlas, attr) 214 else: 215 raise XMLAtlasException("Atlas in " + self.__name__ + " was not read yet")

216 217

218 - def setCoordT(self, coordT):

219 """Set coordT transformation. 220 221 spaceT needs to be adjusted since we glob those two 222 transformations together 223 """ 224 self._coordT = coordT # lets store for debugging etc 225 if self._image is not None: 226 # Combine with the image's qform 227 coordT = Linear(N.linalg.inv(self._image.qform), 228 previous=coordT) 229 self._spaceT = SpaceTransformation( 230 previous=coordT, toRealSpace=False 231 )

232 233

234 - def labelPoint(self, coord, levels=None):

235 """Return labels for the given spatial point at specified levels 236 237 so we first transform point into the voxel space 238 """ 239 coord_ = N.asarray(coord) # or we would alter what should be constant 240 #if not isinstance(coord, N.numpy): 241 #c = self.getVolumeCoordinate(coord) 242 #c = self.spaceT.toVoxelSpace(coord_) 243 #if self.coordT: 244 # coord_t = self.coordT[coord_] 245 #else: 246 # coord_t = coord_ 247 248 c = self.spaceT(coord_) 249 250 result = self.labelVoxel(c, levels) 251 result['coord_queried'] = coord 252 #result['coord_trans'] = coord_t 253 result['voxel_atlas'] = c 254 return result

255 256

257 - def levelsListing(self):

258 lkeys = range(self.Nlevels) 259 return '\n'.join(['%d: ' % k + str(self._levels_dict[k]) 260 for k in lkeys])

261 262

263 - def _getLevels(self, levels=None):

264 """Helper to provide list of levels to operate on 265 266 Depends on given `levels` as well as self.levels 267 """ 268 if levels is None: 269 levels = [ i for i in xrange(self.Nlevels) ] 270 elif (isinstance(levels, slice)): 271 # levels are given as a range 272 if levels.step: step = levels.step 273 else: step = 1 274 275 if levels.start: start = levels.start 276 else: start = 0 277 278 if levels.stop: stop = levels.stop 279 else: stop = self.Nlevels 280 281 levels = [ i for i in xrange(start, stop, step) ] 282 283 elif isinstance(levels, list) or isinstance(levels, tuple): 284 # levels given as list 285 levels = list(levels) 286 287 elif isinstance(levels, int): 288 levels = [ levels ] 289 290 else: 291 raise TypeError('Given levels "%s" are of unsupported type' % `levels`) 292 293 # test given values 294 levels_dict = self.levels_dict 295 for level in levels: 296 if not level in levels_dict: 297 raise ValueError, \ 298 "Levels %s is not known (out of range?). Known levels are:\n%s" \ 299 % (level, self.levelsListing()) 300 301 return levels

302 303

304 - def __getitem__(self, index):

305 """ 306 Accessing the elements via simple indexing. Examples: 307 print atlas[ 0, -7, 20, [1,2,3] ] 308 print atlas[ (0, -7, 20), 1:2 ] 309 print atlas[ (0, -7, 20) ] 310 print atlas[ (0, -7, 20), : ] 311 """ 312 if len(index) in [2, 4]: 313 levels_slice = index[-1] 314 else: 315 if self.levels is None: 316 levels_slice = slice(None,None,None) 317 else: 318 levels_slice = self.levels 319 320 levels = self._getLevels(levels=levels_slice) 321 322 if len(index) in [3, 4]: 323 # we got coordinates 1 by 1 + may be a level 324 coord = index[0:3] 325 326 elif len(index) in [1, 2]: 327 coord = index[0] 328 if isinstance(coord, list) or isinstance(coord, tuple): 329 if len(coord) != 3: 330 raise TypeError("Given coordinates must be in 3D") 331 else: 332 raise TypeError("Given coordinates must be a list or a tuple") 333 334 else: 335 raise TypeError("Unknown shape of parameters `%s`" % `index`) 336 337 if self.query_voxel: 338 return self.labelVoxel(coord, levels) 339 else: 340 return self.labelPoint(coord, levels)

341 342 343 # REDO in some sane fashion so referenceatlas returns levels for the base

344 - def _getLevelsDict(self):

345 return self._getLevelsDict_virtual()

346

347 - def _getLevelsDict_virtual(self):

348 return self._levels_dict

349 350 levels_dict = property(fget=_getLevelsDict) 351 352 353 origin = property(fget=lambda self:self._origin) 354 extent = property(fget=lambda self:self._extent) 355 voxdim = property(fget=lambda self:self._voxdim) 356 spaceT = property(fget=lambda self:self._spaceT) 357 coordT = property(fget=lambda self:self._spaceT, 358 fset=setCoordT)

359

360 -class Label(object):

361 """Represents a label. Just to bring all relevant information together 362 """

363 - def __init__ (self, text, abbr=None, coord=(None, None,None), 364 count=0, index=0):

365 """ 366 :Parameters: 367 text : basestring 368 fullname of the label 369 abbr : basestring 370 abbreviated name (optional) 371 coord : tuple of float 372 coordinates (optional) 373 count : int 374 count of those labels in the atlas (optional) 375 376 """ 377 self.__text = text.strip() 378 if abbr is not None: 379 abbr = abbr.strip() 380 self.__abbr = abbr 381 self.__coord = coord 382 self.__count = count 383 self.__index = int(index)

384 385 386 @property

387 - def index(self):

388 return self.__index

389

390 - def __repr__(self):

391 return "Label(%s%s, coord=(%s, %s, %s), count=%s, index=%s)" % \ 392 ((self.__text, 393 (', abbr=%s' % repr(self.__abbr), '')[int(self.__abbr is None)]) 394 + tuple(self.__coord) + (self.__count, self.__index))

395

396 - def __str__(self):

397 return self.__text

398 399 @staticmethod

400 - def generateFromXML(Elabel):

401 kwargs = {} 402 if Elabel.attrib.has_key('x'): 403 kwargs['coord'] = ( Elabel.attrib.get('x'), 404 Elabel.attrib.get('y'), 405 Elabel.attrib.get('z') ) 406 for l in ('count', 'abbr', 'index'): 407 if Elabel.attrib.has_key(l): 408 kwargs[l] = Elabel.attrib.get(l) 409 return Label(Elabel.text.strip(), **kwargs)

410 411 @property

412 - def count(self): return self.__count

413 @property

414 - def coord(self): return self.__coord

415 @property

416 - def text(self): return self.__text

417 @property

418 - def abbr(self):

419 """Returns abbreviated version if such is available 420 """ 421 if self.__abbr in [None, ""]: 422 return self.__text 423 else: 424 return self.__abbr

425 426

427 -class Level(object):

428 """Represents a level. Just to bring all relevant information together 429 """

430 - def __init__ (self, description):

431 self.description = description 432 self._type = "Base"

433

434 - def __repr__(self):

435 return "%s Level: %s" % \ 436 (self.levelType, self.description)

437

438 - def __str__(self):

439 return self.description

440 441 @staticmethod

442 - def generateFromXML(Elevel, levelType=None):

443 """ 444 Simple factory of levels 445 """ 446 if levelType is None: 447 if not Elevel.attrib.has_key("type"): 448 raise XMLAtlasException("Level must have type specified. Level: " + `Elevel`) 449 levelType = Elevel.get("type") 450 451 levelTypes = { 'label': LabelsLevel, 452 'reference': ReferencesLevel } 453 454 if levelTypes.has_key(levelType): 455 return levelTypes[levelType].generateFromXML(Elevel) 456 else: 457 raise XMLAtlasException("Unknown level type " + levelType)

458 459 levelType = property(lambda self: self._type)

460 461

462 -class LabelsLevel(Level):

463 """Level of labels. 464 465 XXX extend 466 """

467 - def __init__ (self, description, index=None, labels=[]):

468 Level.__init__(self, description) 469 self.__index = index 470 self.__labels = labels 471 self._type = "Labels"

472

473 - def __repr__(self):

474 return Level.__repr__(self) + " [%d] " % \ 475 (self.__index)

476 477 @staticmethod

478 - def generateFromXML(Elevel, levelIndex=[0]):

479 # XXX this is just for label type of level. For distance we need to ... 480 # we need to assure the right indexing 481 482 index = 0 483 if Elevel.attrib.has_key("index"): 484 index = int(Elevel.get("index")) 485 486 maxindex = max([int(i.get('index')) \ 487 for i in Elevel.label[:]]) 488 labels = [ None for i in xrange(maxindex+1) ] 489 for label in Elevel.label[:]: 490 labels[ int(label.get('index')) ] = Label.generateFromXML(label) 491 492 levelIndex[0] = max(levelIndex[0], index) + 1 # assign next one 493 494 return LabelsLevel(Elevel.get('description'), 495 index, 496 labels)

497 498 @property

499 - def index(self): return self.__index

500 501 @property

502 - def labels(self): return self.__labels

503

504 - def __getitem__(self, index):

505 return self.__labels[index]

506 507

508 -class ReferencesLevel(Level):

509 """Level which carries reference points 510 """

511 - def __init__ (self, description, indexes=[]):

512 Level.__init__(self, description) 513 self.__indexes = indexes 514 self._type = "References"

515 516 @staticmethod

517 - def generateFromXML(Elevel):

518 # XXX should probably do the same for the others? 519 requiredAttrs = ['x', 'y', 'z', 'type', 'description'] 520 if not set(requiredAttrs) == set(Elevel.attrib.keys()): 521 raise XMLAtlasException("ReferencesLevel has to have " + 522 "following attributes defined " + 523 `requiredAttrs`) 524 525 indexes = ( int(Elevel.get("x")), int(Elevel.get("y")), 526 int(Elevel.get("z")) ) 527 528 return ReferencesLevel(Elevel.get('description'), 529 indexes)

530 531 @property

532 - def indexes(self): return self.__indexes

533 534

535 -class PyMVPAAtlas(XMLBasedAtlas):

536 """Base class for PyMVPA atlases, such as LabelsAtlas and ReferenceAtlas 537 """ 538 539 source = 'PyMVPA' 540

541 - def __init__(self, *args, **kwargs):

542 XMLBasedAtlas.__init__(self, *args, **kwargs) 543 544 # sanity checks 545 header = self.header 546 headerChildrenTags = XMLBasedAtlas._children_tags(header) 547 if not ('space' in headerChildrenTags) or \ 548 not ('space-flavor' in headerChildrenTags): 549 raise XMLAtlasException("PyMVPA Atlas requires specification of" + 550 " the space in which atlas resides") 551 552 self.__space = header.space.text 553 self.__spaceFlavor = header['space-flavor'].text

554 555

556 - def _loadImages(self):

557 # shortcut 558 imagefile = self.header.images.imagefile 559 #self.Nlevels = len(self._levels_by_id) 560 561 # Set offset if defined in XML file 562 # XXX: should just take one from the qoffset... now that one is 563 # defined... this origin might be misleading actually 564 self._origin = N.array( (0,0,0) ) 565 if imagefile.attrib.has_key('offset'): 566 self._origin = N.array( map(int, 567 imagefile.get('offset').split(',')) ) 568 569 # Load the image file which has labels 570 imagefilename = reuseAbsolutePath(self._filename, imagefile.text) 571 572 try: 573 self._image = NiftiImage(imagefilename) 574 except RuntimeError, e: 575 raise RuntimeError, " Cannot open file " + imagefilename 576 577 self._data = self._image.data 578 579 # remove bogus dimensions on top of 4th 580 if len(self._data.shape[0:-4]) > 0: 581 bogus_dims = self._data.shape[0:-4] 582 if max(bogus_dims)>1: 583 raise RuntimeError, "Atlas %s has more than 4 of non-singular" \ 584 "dimensions" % imagefilename 585 new_shape = self._data.shape[-4:] 586 self._data.reshape(new_shape)

587 588 #if self._image.extent[3] != self.Nlevels: 589 # raise XMLAtlasException("Atlas %s has %d levels defined whenever %s has %d volumes" % \ 590 # ( filename, self.Nlevels, imagefilename, self._image.extent[3] )) 591 592

593 - def _loadData(self):

594 # Load levels 595 self._levels_dict = {} 596 # preprocess labels for different levels 597 self._Nlevels = 0 598 index_incr = 0 599 for index, child in enumerate(self.data.getchildren()): 600 if child.tag == 'level': 601 level = Level.generateFromXML(child) 602 self._levels_dict[level.description] = level 603 if hasattr(level, 'index'): 604 index = level.index 605 else: 606 # to avoid collision if some levels do 607 # have indexes 608 while index_incr in self._levels_dict: 609 index_incr += 1 610 index, index_incr = index_incr, index_incr+1 611 self._levels_dict[index] = level 612 else: 613 raise XMLAtlasException( 614 "Unknown child '%s' within data" % child.tag) 615 self._Nlevels += 1

616 617

618 - def _getNLevelsVirtual(self):

619 return self._Nlevels

620

621 - def _getNLevels(self):

622 return self._getNLevelsVirtual()

623 624 @staticmethod

625 - def _checkVersion(version):

626 # For compatibility lets support "RUMBA" atlases 627 return version.startswith("pymvpa-") or version.startswith("rumba-")

628 629 630 space = property(fget=lambda self:self.__space) 631 spaceFlavor = property(fget=lambda self:self.__spaceFlavor) 632 Nlevels = property(fget=_getNLevels)

633 634

635 -class LabelsAtlas(PyMVPAAtlas):

636 """ 637 Atlas which provides labels for the given coordinate 638 """ 639

640 - def labelVoxel(self, c, levels=None):

641 """ 642 Return labels for the given voxel at specified levels specified by index 643 """ 644 levels = self._getLevels(levels=levels) 645 646 result = {'voxel_queried' : c} 647 648 # check range 649 c = self._checkRange(c) 650 651 resultLevels = [] 652 for level in levels: 653 if self._levels_dict.has_key(level): 654 level_ = self._levels_dict[ level ] 655 else: 656 raise IndexError( 657 "Unknown index or description for level %d" % level) 658 659 resultIndex = int(self._data[ level_.index, \ 660 c[2], c[1], c[0] ]) 661 662 resultLevels += [ {'index': level_.index, 663 'id': level_.description, 664 'label' : level_[ resultIndex ]} ] 665 666 result['labels'] = resultLevels 667 return result

668 669

670 -class ReferencesAtlas(PyMVPAAtlas):

671 """ 672 Atlas which provides references to the other atlases. 673 674 Example: the atlas which has references to the closest points 675 (closest Gray, etc) in another atlas. 676 """ 677

678 - def __init__(self, distance=0, *args, **kwargs):

679 """Initialize `ReferencesAtlas` 680 """ 681 PyMVPAAtlas.__init__(self, *args, **kwargs) 682 # sanity checks 683 if not ('reference-atlas' in XMLBasedAtlas._children_tags(self.header)): 684 raise XMLAtlasException( 685 "ReferencesAtlas must refer to a some other atlas") 686 687 referenceAtlasName = self.header["reference-atlas"].text 688 689 # uff -- another evil import but we better use the factory method 690 from mvpa.atlases.warehouse import Atlas 691 self.__referenceAtlas = Atlas(filename=reuseAbsolutePath( 692 self._filename, referenceAtlasName)) 693 694 if self.__referenceAtlas.space != self.space or \ 695 self.__referenceAtlas.spaceFlavor != self.spaceFlavor: 696 raise XMLAtlasException( 697 "Reference and original atlases should be in the same space") 698 699 self.__referenceLevel = None 700 self.setDistance(distance)

701 702 703 # number of levels must be of the referenced atlas due to 704 # handling of that in __getitem__ 705 #Nlevels = property(fget=lambda self:self.__referenceAtlas.Nlevels)

706 - def _getNLevelsVirtual(self):

707 return self.__referenceAtlas.Nlevels

708 709

710 - def setReferenceLevel(self, level):

711 """ 712 Set the level which will be queried 713 """ 714 if self._levels_dict.has_key(level): 715 self.__referenceLevel = self._levels_dict[level] 716 else: 717 raise IndexError("Unknown reference level " + `level` + 718 ". Known are " + `self._levels_dict.keys()`)

719 720

721 - def labelVoxel(self, c, levels = None):

722 723 if self.__referenceLevel is None: 724 warning("You did not provide what level to use " 725 "for reference. Assigning 0th level -- '%s'" 726 % (self._levels_dict[0],)) 727 self.setReferenceLevel(0) 728 # return self.__referenceAtlas.labelVoxel(c, levels) 729 730 c = self._checkRange(c) 731 732 # obtain coordinates of the closest voxel 733 cref = self._data[ self.__referenceLevel.indexes, c[2], c[1], c[0] ] 734 dist = norm( (cref - c) * self.voxdim ) 735 if __debug__: 736 debug('ATL__', "Closest referenced point for %s is " 737 "%s at distance %3.2f" % (`c`, `cref`, dist)) 738 if (self.distance - dist) >= 1e-3: # neglect everything smaller 739 result = self.__referenceAtlas.labelVoxel(cref, levels) 740 result['voxel_referenced'] = c 741 result['distance'] = dist 742 else: 743 result = self.__referenceAtlas.labelVoxel(c, levels) 744 if __debug__: 745 debug('ATL__', "Closest referenced point is " 746 "further than desired distance %.2f" % self.distance) 747 result['voxel_referenced'] = None 748 result['distance'] = 0 749 return result

750 751

752 - def levelsListing(self):

753 return self.__referenceAtlas.levelsListing()

754

755 - def _getLevelsDict_virtual(self):

756 return self.__referenceAtlas.levels_dict

757

758 - def setDistance(self, distance):

759 """ 760 Set desired maximal distance for the reference 761 """ 762 if distance < 0: 763 raise ValueError("Distance should not be negative. " 764 " Thus '%f' is not a legal value" % distance) 765 if __debug__: 766 debug('ATL__', 767 "Setting maximal distance for queries to be %d" % distance) 768 self.__distance = distance

769 770 distance = property(fget=lambda self:self.__distance, fset=setDistance)

771

Source Code for Module mvpa.atlases.base