Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
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.
 
 JSVD3D (const JMatrix3D &A)
 Constructor.
 
void decompose (const JMatrix3D &A)
 Decompose given matrix.
 
const JMatrix3Dinvert (const double precision=1.0e-12) const
 Get inverted matrix.
 

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.
 

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

◆ JSVD3D() [1/2]

JMATH::JSVD3D::JSVD3D ( )
inline

Default constructor.

Definition at line 32 of file JSVD3D.hh.

33 {}

◆ JSVD3D() [2/2]

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

Constructor.

Parameters
Amatrix

Definition at line 41 of file JSVD3D.hh.

42 {
43 decompose(A);
44 }
void decompose(const JMatrix3D &A)
Decompose given matrix.
Definition JSVD3D.hh:52

Member Function Documentation

◆ decompose()

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

Decompose given matrix.

Parameters
Amatrix

Definition at line 52 of file JSVD3D.hh.

53 {
54 using namespace std;
55
56 // B = A^T * A
57
58 B.a00 = A.a00 * A.a00 + A.a10 * A.a10 + A.a20 * A.a20;
59
60 B.a10 = A.a01 * A.a00 + A.a11 * A.a10 + A.a21 * A.a20;
61 B.a11 = A.a01 * A.a01 + A.a11 * A.a11 + A.a21 * A.a21;
62
63 B.a20 = A.a02 * A.a00 + A.a12 * A.a10 + A.a22 * A.a20;
64 B.a21 = A.a02 * A.a01 + A.a12 * A.a11 + A.a22 * A.a21;
65 B.a22 = A.a02 * A.a02 + A.a12 * A.a12 + A.a22 * A.a22;
66
67 B.a01 = B.a10;
68 B.a02 = B.a20;
69 B.a12 = B.a21;
70
71
72 // symmetric eigen alysis
73
74 const JQuaternion q(B);
75
76 const double w = q[3];
77 const double x = q[0];
78 const double y = q[1];
79 const double z = q[2];
80
81 const double xx = x*x;
82 const double yy = y*y;
83 const double zz = z*z;
84 const double xz = x*z;
85 const double xy = x*y;
86 const double yz = y*z;
87 const double wx = w*x;
88 const double wy = w*y;
89 const double wz = w*z;
90
91 V.a00 = 1.0 - 2.0*(yy + zz); V.a01 = 2.0*(xy - wz); V.a02 = 2.0*(xz + wy);
92 V.a10 = 2.0*(xy + wz); V.a11 = 1.0 - 2.0*(xx + zz); V.a12 = 2.0*(yz - wx);
93 V.a20 = 2.0*(xz - wy); V.a21 = 2.0*(yz + wx); V.a22 = 1.0 - 2.0*(xx + yy);
94
95
96 B.mul(A, V);
97
98
99 // sort singular values.
100
101 double rho1 = getLengthSquared(B.a00, B.a10, B.a20);
102 double rho2 = getLengthSquared(B.a01, B.a11, B.a21);
103 double rho3 = getLengthSquared(B.a02, B.a12, B.a22);
104
105 if (rho1 < rho2) {
106 swap(B.a00,B.a01); B.a01 = -B.a01; swap(V.a00,V.a01); V.a01 = -V.a01;
107 swap(B.a10,B.a11); B.a11 = -B.a11; swap(V.a10,V.a11); V.a11 = -V.a11;
108 swap(B.a20,B.a21); B.a21 = -B.a21; swap(V.a20,V.a21); V.a21 = -V.a21;
109 swap(rho1,rho2);
110 }
111
112 if (rho1 < rho3) {
113 swap(B.a00,B.a02); B.a02 = -B.a02; swap(V.a00,V.a02); V.a02 = -V.a02;
114 swap(B.a10,B.a12); B.a12 = -B.a12; swap(V.a10,V.a12); V.a12 = -V.a12;
115 swap(B.a20,B.a22); B.a22 = -B.a22; swap(V.a20,V.a22); V.a22 = -V.a22;
116 swap(rho1,rho3);
117 }
118
119 if (rho2 < rho3) {
120 swap(B.a01,B.a02); B.a02 = -B.a02; swap(V.a01,V.a02); V.a02 = -V.a02;
121 swap(B.a11,B.a12); B.a12 = -B.a12; swap(V.a11,V.a12); V.a12 = -V.a12;
122 swap(B.a21,B.a22); B.a22 = -B.a22; swap(V.a21,V.a22); V.a22 = -V.a22;
123 }
124
125
126 // QR decomposition
127
128 JMatrix3D& Q = U;
129 JMatrix3D& R = S;
130
131 double a, b;
132
133 const JGivens q1(B.a00, B.a10); // 1st Givens rotation (ch,0,0,sh)
134
135 a = q1.geta();
136 b = q1.getb();
137
138 // apply B = Q' * B
139
140 R.a00 = a*B.a00 + b*B.a10; R.a01 = a*B.a01 + b*B.a11; R.a02 = a*B.a02 + b*B.a12;
141 R.a10 = -b*B.a00 + a*B.a10; R.a11 = -b*B.a01 + a*B.a11; R.a12 = -b*B.a02 + a*B.a12;
142 R.a20 = B.a20; R.a21 = B.a21; R.a22 = B.a22;
143
144 const JGivens q2(R.a00, R.a20); // 2nd Givens rotation (ch,0,-sh,0)
145
146 a = q2.geta();
147 b = q2.getb();
148
149 // apply B = Q' * B;
150
151 B.a00 = a*R.a00 + b*R.a20; B.a01 = a*R.a01 + b*R.a21; B.a02 = a*R.a02 + b*R.a22;
152 B.a10 = R.a10; B.a11 = R.a11; B.a12 = R.a12;
153 B.a20 = -b*R.a00 + a*R.a20; B.a21 = -b*R.a01 + a*R.a21; B.a22 = -b*R.a02 + a*R.a22;
154
155 const JGivens q3(B.a11, B.a21); // 3rd Givens rotation (ch,sh,0,0)
156
157 a = q3.geta();
158 b = q3.getb();
159
160 // R is now set to desired value
161
162 R.a00 = B.a00; R.a01 = B.a01; R.a02 = B.a02;
163 R.a10 = a*B.a10+b*B.a20; R.a11 = a*B.a11 + b*B.a21; R.a12 = a*B.a12 + b*B.a22;
164 R.a20 = -b*B.a10+a*B.a20; R.a21 = -b*B.a11 + a*B.a21; R.a22 = -b*B.a12 + a*B.a22;
165
166 // Q = Q1 * Q2 * Q3
167
168 const double sp1 = -1.0 + 2.0*q1.sh*q1.sh;
169 const double sp2 = -1.0 + 2.0*q2.sh*q2.sh;
170 const double sp3 = -1.0 + 2.0*q3.sh*q3.sh;
171
172 Q.a00 = sp1*sp2;
173 Q.a01 = 4.0*q2.ch*q3.ch*sp1*q2.sh*q3.sh + 2.0*q1.ch*q1.sh*sp3;
174 Q.a02 = 4.0*q1.ch*q3.ch*q1.sh*q3.sh - 2.0*q2.ch*sp1*q2.sh*sp3;
175
176 Q.a10 = -2.0*q1.ch*q1.sh*sp2;
177 Q.a11 = -8.0*q1.ch*q2.ch*q3.ch*q1.sh*q2.sh*q3.sh + sp1*sp3;
178 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);
179
180 Q.a20 = 2.0*q2.ch*q2.sh;
181 Q.a21 = -2.0*q3.ch*sp2*q3.sh;
182 Q.a22 = sp2*sp3;
183 }
JMatrix3D S
Definition JSVD3D.hh:223
static double getLengthSquared(const double x, const double y, const double z)
Get length squared.
Definition JSVD3D.hh:365
JMatrix3D V
Definition JSVD3D.hh:224
JMatrix3D U
Definition JSVD3D.hh:222
const double a

