mvpa.atlases.fsl

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 """FSL atlases interfaces 10 11 """ 12 13 from mvpa.base import warning, externals 14 15 if externals.exists('nifti', raiseException=True): 16 from nifti import NiftiImage 17 18 import os, re 19 import numpy as N 20 21 from mvpa.misc.support import reuseAbsolutePath 22 from mvpa.base.dochelpers import enhancedDocString 23 24 from mvpa.atlases.base import XMLBasedAtlas, LabelsLevel 25 26 if __debug__: 27 from mvpa.base import debug

28 29 # 30 # Atlases from FSL 31 # 32 33 -class FSLAtlas(XMLBasedAtlas):

34 """Base class for FSL atlases 35 36 """ 37 source = 'FSL' 38 39

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

41 """ 42 """ 43 XMLBasedAtlas.__init__(self, *args, **kwargs) 44 self.space = 'MNI'

45 46 47 __doc__ = enhancedDocString('FSLAtlas', locals(), XMLBasedAtlas) 48 49

50 - def _loadImages(self):

51 resolution = self._resolution 52 header = self.header 53 images = header.images 54 # Load present images 55 # XXX might be refactored to avoid duplication of 56 # effort with PyMVPAAtlas 57 ni_image = None 58 resolutions = [] 59 if self._force_image_file is None: 60 imagefile_candidates = [ 61 reuseAbsolutePath(self._filename, i.imagefile.text, force=True) 62 for i in images] 63 else: 64 imagefile_candidates = [self._force_image_file] 65 66 for imagefilename in imagefile_candidates: 67 try: 68 ni_image_ = NiftiImage(imagefilename, load=False) 69 except RuntimeError, e: 70 raise RuntimeError, " Cannot open file " + imagefilename 71 72 resolution_ = ni_image_.pixdim[0] 73 if resolution is None: 74 # select this one if the best 75 if ni_image is None or \ 76 resolution_ < ni_image.pixdim[0]: 77 ni_image = ni_image_ 78 self._image_file = imagefilename 79 else: 80 if resolution_ == resolution: 81 ni_image = ni_image_ 82 self._image_file = imagefilename 83 break 84 else: 85 resolutions += [resolution_] 86 # TODO: also make use of summaryimagefile may be? 87 88 if ni_image is None: 89 msg = "Could not find an appropriate atlas among %d atlases." 90 if resolution is not None: 91 msg += " Atlases had resolutions %s" % \ 92 (resolutions,) 93 raise RuntimeError, msg 94 if __debug__: 95 debug('ATL__', "Loading atlas data from %s" % self._image_file) 96 self._image = ni_image 97 self._resolution = ni_image.pixdim[0] 98 self._origin = N.abs(ni_image.header['qoffset']) * 1.0 # XXX 99 self._data = self._image.data

100 101

102 - def _loadData(self):

103 """ """ 104 # Load levels 105 self._levels_dict = {} 106 # preprocess labels for different levels 107 self.Nlevels = 1 108 #level = Level.generateFromXML(self.data, levelType='label') 109 level = LabelsLevel.generateFromXML(self.data)#, levelType='label') 110 level.description = self.header.name.text 111 self._levels_dict = {0: level}

112 #for index, child in enumerate(self.data.getchildren()): 113 # if child.tag == 'level': 114 # level = Level.generateFromXML(child) 115 # self._levels_dict[level.description] = level 116 # try: 117 # self._levels_dict[level.index] = level 118 # except: 119 # pass 120 # else: 121 # raise XMLAtlasException("Unknown child '%s' within data" % child.tag) 122 # self.Nlevels += 1 123 #pass 124 125 126 @staticmethod

127 - def _checkVersion(version):

128 return re.search('^[0-9]+\.[0-9]', version) is not None

129

130 131 -class FSLProbabilisticAtlas(FSLAtlas):

132 """Probabilistic FSL atlases 133 """ 134

135 - def __init__(self, thr=0.0, strategy='all', sort=True, *args, **kwargs):

136 """ 137 138 :Parameters: 139 thr : float 140 Value to threshold at 141 strategy : basestring 142 Possible values 143 all - all entries above thr 144 max - entry with maximal value 145 sort : bool 146 Either to sort entries for 'all' strategy according to 147 probability 148 """ 149 150 FSLAtlas.__init__(self, *args, **kwargs) 151 self.thr = thr 152 self.strategy = strategy 153 self.sort = sort

154 155 __doc__ = enhancedDocString('FSLProbabilisticAtlas', locals(), FSLAtlas) 156

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

158 """Return labels for the voxel 159 160 :Parameters: 161 - c : tuple of coordinates (xyz) 162 - levels : just for API consistency (heh heh). Must be 0 for FSL atlases 163 """ 164 165 if levels is not None and not (levels in [0, [0], (0,)]): 166 raise ValueError, \ 167 "I guess we don't support levels other than 0 in FSL atlas" 168 169 # check range 170 c = self._checkRange(c) 171 172 # XXX think -- may be we should better assign each map to a 173 # different level 174 level = 0 175 resultLabels = [] 176 for index, area in enumerate(self._levels_dict[level]): 177 prob = int(self._data[index, c[2], c[1], c[0]]) 178 if prob > self.thr: 179 resultLabels += [dict(index=index, 180 #id= 181 label=area.text, 182 prob=prob)] 183 184 if self.sort or self.strategy == 'max': 185 resultLabels.sort(cmp=lambda x,y: cmp(x['prob'], y['prob']), 186 reverse=True) 187 188 if self.strategy == 'max': 189 resultLabels = resultLabels[:1] 190 elif self.strategy == 'all': 191 pass 192 else: 193 raise ValueError, 'Unknown strategy %s' % self.strategy 194 195 result = {'voxel_queried' : c, 196 # in the list since we have only single level but 197 # with multiple entries 198 'labels': [resultLabels]} 199 200 return result

201

202 - def find(self, *args, **kwargs):

203 """Just a shortcut to the only level. 204 205 See :class:`~mvpa.atlases.base.Level.find` for more info 206 """ 207 return self.levels_dict[0].find(*args, **kwargs)

208

209 - def getMap(self, target, strategy='unique'):

210 """Return a probability map 211 212 :Parameters: 213 target : int or str or re._pattern_type 214 If int, map for given index is returned. Otherwise, .find is called 215 with unique=True to find matching area 216 strategy : str in ('unique', 'max') 217 If 'unique', then if multiple areas match, exception would be raised. 218 In case of 'max', each voxel would get maximal value of probabilities 219 from all matching areas 220 """ 221 if isinstance(target, int): 222 return self._data[target] 223 else: 224 lev = self.levels_dict[0] # we have just 1 here 225 if strategy == 'unique': 226 return self.getMap(lev.find(target, unique=True).index) 227 else: 228 maps = N.array(self.getMaps(target)) 229 return N.max(maps, axis=0)

230

231 - def getMaps(self, target):

232 """Return a list of probability maps for the target 233 234 :Parameters: 235 target : str or re._pattern_type 236 .find is called with a target and unique=False to find all matches 237 """ 238 lev = self.levels_dict[0] # we have just 1 here 239 return [self.getMap(l.index) for l in lev.find(target, unique=False)]

240

241 242 -class FSLLabelsAtlas(XMLBasedAtlas):

243 """Not sure what this one was for"""

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

245 """not implemented""" 246 FSLAtlas.__init__(self, *args, **kwargs) 247 raise NotImplementedError

248

Source Code for Module mvpa.atlases.fsl