doc/index.rst

@@ -4,12 +4,57 @@
 
 
 The ``fslpy`` package is a collection of utilities and data abstractions used
-by |fsleyes_apidoc|_.
+within `FSL <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki>`_ and by
+|fsleyes_apidoc|_.
+
+
+The top-level Python package for ``fslpy`` is called :mod:`fsl`.  It is
+broadly split into the following sub-packages:
+
++----------------------+-----------------------------------------------------+
+| :mod:`fsl.data`      | contains data abstractions and I/O routines for a   |
+|                      | range of FSL and neuroimaging file types. Most I/O  |
+|                      | routines use `nibabel <https://nipy.org/nibabel/>`_ |
+|                      | extensively.                                        |
++----------------------+-----------------------------------------------------+
+| :mod:`fsl.utils`     | contains a range of miscellaneous utilities,        |
+|                      | including :mod:`fsl.utils.path`,                    |
+|                      | :mod:`fsl.utils.run`, and :mod:`fsl.utils.bids`     |
++----------------------+-----------------------------------------------------+
+| :mod:`fsl.scripts`   | contains a range of scripts which are installed as  |
+|                      | FSL commands.                                       |
++----------------------+-----------------------------------------------------+
+| :mod:`fsl.transform` | contains functions and classes for working with     |
+|                      | FSL-style linear and non-linear transformations.    |
++----------------------+-----------------------------------------------------+
+| :mod:`fsl.version`   | simply contains the ``fslpy`` version number.       |
++----------------------+-----------------------------------------------------+
+| :mod:`fsl.wrappers`  | contains Python functions which can be used to      |
+|                      | invoke FSL commands.                                |
++----------------------+-----------------------------------------------------+
+
+The :mod:`fsl` package provides the top-level Python package namespace for
+``fslpy``, and for other FSL python libaries. It is a `native namespace
+package <https://packaging.python.org/guides/packaging-namespace-packages/>`_,
+which means that there is no ``fsl/__init__.py`` file.
+
+
+Other libraries can use the ``fsl`` package namepace simply by also omitting a
+``fsl/__init__.py`` file, and by ensuring that there are no naming conflicts
+with any sub-packages of ``fslpy`` or any other projects which use the ``fsl``
+package namespace.
+
+
 
 .. toctree::
    :hidden:
 
    self
-   fsl
+   fsl.data
+   fsl.scripts
+   fsl.transform
+   fsl.utils
+   fsl.wrappers
+   fsl.version
    contributing
    changelog

doc/mock_modules.txt

 deprecation
+dill
+h5py
 nibabel
+nibabel.cifti2
 nibabel.fileslice
 nibabel.freesurfer
 numpy
 numpy.linalg
 scipy
 scipy.ndimage
+scipy.ndimage.interpolation
+six

fsl/__init__.py

-#!/usr/bin/env python
-#
-# __init__.py - The fslpy library.
-#
-# Author: Paul McCarthy <pauldmccarthy@gmail.com>
-#
-"""The :mod:`fsl` package is a library which contains convenience classes
-and functions for use by FSL python tools. It is broadly split into the
-following sub-packages:
-
-.. autosummary::
-
-   fsl.data
-   fsl.utils
-   fsl.scripts
-   fsl.version
-   fsl.wrappers
-
-.. note:: The ``fsl`` namespace is a ``pkgutil``-style *namespace package* -
-          it can be used across different projects - see
-          https://packaging.python.org/guides/packaging-namespace-packages/
-          for details.
-"""
-
-__path__ = __import__('pkgutil').extend_path(__path__, __name__)

fsl/data/bitmap.py

