fsl/data/bitmap.py

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

fsl/data/constants.py

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

fsl/data/dicom.py

@@ -29,22 +29,56 @@ import               os
 import os.path    as op
 import subprocess as sp
 import               re
+import               sys
 import               glob
 import               json
+import               shlex
+import               shutil
 import               logging
+import               binascii
 
+import numpy      as np
 import nibabel    as nib
 
-import fsl.utils.tempdir as tempdir
-import fsl.utils.memoize as memoize
-import fsl.data.image    as fslimage
+import fsl.utils.tempdir  as tempdir
+import fsl.utils.memoize  as memoize
+import fsl.utils.platform as fslplatform
+import fsl.data.image     as fslimage
 
 
 log = logging.getLogger(__name__)
 
 
 MIN_DCM2NIIX_VERSION = (1, 0, 2017, 12, 15)
-"""Minimum version of dcm2niix that is required for this module to work. """
+"""Minimum version of ``dcm2niix`` that is required for this module to work.
+"""
+
+
+CRC_DCM2NIIX_VERSION = (1, 0, 2019, 9, 2)
+"""For versions of ``dcm2niix`` orf this version or newer, the ``-n`` flag,
+used to convert a single DICOM series, requires that a CRC checksum
+identifying the series be passed (see the :func:`seriesCRC`
+function). Versions prior to this require the series number to be passed.
+"""
+
+
+def dcm2niix() -> str:
+    """Tries to find an absolute path to the ``dcm2niix`` command. Returns
+    ``'dcm2niix'`` (unqualified) if a specific executable cannot be found.
+    """
+    fsldir     = fslplatform.platform.fsldir
+    candidates = [
+        shutil.which('dcm2niix')
+    ]
+
+    if fsldir is not None:
+        candidates.insert(0, op.join(fsldir, 'bin', 'dcm2niix'))
+
+    for c in candidates:
+        if c is not None and op.exists(c):
+            return c
+
+    return 'dcm2niix'
 
 
 class DicomImage(fslimage.Image):
@@ -80,12 +114,19 @@ class DicomImage(fslimage.Image):
 
 
 @memoize.memoize
-def enabled():
-    """Returns ``True`` if ``dcm2niix`` is present, and recent enough,
-    ``False`` otherwise.
+def installedVersion():
+    """Return a tuple describing the version of ``dcm2niix`` that is installed,
+    or ``None`` if dcm2niix cannot be found, or its version not parsed.
+
+    The returned tuple contains the following fields, all integers:
+      - Major version number
+      - Minor version number
+      - Year
+      - Month
+      - Day
     """
 
-    cmd            = 'dcm2niix -h'
+    cmd            = f'{dcm2niix()} -h'
     versionPattern = re.compile(r'v'
                                 r'(?P<major>[0-9]+)\.'
                                 r'(?P<minor>[0-9]+)\.'
@@ -102,80 +143,124 @@ def enabled():
 
             match = re.match(versionPattern, word)
 
-            if match is None:
-                continue
+            if match is not None:
+                return (int(match.group('major')),
+                        int(match.group('minor')),
+                        int(match.group('year')),
+                        int(match.group('month')),
+                        int(match.group('day')))
 
-            installedVersion = (
-                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
 
-            # 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
 
-            # if we get here, versions are equal
-            return True
+def compareVersions(v1, v2):
+    """Compares two ``dcm2niix`` versions ``v1`` and ``v2``.  The versions are
+    assumed to be in the format returned by :func:`installedVersion`.
 
-    except Exception as e:
-        log.debug('Error parsing dcm2niix version string: {}'.format(e))
+    :returns: - 1 if ``v1`` is newer than ``v2``
+              - -1 if ``v1`` is older than ``v2``
+              - 0 if ``v1`` the same as ``v2``.
+    """
 
-    return False
+    for iv1, iv2 in zip(v1, v2):
+        if   iv1 > iv2: return  1
+        elif iv1 < iv2: return -1
+    return 0
+
+
+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.
+
+    Some additional metadata is added to each dictionary:
+     - ``DicomDir``: The absolute path to ``dcmdir``
 
-    :arg dcmdir: Directory containing DICOM files.
+    :arg dcmdir: Directory containing DICOM series
 
-    :returns:    A list of dictionaries, each containing metadata about
-                 one DICOM data 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,20 +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()
 
-    with tempdir.tempdir() as td:
+    # Newer versions of dcm2niix
+    # require a CRC to identify
+    # series
+    if compareVersions(version, CRC_DCM2NIIX_VERSION) >= 0:
+        ident = seriesCRC(series)
 
-        with open(os.devnull, 'wb') as devnull:
-            sp.call(cmd.split(), stdout=devnull, stderr=devnull)
+    # Older versions require
+    # the series number
+    else:
+        ident = snum
 
-        files  = glob.glob(op.join(td, '{}*.nii'.format(snum)))
-        images = [nib.load(f) for f in files]
+    cmd = f'{dcm2niix()} -b n -f %s -z n -o . -n "{ident}" "{dcmdir}"'
 
-        # Force-load images into memory
-        [i.get_data() for i in images]
+    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, 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 force-load the
+        # the image data into memory (to
+        # avoid any disk accesses due to
+        # e.g. memmap)
+        images = [nib.Nifti1Image(np.asanyarray(i.dataobj), None, i.header)
+                  for i in images]
 
         return [DicomImage(i, series, dcmdir, name=desc) for i in images]

fsl/data/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`. """