forked from gacevedobolton/myVTKPythonLibrary
-
Notifications
You must be signed in to change notification settings - Fork 0
/
rotateMatrix.py
91 lines (75 loc) · 3.8 KB
/
rotateMatrix.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
#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 myVTKPythonLibrary as myVTK
from mat_vec_tools import *
########################################################################
def rotateMatrix(
old_array,
old_array_storage="vec",
in_vecs=None,
R=None,
out_vecs=None,
verbose=0):
myVTK.myPrint(verbose, "*** rotateMatrix ***")
myVTK.myPrint(min(verbose,1), "*** Warning: in rotateMatrix, the definition of the global rotation is probably the inverse of the definition in previous rotateTensors function. ***")
n_components = old_array.GetNumberOfComponents()
if (old_array_storage == "vec"):
assert (n_components == 6), "Wrong numpber of components (n_components="+str(n_components)+"). Aborting."
elif (old_array_storage == "Cmat"):
assert (n_components == 9), "Wrong numpber of components (n_components="+str(n_components)+"). Aborting."
elif (old_array_storage == "Fmat"):
assert (n_components == 9), "Wrong numpber of components (n_components="+str(n_components)+"). Aborting."
else:
assert (0), "Wrong storage (old_array_storage="+str(old_array_storage)+"). Aborting."
n_tuples = old_array.GetNumberOfTuples()
new_array = myVTK.createFloatArray(old_array.GetName(), n_components, n_tuples)
new_vector = numpy.empty(n_components)
if (old_array_storage == "vec"):
old_vector = numpy.empty(6)
elif (old_array_storage == "Cmat"):
old_vector = numpy.empty(9)
elif (old_array_storage == "Fmat"):
old_vector = numpy.empty(9)
old_matrix = numpy.empty((3,3))
new_matrix = numpy.empty((3,3))
for k_tuple in xrange(n_tuples):
old_array.GetTuple(k_tuple, old_vector)
if (old_array_storage == "vec"):
vec_col6_to_mat_sym33(old_vector, old_matrix)
elif (old_array_storage == "Cmat"):
cvec9_to_mat33(old_vector, old_matrix)
elif (old_array_storage == "Fmat"):
fvec9_to_mat33(old_vector, old_matrix)
if (in_vecs is None):
in_R = numpy.eye(3)
else:
in_R = numpy.transpose(numpy.array([in_vecs[0].GetTuple(k_tuple),
in_vecs[1].GetTuple(k_tuple),
in_vecs[2].GetTuple(k_tuple)]))
if (out_vecs is None):
out_R = numpy.eye(3)
else:
out_R = numpy.transpose(numpy.array([out_vecs[0].GetTuple(k_tuple),
out_vecs[1].GetTuple(k_tuple),
out_vecs[2].GetTuple(k_tuple)]))
if (R is None):
R = numpy.eye(3)
full_R = numpy.dot(numpy.dot(numpy.transpose(in_R), R), out_R)
new_matrix[:] = numpy.dot(numpy.dot(numpy.transpose(full_R), old_matrix), full_R)
if (old_array_storage == "vec"):
mat_sym33_to_vec_col6(new_matrix, new_vector)
elif (old_array_storage == "Cmat"):
mat33_to_cvec9(new_matrix, new_vector)
elif (old_array_storage == "Fmat"):
mat33_to_fvec9(new_matrix, new_vector)
new_array.SetTuple(k_tuple, new_vector)
return new_array