+#!/usr/bin/env python
+#
+# bitmap.py - The Bitmap class
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains the :class:`Bitmap` class, for loading bitmap image
+files. Pillow is required to use the ``Bitmap`` class.
+"""
+
+
+import os.path as op
+import            pathlib
+import            logging
+
+import numpy as np
+
+import fsl.data.image as fslimage
+
+
+log = logging.getLogger(__name__)
+
+
+BITMAP_EXTENSIONS = ['.bmp', '.png',  '.jpg', '.jpeg',
+                     '.tif', '.tiff', '.gif', '.rgba',
+                     '.jp2', '.jpg2', '.jp2k']
+"""File extensions we understand. """
+
+
+BITMAP_DESCRIPTIONS = [
+    'Bitmap',
+    'Portable Network Graphics',
+    'JPEG',
+    'JPEG',
+    'TIFF',
+    'TIFF',
+    'Graphics Interchange Format',
+    'Raw RGBA',
+    'JPEG 2000',
+    'JPEG 2000',
+    'JPEG 2000']
+"""A description for each :attr:`BITMAP_EXTENSION`. """
+
+
+class Bitmap(object):
+    """The ``Bitmap`` class can be used to load a bitmap image. The
+    :meth:`asImage` method will convert the bitmap into an :class:`.Image`
+    instance.
+    """
+
+    def __init__(self, bmp):
+        """Create a ``Bitmap``.
+
+        :arg bmp: File name of an image, or a ``numpy`` array containing image
+                  data.
+        """
+
+        if isinstance(bmp, (pathlib.Path, str)):
+
+            try:
+                # Allow big/truncated images
+                import PIL.Image     as Image
+                import PIL.ImageFile as ImageFile
+                Image    .MAX_IMAGE_PIXELS      = None
+                ImageFile.LOAD_TRUNCATED_IMAGES = True
+
+            except ImportError:
+                raise RuntimeError('Install Pillow to use the Bitmap class')
+
+            src = str(bmp)
+            img = Image.open(src)
+
+            # If this is a palette/LUT
+            # image, convert it into a
+            # regular rgb(a) image.
+            if img.mode == 'P':
+                img = img.convert()
+
+            data = np.array(img)
+
+        elif isinstance(bmp, np.ndarray):
+            src  = 'array'
+            data = np.copy(bmp)
+
+        else:
+            raise ValueError('unknown bitmap: {}'.format(bmp))
+
+        # Make the array (w, h, c). Single channel
+        # (e.g. greyscale) images are returned as
+        # 2D arrays, whereas multi-channel images
+        # are returned as 3D. In either case, the
+        # first two dimensions are (height, width),
+        # but we watn them the other way aruond.
+        data = np.atleast_3d(data)
+        data = np.fliplr(data.transpose((1, 0, 2)))
+        data = np.array(data, dtype=np.uint8, order='C')
+        w, h = data.shape[:2]
+
+        self.__data       = data
+        self.__dataSource = src
+        self.__name       = op.basename(src)
+
+
+    def __hash__(self):
+        """Returns a number which uniquely idenfities this ``Bitmap`` instance
+        (the result of ``id(self)``).
+        """
+        return id(self)
+
+
+    def __str__(self):
+        """Return a string representation of this ``Bitmap`` instance."""
+        return '{}({}, {})'.format(self.__class__.__name__,
+                                   self.dataSource,
+                                   self.shape)
+
+
+    def __repr__(self):
+        """See the :meth:`__str__` method. """
+        return self.__str__()
+
+
+    @property
+    def name(self):
+        """Returns the name of this ``Bitmap``, typically the base name of the
+        file.
+        """
+        return self.__name
+
+
+    @property
+    def dataSource(self):
+        """Returns the bitmap data source - typically the file name. """
+        return self.__dataSource
+
+
+    @property
+    def data(self):
+        """Convenience method which returns the bitmap data as a ``(w, h, c)``
+        array, where ``c`` is either 3 or 4.
+        """
+        return self.__data
+
+
+    @property
+    def shape(self):
+        """Returns the bitmap shape - ``(width, height, nchannels)``. """
+        return self.__data.shape
+
+
+    def asImage(self):
+        """Convert this ``Bitmap`` into an :class:`.Image` instance. """
+
+        width, height, nchannels = self.shape
+
+        if nchannels == 1:
+            dtype = np.uint8
+
+        elif nchannels == 3:
+            dtype = np.dtype([('R', 'uint8'),
+                              ('G', 'uint8'),
+                              ('B', 'uint8')])
+
+        elif nchannels == 4:
+            dtype = np.dtype([('R', 'uint8'),
+                              ('G', 'uint8'),
+                              ('B', 'uint8'),
+                              ('A', 'uint8')])
+
+        else:
+            raise ValueError('Cannot convert bitmap with {} '
+                             'channels into nifti image'.format(nchannels))
+
+        if nchannels == 1:
+            data = self.data.reshape((width, height))
+
+        else:
+            data = np.zeros((width, height), dtype=dtype)
+            for ci, ch in enumerate(dtype.names):
+                data[ch] = self.data[..., ci]
+
+        data = np.asarray(data, order='F')
+
+        return fslimage.Image(data,
+                              name=self.name,
+                              dataSource=self.dataSource)

fsl/data/constants.py

@@ -30,6 +30,7 @@ specification:
    NIFTI_XFORM_ALIGNED_ANAT
    NIFTI_XFORM_TALAIRACH
    NIFTI_XFORM_MNI_152
+   NIFTI_XFORM_TEMPLATE_OTHER
 """
 
 
@@ -81,7 +82,14 @@ NIFTI_XFORM_MNI_152      = 4
 """MNI 152 normalized coordinates."""
 
 
-NIFTI_XFORM_ANALYZE      = 5
+NIFTI_XFORM_TEMPLATE_OTHER = 5
+"""Coordinates aligned to some template that is not MNI152 or Talairach.
+
+See https://www.nitrc.org/forum/message.php?msg_id=26394 for details.
+ """
+
+
+NIFTI_XFORM_ANALYZE      = 6
 """Code which indicates that this is an ANALYZE image, not a NIFTI image. """
 
 