◆ invert()

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

Get inverted matrix.

Parameters
precisionprecision
Returns
inverted matrix

Definition at line 192 of file JSVD3D.hh.

193 {
194 double w = fabs(S.a00);
195
196 if (fabs(S.a11) > w) { w = fabs(S.a11); }
197 if (fabs(S.a22) > w) { w = fabs(S.a22); }
198
199 w *= precision;
200
201 const double s00 = (fabs(S.a00) >= w ? 1.0 / S.a00 : 0.0);
202 const double s11 = (fabs(S.a11) >= w ? 1.0 / S.a11 : 0.0);
203 const double s22 = (fabs(S.a22) >= w ? 1.0 / S.a22 : 0.0);
204
205 // B = V * S^-1
206
207 B.a00 = V.a00*s00; B.a01 = V.a01*s11; B.a02 = V.a02*s22;
208 B.a10 = V.a10*s00; B.a11 = V.a11*s11; B.a12 = V.a12*s22;
209 B.a20 = V.a20*s00; B.a21 = V.a21*s11; B.a22 = V.a22*s22;
210
211 U.transpose();
212
213 // B = V * S^-1 * U^T
214
215 B.mul(U);
216
217 U.transpose();
218
219 return B;
220 }
JMatrix3D & mul(const double factor)
Scale matrix.
JMatrix3D & transpose()
Transpose.
JMatrix3S B
Definition JSVD3D.hh:228

◆ getLengthSquared()

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

Get length squared.

Parameters
xx
yy
zz
Returns
square of length

Definition at line 365 of file JSVD3D.hh.

366 {
367 return x*x + y*y + z*z;
368 }

Member Data Documentation

◆ U

JMatrix3D JMATH::JSVD3D::U
mutable

Definition at line 222 of file JSVD3D.hh.

◆ S

JMatrix3D JMATH::JSVD3D::S
mutable

Definition at line 223 of file JSVD3D.hh.

◆ V

JMatrix3D JMATH::JSVD3D::V
mutable

Definition at line 224 of file JSVD3D.hh.

◆ B

JMatrix3S JMATH::JSVD3D::B
mutableprivate

Definition at line 228 of file JSVD3D.hh.


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