doc/fsl.wrappers.tbss.rst

+``fsl.wrappers.tbss``
+=====================
+
+.. automodule:: fsl.wrappers.tbss
+    :members:
+    :undoc-members:
+    :show-inheritance:

doc/index.rst

@@ -8,24 +8,53 @@ within `FSL <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki>`_ and by
 |fsleyes_apidoc|_.
 
 
-The top-level Python package for ``fslpy`` is called ``fsl``.  It is broadly
-split into the following sub-packages:
+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.
 
 
-.. autosummary::
-
-   fsl.data
-   fsl.utils
-   fsl.scripts
-   fsl.transform
-   fsl.version
-   fsl.wrappers
-
 
 .. 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
@@ -8,3 +10,4 @@ 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.transform
-   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__)  # noqa

fsl/data/atlases.py

@@ -377,7 +377,7 @@ class AtlasLabel(object):
         )
 
 
-class AtlasDescription(object):
+class AtlasDescription:
     """An ``AtlasDescription`` instance parses and stores the information
     stored in the FSL XML file that describes a single FSL atlas.  An XML
     atlas specification file is assumed to have a structure that looks like
@@ -560,7 +560,7 @@ class AtlasDescription(object):
             imagefile        = op.normpath(atlasDir + imagefile)
             summaryimagefile = op.normpath(atlasDir + summaryimagefile)
 
-            i = fslimage.Image(imagefile, loadData=False, calcRange=False)
+            i = fslimage.Image(imagefile)
 
             self.images       .append(imagefile)
             self.summaryImages.append(summaryimagefile)
@@ -880,10 +880,17 @@ class LabelAtlas(Atlas):
                        of each present value. The proportions are returned as
                        values between 0 and 100.
 
+        .. note:: Calling this method will cause the atlas image data to be
+                  loaded into memory.
+
         .. note:: Use the :meth:`find` method to retrieve the ``AtlasLabel``
                   associated with each returned value.
         """
 
+        # Mask-based indexing requires the image
+        # data to be loaded into memory
+        self.data
+
         # Extract the values that are in
         # the mask, and their corresponding
         # mask weights

fsl/data/bitmap.py

@@ -9,20 +9,21 @@ files. Pillow is required to use the ``Bitmap`` class.
 """
 
 
-import os.path      as op
-import                 logging
-import                 six
+import os.path as op
+import            pathlib
+import            logging
 
-import numpy        as np
+import numpy as np
 
-from . import image as fslimage
+import fsl.data.image as fslimage
 
 
 log = logging.getLogger(__name__)
 
 
 BITMAP_EXTENSIONS = ['.bmp', '.png',  '.jpg', '.jpeg',
-                     '.tif', '.tiff', '.gif', '.rgba']
+                     '.tif', '.tiff', '.gif', '.rgba',
+                     '.jp2', '.jpg2', '.jp2k']
 """File extensions we understand. """
 
 
@@ -34,7 +35,10 @@ BITMAP_DESCRIPTIONS = [
     'TIFF',
     'TIFF',
     'Graphics Interchange Format',
-    'Raw RGBA']
+    'Raw RGBA',
+    'JPEG 2000',
+    'JPEG 2000',
+    'JPEG 2000']
 """A description for each :attr:`BITMAP_EXTENSION`. """
 
 
@@ -51,17 +55,19 @@ class Bitmap(object):
                   data.
         """
 
-        if isinstance(bmp, six.string_types):
+        if isinstance(bmp, (pathlib.Path, str)):
 
             try:
-                # Allow big images
-                import PIL.Image as Image
-                Image.MAX_IMAGE_PIXELS = 1e9
+                # 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 = bmp
+            src = str(bmp)
             img = Image.open(src)
 
             # If this is a palette/LUT
@@ -173,7 +179,7 @@ class Bitmap(object):
             for ci, ch in enumerate(dtype.names):
                 data[ch] = self.data[..., ci]
 
-        data = np.array(data, order='F', copy=False)
+        data = np.asarray(data, order='F')
 
         return fslimage.Image(data,
                               name=self.name,

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

@@ -33,15 +33,17 @@ 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__)
@@ -60,6 +62,25 @@ 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):
     """The ``DicomImage`` is a volumetric :class:`.Image` with some associated
     DICOM metadata.
@@ -105,7 +126,7 @@ def installedVersion():
       - Day
     """
 
-    cmd            = 'dcm2niix -h'
+    cmd            = f'{dcm2niix()} -h'
     versionPattern = re.compile(r'v'
                                 r'(?P<major>[0-9]+)\.'
                                 r'(?P<minor>[0-9]+)\.'
@@ -130,7 +151,7 @@ def installedVersion():
                         int(match.group('day')))
 
     except Exception as e:
-        log.debug('Error parsing dcm2niix version string: {}'.format(e))
+        log.debug(f'Error parsing dcm2niix version string: {e}')
     return None
 
 
@@ -177,7 +198,7 @@ def scanDir(dcmdir):
         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)
+    cmd    = f'{dcm2niix()} -b o -ba n -f %s -o . "{dcmdir}"'
     series = []
 
     with tempdir.tempdir() as td:
@@ -194,6 +215,10 @@ def scanDir(dcmdir):
             with open(fn, 'rt') as 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)
 
     # sort by series number
@@ -233,7 +258,7 @@ def seriesCRC(series):
     crc32 = str(binascii.crc32(uid.encode()))
 
     if echo is not None and echo > 1:
-        crc32 = '{}.{}'.format(crc32, echo)
+        crc32 = f'{crc32}.{echo}'
 
     return crc32
 
@@ -268,14 +293,14 @@ def loadSeries(series):
     else:
         ident = snum
 
