Merge pull request #795 from ThoumyreStanislas/correct_scil_convert_s…

…urface_vtk

adding vtk legacy options to scil_convert_surface script

requirements.txt

@@ -38,7 +38,7 @@ scipy==1.9.*
 six==1.16.*
 spams==2.6.*
 statsmodels==0.13.*
-trimeshpy==0.0.2
+trimeshpy==0.0.3
 vtk==9.2.*
 # Dipy requirements
 h5py>=2.8.0

scilpy/io/fetcher.py

@@ -94,7 +94,7 @@ def get_testing_files_dict():
              'a2f982b8d84833f5ccfe709b725307d2'],
             'surface_vtk_fib.zip':
             ['1c9KMNFeSkyYDgu3SH_aMf0kduIlpt7cN',
-             '946beb4271b905a2bd69ad2d80136ca9'],
+             'bf131869a6722778a234869bf585520a'],
             'tracking.zip':
             ['1QSekZYDoMvv-An6FRMSt_s_qPeB3BHfw',
              '6d88910403fb4d9b79604f11100d8915'],

scilpy/surfaces/tests/test_surfaces_utils.py

@@ -0,0 +1,17 @@
+import numpy as np
+
+from scilpy.surfaces.utils import (extract_xform)
+from scilpy.tests.arrays import (xform, xform_matrix_ref)
+
+
+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

@@ -0,0 +1,124 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import vtk
+from nibabel.freesurfer.io import read_geometry
+
+
+def convert_freesurfer_into_polydata(surface_to_polydata, xform):
+    """
+    Convert a freesurfer surface into a polydata surface with vtk.
+
+    Parameters
+    ----------
+    surface_to_vtk: Input a surface from freesurfer.
+        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.
+
+    Returns
+    -------
+    polydata : A polydata surface.
+        A polydata is a mesh structure that can hold data arrays
+        in points, cells, or in the dataset itself.
+    """
+    surface = read_geometry(surface_to_polydata)
+    points = vtk.vtkPoints()
+    triangles = vtk.vtkCellArray()
+
+    flip_LPS = [-1, -1, 1]
+
+    for vertex in surface[0]:
+        id = points.InsertNextPoint((vertex[0:3]+xform)*flip_LPS)
+
+    for vertex_id in surface[1]:
+        triangle = vtk.vtkTriangle()
+        triangle.GetPointIds().SetId(0, vertex_id[0])
+        triangle.GetPointIds().SetId(1, vertex_id[1])
+        triangle.GetPointIds().SetId(2, vertex_id[2])
+        triangles.InsertNextCell(triangle)
+
+    polydata = vtk.vtkPolyData()
+    polydata.SetPoints(points)
+    polydata.SetPolys(triangles)
+    polydata.Modified()
+
+    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):
+    """
+    Apply a flip to the freesurfer surface of the anteroposterior axis.
+
+    Parameters
+    ----------
+    polydata : polydata surface.
+        A surface mesh structure after a transformation in polydata
+        surface with vtk.
+
+    Returns
+    -------
+    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)
+
+    # Apply the transforms
+    transform = vtk.vtkTransform()
+    transform.Concatenate(flip_LPS)
+
+    # Apply the transforms
+    transform = vtk.vtkTransform()
+    transform.Concatenate(flip_LPS)
+
+    # Transform the polydata
+    transform_polydata = vtk.vtkTransformPolyDataFilter()
+    transform_polydata.SetTransform(transform)
+    transform_polydata.SetInputData(polydata)
+    transform_polydata.Update()
+    polydata = transform_polydata.GetOutput()
+
+    return polydata

scilpy/tests/arrays.py

@@ -252,3 +252,33 @@
 ref_out_labels = copy.deepcopy(ref_in_labels)
 for i in range(1, 8, 2):
     ref_out_labels[ref_out_labels == i] = 0
+
+# mri-info log.txt output file from freesurfer input
+xform = ['Volume information for mni__mask.nii.gz', 'type: nii',
+         'dimensions: 193 x 229 x 193',
+         'voxel sizes: 1.000000, 1.000000, 1.000000',
+         'type: SHORT (4)', 'fov: 193.000', 'dof: 1',
+         'xstart: -96.5, xend: 96.5', 'ystart: -114.5, yend: 114.5',
+         'zstart: -96.5, zend: 96.5',
+         'TR: 0.0 msec, TE: 0.0 msec, TI: 0.0 msec, flip angle: 0.0 degrees',
+         'nframes: 1', 'PhEncDir: UNKNOWN', 'FieldStrength: 0.000000',
+         'ras xform present',
+         'xform info: x_r = 1.0000, y_r = 0.0000, z_r = 0.0000, c_r = 0.5000',
+         ': x_a = 0.0000, y_a = 1.0000, z_a = 0.0000, c_a = -17.5000',
+         ': x_s = 0.0000, y_s = 0.0000, z_s = 1.0000, c_s = 18.5000',
+         'Orientation   : RAS', 'Primary Slice Direction: axial',
+         '', 'voxel to ras transform:',
+         '1.0000   0.0000   0.0000   -96.0000',
+         '0.0000   1.0000   0.0000  -132.0000',
+         '0.0000   0.0000   1.0000   -78.0000',
+         '0.0000   0.0000   0.0000     1.0000',
+         '', 'voxel-to-ras determinant 1',
+         '', 'ras to voxel transform:',
+         '1.0000   0.0000   0.0000    96.0000',
+         '0.0000   1.0000   0.0000   132.0000',
+         '0.0000   0.0000   1.0000    78.0000',
+         '0.0000   0.0000   0.0000     1.0000']
+
+# Reference matrix for xform file
+xform_matrix_ref = np.array(((1, 0, 0, 0.5), (0, 1, 0, -17.5),
+                             (0, 0, 1, 18.5), (0, 0, 0, 1)))

scripts/scil_convert_surface.py

@@ -7,10 +7,15 @@
 
 > scil_convert_surface.py surf.vtk converted_surf.ply
 """
