/* Description: Tensor computation Class representing tensor as eigenvalues and eigenvectors; The conversion from tensor to eigensystem is done in tensor.h Basic functions to process eigenvalues and eigenvectors; _lambda1 definies if the system is equal 0 since if _lambda1 is zero the whole tensor has to be equal zero Author: Raimundo Sierra www.rsierra.com Copyright: Copyright (c) 2000 Raimundo Sierra. All rights reserved. LICENSE: This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. Institution: Surgical Planning Laboratory Department of Radiology Brigham and Women's Hospital Harvard Medical School MA 02115 USA Date: 11.2000 - 3.2001 Modifications: */ #ifndef _eigen #define _eigen 1 #include #include #include "tensor.h" const double PI = 3.1415926535897932384626433832795028841971693993751058209; const double GRAD2RAD = PI/180; class eigen{ private: float _lambda1, _lambda2, _lambda3; // eigenvalues 1 to 3, // 1 is largest, 3 smallest eigenvalue // corresponding eigenvectors float _vector1[3], _vector2[3], _vector3[3]; public: // CONSTRUCTOR DESTRUCTOR // constructor which initializes all values eigen(const float eig1, const float eig2, const float eig3, const float *vec1, const float *vec2, const float *vec3); eigen(); // default constructor, all values set to 0 ~eigen(); // destructor eigen(const eigen &sys); // Copy constructor // INPUT OUTPUT // set the eigenvalues and vectors void set(const float eig1, const float eig2, const float eig3, const float *vec1, const float *vec2, const float *vec3); // set the eigenvalues void set(const float eig1, const float eig2, const float eig3); void setVectors(const float *vec1, const float *vec2, const float *vec3); // set eigenvector 2 void setVector2(const float x, const float y, const float z); void resetVector3(); // set eigenvector 3 to vector1 x vector2 void setZero(); // set the eigenvalues and vectors to zero float value1() const; // largest eigenvalue float value2() const; // second eigenvalue float value3() const; // smallest eigenvalue float vector1(const int index) const; // returns element "index" (0|1|2) of eigenvector 1 float vector2(const int index) const; // returns element "index" (0|1|2) of eigenvector 2 float vector3(const int index) const; // returns element "index" (0|1|2) of eigenvector 3 void print() const; // outputs the eigensystem: eigenvalues and // corresponding eigenvectors // OPERATORS eigen& operator= (const eigen &e); // Assignment operator bool isZero() const; // check if eigensystem is zero (true if so) bool operator==(const eigen &e) const; // compares all elements of the eigensystem, no tolerance bool operator!=(const eigen &e) const; // dito bool operator< (const eigen &e) const; // less then, uses measurement specified in measure() bool operator> (const eigen &e) const; // dito bool operator<=(const eigen &e) const; // dito bool operator>=(const eigen &e) const; // dito // FUNCTIONS void scale(const float factor); // scale eigensystem by factor // scale differently in each direction void scale(const float x, const float y, const float z); // Rotate eigensystem angles are in Grad void rotate(const float pitch, const float roll, const float yaw); float det() const; // determinant of the tensor float relation_Lambda1_Lambda2() const; // absolute value of _lambda1/_lambda2 float relation_Lambda1_Lambda3() const; // absolute value of _lambda1/_lambda3 void eigen2tensor(tensor &t) const; // converts an eigensystem into a tensor int tensor2eigen(const tensor &t); // converts a tensor in an eigenvalue/vector system // negative return value is an error, 1 is success // measurement of how close diffusion tensor shape is to float close_line() const; // a line; (lambda1-lambda2)/lambda1 float close_plane() const; // a plane; (lambda2-lambda3)/lambda1 float close_sphere() const; // a sphere; lambda3/lambda1 // anisotropy measure describing deviation float anisotropy() const; // from the spherical case; 1 - (lambda3/lambda1) float diffusion_x() const; // total diffusion in x-direction float diffusion_y() const; // total diffusion in y-direction float diffusion_z() const; // total diffusion in z-direction float magnitude() const; // magnitude of eigensystem = diff_x*diff_y*diff_z // measurement used in functions <, <=, >, >= float measure() const {return this->magnitude(); } // functions between two eigensystems: // angle between the vectors of the maximum eigenvalues float cosinMaxEigen(const eigen &e) const; // returns the position of the closest voxel in the direction of int * closest(); // the largest eigenvector // compare 2 eigensystem, return 1 if identical, 0 if not similar at all bool isSimilar(const eigen &e, const float tolerance, const int similarity) const; bool hasSimilarShape(const eigen &e, const float tolerance) const; bool hasSimilarSize(const eigen &e, const float tolerance) const; bool hasSimilarDirection(const eigen &e, const float tolerance) const; bool threshold(const float threshold=0.8) const; // return true if eigensystem above treshold }; ostream &operator << (ostream &out, const eigen &sys); // output operator #endif