Commit 27db2f14 authored by Max Lyon's avatar Max Lyon
Browse files

removed dependency to gmm by replacing eigenvalue computation with code taken from stack overflow

parent a054dc8c
......@@ -485,19 +485,19 @@ draw_principal_component(const PrincipalComponent& _pc)
// choose color due to eigenvalue sign
if(_pc.sign[i] == true)
{
glColor3f(1.f, 0.f, 0.f);
glColor3f(1.f, 0.f, 0.f);
GLfloat mat_amb_diff[] = {1.0, 0.0, 0.0, 1.0};
glMaterialfv( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_amb_diff);
glColor4fv( mat_amb_diff);
GLfloat mat_amb_diff[] = {1.0, 0.0, 0.0, 1.0};
glMaterialfv( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_amb_diff);
glColor4fv( mat_amb_diff);
}
else
{
glColor3f(0.f, 0.f, 1.f);
glColor3f(0.f, 0.f, 1.f);
GLfloat mat_amb_diff[] = {0.0, 0.0, 1.0, 1.0};
glMaterialfv( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_amb_diff);
glColor4fv( mat_amb_diff);
GLfloat mat_amb_diff[] = {0.0, 0.0, 1.0, 1.0};
glMaterialfv( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_amb_diff);
glColor4fv( mat_amb_diff);
}
}
else
......@@ -516,7 +516,7 @@ draw_principal_component(const PrincipalComponent& _pc)
draw_arrow( a , radius);
if( show_tensor_component_[i] == 2)
draw_arrow( a*(-1), radius);
draw_arrow( a*(-1), radius);
}
else
{
......@@ -524,7 +524,7 @@ draw_principal_component(const PrincipalComponent& _pc)
draw_line( a , width);
if( show_tensor_component_[i] == 2)
draw_line( a*(-1), width);
draw_line( a*(-1), width);
}
}
......@@ -644,6 +644,113 @@ void PrincipalAxisNode::get_axes_colors(Vec4f out_colors[3]) const {
}
//----------------------------------------------------------------------------
// Following code is taken from http://stackoverflow.com/questions/4372224/fast-method-for-computing-3x3-symmetric-matrix-spectral-decomposition:
// Slightly modified version of Stan Melax's code for 3x3 matrix diagonalization (Thanks Stan!)
// source: http://www.melax.com/diag.html?attredirects=0
void PrincipalAxisNode::diagonalize(const double (&A)[3][3], double (&Q)[3][3], double (&D)[3][3])
{
// A must be a symmetric matrix.
// returns Q and D such that
// Diagonal matrix D = QT * A * Q; and A = Q*D*QT
const int maxsteps=24; // certainly wont need that many.
int k0, k1, k2;
double o[3], m[3];
double q [4] = {0.0,0.0,0.0,1.0};
double jr[4];
double sqw, sqx, sqy, sqz;
double tmp1, tmp2, mq;
double AQ[3][3];
double thet, sgn, t, c;
for(int i=0;i < maxsteps;++i)
{
// quat to matrix
sqx = q[0]*q[0];
sqy = q[1]*q[1];
sqz = q[2]*q[2];
sqw = q[3]*q[3];
Q[0][0] = ( sqx - sqy - sqz + sqw);
Q[1][1] = (-sqx + sqy - sqz + sqw);
Q[2][2] = (-sqx - sqy + sqz + sqw);
tmp1 = q[0]*q[1];
tmp2 = q[2]*q[3];
Q[1][0] = 2.0 * (tmp1 + tmp2);
Q[0][1] = 2.0 * (tmp1 - tmp2);
tmp1 = q[0]*q[2];
tmp2 = q[1]*q[3];
Q[2][0] = 2.0 * (tmp1 - tmp2);
Q[0][2] = 2.0 * (tmp1 + tmp2);
tmp1 = q[1]*q[2];
tmp2 = q[0]*q[3];
Q[2][1] = 2.0 * (tmp1 + tmp2);
Q[1][2] = 2.0 * (tmp1 - tmp2);
// AQ = A * Q
AQ[0][0] = Q[0][0]*A[0][0]+Q[1][0]*A[0][1]+Q[2][0]*A[0][2];
AQ[0][1] = Q[0][1]*A[0][0]+Q[1][1]*A[0][1]+Q[2][1]*A[0][2];
AQ[0][2] = Q[0][2]*A[0][0]+Q[1][2]*A[0][1]+Q[2][2]*A[0][2];
AQ[1][0] = Q[0][0]*A[0][1]+Q[1][0]*A[1][1]+Q[2][0]*A[1][2];
AQ[1][1] = Q[0][1]*A[0][1]+Q[1][1]*A[1][1]+Q[2][1]*A[1][2];
AQ[1][2] = Q[0][2]*A[0][1]+Q[1][2]*A[1][1]+Q[2][2]*A[1][2];
AQ[2][0] = Q[0][0]*A[0][2]+Q[1][0]*A[1][2]+Q[2][0]*A[2][2];
AQ[2][1] = Q[0][1]*A[0][2]+Q[1][1]*A[1][2]+Q[2][1]*A[2][2];
AQ[2][2] = Q[0][2]*A[0][2]+Q[1][2]*A[1][2]+Q[2][2]*A[2][2];
// D = Qt * AQ
D[0][0] = AQ[0][0]*Q[0][0]+AQ[1][0]*Q[1][0]+AQ[2][0]*Q[2][0];
D[0][1] = AQ[0][0]*Q[0][1]+AQ[1][0]*Q[1][1]+AQ[2][0]*Q[2][1];
D[0][2] = AQ[0][0]*Q[0][2]+AQ[1][0]*Q[1][2]+AQ[2][0]*Q[2][2];
D[1][0] = AQ[0][1]*Q[0][0]+AQ[1][1]*Q[1][0]+AQ[2][1]*Q[2][0];
D[1][1] = AQ[0][1]*Q[0][1]+AQ[1][1]*Q[1][1]+AQ[2][1]*Q[2][1];
D[1][2] = AQ[0][1]*Q[0][2]+AQ[1][1]*Q[1][2]+AQ[2][1]*Q[2][2];
D[2][0] = AQ[0][2]*Q[0][0]+AQ[1][2]*Q[1][0]+AQ[2][2]*Q[2][0];
D[2][1] = AQ[0][2]*Q[0][1]+AQ[1][2]*Q[1][1]+AQ[2][2]*Q[2][1];
D[2][2] = AQ[0][2]*Q[0][2]+AQ[1][2]*Q[1][2]+AQ[2][2]*Q[2][2];
o[0] = D[1][2];
o[1] = D[0][2];
o[2] = D[0][1];
m[0] = fabs(o[0]);
m[1] = fabs(o[1]);
m[2] = fabs(o[2]);
k0 = (m[0] > m[1] && m[0] > m[2])?0: (m[1] > m[2])? 1 : 2; // index of largest element of offdiag
k1 = (k0+1)%3;
k2 = (k0+2)%3;
if (o[k0]==0.0)
{
break; // diagonal already
}
thet = (D[k2][k2]-D[k1][k1])/(2.0*o[k0]);
sgn = (thet > 0.0)?1.0:-1.0;
thet *= sgn; // make it positive
t = sgn /(thet +((thet < 1.E6)?sqrt(thet*thet+1.0):thet)) ; // sign(T)/(|T|+sqrt(T^2+1))
c = 1.0/sqrt(t*t+1.0); // c= 1/(t^2+1) , t=s/c
if(c==1.0)
{
break; // no room for improvement - reached machine precision.
}
jr[0 ] = jr[1] = jr[2] = jr[3] = 0.0;
jr[k0] = sgn*sqrt((1.0-c)/2.0); // using 1/2 angle identity sin(a/2) = sqrt((1-cos(a))/2)
jr[k0] *= -1.0; // since our quat-to-matrix convention was for v*M instead of M*v
jr[3 ] = sqrt(1.0f - jr[k0] * jr[k0]);
if(jr[3]==1.0)
{
break; // reached limits of floating point precision
}
q[0] = (q[3]*jr[0] + q[0]*jr[3] + q[1]*jr[2] - q[2]*jr[1]);
q[1] = (q[3]*jr[1] - q[0]*jr[2] + q[1]*jr[3] + q[2]*jr[0]);
q[2] = (q[3]*jr[2] + q[0]*jr[1] - q[1]*jr[0] + q[2]*jr[3]);
q[3] = (q[3]*jr[3] - q[0]*jr[0] - q[1]*jr[1] - q[2]*jr[2]);
mq = sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]);
q[0] /= mq;
q[1] /= mq;
q[2] /= mq;
q[3] /= mq;
}
}
//=============================================================================
} // namespace SceneGraph
} // namespace ACG
......
......@@ -233,7 +233,9 @@ public:
void get_axes_colors(Vec4f out_colors[3]) const;
private:
void diagonalize(const double (&A)[3][3], double (&Q)[3][3], double (&D)[3][3]);
// vector of Principal Components
std::vector< PrincipalComponent > pc_;
......
......@@ -52,7 +52,6 @@
#include "PrincipalAxisNode.hh"
#include <ACG/GL/gl.hh>
#include <gmm/gmm.h>
//== NAMESPACES ===============================================================
......@@ -69,7 +68,7 @@ void
PrincipalAxisNode::
set_vector( unsigned int _i, const Vec3d _p, const VectorT& _v)
{
gmm::dense_matrix<double> m(3,3);
Matrix4x4T<double> m;
m(0,0) = _v[0];
m(1,1) = _v[1];
......@@ -86,23 +85,28 @@ set_vector( unsigned int _i, const Vec3d _p, const VectorT& _v)
template<class MatrixT>
void
void
PrincipalAxisNode::
set_matrix( unsigned int _i, const Vec3d _p, const MatrixT& _m)
{
// create gmm matrix
gmm::dense_matrix<double> s(3,3);
// create matrix
double s[3][3];
// copy values
for(unsigned int i=0; i<3; ++i)
for(unsigned int j=0; j<3; ++j)
s(i,j) = _m(i,j);
s[i][j] = _m(i,j);
// compute eigenvalues and eigenvectors
std::vector<double> eigval(3);
gmm::dense_matrix<double> eigvect(3,3);
double Q[3][3];
double D[3][3];
diagonalize(s, Q, D);
std::vector<double> eigval;
for (int i = 0; i < 3; ++i)
eigval.push_back(D[i][i]);
gmm::symmetric_qr_algorithm(s, eigval, eigvect);
unsigned int pstress_idx[3];
if( fabs(eigval[0]) >= fabs(eigval[1]) && fabs(eigval[0]) >= fabs(eigval[2])) pstress_idx[0] = 0;
......@@ -128,23 +132,23 @@ set_matrix( unsigned int _i, const Vec3d _p, const MatrixT& _m)
for(unsigned int j=0; j<3; ++j)
{
// set sign
if( eigval[ pstress_idx[j]] < 0 )
if( eigval[ pstress_idx[j]] < 0 )
sign[j] = false;
else
sign[j] = true;
for(unsigned int k=0; k<3; ++k)
a[j][k] = eigval[ pstress_idx[j]]*eigvect(k, pstress_idx[j]);
a[j][k] = eigval[ pstress_idx[j]]*Q[k][pstress_idx[j]];
}
set(_i,
set(_i,
PrincipalComponent(_p,
a[0],
a[1],
a[2],
sign[0],
sign[1],
sign[2] ) );
a[0],
a[1],
a[2],
sign[0],
sign[1],
sign[2] ) );
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment