Commit d17a1169 authored by Paul McCarthy's avatar Paul McCarthy 🚵
Browse files

Merge branch 'bf/powerspectrum' into 'master'

Bf/powerspectrum

Closes #189

See merge request fsl/fsleyes/fsleyes!224
parents 50edd483 20f371fd
......@@ -41,6 +41,14 @@ Changed
position of the light source.
* FSLeyes no longer ignores the ``LIBGL_ALWAYS_INDIRECT`` environment
variable.
* FSLeyes attempts to determine a suitable value for ``PYOPENGL_PLATFORM``
if it is not already set.
* FSLeyes should now work with both Wayland/EGL and X11/GLX builds of wxPython
on Linux.
* The normalisation method used in the power spectrum panel has been adjusted
so that, instead of the data being normalised to unit variance before the
fourier transform, the fourier-transformed data itself is normalised to the
range [-1, 1].
Fixed
......
......@@ -36,8 +36,7 @@ class DataSeries(props.HasProperties):
- Override the :meth:`redrawProperties` method if necessary
The overlay is accessible as an instance attribute, confusingly called
``overlay``.
The overlay is accessible as an instance attribute called ``overlay``.
.. note:: Some ``DataSeries`` instances may not be associated with
......@@ -212,6 +211,16 @@ class VoxelDataSeries(DataSeries):
It contains a built-in cache which is used to prevent repeated access
to data from the same voxel.
Sub-classes may need to override:
- the :meth:`currentVoxelLocation` method, which
generates the location index for the current location, and which
is used as the unique cache key when the corresponding data is cached
- The :meth:`currentVoxelData` method, which retrieves and/or calculates
the data at the current location. This is the data which is cached,
and which is returned by the :meth:`dataAtCurrentVoxel` method.
"""
......@@ -252,7 +261,10 @@ class VoxelDataSeries(DataSeries):
def getData(self):
"""Returns the data at the current voxel location. """
"""Returns the ``(xdata, ydata)`` at the current voxel location.
This method may be overridden by sub-classes.
"""
xdata = None
ydata = self.dataAtCurrentVoxel()
......@@ -275,19 +287,47 @@ class VoxelDataSeries(DataSeries):
# (which is a good thing).
@idle.mutex
def dataAtCurrentVoxel(self):
"""Returns the data for the current voxel of the overlay. The current
voxel is dictated by the :attr:`.DisplayContext.location` property.
This method is intended to be used by the :meth:`DataSeries.getData`
method of sub-classes.
"""Returns the data for the current voxel of the overlay. This method
is intended to be used within the :meth:`DataSeries.getData` method
of sub-classes.
An internal cache is used to avoid the need to retrieve data for the
same voxel multiple times, as retrieving data from large compressed
4D images can be time consuming.
It may also be overridden by sub-classes, but this implementation
should be called to take advantage of the voxel data cache.
The location for the current voxel is calculated by the
:meth:`currentVoxelLocation` method, and the data lookup is performed
by the :meth:`currentVoxelData` method. These methods may be
overridden by sub-classes.
:returns: A ``numpy`` array containing the data at the current
voxel, or ``None`` if the current location is out of bounds
of the image.
"""
location = self.currentVoxelLocation()
if location is None:
return None
data = self.__cache.get(location, None)
if data is None:
data = self.currentVoxelData(location)
self.__cache.put(location, data)
return data
def currentVoxelLocation(self):
"""Used by :meth:`dataAtCurrentVoxel`. Returns the current voxel
location. This is used as a key for the voxel data cache implemented
within the :meth:`dataAtCurrentVoxel` method, and subsequently passed
to the :meth:`currentVoxelData` method.
This method may be overridden by sub-classes.
"""
opts = self.displayCtx.getOpts(self.overlay)
vdim = opts.volumeDim
voxel = opts.getVoxel()
......@@ -297,10 +337,17 @@ class VoxelDataSeries(DataSeries):
x, y, z = voxel
data = self.__cache.get((x, y, z, vdim), None)
return (x, y, z, vdim)
if data is None:
data = self.overlay[opts.index(voxel, atVolume=False)]
self.__cache.put((x, y, z, vdim), data)
def currentVoxelData(self, location):
"""Used by :meth:`dataAtCurrentVoxel`. Returns the data at the
specified location.
This method may be overridden by sub-classes.
"""
voxel = location[:3]
opts = self.displayCtx.getOpts(self.overlay)
data = self.overlay[opts.index(voxel, atVolume=False)]
return data
......@@ -38,35 +38,16 @@ from . import dataseries
log = logging.getLogger(__name__)
def calcPowerSpectrum(data, varNorm=False):
"""Calculates a power spectrum for the given one-dimensional data array. If
``varNorm is True``, the data is de-meaned and normalised by its standard
deviation before the fourier transformation.
def calcPowerSpectrum(data):
"""Calculates a power spectrum for the given one-dimensional data array.
:arg data: Numpy array containing the time series data
:arg varNorm: Normalise the data before fourier transformation
:returns: If ``data`` contains real values, the magnitude of the power
spectrum is returned. If ``data`` contains complex values,
the complex power spectrum is returned.
"""
# De-mean and normalise
# by standard deviation
if varNorm:
mean = data.mean()
std = data.std()
if not np.isclose(std, 0):
data = data - mean
data = data / std
# If all values in the data
# are the same, it has mean 0
else:
data = np.zeros(data.shape, dtype=data.dtype)
# Fourier transform on complex data
if np.issubdtype(data.dtype, np.complexfloating):
data = fft.fft(data)
......@@ -76,7 +57,7 @@ def calcPowerSpectrum(data, varNorm=False):
# calculate and return the magnitude.
# We also drop the first (zero-frequency)
# term (see the rfft docs) as it is
# useless when varNorm is disabled.
# kind of useless for display purposes
else:
data = fft.rfft(data)[1:]
data = magnitude(data)
......@@ -84,19 +65,19 @@ def calcPowerSpectrum(data, varNorm=False):
return data
def calcFrequencies(data, sampleTime):
def calcFrequencies(nsamples, sampleTime, dtype):
"""Calculates the frequencies of the power spectrum for the given
data.
:arg data: The input time series data
:arg nsamples: Number of samples in the input time series data
:arg sampleTime: Time between each data point
:arg dtype: Data type - the calculation differs depending on
whether the data is real or complex.
:returns: A ``numpy`` array containing the frequencies of the
power spectrum for ``data``
"""
nsamples = len(data)
if np.issubdtype(data.dtype, np.complexfloating):
if np.issubdtype(dtype, np.complexfloating):
xdata = fft.fftfreq(nsamples, sampleTime)
xdata = fft.fftshift(xdata)
else:
......@@ -121,6 +102,17 @@ def phase(data):
return np.arctan2(imag, real)
def normalise(data, dmin=None, dmax=None, nmin=-1, nmax=1):
"""Returns ``data``, rescaled to the range [nmin, nmax].
If dmin and dmax are provided, the data is normalised with respect to them,
rather than being normalised by the data minimum/maximum.
"""
if dmin is None: dmin = data.min()
if dmax is None: dmax = data.max()
return nmin + (nmax - nmin) * (data - dmin) / (dmax - dmin)
def phaseCorrection(spectrum, freqs, p0, p1):
"""Applies phase correction to the given complex power spectrum.
......@@ -134,16 +126,19 @@ def phaseCorrection(spectrum, freqs, p0, p1):
return np.exp(exp) * spectrum
class PowerSpectrumSeries(object):
class PowerSpectrumSeries:
"""The ``PowerSpectrumSeries`` encapsulates a power spectrum data series
from an overlay. The ``PowerSpectrumSeries`` class is a base mixin class
for all other classes in this module.
"""
varNorm = props.Boolean(default=True)
"""If ``True``, the data is normalised to unit variance before the fourier
transformation.
varNorm = props.Boolean(default=False)
"""If ``True``, the fourier-transformed data is normalised to the range
[0, 1] before plotting.
.. note:: The :class:`ComplexPowerSpectrumSeries` applies normalisation
differently.
"""
......@@ -179,16 +174,29 @@ class VoxelPowerSpectrumSeries(dataseries.VoxelDataSeries,
raise ValueError('Overlay is not a 4D image')
def getData(self):
"""Returns the data at the current voxel. """
def currentVoxelData(self, location):
"""Overrides :meth:`.VoxelDataSeries.currentVoxelData`. Retrieves
the data at the specified location, then performs a fourier transform
on it and returnes the result.
"""
data = self.dataAtCurrentVoxel()
data = dataseries.VoxelDataSeries.currentVoxelData(self, location)
data = calcPowerSpectrum(data)
return data
if data is None:
return None, None
xdata = calcFrequencies( data, self.sampleTime)
ydata = calcPowerSpectrum(data, self.varNorm)
def getData(self):
"""Returns the ``(xdata, ydata)`` to be plotted from the current voxel
location.
"""
overlay = self.overlay
ydata = self.dataAtCurrentVoxel()
xdata = calcFrequencies(overlay.shape[3],
self.sampleTime,
overlay.dtype)
if self.varNorm:
ydata = normalise(ydata)
return xdata, ydata
......@@ -224,12 +232,14 @@ class ComplexPowerSpectrumSeries(VoxelPowerSpectrumSeries):
VoxelPowerSpectrumSeries.__init__(
self, overlay, overlayList, displayCtx, plotPanel)
self.__cachedData = (None, None)
self.__imagps = ImaginaryPowerSpectrumSeries(
# Separate DataSeries for the imaginary/
# magnitude/phase signals, returned by
# the extraSeries method
self.__imagps = ImaginaryPowerSpectrumSeries(
self, overlay, overlayList, displayCtx, plotPanel)
self.__magps = MagnitudePowerSpectrumSeries(
self.__magps = MagnitudePowerSpectrumSeries(
self, overlay, overlayList, displayCtx, plotPanel)
self.__phaseps = PhasePowerSpectrumSeries(
self.__phaseps = PhasePowerSpectrumSeries(
self, overlay, overlayList, displayCtx, plotPanel)
for ps in (self.__imagps, self.__magps, self.__phaseps):
......@@ -239,35 +249,60 @@ class ComplexPowerSpectrumSeries(VoxelPowerSpectrumSeries):
ps.bindProps('lineStyle', self)
def makeLabel(self):
"""Returns a string representation of this
``ComplexPowerSpectrumSeries`` instance.
def makeLabelBase(self):
"""Returns a string to be used as the label prefix for this
``ComplexPowerSpectrumSeries`` instance, and for the imaginary,
magnitude, and phase child series.
"""
return '{} ({})'.format(VoxelPowerSpectrumSeries.makeLabel(self),
strings.labels[self])
return VoxelPowerSpectrumSeries.makeLabel(self)
@property
def cachedData(self):
"""Returns the currently cached data (see :meth:`getData`). """
return self.__cachedData
def makeLabel(self):
"""Returns a label to use for this data series. """
return '{} ({})'.format(self.makeLabelBase(), strings.labels[self])
def getData(self):
def getData(self, component='real'):
"""If :attr:`plotReal` is true, returns the real component of the power
spectrum of the data at the current voxel. Otherwise returns ``(None,
None)``.
Every time this method is called, the power spectrum is calculated,
phase correction is applied, and a reference to the resulting complex
power spectrum (and frequencies) is saved; it is accessible via the
:meth:`cachedData` property, for use by the
:class:`ImaginaryPowerSpectrumSeries`,
:class:`MagnitudePowerSpectrumSeries`, and
:class:`PhasePowerSpectrumSeries`.
Every time this method is called, the power spectrum is retrieved (see
the :class:`VoxelPowerSpectrumSeries` class), phase correction is
applied if set, andthe data is normalised, if set. A tuple containing
the ``(xdata, ydata)`` is returned, with ``ydata`` containing the
requested ``component`` ( ``'real'``, ``'imaginary'``,
``'magnitude'``, or ``'phase'``).
This method is called by the :class:`ImaginarySpectrumPowerSeries`,
:class:`MagnitudeSpectrumPowerSeries`, and
:class:`PhasePowerSpectrumPowerSeries` instances that are associated
with this data series.
"""
xdata, ydata = VoxelPowerSpectrumSeries.getData(self)
if ((component == 'real') and (not self.plotReal)) or \
((component == 'imaginary') and (not self.plotImaginary)) or \
((component == 'magnitude') and (not self.plotMagnitude)) or \
((component == 'phase') and (not self.plotPhase)):
return None, None
# See VoxelPowerSpectrumSeries - the data
# is already fourier-transformed
ydata = self.dataAtCurrentVoxel()
if ydata is None:
return None, None
# All of the calculations below are repeated
# for each real/imag/mag/phase series that
# gets plotted. But keeping the code together
# and clean is currently more important than
# performance, as there is not really any
# performance hit.
overlay = self.overlay
xdata = calcFrequencies(overlay.shape[3],
self.sampleTime,
overlay.dtype)
if self.zeroOrderPhaseCorrection != 0 or \
self.firstOrderPhaseCorrection != 0:
......@@ -276,19 +311,21 @@ class ComplexPowerSpectrumSeries(VoxelPowerSpectrumSeries):
self.zeroOrderPhaseCorrection,
self.firstOrderPhaseCorrection)
# Note that we're assuming that this
# ComplexPowerSpectrumSeries.getData
# method will be called before the
# corresponding call(s) to the
# Imaginary/Magnitude/Phase series
# methods.
self.__cachedData = xdata, ydata
if not self.plotReal:
return None, None
if ydata is not None:
ydata = ydata.real
# Normalise magnitude, real, imaginary
# components with respect to magnitude.
# Normalise phase independently.
if self.varNorm:
mag = magnitude(ydata)
mr = mag.min(), mag.max()
if component == 'phase': ydata = normalise(phase(ydata))
elif component == 'magnitude': ydata = normalise(mag)
elif component == 'real': ydata = normalise(ydata.real, *mr)
elif component == 'imaginary': ydata = normalise(ydata.imag, *mr)
elif component == 'real': ydata = ydata.real
elif component == 'imaginary': ydata = ydata.imag
elif component == 'magnitude': ydata = magnitude(ydata)
elif component == 'phase': ydata = phase(ydata)
return xdata, ydata
......@@ -331,19 +368,13 @@ class ImaginaryPowerSpectrumSeries(dataseries.DataSeries):
"""Returns a string representation of this
``ImaginaryPowerSpectrumSeries`` instance.
"""
return '{} ({})'.format(self.__parent.makeLabel(),
return '{} ({})'.format(self.__parent.makeLabelBase(),
strings.labels[self])
def getData(self):
"""Returns the imaginary component of the power spectrum. """
xdata, ydata = self.__parent.cachedData
if ydata is not None:
ydata = ydata.imag
return xdata, ydata
return self.__parent.getData('imaginary')
class MagnitudePowerSpectrumSeries(dataseries.DataSeries):
......@@ -370,16 +401,13 @@ class MagnitudePowerSpectrumSeries(dataseries.DataSeries):
"""Returns a string representation of this
``MagnitudePowerSpectrumSeries`` instance.
"""
return '{} ({})'.format(self.__parent.makeLabel(),
return '{} ({})'.format(self.__parent.makeLabelBase(),
strings.labels[self])
def getData(self):
"""Returns the magnitude of the complex power spectrum. """
xdata, ydata = self.__parent.cachedData
if ydata is not None:
ydata = magnitude(ydata)
return xdata, ydata
return self.__parent.getData('magnitude')
class PhasePowerSpectrumSeries(dataseries.DataSeries):
......@@ -406,16 +434,13 @@ class PhasePowerSpectrumSeries(dataseries.DataSeries):
"""Returns a string representation of this ``PhasePowerSpectrumSeries``
instance.
"""
return '{} ({})'.format(self.__parent.makeLabel(),
return '{} ({})'.format(self.__parent.makeLabelBase(),
strings.labels[self])
def getData(self):
"""Returns the phase of the complex power spectrum. """
xdata, ydata = self.__parent.cachedData
if ydata is not None:
ydata = phase(ydata)
return xdata, ydata
return self.__parent.getData('phase')
class MelodicPowerSpectrumSeries(dataseries.DataSeries,
......@@ -524,7 +549,7 @@ class MeshPowerSpectrumSeries(dataseries.DataSeries,
vd = opts.getVertexData()
data = vd[vidx, :]
xdata = calcFrequencies( data, self.sampleTime)
ydata = calcPowerSpectrum(data, self.varNorm)
xdata = calcFrequencies( len(data), self.sampleTime, data.dtype)
ydata = calcPowerSpectrum(data)
return xdata, ydata
......@@ -1065,7 +1065,7 @@ properties = TypeDict({
'ComplexHistogramSeries.plotMagnitude' : 'Plot magnitude',
'ComplexHistogramSeries.plotPhase' : 'Plot phase',
'PowerSpectrumSeries.varNorm' : 'Normalise to unit variance',
'PowerSpectrumSeries.varNorm' : 'Normalise to [-1, 1]',
'FEATTimeSeries.plotFullModelFit' : 'Plot full model fit',
'FEATTimeSeries.plotEVs' : 'Plot EV{} ({})',
......
......@@ -702,15 +702,14 @@ properties = TypeDict({
'histogram.',
'PowerSpectrumSeries.varNorm' : 'If checked, the data is demeaned and '
'normalised by its standard deviation '
'before its power spectrum is '
'calculated via a fourier transform.',
'PowerSpectrumSeries.varNorm' :
'If checked, the fourier-transformed data is normalised to the range '
'[-1, 1]. Complex valued data are normalised with respect to the '
'absolute value. ',
# Profiles
'OrthoPanel.profile' : 'Switch between view mode '
'and edit mode',
'OrthoPanel.profile' :
'Switch between view mode and edit mode',
'OrthoEditProfile.selectionCursorColour' :
'Colour to use for the selection cursor.',
......
......@@ -37,13 +37,9 @@ def test_calcPowerSpectrum():
]
for data, explen, expdtype in tests:
got1 = psseries.calcPowerSpectrum(data)
got2 = psseries.calcPowerSpectrum(data, True)
assert got1.shape == (explen,)
assert got2.shape == (explen,)
assert np.issubdtype(got1.dtype, expdtype)
assert np.issubdtype(got2.dtype, expdtype)
got = psseries.calcPowerSpectrum(data)
assert got.shape == (explen,)
assert np.issubdtype(got.dtype, expdtype)
def test_calcFrequencies():
......@@ -57,7 +53,7 @@ def test_calcFrequencies():
]
for data, explen in tests:
got = psseries.calcFrequencies(data, 1)
got = psseries.calcFrequencies(len(data), 1, data.dtype)
assert got.shape == (explen,)
......@@ -68,11 +64,14 @@ def test_magnitude():
def test_phase():
assert np.isclose(psseries.phase(4 + 4j), np.pi / 4)
def test_normalise():
data = np.array([1, 2, 3, 4, 5])
assert np.all(psseries.normalise(data) == np.linspace(-1, 1, 5))
def test_phaseCorrection():
data = np.random.random(100) + np.random.random(100) * 1j
spectrum = psseries.calcPowerSpectrum(data)
freqs = psseries.calcFrequencies(data, 1)
freqs = psseries.calcFrequencies(len(data), 1, data.dtype)
corrected = psseries.phaseCorrection(spectrum, freqs, 1, 2)
# TODO write a real test
......@@ -94,14 +93,22 @@ def _test_VoxelPowerSpectrumSeries(panel, overlayList, displayCtx):
displayCtx.location = loc
realYield()
expx = psseries.calcFrequencies( img[x, y, z, :], img.pixdim[3])
expy = psseries.calcPowerSpectrum(img[x, y, z, :], True)
expx = psseries.calcFrequencies(img.shape[3], img.pixdim[3], img.dtype)
expy = psseries.calcPowerSpectrum(img[x, y, z, :])
xdata, ydata = ps.getData()
assert ps.sampleTime == img.pixdim[3]
assert np.all(xdata == expx)
assert np.all(ydata == expy)
ps.varNorm = True
expy = psseries.normalise(expy)
ydata = ps.getData()[1]
assert np.all(ydata == expy)
ps.varNorm = False
x = x - 1
y = y + 1
z = z - 1
......@@ -110,8 +117,8 @@ def _test_VoxelPowerSpectrumSeries(panel, overlayList, displayCtx):
displayCtx.location = loc
realYield()
xdata, ydata = ps.getData()
expx = psseries.calcFrequencies( img[x, y, z, :], img.pixdim[3])
expy = psseries.calcPowerSpectrum(img[x, y, z, :], True)
expx = psseries.calcFrequencies(img.shape[3], img.pixdim[3], img.dtype)
expy = psseries.calcPowerSpectrum(img[x, y, z, :])
assert np.all(xdata == expx)
assert np.all(ydata == expy)
......@@ -132,8 +139,8 @@ def _test_ComplexPowerSpectrumSeries(panel, overlayList, displayCtx):
ps = panel.getDataSeries(img)
displayCtx.location = opts.transformCoords((7, 8, 9), 'voxel', 'display')
expx = psseries.calcFrequencies( img[7, 8, 9, :], 1)
expy = psseries.calcPowerSpectrum(img[7, 8, 9, :], True)
expx = psseries.calcFrequencies( img.shape[3], 1, img.dtype)
expy = psseries.calcPowerSpectrum(img[7, 8, 9, :])
xdata, ydata = ps.getData()
assert np.all(xdata == expx)
......@@ -213,8 +220,8 @@ def _test_MeshPowerSpectrumSeries(panel, overlayList, displayCtx):
mesh.vertices[10, :], 'mesh', 'display')
realYield()
exp = data[10, :]
expx = psseries.calcFrequencies(exp, 1)
expy = psseries.calcPowerSpectrum(exp, True)
expx = psseries.calcFrequencies(len(exp), 1, exp.dtype)
expy = psseries.calcPowerSpectrum(exp)
ps = panel.getDataSeries(mesh)
gotx, goty = ps.getData()
......
Supports Markdown
0%