| |
1
| +#ifndef STIM_TENSOR3_H |
| |
2
| +#define STIM_TENSOR3_H |
| |
3
| + |
| |
4
| +#include "matrix_sym.h" |
| |
5
| +#include <stim/math/constants.h> |
| |
6
| + |
| |
7
| +namespace stim { |
| |
8
| + |
| |
9
| + /*This class represents a symmetric rank-2 2D tensor, useful for structure tensors |
| |
10
| + */ |
| |
11
| + |
| |
12
| + //Matrix ID cheat sheet |
| |
13
| + // | 0 1 2 | |
| |
14
| + // | 1 3 4 | |
| |
15
| + // | 2 4 5 | |
| |
16
| + template<typename T> |
| |
17
| + class tensor3 : public matrix_sym<T, 3> { |
| |
18
| + |
| |
19
| + protected: |
| |
20
| + |
| |
21
| + public: |
| |
22
| + |
| |
23
| + //calculates the determinant of the tensor |
| |
24
| + CUDA_CALLABLE T det() { |
| |
25
| + return M[0] * M[3] * M[5] + 2 * (M[1] * M[4] * M[2]) - M[2] * M[3] * M[2] - M[1] * M[1] * M[5] - M[0] * M[4] * M[4]; |
| |
26
| + } |
| |
27
| + |
| |
28
| + //calculate the eigenvalues for the tensor |
| |
29
| + //adapted from https://en.wikipedia.org/wiki/Eigenvalue_algorithm |
| |
30
| + |
| |
31
| + CUDA_CALLABLE stim::vec3<T> lambda() { |
| |
32
| + stim::vec3<T> lam; |
| |
33
| + T p1 = M[1] * M[1] + M[2] * M[2] + M[4] * M[4]; //calculate the sum of the squared off-diagonal values |
| |
34
| + if (p1 == 0) { //if this value is zero, the matrix is diagonal |
| |
35
| + lam[0] = M[0]; //the eigenvalues are the diagonal values |
| |
36
| + lam[1] = M[3]; |
| |
37
| + lam[2] = M[5]; |
| |
38
| + return lam; //return the eigenvalue vector |
| |
39
| + } |
| |
40
| + |
| |
41
| + T tr = matrix_sym<T, 3>::trace(); //calculate the trace of the matrix |
| |
42
| + T q = tr / 3; |
| |
43
| + T p2 = (M[0] - q) * (M[0] - q) + (M[3] - q) * (M[3] - q) + (M[5] - q) * (M[5] - q) + 2 * p1; |
| |
44
| + T p = sqrt(p2 / 6); |
| |
45
| + tensor3<T> Q; //allocate space for Q (q along the diagonals) |
| |
46
| + Q = (T)0; //initialize Q to zeros |
| |
47
| + Q(0, 0) = Q(1, 1) = Q(2, 2) = q; //set the diagonal values to q |
| |
48
| + tensor3<T> B = *this; // B1 = A |
| |
49
| + B.M[0] = (B.M[0] - q); |
| |
50
| + B.M[3] = (B.M[3] - q); |
| |
51
| + B.M[5] = (B.M[5] - q); |
| |
52
| + matrix_sym<T, 3>::operator_product(B, 1/p); // B = (1/p) * (A - q*I) |
| |
53
| + //B.M[0] = B.M[0] * 1/p; |
| |
54
| + //B.M[1] = B.M[1] * 1/p; |
| |
55
| + //B.M[2] = B.M[2] * 1/p; |
| |
56
| + //B.M[3] = B.M[3] * 1/p; |
| |
57
| + //B.M[4] = B.M[4] * 1/p; |
| |
58
| + //B.M[5] = B.M[5] * 1/p; |
| |
59
| + T r = B.det() / 2; //calculate det(B) / 2 |
| |
60
| + |
| |
61
| + // In exact arithmetic for a symmetric matrix - 1 <= r <= 1 |
| |
62
| + // but computation error can leave it slightly outside this range. |
| |
63
| + T phi; |
| |
64
| + if (r <= -1) phi = stim::PI / 3; |
| |
65
| + else if (r >= 1) phi = 0; |
| |
66
| + else phi = acos(r) / 3; |
| |
67
| + |
| |
68
| + // the eigenvalues satisfy eig3 >= eig2 >= eig1 |
| |
69
| + lam[2] = q + 2 * p * cos(phi); |
| |
70
| + lam[0] = q + 2 * p * cos(phi + (2 * stim::PI / 3)); |
| |
71
| + lam[1] = 3 * q - (lam[2] + lam[0]); |
| |
72
| + |
| |
73
| + return lam; |
| |
74
| + } |
| |
75
| + |
| |
76
| + CUDA_CALLABLE stim::matrix<T, 3> eig(stim::vec3<T>& lambda = stim::vec3<T>()) { |
| |
77
| + stim::matrix<T, 3> V; |
| |
78
| + |
| |
79
| + stim::matrix<T, 3> M1 = matrix_sym<T, 3>::mat(); |
| |
80
| + stim::matrix<T, 3> M2 = matrix_sym<T, 3>::mat(); |
| |
81
| + stim::matrix<T, 3> M3 = matrix_sym<T, 3>::mat(); // fill a tensor with symmetric values |
| |
82
| + |
| |
83
| + M1.operator_minus(M1, lambda[0]); // M1 = A - lambda[0] * I |
| |
84
| + |
| |
85
| + M2.operator_minus(M2, lambda[1]); // M2 = A - lambda[1] * I |
| |
86
| + |
| |
87
| + M3.operator_minus(M3, lambda[2]); // M3 = A - lambda[2] * I |
| |
88
| + |
| |
89
| + T Mod = 0; // module of one column |
| |
90
| + |
| |
91
| + T tmp1[9] = {0}; |
| |
92
| + for(int i = 0; i < 9; i++) { |
| |
93
| + for(int j = 0; j < 3; j++){ |
| |
94
| + tmp1[i] += M2(i%3, j) * M3(j, i/3); |
| |
95
| + } |
| |
96
| + } |
| |
97
| + if(tmp1[0] * tmp1[1] * tmp1[2] != 0) { // test whether it is zero column |
| |
98
| + Mod = sqrt(pow(tmp1[0],2) + pow(tmp1[1],2) + pow(tmp1[2],2)); |
| |
99
| + V(0, 0) = tmp1[0]/Mod; |
| |
100
| + V(1, 0) = tmp1[1]/Mod; |
| |
101
| + V(2, 0) = tmp1[2]/Mod; |
| |
102
| + } |
| |
103
| + else { |
| |
104
| + Mod = sqrt(pow(tmp1[3],2) + pow(tmp1[4],2) + pow(tmp1[5],2)); |
| |
105
| + V(0, 0) = tmp1[3]/Mod; |
| |
106
| + V(1, 0) = tmp1[4]/Mod; |
| |
107
| + V(2, 0) = tmp1[5]/Mod; |
| |
108
| + } |
| |
109
| + |
| |
110
| + T tmp2[9] = {0}; |
| |
111
| + for(int i = 0; i < 9; i++) { |
| |
112
| + for(int j = 0; j < 3; j++){ |
| |
113
| + tmp2[i] += M1(i%3, j) * M3(j, i/3); |
| |
114
| + } |
| |
115
| + } |
| |
116
| + if(tmp2[0] * tmp2[1] * tmp2[2] != 0) { |
| |
117
| + Mod = sqrt(pow(tmp2[0],2) + pow(tmp2[1],2) + pow(tmp2[2],2)); |
| |
118
| + V(0, 1) = tmp2[0]/Mod; |
| |
119
| + V(1, 1) = tmp2[1]/Mod; |
| |
120
| + V(2, 1) = tmp2[2]/Mod; |
| |
121
| + } |
| |
122
| + else { |
| |
123
| + Mod = sqrt(pow(tmp2[3],2) + pow(tmp2[4],2) + pow(tmp2[5],2)); |
| |
124
| + V(0, 1) = tmp2[3]/Mod; |
| |
125
| + V(1, 1) = tmp2[4]/Mod; |
| |
126
| + V(2, 1) = tmp2[5]/Mod; |
| |
127
| + } |
| |
128
| + |
| |
129
| + T tmp3[9] = {0}; |
| |
130
| + for(int i = 0; i < 9; i++) { |
| |
131
| + for(int j = 0; j < 3; j++){ |
| |
132
| + tmp3[i] += M1(i%3, j) * M2(j, i/3); |
| |
133
| + } |
| |
134
| + } |
| |
135
| + if(tmp3[0] * tmp3[1] * tmp3[2] != 0) { |
| |
136
| + Mod = sqrt(pow(tmp3[0],2) + pow(tmp3[1],2) + pow(tmp3[2],2)); |
| |
137
| + V(0, 2) = tmp3[0]/Mod; |
| |
138
| + V(1, 2) = tmp3[1]/Mod; |
| |
139
| + V(2, 2) = tmp3[2]/Mod; |
| |
140
| + } |
| |
141
| + else { |
| |
142
| + Mod = sqrt(pow(tmp3[3],2) + pow(tmp3[4],2) + pow(tmp3[5],2)); |
| |
143
| + V(0, 2) = tmp3[3]/Mod; |
| |
144
| + V(1, 2) = tmp3[4]/Mod; |
| |
145
| + V(2, 2) = tmp3[5]/Mod; |
| |
146
| + } |
| |
147
| + return V; //return the eigenvector matrix |
| |
148
| + } |
| |
149
| + // return one specific eigenvector |
| |
150
| + CUDA_CALLABLE stim::vec3<T> eig(int n, stim::vec3<T>& lambda = stim::vec3<T>()) { |
| |
151
| + stim::matrix<T, 3> V = eig(lambda); |
| |
152
| + stim::vec3<T> v; |
| |
153
| + for(int i = 0; i < 3; i++) |
| |
154
| + v[i] = V(i, n); |
| |
155
| + return v; |
| |
156
| + } |
| |
157
| + |
| |
158
| + |
| |
159
| + CUDA_CALLABLE T linear(stim::vec3<T>& lambda = stim::vec3<T>()) { |
| |
160
| + T cl = (lambda[2] - lambda[1]) / (lambda[0] + lambda[1] + lambda[2]); |
| |
161
| + return cl; |
| |
162
| + } |
| |
163
| + |
| |
164
| + CUDA_CALLABLE T Planar(stim::vec3<T>& lambda = stim::vec3<T>()) { |
| |
165
| + T cp = 2 * (lambda[1] - lambda[0]) / (lambda[0] + lambda[1] + lambda[2]); |
| |
166
| + return cp; |
| |
167
| + } |
| |
168
| + |
| |
169
| + CUDA_CALLABLE T spherical(stim::vec3<T>& lambda = stim::vec3<T>()) { |
| |
170
| + T cs = 3 * lambda[0] / (lambda[0] + lambda[1] + lambda[2]); |
| |
171
| + return cs; |
| |
172
| + } |
| |
173
| + |
| |
174
| + CUDA_CALLABLE T fa(stim::vec3<T>& lambda = stim::vec3<T>()) { |
| |
175
| + T fa = sqrt(1/2) * sqrt(pow(lambda[2] - lambda[1], 2) + pow(lambda[1] - lambda[0], 2) + pow(lambda[0] - lambda[2], 2)) / sqrt(pow(lambda[2], 2) + pow(lambda[1], 2) + pow(lambda[0], 2)); |
| |
176
| + } |
| |
177
| + //JACK 2: write functions to calculate anisotropy |
| |
178
| + //ex: fa(), linear(), planar(), spherical() |
| |
179
| + |
| |
180
| + |
| |
181
| + //calculate the eigenvectors and eigenvalues of the tensor |
| |
182
| + //CUDA_CALLABLE void eig(stim::matrix<T, 3>& v, stim::matrix<T, 3>& lambda){ |
| |
183
| + |
| |
184
| + //} |
| |
185
| + CUDA_CALLABLE tensor3<T> operator=(T rhs) { |
| |
186
| + stim::matrix_sym<T, 3>::operator=(rhs); |
| |
187
| + return *this; |
| |
188
| + } |
| |
189
| + |
| |
190
| + CUDA_CALLABLE tensor3<T> operator=(stim::matrix_sym<T, 3> rhs) { |
| |
191
| + stim::matrix_sym<T, 3>::operator=(rhs); |
| |
192
| + return *this; |
| |
193
| + } |
| |
194
| + }; |
| |
195
| + |
| |
196
| + |
| |
197
| +} //end namespace stim |
| |
198
| + |
| |
199
| + |
| |
200
| +#endif |
0
| \ No newline at end of file |
201
| \ No newline at end of file |