python_for_brain/gii2polydata.py at master · spherik/python_for_brain · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
from vtk import vtkPolyData, vtkFloatArray, vtkPoints, vtkIdTypeArray
from vtk import vtkCellArray, vtkIntArray, vtkSmoothPolyDataFilter, vtkPolyDataNormals
import numpy as np

def gii2polydata(mesh, smooth=False, normals=False):
    """ Convert a Gifti mesh to vtk polydata
        This function converts a Gifti mesh into a VTK polydata.

        :param gii: Gifti mesh
        :type gii: 'nibabel.gifti.gifti.GiftiImage'
        :param smooth: falg to smooth mesh
        :type smooth: bool
        :param normals: flag to compute normals
        :type normals: bool
        :return: vtk mesh
        :rtype: vtk.vtkPolydata
    """
    for da in mesh.darrays:
        if da.intent == 1008:  # nib.nifti1.intent_codes.NIFTI_INTENT_POINTSET
            # Create vtk array from nib data
            points_data = vtkFloatArray()
            points_data.SetNumberOfComponents(3)
            points_data.SetArray(da.data, da.data.shape[0] * 3, 1)

            # Set vtk array as points
            points = vtkPoints()
            points.SetData(points_data)

        elif da.intent == 1009:  # nib.nifti1.intent_codes.NIFTI_INTENT_TRIANGLE
            # To represent triangles, cells need to be [num_idxs, idx0, idx1, idx2], with num_idxs = 3
            triangles_array = np.hstack([np.ones((da.data.shape[0], 1)) * 3, np.array(da.data)]).astype(np.longlong, copy=False)
            # Create vtk array from nib data
            triangles_data = vtkIdTypeArray()
            triangles_data.SetNumberOfComponents(4)
            triangles_data.SetArray(triangles_array, da.data.shape[0] * 4, 1)

            # Set vtk array as cells (triangles)
            triangles = vtkCellArray()
            triangles.SetCells(da.data.shape[0], triangles_data)

    polydata = vtkPolyData()
    polydata.SetPoints(points)
    polydata.SetPolys(triangles)

    if smooth:
        smoothFilter = vtkSmoothPolyDataFilter()
        smoothFilter.SetInputData(polydata)
        smoothFilter.SetNumberOfIterations(15)
        smoothFilter.SetRelaxationFactor(0.1)
        smoothFilter.FeatureEdgeSmoothingOn()
        smoothFilter.BoundarySmoothingOff()
        smoothFilter.Update()
        polydata = smoothFilter.GetOutput()

    if normals:
        normalGenerator = vtkPolyDataNormals()
        normalGenerator.ComputePointNormalsOn()
        normalGenerator.ComputeCellNormalsOn()
        normalGenerator.SetInputData(polydata)
        normalGenerator.AutoOrientNormalsOn()
        normalGenerator.Update()
        polydata = normalGenerator.GetOutput()

    return polydata