-    cmd = 'dcm2niix -b n -f %s -z n -o . -n "{}" "{}"'.format(ident, dcmdir)
+    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(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

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,11 +41,14 @@ The following functions return the names of various files of interest:
    getPEFile
    getCOPEFile
    getZStatFile
+   getZFStatFile
    getClusterMaskFile
+   getFClusterMaskFile
 """
 
 
 import                   collections
+import                   io
 import                   logging
 import os.path        as op
 import numpy          as np
@@ -165,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'])
 
+        if ftests.shape != (nrows, ncols):
+            raise RuntimeError(f'Matrix shape {ftests.shape} does not match '
+                               f'number of EVs/FTests ({ncols}, {nrows})')
 
-def loadSettings(featdir):
-    """Loads the analysis settings from a FEAT directory.
+        ftests = [list(row) for row in ftests]
 
-    Returns a dict containing the settings specified in the ``design.fsf``
-    file within the 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
 
-    :arg featdir: A FEAT directory.
+    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:
 
@@ -251,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.
 
@@ -296,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
@@ -329,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).
                      ============ =========================================
     """
 
@@ -343,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):
 
@@ -353,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,
@@ -386,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
@@ -421,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,
@@ -446,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 -
@@ -466,17 +510,51 @@ def loadClusterResults(featdir, settings, contrast):
         zcog    = [c.zcogx,    c.zcogy,    c.zcogz]
         copemax = [c.copemaxx, c.copemaxy, c.copemaxz]
 
-        zmax    = affine.transform([zmax],    coordXform)[0].round()
-        zcog    = affine.transform([zcog],    coordXform)[0].round()
-        copemax = affine.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``).
@@ -520,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)
 
 
@@ -532,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)
 
 
@@ -544,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.
 
@@ -556,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)

fsl/data/featdesign.py

@@ -160,7 +160,7 @@ class FEATFSFDesign(object):
         # Print a warning if we're
         # using an old version of FEAT
         if version < 6:
-            log.warning('Unsupported FEAT version: {}'.format(version))
+            log.warning('Unsupported FEAT version: %s', version)
 
         # We need to parse the EVS a bit
         # differently depending on whether
@@ -210,8 +210,7 @@ class FEATFSFDesign(object):
                 continue
 
             if (not self.__loadVoxEVs) or (ev.filename is None):
-                log.warning('Voxel EV image missing '
-                            'for ev {}'.format(ev.index))
+                log.warning('Voxel EV image missing for ev %s', ev.index)
                 continue
 
             design[:, ev.index] = ev.getData(x, y, z)
@@ -250,8 +249,7 @@ class VoxelwiseEVMixin(object):
         if op.exists(filename):
             self.__filename = filename
         else:
-            log.warning('Voxelwise EV file does not '
-                        'exist: {}'.format(filename))
+            log.warning('Voxelwise EV file does not exist: %s', filename)
             self.__filename = None
 
         self.__image = None
@@ -502,11 +500,11 @@ def getFirstLevelEVs(featDir, settings, designMat):
     #   - voxelwise EVs
     for origIdx in range(origEVs):
 
-        title    = settings[        'evtitle{}'  .format(origIdx + 1)]
-        shape    = int(settings[    'shape{}'    .format(origIdx + 1)])
-        convolve = int(settings[    'convolve{}' .format(origIdx + 1)])
-        deriv    = int(settings[    'deriv_yn{}' .format(origIdx + 1)])
-        basis    = int(settings.get('basisfnum{}'.format(origIdx + 1), -1))
+        title    = settings[        f'evtitle{origIdx + 1}']
+        shape    = int(settings[    f'shape{origIdx + 1}'])
+        convolve = int(settings[    f'convolve{origIdx + 1}'])
+        deriv    = int(settings[    f'deriv_yn{origIdx + 1}'])
+        basis    = int(settings.get(f'basisfnum{origIdx + 1}', -1))
 
         # Normal EV. This is just a column
         # in the design matrix, defined by
@@ -525,8 +523,7 @@ def getFirstLevelEVs(featDir, settings, designMat):
 
             # The addExt function will raise an
             # error if the file does not exist.
-            filename = op.join(
-                featDir, 'designVoxelwiseEV{}'.format(origIdx + 1))
+            filename = op.join(featDir, f'designVoxelwiseEV{origIdx + 1}')
             filename = fslimage.addExt(filename, True)
 
             evs.append(VoxelwiseEV(len(evs), origIdx, title, filename))
@@ -607,7 +604,7 @@ def getFirstLevelEVs(featDir, settings, designMat):
         startIdx = len(evs) + 1
         if voxConfLocs != list(range(startIdx, startIdx + len(voxConfFiles))):
             raise FSFError('Unsupported voxelwise confound ordering '
-                           '({} -> {})'.format(startIdx, voxConfLocs))
+                           f'({startIdx} -> {voxConfLocs})')
 
         # Create the voxelwise confound EVs.
         # We make a name for the EV from the
@@ -680,7 +677,7 @@ def getHigherLevelEVs(featDir, settings, designMat):
     for origIdx in range(origEVs):
 
         # All we need is the title
-        title = settings['evtitle{}'.format(origIdx + 1)]
+        title = settings[f'evtitle{origIdx + 1}']
         evs.append(NormalEV(len(evs), origIdx, title))
 
     # Only the input file is specified for
@@ -689,7 +686,7 @@ def getHigherLevelEVs(featDir, settings, designMat):
     # name.
     for origIdx in range(voxEVs):
 
-        filename = settings['evs_vox_{}'.format(origIdx + 1)]
+        filename = settings[f'evs_vox_{origIdx + 1}']
         title    = op.basename(fslimage.removeExt(filename))
         evs.append(VoxelwiseEV(len(evs), origIdx, title, filename))
 
@@ -705,12 +702,12 @@ def loadDesignMat(designmat):
     :arg designmat: Path to the ``design.mat`` file.
     """
 
-    log.debug('Loading FEAT design matrix from {}'.format(designmat))
+    log.debug('Loading FEAT design matrix from %s', designmat)
 
     matrix = np.loadtxt(designmat, comments='/', ndmin=2)
 
     if matrix is None or matrix.size == 0 or len(matrix.shape) != 2:
-        raise FSFError('{} does not appear to be a '
-                       'valid design.mat file'.format(designmat))
+        raise FSFError(f'{designmat} does not appear '
+                       'to be a valid design.mat file')
 
     return matrix

fsl/data/featimage.py

@@ -63,8 +63,8 @@ class FEATImage(fslimage.Image):
             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))
+            raise ValueError(f'{path} does not appear to be '
+                             'data from a FEAT analysis')
 
         featDir     = op.dirname(path)
         settings    = featanalysis.loadSettings( featDir)
@@ -72,9 +72,11 @@ class FEATImage(fslimage.Image):
         if featanalysis.hasStats(featDir):
             design      = featanalysis.loadDesign(   featDir, settings)
             names, cons = featanalysis.loadContrasts(featDir)
+            ftests      = featanalysis.loadFTests(   featDir)
         else:
             design      = None
             names, cons = [], []
+            ftests      = []
 
         fslimage.Image.__init__(self, path, **kwargs)
 
@@ -83,26 +85,31 @@ class FEATImage(fslimage.Image):
         self.__design        = design
         self.__contrastNames = names
         self.__contrasts     = cons
+        self.__ftests        = ftests
         self.__settings      = settings
 
         self.__residuals     =  None
         self.__pes           = [None] * self.numEVs()
         self.__copes         = [None] * self.numContrasts()
         self.__zstats        = [None] * self.numContrasts()
+        self.__zfstats       = [None] * self.numFTests()
         self.__clustMasks    = [None] * self.numContrasts()
+        self.__fclustMasks   = [None] * self.numFTests()
 
         if 'name' not in kwargs:
-            self.name = '{}: {}'.format(self.__analysisName, self.name)
+            self.name = f'{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
+        self.__design      = None
+        self.__residuals   = None
+        self.__pes         = None
+        self.__copes       = None
+        self.__zstats      = None
+        self.__zfstats     = None
+        self.__clustMasks  = None
+        self.__fclustMasks = None
 
 
     def getFEATDir(self):
@@ -191,6 +198,11 @@ class FEATImage(fslimage.Image):
         return len(self.__contrasts)
 
 
+    def numFTests(self):
+        """Returns the number of f-tests in the analysis."""
+        return len(self.__ftests)
+
+
     def contrastNames(self):
         """Returns a list containing the name of each contrast in the analysis.
         """
@@ -206,6 +218,15 @@ class FEATImage(fslimage.Image):
         return [list(c) for c in self.__contrasts]
 
 
+    def ftests(self):
+        """Returns a list containing the analysis f-test vectors.
+
+        See :func:`.featanalysis.loadFTests`
+
+        """
+        return [list(f) for f in self.__ftests]
+
+
     def thresholds(self):
         """Returns the statistical thresholds used in the analysis.
 
@@ -214,14 +235,16 @@ class FEATImage(fslimage.Image):
         return featanalysis.getThresholds(self.__settings)
 
 
-    def clusterResults(self, contrast):
-        """Returns the clusters found in the analysis.
+    def clusterResults(self, contrast, ftest=False):
+        """Returns the clusters found in the analysis for the specified
+        contrast or f-test.
 
         See :func:.featanalysis.loadClusterResults`
         """
         return featanalysis.loadClusterResults(self.__featDir,
                                                self.__settings,
-                                               contrast)
+                                               contrast,
+                                               ftest)
 
 
     def getPE(self, ev):
@@ -229,12 +252,10 @@ class FEATImage(fslimage.Image):
 
         if self.__pes[ev] is None:
             pefile = featanalysis.getPEFile(self.__featDir, ev)
