rearranging functions and adding a test

scilus · arnaudbore · Nov 20, 2023 · Oct 25, 2023 · Oct 27, 2023 · Oct 27, 2023
commit 725ecb521302153083d7e933f0ebf46cfb115ad1
diff --git a/scilpy/io/fetcher.py b/scilpy/io/fetcher.py
@@ -93,7  93,7 @@ def get_testing_files_dict():
              'a2f982b8d84833f5ccfe709b725307d2'],
             'surface_vtk_fib.zip':
             ['1c9KMNFeSkyYDgu3SH_aMf0kduIlpt7cN',
-             '946beb4271b905a2bd69ad2d80136ca9'],
              'bf131869a6722778a234869bf585520a'],
             'tracking.zip':
             ['1QSekZYDoMvv-An6FRMSt_s_qPeB3BHfw',
              '6d88910403fb4d9b79604f11100d8915'],

diff --git a/scilpy/surfaces/utils.py b/scilpy/surfaces/utils.py
@@ -0,0  1,117 @@
 # -*- 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: matrix
         Apply a transformation matrix to the surface to align 
         freesurfer surface with T1.
     Returns
     -------
     polydata: A polydata is a surface 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(filename):
     """
     Use the log.txt file from mri_info to generate a transformation
     matrix to align the freesurfer surface with the T1.
     Parameters
     ----------
     filename: .txt file.
         The copy-paste output from mri_info of the surface using:
           mri_info $surface >> log.txt
     Returns
     -------
     Matrix: a transformation matrix to align the surface with the T1.
     """  
     with open(filename) as f:
         content = f.readlines()
     names = [x.strip() for x in content]
 
     raw_xform = []
     for i in names:
         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: The polydata surface 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
diff --git a/scilpy/utils/util.py b/scilpy/utils/util.py
@@ -3,7  3,6 @@
 import collections.abc
 
 from dipy.io.utils import get_reference_info
-import vtk
 from dipy.segment.mask import bounding_box
 import numpy as np
 from numpy.lib.index_tricks import r_ as row
@@ -185,58  184,3 @@ def recursive_print(data):
         recursive_print(data[list(data.keys())[0]])
     else:
         return
-
-
-def extract_xform(filename):
-    """Use the log.txt file of the scil_convert_surface.py script
-    to align the freesurfer surface with T1"""
-
-    with open(filename) as f:
-        content = f.readlines()
-    names = [x.strip() for x in content]
-
-    raw_xform = []
-    for i in names:
-        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
-    in the scil_convert_surface.py script"""
-
-    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
diff --git a/scripts/scil_convert_surface.py b/scripts/scil_convert_surface.py
@@ -9,14  9,13 @@
 """
 import argparse
 import os
-import vtk
-from nibabel.freesurfer.io import read_geometry
 
 from trimeshpy.vtk_util import (load_polydata,
                                 save_polydata)
 
-from scilpy.utils.util import (flip_LPS,
-                               extract_xform)
 from scilpy.surfaces.utils import (convert_freesurfer_into_polydata,
                                    extract_xform,
                                    flip_LPS)
 
 from scilpy.io.utils import (add_overwrite_arg,
                              assert_inputs_exist,
@@ -52,32  51,6 @@ def _build_arg_parser():
     return p
 
 
-def convert_with_vtk_legacy(surface_to_vtk, xform):
-
-    surface = read_geometry(surface_to_vtk)
-    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 main():
 
     parser = _build_arg_parser()
@@ -94,8  67,9 @@ def main():
 
     if not ((os.path.splitext(args.in_surface)[1])
             in ['.vtk', '.vtp', '.fib', '.ply', '.stl', '.xml', '.obj']):
-        polydata = convert_with_vtk_legacy(args.in_surface, xform_translation)
-
         polydata = convert_freesurfer_into_polydata(args.in_surface, 
                                                     xform_translation)
 
     else:
         polydata = load_polydata(args.out_surface)
 

diff --git a/scripts/tests/test_convert_surface.py b/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