myVTKPythonLibrary/mapDataSetToPointData.py at master · gacevedobolton/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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#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 mapDataSetToPointData(
        mesh_from,
        type_of_support,
        mesh_to,
        farray_names,
        radius=1.,
        threshold_min=None,
        threshold_max=None,
        verbose=1):

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

    if (type_of_support == "point"):
        dataset = mesh_from.GetPointData()
        point_locator = getPointLocator(
            mesh_from)
    elif (type_of_support == "cell"):
        dataset = mesh_from.GetPointData()
        pdata_cell_centers_fr = getCellCenters(
            mesh=mesh_from,
            verbose=verbose-1)
        point_locator = getPointLocator(
            pdata_cell_centers_fr)

    n_points = mesh_to.GetNumberOfPoints()

    farrays_avg = {}
    farrays_std = {}
    farrays_rat = {}
    for farray_name in farray_names:
        assert (dataset.HasArray(farray_name)), "mesh has no array named " + farray_name + ". Aborting."

        farray_type = dataset.GetArray(farray_name).GetDataTypeAsString()
        farray_n_components = dataset.GetArray(farray_name).GetNumberOfComponents()
        farrays_avg[farray_name] = createArray(farray_name+"_avg",
                                               farray_n_components,
                                               n_points,
                                               farray_type)
        farrays_std[farray_name] = createArray(farray_name+"_std",
                                               farray_n_components,
                                               n_points,
                                               farray_type)
        farrays_rat[farray_name] = createArray(farray_name+"_rat",
                                               farray_n_components,
                                               n_points,
                                               farray_type)

    points_within_radius = vtk.vtkIdList()

    for k_point in xrange(n_points):
        point_locator.FindClosestNPoints(3,
                                         mesh_to.GetPoints().GetPoint(k_point),
                                         points_within_radius)
        #point_locator.FindPointsWithinRadius(radius,
                                             #mesh_to.GetPoints().GetPoint(k_point),
                                             #points_within_radius)

        for farray_name in farray_names:
            if (points_within_radius.GetNumberOfIds() > 0):
                values = [numpy.array(dataset.GetArray(farray_name).GetTuple(points_within_radius.GetId(k_id))) for k_id in xrange(points_within_radius.GetNumberOfIds())]
                #print "values = " + str(values)
                if (threshold_min != None):
                    values = [value for value in values if (numpy.linalg.norm(value) > threshold_min)]
                if (threshold_max != None):
                    values = [value for value in values if (numpy.linalg.norm(value) < threshold_max)]
                #print "values = " + str(values)
                if (len(values) > 0):
                    avg = numpy.mean(values, 0)
                    std = numpy.std(values, 0)
                    rat = std/avg
                else:
                    avg = [0]*n_components
                    std = [0]*n_components
                    rat = [0]*n_components
            else:
                avg = [0]*n_components
                std = [0]*n_components
                rat = [0]*n_components
            farrays_avg[farray_name].InsertTuple(k_point, avg)
            farrays_std[farray_name].InsertTuple(k_point, std)
            farrays_rat[farray_name].InsertTuple(k_point, rat)

    for farray_name in farray_names:
        mesh_to.GetPointData().AddArray(farrays_avg[farray_name])
        mesh_to.GetPointData().AddArray(farrays_std[farray_name])
        mesh_to.GetPointData().AddArray(farrays_rat[farray_name])