Page MenuHome
Differential D8762 Diff 28333 extern/bullet2/src/LinearMath/btPolarDecomposition.cpp

Changeset View

Standalone View

View Options

extern/bullet2/src/LinearMath/btPolarDecomposition.cpp

#include "btPolarDecomposition.h" #include "btPolarDecomposition.h"
#include "btMinMax.h" #include "btMinMax.h"
namespace namespace
{ {
btScalar abs_column_sum(const btMatrix3x3& a, int i) btScalar abs_column_sum(const btMatrix3x3& a, int i)
{ {
return btFabs(a[0][i]) + btFabs(a[1][i]) + btFabs(a[2][i]); return btFabs(a[0][i]) + btFabs(a[1][i]) + btFabs(a[2][i]);
} }
btScalar abs_row_sum(const btMatrix3x3& a, int i) btScalar abs_row_sum(const btMatrix3x3& a, int i)
{ {
return btFabs(a[i][0]) + btFabs(a[i][1]) + btFabs(a[i][2]); return btFabs(a[i][0]) + btFabs(a[i][1]) + btFabs(a[i][2]);
} }
btScalar p1_norm(const btMatrix3x3& a) btScalar p1_norm(const btMatrix3x3& a)
{ {
const btScalar sum0 = abs_column_sum(a,0); const btScalar sum0 = abs_column_sum(a, 0);
const btScalar sum1 = abs_column_sum(a,1); const btScalar sum1 = abs_column_sum(a, 1);
const btScalar sum2 = abs_column_sum(a,2); const btScalar sum2 = abs_column_sum(a, 2);
return btMax(btMax(sum0, sum1), sum2); return btMax(btMax(sum0, sum1), sum2);
} }
btScalar pinf_norm(const btMatrix3x3& a) btScalar pinf_norm(const btMatrix3x3& a)
{ {
const btScalar sum0 = abs_row_sum(a,0); const btScalar sum0 = abs_row_sum(a, 0);
const btScalar sum1 = abs_row_sum(a,1); const btScalar sum1 = abs_row_sum(a, 1);
const btScalar sum2 = abs_row_sum(a,2); const btScalar sum2 = abs_row_sum(a, 2);
return btMax(btMax(sum0, sum1), sum2); return btMax(btMax(sum0, sum1), sum2);
} }
} } // namespace
const btScalar btPolarDecomposition::DEFAULT_TOLERANCE = btScalar(0.0001);
const unsigned int btPolarDecomposition::DEFAULT_MAX_ITERATIONS = 16;
btPolarDecomposition::btPolarDecomposition(btScalar tolerance, unsigned int maxIterations) btPolarDecomposition::btPolarDecomposition(btScalar tolerance, unsigned int maxIterations)
: m_tolerance(tolerance) : m_tolerance(tolerance), m_maxIterations(maxIterations)
, m_maxIterations(maxIterations)
{ {
} }
unsigned int btPolarDecomposition::decompose(const btMatrix3x3& a, btMatrix3x3& u, btMatrix3x3& h) const unsigned int btPolarDecomposition::decompose(const btMatrix3x3& a, btMatrix3x3& u, btMatrix3x3& h) const
{ {
// Use the 'u' and 'h' matrices for intermediate calculations // Use the 'u' and 'h' matrices for intermediate calculations
u = a; u = a;
h = a.inverse(); h = a.inverse();
for (unsigned int i = 0; i < m_maxIterations; ++i) for (unsigned int i = 0; i < m_maxIterations; ++i)
{ {
const btScalar h_1 = p1_norm(h); const btScalar h_1 = p1_norm(h);
const btScalar h_inf = pinf_norm(h); const btScalar h_inf = pinf_norm(h);
const btScalar u_1 = p1_norm(u); const btScalar u_1 = p1_norm(u);
const btScalar u_inf = pinf_norm(u); const btScalar u_inf = pinf_norm(u);
const btScalar h_norm = h_1 * h_inf; const btScalar h_norm = h_1 * h_inf;
const btScalar u_norm = u_1 * u_inf; const btScalar u_norm = u_1 * u_inf;
// The matrix is effectively singular so we cannot invert it // The matrix is effectively singular so we cannot invert it
if (btFuzzyZero(h_norm) || btFuzzyZero(u_norm)) if (btFuzzyZero(h_norm) || btFuzzyZero(u_norm))
break; break;
const btScalar gamma = btPow(h_norm / u_norm, 0.25f); const btScalar gamma = btPow(h_norm / u_norm, 0.25f);
const btScalar inv_gamma = btScalar(1.0) / gamma; const btScalar inv_gamma = btScalar(1.0) / gamma;
// Determine the delta to 'u' // Determine the delta to 'u'
const btMatrix3x3 delta = (u * (gamma - btScalar(2.0)) + h.transpose() * inv_gamma) * btScalar(0.5); const btMatrix3x3 delta = (u * (gamma - btScalar(2.0)) + h.transpose() * inv_gamma) * btScalar(0.5);
// Update the matrices // Update the matrices
u += delta; u += delta;
h = u.inverse(); h = u.inverse();
// Check for convergence // Check for convergence
if (p1_norm(delta) <= m_tolerance * u_1) if (p1_norm(delta) <= m_tolerance * u_1)
{ {
h = u.transpose() * a; h = u.transpose() * a;
h = (h + h.transpose()) * 0.5; h = (h + h.transpose()) * 0.5;
return i; return i;
} }
} }
// The algorithm has failed to converge to the specified tolerance, but we // The algorithm has failed to converge to the specified tolerance, but we
// want to make sure that the matrices returned are in the right form. // want to make sure that the matrices returned are in the right form.
h = u.transpose() * a; h = u.transpose() * a;
h = (h + h.transpose()) * 0.5; h = (h + h.transpose()) * 0.5;
return m_maxIterations; return m_maxIterations;
} }
unsigned int btPolarDecomposition::maxIterations() const unsigned int btPolarDecomposition::maxIterations() const
{ {
return m_maxIterations; return m_maxIterations;
} }
unsigned int polarDecompose(const btMatrix3x3& a, btMatrix3x3& u, btMatrix3x3& h) unsigned int polarDecompose(const btMatrix3x3& a, btMatrix3x3& u, btMatrix3x3& h)
{ {
static btPolarDecomposition polar; static btPolarDecomposition polar;
return polar.decompose(a, u, h); return polar.decompose(a, u, h);
} }