Package mvpa :: Package misc :: Package plot :: Module mri

[hide private]
[frames] | no frames]

# 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