forked from 541435721/myVTKPythonLibrary
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcomputeFiberDirections.py
128 lines (105 loc) · 4.35 KB
/
computeFiberDirections.py
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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#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 math
import random
import numpy
import myVTKPythonLibrary as myVTK
########################################################################
def computeFiberDirections(
farray_eRR,
farray_eCC,
farray_eLL,
farray_angle_helix,
angles_in_degrees=True,
use_new_definition=False,
shuffle_vectors=False,
verbose=1):
myVTK.myPrint(verbose, "*** computeFiberDirections ***")
n_tuples = farray_angle_helix.GetNumberOfTuples()
farray_eF = myVTK.createFloatArray("eF", 3, n_tuples)
farray_eS = myVTK.createFloatArray("eS", 3, n_tuples)
farray_eN = myVTK.createFloatArray("eN", 3, n_tuples)
eRR = numpy.empty(3)
eCC = numpy.empty(3)
eLL = numpy.empty(3)
for k_tuple in xrange(n_tuples):
farray_eRR.GetTuple(k_tuple, eRR)
farray_eCC.GetTuple(k_tuple, eCC)
farray_eLL.GetTuple(k_tuple, eLL)
assert (round(numpy.linalg.norm(eRR),1) == 1.0),\
"|eRR| = "+str(numpy.linalg.norm(eRR))+"≠ 1. Aborting"
assert (round(numpy.linalg.norm(eCC),1) == 1.0),\
"|eCC| = "+str(numpy.linalg.norm(eCC))+"≠ 1. Aborting"
assert (round(numpy.linalg.norm(eLL),1) == 1.0),\
"|eLL| = "+str(numpy.linalg.norm(eLL))+"≠ 1. Aborting"
angle_helix = farray_angle_helix.GetTuple1(k_tuple)
if (angles_in_degrees): angle_helix = angle_helix*math.pi/180
eF = math.cos(angle_helix) * eCC + math.sin(angle_helix) * eLL
#print "eF = "+str(eF)
if (shuffle_vectors):
eF *= random.choice([-1,+1])
#print "eF = "+str(eF)
if (use_new_definition):
eN = eRR
if (shuffle_vectors):
eN *= random.choice([-1,+1])
assert (abs(numpy.dot(eN,eRR)) > 0.999)
eS = numpy.cross(eN, eF)
else:
eS = eRR
if (shuffle_vectors): eS *= random.choice([-1,+1])
eN = numpy.cross(eF, eS)
assert (round(numpy.linalg.norm(eF),1) == 1.0),\
"|eF| = "+str(numpy.linalg.norm(eF))+"≠ 1. Aborting"
assert (round(numpy.linalg.norm(eS),1) == 1.0),\
"|eS| = "+str(numpy.linalg.norm(eS))+"≠ 1. Aborting"
assert (round(numpy.linalg.norm(eN),1) == 1.0),\
"|eN| = "+str(numpy.linalg.norm(eN))+"≠ 1. Aborting"
farray_eF.SetTuple(k_tuple, eF)
farray_eS.SetTuple(k_tuple, eS)
farray_eN.SetTuple(k_tuple, eN)
return (farray_eF,
farray_eS,
farray_eN)
########################################################################
def addFiberDirections(
ugrid,
type_of_support="cell",
angles_in_degrees=True,
use_new_definition=False,
shuffle_vectors=False,
verbose=1):
myVTK.myPrint(verbose, "*** addFiberDirections ***")
if (type_of_support == "cell"):
ugrid_data = ugrid.GetCellData()
elif (type_of_support == "point"):
ugrid_data = ugrid.GetPointData()
farray_eRR = ugrid_data.GetArray("eRR")
farray_eCC = ugrid_data.GetArray("eCC")
farray_eLL = ugrid_data.GetArray("eLL")
farray_angle_helix = ugrid_data.GetArray("angle_helix")
(farray_eF,
farray_eS,
farray_eN) = computeFiberDirections(
farray_eRR=farray_eRR,
farray_eCC=farray_eCC,
farray_eLL=farray_eLL,
farray_angle_helix=farray_angle_helix,
angles_in_degrees=angles_in_degrees,
use_new_definition=use_new_definition,
shuffle_vectors=shuffle_vectors,
verbose=verbose-1)
ugrid_data.AddArray(farray_eF)
ugrid_data.AddArray(farray_eS)
ugrid_data.AddArray(farray_eN)
return (farray_eF,
farray_eS,
farray_eN)