diff --git a/fsl/fsleyes/controls/timeseriescontrolpanel.py b/fsl/fsleyes/controls/timeseriescontrolpanel.py
index c13d08a6fbedd89753fa51f6639f9c5658ca8e91..d8d71468c483edbacd0973c65dfa083e38ad3e98 100644
--- a/fsl/fsleyes/controls/timeseriescontrolpanel.py
+++ b/fsl/fsleyes/controls/timeseriescontrolpanel.py
@@ -11,12 +11,13 @@ control* which allows the user to configure a :class:`.TimeSeriesPanel`.
 
 import wx
 
-import                         props
-import pwidgets.widgetlist  as widgetlist
+import                                    props
+import pwidgets.widgetlist             as widgetlist
 
-import fsl.fsleyes.panel    as fslpanel
-import fsl.fsleyes.tooltips as fsltooltips
-import fsl.data.strings     as strings
+import fsl.fsleyes.panel               as fslpanel
+import fsl.fsleyes.plotting.timeseries as timeseries
+import fsl.fsleyes.tooltips            as fsltooltips
+import fsl.data.strings                as strings
 
 
 class TimeSeriesControlPanel(fslpanel.FSLEyesPanel):
@@ -262,8 +263,6 @@ class TimeSeriesControlPanel(fslpanel.FSLEyesPanel):
         # We're assuminbg that the TimeSeriesPanel has
         # already updated its current TimeSeries for
         # the newly selected overlay.
-        
-        import fsl.fsleyes.views.timeseriespanel as tsp
 
         if self.__selectedOverlay is not None:
             display = self._displayCtx.getDisplay(self.__selectedOverlay)
@@ -275,7 +274,7 @@ class TimeSeriesControlPanel(fslpanel.FSLEyesPanel):
 
         ts = self.__tsPanel.getCurrent()
 
-        if ts is None or not isinstance(ts, tsp.FEATTimeSeries):
+        if ts is None or not isinstance(ts, timeseries.FEATTimeSeries):
             return
 
         overlay = ts.overlay
diff --git a/fsl/fsleyes/controls/timeserieslistpanel.py b/fsl/fsleyes/controls/timeserieslistpanel.py
index a4da0a51259dfd9eebf70146c138b8794df768ab..fc9dca073a36dcbe0ca531109860589591e2c72f 100644
--- a/fsl/fsleyes/controls/timeserieslistpanel.py
+++ b/fsl/fsleyes/controls/timeserieslistpanel.py
@@ -15,12 +15,13 @@ import          copy
 import          wx
 import numpy as np
 
-import                                      props
-import pwidgets.elistbox                 as elistbox
-import fsl.fsleyes.panel                 as fslpanel
-import fsl.fsleyes.tooltips              as fsltooltips
-import fsl.data.strings                  as strings
-import fsl.fsleyes.colourmaps            as fslcm
+import                                    props
+import fsl.fsleyes.plotting.timeseries as timeseries
+import pwidgets.elistbox               as elistbox
+import fsl.fsleyes.panel               as fslpanel
+import fsl.fsleyes.tooltips            as fsltooltips
+import fsl.data.strings                as strings
+import fsl.fsleyes.colourmaps          as fslcm
 
 
 class TimeSeriesListPanel(fslpanel.FSLEyesPanel):
@@ -114,11 +115,9 @@ class TimeSeriesListPanel(fslpanel.FSLEyesPanel):
         """Creates a label to use for the given :class:`.TimeSeries` instance.
         """
 
-        import fsl.fsleyes.views.timeseriespanel as tsp
-
         display = self._displayCtx.getDisplay(ts.overlay)
 
