doc/fsl.wrappers.fsl_sub.rst

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

doc/fsl.wrappers.fslstats.rst

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

doc/fsl.wrappers.oxford_asl.rst

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

doc/fsl.wrappers.randomise.rst

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

doc/fsl.wrappers.rst

@@ -4,16 +4,29 @@
 .. toctree::
    :hidden:
 
+   fsl.wrappers.avwutils
+   fsl.wrappers.bedpostx
    fsl.wrappers.bet
+   fsl.wrappers.bianca
+   fsl.wrappers.cluster_commands
+   fsl.wrappers.dtifit
    fsl.wrappers.eddy
+   fsl.wrappers.epi_reg
    fsl.wrappers.fast
+   fsl.wrappers.feat
+   fsl.wrappers.first
    fsl.wrappers.flirt
    fsl.wrappers.fnirt
-   fsl.wrappers.fslmaths
    fsl.wrappers.fsl_anat
+   fsl.wrappers.fsl_sub
+   fsl.wrappers.fslmaths
+   fsl.wrappers.fslstats
    fsl.wrappers.fugue
    fsl.wrappers.melodic
    fsl.wrappers.misc
+   fsl.wrappers.oxford_asl
+   fsl.wrappers.randomise
+   fsl.wrappers.tbss
    fsl.wrappers.wrapperutils
 
 .. automodule:: fsl.wrappers

doc/fsl.wrappers.tbss.rst

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

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.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

@@ -36,27 +36,28 @@ load an atlas image, which will be one of the following atlas-specific
    :nosignatures:
 
    LabelAtlas
