JMATH::JSVD3D Class Reference

Singular value decomposition. More...

#include <JSVD3D.hh>

Classes
struct	JGivens
	Givens quaternion. More...
struct	JQuaternion
	Auxiliary class for quaternion computation. More...

Public Member Functions
	JSVD3D ()
	Default constructor. More...
	JSVD3D (const JMatrix3D &A)
	Constructor. More...
void	decompose (const JMatrix3D &A)
	Decompose given matrix. More...
const JMatrix3D &	invert (const double precision=1.0e-12) const
	Get inverted matrix. More...

Public Attributes
JMatrix3D	U
JMatrix3D	S
JMatrix3D	V

Static Private Member Functions
static double	getLengthSquared (const double x, const double y, const double z)
	Get length squared. More...

Private Attributes
JMatrix3S	B

Detailed Description

Singular value decomposition.

See: "Computing the Singular Value Decomposition of 3x3 matrices with minimal branching and elementary floating point operations",
A. McAdams, A. Selle, R. Tamstorf, J. Teran and E. Sifakis,
University of Wisconsin - Madison technical report TR1690, May 2011

Definition at line 27 of file JSVD3D.hh.

Constructor & Destructor Documentation

JMATH::JSVD3D::JSVD3D ( )

inline

Default constructor.

Definition at line 32 of file JSVD3D.hh.

    {}

JMATH::JSVD3D::JSVD3D ( const JMatrix3D &  A )

inline

Constructor.

Parameters

A

matrix

Definition at line 41 of file JSVD3D.hh.

    {
      decompose(A);
    }

Member Function Documentation

void JMATH::JSVD3D::decompose ( const JMatrix3D &  A )

inline

Decompose given matrix.

Parameters

A

matrix

Definition at line 52 of file JSVD3D.hh.

    {
      using namespace std;

      // B = A^T * A

      B.a00  =  A.a00 * A.a00  +  A.a10 * A.a10  +  A.a20 * A.a20;
    
      B.a10  =  A.a01 * A.a00  +  A.a11 * A.a10  +  A.a21 * A.a20;
      B.a11  =  A.a01 * A.a01  +  A.a11 * A.a11  +  A.a21 * A.a21;

      B.a20  =  A.a02 * A.a00  +  A.a12 * A.a10  +  A.a22 * A.a20;
      B.a21  =  A.a02 * A.a01  +  A.a12 * A.a11  +  A.a22 * A.a21;
      B.a22  =  A.a02 * A.a02  +  A.a12 * A.a12  +  A.a22 * A.a22;

      B.a01  =  B.a10;
      B.a02  =  B.a20;
      B.a12  =  B.a21;


      // symmetric eigen alysis

      const JQuaternion q(B);

      const double w  = q[3];
      const double x  = q[0];
      const double y  = q[1];
      const double z  = q[2];

      const double xx = x*x;
      const double yy = y*y;
      const double zz = z*z;
      const double xz = x*z;
      const double xy = x*y;
      const double yz = y*z;
      const double wx = w*x;
      const double wy = w*y;
      const double wz = w*z;

      V.a00 = 1.0 - 2.0*(yy + zz);  V.a01 = 2.0*(xy - wz);        V.a02 = 2.0*(xz + wy);
      V.a10 = 2.0*(xy + wz);        V.a11 = 1.0 - 2.0*(xx + zz);  V.a12 = 2.0*(yz - wx);
      V.a20 = 2.0*(xz - wy);        V.a21 = 2.0*(yz + wx);        V.a22 = 1.0 - 2.0*(xx + yy);


      B.mul(A, V);


      // sort singular values.

      double rho1 = getLengthSquared(B.a00, B.a10, B.a20);
      double rho2 = getLengthSquared(B.a01, B.a11, B.a21);
      double rho3 = getLengthSquared(B.a02, B.a12, B.a22);

      if (rho1 < rho2) {
        swap(B.a00,B.a01);   B.a01 = -B.a01;   swap(V.a00,V.a01);   V.a01 = -V.a01;
        swap(B.a10,B.a11);   B.a11 = -B.a11;   swap(V.a10,V.a11);   V.a11 = -V.a11;
        swap(B.a20,B.a21);   B.a21 = -B.a21;   swap(V.a20,V.a21);   V.a21 = -V.a21;
        swap(rho1,rho2);  
      }

      if (rho1 < rho3) {
        swap(B.a00,B.a02);   B.a02 = -B.a02;   swap(V.a00,V.a02);   V.a02 = -V.a02;
        swap(B.a10,B.a12);   B.a12 = -B.a12;   swap(V.a10,V.a12);   V.a12 = -V.a12;
        swap(B.a20,B.a22);   B.a22 = -B.a22;   swap(V.a20,V.a22);   V.a22 = -V.a22;
        swap(rho1,rho3);  
      }

      if (rho2 < rho3) {
        swap(B.a01,B.a02);   B.a02 = -B.a02;   swap(V.a01,V.a02);   V.a02 = -V.a02;
        swap(B.a11,B.a12);   B.a12 = -B.a12;   swap(V.a11,V.a12);   V.a12 = -V.a12;
        swap(B.a21,B.a22);   B.a22 = -B.a22;   swap(V.a21,V.a22);   V.a22 = -V.a22;
      }


      // QR decomposition

      JMatrix3D& Q = U;
      JMatrix3D& R = S;

      double a, b;
    
      const JGivens q1(B.a00, B.a10);      // 1st Givens rotation (ch,0,0,sh)

      a = q1.geta();
      b = q1.getb();

      // apply B = Q' * B

      R.a00 =  a*B.a00 + b*B.a10;  R.a01 =  a*B.a01 + b*B.a11;  R.a02 =  a*B.a02 + b*B.a12;
      R.a10 = -b*B.a00 + a*B.a10;  R.a11 = -b*B.a01 + a*B.a11;  R.a12 = -b*B.a02 + a*B.a12;
      R.a20 =  B.a20;              R.a21 =  B.a21;              R.a22 =  B.a22;
  
      const JGivens q2(R.a00, R.a20);      // 2nd Givens rotation (ch,0,-sh,0)

      a = q2.geta();
      b = q2.getb();

      // apply B = Q' * B;

      B.a00 =  a*R.a00 + b*R.a20;  B.a01 =  a*R.a01 + b*R.a21;  B.a02 =  a*R.a02 + b*R.a22;
      B.a10 =  R.a10;              B.a11 =  R.a11;              B.a12 =  R.a12;
      B.a20 = -b*R.a00 + a*R.a20;  B.a21 = -b*R.a01 + a*R.a21;  B.a22 = -b*R.a02 + a*R.a22;

      const JGivens q3(B.a11, B.a21);      // 3rd Givens rotation (ch,sh,0,0)

      a = q3.geta();
      b = q3.getb();

      // R is now set to desired value

      R.a00 =  B.a00;              R.a01 =  B.a01;              R.a02 =  B.a02;
      R.a10 =  a*B.a10+b*B.a20;    R.a11 =  a*B.a11 + b*B.a21;  R.a12 =  a*B.a12 + b*B.a22;
      R.a20 = -b*B.a10+a*B.a20;    R.a21 = -b*B.a11 + a*B.a21;  R.a22 = -b*B.a12 + a*B.a22;

      // Q = Q1 * Q2 * Q3

      const double sp1 = -1.0 + 2.0*q1.sh*q1.sh;
      const double sp2 = -1.0 + 2.0*q2.sh*q2.sh;
      const double sp3 = -1.0 + 2.0*q3.sh*q3.sh;

      Q.a00 =  sp1*sp2; 
      Q.a01 =  4.0*q2.ch*q3.ch*sp1*q2.sh*q3.sh  +  2.0*q1.ch*q1.sh*sp3; 
      Q.a02 =  4.0*q1.ch*q3.ch*q1.sh*q3.sh      -  2.0*q2.ch*sp1*q2.sh*sp3;

      Q.a10 = -2.0*q1.ch*q1.sh*sp2;
      Q.a11 = -8.0*q1.ch*q2.ch*q3.ch*q1.sh*q2.sh*q3.sh  +  sp1*sp3; 
      Q.a12 = -2.0*q3.ch*q3.sh                          +  4.0*q1.sh*(q3.ch*q1.sh*q3.sh + q1.ch*q2.ch*q2.sh*sp3);
    
      Q.a20 =  2.0*q2.ch*q2.sh; 
      Q.a21 = -2.0*q3.ch*sp2*q3.sh; 
      Q.a22 =  sp2*sp3;
    }