+            evname = self.evNames()[ev]
             self.__pes[ev] = fslimage.Image(
                 pefile,
-                name='{}: PE{} ({})'.format(
-                    self.__analysisName,
-                    ev + 1,
-                    self.evNames()[ev]))
+                name=f'{self.__analysisName}: PE{ev + 1} ({evname})')
 
         return self.__pes[ev]
 
@@ -246,7 +267,7 @@ class FEATImage(fslimage.Image):
             resfile = featanalysis.getResidualFile(self.__featDir)
             self.__residuals = fslimage.Image(
                 resfile,
-                name='{}: residuals'.format(self.__analysisName))
+                name=f'{self.__analysisName}: residuals')
 
         return self.__residuals
 
@@ -256,12 +277,10 @@ class FEATImage(fslimage.Image):
 
         if self.__copes[con] is None:
             copefile = featanalysis.getCOPEFile(self.__featDir, con)
+            conname  = self.contrastNames()[con]
             self.__copes[con] = fslimage.Image(
                 copefile,
-                name='{}: COPE{} ({})'.format(
-                    self.__analysisName,
-                    con + 1,
-                    self.contrastNames()[con]))
+                name=f'{self.__analysisName}: COPE{con + 1} ({conname})')
         return self.__copes[con]
 
 
@@ -270,35 +289,54 @@ class FEATImage(fslimage.Image):
         """
 
         if self.__zstats[con] is None:
-            zfile = featanalysis.getZStatFile(self.__featDir, con)
-
+            zfile   = featanalysis.getZStatFile(self.__featDir, con)
+            conname = self.contrastNames()[con]
             self.__zstats[con] = fslimage.Image(
                 zfile,
-                name='{}: zstat{} ({})'.format(
-                    self.__analysisName,
-                    con + 1,
-                    self.contrastNames()[con]))
-
+                name=f'{self.__analysisName}: zstat{con + 1} ({conname})')
         return self.__zstats[con]
 
 
+    def getZFStats(self, ftest):
+        """Returns the Z statistic image for the given f-test (0-indexed). """
+
+        if self.__zfstats[ftest] is None:
+            zfile = featanalysis.getZFStatFile(self.__featDir, ftest)
+            self.__zfstats[ftest] = fslimage.Image(
+                zfile,
+                name=f'{self.__analysisName}: zfstat{ftest + 1}')
+        return self.__zfstats[ftest]
+
+
     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)
-
+            mfile   = featanalysis.getClusterMaskFile(self.__featDir, con)
+            conname = self.contrastNames()[con]
             self.__clustMasks[con] = fslimage.Image(
                 mfile,
-                name='{}: cluster mask for zstat{} ({})'.format(
-                    self.__analysisName,
-                    con + 1,
-                    self.contrastNames()[con]))
+                name=f'{self.__analysisName}: cluster mask '
+                     f'for zstat{con + 1} ({conname})')
 
         return self.__clustMasks[con]
 
 
+    def getFClusterMask(self, ftest):
+        """Returns the cluster mask image for the given f-test (0-indexed).
+        """
+
+        if self.__fclustMasks[ftest] is None:
+            mfile = featanalysis.getFClusterMaskFile(self.__featDir, ftest)
+            self.__fclustMasks[ftest] = fslimage.Image(
+                mfile,
+                name=f'{self.__analysisName}: cluster mask '
+                     f'for zfstat{ftest + 1}')
+
+        return self.__fclustMasks[ftest]
+
+
     def fit(self, contrast, xyz):
         """Calculates the model fit for the given contrast vector
         at the given voxel. See the :func:`modelFit` function.

fsl/data/fixlabels.py

@@ -16,18 +16,21 @@
 """
 
 
-import os.path as op
+import itertools as it
+import              math
+import os.path   as op
 
 
 def loadLabelFile(filename,
                   includeLabel=None,
                   excludeLabel=None,
-                  returnIndices=False):
-    """Loads component labels from the specified file. The file is assuemd
+                  returnIndices=False,
+                  missingLabel='Unknown',
+                  returnProbabilities=False):
+    """Loads component labels from the specified file. The file is assumed
     to be of the format generated by FIX, Melview or ICA-AROMA; such a file
     should have a structure resembling the following::
 
-
         filtered_func_data.ica
         1, Signal, False
         2, Unclassified Noise, True
@@ -39,7 +42,6 @@ def loadLabelFile(filename,
         8, Signal, False
         [2, 5, 6, 7]
 
-
     .. note:: This function will also parse files which only contain a
               component list, e.g.::
 
@@ -66,31 +68,46 @@ def loadLabelFile(filename,
       - One or more labels for the component (multiple labels must be
         comma-separated).
 
-      - ``'True'`` if the component has been classified as *bad*,
-        ``'False'`` otherwise. This field is optional - if the last
-        comma-separated token on a line is not equal (case-insensitive)
-        to ``True`` or ``False``, it is interpreted as a component label.
+      - ``'True'`` if the component has been classified as *bad*, ``'False'``
+        otherwise. This field is optional - if the last non-numeric
+        comma-separated token on a line is not equal to ``True`` or ``False``
+        (case-insensitive) , it is interpreted as a component label.
+
+      - A value between 0 and 1, which gives the probability of the component
+        being signal, as generated by an automatic classifier (e.g. FIX). This
+        field is optional - it is output by some versions of FIX.
 
     The last line of the file contains the index (starting from 1) of all
     *bad* components, i.e. those components which are not classified as
     signal or unknown.
 
-    :arg filename:      Name of the label file to load.
+    :arg filename:            Name of the label file to load.
+
+    :arg includeLabel:        If the file contains a single line containing a
+                              list component indices, this label will be used
+                              for the components in the list. Defaults to
+                              ``'Unclassified noise'`` for FIX-like files, and
+                              ``'Movement'`` for ICA-AROMA-like files.
 
-    :arg includeLabel:  If the file contains a single line containing a list
-                        component indices, this label will be used for the
-                        components in the list. Defaults to 'Unclassified
-                        noise' for FIX-like files, and 'Movement' for
-                        ICA-AROMA-like files.
+    :arg excludeLabel:        If the file contains a single line containing
+                              component indices, this label will be used for
+                              the  components that are not in the list.
+                              Defaults to ``'Signal'`` for FIX-like files, and
+                              ``'Unknown'`` for ICA-AROMA-like files.
 
-    :arg excludeLabel:  If the file contains a single line containing component
-                        indices, this label will be used for the components
-                        that are not in the list.  Defaults to 'Signal' for
-                        FIX-like files, and 'Unknown' for ICA-AROMA-like files.
+    :arg returnIndices:       Defaults to ``False``. If ``True``, a list
+                              containing the noisy component numbers that were
+                              listed in the file is returned.
 
-    :arg returnIndices: Defaults to ``False``. If ``True``, a list containing
-                        the noisy component numbers that were listed in the
-                        file is returned.
+    :arg missingLabel:        Label to use for any components which are not
+                              present (only used for label files, not for noise
+                              component files).
+
+    :arg returnProbabilities: Defaults to ``False``. If ``True``, a list
+                              containing the component classification
+                              probabilities is returned. If the file does not
+                              contain probabilities, every value in this list
+                              will be nan.
 
     :returns: A tuple containing:
 
@@ -102,72 +119,55 @@ def loadLabelFile(filename,
 
                - If ``returnIndices is True``, a list of the noisy component
                  indices (starting from 1) that were specified in the file.
+
+               - If ``returnProbabilities is True``, a list of the component
+                 classification probabilities that were specified in the
+                 file (all nan if they are not in the file).
+
+    .. note:: Some label files generated by old versions of FIX/Melview do
+              not contain a line for every component (unknown/unlabelled
+              components may not be listed). For these files, and also for
+              files which only contain a component list, there is no way of
+              knowing how many components were in the data, so the returned
+              list may contain fewer entries than there are components.
     """
 
