Bullet Collision Detection & Physics Library
btPolarDecomposition.cpp
Go to the documentation of this file.
1 #include "btPolarDecomposition.h"
2 #include "btMinMax.h"
3 
4 namespace
5 {
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 
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 
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 
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 
36 : m_tolerance(tolerance)
37 , m_maxIterations(maxIterations)
38 {
39 }
40 
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 
89 {
90  return m_maxIterations;
91 }
92 
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 
btMatrix3x3 inverse() const
Return the inverse of the matrix.
Definition: btMatrix3x3.h:1075
btPolarDecomposition(btScalar tolerance=btScalar(0.0001), unsigned int maxIterations=16)
Creates an instance with optional parameters.
unsigned int maxIterations() const
Returns the maximum number of iterations that this algorithm will perform to achieve convergence...
btScalar btPow(btScalar x, btScalar y)
Definition: btScalar.h:499
unsigned int polarDecompose(const btMatrix3x3 &a, btMatrix3x3 &u, btMatrix3x3 &h)
This functions decomposes the matrix &#39;a&#39; into two parts: an orthogonal matrix &#39;u&#39; and a symmetric...
bool btFuzzyZero(btScalar x)
Definition: btScalar.h:550
btMatrix3x3 transpose() const
Return the transpose of the matrix.
Definition: btMatrix3x3.h:1030
const T & btMax(const T &a, const T &b)
Definition: btMinMax.h:29
The btMatrix3x3 class implements a 3x3 rotation matrix, to perform linear algebra in combination with...
Definition: btMatrix3x3.h:48
This class is used to compute the polar decomposition of a matrix.
unsigned int decompose(const btMatrix3x3 &a, btMatrix3x3 &u, btMatrix3x3 &h) const
Decomposes a matrix into orthogonal and symmetric, positive-definite parts.
int maxIterations
float btScalar
The btScalar type abstracts floating point numbers, to easily switch between double and single floati...
Definition: btScalar.h:292
btScalar btFabs(btScalar x)
Definition: btScalar.h:475