1 #include "btPolarDecomposition.h"
2 #include "btMinMax.h"
3
4 namespace
5 {
abs_column_sum(const btMatrix3x3 & a,int i)6 btScalar abs_column_sum(const btMatrix3x3& a, int i)
7 {
8 return btFabs(a[0][i]) + btFabs(a[1][i]) + btFabs(a[2][i]);
9 }
10
abs_row_sum(const btMatrix3x3 & a,int i)11 btScalar abs_row_sum(const btMatrix3x3& a, int i)
12 {
13 return btFabs(a[i][0]) + btFabs(a[i][1]) + btFabs(a[i][2]);
14 }
15
p1_norm(const btMatrix3x3 & a)16 btScalar p1_norm(const btMatrix3x3& a)
17 {
18 const btScalar sum0 = abs_column_sum(a,0);
19 const btScalar sum1 = abs_column_sum(a,1);
20 const btScalar sum2 = abs_column_sum(a,2);
21 return btMax(btMax(sum0, sum1), sum2);
22 }
23
pinf_norm(const btMatrix3x3 & a)24 btScalar pinf_norm(const btMatrix3x3& a)
25 {
26 const btScalar sum0 = abs_row_sum(a,0);
27 const btScalar sum1 = abs_row_sum(a,1);
28 const btScalar sum2 = abs_row_sum(a,2);
29 return btMax(btMax(sum0, sum1), sum2);
30 }
31 }
32
33
34
btPolarDecomposition(btScalar tolerance,unsigned int maxIterations)35 btPolarDecomposition::btPolarDecomposition(btScalar tolerance, unsigned int maxIterations)
36 : m_tolerance(tolerance)
37 , m_maxIterations(maxIterations)
38 {
39 }
40
decompose(const btMatrix3x3 & a,btMatrix3x3 & u,btMatrix3x3 & h) const41 unsigned int btPolarDecomposition::decompose(const btMatrix3x3& a, btMatrix3x3& u, btMatrix3x3& h) const
42 {
43 // Use the 'u' and 'h' matrices for intermediate calculations
44 u = a;
45 h = a.inverse();
46
47 for (unsigned int i = 0; i < m_maxIterations; ++i)
48 {
49 const btScalar h_1 = p1_norm(h);
50 const btScalar h_inf = pinf_norm(h);
51 const btScalar u_1 = p1_norm(u);
52 const btScalar u_inf = pinf_norm(u);
53
54 const btScalar h_norm = h_1 * h_inf;
55 const btScalar u_norm = u_1 * u_inf;
56
57 // The matrix is effectively singular so we cannot invert it
58 if (btFuzzyZero(h_norm) || btFuzzyZero(u_norm))
59 break;
60
61 const btScalar gamma = btPow(h_norm / u_norm, 0.25f);
62 const btScalar inv_gamma = btScalar(1.0) / gamma;
63
64 // Determine the delta to 'u'
65 const btMatrix3x3 delta = (u * (gamma - btScalar(2.0)) + h.transpose() * inv_gamma) * btScalar(0.5);
66
67 // Update the matrices
68 u += delta;
69 h = u.inverse();
70
71 // Check for convergence
72 if (p1_norm(delta) <= m_tolerance * u_1)
73 {
74 h = u.transpose() * a;
75 h = (h + h.transpose()) * 0.5;
76 return i;
77 }
78 }
79
80 // The algorithm has failed to converge to the specified tolerance, but we
81 // want to make sure that the matrices returned are in the right form.
82 h = u.transpose() * a;
83 h = (h + h.transpose()) * 0.5;
84
85 return m_maxIterations;
86 }
87
maxIterations() const88 unsigned int btPolarDecomposition::maxIterations() const
89 {
90 return m_maxIterations;
91 }
92
polarDecompose(const btMatrix3x3 & a,btMatrix3x3 & u,btMatrix3x3 & h)93 unsigned int polarDecompose(const btMatrix3x3& a, btMatrix3x3& u, btMatrix3x3& h)
94 {
95 static btPolarDecomposition polar;
96 return polar.decompose(a, u, h);
97 }
98
99