Authored by David Mayerich
2 parents 87d0a1d2 54d21d09

### Merge branch 'JACK' into 'master'

```implementing finding eigenvalue and eigenvector for 3x3 matrix.

hi,
I implemented the eigenvalue and eigenvector as well as the variance tensor part.
Fixed some minor error in matrix_sym.h when using in cooperating with matrix.h.

See merge request !11```
Showing 3 changed files with 366 additions and 0 deletions
stim/math/matrix_sym.h 0 → 100644

 1 +#ifndef STIM_MATRIX_SYM_H 2 +#define STIM_MATRIX_SYM_H 3 + 4 +#include 5 +#include 6 + 7 +/* This class represents a rank 2, 3-dimensional tensor viable 8 +for representing tensor fields such as structure and diffusion tensors 9 +*/ 10 +namespace stim{ 11 + 12 +template 13 +class matrix_sym{ 14 + 15 +protected: 16 + //values are stored in column-major order as a lower-triangular matrix 17 + T M[D*(D + 1)/2]; 18 + 19 + static size_t idx(size_t r, size_t c) { 20 + //if the index is in the upper-triangular portion, swap the indices 21 + if(r < c){ 22 + size_t t = r; 23 + r = c; 24 + c = t; 25 + } 26 + 27 + size_t ci = (c + 1) * (D + (D - c))/2 - 1; //index to the end of column c 28 + size_t i = ci - (D - r - 1); 29 + return i; 30 + } 31 + 32 + //calculate the row and column given an index 33 + //static void indices(size_t& r, size_t& c, size_t idx) { 34 + // size_t col = 0; 35 + // for ( ; col < D; col++) 36 + // if(idx <= ((D - col + D) * (col + 1)/2 - 1)) 37 + // break; 38 + 39 + // c = col; 40 + // size_t ci = (D - (col - 1) + D) * col / 2 - 1; //index to the end of last column col -1 41 + // r = idx - ci + c - 1; 42 + //} 43 + static void indices(size_t& r, size_t& c, size_t idx) { 44 + size_t cf = -1/2 * sqrt(4 * D * D + 4 * D - (7 + 8 * idx)) + D - 1/2; 45 + c = ceil(cf); 46 + r = idx - D * c + c * (c + 1) / 2; 47 + } 48 + 49 +public: 50 + //return the symmetric matrix associated with this tensor 51 + stim::matrix mat() { 52 + stim::matrix r; 53 + r.setsym(M); 54 + return r; 55 + } 56 + 57 + CUDA_CALLABLE T& operator()(int r, int c) { 58 + return M[idx(r, c)]; 59 + } 60 + 61 + CUDA_CALLABLE matrix_sym operator=(T rhs) { 62 + int Nsq = D*(D+1)/2; 63 + for(int i=0; i operator=(matrix_sym rhs) { 70 + size_t N = D * (D + 1) / 2; 71 + for (size_t i = 0; i < N; i++) M[i] = rhs.M[i]; 72 + return *this; 73 + } 74 + 75 + CUDA_CALLABLE T trace() { 76 + T tr = 0; 77 + for (size_t i = 0; i < D; i++) //for each diagonal value 78 + tr += M[idx(i, i)]; //add the value on the diagonal 79 + return tr; 80 + } 81 + // overload matrix multiply scalar 82 + CUDA_CALLABLE void operator_product(matrix_sym &B, T rhs) { 83 + int Nsq = D*(D+1)/2; 84 + for(int i=0; i identity() { 105 + matrix_sym I; 106 + I = 0; 107 + for (size_t i = 0; i < D; i++) 108 + I.M[matrix_sym::idx(i, i)] = 1; 109 + return I; 110 + } 111 +}; 112 + 113 + 114 + 115 +} //end namespace stim 116 + 117 + 118 +#endif ... ...
stim/math/tensor2.h 0 → 100644

 1 +#ifndef STIM_TENSOR2_H 2 +#define STIM_TENSOR2_H 3 + 4 +#include "matrix_sym.h" 5 + 6 +namespace stim { 7 + 8 +/*This class represents a symmetric rank-2 2D tensor, useful for structure tensors 9 +*/ 10 +template 11 +class tensor2 : public matrix_sym { 12 + 13 +protected: 14 + 15 +public: 16 + 17 + //calculate the eigenvectors and eigenvalues of the tensor 18 + CUDA_CALLABLE void eig(stim::matrix& v, stim::matrix& lambda) { 19 + 20 + lambda = 0; //initialize the eigenvalue matrix to zero 21 + 22 + T t = M[0] + M[2]; //calculate the trace of the tensor 23 + T d = M[0] * M[2] - M[1] * M[1]; //calculate the determinant of the tensor 24 + 25 + lambda(0, 0) = t / 2 + sqrt(t*t / 4 - d); 26 + lambda(1, 1) = t / 2 - sqrt(t*t / 4 - d); 27 + 28 + if (M[1] == 0) { 29 + v = stim::matrix::identity(); 30 + } 31 + else { 32 + v(0, 0) = lambda(0, 0) - d; 33 + v(0, 1) = lambda(1, 1) - d; 34 + v(1, 0) = v(1, 1) = M[1]; 35 + } 36 + } 37 + 38 + CUDA_CALLABLE tensor2 operator=(stim::matrix_sym rhs){ 39 + stim::matrix_sym::operator=(rhs); 40 + return *this; 41 + } 42 +}; 43 + 44 + 45 +} //end namespace stim 46 + 47 + 48 +#endif 0 49 \ No newline at end of file ... ...
stim/math/tensor3.h 0 → 100644

 1 +#ifndef STIM_TENSOR3_H 2 +#define STIM_TENSOR3_H 3 + 4 +#include "matrix_sym.h" 5 +#include 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 17 + class tensor3 : public matrix_sym { 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 lambda() { 32 + stim::vec3 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::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 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 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::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 eig(stim::vec3& lambda = stim::vec3()) { 77 + stim::matrix V; 78 + 79 + stim::matrix M1 = matrix_sym::mat(); 80 + stim::matrix M2 = matrix_sym::mat(); 81 + stim::matrix M3 = matrix_sym::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 eig(int n, stim::vec3& lambda = stim::vec3()) { 151 + stim::matrix V = eig(lambda); 152 + stim::vec3 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& lambda = stim::vec3()) { 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& lambda = stim::vec3()) { 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& lambda = stim::vec3()) { 170 + T cs = 3 * lambda[0] / (lambda[0] + lambda[1] + lambda[2]); 171 + return cs; 172 + } 173 + 174 + CUDA_CALLABLE T fa(stim::vec3& lambda = stim::vec3()) { 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& v, stim::matrix& lambda){ 183 + 184 + //} 185 + CUDA_CALLABLE tensor3 operator=(T rhs) { 186 + stim::matrix_sym::operator=(rhs); 187 + return *this; 188 + } 189 + 190 + CUDA_CALLABLE tensor3 operator=(stim::matrix_sym rhs) { 191 + stim::matrix_sym::operator=(rhs); 192 + return *this; 193 + } 194 + }; 195 + 196 + 197 +} //end namespace stim 198 + 199 + 200 +#endif 0 201 \ No newline at end of file ... ...