const JMatrix3D& JMATH::JSVD3D::invert ( const double  precision = 1.0e-12 ) const

inline

Get inverted matrix.

Parameters

precision

precision

Returns

inverted matrix

Definition at line 192 of file JSVD3D.hh.

    {
      double w = fabs(S.a00);

      if (fabs(S.a11) > w) { w = fabs(S.a11); }
      if (fabs(S.a22) > w) { w = fabs(S.a22); }

      w *= precision;

      const double s00 = (fabs(S.a00) >= w ? 1.0 / S.a00 : 0.0);
      const double s11 = (fabs(S.a11) >= w ? 1.0 / S.a11 : 0.0);
      const double s22 = (fabs(S.a22) >= w ? 1.0 / S.a22 : 0.0);

      // B = V * S^-1

      B.a00 = V.a00*s00;  B.a01 = V.a01*s11;  B.a02 = V.a02*s22;
      B.a10 = V.a10*s00;  B.a11 = V.a11*s11;  B.a12 = V.a12*s22;
      B.a20 = V.a20*s00;  B.a21 = V.a21*s11;  B.a22 = V.a22*s22;

      U.transpose();

      // B = V * S^-1 * U^T

      B.mul(U);

      U.transpose();

      return B;
    }

static double JMATH::JSVD3D::getLengthSquared ( const double  x, const double  y, const double  z  )

inlinestaticprivate

Get length squared.

Parameters

x	x
y	y
z	z

Returns

square of length

Definition at line 365 of file JSVD3D.hh.

    {
      return x*x + y*y + z*z;
    }

Member Data Documentation

JMatrix3D JMATH::JSVD3D::U

mutable

Definition at line 222 of file JSVD3D.hh.

JMatrix3D JMATH::JSVD3D::S

mutable

Definition at line 223 of file JSVD3D.hh.

JMatrix3D JMATH::JSVD3D::V

mutable

Definition at line 224 of file JSVD3D.hh.

JMatrix3S JMATH::JSVD3D::B

mutableprivate

Definition at line 228 of file JSVD3D.hh.

---

The documentation for this class was generated from the following file:

• JSVD3D.hh