Merge pull request #39 from oesteban/enh/io-fsl-afni

ENH: Refactor of AFNI and FSL i/o with StringStructs

Author: oesteban

nitransforms/__init__.py

@@ -18,7 +18,7 @@
 """
 from .linear import Affine
 from .nonlinear import DisplacementsFieldTransform
-from .io import itk
+from .io import afni, fsl, itk
 
 
-__all__ = ['itk', 'Affine', 'DisplacementsFieldTransform']
+__all__ = ['afni', 'fsl', 'itk', 'Affine', 'DisplacementsFieldTransform']

nitransforms/base.py

@@ -225,13 +225,13 @@ def apply(self, spatialimage, reference=None,
 
         """
         if reference is not None and isinstance(reference, (str, Path)):
-            reference = load(reference)
+            reference = load(str(reference))
 
         _ref = self.reference if reference is None \
             else SpatialReference.factory(reference)
 
         if isinstance(spatialimage, (str, Path)):
-            spatialimage = load(spatialimage)
+            spatialimage = load(str(spatialimage))
 
         data = np.asanyarray(spatialimage.dataobj)
         output_dtype = output_dtype or data.dtype

nitransforms/io/afni.py

@@ -0,0 +1,123 @@
+"""Read/write AFNI's transforms."""
+from math import pi
+import numpy as np
+from nibabel.affines import obliquity, voxel_sizes
+
+from ..patched import shape_zoom_affine
+from .base import BaseLinearTransformList, LinearParameters, TransformFileError
+
+LPS = np.diag([-1, -1, 1, 1])
+OBLIQUITY_THRESHOLD_DEG = 0.01
+
+
+class AFNILinearTransform(LinearParameters):
+    """A string-based structure for AFNI linear transforms."""
+
+    def __str__(self):
+        """Generate a string representation."""
+        param = self.structarr['parameters']
+        return '\t'.join(['%g' % p for p in param[:3, :].reshape(-1)])
+
+    def to_string(self, banner=True):
+        """Convert to a string directly writeable to file."""
+        string = '%s\n' % self
+        if banner:
+            string = '\n'.join(("# 3dvolreg matrices (DICOM-to-DICOM, row-by-row):",
+                                string))
+        return string % self
+
+    @classmethod
+    def from_ras(cls, ras, moving=None, reference=None):
+        """Create an ITK affine from a nitransform's RAS+ matrix."""
+        ras = ras.copy()
+        pre = LPS.copy()
+        post = LPS.copy()
+        if _is_oblique(reference.affine):
+            print('Reference affine axes are oblique.')
+            M = reference.affine
+            A = shape_zoom_affine(reference.shape,
+                                  voxel_sizes(M), x_flip=False, y_flip=False)
+            pre = M.dot(np.linalg.inv(A)).dot(LPS)
+
+        if _is_oblique(moving.affine):
+            print('Moving affine axes are oblique.')
+            M2 = moving.affine
+            A2 = shape_zoom_affine(moving.shape,
+                                   voxel_sizes(M2), x_flip=True, y_flip=True)
+            post = A2.dot(np.linalg.inv(M2))
+
+        # swapaxes is necessary, as axis 0 encodes series of transforms
+        parameters = np.swapaxes(post.dot(ras.dot(pre)), 0, 1)
+
+        tf = cls()
+        tf.structarr['parameters'] = parameters.T
+        return tf
+
+    @classmethod
+    def from_string(cls, string):
+        """Read the struct from string."""
+        tf = cls()
+        sa = tf.structarr
+        lines = [
+            l for l in string.splitlines()
+            if l.strip() and not (l.startswith('#') or '3dvolreg matrices' in l)
+        ]
+
+        if not lines:
+            raise TransformFileError
+
+        parameters = np.vstack((
+            np.genfromtxt([lines[0].encode()],
+                          dtype='f8').reshape((3, 4)),
+            (0., 0., 0., 1.)))
+        sa['parameters'] = parameters
+        return tf
+
+
+class AFNILinearTransformArray(BaseLinearTransformList):
+    """A string-based structure for series of AFNI linear transforms."""
+
+    _inner_type = AFNILinearTransform
+
+    def to_ras(self, moving=None, reference=None):
+        """Return a nitransforms' internal RAS matrix."""
+        return np.stack([xfm.to_ras(moving=moving, reference=reference)
+                         for xfm in self.xforms])
+
+    def to_string(self):
+        """Convert to a string directly writeable to file."""
+        strings = []
+        for i, xfm in enumerate(self.xforms):
+            lines = [
+                l.strip()
+                for l in xfm.to_string(banner=(i == 0)).splitlines()
+                if l.strip()]
+            strings += lines
+        return '\n'.join(strings)
+
+    @classmethod
+    def from_ras(cls, ras, moving=None, reference=None):
+        """Create an ITK affine from a nitransform's RAS+ matrix."""
+        _self = cls()
+        _self.xforms = [cls._inner_type.from_ras(
+            ras[i, ...], moving=moving, reference=reference)
+            for i in range(ras.shape[0])]
+        return _self
+
+    @classmethod
+    def from_string(cls, string):
+        """Read the struct from string."""
+        _self = cls()
+
+        lines = [l.strip() for l in string.splitlines()
+                 if l.strip() and not l.startswith('#')]
+        if not lines:
+            raise TransformFileError('Input string is empty.')
+
+        _self.xforms = [cls._inner_type.from_string(l)
+                        for l in lines]
+        return _self
+
+
+def _is_oblique(affine, thres=OBLIQUITY_THRESHOLD_DEG):
+    return (obliquity(affine).min() * 180 / pi) > thres

