Renamed TensorImage to DTIFitTensor, and tensorimage.py to

dtifit.py. Added new function, decomposeTensorMatrix, which turns a 6
volume tensor matrix image into eigenvalue and eigenvector images.

fsl/data/tensorimage.py → fsl/data/dtifit.py

 #!/usr/bin/env python
 #
-# tensorimage.py - The TensorImage class.
+# dtifit.py - The DTIFitTensor class, and some related utility functions.
 #
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
-"""This module provides the :class:`.TensorImage` class, which encapsulates
-the diffusion tensor data generated by the FSL ``dtifit`` tool.
+"""This module provides the :class:`.DTIFitTensor` class, which encapsulates
+the diffusion tensor data generated by the FSL ``dtifit`` tool. 
+
+The following utility functions are also defined:
+
+  .. autosummary::
+     :nosignatures:
+
+     getDTIFitDataPrefix
+     isDTIFitPath
+     looksLikeTensorImage
+     decomposeTensorMatrix
 """
 
 
@@ -14,7 +24,8 @@ import                 re
 import                 glob
 import os.path      as op
 
-import                 six
+import numpy        as np
+import numpy.linalg as npla
 
 from . import image as fslimage
 
@@ -22,9 +33,10 @@ from . import image as fslimage
 log = logging.getLogger(__name__)
 
 
-def getTensorDataPrefix(path):
-    """Returns the prefix used for the DTI file names in the given
-    directory, or ``None`` if the DTI files could not be identified.
+def getDTIFitDataPrefix(path):
+    """Returns the prefix (a.k,a, base name) used for the ``dtifit`` file 
+    names in the given directory, or ``None`` if the ``dtifit`` files could 
+    not be identified.
     """
     
     v1s   = glob.glob(op.join(path, '*_V1.*'))
@@ -33,14 +45,12 @@ def getTensorDataPrefix(path):
     l1s   = glob.glob(op.join(path, '*_L1.*'))
     l2s   = glob.glob(op.join(path, '*_L2.*'))
     l3s   = glob.glob(op.join(path, '*_L3.*'))
-    fas   = glob.glob(op.join(path, '*_FA.*'))
-    mds   = glob.glob(op.join(path, '*_MD.*'))
-    files = [v1s, v2s, v3s, l1s, l2s, l3s, fas, mds]
+    files = [v1s, v2s, v3s, l1s, l2s, l3s]
 
     # Gather all of the existing file
     # prefixes into a dictionary of
     # prefix : [file list] mappings.
-    pattern  = '^(.*)_(?:V1|V2|V3|L1|L2|L3|FA|MD).*$'
+    pattern  = '^(.*)_(?:V1|V2|V3|L1|L2|L3).*$'
     prefixes = {}
 
     for f in [f for flist in files for f in flist]:
@@ -52,7 +62,7 @@ def getTensorDataPrefix(path):
     # Discard any prefixes which are 
     # not present for every file type.
     for prefix, files in list(prefixes.items()):
-        if len(files) != 8:
+        if len(files) != 6:
             prefixes.pop(prefix)
 
     # Discard any prefixes which
@@ -79,68 +89,110 @@ def getTensorDataPrefix(path):
     return op.basename(prefixes[0])
 
 
-def isPathToTensorData(path):
+def isDTIFitPath(path):
     """Returns ``True`` if the given directory path looks like it contains
     ``dtifit`` data, ``False`` otherwise.
     """
     
-    return getTensorDataPrefix(path) is not None
+    return getDTIFitDataPrefix(path) is not None
 
 
-class TensorImage(fslimage.Nifti):
-    """The ``TensorImage`` class is able to load and encapsulate the diffusion
-    tensor data generated by the FSL ``dtifit`` tool.
+def looksLikeTensorImage(image):
+    """Returns ``True`` if the given :class:`.Image` looks like it could 
+    contain tensor matrix data, ``False`` otherwise.
     """
 
-    
-    def __init__(self, path):
-        """Create a ``TensorImage``.
+    return len(image.shape) == 4 and image.shape[3] == 6
 
-        :arg path: A path to a ``dtifit`` directory. Alternately, the ``path``
-                   may be a dictionary with keys
-                   ``{'v1', 'v2', 'v3', 'l1', 'l2', 'l3'}``, which specify
-                   paths to images containing the tensor eigenvectors and
-                   eigenvalues.
-        """
 
-        dtifitDir = isinstance(path, six.string_types)
+def decomposeTensorMatrix(data):
+    """Decomposes the given ``numpy`` array into six separate arrays,
+    containing the eigenvectors and eigenvalues of the tensor matrix
+    decompositions.
+
+    :arg image: A 4D ``numpy`` array with 6 volumes, which contains
+                the unique elements of diffusion tensor matrices at
+                every voxel.
 