-    signalLabels = None
-    filename     = op.abspath(filename)
+    filename      = op.abspath(filename)
+    probabilities = None
+    signalLabels  = None
 
     with open(filename, 'rt') as f:
         lines = f.readlines()
 
     if len(lines) < 1:
-        raise InvalidLabelFileError('Invalid FIX classification '
+        raise InvalidLabelFileError(f'{filename}: Invalid FIX classification '
                                     'file - not enough lines')
 
     lines = [l.strip() for l in lines]
     lines = [l for l in lines if l != '']
 
-    # If the file contains a single
-    # line, we assume that it is just
-    # a comma-separated list of noise
-    # components.
-    if len(lines) == 1:
-
-        line = lines[0]
-
-        # if the list is contained in
-        # square brackets, we assume
-        # that it is a FIX output file,
-        # where included components have
-        # been classified as noise, and
-        # excluded components as signal.
-        #
-        # Otherwise we assume that it
-        # is an AROMA file, where
-        # included components have
-        # been classified as being due
-        # to motion, and excluded
-        # components unclassified.
-        if includeLabel is None:
-            if line[0] == '[': includeLabel = 'Unclassified noise'
-            else:              includeLabel = 'Movement'
-
-        if excludeLabel is None:
-            if line[0] == '[': excludeLabel = 'Signal'
-            else:              excludeLabel = 'Unknown'
-        else:
-            signalLabels = [excludeLabel]
-
-        # Remove any leading/trailing
-        # whitespace or brackets.
-        line = lines[0].strip(' []')
-
-        melDir     = None
-        noisyComps = [int(i) for i in line.split(',')]
-        allLabels  = []
-
-        for i in range(max(noisyComps)):
-            if (i + 1) in noisyComps: allLabels.append([includeLabel])
-            else:                     allLabels.append([excludeLabel])
-
-    # Otherwise, we assume that
-    # it is a full label file.
-    else:
+    # If the file contains one or two lines, we
+    # assume that it is just a comma-separated list
+    # of noise components (possibly preceeded by
+    # the MELODIC directory path)
+    if len(lines) <= 2:
+        melDir, noisyComps, allLabels, signalLabels = \
+            _parseSingleLineLabelFile(lines, includeLabel, excludeLabel)
+        probabilities = [math.nan] * len(allLabels)
 
-        melDir     = lines[0]
-        noisyComps = lines[-1].strip(' []').split(',')
-        noisyComps = [c      for c in noisyComps if c != '']
-        noisyComps = [int(c) for c in noisyComps]
+    # Otherwise, we assume that it is a full label file.
+    else:
+        melDir, noisyComps, allLabels, probabilities = \
+            _parseFullLabelFile(filename, lines, missingLabel)
+
+    # There's no way to validate
+    # the melodic directory path,
+    # but let's try anyway.
+    if melDir is not None:
+        if len(melDir.split(',')) >= 3:
+            raise InvalidLabelFileError(
+                f'{filename}: First line does not look like '
+                f'a MELODIC directory path: {melDir}')
 
         # The melodic directory path should
         # either be an absolute path, or
@@ -176,38 +176,6 @@ def loadLabelFile(filename,
         if not op.isabs(melDir):
             melDir = op.join(op.dirname(filename), melDir)
 
-        # Parse the labels for every component
-        allLabels = []
-        for i, compLine in enumerate(lines[1:-1]):
-
-            tokens = compLine.split(',')
-            tokens = [t.strip() for t in tokens]
-
-            if len(tokens) < 3:
-                raise InvalidLabelFileError(
-                    'Invalid FIX classification file - '
-                    'line {}: {}'.format(i + 1, compLine))
-
-            try:
-                compIdx = int(tokens[0])
-
-            except ValueError:
-                raise InvalidLabelFileError(
-                    'Invalid FIX classification file - '
-                    'line {}: {}'.format(i + 1, compLine))
-
-            if tokens[-1].lower() in ('true', 'false'):
-                compLabels = tokens[1:-1]
-            else:
-                compLabels = tokens[1:]
-
-            if compIdx != i + 1:
-                raise InvalidLabelFileError(
-                    'Invalid FIX classification file - wrong component '
-                    'number at line {}: {}'.format(i + 1, compLine))
-
-            allLabels.append(compLabels)
-
     # Validate the labels against
     # the noisy list - all components
     # in the noisy list should not
@@ -218,8 +186,8 @@ def loadLabelFile(filename,
         noise = isNoisyComponent(labels, signalLabels)
 
         if noise and (comp not in noisyComps):
-            raise InvalidLabelFileError('Noisy component {} has invalid '
-                                        'labels: {}'.format(comp, labels))
+            raise InvalidLabelFileError(f'{filename}: Noisy component {comp} '
+                                        f'has invalid labels: {labels}')
 
     for comp in noisyComps:
 
@@ -228,44 +196,187 @@ def loadLabelFile(filename,
         noise  = isNoisyComponent(labels, signalLabels)
 
         if not noise:
-            raise InvalidLabelFileError('Noisy component {} is missing '
-                                        'a noise label'.format(comp))
+            raise InvalidLabelFileError(f'{filename}: Noisy component {comp} '
+                                        'is missing a noise label')
+
+    retval = [melDir, allLabels]
 
-    if returnIndices: return melDir, allLabels, noisyComps
-    else:             return melDir, allLabels
+    if returnIndices:       retval.append(noisyComps)
+    if returnProbabilities: retval.append(probabilities)
+
+    return tuple(retval)
+
+
+def _parseSingleLineLabelFile(lines, includeLabel, excludeLabel):
+    """Called by :func:`loadLabelFile`. Parses the contents of an
+    ICA-AROMA-style label file which just contains a list of noise
+    components (and possibly the MELODIC directory path), e.g.::
+
+        filtered_func_data.ica
+        [2, 5, 6, 7]
+    """
+    signalLabels = None
+    noisyComps   = lines[-1]
+
+    if len(lines) == 2: melDir = lines[0]
+    else:               melDir = None
+
+    # if the list is contained in
+    # square brackets, we assume
+    # that it is a FIX output file,
+    # where included components have
+    # been classified as noise, and
+    # excluded components as signal.
+    #
+    # Otherwise we assume that it
+    # is an AROMA file, where
+    # included components have
+    # been classified as being due
+    # to motion, and excluded
+    # components unclassified.
+    if includeLabel is None:
+        if noisyComps[0] == '[': includeLabel = 'Unclassified noise'
+        else:                    includeLabel = 'Movement'
+
+    if excludeLabel is None:
+        if noisyComps[0] == '[': excludeLabel = 'Signal'
+        else:                    excludeLabel = 'Unknown'
+    else:
+        signalLabels = [excludeLabel]
+
+    # Remove any leading/trailing
+    # whitespace or brackets.
+    noisyComps = noisyComps.strip(' []')
+    noisyComps = [int(i) for i in noisyComps.split(',')]
+    allLabels  = []
+
+    for i in range(max(noisyComps)):
+        if (i + 1) in noisyComps: allLabels.append([includeLabel])
+        else:                     allLabels.append([excludeLabel])
+
+    return melDir, noisyComps, allLabels, signalLabels
+
+
+def _parseFullLabelFile(filename, lines, missingLabel):
+    """Called by :func:`loadLabelFile`. Parses the contents of a
+    FIX/Melview-style label file which contains labels for each component,
+    e.g.:
+
+        filtered_func_data.ica
+        1, Signal, False
+        2, Unclassified Noise, True
+        3, Unknown, False
+        4, Signal, False
+        5, Unclassified Noise, True
+        6, Unclassified Noise, True
+        7, Unclassified Noise, True
+        8, Signal, False
+        [2, 5, 6, 7]
+    """
+    melDir     = lines[0]
+    noisyComps = lines[-1].strip(' []').split(',')
+    noisyComps = [c      for c in noisyComps if c != '']
+    noisyComps = [int(c) for c in noisyComps]
+
+    # Parse the labels for every component.
+    # Initially store as a {comp : ([labels], probability)} dict.
+    allLabels = {}
+    for i, compLine in enumerate(lines[1:-1]):
+
+        tokens = compLine.split(',')
+        tokens = [t.strip() for t in tokens]
+
+        if len(tokens) < 3:
+            raise InvalidLabelFileError(
+                f'{filename}: Invalid FIX classification '
+                f'file - line: {i + 1}: {compLine}')
+
+        try:
+            compIdx = int(tokens[0])
+            if compIdx in allLabels:
+                raise ValueError()
+
+        except ValueError:
+            raise InvalidLabelFileError(
+                f'{filename}: Invalid FIX classification '
+                f'file - line {i + 1}: {compLine}')
+
+        tokens      = tokens[1:]
+        probability = math.nan
+
+        # last token could be classification probability
+        if _isfloat(tokens[-1]):
+            probability = float(tokens[-1])
+            tokens      = tokens[:-1]
+
+        # true/false is ignored as it is superfluous
+        if tokens[-1].lower() in ('true', 'false'):
+            tokens = tokens[:-1]
+
+        allLabels[compIdx] = tokens, probability
+
+    # Convert {comp : [labels]} into a list
+    # of lists, filling in missing components
+    allLabelsList = []
+    probabilities = []
+    for i in range(max(it.chain(allLabels.keys(), noisyComps))):
+        labels, prob = allLabels.get(i + 1, ([missingLabel], math.nan))
+        allLabelsList.append(labels)
+        probabilities.append(prob)
+    allLabels = allLabelsList
+
+    return melDir, noisyComps, allLabels, probabilities
+
+
+def _isfloat(s):
+    """Returns True if the given string appears to contain a floating
+    point number, False otherwise.
+    """
+    try:
+        float(s)
+        return True
+    except Exception:
+        return False
 
 
 def saveLabelFile(allLabels,
                   filename,
                   dirname=None,
                   listBad=True,
-                  signalLabels=None):
+                  signalLabels=None,
+                  probabilities=None):
     """Saves the given classification labels to the specified file. The
     classifications are saved in the format described in the
     :func:`loadLabelFile` method.
 
-    :arg allLabels:    A list of lists, one list for each component, where
-                       each list contains the labels for the corresponding
-                       component.
+    :arg allLabels:     A list of lists, one list for each component, where
+                        each list contains the labels for the corresponding
+                        component.
 
-    :arg filename:     Name of the file to which the labels should be saved.
+    :arg filename:      Name of the file to which the labels should be saved.
 
-    :arg dirname:      If provided, is output as the first line of the file.
-                       Intended to be a relative path to the MELODIC analysis
-                       directory with which this label file is associated. If
-                       not provided, a ``'.'`` is output as the first line.
+    :arg dirname:       If provided, is output as the first line of the file.
+                        Intended to be a relative path to the MELODIC analysis
+                        directory with which this label file is associated. If
+                        not provided, a ``'.'`` is output as the first line.
 
-    :arg listBad:      If ``True`` (the default), the last line of the file
-                       will contain a comma separated list of components which
-                       are deemed 'noisy' (see :func:`isNoisyComponent`).
+    :arg listBad:       If ``True`` (the default), the last line of the file
+                        will contain a comma separated list of components which
+                        are deemed 'noisy' (see :func:`isNoisyComponent`).
 
-    :arg signalLabels: Labels which should be deemed 'signal' - see the
-                       :func:`isNoisyComponent` function.
+    :arg signalLabels:  Labels which should be deemed 'signal' - see the
+                        :func:`isNoisyComponent` function.
+
+    :arg probabilities: Classification probabilities. If provided, the
+                        probability for each component is saved to the file.
     """
 
     lines      = []
     noisyComps = []
 
+    if probabilities is not None and len(probabilities) != len(allLabels):
+        raise ValueError('len(probabilities) != len(allLabels)')
+
     # The first line - the melodic directory name
     if dirname is None:
         dirname = '.'
@@ -283,6 +394,9 @@ def saveLabelFile(allLabels,
         labels = [l.replace(',', '_') for l in labels]
         tokens = [str(comp)] + labels + [str(noise)]
 
+        if probabilities is not None:
+            tokens.append(f'{probabilities[i]:0.6f}')
+
         lines.append(', '.join(tokens))
 
         if noise:
@@ -318,4 +432,3 @@ class InvalidLabelFileError(Exception):
     """Exception raised by the :func:`loadLabelFile` function when an attempt
     is made to load an invalid label file.
     """
-    pass

fsl/data/freesurfer.py

@@ -67,7 +67,8 @@ CORE_GEOMETRY_FILES = ['?h.orig',
                        '?h.pial',
                        '?h.white',
                        '?h.inflated',
-                       '?h.sphere']
+                       '?h.sphere',
+                       '?h.pial_semi_inflated']
 """File patterns for identifying the core Freesurfer geometry files. """
 
 
@@ -76,7 +77,8 @@ CORE_GEOMETRY_DESCRIPTIONS = [
     "Freesurfer surface (pial)",
     "Freesurfer surface (white matter)",
     "Freesurfer surface (inflated)",
-    "Freesurfer surface (sphere)"]
+    "Freesurfer surface (sphere)",
+    "Freesurfer surface (pial semi-inflated)"]
 """A description for each extension in :attr:`GEOMETRY_EXTENSIONS`. """
 
 

fsl/data/gifti.py

@@ -101,9 +101,24 @@ class GiftiMesh(fslmesh.Mesh):
 
         for i, v in enumerate(vertices):
             if i == 0: key = infile
-            else:      key = '{}_{}'.format(infile, i)
+            else:      key = f'{infile}_{i}'
             self.addVertices(v, key, select=(i == 0), fixWinding=fixWinding)
-        self.setMeta(infile, surfimg)
+        self.meta[infile] = surfimg
+
+        # Copy all metadata entries for the GIFTI image
+        for k, v in surfimg.meta.items():
+            self.meta[k] = v
+
+        # and also for each GIFTI data array - triangles
+        # are stored under "faces", and pointsets are
+        # stored under "vertices"/[0,1,2...] (as there may
+        # be multiple pointsets in a file)
+        self.meta['vertices'] = {}
+        for i, arr in enumerate(surfimg.darrays):
+            if arr.intent == constants.NIFTI_INTENT_POINTSET:
+                self.meta['vertices'][i] = dict(arr.meta)
+            elif arr.intent == constants.NIFTI_INTENT_TRIANGLE:
+                self.meta['faces'] = dict(arr.meta)
 
         if vdata is not None:
             self.addVertexData(infile, vdata)
@@ -130,7 +145,7 @@ class GiftiMesh(fslmesh.Mesh):
                     continue
 
                 self.addVertices(vertices[0], sfile, select=False)
-                self.setMeta(sfile, surfimg)
+                self.meta[sfile] = surfimg
 
 
     def loadVertices(self, infile, key=None, *args, **kwargs):
@@ -154,10 +169,10 @@ class GiftiMesh(fslmesh.Mesh):
 
         for i, v in enumerate(vertices):
             if i == 0: key = infile
-            else:      key = '{}_{}'.format(infile, i)
+            else:      key = f'{infile}_{i}'
             vertices[i] = self.addVertices(v, key, *args, **kwargs)
 
-        self.setMeta(infile, surfimg)
+        self.meta[infile] = surfimg
 
         return vertices
 
@@ -221,15 +236,15 @@ def loadGiftiMesh(filename):
     vdata     = [d for d in gimg.darrays if d.intent not in (pscode, tricode)]
 
     if len(triangles) != 1:
-        raise ValueError('{}: GIFTI surface files must contain '
-                         'exactly one triangle array'.format(filename))
+        raise ValueError(f'{filename}: GIFTI surface files must '
+                         'contain exactly one triangle array')
 
     if len(pointsets) == 0:
-        raise ValueError('{}: GIFTI surface files must contain '
-                         'at least one pointset array'.format(filename))
+        raise ValueError(f'{filename}: GIFTI surface files must '
+                         'contain at least one pointset array')
 
     vertices = [ps.data for ps in pointsets]
-    indices  = triangles[0].data
+    indices  = np.atleast_2d(triangles[0].data)
 
     if len(vdata) == 0: vdata = None
     else:               vdata = prepareGiftiVertexData(vdata, filename)
@@ -276,14 +291,14 @@ def prepareGiftiVertexData(darrays, filename=None):
     intents = {d.intent for d in darrays}
 
     if len(intents) != 1:
-        raise ValueError('{} contains multiple (or no) intents'
-                         ': {}'.format(filename, intents))
+        raise ValueError(f'{filename} contains multiple '
+                         f'(or no) intents: {intents}')
 
     intent = intents.pop()
 
     if intent in (constants.NIFTI_INTENT_POINTSET,
                   constants.NIFTI_INTENT_TRIANGLE):
-        raise ValueError('{} contains surface data'.format(filename))
+        raise ValueError(f'{filename} contains surface data')
 
     # Just a single array - return it as-is.
     # n.b. Storing (M, N) data in a single
@@ -298,8 +313,8 @@ def prepareGiftiVertexData(darrays, filename=None):
     vdata = [d.data for d in darrays]
 
     if any([len(d.shape) != 1 for d in vdata]):
-        raise ValueError('{} contains one or more non-vector '
-                         'darrays'.format(filename))
+        raise ValueError(f'{filename} contains one or '
+                         'more non-vector darrays')
 
     vdata = np.vstack(vdata).T
     vdata = vdata.reshape(vdata.shape[0], -1)
@@ -314,6 +329,7 @@ def relatedFiles(fname, ftypes=None):
 
     This function assumes that the GIFTI files are named according to a
     standard convention - the following conventions are supported:
+
      - HCP-style, i.e.: ``<subject>.<hemi>.<type>.<space>.<ftype>.gii``
      - BIDS-style, i.e.:
        ``<source_prefix>_hemi-<hemi>[_space-<space>]*_<suffix>.<ftype>.gii``
@@ -373,7 +389,7 @@ def relatedFiles(fname, ftypes=None):
 
     def searchhcp(match, ftype):
         prefix, space = match
-        template      = '{}.*.{}{}'.format(prefix, space, ftype)
+        template      = f'{prefix}.*.{space}{ftype}'
         return glob.glob(op.join(dirname, template))
 
     # BIDS style - extract all entities (kv

fsl/data/imagewrapper.py

@@ -7,6 +7,9 @@
 """This module provides the :class:`ImageWrapper` class, which can be used
 to manage data access to ``nibabel`` NIFTI images.
 
+.. note:: This module is deprecated - it is being moved to the FSLeyes project,
+          and will be removed in a future version of fslpy.
+
 
 Terminology
 -----------
@@ -42,9 +45,10 @@ import                    collections
 import collections.abc as abc
 import itertools       as it
 
-import numpy     as np
-import nibabel   as nib
+import numpy           as np
 
+import fsl.data.image        as fslimage
+import fsl.utils.deprecated  as deprecated
 import fsl.utils.notifier    as notifier
 import fsl.utils.naninfrange as nir
 import fsl.utils.idle        as idle
@@ -148,6 +152,8 @@ class ImageWrapper(notifier.Notifier):
     """
 
 
+    @deprecated.deprecated('3.9.0', '4.0.0',
+                           'The ImageWrapper has been migrated to FSLeyes')
     def __init__(self,
                  image,
                  name=None,
@@ -175,8 +181,6 @@ class ImageWrapper(notifier.Notifier):
                         data range is updated directly on reads/writes.
         """
 
-        import fsl.data.image as fslimage
-
         self.__image      = image
         self.__name       = name
         self.__taskThread = None
@@ -388,6 +392,14 @@ class ImageWrapper(notifier.Notifier):
             self.__data = np.asanyarray(self.__image.dataobj)
 
 
+    @property
+    def dataIsLoaded(self):
+        """Return true if the image data has been loaded into memory, ``False``
+        otherwise.
+        """
+        return self.__data is not None
+
+
     def __getData(self, sliceobj, isTuple=False):
         """Retrieves the image data at the location specified by ``sliceobj``.
 
@@ -717,104 +729,29 @@ class ImageWrapper(notifier.Notifier):
         self.__updateDataRangeOnWrite(slices, values)
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to fsl.data.image')
 def isValidFancySliceObj(sliceobj, shape):
-    """Returns ``True`` if the given ``sliceobj`` is a valid and fancy slice
-    object.
-
-    ``nibabel`` refers to slice objects as "fancy" if they comprise anything
-    but tuples of integers and simple ``slice`` objects. The ``ImageWrapper``
-    class supports one type of "fancy" slicing, where the ``sliceobj`` is a
-    boolean ``numpy`` array of the same shape as the image.
-
-    This function returns ``True`` if the given ``sliceobj`` adheres to these
-    requirements, ``False`` otherwise.
-    """
-
-    # We only support boolean numpy arrays
-    # which have the same shape as the image
-    return (isinstance(sliceobj, np.ndarray) and
-            sliceobj.dtype == np.bool        and
-            np.prod(sliceobj.shape) == np.prod(shape))
+    """Deprecated - moved to :mod:`fsl.data.image`."""
+    return fslimage.isValidFancySliceObj(sliceobj, shape)
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to fsl.data.image')
 def canonicalSliceObj(sliceobj, shape):
-    """Returns a canonical version of the given ``sliceobj``. See the
-    ``nibabel.fileslice.canonical_slicers`` function.
-    """
-
-    # Fancy slice objects must have
-    # the same shape as the data
-    if isValidFancySliceObj(sliceobj, shape):
-        return sliceobj.reshape(shape)
-
-    else:
-
-        if not isinstance(sliceobj, tuple):
-            sliceobj = (sliceobj,)
-
-        if len(sliceobj) > len(shape):
-            sliceobj = sliceobj[:len(shape)]
-
-        return nib.fileslice.canonical_slicers(sliceobj, shape)
+    """Deprecated - moved to :mod:`fsl.data.image`."""
+    return fslimage.canonicalSliceObj(sliceobj, shape)
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to fsl.data.image')
 def expectedShape(sliceobj, shape):
-    """Given a slice object, and the shape of an array to which
-    that slice object is going to be applied, returns the expected
-    shape of the result.
-
-    .. note:: It is assumed that the ``sliceobj`` has been passed through
-              the :func:`canonicalSliceObj` function.
-
-    :arg sliceobj: Something which can be used to slice an array
-                   of shape ``shape``.
-
-    :arg shape:    Shape of the array being sliced.
-
-    :returns:      A tuple containing:
-
-                     - Expected number of dimensions of the result
-
-                     - Expected shape of the result (or ``None`` if
-                       ``sliceobj`` is fancy).
-    """
-
-    if isValidFancySliceObj(sliceobj, shape):
-        return 1, None
-
-    # Truncate some dimensions from the
-    # slice object if it has too many
-    # (e.g. trailing dims of length 1).
-    elif len(sliceobj) > len(shape):
-        sliceobj = sliceobj[:len(shape)]
-
-    # Figure out the number of dimensions
-    # that the result should have, given
-    # this slice object.
-    expShape = []
-
-    for i in range(len(sliceobj)):
-
-        # Each dimension which has an
-        # int slice will be collapsed
-        if isinstance(sliceobj[i], int):
-            continue
-
-        start = sliceobj[i].start
-        stop  = sliceobj[i].stop
-
-        if start is None: start = 0
-        if stop  is None: stop  = shape[i]
-
-        stop = min(stop, shape[i])
-
-        expShape.append(stop - start)
-
-    return len(expShape), expShape
+    """Deprecated - moved to :mod:`fsl.data.image`."""
+    return fslimage.expectedShape(sliceobj, shape)
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to FSLeyes')
 def sliceObjToSliceTuple(sliceobj, shape):
-    """Turns an array slice object into a tuple of (low, high) index
+    """Deprecated - the imagewrapper has been moved to FSLeyes.
+
+    Turns an array slice object into a tuple of (low, high) index
     pairs, one pair for each dimension in the given shape
 
     :arg sliceobj: Something which can be used to slice an array of shape
@@ -853,8 +790,11 @@ def sliceObjToSliceTuple(sliceobj, shape):
     return tuple(indices)
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to FSLeyes')
 def sliceTupleToSliceObj(slices):
-    """Turns a sequence of (low, high) index pairs into a tuple of array
+    """Deprecated - the imagewrapper has been moved to FSLeyes.
+
+    Turns a sequence of (low, high) index pairs into a tuple of array
     ``slice`` objects.
 
     :arg slices: A sequence of (low, high) index pairs.
@@ -868,8 +808,11 @@ def sliceTupleToSliceObj(slices):
     return tuple(sliceobj)
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to FSLeyes')
 def adjustCoverage(oldCoverage, slices):
-    """Adjusts/expands the given ``oldCoverage`` so that it covers the
+    """Deprecated - the imagewrapper has been moved to FSLeyes.
+
+    Adjusts/expands the given ``oldCoverage`` so that it covers the
     given set of ``slices``.
 
     :arg oldCoverage: A ``numpy`` array of shape ``(2, n)`` containing
@@ -916,8 +859,11 @@ return code for the :func:`sliceOverlap` function.
 """
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to FSLeyes')
 def sliceOverlap(slices, coverage):
-    """Determines whether the given ``slices`` overlap with the given
+    """Deprecated - the imagewrapper has been moved to FSLeyes.
+
+    Determines whether the given ``slices`` overlap with the given
     ``coverage``.
 
     :arg slices:    A sequence of (low, high) index pairs, assumed to cover
@@ -983,8 +929,11 @@ def sliceOverlap(slices, coverage):
     elif np.all(overlapStates == OVERLAP_ALL):  return OVERLAP_ALL
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to FSLeyes')
 def sliceCovered(slices, coverage):
-    """Returns ``True`` if the portion of the image data calculated by
+    """Deprecated - the imagewrapper has been moved to FSLeyes.
+
+    Returns ``True`` if the portion of the image data calculated by
     the given ``slices` has already been calculated, ``False`` otherwise.
 
     :arg slices:    A sequence of (low, high) index pairs, assumed to cover
@@ -1015,8 +964,11 @@ def sliceCovered(slices, coverage):
     return True
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to FSLeyes')
 def calcExpansion(slices, coverage):
-    """Calculates a series of *expansion* slices, which can be used to expand
+    """Deprecated - the imagewrapper has been moved to FSLeyes.
+
+    Calculates a series of *expansion* slices, which can be used to expand
     the given ``coverage`` so that it includes the given ``slices``.
 
     :arg slices:   Slices that the coverage needs to be expanded to cover.
@@ -1185,8 +1137,11 @@ def calcExpansion(slices, coverage):
     return volumes, expansions
 
 
+@deprecated.deprecated('3.9.0', '4.0.0', 'Moved to FSLeyes')
 def collapseExpansions(expansions, numDims):
-    """Scans through the given list of expansions (each assumed to pertain
+    """Deprecated - the imagewrapper has been moved to FSLeyes.
+
+    Scans through the given list of expansions (each assumed to pertain
     to a single 3D image), and combines any which cover the same
     image area, and cover adjacent volumes.
 

fsl/data/melodicanalysis.py

@@ -33,7 +33,6 @@ import logging
 import os.path as op
 import numpy   as np
 
-import fsl.utils.path        as fslpath
 import fsl.data.image        as fslimage
 import fsl.data.featanalysis as featanalysis
 
@@ -63,10 +62,9 @@ def isMelodicImage(path):
 
 
 def isMelodicDir(path):
-    """Returns ``True`` if the given path looks like it is contained within
-    a MELODIC directory, ``False`` otherwise. A melodic directory:
+    """Returns ``True`` if the given path looks like it is a MELODIC directory,
+    ``False`` otherwise. A MELODIC directory:
 
-      - Must be named ``*.ica``.
       - Must contain a file called ``melodic_IC.nii.gz`` or
         ``melodic_oIC.nii.gz``.
       - Must contain a file called ``melodic_mix``.
@@ -75,12 +73,7 @@ def isMelodicDir(path):
 
     path = op.abspath(path)
 
-    if op.isdir(path): dirname = path
-    else:              dirname = op.dirname(path)
-
-    sufs = ['.ica']
-
-    if not any([dirname.endswith(suf) for suf in sufs]):
+    if not op.isdir(path):
         return False
 
     # Must contain an image file called
@@ -88,7 +81,7 @@ def isMelodicDir(path):
     prefixes = ['melodic_IC', 'melodic_oIC']
     for p in prefixes:
         try:
-            fslimage.addExt(op.join(dirname, p))
+            fslimage.addExt(op.join(path, p))
             break
         except fslimage.PathError:
             pass
@@ -97,8 +90,8 @@ def isMelodicDir(path):
 
     # Must contain files called
     # melodic_mix and melodic_FTmix
-    if not op.exists(op.join(dirname, 'melodic_mix')):   return False
-    if not op.exists(op.join(dirname, 'melodic_FTmix')): return False
+    if not op.exists(op.join(path, 'melodic_mix')):   return False
+    if not op.exists(op.join(path, 'melodic_FTmix')): return False
 
     return True
 
@@ -108,10 +101,13 @@ def getAnalysisDir(path):
     to that MELODIC directory is returned. Otherwise, ``None`` is returned.
     """
 
-    meldir = fslpath.deepest(path, ['.ica'])
+    if not op.isdir(path):
+        path = op.dirname(path)
 
-    if meldir is not None and isMelodicDir(meldir):
-        return meldir
+    while path not in (op.sep, ''):
+        if isMelodicDir(path):
+            return path
+        path = op.dirname(path)
 
     return None
 
@@ -137,10 +133,18 @@ def getDataFile(meldir):
     if topDir is None:
         return None
 
-    dataFile = op.join(topDir, 'filtered_func_data')
+    # People often rename filtered_func_data.nii.gz
+    # to something like filtered_func_data_clean.nii.gz,
+    # because that is the recommended approach when
+    # performing ICA-based denoising). So we try both.
+    candidates = ['filtered_func_data', 'filtered_func_data_clean']
 
-    try:                       return fslimage.addExt(dataFile)
-    except fslimage.PathError: return None
+    for candidate in candidates:
+        dataFile = op.join(topDir, candidate)
+        try:                       return fslimage.addExt(dataFile)
+        except fslimage.PathError: continue
+
+    return None
 
 
 def getMeanFile(meldir):
@@ -187,7 +191,7 @@ def getNumComponents(meldir):
     contained in the given directrory.
     """
 
-    icImg = fslimage.Image(getICFile(meldir), loadData=False, calcRange=False)
+    icImg = fslimage.Image(getICFile(meldir))
     return icImg.shape[3]
 
 

fsl/data/melodicimage.py

@@ -74,9 +74,7 @@ class MelodicImage(fslimage.Image):
         dataFile = self.getDataFile()
 
         if dataFile is not None:
-            dataImage = fslimage.Image(dataFile,
-                                       loadData=False,
-                                       calcRange=False)
+            dataImage = fslimage.Image(dataFile)
             if dataImage.ndim >= 4:
                 self.__tr = dataImage.pixdim[3]
 

fsl/data/mesh.py

@@ -41,6 +41,13 @@ import fsl.transform.affine as affine
 log = logging.getLogger(__name__)
 
 
+class IncompatibleVerticesError(ValueError):
+    """``ValueError`` raised by the :meth:`Mesh.addVertices` method if
+    an attempt is made to add a vertex set with the wrong number of
+    vertices.
+    """
+
+
 class Mesh(notifier.Notifier, meta.Meta):
     """The ``Mesh`` class represents a 3D model. A mesh is defined by a
     collection of ``N`` vertices, and ``M`` triangles.  The triangles are
@@ -155,7 +162,7 @@ class Mesh(notifier.Notifier, meta.Meta):
 
 
     def __init__(self,
-                 indices,
+                 indices=None,
                  name='mesh',
                  dataSource=None,
                  vertices=None,
@@ -166,7 +173,8 @@ class Mesh(notifier.Notifier, meta.Meta):
         :meth:`addVertices` method.
 
         :arg indices:    A list of indices into the vertex data, defining the
-                         mesh triangles.
+                         mesh triangles. If not provided, must be provided
+                         after creation via the :meth:`indices` setter method.
 
         :arg name:       A name for this ``Mesh``.
 
@@ -179,22 +187,31 @@ class Mesh(notifier.Notifier, meta.Meta):
                          :meth:`addVertices` method along with ``vertices``.
         """
 
+        if indices is None and vertices is not None:
+            raise ValueError('Indices must be provided '
+                             'if vertices are provided')
+
         self.__name       = name
         self.__dataSource = dataSource
-        self.__nvertices  = int(indices.max()) + 1
-        self.__selected   = None
+
+        # nvertices/indices are assigned in the
+        # indices setter method.
 
         # We potentially store two copies of
-        # the indices, with opposite unwinding
-        # orders. The vindices dict stores refs
-        # to one or the other for each vertex
-        # set.
-        self.__indices      = np.asarray(indices, dtype=np.int32).reshape((-1, 3))
+        # the indices, - one set (__indices)
+        # as provided, and the other
+        # (__fixedIndices) with opposite
+        # unwinding orders. The vindices dict
+        # stores refs to one or the other for
+        # each vertex set.
+        self.__nvertices    = None
+        self.__indices      = None
         self.__fixedIndices = None
         self.__vindices     = collections.OrderedDict()
 
         # All of these are populated
         # in the addVertices method
+        self.__selected = None
         self.__vertices = collections.OrderedDict()
         self.__loBounds = collections.OrderedDict()
         self.__hiBounds = collections.OrderedDict()
@@ -212,8 +229,9 @@ class Mesh(notifier.Notifier, meta.Meta):
         # in the trimesh method
         self.__trimesh = collections.OrderedDict()
 
-        # Add initial vertex
-        # set if provided
+        # Add initial indices/vertices if provided
+        if indices is not None:
+            self.indices = indices
         if vertices is not None:
             self.addVertices(vertices, fixWinding=fixWinding)
 
@@ -282,10 +300,25 @@ class Mesh(notifier.Notifier, meta.Meta):
 
     @property
     def indices(self):
-        """The ``(M, 3)`` triangles of this mesh. """
+        """The ``(M, 3)`` triangles of this mesh. Returns ``None`` if
+        indices have not yet been assigned.
+        """
+        if self.__indices is None:
+            return None
         return self.__vindices[self.__selected]
 
 
+    @indices.setter
+    def indices(self, indices):
+        """Set the indices for this mesh. """
+        if self.__indices is not None:
+            raise ValueError('Indices are already set')
+
+        indices          = np.asarray(indices, dtype=np.int32)
+        self.__nvertices = int(indices.max()) + 1
+        self.__indices   = indices.reshape((-1, 3))
+
+
     @property
     def normals(self):
         """A ``(M, 3)`` array containing surface normals for every
@@ -329,7 +362,13 @@ class Mesh(notifier.Notifier, meta.Meta):
         ``Mesh`` instance. The bounding box is arranged like so:
 
             ``((xlow, ylow, zlow), (xhigh, yhigh, zhigh))``
+
+        Returns ``None`` if indices or vertices have not yet been assigned.
         """
+
+        if self.__indices is None or len(self.__vertices) == 0:
+            return None
+
         lo = self.__loBounds[self.__selected]
         hi = self.__hiBounds[self.__selected]
         return lo, hi
@@ -380,10 +419,13 @@ class Mesh(notifier.Notifier, meta.Meta):
 
         :returns:        The vertices, possibly reshaped
 
-        :raises:         ``ValueError`` if the provided ``vertices`` array
-                         has the wrong number of vertices.
+        :raises:         ``IncompatibleVerticesError`` if the provided
+                         ``vertices`` array has the wrong number of vertices.
         """
 
+        if self.__indices is None:
+            raise ValueError('Mesh indices have not yet been set')
+
         if key is None:
             key = 'default'
 
@@ -399,10 +441,9 @@ class Mesh(notifier.Notifier, meta.Meta):
         # reshape raised an error -
         # wrong number of vertices
         except ValueError:
-            raise ValueError('{}: invalid number of vertices: '
-                             '{} != ({}, 3)'.format(key,
-                                                    vertices.shape,
-                                                    self.nvertices))
+            raise IncompatibleVerticesError(
+                f'{key}: invalid number of vertices: '
+                f'{vertices.shape} != ({self.nvertices}, 3)')
 
         self.__vertices[key] = vertices
         self.__vindices[key] = self.__indices
@@ -573,7 +614,7 @@ class Mesh(notifier.Notifier, meta.Meta):
 
         # sort by ray. I'm Not sure if this is
         # needed - does trimesh do it for us?
-        rayIdxs = np.asarray(np.argsort(rays), np.int)
+        rayIdxs = np.asarray(np.argsort(rays))
         locs    = locs[rayIdxs]
         tris    = tris[rayIdxs]
 
@@ -687,7 +728,7 @@ def calcFaceNormals(vertices, indices):
     fnormals = np.cross((v1 - v0), (v2 - v0))
     fnormals = affine.normalise(fnormals)
 
-    return fnormals
+    return np.atleast_2d(fnormals)
 
 
 def calcVertexNormals(vertices, indices, fnormals):
@@ -701,7 +742,7 @@ def calcVertexNormals(vertices, indices, fnormals):
                    the mesh.
     """
 
-    vnormals = np.zeros((vertices.shape[0], 3), dtype=np.float)
+    vnormals = np.zeros((vertices.shape[0], 3), dtype=float)
 
     # TODO make fast. I can't figure
     # out how to use np.add.at to