Merge pull request #1015 from arnaudbore/add_surface_apply_transfo

[ENH] Uniformize surface apply transformation

scilpy/surfaces/surface_operations.py

@@ -1,11 +1,12 @@
 # -*- coding: utf-8 -*-
 
+from nibabel.affines import apply_affine
 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):
+def apply_transform(mesh, linear_transfo, deformation_img=None, inverse=False):
     """
     Apply transformation to a surface
         - Apply linear transformation with affine
@@ -15,37 +16,48 @@ def apply_transform(mesh, ants_affine, ants_warp=None):
     ----------
     mesh: trimeshpy - Triangle Mesh VTK class
         Moving surface
-    ants_affine: numpy.ndarray
+    linear_transfo: numpy.ndarray
         Transformation matrix to be applied
-    ants_warp: nib.Nifti1Image
+    deformation_img: nib.Nifti1Image
         Warp image from ANTs
+    inverse: boolean
+        Apply the inverse linear transformation.
 
     Returns
     -------
     mesh: trimeshpy - Triangle Mesh VTK class
         Surface moved
     """
     # Affine transformation
-    inv_affine = np.linalg.inv(ants_affine)
+    affine = linear_transfo
+    if inverse:
+        affine = np.linalg.inv(linear_transfo)
 
-    # Transform mesh vertices
-    mesh.set_vertices(mesh.vertices_affine(inv_affine))
+    # Apply affine transformation
+    flip = np.diag([-1, -1, 1, 1])
+    mesh.set_vertices(apply_affine(flip, mesh.get_vertices()))
+    mesh.set_vertices(apply_affine(affine, mesh.get_vertices()))
 
     # Flip triangle face, if needed
-    if mesh.is_transformation_flip(inv_affine):
+    if mesh.is_transformation_flip(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))
+    if deformation_img is not None:
+        deformation_data = np.squeeze(deformation_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_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)
+        inv_affine = np.linalg.inv(deformation_img.affine)
+        vts_vox = apply_affine(inv_affine, mesh.get_vertices())
+        tx = map_coordinates(deformation_data[..., 0], vts_vox.T, order=1)
+        ty = map_coordinates(deformation_data[..., 1], vts_vox.T, order=1)
+        tz = map_coordinates(deformation_data[..., 2], vts_vox.T, order=1)
 
         # Apply vertices translation in world coordinates
-        mesh.set_vertices(mesh.get_vertices() + np.array([tx, ty, tz]).T)
+        # LPS versus RAS (-tx, -ty)
+        mesh.set_vertices(mesh.get_vertices() + np.array([-tx, -ty, tz]).T)
+
+    mesh.set_vertices(apply_affine(flip, mesh.get_vertices()))
+    mesh.update_polydata()
 
     return mesh
 

scripts/scil_surface_apply_transform.py

@@ -3,12 +3,12 @@
 
 """
 Script to apply a transform to a surface (FreeSurfer or VTK supported),
-using output from ANTs registration tools (i.e. affine.txt, warp.nii.gz).
+using output from ANTs registration tools (i.e. vtk_transfo.txt, output1InverseWarp.nii.gz).
 
 Example usage from T1 to b0 using ANTs transforms:
 > ConvertTransformFile 3 output0GenericAffine.mat vtk_transfo.txt --hm
-> scil_surface_apply_transform.py lh_white_lps.vtk affine.txt lh_white_b0.vtk\\
-    --ants_warp warp.nii.gz
+> scil_surface_apply_transform.py lh_white_lps.vtk vtk_transfo.txt lh_white_b0.vtk\\
+    --in_deformation output1InverseWarp.nii.gz --inverse
 
 Important: The input surface needs to be in *T1 world LPS* coordinates
 (aligned over the T1 in MI-Brain).
@@ -24,13 +24,14 @@
 import logging
 
 import nibabel as nib
-import numpy as np
 from trimeshpy.io import load_mesh_from_file
+from trimeshpy.vtk_util import save_polydata
 
 from scilpy.io.utils import (add_overwrite_arg,
                              add_verbose_arg,
                              assert_inputs_exist,
-                             assert_outputs_exist)
+                             assert_outputs_exist,
+                             load_matrix_in_any_format)
 from scilpy.surfaces.surface_operations import apply_transform
 
 
@@ -45,17 +46,19 @@ def _build_arg_parser():
     p = argparse.ArgumentParser(description=__doc__, epilog=EPILOG,
                                 formatter_class=argparse.RawTextHelpFormatter)
 
-    p.add_argument('in_surface',
+    p.add_argument('in_moving_surface',
                    help='Input surface (.vtk).')
-
-    p.add_argument('ants_affine',
-                   help='Affine transform from ANTs (.txt or .mat).')
-
+    p.add_argument('in_transfo',
+                   help='Path of the file containing the 4x4 \n'
+                        'transformation, matrix (.txt, .npy or .mat).')
     p.add_argument('out_surface',
                    help='Output surface (.vtk).')
 
-    p.add_argument('--ants_warp',
-                   help='Warp image from ANTs (Nifti image).')
+    g = p.add_argument_group("Transformation options")
+    g.add_argument('--inverse', action='store_true',
+                   help='Apply the inverse linear transformation.')
+    g.add_argument('--in_deformation', metavar='file',
+                   help='Path to the file containing a deformation field.')
 
     add_verbose_arg(p)
     add_overwrite_arg(p)
@@ -68,26 +71,25 @@ def main():
     args = parser.parse_args()
     logging.getLogger().setLevel(logging.getLevelName(args.verbose))
 
-    assert_inputs_exist(parser, [args.in_surface, args.ants_affine],
-                        args.ants_warp)
+    assert_inputs_exist(parser, [args.in_moving_surface, args.in_transfo],
+                        args.in_deformation)
     assert_outputs_exist(parser, args, args.out_surface)
 
     # Load mesh
-    mesh = load_mesh_from_file(args.in_surface)
+    mesh = load_mesh_from_file(args.in_moving_surface)
 
-    warp_img = None
-    # Load affine
-    affine = np.loadtxt(args.ants_affine)
+    # Load transformation
+    transfo = load_matrix_in_any_format(args.in_transfo)
 
-    # Non linear transformation
-    if args.ants_warp:
-        # Load warp
-        warp_img = nib.load(args.ants_warp)
+    deformation_data = None
+    if args.in_deformation is not None:
+        deformation_data = nib.load(args.in_deformation)
 
-    mesh = apply_transform(mesh, affine, warp_img)
+    mesh = apply_transform(mesh, transfo, deformation_data,
+                           inverse=args.inverse)
 
     # Save mesh
-    mesh.save(args.out_surface)
+    save_polydata(mesh.__polydata__, args.out_surface, legacy_vtk_format=True)
 
 
 if __name__ == "__main__":

scripts/tests/test_surface_apply_transform.py

@@ -24,5 +24,5 @@ def test_execution_surface_vtk_fib(script_runner, monkeypatch):
     in_aff = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
                           'affine.txt')
     ret = script_runner.run('scil_surface_apply_transform.py', in_surf,
-                            in_aff, 'lhpialt_lin.vtk')
+                            in_aff, 'lhpialt_lin.vtk', '--inverse')
     assert ret.success