@@ -98,6 +106,36 @@ NIFTI_UNITS_PPM     = 40
 NIFTI_UNITS_RADS    = 48
 
 
+# NIFTI datatype codes
+NIFTI_DT_NONE          = 0
+NIFTI_DT_UNKNOWN       = 0
+NIFTI_DT_BINARY        = 1
+NIFTI_DT_UNSIGNED_CHAR = 2
+NIFTI_DT_SIGNED_SHORT  = 4
+NIFTI_DT_SIGNED_INT    = 8
+NIFTI_DT_FLOAT         = 16
+NIFTI_DT_COMPLEX       = 32
+NIFTI_DT_DOUBLE        = 64
+NIFTI_DT_RGB           = 128
+NIFTI_DT_ALL           = 255
+NIFTI_DT_UINT8         = 2
+NIFTI_DT_INT16         = 4
+NIFTI_DT_INT32         = 8
+NIFTI_DT_FLOAT32       = 16
+NIFTI_DT_COMPLEX64     = 32
+NIFTI_DT_FLOAT64       = 64
+NIFTI_DT_RGB24         = 128
+NIFTI_DT_INT8          = 256
+NIFTI_DT_UINT16        = 512
+NIFTI_DT_UINT32        = 768
+NIFTI_DT_INT64         = 1024
+NIFTI_DT_UINT64        = 1280
+NIFTI_DT_FLOAT128      = 1536
+NIFTI_DT_COMPLEX128    = 1792
+NIFTI_DT_COMPLEX256    = 2048
+NIFTI_DT_RGBA32        = 2304
+
+
 # NIFTI file intent codes
 NIFTI_INTENT_NONE          = 0
 NIFTI_INTENT_CORREL        = 2

fsl/data/dicom.py

@@ -29,22 +29,56 @@ import               os
 import os.path    as op
 import subprocess as sp
 import               re
+import               sys
 import               glob
 import               json
+import               shlex
+import               shutil
 import               logging
+import               binascii
 
+import numpy      as np
 import nibabel    as nib
 
-import fsl.utils.tempdir as tempdir
-import fsl.utils.memoize as memoize
-import fsl.data.image    as fslimage
+import fsl.utils.tempdir  as tempdir
+import fsl.utils.memoize  as memoize
+import fsl.utils.platform as fslplatform
+import fsl.data.image     as fslimage
 
 
 log = logging.getLogger(__name__)
 
 
 MIN_DCM2NIIX_VERSION = (1, 0, 2017, 12, 15)
-"""Minimum version of dcm2niix that is required for this module to work. """
+"""Minimum version of ``dcm2niix`` that is required for this module to work.
+"""
+
+
+CRC_DCM2NIIX_VERSION = (1, 0, 2019, 9, 2)
+"""For versions of ``dcm2niix`` orf this version or newer, the ``-n`` flag,
+used to convert a single DICOM series, requires that a CRC checksum
+identifying the series be passed (see the :func:`seriesCRC`
+function). Versions prior to this require the series number to be passed.
+"""
+
+
+def dcm2niix() -> str:
+    """Tries to find an absolute path to the ``dcm2niix`` command. Returns
+    ``'dcm2niix'`` (unqualified) if a specific executable cannot be found.
+    """
+    fsldir     = fslplatform.platform.fsldir
+    candidates = [
+        shutil.which('dcm2niix')
+    ]
+
+    if fsldir is not None:
+        candidates.insert(0, op.join(fsldir, 'bin', 'dcm2niix'))
+
+    for c in candidates:
+        if c is not None and op.exists(c):
+            return c
+
+    return 'dcm2niix'
 
 
 class DicomImage(fslimage.Image):
@@ -80,12 +114,19 @@ class DicomImage(fslimage.Image):
 
 
 @memoize.memoize
-def enabled():
-    """Returns ``True`` if ``dcm2niix`` is present, and recent enough,
-    ``False`` otherwise.
+def installedVersion():
+    """Return a tuple describing the version of ``dcm2niix`` that is installed,
+    or ``None`` if dcm2niix cannot be found, or its version not parsed.
+
+    The returned tuple contains the following fields, all integers:
+      - Major version number
+      - Minor version number
+      - Year
+      - Month
+      - Day
     """
 
-    cmd            = 'dcm2niix -h'
+    cmd            = f'{dcm2niix()} -h'
     versionPattern = re.compile(r'v'
                                 r'(?P<major>[0-9]+)\.'
                                 r'(?P<minor>[0-9]+)\.'
@@ -102,80 +143,124 @@ def enabled():
 
             match = re.match(versionPattern, word)
 
-            if match is None:
-                continue
+            if match is not None:
+                return (int(match.group('major')),
+                        int(match.group('minor')),
+                        int(match.group('year')),
+                        int(match.group('month')),
+                        int(match.group('day')))
+
+    except Exception as e:
+        log.debug(f'Error parsing dcm2niix version string: {e}')
+    return None
 
-            installedVersion = (
-                int(match.group('major')),
-                int(match.group('minor')),
-                int(match.group('year')),
-                int(match.group('month')),
-                int(match.group('day')))
 
-            # make sure installed version
-            # is equal to or newer than
-            # minimum required version
-            for iv, mv in zip(installedVersion, MIN_DCM2NIIX_VERSION):
-                if   iv > mv: return True
-                elif iv < mv: return False
+def compareVersions(v1, v2):
+    """Compares two ``dcm2niix`` versions ``v1`` and ``v2``.  The versions are
+    assumed to be in the format returned by :func:`installedVersion`.
 
