1
2
3
4
5
6
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
36 """
37 Check if coordinates are within range (0,0,0) - (range)
38 Return True on success
39 """
40
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
53 """Base class for the atlases.
54 """
55
57 """
58 Create an atlas object based on the... XXX
59 """
60 self.__name = "blank"
61
64 """ Exception to be thrown if smth goes wrong dealing with XML based atlas
65 """
70
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
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
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
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
127
128 self.setCoordT(coordT)
129 self._loadData()
130
131
133 """ check and adjust the voxel coordinates"""
134
135
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
142 c = [0]*3;
143 return c
144
145
146 @staticmethod
148 """To be overriden in the derived classes. By default anything is good"""
149 return True
150
151
153 """To be overriden in the derived classes. By default does nothing"""
154 pass
155
156
158 """To be overriden in the derived classes. By default does nothing"""
159 pass
160
161
163 if __debug__: debug('ATL_', "Loading atlas definition xml file " + filename)
164
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
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
183 raise RuntimeError, "DEPRECATED _compare_lists"
184 checkitems.sort()
185 neededitems.sort()
186 return (checkitems == neededitems)
187
188
189 @staticmethod
192
193
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
205 """Set coordT transformation.
206
207 spaceT needs to be adjusted since we glob those two
208 transformations together
209 """
210 self._coordT = coordT
211 if self._image is not None:
212
213 coordT = Linear(N.linalg.inv(self._image.qform),
214 previous=coordT)
215 self._spaceT = SpaceTransformation(
216 previous=coordT, toRealSpace=False
217 )
218
219
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)
226
227
228
229
230
231
232
233
234 c = self.spaceT(coord_)
235
236 result = self.labelVoxel(c, levels)
237 result['coord_queried'] = coord
238
239 result['voxel_atlas'] = c
240 return result
241
242
244 lkeys = range(self.Nlevels)
245 return '\n'.join(['%d: ' % k + str(self._levels_dict[k])
246 for k in lkeys])
247
248
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
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
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
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
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
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
332
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
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
375
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
384
385 @staticmethod
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
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
414 """Represents a level. Just to bring all relevant information together
415 """
417 self.description = description
418 self._type = "Base"
419
421 return "%s Level: %s" % \
422 (self.levelType, self.description)
423
425 return self.description
426
427 @staticmethod
444
445 levelType = property(lambda self: self._type)
446
449 """Level of labels.
450
451 XXX extend
452 """
453 - def __init__ (self, description, index=None, labels=[]):
458
460 return Level.__repr__(self) + " [%d] " % \
461 (self.__index)
462
463 @staticmethod
465
466
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
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
491 return self.__labels[index]
492
495 """Level which carries reference points
496 """
497 - def __init__ (self, description, indexes=[]):
501
502 @staticmethod
504
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
522 """Base class for PyMVPA atlases, such as LabelsAtlas and ReferenceAtlas
523 """
524
525 source = 'PyMVPA'
526
528 XMLBasedAtlas.__init__(self, *args, **kwargs)
529
530
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
543
544 imagefile = self.header.images.imagefile
545
546
547
548
549
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
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
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
575
576
577
578
580
581 self._levels_dict = {}
582
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
593
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
606
609
610 @staticmethod
612
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
622 """
623 Atlas which provides labels for the given coordinate
624 """
625
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
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
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):
687
688
689
690
691
693 return self.__referenceAtlas.Nlevels
694
695
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
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
715
716 c = self._checkRange(c)
717
718
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:
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
740
743
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