+   StatisticAtlas
    ProbabilisticAtlas
 """
 
 
 from __future__ import division
 
-import xml.etree.ElementTree              as et
-import os.path                            as op
-import                                       glob
-import                                       bisect
-import                                       logging
+import xml.etree.ElementTree    as et
+import os.path                  as op
+import                             glob
+import                             bisect
+import                             logging
 
-import numpy                              as np
+import numpy                    as np
 
-import fsl.data.image                     as fslimage
-import fsl.data.constants                 as constants
-from   fsl.utils.platform import platform as platform
-import fsl.utils.image.resample           as resample
-import fsl.transform.affine               as affine
-import fsl.utils.notifier                 as notifier
-import fsl.utils.settings                 as fslsettings
+import fsl.data.image           as fslimage
+import fsl.data.constants       as constants
+from   fsl.utils.platform import   platform
+import fsl.utils.image.resample as resample
+import fsl.transform.affine     as affine
+import fsl.utils.notifier       as notifier
+import fsl.utils.settings       as fslsettings
 
 
 log = logging.getLogger(__name__)
@@ -322,9 +323,9 @@ class AtlasLabel(object):
     ========= ================================================================
     ``name``  Region name
     ``index`` The index of this label into the list of all labels in the
-              ``AtlasDescription`` that owns it. For probabilistic atlases,
-              this is also the index into the 4D atlas image of the volume
-              that corresponds to this region.
+              ``AtlasDescription`` that owns it. For statistic/probabilistic
+              atlases, this is also the index into the 4D atlas image of the
+              volume that corresponds to this region.
     ``value`` For label atlases and summary images, the value of voxels that
               are in this region.
     ``x``     X coordinate of the region in world space
@@ -376,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
@@ -386,8 +387,13 @@ class AtlasDescription(object):
 
        <atlas>
          <header>
-           <name></name>        # Atlas name
-           <type></type>        # 'Probabilistic' or 'Label'
+           <name></name>           # Atlas name
+           <type></type>           # 'Statistic', 'Probabilistic' or 'Label'
+           <statistic></statistic> # Optional. Type of statistic
+           <units></units>         # Optional. Units of measurement
+           <precision></precision> # Optional. Decimal precision to report
+           <upper></upper>         # Optional. Upper threshold
+           <lower></lower>         # Optional. Lower threshold
            <images>
             <imagefile>
             </imagefile>        # If type is Probabilistic, path
@@ -412,11 +418,12 @@ class AtlasDescription(object):
          </header>
         <data>
 
-         # index - For probabilistic atlases, index of corresponding volume in
-         #         4D image file. For label images, the value of voxels which
-         #         are in the corresponding region. For probabilistic atlases,
-         #         it is assumed that the value for each region in the  summary
-         #         image(s) are equal to ``index + 1``.
+         # index - For statistic/probabilistic atlases, index of corresponding
+         #         volume in 4D image file. For label images, the value of
+         #         voxels which are in the corresponding region. For
+         #         statistic/probabilistic atlases, it is assumed that the
+         #         value for each region in the summary image(s) are equal to
+         #         ``index + 1``.
          #
          #
          # x    |
@@ -452,7 +459,18 @@ class AtlasDescription(object):
 
     ``specPath``      Path to the atlas XML specification file.
 
-    ``atlasType``     Atlas type - either *probabilistic* or *label*.
+    ``atlasType``     Atlas type - either *statistic*, *probabilistic* or
+                      *label*.
+
+    ``statistic``     Type of statistic, for statistic atlases.
+
+    ``units``         Unit of measurement, for statistic atlases.
+
+    ``precision``     Reporting precision, for statistic atlases.
+
+    ``upper``         Upper threshold, for statistic atlases.
+
+    ``lower``         Lower threshold, for statistic atlases.
 
     ``images``        A list of images available for this atlas - usually
                       :math:`1mm^3` and :math:`2mm^3` images are present.
@@ -500,6 +518,29 @@ class AtlasDescription(object):
         if self.atlasType == 'probabalistic':
             self.atlasType = 'probabilistic'
 
+        if self.atlasType == 'statistic':
+
+            fields = ['statistic', 'units', 'lower', 'upper', 'precision']
+            values = {}
+
+            for field in fields:
+                elem = header.find(field)
+                if elem is not None and elem.text is not None:
+                    values[field] = elem.text.strip()
+
+            self.statistic =       values.get('statistic', '')
+            self.units     =       values.get('units',     '')
+            self.lower     = float(values.get('lower',     0))
+            self.upper     = float(values.get('upper',     100))
+            self.precision = int(  values.get('precision', 2))
+
+        elif self.atlasType == 'probabilistic':
+            self.statistic = ''
+            self.units     = '%'
+            self.lower     = 5
+            self.upper     = 100
+            self.precision = 0
+
         images             = header.findall('images')
         self.images        = []
         self.summaryImages = []
@@ -519,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)
@@ -661,7 +702,7 @@ class Atlas(fslimage.Image):
         :arg resolution: Desired isotropic resolution in millimetres.
 
         :arg isLabel:    Pass in ``True`` for label atlases, ``False`` for
-                         probabilistic atlases.
+                         statistic/probabilistic atlases.
 
         All other arguments are passed to :meth:`.Image.__init__`.
         """
@@ -708,7 +749,7 @@ class Atlas(fslimage.Image):
         """Makes sure that the given mask has the same resolution as this
         atlas, so it can be used for querying. Used by the
         :meth:`.LabelAtlas.maskLabels` and
-        :meth:`.ProbabilisticAtlas.maskProportions` methods.
+        :meth:`.StatisticAtlas.maskValues` methods.
 
         :arg mask: A :class:`.Image`
 
@@ -738,13 +779,11 @@ class Atlas(fslimage.Image):
         return mask
 
 
-
 class MaskError(Exception):
     """Exception raised by the :meth:`LabelAtlas.maskLabel` and
