Merge pull request #36 from mgxd/tst/lta

TST+FIX: LTA conversions

Author: oesteban

nitransforms/io/lta.py

@@ -3,7 +3,7 @@
 from nibabel.volumeutils import Recoder
 from nibabel.affines import voxel_sizes
 
-from .base import StringBasedStruct
+from .base import StringBasedStruct, TransformFileError
 
 
 transform_codes = Recoder((
@@ -20,10 +20,10 @@ class VolumeGeometry(StringBasedStruct):
         ('valid', 'i4'),              # Valid values: 0, 1
         ('volume', 'i4', (3, 1)),     # width, height, depth
         ('voxelsize', 'f4', (3, 1)),  # xsize, ysize, zsize
-        ('xras', 'f4', (3, 1)),       # x_r, x_a, x_s
-        ('yras', 'f4', (3, 1)),       # y_r, y_a, y_s
-        ('zras', 'f4', (3, 1)),       # z_r, z_a, z_s
-        ('cras', 'f4', (3, 1)),       # c_r, c_a, c_s
+        ('xras', 'f8', (3, 1)),       # x_r, x_a, x_s
+        ('yras', 'f8', (3, 1)),       # y_r, y_a, y_s
+        ('zras', 'f8', (3, 1)),       # z_r, z_a, z_s
+        ('cras', 'f8', (3, 1)),       # c_r, c_a, c_s
         ('filename', 'U1024')])       # Not conformant (may be >1024 bytes)
     dtype = template_dtype
 
@@ -80,7 +80,7 @@ def from_string(klass, string):
         lines = string.splitlines()
         for key in ('valid', 'filename', 'volume', 'voxelsize',
                     'xras', 'yras', 'zras', 'cras'):
