Merge pull request #1001 from frheault/xform_in_convert

Removing xform from convert script

scilpy/surfaces/tests/test_surfaces_utils.py

@@ -1,17 +1,8 @@
-import numpy as np
-
-from scilpy.surfaces.utils import (extract_xform)
-from scilpy.tests.arrays import (xform, xform_matrix_ref)
-
+# -*- coding: utf-8 -*-
 
 def test_convert_freesurfer_into_polydata():
     pass
 
 
 def test_flip_LPS():
     pass
-
-
-def test_extract_xform():
-    out = extract_xform(xform)
-    assert np.allclose(out, xform_matrix_ref)

scilpy/surfaces/utils.py

@@ -2,10 +2,11 @@
 
 import numpy as np
 import vtk
+import nibabel as nib
 from nibabel.freesurfer.io import read_geometry
 
 
-def convert_freesurfer_into_polydata(surface_to_polydata, xform):
+def convert_freesurfer_into_polydata(surface_to_polydata, reference):
     """
     Convert a freesurfer surface into a polydata surface with vtk.
 
@@ -15,9 +16,8 @@ def convert_freesurfer_into_polydata(surface_to_polydata, xform):
         The header must not contain any of these suffixes:
         '.vtk', '.vtp', '.fib', '.ply', '.stl', '.xml', '.obj'.
 
-    xform: array [float]
-        Apply a transformation matrix to the surface to align
-        freesurfer surface with T1.
+    reference: Reference image to extract the transformation matrix.
+        The reference image is used to extract the transformation matrix
 
     Returns
     -------
@@ -29,10 +29,14 @@ def convert_freesurfer_into_polydata(surface_to_polydata, xform):
     points = vtk.vtkPoints()
     triangles = vtk.vtkCellArray()
 
-    flip_LPS = [-1, -1, 1]
+    img = nib.load(reference)
+    affine = img.affine
+    center_volume = (np.array(img.shape) / 2)
+
+    xform_translation = np.dot(affine[0:3, 0:3], center_volume) + affine[0:3, 3]
 
     for vertex in surface[0]:
-        id = points.InsertNextPoint((vertex[0:3]+xform)*flip_LPS)
+        _ = points.InsertNextPoint((vertex[0:3] + xform_translation))
 
     for vertex_id in surface[1]:
         triangle = vtk.vtkTriangle()
@@ -49,44 +53,7 @@ def convert_freesurfer_into_polydata(surface_to_polydata, xform):
     return polydata
 
 
-def extract_xform(xform):
-    """
-    Use the log.txt file from mri_info to generate a transformation
-    matrix to align the freesurfer surface with the T1.
-
-    Parameters
-    ----------
-    filename : list
-        The copy-paste output from mri_info of the surface using:
-        mri_info $surface >> log.txt
-
-    Returns
-    -------
-    Matrix : np.array
-        a transformation matrix to align the surface with the T1.
-    """
-
-    raw_xform = []
-    for i in xform:
-        raw_xform.extend(i.split())
-
-    start_read = 0
-    for i, value in enumerate(raw_xform):
-        if value == 'xform':
-            start_read = int(i)
-            break
-
-    if start_read == 0:
-        raise ValueError('No xform in that file...')
-
-    matrix = np.eye(4)
-    for i in range(3):
-        for j in range(4):
-            matrix[i, j] = float(raw_xform[13*i + (j*3) + 4+2+start_read][:-1])
-    return matrix
-
-
-def flip_LPS(polydata):
+def flip_surfaces_axes(polydata, axes=[-1, -1, 1]):
     """
     Apply a flip to the freesurfer surface of the anteroposterior axis.
 