-        if dtifitDir:
+    :returns:   
+    """
 
-            prefix = getTensorDataPrefix(path)
+    # The image contains 6 volumes, corresponding
+    # to the Dxx, Dxy, Dxz, Dyy, Dyz, Dzz elements
+    # of the tensor matrix, at each voxel.
+    #
+    # We need to re-organise this into a series of
+    # complete 3x3 tensor matrices, one for each
+    # voxel.
+
+    shape    = data.shape[:3]
+    nvoxels  = np.prod(shape)
+    matrices = np.zeros((nvoxels, 3, 3), dtype=np.float32)
+
+    # Copy the tensor matrix elements
+    # into their respective locations
+    matrices[:, 0, 0] = data[..., 0].flat
+    matrices[:, 0, 1] = data[..., 1].flat
+    matrices[:, 0, 2] = data[..., 2].flat
+    matrices[:, 1, 0] = data[..., 1].flat
+    matrices[:, 1, 1] = data[..., 3].flat
+    matrices[:, 1, 2] = data[..., 4].flat
+    matrices[:, 2, 0] = data[..., 2].flat
+    matrices[:, 2, 1] = data[..., 4].flat
+    matrices[:, 2, 2] = data[..., 5].flat
+
+    # Calculate the eigenvectors and 
+    # values on all of those matrices
+    vals, vecs = npla.eig(matrices)
+    vecShape   = list(shape) + [3]
+
+    l1 = vals[:, 0]   .reshape(shape)
+    l2 = vals[:, 1]   .reshape(shape)
+    l3 = vals[:, 2]   .reshape(shape)
+    v1 = vecs[:, 0, :].reshape(vecShape)
+    v2 = vecs[:, 1, :].reshape(vecShape)
+    v3 = vecs[:, 2, :].reshape(vecShape)
+
+    return v1, v2, v3, l1, l2, l3
+
+
+class DTIFitTensor(fslimage.Nifti):
+    """The ``DTIFitTensor`` class is able to load and encapsulate the diffusion
+    tensor data generated by the FSL ``dtifit`` tool.  The ``DtiFitTensor``
+    class supports tensor model data generated by ``dtifit``, where the
+    eigenvectors and eigenvalues of the tensor matrices have been saved as six
+    separate NIFTI images.
+    """
 
-            if prefix is None:
-                raise ValueError('Invalid path: {}'.format(path))
+    
+    def __init__(self, path):
+        """Create a ``DTIFitTensor``.
 
-            v1 = op.join(path, '{}_V1'.format(prefix))
-            v2 = op.join(path, '{}_V2'.format(prefix))
-            v3 = op.join(path, '{}_V3'.format(prefix))
-            l1 = op.join(path, '{}_L1'.format(prefix))
-            l2 = op.join(path, '{}_L2'.format(prefix))
-            l3 = op.join(path, '{}_L3'.format(prefix))
+        :arg path: A path to a ``dtifit`` directory.
+        """
 
-            paths = {'v1' : v1, 'v2' : v2, 'v3' : v3,
-                     'l1' : l1, 'l2' : l2, 'l3' : l3}
-            
-        else:
-            paths = path
+        prefix      = getDTIFitDataPrefix(path)
+        isDTIfitDir = prefix is not None
 
+        if not isDTIfitDir:
+            raise ValueError('{} does not look like a dtifit '
+                             'output directory!'.format(path))
 
-        self.__v1 = fslimage.Image(paths['v1'])
-        self.__v2 = fslimage.Image(paths['v2'])
-        self.__v3 = fslimage.Image(paths['v3'])
-        self.__l1 = fslimage.Image(paths['l1'])
-        self.__l2 = fslimage.Image(paths['l2'])
-        self.__l3 = fslimage.Image(paths['l3']) 
+        # DTIFit output directory with separate
+        # eigenvector/eigenvalue images
 
-        fslimage.Nifti.__init__(self, self.__l1.header)
+        self.__v1 = fslimage.Image(op.join(path, '{}_V1'.format(prefix)))
+        self.__v2 = fslimage.Image(op.join(path, '{}_V2'.format(prefix)))
+        self.__v3 = fslimage.Image(op.join(path, '{}_V3'.format(prefix)))
+        self.__l1 = fslimage.Image(op.join(path, '{}_L1'.format(prefix)))
+        self.__l2 = fslimage.Image(op.join(path, '{}_L2'.format(prefix)))
+        self.__l3 = fslimage.Image(op.join(path, '{}_L3'.format(prefix)))
 
-        l1dir = op.abspath(op.dirname(paths['l1']))
+        fslimage.Nifti.__init__(self, self.__l1.header)
 
-        self.dataSource = l1dir
-        self.name       = '{}[tensor]'.format(op.basename(l1dir))
+        self.dataSource = op.abspath(path)
+        self.name       = '{}[tensor]'.format(op.basename(path))
 
-        
     def V1(self): return self.__v1
     def V2(self): return self.__v2
     def V3(self): return self.__v3