-            # if we get here, versions are equal
-            return True
+    :returns: - 1 if ``v1`` is newer than ``v2``
+              - -1 if ``v1`` is older than ``v2``
+              - 0 if ``v1`` the same as ``v2``.
+    """
 
-    except Exception as e:
-        log.debug('Error parsing dcm2niix version string: {}'.format(e))
+    for iv1, iv2 in zip(v1, v2):
+        if   iv1 > iv2: return  1
+        elif iv1 < iv2: return -1
+    return 0
 
-    return False
+
+def enabled():
+    """Returns ``True`` if ``dcm2niix`` is present, and recent enough,
+    ``False`` otherwise.
+    """
+    installed = installedVersion()
+    required  = MIN_DCM2NIIX_VERSION
+    return ((installed is not None) and
+            (compareVersions(installed, required) >= 0))
 
 
 def scanDir(dcmdir):
-    """Uses ``dcm2niix`` to scans the given DICOM directory, and returns a
-    list of dictionaries, one for each data series that was identified.
-    Each dictionary is populated with some basic metadata about the series.
+    """Uses the ``dcm2niix -b o`` option to generate a BIDS sidecar JSON
+    file for each series in the given DICOM directory. Reads them all in,
+    and returns them as a sequence of dicts.
 
-    :arg dcmdir: Directory containing DICOM files.
+    Some additional metadata is added to each dictionary:
+     - ``DicomDir``: The absolute path to ``dcmdir``
 
-    :returns:    A list of dictionaries, each containing metadata about
-                 one DICOM data series.
+    :arg dcmdir: Directory containing DICOM series
+
+    :returns:    A list of dicts, each containing the BIDS sidecar JSON
+                 metadata for one DICOM series.
     """
 
     if not enabled():
         raise RuntimeError('dcm2niix is not available or is too old')
 
-    dcmdir      = op.abspath(dcmdir)
-    cmd         = 'dcm2niix -b o -ba n -f %s -o . {}'.format(dcmdir)
-    snumPattern = re.compile('^[0-9]+')
+    dcmdir = op.abspath(dcmdir)
+    cmd    = f'{dcm2niix()} -b o -ba n -f %s -o . "{dcmdir}"'
+    series = []
 
     with tempdir.tempdir() as td:
 
         with open(os.devnull, 'wb') as devnull:
-            sp.call(cmd.split(), stdout=devnull, stderr=devnull)
+            sp.call(shlex.split(cmd), stdout=devnull, stderr=devnull)
 
         files = glob.glob(op.join(td, '*.json'))
 
         if len(files) == 0:
             return []
 
-        # sort numerically by series number if possible
-        try:
-            def sortkey(f):
-                match = re.match(snumPattern, f)
-                snum  = int(match.group(0))
-                return snum
-
-            files = sorted(files, key=sortkey)
-
-        except Exception:
-            files = sorted(files)
-
-        series = []
         for fn in files:
             with open(fn, 'rt') as f:
-                meta = json.load(f)
+                meta             = json.load(f)
                 meta['DicomDir'] = dcmdir
+                # SeriesDescription is not
+                # guaranteed to be present
+                if 'SeriesDescription' not in meta:
+                    meta['SeriesDescription'] = meta['SeriesNumber']
                 series.append(meta)
 
-        return series
+    # sort by series number
+    def key(s):
+        return s.get('SeriesNumber', sys.maxsize)
+
+    series = list(sorted(series, key=key))
+
+    return series
+
+
+def seriesCRC(series):
+    """Calculate a checksum string of the given DICOM series.
+
+    The returned string is of the form::
+
+         SeriesCRC[.echonumber]
+
+    Where ``SeriesCRC`` is an unsigned integer which is the CRC32
+    checksum of the ``SeriesInstanceUID``, and ``echonumber`` is
+    the ``EchoNumber`` of the series - this is only present for
+    multi-echo data, where the series is from the second or subsequent
+    echos.
+
+    :arg series: Dict containing BIDS metadata about a DICOM series,
+                 as returned by :func:`scanDir`.
+
+    :returns:    String containing a CRC32 checksum for the series.
+    """
+
+    uid  = series.get('SeriesInstanceUID', None)
+    echo = series.get('EchoNumber',        None)
+
+    if uid is None:
+        return None
+
+    crc32 = str(binascii.crc32(uid.encode()))
+
+    if echo is not None and echo > 1:
+        crc32 = f'{crc32}.{echo}'
+
+    return crc32
 
 
 def loadSeries(series):
