myVTKPythonLibrary/addPDataNormals.py at master · 541435721/myVTKPythonLibrary · 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
65
66
#coding=utf8

########################################################################
###                                                                  ###
### Created by Martin Genet, 2012-2016                               ###
###                                                                  ###
### University of California at San Francisco (UCSF), USA            ###
### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
### École Polytechnique, Palaiseau, France                           ###
###                                                                  ###
########################################################################

import numpy
import vtk

import myVTKPythonLibrary as myVTK

########################################################################

def addPDataNormals(
        pdata,
        orient_outward=1,
        verbose=1):

    myVTK.myPrint(verbose, "*** addPDataNormals ***")

    poly_data_normals = vtk.vtkPolyDataNormals()
    poly_data_normals.ComputePointNormalsOff()
    poly_data_normals.ComputeCellNormalsOn()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        poly_data_normals.SetInputData(pdata)
    else:
        poly_data_normals.SetInput(pdata)
    poly_data_normals.Update()
    pdata = poly_data_normals.GetOutput()

    if (orient_outward):
        cell_centers = myVTK.getCellCenters(
            mesh=pdata,
            verbose=verbose-1)
        cell_center = numpy.empty(3)

        mesh_center = numpy.array(pdata.GetCenter())

        normals = pdata.GetCellData().GetNormals()
        normal = numpy.empty(3)

        cnt_pos = 0
        cnt_neg = 0
        for k_cell in xrange(pdata.GetNumberOfCells()):
            cell_centers.GetPoint(k_cell, cell_center)
            outward  = cell_center-mesh_center
            outward /= numpy.linalg.norm(outward)
            normals.GetTuple(k_cell, normal)
            proj = numpy.dot(outward, normal)
            if (proj > 0): cnt_pos += 1
            else:          cnt_neg += 1
        #print cnt_pos
        #print cnt_neg

        if (cnt_neg > cnt_pos):
            poly_data_normals.FlipNormalsOn()
            poly_data_normals.Update()
            pdata = poly_data_normals.GetOutput()

    return pdata