mvpa.misc.plot.mri

Source Code for Module mvpa.misc.plot.mri

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 """Basic (f)MRI plotting with ability to interactively perform thresholding 10 11 """ 12 13 import pylab as P 14 import numpy as N 15 import matplotlib as mpl 16 17 from mvpa.base import warning, externals 18 19 if externals.exists('nifti', raiseException=True): 20 from nifti import NiftiImage 21 22 _interactive_backends = ['GTKAgg', 'TkAgg'] 23

24 -def plotMRI(background=None, background_mask=None, cmap_bg='gray', 25 overlay=None, overlay_mask=None, cmap_overlay='autumn', 26 vlim=(0.0, None), vlim_type=None, 27 do_stretch_colors=False, 28 add_info=True, add_hist=True, add_colorbar=True, 29 fig=None, interactive=None, 30 nrows=None, ncolumns=None 31 ):

32 """Very basic plotting of 3D data with interactive thresholding. 33 34 Background/overlay could be nifti files names or NiftiImage 35 objects, or 3D ndarrays. if no mask provided, only non-0 elements 36 are plotted 37 38 :Parameters: 39 do_stretch_colors : bool 40 Stratch color range to the data (not just to visible data) 41 vlim 42 2 element tuple of low/upper bounds of values to plot 43 vlim_type : None or 'symneg_z' 44 If not None, then vlim would be treated accordingly: 45 symneg_z 46 z-score values of symmetric normal around 0, estimated 47 by symmetrizing negative part of the distribution, which 48 often could be assumed when total distribution is a mixture of 49 by-chance performance normal around 0, and some other in the 50 positive tail 51 ncolumns : int or None 52 Explicit starting number of columns into which position the 53 slice renderings. 54 If None, square arrangement would be used 55 nrows : int or None 56 Explicit starting number of rows into which position the 57 slice renderings. 58 If None, square arrangement would be used 59 add_hist : bool or tuple (int, int) 60 If True, add histogram and position automagically. 61 If a tuple -- use as (row, column) 62 add_info : bool or tuple (int, int) 63 If True, add information and position automagically. 64 If a tuple -- use as (row, column). 65 66 Available colormaps are presented nicely on 67 http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps 68 69 TODO: 70 * Make interface more attractive/usable 71 * allow multiple overlays... or just unify for them all to be just a list of entries 72 * handle cases properly when there is only one - background/overlay 73 """ 74 # 75 if False: # for easy debugging 76 impath = '/research/fusion/herrman/be37/fMRI' 77 background = NiftiImage('%s/anat_slices_brain_inbold.nii.gz' % impath) 78 background_mask = None 79 overlay = NiftiImage('/research/fusion/herrman/code/CCe-1.nii.gz') 80 overlay_mask = NiftiImage('%s/masks/example_func_brain_mask.nii.gz' % impath) 81 82 do_stretch_colors = False 83 add_info = True 84 add_hist = True 85 add_colorbar = True 86 cmap_bg = 'gray' 87 cmap_overlay = 'hot' # YlOrRd_r # P.cm.autumn 88 89 fig = None 90 # vlim describes value limits 91 # clim color limits (same by default) 92 vlim = [2.3, None] 93 vlim_type = 'symneg_z' 94 interactive = False 95 96 # 97 # process data arguments 98 99 def handle_arg(arg): 100 """Helper which would read in NiftiImage if necessary 101 """ 102 if isinstance(arg, basestring): 103 arg = NiftiImage(arg) 104 argshape = arg.data.shape 105 # Assure that we have 3D (at least) 106 if len(argshape)<3: 107 arg.data = arg.data.reshape((1,)*(3-len(argshape)) + argshape) 108 if isinstance(arg, N.ndarray): 109 if len(arg.shape) != 3: 110 raise ValueError, "For now just handling 3D volumes" 111 return arg