-            label, valstring = lines.pop(0).split(' = ')
+            label, valstring = lines.pop(0).split(' =')
             assert label.strip() == key
 
             val = np.genfromtxt([valstring.encode()],
@@ -92,11 +92,11 @@ def from_string(klass, string):
 
 class LinearTransform(StringBasedStruct):
     template_dtype = np.dtype([
-        ('mean', 'f4', (3, 1)),       # x0, y0, z0
+        ('mean', 'f4', (3, 1)),  # x0, y0, z0
         ('sigma', 'f4'),
-        ('m_L', 'f4', (4, 4)),
-        ('m_dL', 'f4', (4, 4)),
-        ('m_last_dL', 'f4', (4, 4)),
+        ('m_L', 'f8', (4, 4)),
+        ('m_dL', 'f8', (4, 4)),
+        ('m_last_dL', 'f8', (4, 4)),
         ('src', VolumeGeometry),
         ('dst', VolumeGeometry),
         ('label', 'i4')])
@@ -190,7 +190,10 @@ def to_string(self):
     def from_string(klass, string):
         lta = klass()
         sa = lta.structarr
-        lines = string.splitlines()
+        lines = [l.strip() for l in string.splitlines()
+                 if l.strip() and not l.strip().startswith('#')]
+        if not lines or not lines[0].startswith('type'):
+            raise TransformFileError("Invalid LTA format")
         for key in ('type', 'nxforms'):
             label, valstring = lines.pop(0).split(' = ')
             assert label.strip() == key
@@ -207,12 +210,16 @@ def from_string(klass, string):
                 # Optional keys
                 if not lines[0].startswith(key):
                     continue
-                label, valstring = lines.pop(0).split(' ')
-                assert label.strip() == key
+                try:
+                    label, valstring = lines.pop(0).split(' ')
+                except ValueError:
+                    sa[key] = ''
+                else:
+                    assert label.strip() == key
 
-                val = np.genfromtxt([valstring.encode()],
-                                    dtype=klass.dtype[key])
-                sa[key] = val.reshape(sa[key].shape) if val.size else ''
+                    val = np.genfromtxt([valstring.encode()],
+                                        dtype=klass.dtype[key])
+                    sa[key] = val.reshape(sa[key].shape) if val.size else ''
 
         assert len(lta._xforms) == sa['nxforms']
         return lta
@@ -240,7 +247,7 @@ def set_type(self, target):
 
         # VOX2VOX -> RAS2RAS
         if current == 0 and target == 1:
-            M = np.linalg.inv(src.as_affine()).dot(xform['m_L']).dot(dst.as_affine())
+            M = dst.as_affine().dot(xform['m_L'].dot(np.linalg.inv(src.as_affine())))
         else:
             raise NotImplementedError(
                 "Converting {0} to {1} is not yet available".format(

nitransforms/linear.py

@@ -247,6 +247,8 @@ def load(filename, fmt='X5', reference=None):
         if lta['nxforms'] > 1:
             raise NotImplementedError("Multiple transforms are not yet supported.")
         if lta['type'] != 1:
+            # To make transforms generalize across use-cases, LTA transforms
+            # are converted to RAS-to-RAS.
             lta.set_type(1)
         matrix = lta['xforms'][0]['m_L']
     elif fmt.lower() in ('x5', 'bids'):

nitransforms/tests/data/affine-RAS.fs.v2v.lta

@@ -0,0 +1,31 @@
+# transform file affine-RAS.fs.v2v.lta
+# created by a user on Mon Oct 28 10:09:56 2019
+
+type      = 0 # LINEAR_VOX_TO_VOX
+nxforms   = 1
+mean      = 0.0000 0.0000 0.0000
+sigma     = 1.0000
+1 4 4
+9.999989867210388e-01 -9.999993490055203e-04 9.999998146668077e-04 1.454572677612305e+00
+1.404936425387859e-03 6.216088533401489e-01 -7.833265066146851e-01 3.912971496582031e+01
+1.617172238184139e-04 7.833272218704224e-01 6.216096878051758e-01 -1.376794815063477e+01
+0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00
+src volume info
+valid = 1  # volume info valid
+filename = unknown
+volume = 57 67 56
+voxelsize = 2.750000000000000e+00 2.750000000000000e+00 2.750000000000000e+00
+xras   = 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
+yras   = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00
+zras   = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00
+cras   = 3.750000000000000e-01 1.125000000000000e+00 -1.400000000000000e+01
+dst volume info
+valid = 1  # volume info valid
+filename = unknown
+volume = 57 67 56
+voxelsize = 2.750000000000000e+00 2.750000000000000e+00 2.750000000000000e+00
+xras   = 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
+yras   = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00
+zras   = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00
+cras   = 3.750000000000000e-01 1.125000000000000e+00 -1.400000000000000e+01
+fscale 0.100000

nitransforms/tests/data/inv.lta

@@ -45,7 +45,13 @@ def test_LinearTransformArray(tmpdir, data_path):
     assert lta['nxforms'] == 0
     assert len(lta['xforms']) == 0
 
-    test_lta = str(data_path / 'inv.lta')
+    # read invalid LTA file
+    test_lta = str(data_path / 'affine-RAS.fsl')
+    with pytest.raises(TransformFileError):
+        with open(test_lta) as fp:
+            LTA.from_fileobj(fp)
+
+    test_lta = str(data_path / 'affine-RAS.fs.lta')
     with open(test_lta) as fp:
         lta = LTA.from_fileobj(fp)
 
@@ -66,6 +72,24 @@ def test_LinearTransformArray(tmpdir, data_path):
     assert np.allclose(lta['xforms'][0]['m_L'], lta2['xforms'][0]['m_L'])
 
 
+def test_LT_conversions(data_path):
+    r = str(data_path / 'affine-RAS.fs.lta')
+    v = str(data_path / 'affine-RAS.fs.v2v.lta')
+    with open(r) as fa, open(v) as fb:
+        r2r = LTA.from_fileobj(fa)
+        v2v = LTA.from_fileobj(fb)
+    assert r2r['type'] == 1
+    assert v2v['type'] == 0
+
+    r2r_m = r2r['xforms'][0]['m_L']
+    v2v_m = v2v['xforms'][0]['m_L']
+    assert np.any(r2r_m != v2v_m)
+    # convert vox2vox LTA to ras2ras
+    v2v.set_type('LINEAR_RAS_TO_RAS')
+    assert v2v['type'] == 1
+    assert np.allclose(r2r_m, v2v_m, atol=1e-05)
+
+
 def test_ITKLinearTransform(tmpdir, data_path):
     tmpdir.chdir()