From 54d21d09ab5ad045c331fc17f33c61d6764c3652 Mon Sep 17 00:00:00 2001 From: Jiaming Guo Date: Mon, 21 Nov 2016 16:57:03 -0600 Subject: [PATCH] implementing finding eigenvalue and eigenvector for 3x3 matrix. --- stim/math/matrix_sym.h | 118 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ stim/math/tensor2.h | 48 ++++++++++++++++++++++++++++++++++++++++++++++++ stim/math/tensor3.h | 200 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 366 insertions(+), 0 deletions(-) create mode 100644 stim/math/matrix_sym.h create mode 100644 stim/math/tensor2.h create mode 100644 stim/math/tensor3.h diff --git a/stim/math/matrix_sym.h b/stim/math/matrix_sym.h new file mode 100644 index 0000000..4cac755 --- /dev/null +++ b/stim/math/matrix_sym.h @@ -0,0 +1,118 @@ +#ifndef STIM_MATRIX_SYM_H +#define STIM_MATRIX_SYM_H + +#include +#include + +/* This class represents a rank 2, 3-dimensional tensor viable +for representing tensor fields such as structure and diffusion tensors +*/ +namespace stim{ + +template +class matrix_sym{ + +protected: + //values are stored in column-major order as a lower-triangular matrix + T M[D*(D + 1)/2]; + + static size_t idx(size_t r, size_t c) { + //if the index is in the upper-triangular portion, swap the indices + if(r < c){ + size_t t = r; + r = c; + c = t; + } + + size_t ci = (c + 1) * (D + (D - c))/2 - 1; //index to the end of column c + size_t i = ci - (D - r - 1); + return i; + } + + //calculate the row and column given an index + //static void indices(size_t& r, size_t& c, size_t idx) { + // size_t col = 0; + // for ( ; col < D; col++) + // if(idx <= ((D - col + D) * (col + 1)/2 - 1)) + // break; + + // c = col; + // size_t ci = (D - (col - 1) + D) * col / 2 - 1; //index to the end of last column col -1 + // r = idx - ci + c - 1; + //} + static void indices(size_t& r, size_t& c, size_t idx) { + size_t cf = -1/2 * sqrt(4 * D * D + 4 * D - (7 + 8 * idx)) + D - 1/2; + c = ceil(cf); + r = idx - D * c + c * (c + 1) / 2; + } + +public: + //return the symmetric matrix associated with this tensor + stim::matrix mat() { + stim::matrix r; + r.setsym(M); + return r; + } + + CUDA_CALLABLE T& operator()(int r, int c) { + return M[idx(r, c)]; + } + + CUDA_CALLABLE matrix_sym operator=(T rhs) { + int Nsq = D*(D+1)/2; + for(int i=0; i operator=(matrix_sym rhs) { + size_t N = D * (D + 1) / 2; + for (size_t i = 0; i < N; i++) M[i] = rhs.M[i]; + return *this; + } + + CUDA_CALLABLE T trace() { + T tr = 0; + for (size_t i = 0; i < D; i++) //for each diagonal value + tr += M[idx(i, i)]; //add the value on the diagonal + return tr; + } + // overload matrix multiply scalar + CUDA_CALLABLE void operator_product(matrix_sym &B, T rhs) { + int Nsq = D*(D+1)/2; + for(int i=0; i identity() { + matrix_sym I; + I = 0; + for (size_t i = 0; i < D; i++) + I.M[matrix_sym::idx(i, i)] = 1; + return I; + } +}; + + + +} //end namespace stim + + +#endif diff --git a/stim/math/tensor2.h b/stim/math/tensor2.h new file mode 100644 index 0000000..fc0f639 --- /dev/null +++ b/stim/math/tensor2.h @@ -0,0 +1,48 @@ +#ifndef STIM_TENSOR2_H +#define STIM_TENSOR2_H + +#include "matrix_sym.h" + +namespace stim { + +/*This class represents a symmetric rank-2 2D tensor, useful for structure tensors +*/ +template +class tensor2 : public matrix_sym { + +protected: + +public: + + //calculate the eigenvectors and eigenvalues of the tensor + CUDA_CALLABLE void eig(stim::matrix& v, stim::matrix& lambda) { + + lambda = 0; //initialize the eigenvalue matrix to zero + + T t = M[0] + M[2]; //calculate the trace of the tensor + T d = M[0] * M[2] - M[1] * M[1]; //calculate the determinant of the tensor + + lambda(0, 0) = t / 2 + sqrt(t*t / 4 - d); + lambda(1, 1) = t / 2 - sqrt(t*t / 4 - d); + + if (M[1] == 0) { + v = stim::matrix::identity(); + } + else { + v(0, 0) = lambda(0, 0) - d; + v(0, 1) = lambda(1, 1) - d; + v(1, 0) = v(1, 1) = M[1]; + } + } + + CUDA_CALLABLE tensor2 operator=(stim::matrix_sym rhs){ + stim::matrix_sym::operator=(rhs); + return *this; + } +}; + + +} //end namespace stim + + +#endif \ No newline at end of file diff --git a/stim/math/tensor3.h b/stim/math/tensor3.h new file mode 100644 index 0000000..f0ed778 --- /dev/null +++ b/stim/math/tensor3.h @@ -0,0 +1,200 @@ +#ifndef STIM_TENSOR3_H +#define STIM_TENSOR3_H + +#include "matrix_sym.h" +#include + +namespace stim { + + /*This class represents a symmetric rank-2 2D tensor, useful for structure tensors + */ + + //Matrix ID cheat sheet + // | 0 1 2 | + // | 1 3 4 | + // | 2 4 5 | + template + class tensor3 : public matrix_sym { + + protected: + + public: + + //calculates the determinant of the tensor + CUDA_CALLABLE T det() { + 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]; + } + + //calculate the eigenvalues for the tensor + //adapted from https://en.wikipedia.org/wiki/Eigenvalue_algorithm + + CUDA_CALLABLE stim::vec3 lambda() { + stim::vec3 lam; + T p1 = M[1] * M[1] + M[2] * M[2] + M[4] * M[4]; //calculate the sum of the squared off-diagonal values + if (p1 == 0) { //if this value is zero, the matrix is diagonal + lam[0] = M[0]; //the eigenvalues are the diagonal values + lam[1] = M[3]; + lam[2] = M[5]; + return lam; //return the eigenvalue vector + } + + T tr = matrix_sym::trace(); //calculate the trace of the matrix + T q = tr / 3; + T p2 = (M[0] - q) * (M[0] - q) + (M[3] - q) * (M[3] - q) + (M[5] - q) * (M[5] - q) + 2 * p1; + T p = sqrt(p2 / 6); + tensor3 Q; //allocate space for Q (q along the diagonals) + Q = (T)0; //initialize Q to zeros + Q(0, 0) = Q(1, 1) = Q(2, 2) = q; //set the diagonal values to q + tensor3 B = *this; // B1 = A + B.M[0] = (B.M[0] - q); + B.M[3] = (B.M[3] - q); + B.M[5] = (B.M[5] - q); + matrix_sym::operator_product(B, 1/p); // B = (1/p) * (A - q*I) + //B.M[0] = B.M[0] * 1/p; + //B.M[1] = B.M[1] * 1/p; + //B.M[2] = B.M[2] * 1/p; + //B.M[3] = B.M[3] * 1/p; + //B.M[4] = B.M[4] * 1/p; + //B.M[5] = B.M[5] * 1/p; + T r = B.det() / 2; //calculate det(B) / 2 + + // In exact arithmetic for a symmetric matrix - 1 <= r <= 1 + // but computation error can leave it slightly outside this range. + T phi; + if (r <= -1) phi = stim::PI / 3; + else if (r >= 1) phi = 0; + else phi = acos(r) / 3; + + // the eigenvalues satisfy eig3 >= eig2 >= eig1 + lam[2] = q + 2 * p * cos(phi); + lam[0] = q + 2 * p * cos(phi + (2 * stim::PI / 3)); + lam[1] = 3 * q - (lam[2] + lam[0]); + + return lam; + } + + CUDA_CALLABLE stim::matrix eig(stim::vec3& lambda = stim::vec3()) { + stim::matrix V; + + stim::matrix M1 = matrix_sym::mat(); + stim::matrix M2 = matrix_sym::mat(); + stim::matrix M3 = matrix_sym::mat(); // fill a tensor with symmetric values + + M1.operator_minus(M1, lambda[0]); // M1 = A - lambda[0] * I + + M2.operator_minus(M2, lambda[1]); // M2 = A - lambda[1] * I + + M3.operator_minus(M3, lambda[2]); // M3 = A - lambda[2] * I + + T Mod = 0; // module of one column + + T tmp1[9] = {0}; + for(int i = 0; i < 9; i++) { + for(int j = 0; j < 3; j++){ + tmp1[i] += M2(i%3, j) * M3(j, i/3); + } + } + if(tmp1[0] * tmp1[1] * tmp1[2] != 0) { // test whether it is zero column + Mod = sqrt(pow(tmp1[0],2) + pow(tmp1[1],2) + pow(tmp1[2],2)); + V(0, 0) = tmp1[0]/Mod; + V(1, 0) = tmp1[1]/Mod; + V(2, 0) = tmp1[2]/Mod; + } + else { + Mod = sqrt(pow(tmp1[3],2) + pow(tmp1[4],2) + pow(tmp1[5],2)); + V(0, 0) = tmp1[3]/Mod; + V(1, 0) = tmp1[4]/Mod; + V(2, 0) = tmp1[5]/Mod; + } + + T tmp2[9] = {0}; + for(int i = 0; i < 9; i++) { + for(int j = 0; j < 3; j++){ + tmp2[i] += M1(i%3, j) * M3(j, i/3); + } + } + if(tmp2[0] * tmp2[1] * tmp2[2] != 0) { + Mod = sqrt(pow(tmp2[0],2) + pow(tmp2[1],2) + pow(tmp2[2],2)); + V(0, 1) = tmp2[0]/Mod; + V(1, 1) = tmp2[1]/Mod; + V(2, 1) = tmp2[2]/Mod; + } + else { + Mod = sqrt(pow(tmp2[3],2) + pow(tmp2[4],2) + pow(tmp2[5],2)); + V(0, 1) = tmp2[3]/Mod; + V(1, 1) = tmp2[4]/Mod; + V(2, 1) = tmp2[5]/Mod; + } + + T tmp3[9] = {0}; + for(int i = 0; i < 9; i++) { + for(int j = 0; j < 3; j++){ + tmp3[i] += M1(i%3, j) * M2(j, i/3); + } + } + if(tmp3[0] * tmp3[1] * tmp3[2] != 0) { + Mod = sqrt(pow(tmp3[0],2) + pow(tmp3[1],2) + pow(tmp3[2],2)); + V(0, 2) = tmp3[0]/Mod; + V(1, 2) = tmp3[1]/Mod; + V(2, 2) = tmp3[2]/Mod; + } + else { + Mod = sqrt(pow(tmp3[3],2) + pow(tmp3[4],2) + pow(tmp3[5],2)); + V(0, 2) = tmp3[3]/Mod; + V(1, 2) = tmp3[4]/Mod; + V(2, 2) = tmp3[5]/Mod; + } + return V; //return the eigenvector matrix + } + // return one specific eigenvector + CUDA_CALLABLE stim::vec3 eig(int n, stim::vec3& lambda = stim::vec3()) { + stim::matrix V = eig(lambda); + stim::vec3 v; + for(int i = 0; i < 3; i++) + v[i] = V(i, n); + return v; + } + + + CUDA_CALLABLE T linear(stim::vec3& lambda = stim::vec3()) { + T cl = (lambda[2] - lambda[1]) / (lambda[0] + lambda[1] + lambda[2]); + return cl; + } + + CUDA_CALLABLE T Planar(stim::vec3& lambda = stim::vec3()) { + T cp = 2 * (lambda[1] - lambda[0]) / (lambda[0] + lambda[1] + lambda[2]); + return cp; + } + + CUDA_CALLABLE T spherical(stim::vec3& lambda = stim::vec3()) { + T cs = 3 * lambda[0] / (lambda[0] + lambda[1] + lambda[2]); + return cs; + } + + CUDA_CALLABLE T fa(stim::vec3& lambda = stim::vec3()) { + 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)); + } + //JACK 2: write functions to calculate anisotropy + //ex: fa(), linear(), planar(), spherical() + + + //calculate the eigenvectors and eigenvalues of the tensor + //CUDA_CALLABLE void eig(stim::matrix& v, stim::matrix& lambda){ + + //} + CUDA_CALLABLE tensor3 operator=(T rhs) { + stim::matrix_sym::operator=(rhs); + return *this; + } + + CUDA_CALLABLE tensor3 operator=(stim::matrix_sym rhs) { + stim::matrix_sym::operator=(rhs); + return *this; + } + }; + + +} //end namespace stim + + +#endif \ No newline at end of file -- libgit2 0.21.4