codebase / stimlib

Commit 54d21d09ab5ad045c331fc17f33c61d6764c3652

Authored by Jiaming Guo
1 parent 87d0a1d2

implementing finding eigenvalue and eigenvector for 3x3 matrix.

Showing 3 changed files with 366 additions and 0 deletions   Show diff stats
Inline Side-by-side
stim/math/matrix_sym.h 0 → 100644
  Show/Hide comments   View file @54d21d0
  1 +#ifndef STIM_MATRIX_SYM_H
  2 +#define STIM_MATRIX_SYM_H
  3 +
  4 +#include <stim/cuda/cudatools/callable.h>
  5 +#include <stim/math/matrix.h>
  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 <typename T, int D>
  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<T, D> mat() {
  52 + stim::matrix<T, D> 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<T, D> operator=(T rhs) {
  62 + int Nsq = D*(D+1)/2;
  63 + for(int i=0; i<Nsq; i++)
  64 + M[i] = rhs;
  65 +
  66 + return *this;
  67 + }
  68 +
  69 + CUDA_CALLABLE matrix_sym<T, D> operator=(matrix_sym<T, D> 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<T, D> &B, T rhs) {
  83 + int Nsq = D*(D+1)/2;
  84 + for(int i=0; i<Nsq; i++)
  85 + B.M[i] *= rhs;
  86 + }
  87 +
  88 + //return the tensor as a string
  89 + std::string str() {
  90 + std::stringstream ss;
  91 + for(int r = 0; r < D; r++){
  92 + ss << "| ";
  93 + for(int c=0; c<D; c++)
  94 + {
  95 + ss << (*this)(r, c) << " ";
  96 + }
  97 + ss << "|" << std::endl;
  98 + }
  99 +
  100 + return ss.str();
  101 + }
  102 +
  103 + //returns an identity matrix
  104 + static matrix_sym<T, D> identity() {
  105 + matrix_sym<T, D> I;
  106 + I = 0;
  107 + for (size_t i = 0; i < D; i++)
  108 + I.M[matrix_sym<T, D>::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
  Show/Hide comments   View file @54d21d0
  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<typename T>
  11 +class tensor2 : public matrix_sym<T, 2> {
  12 +
  13 +protected:
  14 +
  15 +public:
  16 +
  17 + //calculate the eigenvectors and eigenvalues of the tensor
  18 + CUDA_CALLABLE void eig(stim::matrix<T, 2>& v, stim::matrix<T, 2>& 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<T, 2>::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<T> operator=(stim::matrix_sym<T, 2> rhs){
  39 + stim::matrix_sym<T, 2>::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
  Show/Hide comments   View file @54d21d0
  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 201 \ No newline at end of file
... ...