@@ -192,24 +277,39 @@ def loadSeries(series):
     if not enabled():
         raise RuntimeError('dcm2niix is not available or is too old')
 
-    dcmdir = series['DicomDir']
-    snum   = series['SeriesNumber']
-    desc   = series['SeriesDescription']
-    cmd    = 'dcm2niix -b n -f %s -z n -o . -n {} {}'.format(snum, dcmdir)
+    dcmdir  = series['DicomDir']
+    snum    = series['SeriesNumber']
+    desc    = series['SeriesDescription']
+    version = installedVersion()
+
+    # Newer versions of dcm2niix
+    # require a CRC to identify
+    # series
+    if compareVersions(version, CRC_DCM2NIIX_VERSION) >= 0:
+        ident = seriesCRC(series)
+
+    # Older versions require
+    # the series number
+    else:
+        ident = snum
+
+    cmd = f'{dcm2niix()} -b n -f %s -z n -o . -n "{ident}" "{dcmdir}"'
 
     with tempdir.tempdir() as td:
 
         with open(os.devnull, 'wb') as devnull:
-            sp.call(cmd.split(), stdout=devnull, stderr=devnull)
+            sp.call(shlex.split(cmd), stdout=devnull, stderr=devnull)
 
-        files  = glob.glob(op.join(td, '{}*.nii'.format(snum)))
+        files  = glob.glob(op.join(td, f'{snum}*.nii'))
         images = [nib.load(f, mmap=False) for f in files]
 
         # copy images so nibabel no longer
         # refs to the files (as they will
