Package mvpa :: Package misc :: Package plot :: Module topo

[hide private]
[frames] | no frames]

# Source Code for Module mvpa.misc.plot.topo

      1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil; encoding: utf-8 -*-
      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  #   The initial version of the code was contributed by Ingo Fründ and is
      9  #   Coypright (c) 2008 by Ingo Fründ ingo.fruend@googlemail.com
     10  #
     11  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
     12  """Plot parameter distributions on a head surface (topography plots)."""
     13
     14  __docformat__ = 'restructuredtext'
     15
     16  import numpy as N
     17
     18  from mvpa.base import externals
     19
     20  if externals.exists("pylab", raiseException=True):
     21      import pylab as P
     22      import numpy.ma as M
     23
     24  if externals.exists("griddata", raiseException=True):
     25      from mvpa.support.griddata import griddata
     26
     27  if externals.exists("scipy", raiseException=True):
     28      from scipy.optimize import leastsq
     29
     30
     31  # TODO : add optional plotting labels for the sensors
     32 -def plotHeadTopography(topography, sensorlocations, plotsensors=False,
     33                         resolution=51, masked=True, plothead=True,
     34                         plothead_kwargs=None, **kwargs):
     35      """Plot distribution to a head surface, derived from some sensor locations.
     36
     37      The sensor locations are first projected onto the best fitting sphere and
     38      finally projected onto a circle (by simply ignoring the z-axis).
     39
     40      :Parameters:
     41        topography: array
     42          A vector of some values corresponding to each sensor.
     43        sensorlocations: (nsensors x 3) array
     44          3D coordinates of each sensor. The order of the sensors has to match
     45          with the `topography` vector.
     46        plotsensors: bool
     47          If True, sensor will be plotted on their projected coordinates.
     48          No sensor are shown otherwise.
     49        plothead: bool
     50          If True, a head outline is plotted.
     51        plothead_kwargs: dict
     52          Additional keyword arguments passed to `plotHeadOutline()`.
     53        resolution: int
     54          Number of surface samples along both x and y-axis.
     55        masked: bool
     56          If True, all surface sample extending to head outline will be
     57          masked.
     58        **kwargs:
     59          All additional arguments will be passed to `pylab.imshow()`.
     60
     61
     62      :Returns:
     63        (map, head, sensors)
     64          The corresponding matplotlib objects are returned if plotted, ie.
     65          if plothead is set to `False`, `head` will be `None`.
     66
     67            map
     68              The colormap that makes the actual plot, a
     69              matplotlib.image.AxesImage instance.
     70            head
     71              What is returned by `plotHeadOutline()`.
     72            sensors
     73              The dots marking the electrodes, a matplotlib.lines.Line2d
     74              instance.
     75      """
     76      # give sane defaults
     77      if plothead_kwargs is None:
     78          plothead_kwargs = {}
     79
     80      # error function to fit the sensor locations to a sphere
     81      def err(params):
     82          r, cx, cy, cz = params
     83          return (sensorlocations[:, 0] - cx) ** 2 \
     84                 + (sensorlocations[:, 1] - cy) ** 2 \
     85                 + (sensorlocations[:, 2] - cz) ** 2 \
     86                 - r ** 2
     87
     88      # initial guess of sphere parameters (radius and center)
     89      params = (1, 0, 0, 0)
     90
     91      # do fit
     92      (r, cx, cy, cz), stuff = leastsq(err, params)
     93
     94      # size of each square
     95      ssh = float(r) / resolution         # half-size
     96      ss = ssh * 2.0                      # full-size
     97
     98      # Generate a grid and interpolate using the griddata module
     99      x = N.arange(cx - r, cx + r, ss) + ssh
    100      y = N.arange(cy - r, cy + r, ss) + ssh
    101      x, y = P.meshgrid(x, y)
    102
    103      # project the sensor locations onto the sphere
    104      sphere_center = N.array((cx, cy, cz))
    105      sproj = sensorlocations - sphere_center
    106      sproj = r * sproj / N.c_[N.sqrt(N.sum(sproj ** 2, axis=1))]
    107      sproj += sphere_center
    108
    109      # fit topology onto xy projection of sphere
    110      topo = griddata(sproj[:, 0], sproj[:, 1],
    111              N.ravel(N.array(topography)), x, y)
    112
    113      # mask values outside the head
    114      if masked:
    115          notinhead = N.greater_equal((x - cx) ** 2 + (y - cy) ** 2,
    116                                      (1.0 * r) ** 2)
    117          topo = M.masked_where(notinhead, topo)
    118
    119      # show surface
    120      map = P.imshow(topo, origin="lower", extent=(-r, r, -r, r), **kwargs)
    121      P.axis('off')
    122
    123      if plothead:
    124          # plot scaled head outline
    125          head = plotHeadOutline(scale=r, shift=(cx/2.0, cy/2.0), **plothead_kwargs)
    126      else:
    127          head = None
    128
    129      if plotsensors:
    130          # plot projected sensor locations
    131
    132          # reorder sensors so the ones below plotted first
    133          # TODO: please fix with more elegant solution
    134          zenum = [x[::-1] for x in enumerate(sproj[:, 2].tolist())]
    135          zenum.sort()
    136          indx = [ x[1] for x in zenum ]
    137          sensors = P.plot(sproj[indx, 0] - cx/2.0, sproj[indx, 1] - cy/2.0, 'wo')
    138      else:
    139          sensors = None
    140
    141      return map, head, sensors
    142
    143
    144 -def plotHeadOutline(scale=1, shift=(0, 0), color='k', linewidth='5', **kwargs):
    145      """Plots a simple outline of a head viewed from the top.
    146
    147      The plot contains schematic representations of the nose and ears. The
    148      size of the head is basically a unit circle for nose and ears attached
    149      to it.
    150
    151      :Parameters:
    152        scale: float
    153          Factor to scale the size of the head.
    154        shift: 2-tuple of floats
    155          Shift the center of the head circle by these values.
    156        color: matplotlib color spec
    157          The color the outline should be plotted in.
    158        linewidth: int
    159          Linewidth of the head outline.
    160        **kwargs:
    161          All additional arguments are passed to `pylab.plot()`.
    162
    163      :Returns:
    164        Matplotlib lines2D object
    165          can be used to tweak the look of the head outline.
    166      """
    167
    168      rmax = 0.5
    169      # factor used all the time
    170      fac = 2 * N.pi * 0.01
    171
    172      # Koordinates for the ears
    173      EarX1 =  -1 * N.array(
    174              [.497, .510, .518, .5299,
    175              .5419, .54, .547, .532, .510,
    176              rmax * N.cos(fac * (54 + 42))])
    177      EarY1 = N.array(
    178              [.0655, .0775, .0783, .0746, .0555,
    179              -.0055, -.0932, -.1313, -.1384,
    180              rmax * N.sin(fac * (54 + 42))])
    181      EarX2 = N.array(
    182              [rmax * N.cos(fac * (54 + 42)),
    183              .510, .532, .547, .54, .5419,
    184              .5299, .518, .510, .497] )
    185      EarY2 = N.array(
    186              [rmax * N.sin(fac * (54 + 42)),
    187              -.1384, -.1313, -.0932, -.0055,
    188              .0555, .0746, .0783, .0775, .0655] )
    189
    190      # Coordinates for the Head
    191      HeadX1 = N.fromfunction(
    192              lambda x: rmax * N.cos(fac * (x + 2)), (21,))
    193      HeadY1 = N.fromfunction(
    194              lambda y: rmax * N.sin(fac * (y + 2)), (21,))
    195      HeadX2 = N.fromfunction(
    196              lambda x: rmax * N.cos(fac * (x + 28)), (21,))
    197      HeadY2 = N.fromfunction(
    198              lambda y: rmax * N.sin(fac * (y + 28)), (21,))
    199      HeadX3 = N.fromfunction(
    200              lambda x: rmax * N.cos(fac * (x + 54)), (43,))
    201      HeadY3 = N.fromfunction(
    202              lambda y: rmax * N.sin(fac * (y + 54)), (43,))
    203
    204      # Coordinates for the Nose
    205      NoseX = N.array([.18 * rmax, 0, -.18 * rmax])
    206      NoseY = N.array([rmax - 0.004, rmax * 1.15, rmax - 0.004])
    207
    208      # Combine to one
    209      X = N.concatenate((EarX2, HeadX1, NoseX, HeadX2, EarX1, HeadX3))
    210      Y = N.concatenate((EarY2, HeadY1, NoseY, HeadY2, EarY1, HeadY3))
    211
    212      X *= 2 * scale
    213      Y *= 2 * scale
    214      X += shift[0]
    215      Y += shift[1]
    216
    217      return P.plot(X, Y, color=color, linewidth=linewidth)
    218