-
 import argparse
+import os
+
+from trimeshpy.vtk_util import (load_polydata,
+                                save_polydata)
 
-from trimeshpy.io import load_mesh_from_file
+from scilpy.surfaces.utils import (convert_freesurfer_into_polydata,
+                                   extract_xform,
+                                   flip_LPS)
 
 from scilpy.io.utils import (add_overwrite_arg,
                              assert_inputs_exist,
@@ -33,19 +38,47 @@ def _build_arg_parser():
     p.add_argument('out_surface',
                    help='Output flipped 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')
+
     add_overwrite_arg(p)
+
     return p
 
 
 def main():
+
     parser = _build_arg_parser()
     args = parser.parse_args()
 
     assert_inputs_exist(parser, args.in_surface)
     assert_outputs_exist(parser, args, args.out_surface)
 
-    mesh = load_mesh_from_file(args.in_surface)
-    mesh.save(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]
+
+    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)
+    else:
+        polydata = load_polydata(args.out_surface)
+
+    if args.to_lps:
+        polydata = flip_LPS(polydata)
+
+    save_polydata(polydata, args.out_surface, legacy_vtk_format=True)
 
 
 if __name__ == "__main__":

scripts/tests/test_convert_surface.py

@@ -24,3 +24,15 @@ def test_execution_surface_vtk_fib(script_runner):
     ret = script_runner.run('scil_convert_surface.py', in_surf,
                             'rhpialt.ply')
     assert ret.success
+
+
+def test_execution_surface_vtk_xfrom(script_runner):
+    os.chdir(os.path.expanduser(tmp_dir.name))
+    in_surf = os.path.join(get_home(), 'surface_vtk_fib',
+                           'lh.pialt_xform')
+    x_form = os.path.join(get_home(), 'surface_vtk_fib',
+                          'log.txt')
+    ret = script_runner.run('scil_convert_surface.py', in_surf,
+                            'lh.pialt_xform.vtk', '--xform', x_form,
+                            '--to_lps')
+    assert ret.success