-        # be deleted), and use get_data()
-        # to force-load the image data.
-        images = [nib.Nifti1Image(i.get_data(), None, i.header)
+        # be deleted), and force-load the
+        # the image data into memory (to
+        # avoid any disk accesses due to
+        # e.g. memmap)
+        images = [nib.Nifti1Image(np.asanyarray(i.dataobj), None, i.header)
                   for i in images]
 
         return [DicomImage(i, series, dcmdir, name=desc) for i in images]

fsl/data/featanalysis.py

@@ -22,9 +22,12 @@ following functions are provided:
    isFirstLevelAnalysis
    loadDesign
    loadContrasts
+   loadFTests
+   loadFsf
    loadSettings
    getThresholds
    loadClusterResults
+   loadFEATDesignFile
 
 
 The following functions return the names of various files of interest:
@@ -38,20 +41,22 @@ The following functions return the names of various files of interest:
    getPEFile
    getCOPEFile
    getZStatFile
+   getZFStatFile
    getClusterMaskFile
+   getFClusterMaskFile
 """
 
 
-import                        collections
-import                        logging
-import os.path             as op
-import numpy               as np
+import                   collections
+import                   io
+import                   logging
+import os.path        as op
+import numpy          as np
 
-import fsl.utils.path      as fslpath
-import fsl.utils.transform as transform
-
-from . import image as fslimage
-from . import          featdesign
+import fsl.utils.path       as fslpath
+import fsl.transform.affine as affine
+from . import image         as fslimage
+from . import                  featdesign
 
 
 log = logging.getLogger(__name__)
@@ -166,70 +171,83 @@ def loadContrasts(featdir):
     :arg featdir: A FEAT directory.
     """
 
-    matrix       = None
-    numContrasts = 0
-    names        = {}
-    designcon    = op.join(featdir, 'design.con')
+    filename = op.join(featdir, 'design.con')
 
-    log.debug('Loading FEAT contrasts from {}'.format(designcon))
+    log.debug('Loading FEAT contrasts from %s', filename)
 
-    with open(designcon, 'rt') as f:
+    try:
+        designcon    = loadFEATDesignFile(filename)
+        contrasts    = np.genfromtxt(io.StringIO(designcon['Matrix']), ndmin=2)
+        numContrasts = int(designcon['NumContrasts'])
+        names        = []
 
-        while True:
-            line = f.readline().strip()
+        if numContrasts != contrasts.shape[0]:
+            raise RuntimeError(f'Matrix shape {contrasts.shape} does not '
+                               f'match number of contrasts {numContrasts}')
 
-            if line.startswith('/ContrastName'):
-                tkns       = line.split(None, 1)
-                num        = [c for c in tkns[0] if c.isdigit()]
-                num        = int(''.join(num))
+        contrasts = [list(row) for row in contrasts]
 
-                # The /ContrastName field may not
-                # actually have a name specified
-                if len(tkns) > 1:
-                    name       = tkns[1].strip()
-                    names[num] = name
+        for i in range(numContrasts):
+            cname = designcon.get(f'ContrastName{i + 1}', '')
+            if cname == '':
+                cname = f'{i + 1}'
+            names.append(cname)
 
-            elif line.startswith('/NumContrasts'):
-                numContrasts = int(line.split()[1])
+    except Exception as e:
+        log.debug('Error reading %s: %s', filename, e, exc_info=True)
+        raise RuntimeError(f'{filename} does not appear '
+                           'to be a valid design.con file') from e
 
-            elif line == '/Matrix':
-                break
+    return names, contrasts
 
-        matrix = np.loadtxt(f, ndmin=2)
 
-    if matrix       is None             or \
-       numContrasts != matrix.shape[0]:
-        raise RuntimeError('{} does not appear to be a '
-                           'valid design.con file'.format(designcon))
+def loadFTests(featdir):
+    """Loads F-tests from a FEAT directory. Returns a list of f-test vectors
+    (each of which is a list itself), where each vector contains a 1 or a 0
+    denoting the contrasts included in the F-test.
 
-    # Fill in any missing contrast names
-    if len(names) != numContrasts:
-        for i in range(numContrasts):
-            if i + 1 not in names:
-                names[i + 1] = str(i + 1)
+    :arg featdir: A FEAT directory.
+    """
 
-    names     = [names[c + 1] for c in range(numContrasts)]
-    contrasts = []
+    filename = op.join(featdir, 'design.fts')
 
-    for row in matrix:
-        contrasts.append(list(row))
+    if not op.exists(filename):
+        return []
 
-    return names, contrasts
+    log.debug('Loading FEAT F-tests from %s', filename)
 
+    try:
+        desfts = loadFEATDesignFile(filename)
+        ftests = np.genfromtxt(io.StringIO(desfts['Matrix']), ndmin=2)
+        ncols  = int(desfts['NumWaves'])
+        nrows  = int(desfts['NumContrasts'])
 
-def loadSettings(featdir):
-    """Loads the analysis settings from a FEAT directory.
+        if ftests.shape != (nrows, ncols):
+            raise RuntimeError(f'Matrix shape {ftests.shape} does not match '
+                               f'number of EVs/FTests ({ncols}, {nrows})')
 
-    Returns a dict containing the settings specified in the ``design.fsf``
-    file within the directory
+        ftests = [list(row) for row in ftests]
 
-    :arg featdir: A FEAT directory.
+    except Exception as e:
+        log.debug('Error reading %s: %s', filename, e, exc_info=True)
+        raise RuntimeError(f'{filename} does not appear '
+                           'to be a valid design.fts file') from e
+
+    return ftests
+
+
+def loadFsf(designfsf):
+    """Loads the analysis settings from a text file (.fsf) used to configure
+    FEAT.
+
+    Returns a dict containing the settings specified in the file
+
+    :arg designfsf: A .fsf file.
     """
 
     settings  = collections.OrderedDict()
-    designfsf = op.join(featdir, 'design.fsf')
 
-    log.debug('Loading FEAT settings from {}'.format(designfsf))
+    log.debug('Loading FEAT settings from %s', designfsf)
 
     with open(designfsf, 'rt') as f:
 
@@ -252,6 +270,20 @@ def loadSettings(featdir):
     return settings
 
 
+def loadSettings(featdir):
+    """Loads the analysis settings from a FEAT directory.
+
+    Returns a dict containing the settings specified in the ``design.fsf``
+    file within the directory
+
+    :arg featdir: A FEAT directory.
+    """
+
+    designfsf = op.join(featdir, 'design.fsf')
+
+    return loadFsf(designfsf)
+
+
 def loadDesign(featdir, settings):
     """Loads the design matrix from a FEAT directory.
 
@@ -297,19 +329,22 @@ def isFirstLevelAnalysis(settings):
     return settings['level'] == '1'
 
 
-def loadClusterResults(featdir, settings, contrast):
+def loadClusterResults(featdir, settings, contrast, ftest=False):
     """If cluster thresholding was used in the FEAT analysis, this function
     will load and return the cluster results for the specified (0-indexed)
-    contrast number.
+    contrast or f-test.
 
-    If there are no cluster results for the given contrast, ``None`` is
-    returned.
+    If there are no cluster results for the given contrast/f-test, ``None``
+    is returned.
 
     An error will be raised if the cluster file cannot be parsed.
 
     :arg featdir:  A FEAT directory.
     :arg settings: A FEAT settings dictionary.
-    :arg contrast: 0-indexed contrast number.
+    :arg contrast: 0-indexed contrast or f-test number.
+    :arg ftest:    If ``False`` (default), return cluster results for
+                   the contrast numbered ``contrast``. Otherwise, return
+                   cluster results for the f-test numbered ``contrast``.
 
     :returns:      A list of ``Cluster`` instances, each of which contains
                    information about one cluster. A ``Cluster`` object has the
@@ -330,11 +365,16 @@ def loadClusterResults(featdir, settings, contrast):
                                   gravity.
                      ``zcogz``    Z voxel coordinate of cluster centre of
                                   gravity.
-                     ``copemax``  Maximum COPE value in cluster.
-                     ``copemaxx`` X voxel coordinate of maximum COPE value.
+                     ``copemax``  Maximum COPE value in cluster (not
+                                  present for f-tests).
+                     ``copemaxx`` X voxel coordinate of maximum COPE value
+                                  (not present for f-tests).
                      ``copemaxy`` Y voxel coordinate of maximum COPE value.
+                                  (not present for f-tests).
                      ``copemaxz`` Z voxel coordinate of maximum COPE value.
+                                  (not present for f-tests).
                      ``copemean`` Mean COPE of all voxels in the cluster.
+                                  (not present for f-tests).
                      ============ =========================================
     """
 
@@ -344,8 +384,11 @@ def loadClusterResults(featdir, settings, contrast):
     # the ZMax/COG etc coordinates
     # are usually in voxel coordinates
     coordXform  = np.eye(4)
-    clusterFile = op.join(
-        featdir, 'cluster_zstat{}.txt'.format(contrast + 1))
+
+    if ftest: prefix = 'cluster_zfstat'
+    else:     prefix = 'cluster_zstat'
+
+    clusterFile = op.join(featdir, f'{prefix}{contrast + 1}.txt')
 
     if not op.exists(clusterFile):
 
@@ -354,22 +397,16 @@ def loadClusterResults(featdir, settings, contrast):
         # the cluster file will instead be called
         # 'cluster_zstatX_std.txt', so we'd better
         # check for that too.
-        clusterFile = op.join(
-            featdir, 'cluster_zstat{}_std.txt'.format(contrast + 1))
+        clusterFile = op.join(featdir, f'{prefix}{contrast + 1}_std.txt')
 
         if not op.exists(clusterFile):
             return None
 
-        # In higher levle analysis run in some standard
+        # In higher level analysis run in some standard
         # space, the cluster coordinates are in standard
         # space. We transform them to voxel coordinates.
         # later on.
-        coordXform = fslimage.Image(
-            getDataFile(featdir),
-            loadData=False).worldToVoxMat
-
-    log.debug('Loading cluster results for contrast {} from {}'.format(
-        contrast, clusterFile))
+        coordXform = fslimage.Image(getDataFile(featdir)).worldToVoxMat
 
     # The cluster.txt file is converted
     # into a list of Cluster objects,
@@ -387,10 +424,18 @@ def loadClusterResults(featdir, settings, contrast):
 
             # if cluster thresholding was not used,
             # the cluster table will not contain
-            # P valuse.
+            # P values.
             if not hasattr(self, 'p'):    self.p    = 1.0
             if not hasattr(self, 'logp'): self.logp = 0.0
 
+            # F-test cluster results will not have
+            # COPE-* results
+            if not hasattr(self, 'copemax'):  self.copemax  = np.nan
+            if not hasattr(self, 'copemaxx'): self.copemaxx = np.nan
+            if not hasattr(self, 'copemaxy'): self.copemaxy = np.nan
+            if not hasattr(self, 'copemaxz'): self.copemaxz = np.nan
+            if not hasattr(self, 'copemean'): self.copemean = np.nan
+
     # This dict provides a mapping between
     # Cluster object attribute names, and
     # the corresponding column name in the
@@ -422,10 +467,9 @@ def loadClusterResults(featdir, settings, contrast):
         'COPE-MAX Z (mm)'  : 'copemaxz',
         'COPE-MEAN'        : 'copemean'}
 
