featimage.py

#!/usr/bin/env python
#
# featimage.py - An Image subclass which has some FEAT-specific functionality.
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module provides the :class:`FEATImage` class, a subclass of
:class:`.Image` designed to encapsulate data from a FEAT analysis.
This module also provides the :func:`modelFit` function.
"""


import os.path      as op

import numpy        as np

from . import image as fslimage
from . import          featanalysis


class FEATImage(fslimage.Image):
    """An ``Image`` which contains the input data from a FEAT analysis.

    The :class:`FEATImage` class makes use of the functions defined in the
    :mod:`.featanalysis` module.


    An example of using the ``FEATImage`` class::

        import fsl.data.featimage as featimage

        # You can pass in the name of the
        # .feat directory, or the filtered_func_data
        # file contained within that directory.
        img = featimage.FEATImage('myanalysis.feat/filtered_func_data.nii.gz')

        # Query information about the FEAT analysis
        print(img.numEVs())
        print(img.contrastNames())
        print(img.numPoints())

        # Get the model fit residuals
        res4d = img.getResiduals()

        # Get the full model fit for voxel
        # [23, 30, 42] (in this example, we
        # have 4 EVs - the first argument
        # is a contrast vector).
        img.fit([1, 1, 1, 1], [23, 30, 42], fullModel=True)
    """


    def __init__(self, path, **kwargs):
        """Create a ``FEATImage`` instance.

        :arg path:   A FEAT analysis directory, or the input data image file
                     contained within such a directory.

        :arg kwargs: Passed to the :class:`.Image` constructor.
        """

        if op.isdir(path):
            path = op.join(path, 'filtered_func_data')

        if not featanalysis.isFEATImage(path):
            raise ValueError('{} does not appear to be data '
                             'from a FEAT analysis'.format(path))

        featDir     = op.dirname(path)
        settings    = featanalysis.loadSettings( featDir)

        if featanalysis.hasStats(featDir):
            design      = featanalysis.loadDesign(   featDir, settings)
            names, cons = featanalysis.loadContrasts(featDir)
        else:
            design      = None
            names, cons = [], []

        fslimage.Image.__init__(self, path, **kwargs)

        self.__analysisName  = op.splitext(op.basename(featDir))[0]
        self.__featDir       = featDir
        self.__design        = design
        self.__contrastNames = names
        self.__contrasts     = cons
        self.__settings      = settings

        self.__residuals     =  None
        self.__pes           = [None] * self.numEVs()
        self.__copes         = [None] * self.numContrasts()
        self.__zstats        = [None] * self.numContrasts()
        self.__clustMasks    = [None] * self.numContrasts()

        if 'name' not in kwargs:
            self.name = '{}: {}'.format(self.__analysisName, self.name)


    def __del__(self):
        """Clears references to any loaded images."""
        self.__design     = None
        self.__residuals  = None
        self.__pes        = None
        self.__copes      = None
        self.__zstats     = None
        self.__clustMasks = None


    def getFEATDir(self):
        """Returns the FEAT directory this image is contained in."""
        return self.__featDir


    def getAnalysisName(self):
        """Returns the FEAT analysis name, which is the FEAT directory
        name, minus the ``.feat`` / ``.gfeat`` suffix.
        """
        return self.__analysisName


    def isFirstLevelAnalysis(self):
        """Returns ``True`` if the FEAT analysis described by ``settings``
        is a first level analysis, ``False`` otherwise.
        """
        return featanalysis.isFirstLevelAnalysis(self.__settings)


    def getTopLevelAnalysisDir(self):
        """Returns the path to the higher level analysis directory of
        which this FEAT analysis is a part, or ``None`` if this analysis
        is not part of another analysis.
        """
        return featanalysis.getTopLevelAnalysisDir(self.__featDir)


    def getReportFile(self):
        """Returns the path to the FEAT report - see
        :func:`.featanalysis.getReportFile`.
        """
        return featanalysis.getReportFile(self.__featDir)


    def hasStats(self):
        """Returns ``True`` if the analysis for this ``FEATImage`` contains
        a statistical analysis.
        """
        return self.__design is not None


    def getDesign(self, voxel=None):
        """Returns the analysis design matrix as a :mod:`numpy` array
        with shape :math:`numPoints\\times numEVs`.
        See :meth:`.FEATFSFDesign.getDesign`.
        """

        if self.__design is None:
            return None

        return self.__design.getDesign(voxel)


    def numPoints(self):
        """Returns the number of points (e.g. time points, number of
        subjects, etc) in the analysis.
        """
        if self.__design is None:
            return 0

        return self.__design.getDesign().shape[0]


    def numEVs(self):
        """Returns the number of explanatory variables (EVs) in the analysis.
        """
        if self.__design is None:
            return 0

        return len(self.__design.getEVs())


    def evNames(self):
        """Returns a list containing the name of each EV in the analysis."""

        if self.__design is None:
            return []

        return [ev.title for ev in self.__design.getEVs()]


    def numContrasts(self):
        """Returns the number of contrasts in the analysis."""
        return len(self.__contrasts)


    def contrastNames(self):
        """Returns a list containing the name of each contrast in the analysis.
        """
        return list(self.__contrastNames)


    def contrasts(self):
        """Returns a list containing the analysis contrast vectors.

        See :func:`.featanalysis.loadContrasts`

        """
        return [list(c) for c in self.__contrasts]


    def thresholds(self):
        """Returns the statistical thresholds used in the analysis.

        See :func:`.featanalysis.getThresholds`
        """
        return featanalysis.getThresholds(self.__settings)


    def clusterResults(self, contrast):
        """Returns the clusters found in the analysis.

        See :func:.featanalysis.loadClusterResults`
        """
        return featanalysis.loadClusterResults(self.__featDir,
                                               self.__settings,
                                               contrast)


    def getPE(self, ev):
        """Returns the PE image for the given EV (0-indexed). """

        if self.__pes[ev] is None:
            pefile = featanalysis.getPEFile(self.__featDir, ev)
            self.__pes[ev] = fslimage.Image(
                pefile,
                name='{}: PE{} ({})'.format(
                    self.__analysisName,
                    ev + 1,
                    self.evNames()[ev]))

        return self.__pes[ev]


    def getResiduals(self):
        """Returns the residuals of the full model fit. """

        if self.__residuals is None:
            resfile = featanalysis.getResidualFile(self.__featDir)
            self.__residuals = fslimage.Image(
                resfile,
                name='{}: residuals'.format(self.__analysisName))

        return self.__residuals


    def getCOPE(self, con):
        """Returns the COPE image for the given contrast (0-indexed). """

        if self.__copes[con] is None:
            copefile = featanalysis.getCOPEFile(self.__featDir, con)
            self.__copes[con] = fslimage.Image(
                copefile,
                name='{}: COPE{} ({})'.format(
                    self.__analysisName,
                    con + 1,
                    self.contrastNames()[con]))
        return self.__copes[con]


    def getZStats(self, con):
        """Returns the Z statistic image for the given contrast (0-indexed).
        """

        if self.__zstats[con] is None:
            zfile = featanalysis.getZStatFile(self.__featDir, con)

            self.__zstats[con] = fslimage.Image(
                zfile,
                name='{}: zstat{} ({})'.format(
                    self.__analysisName,
                    con + 1,
                    self.contrastNames()[con]))

        return self.__zstats[con]


    def getClusterMask(self, con):
        """Returns the cluster mask image for the given contrast (0-indexed).
        """

        if self.__clustMasks[con] is None:
            mfile = featanalysis.getClusterMaskFile(self.__featDir, con)

            self.__clustMasks[con] = fslimage.Image(
                mfile,
                name='{}: cluster mask for zstat{} ({})'.format(
                    self.__analysisName,
                    con + 1,
                    self.contrastNames()[con]))

        return self.__clustMasks[con]


    def fit(self, contrast, xyz):
        """Calculates the model fit for the given contrast vector
        at the given voxel. See the :func:`modelFit` function.

        :arg contrast:  The contrast vector (pass all 1s for a full model
                        fit).

        :arg xyz:       Coordinates of the voxel to calculate the model fit
                        for.
        """

        if self.__design is None:
            raise RuntimeError('No design')

        x, y, z    = xyz
        firstLevel = self.isFirstLevelAnalysis()
        numEVs     = self.numEVs()

        if len(contrast) != numEVs:
            raise ValueError('Contrast is wrong length')

        design = self.getDesign(xyz)
        data   = self[x, y, z, :]
        pes    = [self.getPE(i)[x, y, z] for i in range(numEVs)]

        return modelFit(data, design, contrast, pes, firstLevel)


    def partialFit(self, contrast, xyz):
        """Calculates and returns the partial model fit for the specified
        contrast vector at the specified voxel.

        See :meth:`fit` for details on the arguments.
        """

        x, y, z   = xyz
        residuals = self.getResiduals()[x, y, z, :]
        modelfit  = self.fit(contrast, xyz)

        return residuals + modelfit


def modelFit(data, design, contrast, pes, firstLevel=True):
    """Calculates the model fit to the given data for the given contrast
    vector.

    :arg data:       The input data

    :arg design:     The design matrix

    :arg contrast:   The contrast vector (pass all 1s for a full model
                     fit)

    :arg pes:        Parameter estimates for each EV in the design matrix

    :arg firstLevel: If ``True`` (the default), the mean of the input
                     data is added to the result.

    :returns: The best fit of the model to the data.
    """

    # Here we are basically trying to
    # replicate the behaviour of tsplot.
    # There are some differences though -
    # by default, tsplot weights the
    # data by Z statistics. We're not
    # doing that here.

    # Normalise the contrast vector.
    # The scaling factor is arbitrary,
    # but should result in a visually
    # sensible scaling of the model fit.
    # For a vector of all 1's (i.e. a
    # full model fit) this is a no-op.
    #
    # We also take the absolute value
    # of the values in the contrast
    # vector, as the parameter estimates
    # should be appropriately signed,
    # so we don't want negative contrast
    # vector values to invert them.
    contrast = np.array(contrast)
    nevs     = len(contrast)
    nonzero  = sum(~np.isclose(contrast, 0))
    contrast = contrast / np.sqrt((contrast ** 2).sum())
    contrast = np.abs(contrast * np.sqrt(nonzero))

    modelfit = np.zeros(len(data))

    for i in range(nevs):

        ev        = design[:, i]
        pe        = pes[i]
        modelfit += ev * pe * contrast[i]

    # Make sure the model fit has an
    # appropriate mean.  The data in
    # first level analyses is demeaned
    # before model fitting, so we need
    # to add it back in.
    if firstLevel: return modelfit + data.mean()
    else:          return modelfit