112 113 bg = handle_arg(background) 114 if isinstance(bg, NiftiImage): 115 # figure out aspect 116 fov = (N.array(bg.header['pixdim']) * bg.header['dim'])[3:0:-1] 117 # XXX might be vise-verse ;-) 118 aspect = fov[2]/fov[1] 119 120 bg = bg.data[...,::-1,::-1] # XXX custom for now 121 else: 122 aspect = 1.0 123 124 if bg is not None: 125 bg_mask = handle_arg(background_mask) 126 if isinstance(bg_mask, NiftiImage): 127 bg_mask = bg_mask.data[...,::-1,::-1] # XXX 128 if bg_mask is not None: 129 bg_mask = bg_mask != 0 130 else: 131 bg_mask = bg != 0 132 133 func = handle_arg(overlay) 134 135 if func is not None: 136 if isinstance(func, NiftiImage): 137 func = func.data[..., ::-1, :] # XXX 138 139 func_mask = handle_arg(overlay_mask) 140 if isinstance(func_mask, NiftiImage): 141 func_mask = func_mask.data[..., ::-1, :] # XXX 142 if func_mask is not None: 143 func_mask = func_mask != 0 144 else: 145 func_mask = func != 0 146 147 148 # process vlim 149 vlim = list(vlim) 150 vlim_orig = vlim[:] 151 add_dist2hist = [] 152 if isinstance(vlim_type, basestring): 153 if vlim_type == 'symneg_z': 154 func_masked = func[func_mask] 155 fnonpos = func_masked[func_masked<=0] 156 fneg = func_masked[func_masked<0] 157 # take together with sign-reverted negative values 158 fsym = N.hstack((-fneg, fnonpos)) 159 nfsym = len(fsym) 160 # Estimate normal std under assumption of mean=0 161 std = N.sqrt(N.mean(abs(fsym)**2)) 162 # convert vlim assuming it is z-scores 163 for i,v in enumerate(vlim): 164 if v is not None: 165 vlim[i] = std * v 166 # add a plot to histogram 167 add_dist2hist = [(lambda x: nfsym/(N.sqrt(2*N.pi)*std)*N.exp(-(x**2)/(2*std**2)), 168 {})] 169 else: 170 raise ValueError, 'Unknown specification of vlim=%s' % vlim + \ 171 ' Known is: symneg' 172 173 174 class Plotter(object): 175 """ 176 TODO 177 """ 178 179 #_store_attribs = ('vlim', 'fig', 'bg', 'bg_mask') 180 181 def __init__(self, _locals): 182 """TODO""" 183 self._locals = _locals 184 self.fig = _locals['fig'] 185 186 def do_plot(self): 187 """TODO""" 188 # silly yarik didn't find proper way 189 vlim = self._locals['vlim'] 190 bg = self._locals['bg'] 191 bg_mask = self._locals['bg_mask'] 192 ncolumns = self._locals['ncolumns'] 193 nrows = self._locals['nrows'] 194 add_info = self._locals['add_info'] 195 add_hist = self._locals['add_hist'] 196 #print locals() 197 if N.isscalar(vlim): vlim = (vlim, None) 198 if vlim[0] is None: vlim = (N.min(func), vlim[1]) 199 if vlim[1] is None: vlim = (vlim[0], N.max(func)) 200 invert = vlim[1] < vlim[0] 201 if invert: 202 vlim = (vlim[1], vlim[0]) 203 print "Not yet fully supported" 204 205 # adjust lower bound if it is too low 206 if vlim[0] < N.min(func[func_mask]): 207 vlim = list(vlim) 208 vlim[0] = N.min(func[func_mask]) 209 vlim = tuple(vlim) 210 211 bound_above = (max(vlim) < N.max(func)) 212 bound_below = (min(vlim) > N.min(func)) 213 214 # 215 # reverse the map if needed 216 cmap_ = cmap_overlay 217 if not bound_below and bound_above: 218 if cmap_.endswith('_r'): 219 cmap_ = cmap_[:-2] 220 else: 221 cmap_ += '_r' 222 223 func_cmap = eval("P.cm.%s" % cmap_) 224 bg_cmap = eval("P.cm.%s" % cmap_bg) 225 226 227 if do_stretch_colors: 228 clim = (N.min(func), N.max(func))#vlim 229 else: 230 clim = vlim 231 232 # 233 # figure out 'extend' for colorbar and threshold string 234 extend, thresh_str = { 235 (True, True) : ('both', 'x in [%.3g, %.3g]' % tuple(vlim)), 236 (True, False): ('min', 'x in [%.3g, +inf]' % vlim[0]), 237 (False, True): ('max', 'x in (-inf, %.3g]' % vlim[1]), 238 (False, False): ('neither', 'none') }[(bound_below, 239 bound_above)] 240 241 # 242 # Figure out subplots 243 dshape = func.shape 244 nslices = func.shape[0] 245 246 # Check if additional column/row information was provided and extend nrows/ncolumns 247 for v in (add_hist, add_info): 248 if v and not isinstance(v, bool): 249 ncolumns = max(ncolumns, v[1]+1) 250 nrows = max(nrows, v[0]+1) 251 252 # more or less square alignment ;-) 253 if ncolumns is None: 254 ncolumns = int(N.sqrt(nslices)) 255 ndcolumns = ncolumns 256 nrows = max(nrows, int(N.ceil(nslices*1.0/ncolumns))) 257 258 259 # Decide either we need more cells where to add hist and/or info 260 nadd = bool(add_info) + bool(add_hist) 261 while ncolumns*nrows - (nslices + nadd) < 0: 262 ncolumns += 1 263 264 locs = ['' for i in xrange(ncolumns*nrows)] 265 266 # Fill in predefined locations 267 for v,vl in ((add_hist, 'hist'), 268 (add_info, 'info')): 269 if v and not isinstance(v, bool): 270 locs[ncolumns*v[0] + v[1]] = vl 271 272 # Fill in slices 273 for islice in xrange(nslices): 274 locs[locs.index('')] = islice 275 276 # Fill the last available if necessary 277 if add_hist and isinstance(add_hist, bool): 278 locs[locs.index('')] = 'hist' 279 if add_info and isinstance(add_info, bool): 280 locs[locs.index('')] = 'info' 281 282 print ncolumns, nrows 283 print locs 284 285 # should compare by backend? 286 if P.matplotlib.get_backend() in _interactive_backends: 287 P.ioff() 288 289 if self.fig is None: 290 self.fig = P.figure(facecolor='white', figsize=(4*ncolumns, 4*nrows)) 291 else: 292 self.fig.clf() 293 fig = self.fig 294 # fig.clf() 295 296 # 297 # how to threshold images 298 thresholder = lambda x: N.logical_and(x>=vlim[0], x<=vlim[1]) ^ invert 299 300 # 301 # Draw all slices 302 self.slices = [] 303 for si in range(nslices)[::-1]: 304 ax = fig.add_subplot(nrows, ncolumns, 305 locs.index(si) + 1, 306 frame_on=False) 307 self.slices.append(ax) 308 ax.axison = False 309 slice_bg = bg[si] 310 slice_bg_ = N.ma.masked_array(slice_bg, 311 mask=N.logical_not(bg_mask[si]))#slice_bg<=0) 312 313 slice_sl = func[si] 314 315 in_thresh = thresholder(slice_sl) 316 out_thresh = N.logical_not(in_thresh) 317 slice_sl_ = N.ma.masked_array(slice_sl, 318 mask=N.logical_or(out_thresh, 319 N.logical_not(func_mask[si]))) 320 321 kwargs = dict(aspect=aspect, origin='lower') 322 323 # paste a blank white background first, since otherwise 324 # recent matplotlib screws up those masked imshows 325 im = ax.imshow(N.ones(slice_sl_.shape), 326 cmap=bg_cmap, 327 extent=(0, slice_bg.shape[0], 328 0, slice_bg.shape[1]), 329 **kwargs) 330 im.set_clim((0,1)) 331 332 # ax.clim((0,1)) 333 ax.imshow(slice_bg_, 334 interpolation='bilinear', 335 cmap=bg_cmap, 336 **kwargs) 337 338 im = ax.imshow(slice_sl_, 339 interpolation='nearest', 340 cmap=func_cmap, 341 alpha=0.8, 342 extent=(0, slice_bg.shape[0], 343 0, slice_bg.shape[1]), 344 **kwargs) 345 im.set_clim(*clim) 346 347 if si == 0: 348 im0 = im 349 350 func_masked = func[func_mask] 351 352 # 353 # Add summary information 354 func_thr = func[N.logical_and(func_mask, thresholder(func))] 355 if add_info and len(func_thr): 356 self.info = ax = fig.add_subplot(nrows, ncolumns, 357 locs.index('info')+1, 358 frame_on=False) 359 # cb = P.colorbar(shrink=0.8) 360 # #cb.set_clim(clim[0], clim[1]) 361 ax.axison = False 362 #if add_colorbar: 363 # cb = P.colorbar(im, shrink=0.8, pad=0.0, drawedges=False, 364 # extend=extend, cmap=func_cmap) 365 366 stats = {'v':len(func_masked), 367 'vt': len(func_thr), 368 'm': N.mean(func_masked), 369 'mt': N.mean(func_thr), 370 'min': N.min(func_masked), 371 'mint': N.min(func_thr), 372 'max': N.max(func_masked), 373 'maxt': N.max(func_thr), 374 'mm': N.median(func_masked), 375 'mmt': N.median(func_thr), 376 'std': N.std(func_masked), 377 'stdt': N.std(func_thr), 378 'sthr': thresh_str} 379 P.text(0, 0.5, """ 380 Original: 381 voxels = %(v)d 382 range = [%(min).3g, %(max).3g] 383 mean = %(m).3g 384 median = %(mm).3g 385 std = %(std).3g 386 387 Thresholded: %(sthr)s: 388 voxels = %(vt)d 389 range = [%(mint).3g, %(maxt).3g] 390 median = %(mt).3g 391 mean = %(mmt).3g 392 std = %(stdt).3g 393 """ % stats, 394 horizontalalignment='left', 395 verticalalignment='center', 396 transform = ax.transAxes, 397 fontsize=14) 398 399 if add_colorbar: 400 kwargs_cb = {} 401 #if add_hist: 402 # kwargs_cb['cax'] = self.hist 403 self.cb = cb = P.colorbar( 404 im0, #self.hist, 405 shrink=0.8, pad=0.0, drawedges=False, 406 extend=extend, cmap=func_cmap, **kwargs_cb) 407 cb.set_clim(*clim) 408 409 # Add histogram 410 if add_hist: 411 self.hist = fig.add_subplot(nrows, ncolumns, 412 locs.index('hist') + 1, 413 frame_on=True) 414 415 minv, maxv = N.min(func_masked), N.max(func_masked) 416 if minv<0 and maxv>0: # then make it centered on 0 417 maxx = max(-minv, maxv) 418 range_ = (-maxx, maxx) 419 else: 420 range_ = (minv, maxv) 421 H = N.histogram(func_masked, range=range_, bins=31) 422 H2 = P.hist(func_masked, bins=H[1], align='center', 423 facecolor='#FFFFFF', hold=True) 424 for a, kwparams in add_dist2hist: 425 dbin = (H[1][1] - H[1][0]) 426 P.plot(H2[1], [a(x) * dbin for x in H2[1]], **kwparams) 427 if add_colorbar: 428 cbrgba = cb.to_rgba(H2[1]) 429 for face, facecolor, value in zip(H2[2], cbrgba, H2[1]): 430 if not thresholder(value): 431 color = '#FFFFFF' 432 else: 433 color = facecolor 434 face.set_facecolor(color) 435 436 437 fig.subplots_adjust(left=0.01, right=0.95, hspace=0.01) # , bottom=0.01 438 if ncolumns - int(bool(add_info) or bool(add_hist)) < 2: 439 fig.subplots_adjust(wspace=0.4) 440 else: 441 fig.subplots_adjust(wspace=0.1) 442 443 if P.matplotlib.get_backend() in _interactive_backends: 444 P.draw() 445 P.ion() 446 447 def on_click(self, event): 448 """Actions to perform on click 449 """ 450 if id(event.inaxes) != id(plotter.hist): 451 return 452 xdata, ydata, button = event.xdata, event.ydata, event.button 453 print xdata, ydata, button 454 vlim = self._locals['vlim'] 455 if button == 1: 456 vlim[0] = xdata 457 elif button == 3: 458 vlim[1] = xdata 459 elif button == 2: 460 vlim[0], vlim[1] = vlim[1], vlim[0] 461 self.do_plot() 462 463 plotter = Plotter(locals()) 464 plotter.do_plot() 465 466 if interactive is None: 467 interactive = P.matplotlib.get_backend() in _interactive_backends 468 469 # Global adjustments 470 if interactive: 471 # if P.matplotlib.is_interactive(): 472 P.connect('button_press_event', plotter.on_click) 473 P.show() 474 475 plotter.fig.plotter = plotter 476 return plotter.fig 477