nitransforms/io/base.py

@@ -1,4 +1,5 @@
 """Read/write linear transforms."""
+import numpy as np
 from scipy.io.matlab.miobase import get_matfile_version
 from scipy.io.matlab.mio4 import MatFile4Reader
 from scipy.io.matlab.mio5 import MatFile5Reader
@@ -29,12 +30,141 @@ def __array__(self):
         return self._structarr
 
 
+class LinearParameters(StringBasedStruct):
+    """
+    A string-based structure for linear transforms.
+
+    Examples
+    --------
+    >>> lp = LinearParameters()
+    >>> np.all(lp.structarr['parameters'] == np.eye(4))
+    True
+
+    >>> p = np.diag([2., 2., 2., 1.])
+    >>> lp = LinearParameters(p)
+    >>> np.all(lp.structarr['parameters'] == p)
+    True
+
+    """
+
+    template_dtype = np.dtype([
+        ('parameters', 'f8', (4, 4)),
+    ])
+    dtype = template_dtype
+
+    def __init__(self, parameters=None):
+        """Initialize with default parameters."""
+        super().__init__()
+        self.structarr['parameters'] = np.eye(4)
+        if parameters is not None:
+            self.structarr['parameters'] = parameters
+
+    def to_filename(self, filename):
+        """Store this transform to a file with the appropriate format."""
+        with open(str(filename), 'w') as f:
+            f.write(self.to_string())
+
+    def to_ras(self, moving=None, reference=None):
+        """Return a nitransforms internal RAS+ matrix."""
+        raise NotImplementedError
+
+    @classmethod
+    def from_filename(cls, filename):
+        """Read the struct from a file given its path."""
+        with open(str(filename)) as f:
+            string = f.read()
+        return cls.from_string(string)
+
+    @classmethod
+    def from_fileobj(cls, fileobj, check=True):
+        """Read the struct from a file object."""
+        return cls.from_string(fileobj.read())
+
+    @classmethod
+    def from_string(cls, string):
+        """Read the struct from string."""
+        raise NotImplementedError
+
+
+class BaseLinearTransformList(StringBasedStruct):
+    """A string-based structure for series of linear transforms."""
+
+    template_dtype = np.dtype([('nxforms', 'i4')])
+    dtype = template_dtype
+    _xforms = None
+    _inner_type = LinearParameters
+
+    def __init__(self,
+                 xforms=None,
+                 binaryblock=None,
+                 endianness=None,
+                 check=True):
+        """Initialize with (optionally) a list of transforms."""
+        super().__init__(binaryblock, endianness, check)
+        self.xforms = [self._inner_type(parameters=mat)
+                       for mat in xforms or []]
+
+    @property
+    def xforms(self):
+        """Get the list of internal transforms."""
+        return self._xforms
+
+    @xforms.setter
+    def xforms(self, value):
+        self._xforms = list(value)
+
+    def __getitem__(self, idx):
+        """Allow dictionary access to the transforms."""
+        if idx == 'xforms':
+            return self._xforms
+        if idx == 'nxforms':
+            return len(self._xforms)
+        raise KeyError(idx)
+
+    def to_filename(self, filename):
+        """Store this transform to a file with the appropriate format."""
+        with open(str(filename), 'w') as f:
+            f.write(self.to_string())
+
+    def to_ras(self, moving=None, reference=None):
+        """Return a nitransforms' internal RAS matrix."""
+        raise NotImplementedError
+
+    def to_string(self):
+        """Convert to a string directly writeable to file."""
+        raise NotImplementedError
+
+    @classmethod
+    def from_filename(cls, filename):
+        """Read the struct from a file given its path."""
+        with open(str(filename)) as f:
+            string = f.read()
+        return cls.from_string(string)
+
+    @classmethod
+    def from_fileobj(cls, fileobj, check=True):
+        """Read the struct from a file object."""
+        return cls.from_string(fileobj.read())
+
+    @classmethod
+    def from_ras(cls, ras, moving=None, reference=None):
+        """Create an ITK affine from a nitransform's RAS+ matrix."""
+        raise NotImplementedError
+
+    @classmethod
+    def from_string(cls, string):
+        """Read the struct from string."""
+        raise NotImplementedError
+
+
 def _read_mat(byte_stream):
     mjv, _ = get_matfile_version(byte_stream)
     if mjv == 0:
         reader = MatFile4Reader(byte_stream)
     elif mjv == 1:
         reader = MatFile5Reader(byte_stream)
     elif mjv == 2:
-        raise TransformFileError('Please use HDF reader for matlab v7.3 files')
+        raise TransformFileError('Please use HDF reader for Matlab v7.3 files')
+    else:
+        raise TransformFileError('Not a Matlab file.')
     return reader.get_variables()