-    :meth:`ProbabilisticAtlas.maskProportions` when a mask is provided which
+    :meth:`StatisticAtlas.maskValues` when a mask is provided which
     does not match the atlas space.
     """
-    pass
 
 
 class LabelAtlas(Atlas):
@@ -841,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
@@ -877,17 +923,20 @@ class LabelAtlas(Atlas):
 
         return values, props
 
-    def get(self, label=None, index=None, value=None, name=None):
-        """
-        Returns the binary image for given label
+
+    def get(self, label=None, index=None, value=None, name=None, binary=True):
+        """Returns the binary image for the given label.
 
         Only one of the arguments should be used to define the label
 
-        :arg label: AtlasLabel contained within this atlas
-        :arg index: index of the label
-        :arg value: value of the label
-        :arg name: string of the label
-        :return: image.Image with the mask
+        :arg label:  :class:`AtlasLabel` contained within this atlas
+        :arg index:  index of the label
+        :arg value:  value of the label
+        :arg name:   string of the label
+        :arg binary: If ``True`` (the default), the image will contain 1s in
+                     the label region. Otherwise the image will contain the
+                     label value.
+        :return:     :class:`.Image` with the mask
         """
         if ((label is not None) + (index is not None) +
             (value is not None) + (name is not None)) != 1:
@@ -896,19 +945,27 @@ class LabelAtlas(Atlas):
             label = self.find(index=index, name=name, value=value)
         elif label not in self.desc.labels:
             raise ValueError("Unknown label provided")
-        arr = (self.data == label.value).astype(int)
+
+        arr = (self.data == label.value).astype(np.int32)
+
+        if not binary:
+            arr[arr > 0] = label.value
+
         return fslimage.Image(arr, name=label.name, header=self.header)
 
 
-class ProbabilisticAtlas(Atlas):
-    """A 4D atlas which contains one volume for each region.
+class StatisticAtlas(Atlas):
+    """A ``StatisticAtlas`` is a 4D image which contains one volume for
+    each region in the atlas; each volume contains some statistic value
+    for the corresponding region.
 
-    The ``ProbabilisticAtlas`` provides the :meth`proportions` method,
-    which makes looking up region probabilities easy.
+    The :class:`ProbabilisticAtlas` is a specialisation of the
+    ``StatisticAtlas``
     """
 
+
     def __init__(self, atlasDesc, resolution=None, **kwargs):
-        """Create a ``ProbabilisticAtlas`` instance.
+        """Create a ``StatisticAtlas`` instance.
 
         :arg atlasDesc:  The :class:`AtlasDescription` instance describing
                          the atlas.
@@ -917,17 +974,18 @@ class ProbabilisticAtlas(Atlas):
         """
         Atlas.__init__(self, atlasDesc, resolution, False, **kwargs)
 
+
     def get(self, label=None, index=None, value=None, name=None):
-        """
-        Returns the probabilistic image for given label
+        """Returns the statistic image at the given label.
 
         Only one of the arguments should be used to define the label
 
-        :arg label: AtlasLabel contained within this atlas
+        :arg label: :class:`AtlasLabel` contained within this atlas
         :arg index: index of the label
         :arg value: value of the label
-        :arg name: string of the label
-        :return: image.Image with the probabilistic mask
+        :arg name:  string of the label
+        :return:    :class:`.Image` with the statistic values for the
+                    specified label.
         """
         if ((label is not None) + (index is not None) +
             (value is not None) + (name is not None)) != 1:
@@ -939,36 +997,37 @@ class ProbabilisticAtlas(Atlas):
         arr = self[..., label.index]
         return fslimage.Image(arr, name=label.name, header=self.header)
 
-    def proportions(self, location, *args, **kwargs):
-        """Looks up and returns the proportions of of all regions at the given
+
+    def values(self, location, *args, **kwargs):
+        """Looks up and returns the values of of all regions at the given
         location.
 
         :arg location: Can be one of the following:
 
                         - A sequence of three values, interpreted as atlas
-                          coordinates. In this case, :meth:`coordProportions`
+                          coordinates. In this case, :meth:`coordValues`
                           is called.
 
                         - An :class:`.Image` which is interpreted as a
-                          weighted mask. In this case, :meth:`maskProportions`
+                          weighted mask. In this case, :meth:`maskValues`
                           is called.
 
-        All other arguments are passed through to the :meth:`coordProportions`
-        or :meth:`maskProportions` methods.
+        All other arguments are passed through to the :meth:`coordValues`
+        or :meth:`maskValues` methods.
 
 
-        :returns: The return value of either :meth:`coordProportions` or
-                  :meth:`maskProportions`.
+        :returns: The return value of either :meth:`coordValues` or
+                  :meth:`maskValues`.
         """
 
         if isinstance(location, fslimage.Image):