-    # An error will be raised if the
-    # cluster file does not exist (e.g.
-    # if the specified contrast index
-    # is invalid)
+    log.debug('Loading cluster results for contrast %s from %s',
+              contrast, clusterFile)
+
     with open(clusterFile, 'rt') as f:
 
         # Get every line in the file,
@@ -447,12 +491,11 @@ def loadClusterResults(featdir, settings, contrast):
     colNames     =  colNames.split('\t')
     clusterLines = [cl      .split('\t') for cl in clusterLines]
 
-    # Turn each cluster line into a
-    # Cluster instance. An error will
-    # be raised if the columm names
-    # are unrecognised (i.e. not in
-    # the colmap above), or if the
-    # file is poorly formed.
+    # Turn each cluster line into a Cluster
+    # instance. An error will be raised if the
+    # columm names are unrecognised (i.e. not
+    # in the colmap above), or if the file is
+    # poorly formed.
     clusters = [Cluster(**dict(zip(colNames, cl))) for cl in clusterLines]
 
     # Make sure all coordinates are in voxels -
@@ -467,17 +510,51 @@ def loadClusterResults(featdir, settings, contrast):
         zcog    = [c.zcogx,    c.zcogy,    c.zcogz]
         copemax = [c.copemaxx, c.copemaxy, c.copemaxz]
 
