diff --git a/scilpy/surfaces/surface_operations.py b/scilpy/surfaces/surface_operations.py new file mode 100644 index 000000000..32a59b95c --- /dev/null +++ b/scilpy/surfaces/surface_operations.py @@ -0,0 +1,82 @@ +# -*- coding: utf-8 -*- + +import numpy as np +from scipy.ndimage import map_coordinates +import trimeshpy.vtk_util as vtk_u + + +def apply_transform(mesh, ants_affine, ants_warp=None): + """ + Apply transformation to a surface + - Apply linear transformation with affine + - Apply non linear transformation with ants_warp + + Parameters + ---------- + mesh: trimeshpy - Triangle Mesh VTK class + Moving surface + ants_affine: numpy.ndarray + Transformation matrix to be applied + ants_warp: nib.Nifti1Image + Warp image from ANTs + + Returns + ------- + mesh: trimeshpy - Triangle Mesh VTK class + Surface moved + """ + # Affine transformation + inv_affine = np.linalg.inv(ants_affine) + + # Transform mesh vertices + mesh.set_vertices(mesh.vertices_affine(inv_affine)) + + # Flip triangle face, if needed + if mesh.is_transformation_flip(inv_affine): + mesh.set_triangles(mesh.triangles_face_flip()) + + if ants_warp is not None: + warp_img = np.squeeze(ants_warp.get_fdata(dtype=np.float32)) + + # Get vertices translation in voxel space, from the warp image + vts_vox = vtk_u.vtk_to_vox(mesh.get_vertices(), warp_img) + tx = map_coordinates(warp_img[..., 0], vts_vox.T, order=1) + ty = map_coordinates(warp_img[..., 1], vts_vox.T, order=1) + tz = map_coordinates(warp_img[..., 2], vts_vox.T, order=1) + + # Apply vertices translation in world coordinates + mesh.set_vertices(mesh.get_vertices() + np.array([tx, ty, tz]).T) + + return mesh + + +def flip(mesh, axes): + """ + Apply flip to a surface + + Parameters + ---------- + mesh: trimeshpy - Triangle Mesh VTK class + Moving surface + axes: list + Axes (or normal orientation) you want to flip + + Returns + ------- + mesh: trimeshpy - Triangle Mesh VTK class + Surface flipped + """ + # Flip axes + flip = (-1 if 'x' in axes else 1, + -1 if 'y' in axes else 1, + -1 if 'z' in axes else 1) + tris, vts = mesh.flip_triangle_and_vertices(flip) + mesh.set_vertices(vts) + mesh.set_triangles(tris) + + # Reverse surface orientation + if 'n' in axes: + tris = mesh.triangles_face_flip() + mesh.set_triangles(tris) + + return mesh diff --git a/scilpy/surfaces/tests/test_surface_operations.py b/scilpy/surfaces/tests/test_surface_operations.py new file mode 100644 index 000000000..41b2aa74e --- /dev/null +++ b/scilpy/surfaces/tests/test_surface_operations.py @@ -0,0 +1,10 @@ +# -*- coding: utf-8 -*- + +def test_surface_apply_transfo(): + # toDO + pass + + +def test_surface_flip(): + # toDo + pass diff --git a/scripts/legacy/scil_apply_transform_to_surface.py b/scripts/legacy/scil_apply_transform_to_surface.py new file mode 100755 index 000000000..684466a88 --- /dev/null +++ b/scripts/legacy/scil_apply_transform_to_surface.py @@ -0,0 +1,22 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from scilpy.io.deprecator import deprecate_script +from scripts.scil_surface_apply_transform import main as new_main + + +DEPRECATION_MSG = """ +This script has been renamed scil_surface_apply_transform.py. Please change +your existing pipelines accordingly. + +""" + + +@deprecate_script("scil_apply_transform_to_surface.py", + DEPRECATION_MSG, '1.7.0') +def main(): + new_main() + + +if __name__ == "__main__": + main() diff --git a/scripts/legacy/scil_convert_surface.py b/scripts/legacy/scil_convert_surface.py new file mode 100755 index 000000000..3842baaf0 --- /dev/null +++ b/scripts/legacy/scil_convert_surface.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from scilpy.io.deprecator import deprecate_script +from scripts.scil_surface_convert import main as new_main + + +DEPRECATION_MSG = """ +This script has been renamed scil_surface_convert.py. +Please change your existing pipelines accordingly. +""" + + +@deprecate_script("scil_convert_surface.py", + DEPRECATION_MSG, '1.7.0') +def main(): + new_main() + + +if __name__ == "__main__": + main() diff --git a/scripts/legacy/scil_flip_surface.py b/scripts/legacy/scil_flip_surface.py new file mode 100755 index 000000000..9304335db --- /dev/null +++ b/scripts/legacy/scil_flip_surface.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from scilpy.io.deprecator import deprecate_script +from scripts.scil_surface_flip import main as new_main + + +DEPRECATION_MSG = """ +This script has been renamed scil_surface_flip.py. +Please change your existing pipelines accordingly. +""" + + +@deprecate_script("scil_flip_surface.py", + DEPRECATION_MSG, '1.7.0') +def main(): + new_main() + + +if __name__ == "__main__": + main() diff --git a/scripts/legacy/scil_smooth_surface.py b/scripts/legacy/scil_smooth_surface.py new file mode 100755 index 000000000..2030e083a --- /dev/null +++ b/scripts/legacy/scil_smooth_surface.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from scilpy.io.deprecator import deprecate_script +from scripts.scil_surface_smooth import main as new_main + + +DEPRECATION_MSG = """ +This script has been renamed scil_surface_smooth.py. +Please change your existing pipelines accordingly. +""" + + +@deprecate_script("scil_smooth_surface.py", + DEPRECATION_MSG, '1.7.0') +def main(): + new_main() + + +if __name__ == "__main__": + main() diff --git a/scripts/scil_apply_transform_to_surface.py b/scripts/scil_surface_apply_transform.py similarity index 64% rename from scripts/scil_apply_transform_to_surface.py rename to scripts/scil_surface_apply_transform.py index 1a3b6ede1..6a635d4ac 100755 --- a/scripts/scil_apply_transform_to_surface.py +++ b/scripts/scil_surface_apply_transform.py @@ -7,11 +7,12 @@ Best usage with ANTs from T1 to b0: > ConvertTransformFile 3 output0GenericAffine.mat vtk_transfo.txt --hm -> scil_transform_surface.py lh_white_lps.vtk affine.txt lh_white_b0.vtk\\ +> scil_surface_apply_transform.py lh_white_lps.vtk affine.txt lh_white_b0.vtk\\ --ants_warp warp.nii.gz The input surface needs to be in *T1 world LPS* coordinates (aligned over the T1 in MI-Brain). +The script will use the linear affine first and then the warp image from ANTs. The resulting surface should be aligned *b0 world LPS* coordinates (aligned over the b0 in MI-Brain). """ @@ -20,13 +21,12 @@ import nibabel as nib import numpy as np -from scipy.ndimage import map_coordinates from trimeshpy.io import load_mesh_from_file -import trimeshpy.vtk_util as vtk_u from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist, assert_outputs_exist) +from scilpy.surfaces.surface_operations import apply_transform EPILOG = """ @@ -67,32 +67,16 @@ def main(): # Load mesh mesh = load_mesh_from_file(args.in_surface) - # Affine transformation - if args.ants_affine: - # Load affine - affine = np.loadtxt(args.ants_affine) - inv_affine = np.linalg.inv(affine) - - # Transform mesh vertices - mesh.set_vertices(mesh.vertices_affine(inv_affine)) - - # Flip triangle face, if needed - if mesh.is_transformation_flip(inv_affine): - mesh.set_triangles(mesh.triangles_face_flip()) + warp_img = None + # Load affine + affine = np.loadtxt(args.ants_affine) + # Non linear transformation if args.ants_warp: # Load warp warp_img = nib.load(args.ants_warp) - warp = np.squeeze(warp_img.get_fdata(dtype=np.float32)) - - # Get vertices translation in voxel space, from the warp image - vts_vox = vtk_u.vtk_to_vox(mesh.get_vertices(), warp_img) - tx = map_coordinates(warp[..., 0], vts_vox.T, order=1) - ty = map_coordinates(warp[..., 1], vts_vox.T, order=1) - tz = map_coordinates(warp[..., 2], vts_vox.T, order=1) - # Apply vertices translation in world coordinates - mesh.set_vertices(mesh.get_vertices() + np.array([tx, ty, tz]).T) + mesh = apply_transform(mesh, affine, warp_img) # Save mesh mesh.save(args.out_surface) diff --git a/scripts/scil_convert_surface.py b/scripts/scil_surface_convert.py similarity index 97% rename from scripts/scil_convert_surface.py rename to scripts/scil_surface_convert.py index 1a99b3fba..625f4edc2 100755 --- a/scripts/scil_convert_surface.py +++ b/scripts/scil_surface_convert.py @@ -5,7 +5,7 @@ Script to convert a surface (FreeSurfer or VTK supported). ".vtk", ".vtp", ".ply", ".stl", ".xml", ".obj" -> scil_convert_surface.py surf.vtk converted_surf.ply +> scil_surface_convert.py surf.vtk converted_surf.ply """ import argparse import os diff --git a/scripts/scil_flip_surface.py b/scripts/scil_surface_flip.py similarity index 79% rename from scripts/scil_flip_surface.py rename to scripts/scil_surface_flip.py index fc8374c5f..793929e28 100755 --- a/scripts/scil_flip_surface.py +++ b/scripts/scil_surface_flip.py @@ -10,7 +10,7 @@ Best usage for FreeSurfer to LPS vtk (for MI-Brain): !!! important FreeSurfer surfaces must be in their respective folder !!! > mris_convert --to-scanner lh.white lh.white.vtk -> scil_flip_surface.py lh.white.vtk lh_white_lps.vtk x y +> scil_surface_flip.py lh.white.vtk lh_white_lps.vtk x y """ import argparse @@ -20,6 +20,7 @@ from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist, assert_outputs_exist) +from scilpy.surfaces.surface_operations import flip EPILOG = """ References: @@ -57,19 +58,7 @@ def main(): # Load mesh mesh = load_mesh_from_file(args.in_surface) - # Flip axes - flip = (-1 if 'x' in args.axes else 1, - -1 if 'y' in args.axes else 1, - -1 if 'z' in args.axes else 1) - tris, vts = mesh.flip_triangle_and_vertices(flip) - mesh.set_vertices(vts) - mesh.set_triangles(tris) - - # Reverse surface orientation - if 'n' in args.axes: - tris = mesh.triangles_face_flip() - mesh.set_triangles(tris) - + mesh = flip(mesh, args.axes) # Save mesh.save(args.out_surface) diff --git a/scripts/scil_smooth_surface.py b/scripts/scil_surface_smooth.py similarity index 100% rename from scripts/scil_smooth_surface.py rename to scripts/scil_surface_smooth.py diff --git a/scripts/tests/test_apply_transform_to_surface.py b/scripts/tests/test_surface_apply_transform.py similarity index 83% rename from scripts/tests/test_apply_transform_to_surface.py rename to scripts/tests/test_surface_apply_transform.py index 9d45dff51..4ebc278de 100644 --- a/scripts/tests/test_apply_transform_to_surface.py +++ b/scripts/tests/test_surface_apply_transform.py @@ -13,7 +13,7 @@ def test_help_option(script_runner): - ret = script_runner.run('scil_apply_transform_to_surface.py', '--help') + ret = script_runner.run('scil_surface_apply_transform.py', '--help') assert ret.success @@ -23,6 +23,6 @@ def test_execution_surface_vtk_fib(script_runner): 'lhpialt.vtk') in_aff = os.path.join(get_home(), 'surface_vtk_fib', 'affine.txt') - ret = script_runner.run('scil_apply_transform_to_surface.py', in_surf, + ret = script_runner.run('scil_surface_apply_transform.py', in_surf, in_aff, 'lhpialt_lin.vtk') assert ret.success diff --git a/scripts/tests/test_flip_surface.py b/scripts/tests/test_surface_flip.py similarity index 85% rename from scripts/tests/test_flip_surface.py rename to scripts/tests/test_surface_flip.py index e780a37d4..44efd7137 100644 --- a/scripts/tests/test_flip_surface.py +++ b/scripts/tests/test_surface_flip.py @@ -13,7 +13,7 @@ def test_help_option(script_runner): - ret = script_runner.run('scil_flip_surface.py', '--help') + ret = script_runner.run('scil_surface_flip.py', '--help') assert ret.success @@ -23,6 +23,6 @@ def test_execution_surface_vtk_fib(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_surf = os.path.join(get_home(), 'surface_vtk_fib', 'lhpialt.vtk') - ret = script_runner.run('scil_flip_surface.py', in_surf, 'rhpialt.vtk', + ret = script_runner.run('scil_surface_flip.py', in_surf, 'rhpialt.vtk', 'x') assert ret.success diff --git a/scripts/tests/test_smooth_surface.py b/scripts/tests/test_surface_smooth.py similarity index 84% rename from scripts/tests/test_smooth_surface.py rename to scripts/tests/test_surface_smooth.py index 9a31cc161..4b55e1c2a 100644 --- a/scripts/tests/test_smooth_surface.py +++ b/scripts/tests/test_surface_smooth.py @@ -13,7 +13,7 @@ def test_help_option(script_runner): - ret = script_runner.run('scil_smooth_surface.py', '--help') + ret = script_runner.run('scil_surface_smooth.py', '--help') assert ret.success @@ -21,6 +21,6 @@ def test_execution_surface_vtk_fib(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_surf = os.path.join(get_home(), 'surface_vtk_fib', 'lhpialt.vtk') - ret = script_runner.run('scil_smooth_surface.py', in_surf, + ret = script_runner.run('scil_surface_smooth.py', in_surf, 'lhpialt_smooth.vtk', '-n', '5', '-s', '1') assert ret.success