Merge pull request #31 from oesteban/enh/19

ENH: Rewrite load/save utilities for ITK's MatrixOffsetBased transforms in ``io``

Author: oesteban

nitransforms/__init__.py

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

nitransforms/io/__init__.py

@@ -0,0 +1,22 @@
+# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+#
+#   See COPYING file distributed along with the NiBabel package for the
+#   copyright and license terms.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+"""
+Read and write transforms.
+
+.. currentmodule:: nitransforms.io
+
+.. autosummary::
+   :toctree: ../generated
+
+   io
+"""
+from .lta import LinearTransform, LinearTransformArray, VolumeGeometry
+
+
+__all__ = ['LinearTransform', 'LinearTransformArray', 'VolumeGeometry']

nitransforms/io/base.py

@@ -0,0 +1,23 @@
+"""Read/write linear transforms."""
+import numpy as np
+from nibabel.wrapstruct import LabeledWrapStruct as LWS
+
+
+class LabeledWrapStruct(LWS):
+    def __setitem__(self, item, value):
+        self._structarr[item] = np.asanyarray(value)
+
+
+class StringBasedStruct(LabeledWrapStruct):
+    def __init__(self,
+                 binaryblock=None,
+                 endianness=None,
+                 check=True):
+        if binaryblock is not None and getattr(binaryblock, 'dtype',
+                                               None) == self.dtype:
+            self._structarr = binaryblock.copy()
+            return
+        super(StringBasedStruct, self).__init__(binaryblock, endianness, check)
+
+    def __array__(self):
+        return self._structarr

nitransforms/io/itk.py

@@ -0,0 +1,134 @@
+"""Read/write ITK transforms."""
+import numpy as np
+from .base import StringBasedStruct
+
+
+class ITKLinearTransform(StringBasedStruct):
+    """A string-based structure for ITK linear transforms."""
+
+    template_dtype = np.dtype([
+        ('type', 'i4'),
+        ('index', 'i4'),
+        ('parameters', 'f4', (4, 4)),
+        ('offset', 'f4', 3),  # Center of rotation
+    ])
+    dtype = template_dtype
+
+    def __init__(self):
+        """Initialize with default offset and index."""
+        super().__init__()
+        self.structarr['offset'] = [0, 0, 0]
+        self.structarr['index'] = 1
+        self.structarr['parameters'] = np.eye(4)
+
+    def __str__(self):
+        """Generate a string representation."""
+        sa = self.structarr
+        lines = [
+            '#Transform {:d}'.format(sa['index']),
+            'Transform: MatrixOffsetTransformBase_double_3_3',
+            'Parameters: {}'.format(' '.join(
+                ['%g' % p
+                 for p in sa['parameters'][:3, :3].reshape(-1).tolist() +
+                 sa['parameters'][:3, 3].tolist()])),
+            'FixedParameters: {:g} {:g} {:g}'.format(*sa['offset']),
+            '',
+        ]
+        return '\n'.join(lines)
+
+    def to_string(self, banner=None):
+        """Convert to a string directly writeable to file."""
+        string = '%s'
+
+        if banner is None:
+            banner = self.structarr['index'] == 1
+
+        if banner:
+            string = '#Insight Transform File V1.0\n%s'
+        return string % self
+
+    @classmethod
+    def from_string(klass, string):
+        """Read the struct from string."""
+        tf = klass()
+        sa = tf.structarr
+        lines = [l for l in string.splitlines()
+                 if l.strip()]
+        assert lines[0][0] == '#'
+        if lines[1][0] == '#':
+            lines = lines[1:]  # Drop banner with version
+
+        parameters = np.eye(4, dtype='f4')
+        sa['index'] = int(lines[0][lines[0].index('T'):].split()[1])
+        sa['offset'] = np.genfromtxt([lines[3].split(':')[-1].encode()],
+                                     dtype=klass.dtype['offset'])
+        vals = np.genfromtxt([lines[2].split(':')[-1].encode()],
+                             dtype='f4')
+        parameters[:3, :3] = vals[:-3].reshape((3, 3))
+        parameters[:3, 3] = vals[-3:]
+        sa['parameters'] = parameters
+        return tf
+
+    @classmethod
+    def from_fileobj(klass, fileobj, check=True):
+        """Read the struct from a file object."""
+        return klass.from_string(fileobj.read())
+
+
+class ITKLinearTransformArray(StringBasedStruct):
+    """A string-based structure for series of ITK linear transforms."""
+
+    template_dtype = np.dtype([('nxforms', 'i4')])
+    dtype = template_dtype
+    _xforms = None
+
+    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 = []
+        for i, mat in enumerate(xforms or []):
+            xfm = ITKLinearTransform()
+            xfm['parameters'] = mat
+            xfm['index'] = i + 1
+            self._xforms.append(xfm)
+
+    def __getitem__(self, idx):
+        """Allow dictionary access to the transforms."""
+        if idx == 'xforms':
+            return self._xforms
+        if idx == 'nxforms':
+            return len(self._xforms)
+        return super().__getitem__(idx)
+
+    def to_string(self):
+        """Convert to a string directly writeable to file."""
+        strings = []
+        for i, xfm in enumerate(self._xforms):
+            xfm.structarr['index'] = i + 1
+            strings.append(xfm.to_string())
+        return '\n'.join(strings)
+
+    @classmethod
+    def from_string(klass, string):
+        """Read the struct from string."""
+        _self = klass()
+        lines = [l.strip() for l in string.splitlines()
+                 if l.strip()]
+
+        if lines[0][0] != '#' or 'Insight Transform File V1.0' not in lines[0]:
+            raise ValueError('Unknown Insight Transform File format.')
+
+        string = '\n'.join(lines[1:])
+        for xfm in string.split('#')[1:]:
+            _self._xforms.append(ITKLinearTransform.from_string(
+                '#%s' % xfm))
+        return _self
+
+    @classmethod
+    def from_fileobj(klass, fileobj, check=True):
+        """Read the struct from a file object."""
+        return klass.from_string(fileobj.read())

nitransforms/io.py renamed to nitransforms/io/lta.py

@@ -1,8 +1,10 @@
+"""Read/write linear transforms."""
 import numpy as np
 from nibabel.volumeutils import Recoder
 from nibabel.affines import voxel_sizes
 
-from .patched import LabeledWrapStruct
+from .base import StringBasedStruct
+
 
 transform_codes = Recoder((
     (0, 'LINEAR_VOX_TO_VOX'),
@@ -13,21 +15,6 @@
     fields=('code', 'label'))
 
 
-class StringBasedStruct(LabeledWrapStruct):
-    def __init__(self,
-                 binaryblock=None,
-                 endianness=None,
-                 check=True):
-        if binaryblock is not None and getattr(binaryblock, 'dtype',
-                                               None) == self.dtype:
-            self._structarr = binaryblock.copy()
-            return
-        super(StringBasedStruct, self).__init__(binaryblock, endianness, check)
-
-    def __array__(self):
-        return self._structarr
-
-
 class VolumeGeometry(StringBasedStruct):
     template_dtype = np.dtype([
         ('valid', 'i4'),              # Valid values: 0, 1

nitransforms/linear.py

@@ -133,19 +133,10 @@ def _to_hdf5(self, x5_root):
     def to_filename(self, filename, fmt='X5', moving=None):
         """Store the transform in BIDS-Transforms HDF5 file format (.x5)."""
         if fmt.lower() in ['itk', 'ants', 'elastix']:
+            itkobj = io.itk.ITKLinearTransformArray(
+                xforms=[LPS.dot(m.dot(LPS)) for m in self.matrix])
             with open(filename, 'w') as f:
-                f.write('#Insight Transform File V1.0\n')
-
-                for i in range(self.matrix.shape[0]):
-                    parameters = LPS.dot(self.matrix[i].dot(LPS))
-                    parameters = parameters[:3, :3].reshape(-1).tolist() + \
-                        parameters[:3, 3].tolist()
-                    itkfmt = """\
-#Transform {0}
-Transform: MatrixOffsetTransformBase_double_3_3
-Parameters: {1}
-FixedParameters: 0 0 0\n""".format
-                    f.write(itkfmt(i, ' '.join(['%g' % p for p in parameters])))
+                f.write(itkobj.to_string())
             return filename
 
         if fmt.lower() == 'afni':
@@ -239,23 +230,16 @@ def load(filename, fmt='X5', reference=None):
     """Load a linear transform."""
     if fmt.lower() in ['itk', 'ants', 'elastix', 'nifty']:
         with open(filename) as itkfile:
-            itkxfm = itkfile.read().splitlines()
-
-        if 'Insight Transform File' not in itkxfm[0]:
-            raise ValueError(
-                "File %s does not look like an ITK transform text file." % filename)
+            itkxfm = io.itk.ITKLinearTransformArray.from_fileobj(
+                itkfile)
 
         matlist = []
-        for nxfm in range(len(itkxfm[1:]) // 4):
-            parameters = np.fromstring(itkxfm[nxfm * 4 + 3].split(':')[-1].strip(),
-                                       dtype=float, sep=' ')
-            offset = np.fromstring(itkxfm[nxfm * 4 + 4].split(':')[-1].strip(),
-                                   dtype=float, sep=' ')
-            if len(parameters) == 12:
-                matrix = from_matvec(parameters[:9].reshape((3, 3)), parameters[9:])
-                c_neg = from_matvec(np.eye(3), offset * -1.0)
-                c_pos = from_matvec(np.eye(3), offset)
-                matlist.append(LPS.dot(c_pos.dot(matrix.dot(c_neg.dot(LPS)))))
+        for xfm in itkxfm['xforms']:
+            matrix = xfm['parameters']
+            offset = xfm['offset']
+            c_neg = from_matvec(np.eye(3), offset * -1.0)
+            c_pos = from_matvec(np.eye(3), offset)
+            matlist.append(LPS.dot(c_pos.dot(matrix.dot(c_neg.dot(LPS)))))
         matrix = np.stack(matlist)
     # elif fmt.lower() == 'afni':
     #     parameters = LPS.dot(self.matrix.dot(LPS))

nitransforms/patched.py

@@ -1,5 +1,4 @@
 import numpy as np
-from nibabel.wrapstruct import LabeledWrapStruct as LWS
 
 
 def shape_zoom_affine(shape, zooms, x_flip=True, y_flip=False):
@@ -64,8 +63,3 @@ def shape_zoom_affine(shape, zooms, x_flip=True, y_flip=False):
     aff[:3, :3] = np.diag(zooms)
     aff[:3, -1] = -origin * zooms
     return aff
-
-
-class LabeledWrapStruct(LWS):
-    def __setitem__(self, item, value):
-        self._structarr[item] = np.asanyarray(value)

nitransforms/tests/data/itktflist.tfm

@@ -0,0 +1,45 @@
+#Insight Transform File V1.0
+#Transform 1
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: 1 0 0 0 1 0 0 0 1 0 0 0
+FixedParameters: 10 10 10
+
+#Transform 2
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: 1 0 0 0 1 0 0 0 1 0 0 0
+FixedParameters: 10 10 10
+
+#Transform 3
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: 1 0 0 0 1 0 0 0 1 0 0 0
+FixedParameters: 10 10 10
+
+#Transform 4
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: 1 0 0 0 1 0 0 0 1 0 0 0
+FixedParameters: 10 10 10
+
+#Transform 5
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: -1.53626 0.71973 -0.639856 -0.190759 -1.80082 -0.915885 0.502537 1.12532 0.275748 0.393413 1.13855 0.761131
+FixedParameters: -0.0993171 0.364984 1.99264
+
+#Transform 6
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: -0.130507 -1.03017 2.08189 -1.51723 1.37849 -0.0890962 -0.656323 0.242694 2.15801 -1.26689 0.367131 1.23616
+FixedParameters: 0.626607 0.15351 1.24982
+
+#Transform 7
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: -1.55395 -0.36383 -0.17749 1.3387 -0.384534 -0.901462 -1.06598 -0.448228 -1.07535 1.92599 0.454696 0.576697
+FixedParameters: -0.425602 0.333406 -1.14957
+
+#Transform 8
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: 0.723719 -1.05617 -0.800562 -2.47048 -1.76301 -1.4447 -0.749896 1.29774 -1.48893 1.02789 0.65017 -1.48326
+FixedParameters: 0.800882 -1.20202 1.25495
+
+#Transform 9
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: 1.24025 -0.77628 0.618013 -0.523829 1.09471 1.66921 0.73753 -1.33588 -0.627659 -0.449913 -0.00124181 0.21433
+FixedParameters: -0.226504 -0.877893 0.2608

nitransforms/tests/test_io.py

@@ -4,12 +4,12 @@
 import numpy as np
 
 from ..io import (
+    itk,
     VolumeGeometry as VG,
     LinearTransform as LT,
     LinearTransformArray as LTA,
 )
 
-
 def test_VolumeGeometry(tmpdir, get_testdata):
     vg = VG()
     assert vg['valid'] == 0
@@ -55,3 +55,32 @@ def test_LinearTransformArray(tmpdir, data_path):
     with open(outlta) as fp:
         lta2 = LTA.from_fileobj(fp)
     assert np.allclose(lta['xforms'][0]['m_L'], lta2['xforms'][0]['m_L'])
+
+
+def test_ITKLinearTransformArray(tmpdir, data_path):
+    tmpdir.chdir()
+
+    with open(str(data_path / 'itktflist.tfm')) as f:
+        text = f.read()
+        f.seek(0)
+        itklist = itk.ITKLinearTransformArray.from_fileobj(f)
+
+    assert itklist['nxforms'] == 9
+    assert text == itklist.to_string()
+
+    itklist = itk.ITKLinearTransformArray(
+        xforms=[np.around(np.random.normal(size=(4, 4)), decimals=5)
+                for _ in range(4)])
+
+    assert itklist['nxforms'] == 4
+    assert itklist['xforms'][1].structarr['index'] == 2
+
+    xfm = itklist['xforms'][1]
+    xfm['index'] = 1
+    with open('extracted.tfm', 'w') as f:
+        f.write(xfm.to_string())
+
+    with open('extracted.tfm') as f:
+        xfm2 = itk.ITKLinearTransform.from_fileobj(f)
+    assert np.allclose(xfm.structarr['parameters'][:3, ...],
+                       xfm2.structarr['parameters'][:3, ...])