-        zmax    = transform.transform([zmax],    coordXform)[0].round()
-        zcog    = transform.transform([zcog],    coordXform)[0].round()
-        copemax = transform.transform([copemax], coordXform)[0].round()
+        zmax    = affine.transform([zmax],    coordXform)[0]
+        zcog    = affine.transform([zcog],    coordXform)[0]
+        copemax = affine.transform([copemax], coordXform)[0]
 
-        c.zmaxx,   c.zmaxy,    c.zmaxz    = zmax
-        c.zcogx,   c.zcogy,    c.zcogz    = zcog
-        c.copemax, c.copemaxy, c.copemaxz = copemax
+        c.zmaxx,    c.zmaxy,    c.zmaxz    = zmax
+        c.zcogx,    c.zcogy,    c.zcogz    = zcog
+        c.copemaxx, c.copemaxy, c.copemaxz = copemax
 
     return clusters
 
 
+def loadFEATDesignFile(filename):
+    """Load a FEAT design file, e.g. ``design.mat``, ``design.con``, ``design.fts``.
+
+    These files contain key-value pairs, and are formatted according to an
+    undocumented structure where each key is of the form "/KeyName", and is
+    followed immediately by a whitespace character, and then the value.
+
+    :arg filename: File to load
+    :returns:      A dictionary of key-value pairs. The values are all left
+                   as strings.
+    """
+
+    fields = {}
+
+    with open(filename, 'rt') as f:
+        content = f.read()
+
+    content = content.split('/')
+    for line in content:
+        line = line.strip()
+        if line == '':
+            continue
+
+        tokens = line.split(maxsplit=1)
+        if len(tokens) == 1:
+            name, value = tokens[0], ''
+        else:
+            name, value = tokens
+
+        fields[name] = value
+
+    return fields
+
+
 def getDataFile(featdir):
     """Returns the name of the file in the FEAT directory which contains
     the model input data (typically called ``filtered_func_data.nii.gz``).
@@ -521,7 +598,7 @@ def getPEFile(featdir, ev):
     :arg featdir: A FEAT directory.
     :arg ev:      The EV number (0-indexed).
     """
-    pefile = op.join(featdir, 'stats', 'pe{}'.format(ev + 1))
+    pefile = op.join(featdir, 'stats', f'pe{ev + 1}')
     return fslimage.addExt(pefile, mustExist=True)
 
 
@@ -533,7 +610,7 @@ def getCOPEFile(featdir, contrast):
     :arg featdir:  A FEAT directory.
     :arg contrast: The contrast number (0-indexed).
     """
-    copefile = op.join(featdir, 'stats', 'cope{}'.format(contrast + 1))
+    copefile = op.join(featdir, 'stats', f'cope{contrast + 1}')
     return fslimage.addExt(copefile, mustExist=True)
 
 
@@ -545,10 +622,22 @@ def getZStatFile(featdir, contrast):
     :arg featdir:  A FEAT directory.
     :arg contrast: The contrast number (0-indexed).
     """
-    zfile = op.join(featdir, 'stats', 'zstat{}'.format(contrast + 1))
+    zfile = op.join(featdir, 'stats', f'zstat{contrast + 1}')
     return fslimage.addExt(zfile, mustExist=True)
 
 
+def getZFStatFile(featdir, ftest):
+    """Returns the path of the Z-statistic file for the specified F-test.
+
+    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
+
+    :arg featdir: A FEAT directory.
+    :arg ftest:   The F-test number (0-indexed).
+    """
+    zffile = op.join(featdir, 'stats', f'zfstat{ftest + 1}')
+    return fslimage.addExt(zffile, mustExist=True)
+
+
 def getClusterMaskFile(featdir, contrast):
     """Returns the path of the cluster mask file for the specified contrast.
 
@@ -557,5 +646,17 @@ def getClusterMaskFile(featdir, contrast):
     :arg featdir:  A FEAT directory.
     :arg contrast: The contrast number (0-indexed).
     """
-    mfile = op.join(featdir, 'cluster_mask_zstat{}'.format(contrast + 1))
+    mfile = op.join(featdir, f'cluster_mask_zstat{contrast + 1}')
+    return fslimage.addExt(mfile, mustExist=True)
+
+
+def getFClusterMaskFile(featdir, ftest):
+    """Returns the path of the cluster mask file for the specified f-test.
+
+    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
+
+    :arg featdir:  A FEAT directory.
+    :arg contrast: The f-test number (0-indexed).
+    """
+    mfile = op.join(featdir, f'cluster_mask_zfstat{ftest + 1}')
     return fslimage.addExt(mfile, mustExist=True)