Package mvpa :: Package atlases :: Module base

[hide private]
[frames] | no frames]

# Source Code for Module mvpa.atlases.base

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