forked from 541435721/myVTKPythonLibrary
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcomputeFractionalAnisotropy.py
84 lines (65 loc) · 2.75 KB
/
computeFractionalAnisotropy.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
#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 myVTKPythonLibrary as myVTK
########################################################################
def computeFractionalAnisotropy(
farray_e1,
farray_e2,
farray_e3,
verbose=1):
myVTK.myPrint(verbose, "*** computeFractionalAnisotropy ***")
n_tuples = farray_e1.GetNumberOfTuples()
farray_FA = myVTK.createFloatArray("FA" , 1, n_tuples)
farray_FA12 = myVTK.createFloatArray("FA_12", 1, n_tuples)
farray_FA23 = myVTK.createFloatArray("FA_23", 1, n_tuples)
for k_tuple in xrange(n_tuples):
e1 = farray_e1.GetTuple1(k_tuple)
e2 = farray_e2.GetTuple1(k_tuple)
e3 = farray_e3.GetTuple1(k_tuple)
FA = ((e1-e2)**2+(e1-e3)**2+(e2-e3)**2)**(0.5) / (2*(e1**2+e2**2+e3**2))**(0.5)
FA12 = ((e1-e2)**2)**(0.5) / (e1**2+e2**2)**(0.5)
FA23 = ((e2-e3)**2)**(0.5) / (e2**2+e3**2)**(0.5)
farray_FA.SetTuple1(k_tuple, FA)
farray_FA12.SetTuple1(k_tuple, FA12)
farray_FA23.SetTuple1(k_tuple, FA23)
return (farray_FA,
farray_FA12,
farray_FA23)
########################################################################
def addFractionalAnisotropy(
ugrid,
field_name,
type_of_support="cell",
verbose=1):
myVTK.myPrint(verbose, "*** addFractionalAnisotropy ***")
if (type_of_support == "cell"):
data = ugrid.GetCellData()
elif (type_of_support == "point"):
data = ugrid.GetPointData()
farray_e1 = data.GetArray(field_name+"_Lmax")
farray_e2 = data.GetArray(field_name+"_Lmid")
farray_e3 = data.GetArray(field_name+"_Lmin")
(farray_FA,
farray_FA12,
farray_FA23) = computeFractionalAnisotropy(
farray_e1=farray_e1,
farray_e2=farray_e2,
farray_e3=farray_e3,
verbose=verbose-1)
farray_FA.SetName(field_name+"_FA")
farray_FA12.SetName(field_name+"_FA_12")
farray_FA23.SetName(field_name+"_FA_23")
data.AddArray(farray_FA)
data.AddArray(farray_FA12)
data.AddArray(farray_FA23)
return (farray_FA,
farray_FA12,
farray_FA23)