@@ -101,18 +68,15 @@ def flip_LPS(polydata):
     polydata : polydata surface.
         return the polydata turned over.
     """
-    flip_LPS = vtk.vtkMatrix4x4()
-    flip_LPS.Identity()
-    flip_LPS.SetElement(0, 0, -1)
-    flip_LPS.SetElement(1, 1, -1)
+    flip_matrix = vtk.vtkMatrix4x4()
+    flip_matrix.Identity()
 
-    # Apply the transforms
-    transform = vtk.vtkTransform()
-    transform.Concatenate(flip_LPS)
+    for i in range(3):
+        flip_matrix.SetElement(i, i, axes[i])
 
     # Apply the transforms
     transform = vtk.vtkTransform()
-    transform.Concatenate(flip_LPS)
+    transform.Concatenate(flip_matrix)
 
     # Transform the polydata
     transform_polydata = vtk.vtkTransformPolyDataFilter()

scripts/scil_surface_convert.py

@@ -20,8 +20,7 @@
                                 save_polydata)
 
 from scilpy.surfaces.utils import (convert_freesurfer_into_polydata,
-                                   extract_xform,
-                                   flip_LPS)
+                                   flip_surfaces_axes)
 
 from scilpy.io.utils import (add_overwrite_arg,
                              add_verbose_arg,
@@ -41,18 +40,16 @@ def _build_arg_parser():
 
     p.add_argument('in_surface',
                    help='Input a surface (FreeSurfer or supported by VTK).')
-
     p.add_argument('out_surface',
-                   help='Output surface (formats supported by VTK).')
-
-    p.add_argument('--xform',
-                   help='Path of the copy-paste output from mri_info \n'
-                        'Using: mri_info $input >> log.txt, \n'
-                        'The file log.txt would be this parameter')
-
-    p.add_argument('--to_lps', action='store_true',
-                   help='Flip for Surface/MI-Brain LPS')
-
+                   help='Output surface (formats supported by VTK).\n'
+                        'Recommended extension: .vtk or .ply')
+
+    p.add_argument('--flip_axes', default=[-1, -1, 1], type=int, nargs=3,
+                   help='Flip axes for RAS or LPS convention. '
+                        'Default is LPS convention (MI-Brain) %(default)s.')
+    p.add_argument('--reference',
+                   help='Reference image to extract the transformation matrix '
+                        'to align the freesurfer surface with the T1.')
     add_verbose_arg(p)
     add_overwrite_arg(p)
 
@@ -64,27 +61,24 @@ def main():
     args = parser.parse_args()
     logging.getLogger().setLevel(logging.getLevelName(args.verbose))
 
-    assert_inputs_exist(parser, args.in_surface)
+    assert_inputs_exist(parser, args.in_surface, optional=args.reference)
     assert_outputs_exist(parser, args, args.out_surface)
 
-    if args.xform:
-        with open(args.xform) as f:
-            content = f.readlines()
-        xform = [x.strip() for x in content]
-        xform_matrix = extract_xform(xform)
-        xform_translation = xform_matrix[0:3, 3]
-    else:
-        xform_translation = [0, 0, 0]
+    _, ext = os.path.splitext(args.in_surface)
+    # FreeSurfer surfaces have no extension, verify if the input has one of the
+    # many supported extensions, otherwise it is (likely) a FreeSurfer surface
+    if ext not in ['.vtk', '.vtp', '.fib', '.ply', '.stl', '.xml', '.obj']:
+        if args.reference is None:
+            parser.error('The reference image is required for FreeSurfer '
+                         'surfaces.')
 
-    if not ((os.path.splitext(args.in_surface)[1])
-            in ['.vtk', '.vtp', '.fib', '.ply', '.stl', '.xml', '.obj']):
         polydata = convert_freesurfer_into_polydata(args.in_surface,
-                                                    xform_translation)
+                                                    args.reference)
     else:
         polydata = load_polydata(args.in_surface)
 
-    if args.to_lps:
-        polydata = flip_LPS(polydata)
+    if args.flip_axes:
+        polydata = flip_surfaces_axes(polydata, args.flip_axes)
 
     save_polydata(polydata, args.out_surface, legacy_vtk_format=True)
 

scripts/tests/test_surface_convert.py

@@ -30,9 +30,8 @@ def test_execution_surface_vtk_xfrom(script_runner, monkeypatch):
     monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
     in_surf = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
                            'lh.pialt_xform')
-    x_form = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
-                          'log.txt')
+    ref = os.path.join(SCILPY_HOME, 'surface_vtk_fib', 'fa.nii.gz')
     ret = script_runner.run('scil_surface_convert.py', in_surf,
-                            'lh.pialt_xform.vtk', '--xform', x_form,
-                            '--to_lps')
+                            'lh.pialt_xform.vtk', '--reference', ref,
+                            '--flip_axes', '-1', '-1', '1')
     assert ret.success

scripts/tests/test_tractogram_assign_uniform_color.py

@@ -40,5 +40,5 @@ def test_execution_dict(script_runner, monkeypatch):
 
     ret = script_runner.run('scil_tractogram_assign_uniform_color.py',
                             in_bundle, '--dict_colors', json_file,
-                            '--out_suffix', 'colored')
+                            '--out_suffix', 'colored', '-f')
     assert ret.success