-        if isinstance(ts, tsp.MelodicTimeSeries):
+        if isinstance(ts, timeseries.MelodicTimeSeries):
             return '{} [component {}]'.format(display.name, ts.coords)
         else:
             return '{} [{} {} {}]'.format(display.name,
diff --git a/fsl/fsleyes/plotting/__init__.py b/fsl/fsleyes/plotting/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/fsl/fsleyes/plotting/dataseries.py b/fsl/fsleyes/plotting/dataseries.py
new file mode 100644
index 0000000000000000000000000000000000000000..cb153735abe31143045716c857590c17f7528619
--- /dev/null
+++ b/fsl/fsleyes/plotting/dataseries.py
@@ -0,0 +1,74 @@
+#!/usr/bin/env python
+#
+# dataseries.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import props
+
+
+class DataSeries(props.HasProperties):
+    """A ``DataSeries`` instance encapsulates some data to be plotted by
+    a :class:`PlotPanel`, with the data extracted from an overlay in the
+    :class:`.OverlayList`. 
+
+    Sub-class implementations must accept an overlay object, pass this
+    overlay to the ``DataSeries`` constructor, and override the
+    :meth:`getData` method. The overlay is accessible as an instance
+    attribute, confusingly called ``overlay``.
+
+    Each``DataSeries`` instance is plotted as a line, with the line
+    style defined by properties on the ``DataSeries`` instance,
+    such as :attr:`colour`, :attr:`lineWidth` etc.
+    """
+
+    colour = props.Colour()
+    """Line colour. """
+
+    
+    alpha = props.Real(minval=0.0, maxval=1.0, default=1.0, clamped=True)
+    """Line transparency."""
+
+    
+    label = props.String()
+    """Line label (used in the plot legend)."""
+
+    
+    lineWidth = props.Choice((0.5, 1, 2, 3, 4, 5))
+    """Line width. """
+
+    
+    lineStyle = props.Choice(('-', '--', '-.', ':'))
+    """Line style. """
+
+    
+    def __init__(self, overlay):
+        """Create a ``DataSeries``.
+
+        :arg overlay: The overlay from which the data to be plotted is
+                      derived. 
+        """
+        
+        self.overlay = overlay
+
+
+    def __copy__(self):
+        """``DataSeries`` copy operator. Sub-classes with constructors
+        that require more than just the overlay object will need to
+        implement their own copy operator.
+        """
+        return type(self)(self.overlay)
+
+
+    def getData(self):
+        """This method must be implemented by sub-classes. It must return
+        the data to be plotted, as a tuple of the form:
+        
+            ``(xdata, ydata)``
+
+        where ``xdata`` and ``ydata`` are sequences containing the x/y data
+        to be plotted.
+        """
+        raise NotImplementedError('The getData method must be '
+                                  'implemented by subclasses')
diff --git a/fsl/fsleyes/plotting/histogramseries.py b/fsl/fsleyes/plotting/histogramseries.py
new file mode 100644
index 0000000000000000000000000000000000000000..ad45a9a016d3dead3858278bc453d7200e06793f
--- /dev/null
+++ b/fsl/fsleyes/plotting/histogramseries.py
@@ -0,0 +1,478 @@
+#!/usr/bin/env python
+#
+# histogramseries.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import logging
+
+import numpy as np
+
+import props
+
+import fsl.data.image as fslimage
+import                   dataseries
+
+
+log = logging.getLogger(__name__)
+
+
+class HistogramSeries(dataseries.DataSeries):
+    """A ``HistogramSeries`` generates histogram data from an :class:`.Image`
+    instance.
+    """
+
+    nbins = props.Int(minval=10, maxval=500, default=100, clamped=True)
+    """Number of bins to use in the histogram. This value is overridden
+    by the :attr:`HistogramPanel.autoBin` setting.
+
+    .. note:: I'm not sure why ``autoBin`` is a :class:`HistogramPanel`
+              setting, rather than a ``HistogramSeries`` setting. I might 
+              change this some time.
+    """
+
+    
+    ignoreZeros = props.Boolean(default=True)
+    """If ``True``, zeros are excluded from the calculated histogram. """
+
+    
+    showOverlay = props.Boolean(default=False)
+    """If ``True``, a 3D mask :class:`.Image` overlay is added to the
+    :class:`.OverlayList`, which highlights the voxels that have been included
+    in the histogram.
+    """
+
+
+    includeOutliers = props.Boolean(default=False)
+    """If ``True``, values which are outside of the :attr:`dataRange` are
+    included in the histogram end bins.
+    """
+
+    
+    volume = props.Int(minval=0, maxval=0, clamped=True)
+    """If the :class:`.Image` overlay associated with this ``HistogramSeries`` 
+    is 4D, this settings specifies the index of the volume that the histogram
+    is calculated upon.
+
+    .. note:: Calculating the histogram over an entire 4D :class:`.Image` is
+              not yet supported.
+    """
+
+    
+    dataRange = props.Bounds(ndims=1)
+    """Specifies the range of data which should be included in the histogram.
+    See the :attr:`includeOutliers` property.
+    """
+
+
+    def __init__(self,
+                 overlay,
+                 hsPanel,
+                 displayCtx,
+                 overlayList,
+                 volume=0,
+                 baseHs=None):
+        """Create a ``HistogramSeries``.
+
+        :arg overlay:     The :class:`.Image` overlay to calculate a histogram
+                          for.
+        
+        :arg hsPanel:     The :class:`HistogramPanel` that is displaying this
+                          ``HistogramSeries``.
+        
+        :arg displayCtx:  The :class:`.DisplayContext` instance.
+        
+        :arg overlayList: The :class:`.OverlayList` instance.
+        
+        :arg volume:      If the ``overlay`` is 4D, the initial value for the
+                          :attr:`volume` property.
+        
+        :arg baseHs:      If a ``HistogramSeries`` has already been created
+                          for the ``overlay``, it may be passed in here, so
+                          that the histogram data can be copied instead of
+                          having to be re-calculated.
+        """
+
+        log.debug('New HistogramSeries instance for {} '
+                  '(based on existing instance: {})'.format(
+                      overlay.name, baseHs is not None)) 
+
+        dataseries.DataSeries.__init__(self, overlay)
+
+        self.volume        = volume
+        self.__hsPanel     = hsPanel
+        self.__name        = '{}_{}'.format(type(self).__name__, id(self))
+
+        self.__displayCtx  = displayCtx
+        self.__overlayList = overlayList
+        self.__overlay3D   = None
+
+        if overlay.is4DImage():
+            self.setConstraint('volume', 'maxval', overlay.shape[3] - 1)
+
+        # If we have a baseHS, we 
+        # can copy all its data
+        if baseHs is not None:
+            self.dataRange.xmin     = baseHs.dataRange.xmin
+            self.dataRange.xmax     = baseHs.dataRange.xmax
+            self.dataRange.x        = baseHs.dataRange.x
+            self.nbins              = baseHs.nbins
+            self.volume             = baseHs.volume
+            self.ignoreZeros        = baseHs.ignoreZeros
+            self.includeOutliers    = baseHs.includeOutliers
+            
+            self.__nvals              =          baseHs.__nvals
+            self.__xdata              = np.array(baseHs.__xdata)
+            self.__ydata              = np.array(baseHs.__ydata)
+            self.__finiteData         = np.array(baseHs.__finiteData)
+            self.__nonZeroData        = np.array(baseHs.__nonZeroData)
+            self.__clippedFiniteData  = np.array(baseHs.__finiteData)
+            self.__clippedNonZeroData = np.array(baseHs.__nonZeroData) 
+
+        # Otherwise we need to calculate
+        # it all for ourselves
+        else:
+            self.__initProperties()
+        
+        overlayList.addListener('overlays',
+                                self.__name,
+                                self.__overlayListChanged)
+        self       .addListener('volume',
+                                self.__name,
+                                self.__volumeChanged)
+        self       .addListener('dataRange',
+                                self.__name,
+                                self.__dataRangeChanged)
+        self       .addListener('nbins',
+                                self.__name,
+                                self.__histPropsChanged)
+        self       .addListener('ignoreZeros',
+                                self.__name,
+                                self.__histPropsChanged)
+        self       .addListener('includeOutliers',
+                                self.__name,
+                                self.__histPropsChanged)
+        self       .addListener('showOverlay',
+                                self.__name,
+                                self.__showOverlayChanged)
+
+        
+    def update(self):
+        """This method may be called to force re-calculation of the
+        histogram data.
+        """
+        self.__histPropsChanged()
+
+        
+    def destroy(self):
+        """This needs to be called when this ``HistogramSeries`` instance
+        is no longer being used.
+
+        It removes several property listeners and, if the :attr:`overlay3D`
+        property is ``True``, removes the mask overlay from the
+        :class:`.OverlayList`.
+        """
+        self              .removeListener('nbins',           self.__name)
+        self              .removeListener('ignoreZeros',     self.__name)
+        self              .removeListener('includeOutliers', self.__name)
+        self              .removeListener('volume',          self.__name)
+        self              .removeListener('dataRange',       self.__name)
+        self              .removeListener('nbins',           self.__name)
+        self.__overlayList.removeListener('overlays',        self.__name)
+        
+        if self.__overlay3D is not None:
+            self.__overlayList.remove(self.__overlay3D)
+            self.__overlay3D = None
+
+        
+    def __initProperties(self):
+        """Called by :meth:`__init__`. Calculates and caches some things which
+        are needed for the histogram calculation.
+
+        .. note:: This method is never called if a ``baseHs`` is provided to
+                 :meth:`__init__`.
+        """
+
+        log.debug('Performining initial histogram '
+                  'calculations for overlay {}'.format(
+                      self.overlay.name))
+
+        data  = self.overlay.data[:]
+        
+        finData = data[np.isfinite(data)]
+        dmin    = finData.min()
+        dmax    = finData.max()
+        dist    = (dmax - dmin) / 10000.0
+        
+        nzData = finData[finData != 0]
+        nzmin  = nzData.min()
+        nzmax  = nzData.max()
+
+        self.dataRange.xmin = dmin
+        self.dataRange.xmax = dmax  + dist
+        self.dataRange.xlo  = nzmin
+        self.dataRange.xhi  = nzmax + dist
+
+        self.nbins = self.__autoBin(nzData, self.dataRange.x)
+
+        if not self.overlay.is4DImage():
+            
+            self.__finiteData  = finData
+            self.__nonZeroData = nzData
+            self.__dataRangeChanged(callHistPropsChanged=False)
+            
+        else:
+            self.__volumeChanged(callHistPropsChanged=False)
+
+        self.__histPropsChanged()
+
+        
+    def __volumeChanged(
+            self,
+            ctx=None,
+            value=None,
+            valid=None,
+            name=None,
+            callHistPropsChanged=True):
+        """Called when the :attr:`volume` property changes, and also by the
+        :meth:`__initProperties` method.
+
+        Re-calculates some things for the new overlay volume.
+
+        :arg callHistPropsChanged: If ``True`` (the default), the
+                                   :meth:`__histPropsChanged` method will be
+                                   called.
+
+        All other arguments are ignored, but are passed in when this method is
+        called due to a property change (see the
+        :meth:`.HasProperties.addListener` method).        
+        """
+
+        if self.overlay.is4DImage(): data = self.overlay.data[..., self.volume]
+        else:                        data = self.overlay.data[:]
+
+        data = data[np.isfinite(data)]
+
+        self.__finiteData  = data
+        self.__nonZeroData = data[data != 0]
+
+        self.__dataRangeChanged(callHistPropsChanged=False)
+
+        if callHistPropsChanged:
+            self.__histPropsChanged()
+
+
+    def __dataRangeChanged(
+            self,
+            ctx=None,
+            value=None,
+            valid=None,
+            name=None,
+            callHistPropsChanged=True):
+        """Called when the :attr:`dataRange` property changes, and also by the
+        :meth:`__initProperties` and :meth:`__volumeChanged` methods.
+
+        :arg callHistPropsChanged: If ``True`` (the default), the
+                                   :meth:`__histPropsChanged` method will be
+                                   called.
+
+        All other arguments are ignored, but are passed in when this method is
+        called due to a property change (see the
+        :meth:`.HasProperties.addListener` method).
+        """
+        finData = self.__finiteData
+        nzData  = self.__nonZeroData
+        
+        self.__clippedFiniteData  = finData[(finData >= self.dataRange.xlo) &
+                                            (finData <  self.dataRange.xhi)]
+        self.__clippedNonZeroData = nzData[ (nzData  >= self.dataRange.xlo) &
+                                            (nzData  <  self.dataRange.xhi)]
+
+        if callHistPropsChanged:
+            self.__histPropsChanged()
+
+            
+    def __histPropsChanged(self, *a):
+        """Called internally, and when any histogram settings change.
+        Re-calculates the histogram data.
+        """
+
+        log.debug('Calculating histogram for '
+                  'overlay {}'.format(self.overlay.name))
+
+        if self.dataRange.xhi - self.dataRange.xlo < 0.00000001:
+            self.__xdata = np.array([])
+            self.__ydata = np.array([])
+            self.__nvals = 0
+            return
+
+        if self.ignoreZeros:
+            if self.includeOutliers: data = self.__nonZeroData
+            else:                    data = self.__clippedNonZeroData
+        else:
+            if self.includeOutliers: data = self.__finiteData
+            else:                    data = self.__clippedFiniteData 
+        
+        if self.__hsPanel.autoBin:
+            nbins = self.__autoBin(data, self.dataRange.x)
+
+            self.disableListener('nbins', self.__name)
+            self.nbins = nbins
+            self.enableListener('nbins', self.__name)
+
+        # Calculate bin edges
+        bins = np.linspace(self.dataRange.xlo,
+                           self.dataRange.xhi,
+                           self.nbins + 1)
+
+        if self.includeOutliers:
+            bins[ 0] = self.dataRange.xmin
+            bins[-1] = self.dataRange.xmax
+            
+        # Calculate the histogram
+        histX    = bins
+        histY, _ = np.histogram(data.flat, bins=bins)
+            
+        self.__xdata = histX
+        self.__ydata = histY
+        self.__nvals = histY.sum()
+
+        log.debug('Calculated histogram for overlay '
+                  '{} (number of values: {}, number '
+                  'of bins: {})'.format(
+                      self.overlay.name,
+                      self.__nvals,
+                      self.nbins))
+
+
+    def __showOverlayChanged(self, *a):
+        """Called when the :attr:`showOverlay` property changes.
+
+        Adds/removes a 3D mask :class:`.Image` to the :class:`.OverlayList`,
+        which highlights the voxels that have been included in the histogram.
+        The :class:`.MaskOpts.threshold` property is bound to the
+        :attr:`dataRange` property, so the masked voxels are updated whenever
+        the histogram data range changes, and vice versa.
+        """
+
+        if not self.showOverlay:
+            if self.__overlay3D is not None:
+
+                log.debug('Removing 3D histogram overlay mask for {}'.format(
+                    self.overlay.name))
+                self.__overlayList.remove(self.__overlay3D)
+                self.__overlay3D = None
+
+        else:
+
+            log.debug('Creating 3D histogram overlay mask for {}'.format(
+                self.overlay.name))
+            
+            self.__overlay3D = fslimage.Image(
+                self.overlay.data,
+                name='{}/histogram/mask'.format(self.overlay.name),
+                header=self.overlay.nibImage.get_header())
+
+            self.__overlayList.append(self.__overlay3D)
+
+            opts = self.__displayCtx.getOpts(self.__overlay3D,
+                                             overlayType='mask')
+
+            opts.bindProps('volume',    self)
+            opts.bindProps('colour',    self)
+            opts.bindProps('threshold', self, 'dataRange')
+
+
+    def __overlayListChanged(self, *a):
+        """Called when the :class:`.OverlayList` changes.
+
+        If a 3D mask overlay was being shown, and it has been removed from the
+        ``OverlayList``, the :attr:`showOverlay` property is updated
+        accordingly.
+        """
+        
+        if self.__overlay3D is None:
+            return
+
+        # If a 3D overlay was being shown, and it
+        # has been removed from the overlay list
+        # by the user, turn the showOverlay property
+        # off
+        if self.__overlay3D not in self.__overlayList:
+            
+            self.disableListener('showOverlay', self.__name)
+            self.showOverlay = False
+            self.__showOverlayChanged()
+            self.enableListener('showOverlay', self.__name)
+
+        
+    def __autoBin(self, data, dataRange):
+        """Calculates the number of bins which should be used for a histogram
+        of the given data. The calculation is identical to that implemented
+        in the original FSLView.
+
+        :arg data:      The data that the histogram is to be calculated on.
+
+        :arg dataRange: A tuple containing the ``(min, max)`` histogram range.
+        """
+
+        dMin, dMax = dataRange
+        dRange     = dMax - dMin
+
+        binSize = np.power(10, np.ceil(np.log10(dRange) - 1) - 1)
+
+        nbins = dRange / binSize
+
+        while nbins < 100:
+            binSize /= 2
+            nbins    = dRange / binSize
+
+        if issubclass(data.dtype.type, np.integer):
+            binSize = max(1, np.ceil(binSize))
+
+        adjMin = np.floor(dMin / binSize) * binSize
+        adjMax = np.ceil( dMax / binSize) * binSize
+
+        nbins = int((adjMax - adjMin) / binSize) + 1
+
+        return nbins
+            
+
+    def getData(self):
+        """Overrides :meth:`.DataSeries.getData`.
+
+        Returns  a tuple containing the ``(x, y)`` histogram data.
+        """
+
+        if len(self.__xdata) == 0 or \
+           len(self.__ydata) == 0:
+            return self.__xdata, self.__ydata
+
+        # If smoothing is not enabled, we'll
+        # munge the histogram data a bit so
+        # that plt.plot(drawstyle='steps-pre')
+        # plots it nicely.
+        if not self.__hsPanel.smooth:
+
+            xdata = np.zeros(len(self.__xdata) + 1, dtype=np.float32)
+            ydata = np.zeros(len(self.__ydata) + 2, dtype=np.float32)
+
+            xdata[ :-1] = self.__xdata
+            xdata[  -1] = self.__xdata[-1]
+            ydata[1:-1] = self.__ydata
+
+
+        # If smoothing is enabled, the above munge
+        # is not necessary, and will probably cause
+        # the spline interpolation (performed by 
+        # the PlotPanel) to fail.
+        else:
+            xdata = np.array(self.__xdata[:-1], dtype=np.float32)
+            ydata = np.array(self.__ydata,      dtype=np.float32)
+
+        nvals    = self.__nvals
+        histType = self.__hsPanel.histType
+            
+        if   histType == 'count':       return xdata, ydata
+        elif histType == 'probability': return xdata, ydata / nvals
diff --git a/fsl/fsleyes/plotting/timeseries.py b/fsl/fsleyes/plotting/timeseries.py
new file mode 100644
index 0000000000000000000000000000000000000000..548bd34a2a78ebb5909d3a02f0e10b6d58b76029
--- /dev/null
+++ b/fsl/fsleyes/plotting/timeseries.py
@@ -0,0 +1,632 @@
+#!/usr/bin/env python
+#
+# timeseries.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import logging
+
+import numpy as np
+
+import props
+
+import dataseries
+
+
+
+log = logging.getLogger(__name__)
+
+
+
+class TimeSeries(dataseries.DataSeries):
+    """Encapsulates time series data from a specific voxel in an
+    :class:`.Image` overlay.
+    """
+
+    
+    def __init__(self, tsPanel, overlay, coords):
+        """Create a ``TimeSeries`` instane.
+
+        :arg tsPanel: The :class:`TimeSeriesPanel` which owns this
+                      ``TimeSeries``.
+
+        :arg overlay: The :class:`.Image` instance to extract the data from.
+
+        :arg coords:  The voxel coordinates of the time series.
+        """
+        dataseries.DataSeries.__init__(self, overlay)
+
+        self.tsPanel = tsPanel
+        self.coords  = None
+        self.data    = None
+        self.update(coords)
+
+
+    def __copy__(self):
+        """Copy operator, creates and returns a copy of this ``TimeSeries``
+        instance.
+        """
+        return type(self)(self.tsPanel, self.overlay, self.coords)
+
+        
+    def update(self, coords):
+        """Updates the voxel coordinates and time series data encapsulated
+        by this ``TimeSeries``.
+
+        .. warning:: This method is only intended for use on the *current*
+                     ``TimeSeriesPanel`` series, not for time series instances
+                     which have been added to the :attr:`.PlotPanel.dataSeries`
+                     list.
+        
+                     It is used by the :class:`TimeSeriesPanel` so that
+                     ``TimeSeries`` instances do not have to constantly be
+                     destroyed/recreated whenever the
+                     :attr:`.DisplayContext.location` changes.
+        """
+
+        self.coords = coords
+        self.data = self._getData(coords)
+        return True
+
+
+    def _getData(self, coords):
+        return self.overlay.data[coords[0], coords[1], coords[2], :]
+    
+        
+    def getData(self, xdata=None, ydata=None):
+        """Overrides :meth:`.DataSeries.getData` Returns the data associated
+        with this ``TimeSeries`` instance.
+
+        The ``xdata`` and ``ydata`` arguments may be used by subclasses to
+        override the x/y data in the event that they have already performed
+        some processing on the data.
+        """
+
+        if xdata is None: xdata = np.arange(len(self.data), dtype=np.float32)
+        if ydata is None: ydata = np.array(     self.data,  dtype=np.float32)
+
+        if self.tsPanel.usePixdim:
+            xdata *= self.overlay.pixdim[3]
+        
+        if self.tsPanel.plotMode == 'demean':
+            ydata = ydata - ydata.mean()
+
+        elif self.tsPanel.plotMode == 'normalise':
+            ymin  = ydata.min()
+            ymax  = ydata.max()
+            ydata = 2 * (ydata - ymin) / (ymax - ymin) - 1
+            
+        elif self.tsPanel.plotMode == 'percentChange':
+            mean  = ydata.mean()
+            ydata =  100 * (ydata / mean) - 100
+            
+        return xdata, ydata
+
+ 
+class FEATTimeSeries(TimeSeries):
+    """A :class:`TimeSeries` class for use with :class:`FEATImage` instances,
+    containing some extra FEAT specific options.
+
+    The ``FEATTimeSeries`` class acts as a container for several
+    ``TimeSeries`` instances, each of which represent some part of a FEAT
+    analysis. Therefore, the data returned by a call to
+    :meth:`.TimeSeries.getData` on a ``FEATTimeSeries`` instance should not
+    be plotted.
+
+    Instead, the :meth:`getModelTimeSeries` method should be used to retrieve
+    a list of all the ``TimeSeries`` instances which are associated with the
+    ``FEATTimeSeries`` instance - all of these ``TimeSeries`` instances should
+    be plotted instead.
+
+    For example, if the :attr:`plotData` and :attr:`plotFullModelFit` settings
+    are ``True``, the :meth:`getModelTimeSeries` method will return a list
+    containing two ``TimeSeries`` instances - one which will return the FEAT
+    analysis input data, and another which will return the full model fit, for
+    the voxel in question.
+
+
+    .. note:: The ``getData`` method of a ``FEATTimeSeries`` instance will
+              return the FEAT input data; therefore, when :attr:`plotData` is
+              ``True``, the ``FEATTimeSeries`` instance will itself be included
+              in the list returned by :meth:`getModelTimeSeries`.
+
+    
+    The following classes are used to represent the various parts of a FEAT
+    analysis:
+
+    .. autosummary::
+       :nosignatures:
+    
+       FEATEVTimeSeries
+       FEATResidualTimeSeries
+       FEATPartialFitTimeSeries
+       FEATModelFitTimeSeries
+    """
+
+
+    plotData = props.Boolean(default=True)
+    """If ``True``, the FEAT input data is plotted. """
+
+    
+    plotFullModelFit = props.Boolean(default=True)
+    """If ``True``, the FEAT full model fit is plotted. """
+
+    
+    plotResiduals = props.Boolean(default=False)
+    """If ``True``, the FEAT model residuals are plotted. """
+
+    
+    plotEVs = props.List(props.Boolean(default=False))
+    """A list of ``Boolean`` properties, one for each EV in the FEAT analysis.
+    For elements that are ``True``, the corresponding FEAT EV time course is
+    plotted.
+    """
+
+    
+    plotPEFits = props.List(props.Boolean(default=False))
+    """A list of ``Boolean`` properties, one for each EV in the FEAT analysis.
+    For elements that are ``True``, the model fit for the corresponding FEAT
+    EV is plotted.    
+    """
+
+    
+    plotCOPEFits = props.List(props.Boolean(default=False))
+    """A list of ``Boolean`` properties, one for each EV in the FEAT analysis.
+    For elements that are ``True``, the model fit for the corresponding FEAT
+    contrast is plotted.    
+    """ 
+
+    
+    plotPartial = props.Choice()
+    """Plot the raw data, after regression against a chosen EV or contrast.
+    The options are populated in the :meth:`__init__` method.
+    """
+    
+
+    def __init__(self, *args, **kwargs):
+        """Create a ``FEATTimeSeries``.
+
+        All arguments are passed through to the :class:`TimeSeries`
+        constructor.
+        """
+        
+        TimeSeries.__init__(self, *args, **kwargs)
+        self.name = '{}_{}'.format(type(self).__name__, id(self))
+
+        numEVs    = self.overlay.numEVs()
+        numCOPEs  = self.overlay.numContrasts()
+        copeNames = self.overlay.contrastNames()
+        
+        reduceOpts = ['none'] + \
+                     ['PE{}'.format(i + 1) for i in range(numEVs)]
+
+        for i in range(numCOPEs):
+            name = 'COPE{} ({})'.format(i + 1, copeNames[i])
+            reduceOpts.append(name)
+
+        self.getProp('plotPartial').setChoices(reduceOpts, instance=self)
+
+        for i in range(numEVs):
+            self.plotPEFits.append(False)
+            self.plotEVs   .append(False)
+
+        for i in range(numCOPEs):
+            self.plotCOPEFits.append(False) 
+
+        self.__fullModelTs =  None
+        self.__partialTs   =  None
+        self.__resTs       =  None
+        self.__evTs        = [None] * numEVs
+        self.__peTs        = [None] * numEVs
+        self.__copeTs      = [None] * numCOPEs
+        
+        self.addListener('plotFullModelFit',
+                         self.name,
+                         self.__plotFullModelFitChanged)
+        self.addListener('plotResiduals',
+                         self.name,
+                         self.__plotResidualsChanged)
+        self.addListener('plotPartial',
+                         self.name,
+                         self.__plotPartialChanged)
+
+        self.addListener('plotEVs',      self.name, self.__plotEVChanged)
+        self.addListener('plotPEFits',   self.name, self.__plotPEFitChanged)
+        self.addListener('plotCOPEFits', self.name, self.__plotCOPEFitChanged)
+
+        # plotFullModelFit defaults to True, so
+        # force the model fit ts creation here
+        self.__plotFullModelFitChanged()
+
+
+    def __copy__(self):
+        """Copy operator for a ``FEATTimeSeries`` instance."""
+        
+        copy = type(self)(self.tsPanel, self.overlay, self.coords)
+
+        copy.colour           = self.colour
+        copy.alpha            = self.alpha 
+        copy.label            = self.label 
+        copy.lineWidth        = self.lineWidth
+        copy.lineStyle        = self.lineStyle
+
+        # When these properties are changed 
+        # on the copy instance, it will create 
+        # its own FEATModelFitTimeSeries 
+        # instances accordingly
+        copy.plotFullModelFit = self.plotFullModelFit
+        copy.plotEVs[     :]  = self.plotEVs[     :]
+        copy.plotPEFits[  :]  = self.plotPEFits[  :]
+        copy.plotCOPEFits[:]  = self.plotCOPEFits[:]
+        copy.plotPartial      = self.plotPartial
+        copy.plotResiduals    = self.plotResiduals
+
+        return copy
+ 
+
+    def getModelTimeSeries(self):
+        """Returns a list containing all of the ``TimeSeries`` instances
+        which should be plotted in place of this ``FEATTimeSeries``.
+        """
+        
+        
+        modelts = []
+
+        if self.plotData:              modelts.append(self)
+        if self.plotFullModelFit:      modelts.append(self.__fullModelTs)
+        if self.plotResiduals:         modelts.append(self.__resTs)
+        if self.plotPartial != 'none': modelts.append(self.__partialTs)
+        
+        for i in range(self.overlay.numEVs()):
+            if self.plotPEFits[i]:
+                modelts.append(self.__peTs[i])
+
+        for i in range(self.overlay.numEVs()):
+            if self.plotEVs[i]:
+                modelts.append(self.__evTs[i]) 
+
+        for i in range(self.overlay.numContrasts()):
+            if self.plotCOPEFits[i]:
+                modelts.append(self.__copeTs[i])
+
+        return modelts
+
+
+    def update(self, coords):
+        """Overrides :meth:`TimeSeries.update`.
+
+        Updates the coordinates and data associated wsith this
+        ``FEATTimeSeries`` instance.
+        """
+        
+        if not TimeSeries.update(self, coords):
+            return False
+            
+        for modelTs in self.getModelTimeSeries():
+            if modelTs is self:
+                continue
+            modelTs.update(coords)
+
+        return True
+
+
+    def __getContrast(self, fitType, idx):
+        """Returns a contrast vector for the given model fit type, and index.
+
+        :arg fitType: either ``'full'``, ``'pe'``, or ``'cope'``. If
+                      ``'full'``, the ``idx`` argument is ignored.
+
+        :arg idx:     The EV or contrast index for ``'pe'`` or ``'cope'`` model
+                      fits.
+        """
+
+        if fitType == 'full':
+            return [1] * self.overlay.numEVs()
+        elif fitType == 'pe':
+            con      = [0] * self.overlay.numEVs()
+            con[idx] = 1
+            return con
+        elif fitType == 'cope':
+            return self.overlay.contrasts()[idx]
+
+        
+    def __createModelTs(self, tsType, *args, **kwargs):
+        """Creates a ``TimeSeries`` instance of the given ``tsType``, and
+        sets its display settings  according to those of this
+        ``FEATTimeSeries``.
+
+        :arg tsType: The type to create, e.g. :class:`FEATModelFitTimeSeries`,
+                     :class:`FEATEVTimeSeries`, etc.
+
+        :arg args:   Passed to the ``tsType`` constructor.
+
+        :arg kwargs: Passed to the ``tsType`` constructor.
+        """
+
+        ts = tsType(self.tsPanel, self.overlay, self.coords, *args, **kwargs)
+
+        ts.alpha     = self.alpha
+        ts.label     = self.label
+        ts.lineWidth = self.lineWidth
+        ts.lineStyle = self.lineStyle
+
+        if   isinstance(ts, FEATPartialFitTimeSeries):
+            ts.colour = (0, 0.6, 0.6)
+        elif isinstance(ts, FEATResidualTimeSeries):
+            ts.colour = (0.8, 0.4, 0)
+        elif isinstance(ts, FEATEVTimeSeries):
+            ts.colour = (0, 0.7, 0.35)
+        elif isinstance(ts, FEATModelFitTimeSeries):
+            if   ts.fitType == 'full': ts.colour = (0,   0, 1)
+            elif ts.fitType == 'cope': ts.colour = (0,   1, 0)
+            elif ts.fitType == 'pe':   ts.colour = (0.7, 0, 0)
+
+        return ts
+
+
+    def __plotPartialChanged(self, *a):
+        """Called when the :attr:`plotPartial` setting changes.
+
+        If necessary, creates and caches a :class:`FEATPartialFitTimeSeries`
+        instance.
+        """
+            
+        partial = self.plotPartial
+
+        if partial == 'none' and self.__partialTs is not None:
+            self.__partialTs = None
+            return
+
+        partial = partial.split()[0]
+
+        # fitType is either 'cope' or 'pe'
+        fitType = partial[:-1].lower()
+        idx     = int(partial[-1]) - 1
+
+        self.__partialTs = self.__createModelTs(
+            FEATPartialFitTimeSeries,
+            self.__getContrast(fitType, idx),
+            fitType,
+            idx) 
+
+
+    def __plotResidualsChanged(self, *a):
+        """Called when the :attr:`plotResiduals` setting changes.
+
+        If necessary, creates and caches a :class:`FEATResidualTimeSeries`
+        instance.
+        """ 
+        
+        if not self.plotResiduals:
+            self.__resTs = None
+            return
+
+        self.__resTs = self.__createModelTs(FEATResidualTimeSeries)
+
+
+    def __plotEVChanged(self, *a):
+        """Called when the :attr:`plotEVs` setting changes.
+
+        If necessary, creates and caches one or more :class:`FEATEVTimeSeries`
+        instances.
+        """ 
+
+        for evnum, plotEV in enumerate(self.plotEVs):
+
+            if not self.plotEVs[evnum]:
+                self.__evTs[evnum] = None
+                
+            elif self.__evTs[evnum] is None:
+                self.__evTs[evnum] = self.__createModelTs(
+                    FEATEVTimeSeries, evnum)
+            
+    
+    def __plotCOPEFitChanged(self, *a):
+        """Called when the :attr:`plotCOPEFits` setting changes.
+
+        If necessary, creates and caches one or more
+        :class:`FEATModelFitTimeSeries` instances.
+        """ 
+
+        for copenum, plotCOPE in enumerate(self.plotCOPEFits):
+
+            if not self.plotCOPEFits[copenum]:
+                self.__copeTs[copenum] = None
+            
+            elif self.__copeTs[copenum] is None:
+                self.__copeTs[copenum] = self.__createModelTs(
+                    FEATModelFitTimeSeries,
+                    self.__getContrast('cope', copenum),
+                    'cope',
+                    copenum)
+
+
+    def __plotPEFitChanged(self, evnum):
+        """Called when the :attr:`plotPEFits` setting changes.
+
+        If necessary, creates and caches one or more
+        :class:`FEATModelFitTimeSeries` instances.
+        """         
+
+        for evnum, plotPE in enumerate(self.plotPEFits):
+
+            if not self.plotPEFits[evnum]:
+                self.__peTs[evnum] = None
+
+            elif self.__peTs[evnum] is None:
+                self.__peTs[evnum] = self.__createModelTs(
+                    FEATModelFitTimeSeries,
+                    self.__getContrast('pe', evnum),
+                    'pe',
+                    evnum)
+
+
+    def __plotFullModelFitChanged(self, *a):
+        """Called when the :attr:`plotFullModelFit` setting changes.
+
+        If necessary, creates and caches a
+        :class:`FEATModelFitTimeSeries` instance.
+        """
+        
+        if not self.plotFullModelFit:
+            self.__fullModelTs = None
+            return
+
+        self.__fullModelTs = self.__createModelTs(
+            FEATModelFitTimeSeries, self.__getContrast('full', -1), 'full', -1)
+
+
+class FEATPartialFitTimeSeries(TimeSeries):
+    """A :class:`TimeSeries` class which represents the partial model fit
+    of an EV or contrast from a FEAT analysis. Instances of this class
+    are created by the :class:`FEATTimeSeries` class.
+    """
+    def __init__(self, tsPanel, overlay, coords, contrast, fitType, idx):
+        """Create a ``FEATPartialFitTimeSeries``.
+
+        :arg tsPanel:  The :class:`TimeSeriesPanel` that owns this
+                       ``FEATPartialFitTimeSeries`` instance.
+        
+        :arg overlay:  A :class:`.FEATImage` overlay.
+        
+        :arg coords:   Voxel coordinates.
+        
+        :arg contrast: The contrast vector to calculate the partial model
+                       fit for.
+        
+        :arg fitType:  The model fit type, either ``'full'``, ``'pe'`` or
+                       ``'cope'``.
+        
+        :arg idx:      If the model fit type is ``'pe'`` or ``'cope'``,
+                       the EV/contrast index.
+        """
+        TimeSeries.__init__(self, tsPanel, overlay, coords)
+
+        self.contrast = contrast
+        self.fitType  = fitType
+        self.idx      = idx
+
+        
+    def getData(self):
+        """Returns the partial model fit for the voxel and model fit type
+        specified in the constructop.
+
+        See the :meth:`.FEATImage.partialFit` method.
+        """
+        
+        data = self.overlay.partialFit(self.contrast, self.coords, False)
+        return TimeSeries.getData(self, ydata=data)
+
+    
+class FEATEVTimeSeries(TimeSeries):
+    """A :class:`TimeSeries` class which represents the time course of an
+    EV from a FEAT analysis. Instances of this class are created by the
+    :class:`FEATTimeSeries` class.
+    """
+    
+    def __init__(self, tsPanel, overlay, coords, idx):
+        """Create a ``FEATEVTimeSeries``.
+
+        :arg tsPanel:  The :class:`TimeSeriesPanel` that owns this
+                       ``FEATEVTimeSeries`` instance.
+        
+        :arg overlay:  A :class:`.FEATImage` overlay.
+        
+        :arg coords:   Voxel coordinates.
+
+        :arg idx:      The EV index.
+        """
+        TimeSeries.__init__(self, tsPanel, overlay, coords)
+        self.idx = idx
+
+        
+    def getData(self):
+        """Returns the time course of the EV specified in the constructor. """
+        data = self.overlay.getDesign()[:, self.idx]
+        return TimeSeries.getData(self, ydata=data)
+    
+
+class FEATResidualTimeSeries(TimeSeries):
+    """A :class:`TimeSeries` class which represents the time course of the
+    residuals from a FEAT analysis. Instances of this class are created by
+    the :class:`FEATTimeSeries` class.
+    """
+
+    
+    def getData(self):
+        """Returns the residuals for the voxel specified in the constructor.
+        """
+        x, y, z = self.coords
+        data    = self.overlay.getResiduals().data[x, y, z, :]
+        
+        return TimeSeries.getData(self, ydata=np.array(data))
+            
+
+class FEATModelFitTimeSeries(TimeSeries):
+    """A :class:`TimeSeries` class which represents the time course for 
+    a model fit from a FEAT analysis. Instances of this class are created by
+    the :class:`FEATTimeSeries` class.
+    """ 
+
+    def __init__(self, tsPanel, overlay, coords, contrast, fitType, idx):
+        """Create a ``FEATModelFitTimeSeries``.
+        
+        :arg tsPanel:  The :class:`TimeSeriesPanel` that owns this
+                       ``FEATModelFitTimeSeries`` instance.
+        
+        :arg overlay:  A :class:`.FEATImage` overlay.
+        
+        :arg coords:   Voxel coordinates.
+        
+        :arg contrast: The contrast vector to calculate the partial model
+                       fit for.
+        
+        :arg fitType:  The model fit type, either ``'full'``, ``'pe'`` or
+                       ``'cope'``.
+        
+        :arg idx:      If the model fit type is ``'pe'`` or ``'cope'``,
+                       the EV/contrast index.
+        """
+        
+        if fitType not in ('full', 'cope', 'pe'):
+            raise ValueError('Unknown model fit type {}'.format(fitType))
+        
+        TimeSeries.__init__(self, tsPanel, overlay, coords)
+        self.fitType  = fitType
+        self.idx      = idx
+        self.contrast = contrast
+        self.__updateModelFit()
+
+        
+    def update(self, coords):
+        """Overrides :meth:`TimeSeries.update`.
+
+        Updates the coordinates and the data encapsulated by this
+        ``FEATModelFitTimeSeries``.
+        """
+        if not TimeSeries.update(self, coords):
+            return
+        self.__updateModelFit()
+        
+
+    def __updateModelFit(self):
+        """Called by :meth:`update`, and in the constructor.  Updates the model
+        fit. See the :meth:`.FEATImage.fit` method.
+        """
+
+        fitType   = self.fitType
+        contrast  = self.contrast
+        xyz       = self.coords
+        self.data = self.overlay.fit(contrast, xyz, fitType == 'full')
+
+
+
+class MelodicTimeSeries(TimeSeries):
+
+    def __init__(self, tsPanel, overlay, component):
+        TimeSeries.__init__(self, tsPanel, overlay, component)
+
+
+    def _getData(self, component):
+        return self.overlay.getComponentTimeSeries(component)
diff --git a/fsl/fsleyes/views/histogrampanel.py b/fsl/fsleyes/views/histogrampanel.py
index 43c08686d65ff956d4ba4ab72b1f0a9a4a0c06c2..b035c9fdbd87b7253383c3e5ebabaa59e13e8d2b 100644
--- a/fsl/fsleyes/views/histogrampanel.py
+++ b/fsl/fsleyes/views/histogrampanel.py
@@ -22,13 +22,13 @@ import logging
 
 import wx
 
-import numpy as np
 
 import props
 
 import fsl.data.image                             as fslimage
 import fsl.data.strings                           as strings
 import fsl.utils.dialog                           as fsldlg
+import fsl.fsleyes.plotting.histogramseries       as histogramseries
 import fsl.fsleyes.controls.histogramcontrolpanel as histogramcontrolpanel
 import fsl.fsleyes.controls.histogramlistpanel    as histogramlistpanel
 import                                               plotpanel
@@ -198,11 +198,11 @@ class HistogramPanel(plotpanel.PlotPanel):
         if self.__current is None:
             return None
 
-        return HistogramSeries(self.__current.overlay,
-                               self,
-                               self._displayCtx,
-                               self._overlayList,
-                               baseHs=self.__current) 
+        return histogramseries.HistogramSeries(self.__current.overlay,
+                                               self,
+                                               self._displayCtx,
+                                               self._overlayList,
+                                               baseHs=self.__current) 
 
 
     def draw(self, *a):
@@ -331,11 +331,11 @@ class HistogramPanel(plotpanel.PlotPanel):
             baseHs = None
 
         def loadHs():
-            return HistogramSeries(overlay,
-                                   self,
-                                   self._displayCtx,
-                                   self._overlayList,
-                                   baseHs=baseHs)
+            return histogramseries.HistogramSeries(overlay,
+                                                   self,
+                                                   self._displayCtx,
+                                                   self._overlayList,
+                                                   baseHs=baseHs)
 
         # We are creating a new HS instance, so it
         # needs to do some initla data range/histogram
@@ -368,463 +368,3 @@ class HistogramPanel(plotpanel.PlotPanel):
         hs.label       = None
 
         self.__current = hs
-
-
-class HistogramSeries(plotpanel.DataSeries):
-    """A ``HistogramSeries`` generates histogram data from an :class:`.Image`
-    instance.
-    """
-
-    nbins = props.Int(minval=10, maxval=500, default=100, clamped=True)
-    """Number of bins to use in the histogram. This value is overridden
-    by the :attr:`HistogramPanel.autoBin` setting.
-
-    .. note:: I'm not sure why ``autoBin`` is a :class:`HistogramPanel`
-              setting, rather than a ``HistogramSeries`` setting. I might 
-              change this some time.
-    """
-
-    
-    ignoreZeros = props.Boolean(default=True)
-    """If ``True``, zeros are excluded from the calculated histogram. """
-
-    
-    showOverlay = props.Boolean(default=False)
-    """If ``True``, a 3D mask :class:`.Image` overlay is added to the
-    :class:`.OverlayList`, which highlights the voxels that have been included
-    in the histogram.
-    """
-
-
-    includeOutliers = props.Boolean(default=False)
-    """If ``True``, values which are outside of the :attr:`dataRange` are
-    included in the histogram end bins.
-    """
-
-    
-    volume = props.Int(minval=0, maxval=0, clamped=True)
-    """If the :class:`.Image` overlay associated with this ``HistogramSeries`` 
-    is 4D, this settings specifies the index of the volume that the histogram
-    is calculated upon.
-
-    .. note:: Calculating the histogram over an entire 4D :class:`.Image` is
-              not yet supported.
-    """
-
-    
-    dataRange = props.Bounds(ndims=1)
-    """Specifies the range of data which should be included in the histogram.
-    See the :attr:`includeOutliers` property.
-    """
-
-
-    def __init__(self,
-                 overlay,
-                 hsPanel,
-                 displayCtx,
-                 overlayList,
-                 volume=0,
-                 baseHs=None):
-        """Create a ``HistogramSeries``.
-
-        :arg overlay:     The :class:`.Image` overlay to calculate a histogram
-                          for.
-        
-        :arg hsPanel:     The :class:`HistogramPanel` that is displaying this
-                          ``HistogramSeries``.
-        
-        :arg displayCtx:  The :class:`.DisplayContext` instance.
-        
-        :arg overlayList: The :class:`.OverlayList` instance.
-        
-        :arg volume:      If the ``overlay`` is 4D, the initial value for the
-                          :attr:`volume` property.
-        
-        :arg baseHs:      If a ``HistogramSeries`` has already been created
-                          for the ``overlay``, it may be passed in here, so
-                          that the histogram data can be copied instead of
-                          having to be re-calculated.
-        """
-
-        log.debug('New HistogramSeries instance for {} '
-                  '(based on existing instance: {})'.format(
-                      overlay.name, baseHs is not None)) 
-
-        plotpanel.DataSeries.__init__(self, overlay)
-
-        self.volume        = volume
-        self.__hsPanel     = hsPanel
-        self.__name        = '{}_{}'.format(type(self).__name__, id(self))
-
-        self.__displayCtx  = displayCtx
-        self.__overlayList = overlayList
-        self.__overlay3D   = None
-
-        if overlay.is4DImage():
-            self.setConstraint('volume', 'maxval', overlay.shape[3] - 1)
-
-        # If we have a baseHS, we 
-        # can copy all its data
-        if baseHs is not None:
-            self.dataRange.xmin     = baseHs.dataRange.xmin
-            self.dataRange.xmax     = baseHs.dataRange.xmax
-            self.dataRange.x        = baseHs.dataRange.x
-            self.nbins              = baseHs.nbins
-            self.volume             = baseHs.volume
-            self.ignoreZeros        = baseHs.ignoreZeros
-            self.includeOutliers    = baseHs.includeOutliers
-            
-            self.__nvals              =          baseHs.__nvals
-            self.__xdata              = np.array(baseHs.__xdata)
-            self.__ydata              = np.array(baseHs.__ydata)
-            self.__finiteData         = np.array(baseHs.__finiteData)
-            self.__nonZeroData        = np.array(baseHs.__nonZeroData)
-            self.__clippedFiniteData  = np.array(baseHs.__finiteData)
-            self.__clippedNonZeroData = np.array(baseHs.__nonZeroData) 
-
-        # Otherwise we need to calculate
-        # it all for ourselves
-        else:
-            self.__initProperties()
-        
-        overlayList.addListener('overlays',
-                                self.__name,
-                                self.__overlayListChanged)
-        self       .addListener('volume',
-                                self.__name,
-                                self.__volumeChanged)
-        self       .addListener('dataRange',
-                                self.__name,
-                                self.__dataRangeChanged)
-        self       .addListener('nbins',
-                                self.__name,
-                                self.__histPropsChanged)
-        self       .addListener('ignoreZeros',
-                                self.__name,
-                                self.__histPropsChanged)
-        self       .addListener('includeOutliers',
-                                self.__name,
-                                self.__histPropsChanged)
-        self       .addListener('showOverlay',
-                                self.__name,
-                                self.__showOverlayChanged)
-
-        
-    def update(self):
-        """This method may be called to force re-calculation of the
-        histogram data.
-        """
-        self.__histPropsChanged()
-
-        
-    def destroy(self):
-        """This needs to be called when this ``HistogramSeries`` instance
-        is no longer being used.
-
-        It removes several property listeners and, if the :attr:`overlay3D`
-        property is ``True``, removes the mask overlay from the
-        :class:`.OverlayList`.
-        """
-        self              .removeListener('nbins',           self.__name)
-        self              .removeListener('ignoreZeros',     self.__name)
-        self              .removeListener('includeOutliers', self.__name)
-        self              .removeListener('volume',          self.__name)
-        self              .removeListener('dataRange',       self.__name)
-        self              .removeListener('nbins',           self.__name)
-        self.__overlayList.removeListener('overlays',        self.__name)
-        
-        if self.__overlay3D is not None:
-            self.__overlayList.remove(self.__overlay3D)
-            self.__overlay3D = None
-
-        
-    def __initProperties(self):
-        """Called by :meth:`__init__`. Calculates and caches some things which
-        are needed for the histogram calculation.
-
-        .. note:: This method is never called if a ``baseHs`` is provided to
-                 :meth:`__init__`.
-        """
-
-        log.debug('Performining initial histogram '
-                  'calculations for overlay {}'.format(
-                      self.overlay.name))
-
-        data  = self.overlay.data[:]
-        
-        finData = data[np.isfinite(data)]
-        dmin    = finData.min()
-        dmax    = finData.max()
-        dist    = (dmax - dmin) / 10000.0
-        
-        nzData = finData[finData != 0]
-        nzmin  = nzData.min()
-        nzmax  = nzData.max()
-
-        self.dataRange.xmin = dmin
-        self.dataRange.xmax = dmax  + dist
-        self.dataRange.xlo  = nzmin
-        self.dataRange.xhi  = nzmax + dist
-
-        self.nbins = self.__autoBin(nzData, self.dataRange.x)
-
-        if not self.overlay.is4DImage():
-            
-            self.__finiteData  = finData
-            self.__nonZeroData = nzData
-            self.__dataRangeChanged(callHistPropsChanged=False)
-            
-        else:
-            self.__volumeChanged(callHistPropsChanged=False)
-
-        self.__histPropsChanged()
-
-        
-    def __volumeChanged(
-            self,
-            ctx=None,
-            value=None,
-            valid=None,
-            name=None,
-            callHistPropsChanged=True):
-        """Called when the :attr:`volume` property changes, and also by the
-        :meth:`__initProperties` method.
-
-        Re-calculates some things for the new overlay volume.
-
-        :arg callHistPropsChanged: If ``True`` (the default), the
-                                   :meth:`__histPropsChanged` method will be
-                                   called.
-
-        All other arguments are ignored, but are passed in when this method is
-        called due to a property change (see the
-        :meth:`.HasProperties.addListener` method).        
-        """
-
-        if self.overlay.is4DImage(): data = self.overlay.data[..., self.volume]
-        else:                        data = self.overlay.data[:]
-
-        data = data[np.isfinite(data)]
-
-        self.__finiteData  = data
-        self.__nonZeroData = data[data != 0]
-
-        self.__dataRangeChanged(callHistPropsChanged=False)
-
-        if callHistPropsChanged:
-            self.__histPropsChanged()
-
-
-    def __dataRangeChanged(
-            self,
-            ctx=None,
-            value=None,
-            valid=None,
-            name=None,
-            callHistPropsChanged=True):
-        """Called when the :attr:`dataRange` property changes, and also by the
-        :meth:`__initProperties` and :meth:`__volumeChanged` methods.
-
-        :arg callHistPropsChanged: If ``True`` (the default), the
-                                   :meth:`__histPropsChanged` method will be
-                                   called.
-
-        All other arguments are ignored, but are passed in when this method is
-        called due to a property change (see the
-        :meth:`.HasProperties.addListener` method).
-        """
-        finData = self.__finiteData
-        nzData  = self.__nonZeroData
-        
-        self.__clippedFiniteData  = finData[(finData >= self.dataRange.xlo) &
-                                            (finData <  self.dataRange.xhi)]
-        self.__clippedNonZeroData = nzData[ (nzData  >= self.dataRange.xlo) &
-                                            (nzData  <  self.dataRange.xhi)]
-
-        if callHistPropsChanged:
-            self.__histPropsChanged()
-
-            
-    def __histPropsChanged(self, *a):
-        """Called internally, and when any histogram settings change.
-        Re-calculates the histogram data.
-        """
-
-        log.debug('Calculating histogram for '
-                  'overlay {}'.format(self.overlay.name))
-
-        if self.dataRange.xhi - self.dataRange.xlo < 0.00000001:
-            self.__xdata = np.array([])
-            self.__ydata = np.array([])
-            self.__nvals = 0
-            return
-
-        if self.ignoreZeros:
-            if self.includeOutliers: data = self.__nonZeroData
-            else:                    data = self.__clippedNonZeroData
-        else:
-            if self.includeOutliers: data = self.__finiteData
-            else:                    data = self.__clippedFiniteData 
-        
-        if self.__hsPanel.autoBin:
-            nbins = self.__autoBin(data, self.dataRange.x)
-
-            self.disableListener('nbins', self.__name)
-            self.nbins = nbins
-            self.enableListener('nbins', self.__name)
-
-        # Calculate bin edges
-        bins = np.linspace(self.dataRange.xlo,
-                           self.dataRange.xhi,
-                           self.nbins + 1)
-
-        if self.includeOutliers:
-            bins[ 0] = self.dataRange.xmin
-            bins[-1] = self.dataRange.xmax
-            
-        # Calculate the histogram
-        histX    = bins
-        histY, _ = np.histogram(data.flat, bins=bins)
-            
-        self.__xdata = histX
-        self.__ydata = histY
-        self.__nvals = histY.sum()
-
-        log.debug('Calculated histogram for overlay '
-                  '{} (number of values: {}, number '
-                  'of bins: {})'.format(
-                      self.overlay.name,
-                      self.__nvals,
-                      self.nbins))
-
-
-    def __showOverlayChanged(self, *a):
-        """Called when the :attr:`showOverlay` property changes.
-
-        Adds/removes a 3D mask :class:`.Image` to the :class:`.OverlayList`,
-        which highlights the voxels that have been included in the histogram.
-        The :class:`.MaskOpts.threshold` property is bound to the
-        :attr:`dataRange` property, so the masked voxels are updated whenever
-        the histogram data range changes, and vice versa.
-        """
-
-        if not self.showOverlay:
-            if self.__overlay3D is not None:
-
-                log.debug('Removing 3D histogram overlay mask for {}'.format(
-                    self.overlay.name))
-                self.__overlayList.remove(self.__overlay3D)
-                self.__overlay3D = None
-
-        else:
-
-            log.debug('Creating 3D histogram overlay mask for {}'.format(
-                self.overlay.name))
-            
-            self.__overlay3D = fslimage.Image(
-                self.overlay.data,
-                name='{}/histogram/mask'.format(self.overlay.name),
-                header=self.overlay.nibImage.get_header())
-
-            self.__overlayList.append(self.__overlay3D)
-
-            opts = self.__displayCtx.getOpts(self.__overlay3D,
-                                             overlayType='mask')
-
-            opts.bindProps('volume',    self)
-            opts.bindProps('colour',    self)
-            opts.bindProps('threshold', self, 'dataRange')
-
-
-    def __overlayListChanged(self, *a):
-        """Called when the :class:`.OverlayList` changes.
-
-        If a 3D mask overlay was being shown, and it has been removed from the
-        ``OverlayList``, the :attr:`showOverlay` property is updated
-        accordingly.
-        """
-        
-        if self.__overlay3D is None:
-            return
-
-        # If a 3D overlay was being shown, and it
-        # has been removed from the overlay list
-        # by the user, turn the showOverlay property
-        # off
-        if self.__overlay3D not in self.__overlayList:
-            
-            self.disableListener('showOverlay', self.__name)
-            self.showOverlay = False
-            self.__showOverlayChanged()
-            self.enableListener('showOverlay', self.__name)
-
-        
-    def __autoBin(self, data, dataRange):
-        """Calculates the number of bins which should be used for a histogram
-        of the given data. The calculation is identical to that implemented
-        in the original FSLView.
-
-        :arg data:      The data that the histogram is to be calculated on.
-
-        :arg dataRange: A tuple containing the ``(min, max)`` histogram range.
-        """
-
-        dMin, dMax = dataRange
-        dRange     = dMax - dMin
-
-        binSize = np.power(10, np.ceil(np.log10(dRange) - 1) - 1)
-
-        nbins = dRange / binSize
-
-        while nbins < 100:
-            binSize /= 2
-            nbins    = dRange / binSize
-
-        if issubclass(data.dtype.type, np.integer):
-            binSize = max(1, np.ceil(binSize))
-
-        adjMin = np.floor(dMin / binSize) * binSize
-        adjMax = np.ceil( dMax / binSize) * binSize
-
-        nbins = int((adjMax - adjMin) / binSize) + 1
-
-        return nbins
-            
-
-    def getData(self):
-        """Overrides :meth:`.DataSeries.getData`.
-
-        Returns  a tuple containing the ``(x, y)`` histogram data.
-        """
-
-        if len(self.__xdata) == 0 or \
-           len(self.__ydata) == 0:
-            return self.__xdata, self.__ydata
-
-        # If smoothing is not enabled, we'll
-        # munge the histogram data a bit so
-        # that plt.plot(drawstyle='steps-pre')
-        # plots it nicely.
-        if not self.__hsPanel.smooth:
-
-            xdata = np.zeros(len(self.__xdata) + 1, dtype=np.float32)
-            ydata = np.zeros(len(self.__ydata) + 2, dtype=np.float32)
-
-            xdata[ :-1] = self.__xdata
-            xdata[  -1] = self.__xdata[-1]
-            ydata[1:-1] = self.__ydata
-
-
-        # If smoothing is enabled, the above munge
-        # is not necessary, and will probably cause
-        # the spline interpolation (performed by 
-        # the PlotPanel) to fail.
-        else:
-            xdata = np.array(self.__xdata[:-1], dtype=np.float32)
-            ydata = np.array(self.__ydata,      dtype=np.float32)
-
-        nvals    = self.__nvals
-        histType = self.__hsPanel.histType
-            
-        if   histType == 'count':       return xdata, ydata
-        elif histType == 'probability': return xdata, ydata / nvals
diff --git a/fsl/fsleyes/views/plotpanel.py b/fsl/fsleyes/views/plotpanel.py
index c69b7c484400616d292717fd0d4adf1a9a68db14..0a64d0e59dd5dcafbdf053cb7f0783f3806eb696 100644
--- a/fsl/fsleyes/views/plotpanel.py
+++ b/fsl/fsleyes/views/plotpanel.py
@@ -4,10 +4,9 @@
 #
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
-"""This module provides the :class:`PlotPanel` and :class:`DataSeries`
-classes. The ``PlotPanel`` class is the base class for all *FSLeyes views*
-which display some sort of data plot.  See the :mod:`~fsl.fsleyes` package
-documentation for more details..
+"""This module provides the :class:`PlotPanel` class. The ``PlotPanel`` class
+is the base class for all *FSLeyes views* which display some sort of data
+plot.  See the :mod:`~fsl.fsleyes` package documentation for more details..
 """
 
 
@@ -54,7 +53,7 @@ class PlotPanel(viewpanel.ViewPanel):
 
       1. Call the ``PlotPanel`` constructor.
 
-      2. Define a :class:`DataSeries` sub-class.
+      2. Define a :class:`.DataSeries` sub-class.
 
       3. Override the :meth:`draw` method, so it calls the
          :meth:`drawDataSeries` method.
@@ -66,8 +65,9 @@ class PlotPanel(viewpanel.ViewPanel):
     **Data series**
 
     A ``PlotPanel`` instance plots data contained in one or more
-    :class:`DataSeries` instances.  Therefore, ``PlotPanel`` sub-classes also
-    need to define a sub-class of the :class:`DataSeries` base class.
+    :class:`.DataSeries` instances; all ``DataSeries`` classes are defined in
+    the :mod:`.plotting` sub-package.  Therefore, ``PlotPanel`` sub-classes
+    also need to define a sub-class of the :class:`.DataSeries` base class.
 
     ``DataSeries`` objects can be plotted by passing them to the
     :meth:`drawDataSeries` method.
@@ -105,7 +105,7 @@ class PlotPanel(viewpanel.ViewPanel):
 
 
     dataSeries = props.List()
-    """This list contains :class:`DataSeries` instances which are plotted
+    """This list contains :class:`.DataSeries` instances which are plotted
     on every call to :meth:`drawDataSeries`. ``DataSeries`` instances can
     be added/removed directly to/from this list.
     """
@@ -266,7 +266,7 @@ class PlotPanel(viewpanel.ViewPanel):
     def draw(self, *a):
         """This method must be overridden by ``PlotPanel`` sub-classes.
 
-        It is called whenever a :class:`DataSeries` is added to the
+        It is called whenever a :class:`.DataSeries` is added to the
         :attr:`dataSeries` list, or when any plot display properties change.
 
         Sub-class implementations should call the :meth:`drawDataSeries`
@@ -356,7 +356,7 @@ class PlotPanel(viewpanel.ViewPanel):
 
 
     def drawDataSeries(self, extraSeries=None, **plotArgs):
-        """Plots all of the :class:`DataSeries` instances in the
+        """Plots all of the :class:`.DataSeries` instances in the
         :attr:`dataSeries` list
 
         :arg extraSeries: A sequence of additional ``DataSeries`` to be
@@ -475,7 +475,7 @@ class PlotPanel(viewpanel.ViewPanel):
 
         
     def __drawOneDataSeries(self, ds, **plotArgs):
-        """Plots a single :class:`DataSeries` instance. This method is called
+        """Plots a single :class:`.DataSeries` instance. This method is called
         by the :meth:`drawDataSeries` method.
 
         :arg ds:       The ``DataSeries`` instance.
@@ -586,7 +586,7 @@ class PlotPanel(viewpanel.ViewPanel):
 
     def __dataSeriesChanged(self, *a):
         """Called when the :attr:`dataSeries` list changes. Adds listeners
-        to any new :class:`DataSeries` instances, and then calls :meth:`draw`.
+        to any new :class:`.DataSeries` instances, and then calls :meth:`draw`.
         """
         
         for ds in self.dataSeries:
@@ -666,69 +666,3 @@ class PlotPanel(viewpanel.ViewPanel):
         self.enableListener('limits', self.__name)            
  
         return (xmin, xmax), (ymin, ymax)
-
-
-class DataSeries(props.HasProperties):
-    """A ``DataSeries`` instance encapsulates some data to be plotted by
-    a :class:`PlotPanel`, with the data extracted from an overlay in the
-    :class:`.OverlayList`. 
-
-    Sub-class implementations must accept an overlay object, pass this
-    overlay to the ``DataSeries`` constructor, and override the
-    :meth:`getData` method. The overlay is accessible as an instance
-    attribute, confusingly called ``overlay``.
-
-    Each``DataSeries`` instance is plotted as a line, with the line
-    style defined by properties on the ``DataSeries`` instance,
-    such as :attr:`colour`, :attr:`lineWidth` etc.
-    """
-
-    colour = props.Colour()
-    """Line colour. """
-
-    
-    alpha = props.Real(minval=0.0, maxval=1.0, default=1.0, clamped=True)
-    """Line transparency."""
-
-    
-    label = props.String()
-    """Line label (used in the plot legend)."""
-
-    
-    lineWidth = props.Choice((0.5, 1, 2, 3, 4, 5))
-    """Line width. """
-
-    
-    lineStyle = props.Choice(('-', '--', '-.', ':'))
-    """Line style. """
-
-    
-    def __init__(self, overlay):
-        """Create a ``DataSeries``.
-
-        :arg overlay: The overlay from which the data to be plotted is
-                      derived. 
-        """
-        
-        self.overlay = overlay
-
-
-    def __copy__(self):
-        """``DataSeries`` copy operator. Sub-classes with constructors
-        that require more than just the overlay object will need to
-        implement their own copy operator.
-        """
-        return type(self)(self.overlay)
-
-
-    def getData(self):
-        """This method must be implemented by sub-classes. It must return
-        the data to be plotted, as a tuple of the form:
-        
-            ``(xdata, ydata)``
-
-        where ``xdata`` and ``ydata`` are sequences containing the x/y data
-        to be plotted.
-        """
-        raise NotImplementedError('The getData method must be '
-                                  'implemented by subclasses')
diff --git a/fsl/fsleyes/views/timeseriespanel.py b/fsl/fsleyes/views/timeseriespanel.py
index 79bfcf41d2ff15f9b7fd1a8ed6a28b8d6223adbf..a087d025598c3ccea0955f88866d581b85d8f6de 100644
--- a/fsl/fsleyes/views/timeseriespanel.py
+++ b/fsl/fsleyes/views/timeseriespanel.py
@@ -6,17 +6,11 @@
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
 """This module provides the :class:`TimeSeriesPanel`, which is a *FSLeyes
-view* for displaying time series data from :class:`.Image` overlays. A
-``TimeSeriesPanel`` looks something like the following:
+view* for displaying time series data from :class:`.Image` overlays.
 
-.. image:: images/timeseriespanel.png
-   :scale: 50%
-   :align: center
-
-
-The following :class:`.DataSeries` sub-classes are used by the
-:class:`TimeSeriesPanel` (see the :class:`.PlotPanel` documentation for more
-details):
+The following :class:`.DataSeries` sub-classes, defined in the
+:mod:`.plotting.timeseries` module, are used by the :class:`TimeSeriesPanel`
+(see the :class:`.PlotPanel` documentation for more details):
 
 .. autosummary::
    :nosignatures:
@@ -42,6 +36,7 @@ import fsl.data.melodicimage                       as fslmelimage
 import fsl.data.image                              as fslimage
 import fsl.fsleyes.displaycontext                  as fsldisplay
 import fsl.fsleyes.colourmaps                      as fslcmaps
+import fsl.fsleyes.plotting.timeseries             as timeseries
 import fsl.fsleyes.controls.timeseriescontrolpanel as timeseriescontrolpanel
 import fsl.fsleyes.controls.timeserieslistpanel    as timeserieslistpanel
 
@@ -49,622 +44,15 @@ import fsl.fsleyes.controls.timeserieslistpanel    as timeserieslistpanel
 log = logging.getLogger(__name__)
 
 
-class TimeSeries(plotpanel.DataSeries):
-    """Encapsulates time series data from a specific voxel in an
-    :class:`.Image` overlay.
-    """
-
-    
-    def __init__(self, tsPanel, overlay, coords):
-        """Create a ``TimeSeries`` instane.
-
-        :arg tsPanel: The :class:`TimeSeriesPanel` which owns this
-                      ``TimeSeries``.
-
-        :arg overlay: The :class:`.Image` instance to extract the data from.
-
-        :arg coords:  The voxel coordinates of the time series.
-        """
-        plotpanel.DataSeries.__init__(self, overlay)
-
-        self.tsPanel = tsPanel
-        self.coords  = None
-        self.data    = None
-        self.update(coords)
-
-
-    def __copy__(self):
-        """Copy operator, creates and returns a copy of this ``TimeSeries``
-        instance.
-        """
-        return type(self)(self.tsPanel, self.overlay, self.coords)
-
-        
-    def update(self, coords):
-        """Updates the voxel coordinates and time series data encapsulated
-        by this ``TimeSeries``.
-
-        .. warning:: This method is only intended for use on the *current*
-                     ``TimeSeriesPanel`` series, not for time series instances
-                     which have been added to the :attr:`.PlotPanel.dataSeries`
-                     list.
-        
-                     It is used by the :class:`TimeSeriesPanel` so that
-                     ``TimeSeries`` instances do not have to constantly be
-                     destroyed/recreated whenever the
-                     :attr:`.DisplayContext.location` changes.
-        """
-
-        self.coords = coords
-        self.data = self._getData(coords)
-        return True
-
-
-    def _getData(self, coords):
-        return self.overlay.data[coords[0], coords[1], coords[2], :]
-    
-        
-    def getData(self, xdata=None, ydata=None):
-        """Overrides :meth:`.DataSeries.getData` Returns the data associated
-        with this ``TimeSeries`` instance.
-
-        The ``xdata`` and ``ydata`` arguments may be used by subclasses to
-        override the x/y data in the event that they have already performed
-        some processing on the data.
-        """
-
-        if xdata is None: xdata = np.arange(len(self.data), dtype=np.float32)
-        if ydata is None: ydata = np.array(     self.data,  dtype=np.float32)
-
-        if self.tsPanel.usePixdim:
-            xdata *= self.overlay.pixdim[3]
-        
-        if self.tsPanel.plotMode == 'demean':
-            ydata = ydata - ydata.mean()
-
-        elif self.tsPanel.plotMode == 'normalise':
-            ymin  = ydata.min()
-            ymax  = ydata.max()
-            ydata = 2 * (ydata - ymin) / (ymax - ymin) - 1
-            
-        elif self.tsPanel.plotMode == 'percentChange':
-            mean  = ydata.mean()
-            ydata =  100 * (ydata / mean) - 100
-            
-        return xdata, ydata
-
- 
-class FEATTimeSeries(TimeSeries):
-    """A :class:`TimeSeries` class for use with :class:`FEATImage` instances,
-    containing some extra FEAT specific options.
-
-    The ``FEATTimeSeries`` class acts as a container for several
-    ``TimeSeries`` instances, each of which represent some part of a FEAT
-    analysis. Therefore, the data returned by a call to
-    :meth:`.TimeSeries.getData` on a ``FEATTimeSeries`` instance should not
-    be plotted.
-
-    Instead, the :meth:`getModelTimeSeries` method should be used to retrieve
-    a list of all the ``TimeSeries`` instances which are associated with the
-    ``FEATTimeSeries`` instance - all of these ``TimeSeries`` instances should
-    be plotted instead.
-
-    For example, if the :attr:`plotData` and :attr:`plotFullModelFit` settings
-    are ``True``, the :meth:`getModelTimeSeries` method will return a list
-    containing two ``TimeSeries`` instances - one which will return the FEAT
-    analysis input data, and another which will return the full model fit, for
-    the voxel in question.
-
-
-    .. note:: The ``getData`` method of a ``FEATTimeSeries`` instance will
-              return the FEAT input data; therefore, when :attr:`plotData` is
-              ``True``, the ``FEATTimeSeries`` instance will itself be included
-              in the list returned by :meth:`getModelTimeSeries`.
-
-    
-    The following classes are used to represent the various parts of a FEAT
-    analysis:
-
-    .. autosummary::
-       :nosignatures:
-    
-       FEATEVTimeSeries
-       FEATResidualTimeSeries
-       FEATPartialFitTimeSeries
-       FEATModelFitTimeSeries
-    """
-
-
-    plotData = props.Boolean(default=True)
-    """If ``True``, the FEAT input data is plotted. """
-
-    
-    plotFullModelFit = props.Boolean(default=True)
-    """If ``True``, the FEAT full model fit is plotted. """
-
-    
-    plotResiduals = props.Boolean(default=False)
-    """If ``True``, the FEAT model residuals are plotted. """
-
-    
-    plotEVs = props.List(props.Boolean(default=False))
-    """A list of ``Boolean`` properties, one for each EV in the FEAT analysis.
-    For elements that are ``True``, the corresponding FEAT EV time course is
-    plotted.
-    """
-
-    
-    plotPEFits = props.List(props.Boolean(default=False))
-    """A list of ``Boolean`` properties, one for each EV in the FEAT analysis.
-    For elements that are ``True``, the model fit for the corresponding FEAT
-    EV is plotted.    
-    """
-
-    
-    plotCOPEFits = props.List(props.Boolean(default=False))
-    """A list of ``Boolean`` properties, one for each EV in the FEAT analysis.
-    For elements that are ``True``, the model fit for the corresponding FEAT
-    contrast is plotted.    
-    """ 
-
-    
-    plotPartial = props.Choice()
-    """Plot the raw data, after regression against a chosen EV or contrast.
-    The options are populated in the :meth:`__init__` method.
-    """
-    
-
-    def __init__(self, *args, **kwargs):
-        """Create a ``FEATTimeSeries``.
-
-        All arguments are passed through to the :class:`TimeSeries`
-        constructor.
-        """
-        
-        TimeSeries.__init__(self, *args, **kwargs)
-        self.name = '{}_{}'.format(type(self).__name__, id(self))
-
-        numEVs    = self.overlay.numEVs()
-        numCOPEs  = self.overlay.numContrasts()
-        copeNames = self.overlay.contrastNames()
-        
-        reduceOpts = ['none'] + \
-                     ['PE{}'.format(i + 1) for i in range(numEVs)]
-
-        for i in range(numCOPEs):
-            name = 'COPE{} ({})'.format(i + 1, copeNames[i])
-            reduceOpts.append(name)
-
-        self.getProp('plotPartial').setChoices(reduceOpts, instance=self)
-
-        for i in range(numEVs):
-            self.plotPEFits.append(False)
-            self.plotEVs   .append(False)
-
-        for i in range(numCOPEs):
-            self.plotCOPEFits.append(False) 
-
-        self.__fullModelTs =  None
-        self.__partialTs   =  None
-        self.__resTs       =  None
-        self.__evTs        = [None] * numEVs
-        self.__peTs        = [None] * numEVs
-        self.__copeTs      = [None] * numCOPEs
-        
-        self.addListener('plotFullModelFit',
-                         self.name,
-                         self.__plotFullModelFitChanged)
-        self.addListener('plotResiduals',
-                         self.name,
-                         self.__plotResidualsChanged)
-        self.addListener('plotPartial',
-                         self.name,
-                         self.__plotPartialChanged)
-
-        self.addListener('plotEVs',      self.name, self.__plotEVChanged)
-        self.addListener('plotPEFits',   self.name, self.__plotPEFitChanged)
-        self.addListener('plotCOPEFits', self.name, self.__plotCOPEFitChanged)
-
-        # plotFullModelFit defaults to True, so
-        # force the model fit ts creation here
-        self.__plotFullModelFitChanged()
-
-
-    def __copy__(self):
-        """Copy operator for a ``FEATTimeSeries`` instance."""
-        
-        copy = type(self)(self.tsPanel, self.overlay, self.coords)
-
-        copy.colour           = self.colour
-        copy.alpha            = self.alpha 
-        copy.label            = self.label 
-        copy.lineWidth        = self.lineWidth
-        copy.lineStyle        = self.lineStyle
-
-        # When these properties are changed 
-        # on the copy instance, it will create 
-        # its own FEATModelFitTimeSeries 
-        # instances accordingly
-        copy.plotFullModelFit = self.plotFullModelFit
-        copy.plotEVs[     :]  = self.plotEVs[     :]
-        copy.plotPEFits[  :]  = self.plotPEFits[  :]
-        copy.plotCOPEFits[:]  = self.plotCOPEFits[:]
-        copy.plotPartial      = self.plotPartial
-        copy.plotResiduals    = self.plotResiduals
-
-        return copy
- 
-
-    def getModelTimeSeries(self):
-        """Returns a list containing all of the ``TimeSeries`` instances
-        which should be plotted in place of this ``FEATTimeSeries``.
-        """
-        
-        
-        modelts = []
-
-        if self.plotData:              modelts.append(self)
-        if self.plotFullModelFit:      modelts.append(self.__fullModelTs)
-        if self.plotResiduals:         modelts.append(self.__resTs)
-        if self.plotPartial != 'none': modelts.append(self.__partialTs)
-        
-        for i in range(self.overlay.numEVs()):
-            if self.plotPEFits[i]:
-                modelts.append(self.__peTs[i])
-
-        for i in range(self.overlay.numEVs()):
-            if self.plotEVs[i]:
-                modelts.append(self.__evTs[i]) 
-
-        for i in range(self.overlay.numContrasts()):
-            if self.plotCOPEFits[i]:
-                modelts.append(self.__copeTs[i])
-
-        return modelts
-
-
-    def update(self, coords):
-        """Overrides :meth:`TimeSeries.update`.
-
-        Updates the coordinates and data associated wsith this
-        ``FEATTimeSeries`` instance.
-        """
-        
-        if not TimeSeries.update(self, coords):
-            return False
-            
-        for modelTs in self.getModelTimeSeries():
-            if modelTs is self:
-                continue
-            modelTs.update(coords)
-
-        return True
-
-
-    def __getContrast(self, fitType, idx):
-        """Returns a contrast vector for the given model fit type, and index.
-
-        :arg fitType: either ``'full'``, ``'pe'``, or ``'cope'``. If
-                      ``'full'``, the ``idx`` argument is ignored.
-
-        :arg idx:     The EV or contrast index for ``'pe'`` or ``'cope'`` model
-                      fits.
-        """
-
-        if fitType == 'full':
-            return [1] * self.overlay.numEVs()
-        elif fitType == 'pe':
-            con      = [0] * self.overlay.numEVs()
-            con[idx] = 1
-            return con
-        elif fitType == 'cope':
-            return self.overlay.contrasts()[idx]
-
-        
-    def __createModelTs(self, tsType, *args, **kwargs):
-        """Creates a ``TimeSeries`` instance of the given ``tsType``, and
-        sets its display settings  according to those of this
-        ``FEATTimeSeries``.
-
-        :arg tsType: The type to create, e.g. :class:`FEATModelFitTimeSeries`,
-                     :class:`FEATEVTimeSeries`, etc.
-
-        :arg args:   Passed to the ``tsType`` constructor.
-
-        :arg kwargs: Passed to the ``tsType`` constructor.
-        """
-
-        ts = tsType(self.tsPanel, self.overlay, self.coords, *args, **kwargs)
-
-        ts.alpha     = self.alpha
-        ts.label     = self.label
-        ts.lineWidth = self.lineWidth
-        ts.lineStyle = self.lineStyle
 
-        if   isinstance(ts, FEATPartialFitTimeSeries):
-            ts.colour = (0, 0.6, 0.6)
-        elif isinstance(ts, FEATResidualTimeSeries):
-            ts.colour = (0.8, 0.4, 0)
-        elif isinstance(ts, FEATEVTimeSeries):
-            ts.colour = (0, 0.7, 0.35)
-        elif isinstance(ts, FEATModelFitTimeSeries):
-            if   ts.fitType == 'full': ts.colour = (0,   0, 1)
-            elif ts.fitType == 'cope': ts.colour = (0,   1, 0)
-            elif ts.fitType == 'pe':   ts.colour = (0.7, 0, 0)
-
-        return ts
-
-
-    def __plotPartialChanged(self, *a):
-        """Called when the :attr:`plotPartial` setting changes.
-
-        If necessary, creates and caches a :class:`FEATPartialFitTimeSeries`
-        instance.
-        """
-            
-        partial = self.plotPartial
-
-        if partial == 'none' and self.__partialTs is not None:
-            self.__partialTs = None
-            return
-
-        partial = partial.split()[0]
-
-        # fitType is either 'cope' or 'pe'
-        fitType = partial[:-1].lower()
-        idx     = int(partial[-1]) - 1
-
-        self.__partialTs = self.__createModelTs(
-            FEATPartialFitTimeSeries,
-            self.__getContrast(fitType, idx),
-            fitType,
-            idx) 
-
-
-    def __plotResidualsChanged(self, *a):
-        """Called when the :attr:`plotResiduals` setting changes.
-
-        If necessary, creates and caches a :class:`FEATResidualTimeSeries`
-        instance.
-        """ 
-        
-        if not self.plotResiduals:
-            self.__resTs = None
-            return
-
-        self.__resTs = self.__createModelTs(FEATResidualTimeSeries)
-
-
-    def __plotEVChanged(self, *a):
-        """Called when the :attr:`plotEVs` setting changes.
-
-        If necessary, creates and caches one or more :class:`FEATEVTimeSeries`
-        instances.
-        """ 
-
-        for evnum, plotEV in enumerate(self.plotEVs):
-
-            if not self.plotEVs[evnum]:
-                self.__evTs[evnum] = None
-                
-            elif self.__evTs[evnum] is None:
-                self.__evTs[evnum] = self.__createModelTs(
-                    FEATEVTimeSeries, evnum)
-            
-    
-    def __plotCOPEFitChanged(self, *a):
-        """Called when the :attr:`plotCOPEFits` setting changes.
-
-        If necessary, creates and caches one or more
-        :class:`FEATModelFitTimeSeries` instances.
-        """ 
-
-        for copenum, plotCOPE in enumerate(self.plotCOPEFits):
-
-            if not self.plotCOPEFits[copenum]:
-                self.__copeTs[copenum] = None
-            
-            elif self.__copeTs[copenum] is None:
-                self.__copeTs[copenum] = self.__createModelTs(
-                    FEATModelFitTimeSeries,
-                    self.__getContrast('cope', copenum),
-                    'cope',
-                    copenum)
-
-
-    def __plotPEFitChanged(self, evnum):
-        """Called when the :attr:`plotPEFits` setting changes.
-
-        If necessary, creates and caches one or more
-        :class:`FEATModelFitTimeSeries` instances.
-        """         
-
-        for evnum, plotPE in enumerate(self.plotPEFits):
-
-            if not self.plotPEFits[evnum]:
-                self.__peTs[evnum] = None
-
-            elif self.__peTs[evnum] is None:
-                self.__peTs[evnum] = self.__createModelTs(
-                    FEATModelFitTimeSeries,
-                    self.__getContrast('pe', evnum),
-                    'pe',
-                    evnum)
-
-
-    def __plotFullModelFitChanged(self, *a):
-        """Called when the :attr:`plotFullModelFit` setting changes.
-
-        If necessary, creates and caches a
-        :class:`FEATModelFitTimeSeries` instance.
-        """
-        
-        if not self.plotFullModelFit:
-            self.__fullModelTs = None
-            return
-
-        self.__fullModelTs = self.__createModelTs(
-            FEATModelFitTimeSeries, self.__getContrast('full', -1), 'full', -1)
-
-
-class FEATPartialFitTimeSeries(TimeSeries):
-    """A :class:`TimeSeries` class which represents the partial model fit
-    of an EV or contrast from a FEAT analysis. Instances of this class
-    are created by the :class:`FEATTimeSeries` class.
-    """
-    def __init__(self, tsPanel, overlay, coords, contrast, fitType, idx):
-        """Create a ``FEATPartialFitTimeSeries``.
-
-        :arg tsPanel:  The :class:`TimeSeriesPanel` that owns this
-                       ``FEATPartialFitTimeSeries`` instance.
-        
-        :arg overlay:  A :class:`.FEATImage` overlay.
-        
-        :arg coords:   Voxel coordinates.
-        
-        :arg contrast: The contrast vector to calculate the partial model
-                       fit for.
-        
-        :arg fitType:  The model fit type, either ``'full'``, ``'pe'`` or
-                       ``'cope'``.
-        
-        :arg idx:      If the model fit type is ``'pe'`` or ``'cope'``,
-                       the EV/contrast index.
-        """
-        TimeSeries.__init__(self, tsPanel, overlay, coords)
-
-        self.contrast = contrast
-        self.fitType  = fitType
-        self.idx      = idx
-
-        
-    def getData(self):
-        """Returns the partial model fit for the voxel and model fit type
-        specified in the constructop.
-
-        See the :meth:`.FEATImage.partialFit` method.
-        """
-        
-        data = self.overlay.partialFit(self.contrast, self.coords, False)
-        return TimeSeries.getData(self, ydata=data)
-
-    
-class FEATEVTimeSeries(TimeSeries):
-    """A :class:`TimeSeries` class which represents the time course of an
-    EV from a FEAT analysis. Instances of this class are created by the
-    :class:`FEATTimeSeries` class.
-    """
-    
-    def __init__(self, tsPanel, overlay, coords, idx):
-        """Create a ``FEATEVTimeSeries``.
-
-        :arg tsPanel:  The :class:`TimeSeriesPanel` that owns this
-                       ``FEATEVTimeSeries`` instance.
-        
-        :arg overlay:  A :class:`.FEATImage` overlay.
-        
-        :arg coords:   Voxel coordinates.
-
-        :arg idx:      The EV index.
-        """
-        TimeSeries.__init__(self, tsPanel, overlay, coords)
-        self.idx = idx
-
-        
-    def getData(self):
-        """Returns the time course of the EV specified in the constructor. """
-        data = self.overlay.getDesign()[:, self.idx]
-        return TimeSeries.getData(self, ydata=data)
-    
-
-class FEATResidualTimeSeries(TimeSeries):
-    """A :class:`TimeSeries` class which represents the time course of the
-    residuals from a FEAT analysis. Instances of this class are created by
-    the :class:`FEATTimeSeries` class.
-    """
-
-    
-    def getData(self):
-        """Returns the residuals for the voxel specified in the constructor.
-        """
-        x, y, z = self.coords
-        data    = self.overlay.getResiduals().data[x, y, z, :]
-        
-        return TimeSeries.getData(self, ydata=np.array(data))
-            
-
-class FEATModelFitTimeSeries(TimeSeries):
-    """A :class:`TimeSeries` class which represents the time course for 
-    a model fit from a FEAT analysis. Instances of this class are created by
-    the :class:`FEATTimeSeries` class.
-    """ 
-
-    def __init__(self, tsPanel, overlay, coords, contrast, fitType, idx):
-        """Create a ``FEATModelFitTimeSeries``.
-        
-        :arg tsPanel:  The :class:`TimeSeriesPanel` that owns this
-                       ``FEATModelFitTimeSeries`` instance.
-        
-        :arg overlay:  A :class:`.FEATImage` overlay.
-        
-        :arg coords:   Voxel coordinates.
-        
-        :arg contrast: The contrast vector to calculate the partial model
-                       fit for.
-        
-        :arg fitType:  The model fit type, either ``'full'``, ``'pe'`` or
-                       ``'cope'``.
-        
-        :arg idx:      If the model fit type is ``'pe'`` or ``'cope'``,
-                       the EV/contrast index.
-        """
-        
-        if fitType not in ('full', 'cope', 'pe'):
-            raise ValueError('Unknown model fit type {}'.format(fitType))
-        
-        TimeSeries.__init__(self, tsPanel, overlay, coords)
-        self.fitType  = fitType
-        self.idx      = idx
-        self.contrast = contrast
-        self.__updateModelFit()
-
-        
-    def update(self, coords):
-        """Overrides :meth:`TimeSeries.update`.
-
-        Updates the coordinates and the data encapsulated by this
-        ``FEATModelFitTimeSeries``.
-        """
-        if not TimeSeries.update(self, coords):
-            return
-        self.__updateModelFit()
-        
-
-    def __updateModelFit(self):
-        """Called by :meth:`update`, and in the constructor.  Updates the model
-        fit. See the :meth:`.FEATImage.fit` method.
-        """
-
-        fitType   = self.fitType
-        contrast  = self.contrast
-        xyz       = self.coords
-        self.data = self.overlay.fit(contrast, xyz, fitType == 'full')
-
-
-
-class MelodicTimeSeries(TimeSeries):
-
-    def __init__(self, tsPanel, overlay, component):
-        TimeSeries.__init__(self, tsPanel, overlay, component)
-
-
-    def _getData(self, component):
-        return self.overlay.getComponentTimeSeries(component)
+class TimeSeriesPanel(plotpanel.PlotPanel):
+    """A :class:`.PlotPanel` which plots time series data from :class:`.Image`
+    overlays. A ``TimeSeriesPanel`` looks something like the following:
 
+    .. image:: images/timeseriespanel.png
+       :scale: 50%
+       :align: center
 
-class TimeSeriesPanel(plotpanel.PlotPanel):
-    """A :class:`.PlotPanel` which plots time series data from
-    :class:`.Image` overlays.
 
     A ``TimeSeriesPanel`` plots one or more :class:`TimeSeries` instances,
     which encapsulate time series data for a voxel from a :class:`.Image`
@@ -761,19 +149,19 @@ class TimeSeriesPanel(plotpanel.PlotPanel):
     ================= =======================================================
     """
 
-    currentColour = copy.copy(TimeSeries.colour)
+    currentColour = copy.copy(timeseries.TimeSeries.colour)
     """Colour to use for the current time course. """
 
     
-    currentAlpha = copy.copy(TimeSeries.alpha)
+    currentAlpha = copy.copy(timeseries.TimeSeries.alpha)
     """Transparency to use for the current time course. """
 
     
-    currentLineWidth = copy.copy(TimeSeries.lineWidth)
+    currentLineWidth = copy.copy(timeseries.TimeSeries.lineWidth)
     """Line width to use for the current time course. """
 
     
-    currentLineStyle = copy.copy(TimeSeries.lineStyle)
+    currentLineStyle = copy.copy(timeseries.TimeSeries.lineStyle)
     """Line style to use for the current time course. """
 
 
@@ -883,7 +271,7 @@ class TimeSeriesPanel(plotpanel.PlotPanel):
             currOverlay = None
 
             if current is not None:
-                if isinstance(current, FEATTimeSeries):
+                if isinstance(current, timeseries.FEATTimeSeries):
                     extras = current.getModelTimeSeries()
                 else:
                     extras = [current]
@@ -912,7 +300,7 @@ class TimeSeriesPanel(plotpanel.PlotPanel):
                     ts.colour = colour
                     self.__overlayColours[ts.overlay] = colour
                         
-                    if isinstance(ts, FEATTimeSeries):
+                    if isinstance(ts, timeseries.FEATTimeSeries):
                         extras.extend(ts.getModelTimeSeries())
                 
             self.drawDataSeries(extras)
@@ -932,7 +320,7 @@ class TimeSeriesPanel(plotpanel.PlotPanel):
 
         tss = [self.__currentTs]
         
-        if isinstance(self.__currentTs, FEATTimeSeries):
+        if isinstance(self.__currentTs, timeseries.FEATTimeSeries):
             tss = self.__currentTs.getModelTimeSeries()
 
             for ts in tss:
@@ -1026,13 +414,13 @@ class TimeSeriesPanel(plotpanel.PlotPanel):
             return None
 
         if isinstance(overlay, fslfeatimage.FEATImage):
-            ts = FEATTimeSeries(self, overlay, loc)
+            ts = timeseries.FEATTimeSeries(self, overlay, loc)
             
         elif isinstance(overlay, fslmelimage.MelodicImage):
-            ts = MelodicTimeSeries(self, overlay, loc)
+            ts = timeseries.MelodicTimeSeries(self, overlay, loc)
             
         else:
-            ts = TimeSeries(self, overlay, loc)
+            ts = timeseries.TimeSeries(self, overlay, loc)
 
         ts.colour    = self.currentColour
         ts.alpha     = self.currentAlpha