-            return self.maskProportions(location, *args, **kwargs)
+            return self.maskValues(location, *args, **kwargs)
         else:
-            return self.coordProportions(location, *args, **kwargs)
+            return self.coordValues(location, *args, **kwargs)
 
 
-    def coordProportions(self, loc, voxel=False):
-        """Looks up the region probabilities for the given location.
+    def coordValues(self, loc, voxel=False):
+        """Looks up the region values for the given location.
 
         :arg loc:   A sequence of three values, interpreted as atlas
                     world or voxel coordinates.
@@ -976,10 +1035,8 @@ class ProbabilisticAtlas(Atlas):
         :arg voxel: Defaults to ``False``. If ``True``, the ``loc``
                     argument is interpreted as voxel coordinates.
 
-        :returns: a list of values, one per region, which represent
-                  the probability of each region for the specified
-                  location. Returns an empty list if the given
-                  location is out of bounds.
+        :returns: a list of values, one per region.  Returns an empty
+                  list if the given location is out of bounds.
         """
 
         if not voxel:
@@ -994,30 +1051,27 @@ class ProbabilisticAtlas(Atlas):
            loc[2] >= self.shape[2]:
             return []
 
-        props = self[loc[0], loc[1], loc[2], :]
+        vals = self[loc[0], loc[1], loc[2], :]
 
         # We only return labels for this atlas -
         # the underlying image may have more
         # volumes than this atlas has labels.
-        return [props[l.index] for l in self.desc.labels]
+        return [vals[l.index] for l in self.desc.labels]
 
 
-    def maskProportions(self, mask):
-        """Looks up the probabilities of all regions in the given ``mask``.
+    def maskValues(self, mask):
+        """Looks up the average values of all regions in the given ``mask``.
 
         :arg mask: A 3D :class:`.Image`` which is interpreted as a weighted
                    mask. If the ``mask`` shape does not match that of this
-                   ``ProbabilisticAtlas``, it is resampled using
-                   :meth:`.Image.resample`, with nearest-neighbour
-                   interpolation.
+                   ``StatisticAtlas``, it is resampled using
+                   :meth:`Atlas.prepareMask`.
 
-        :returns:  A sequence containing the proportion, within the mask,
-                   of all regions in the atlas. The proportions are returned as
-                   values between 0 and 100.
+        :returns:  A sequence containing the average value, within the mask,
+                   of all regions in the atlas.
         """
 
-        props = []
-
+        avgvals   = []
         mask      = self.prepareMask(mask)
         boolmask  = mask > 0
         weights   = mask[boolmask]
@@ -1030,11 +1084,17 @@ class ProbabilisticAtlas(Atlas):
 
             vals  = self[..., label.index]
             vals  = vals[boolmask] * weights
-            prop  = vals.sum() / weightsum
+            val   = vals.sum() / weightsum
 
-            props.append(prop)
+            avgvals.append(val)
 
-        return props
+        return avgvals
+
+
+class ProbabilisticAtlas(StatisticAtlas):
+    """A 4D atlas which contains one volume for each region. Each volume
+    contains probabiliy values for one region, between 0 and 100.
+    """
 
 
 registry            = AtlasRegistry()

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

@@ -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]