Skip to content

Commit

Permalink
Matrix4 now uses an Eigen::Transform internally
Browse files Browse the repository at this point in the history 
Eigen::Transform is a wrapper for an Eigen::Matrix structure which adds
functionality specific to geometric transformations, e.g. extracting the affine
part from an arbitrary (projective) transformation. Our own Matrix4 class is
now just interface glue for this Eigen::Transformation object, which means we
no longer need to implement linear algebra ourselves.

Since Matrix4 is well covered by unit tests (which still pass), there should be
no mathematical problems with this change, however the code is not yet optimal
since some functions are still performing maths manually using xx(), xy() etc
when they could be implemented using functionality exposed directly by Eigen
(some of which may be optimised using SIMD instructions).
  • Loading branch information
Matthew Mott committed Mar 27, 2021
1 parent f03e405 commit c137351
Show file tree
Hide file tree
Showing 2 changed files with 67 additions and 126 deletions.
56 changes: 5 additions & 51 deletions libs/math/Matrix4.cpp
View file
Expand Up @@ -252,67 +252,21 @@ Matrix4 Matrix4::getScale(const Vector3& scale)
// Transpose the matrix in-place
void Matrix4::transpose()
{
std::swap(_m[1], _m[4]); // xy <=> yx
std::swap(_m[2], _m[8]); // xz <=> zx
std::swap(_m[3], _m[12]); // xw <=> tx
std::swap(_m[6], _m[9]); // yz <=> zy
std::swap(_m[7], _m[13]); // yw <=> ty
std::swap(_m[11], _m[14]); // zw <=> tz
_transform.matrix().transposeInPlace();
}

// Return transposed copy
Matrix4 Matrix4::getTransposed() const
{
return Matrix4(
xx(), yx(), zx(), tx(),
xy(), yy(), zy(), ty(),
xz(), yz(), zz(), tz(),
xw(), yw(), zw(), tw()
);
Matrix4 copy = *this;
copy.transpose();
return copy;
}

// Return affine inverse
Matrix4 Matrix4::getInverse() const
{
Matrix4 result;

// determinant of rotation submatrix
double det
= _m[0] * ( _m[5]*_m[10] - _m[9]*_m[6] )
- _m[1] * ( _m[4]*_m[10] - _m[8]*_m[6] )
+ _m[2] * ( _m[4]*_m[9] - _m[8]*_m[5] );

// throw exception here if (det*det < 1e-25)

// invert rotation submatrix
det = 1.0f / det;

result[0] = ( (_m[5]*_m[10]- _m[6]*_m[9] )*det);
result[1] = (- (_m[1]*_m[10]- _m[2]*_m[9] )*det);
result[2] = ( (_m[1]*_m[6] - _m[2]*_m[5] )*det);
result[3] = 0;
result[4] = (- (_m[4]*_m[10]- _m[6]*_m[8] )*det);
result[5] = ( (_m[0]*_m[10]- _m[2]*_m[8] )*det);
result[6] = (- (_m[0]*_m[6] - _m[2]*_m[4] )*det);
result[7] = 0;
result[8] = ( (_m[4]*_m[9] - _m[5]*_m[8] )*det);
result[9] = (- (_m[0]*_m[9] - _m[1]*_m[8] )*det);
result[10]= ( (_m[0]*_m[5] - _m[1]*_m[4] )*det);
result[11] = 0;

// multiply translation part by rotation
result[12] = - (_m[12] * result[0] +
_m[13] * result[4] +
_m[14] * result[8]);
result[13] = - (_m[12] * result[1] +
_m[13] * result[5] +
_m[14] * result[9]);
result[14] = - (_m[12] * result[2] +
_m[13] * result[6] +
_m[14] * result[10]);
result[15] = 1;

return result;
return Matrix4(_transform.inverse(Eigen::Affine));
}

Matrix4 Matrix4::getFullInverse() const
Expand Down
137 changes: 62 additions & 75 deletions libs/math/Matrix4.h
View file
Expand Up @@ -8,7 +8,7 @@
#include "math/pi.h"

#undef Success // get rid of fuckwit X.h macro
#include <Eigen/Dense>
#include <Eigen/Geometry>

class Quaternion;

Expand All @@ -31,10 +31,10 @@ class Quaternion;
*/
class Matrix4
{
// Elements of the 4x4 matrix. These appear to be treated COLUMNWISE, i.e.
// elements [0] through [3] are the first column, [4] through [7] are the
// second column, etc.
double _m[16];
// Underlying Eigen transform object (which can in turn be reduced to a 4x4
// Matrix)
using Transform = Eigen::Projective3d;
Transform _transform;

private:

Expand All @@ -61,6 +61,16 @@ class Matrix4
/// Construct a matrix with uninitialised values.
Matrix4() { }

/// Construct from Eigen transform
explicit Matrix4(const Eigen::Projective3d& t): _transform(t)
{}

/// Get the underlying Eigen transform
Eigen::Projective3d& eigen() { return _transform; }

/// Get the underlying const Eigen transform
const Eigen::Projective3d& eigen() const { return _transform; }

/// Obtain the identity matrix.
static const Matrix4& getIdentity();

Expand Down Expand Up @@ -152,38 +162,41 @@ class Matrix4
* Return matrix elements
* \{
*/
double& xx() { return _m[0]; }
const double& xx() const { return _m[0]; }
double& xy() { return _m[1]; }
const double& xy() const { return _m[1]; }
double& xz() { return _m[2]; }
const double& xz() const { return _m[2]; }
double& xw() { return _m[3]; }
const double& xw() const { return _m[3]; }
double& yx() { return _m[4]; }
const double& yx() const { return _m[4]; }
double& yy() { return _m[5]; }
const double& yy() const { return _m[5]; }
double& yz() { return _m[6]; }
const double& yz() const { return _m[6]; }
double& yw() { return _m[7]; }
const double& yw() const { return _m[7]; }
double& zx() { return _m[8]; }
const double& zx() const { return _m[8]; }
double& zy() { return _m[9]; }
const double& zy() const { return _m[9]; }
double& zz() { return _m[10]; }
const double& zz() const { return _m[10]; }
double& zw() { return _m[11]; }
const double& zw() const { return _m[11]; }
double& tx() { return _m[12]; }
const double& tx() const { return _m[12]; }
double& ty() { return _m[13]; }
const double& ty() const { return _m[13]; }
double& tz() { return _m[14]; }
const double& tz() const { return _m[14]; }
double& tw() { return _m[15]; }
const double& tw() const { return _m[15]; }
double& xx() { return _transform.matrix()(0, 0); }
const double& xx() const { return _transform.matrix()(0, 0); }
double& xy() { return _transform.matrix()(1, 0); }
const double& xy() const { return _transform.matrix()(1, 0); }
double& xz() { return _transform.matrix()(2, 0); }
const double& xz() const { return _transform.matrix()(2, 0); }
double& xw() { return _transform.matrix()(3, 0); }
const double& xw() const { return _transform.matrix()(3, 0); }

double& yx() { return _transform.matrix()(0, 1); }
const double& yx() const { return _transform.matrix()(0, 1); }
double& yy() { return _transform.matrix()(1, 1); }
const double& yy() const { return _transform.matrix()(1, 1); }
double& yz() { return _transform.matrix()(2, 1); }
const double& yz() const { return _transform.matrix()(2, 1); }
double& yw() { return _transform.matrix()(3, 1); }
const double& yw() const { return _transform.matrix()(3, 1); }

double& zx() { return _transform.matrix()(0, 2); }
const double& zx() const { return _transform.matrix()(0, 2); }
double& zy() { return _transform.matrix()(1, 2); }
const double& zy() const { return _transform.matrix()(1, 2); }
double& zz() { return _transform.matrix()(2, 2); }
const double& zz() const { return _transform.matrix()(2, 2); }
double& zw() { return _transform.matrix()(3, 2); }
const double& zw() const { return _transform.matrix()(3, 2); }

double& tx() { return _transform.matrix()(0, 3); }
const double& tx() const { return _transform.matrix()(0, 3); }
double& ty() { return _transform.matrix()(1, 3); }
const double& ty() const { return _transform.matrix()(1, 3); }
double& tz() { return _transform.matrix()(2, 3); }
const double& tz() const { return _transform.matrix()(2, 3); }
double& tw() { return _transform.matrix()(3, 3); }
const double& tw() const { return _transform.matrix()(3, 3); }
/**
* \}
*/
Expand Down Expand Up @@ -234,15 +247,15 @@ class Matrix4
*/
operator double* ()
{
return _m;
return _transform.matrix().data();
}

/**
* Cast to const double* to provide operator[] for const objects.
*/
operator const double* () const
{
return _m;
return _transform.matrix().data();
}

/// Transpose this matrix in-place.
Expand Down Expand Up @@ -296,11 +309,11 @@ class Matrix4
template<typename Element>
BasicVector4<Element> transform(const BasicVector4<Element>& vector4) const;

/**
* \brief
* Return the result of this matrix post-multiplied by another matrix.
*/
Matrix4 getMultipliedBy(const Matrix4& other) const;
/// Return the result of this matrix post-multiplied by another matrix.
Matrix4 getMultipliedBy(const Matrix4& other) const
{
return Matrix4(_transform * other.eigen());
}

/**
* \brief
Expand Down Expand Up @@ -474,29 +487,6 @@ inline Matrix4 Matrix4::byRows(double xx, double yx, double zx, double tx,
tx, ty, tz, tw);
}

// Post-multiply this with other
inline Matrix4 Matrix4::getMultipliedBy(const Matrix4& other) const
{
return Matrix4::byColumns(
other[0] * _m[0] + other[1] * _m[4] + other[2] * _m[8] + other[3] * _m[12],
other[0] * _m[1] + other[1] * _m[5] + other[2] * _m[9] + other[3] * _m[13],
other[0] * _m[2] + other[1] * _m[6] + other[2] * _m[10]+ other[3] * _m[14],
other[0] * _m[3] + other[1] * _m[7] + other[2] * _m[11]+ other[3] * _m[15],
other[4] * _m[0] + other[5] * _m[4] + other[6] * _m[8] + other[7] * _m[12],
other[4] * _m[1] + other[5] * _m[5] + other[6] * _m[9] + other[7] * _m[13],
other[4] * _m[2] + other[5] * _m[6] + other[6] * _m[10]+ other[7] * _m[14],
other[4] * _m[3] + other[5] * _m[7] + other[6] * _m[11]+ other[7] * _m[15],
other[8] * _m[0] + other[9] * _m[4] + other[10]* _m[8] + other[11]* _m[12],
other[8] * _m[1] + other[9] * _m[5] + other[10]* _m[9] + other[11]* _m[13],
other[8] * _m[2] + other[9] * _m[6] + other[10]* _m[10]+ other[11]* _m[14],
other[8] * _m[3] + other[9] * _m[7] + other[10]* _m[11]+ other[11]* _m[15],
other[12]* _m[0] + other[13]* _m[4] + other[14]* _m[8] + other[15]* _m[12],
other[12]* _m[1] + other[13]* _m[5] + other[14]* _m[9] + other[15]* _m[13],
other[12]* _m[2] + other[13]* _m[6] + other[14]* _m[10]+ other[15]* _m[14],
other[12]* _m[3] + other[13]* _m[7] + other[14]* _m[11]+ other[15]* _m[15]
);
}

/// Compare two matrices elementwise for equality
inline bool operator==(const Matrix4& l, const Matrix4& r)
{
Expand Down Expand Up @@ -553,15 +543,12 @@ BasicVector3<Element> Matrix4::transformDirection(const BasicVector3<Element>& d
);
}

template<typename Element>
BasicVector4<Element> Matrix4::transform(const BasicVector4<Element>& vector4) const
template<typename T>
BasicVector4<T> Matrix4::transform(const BasicVector4<T>& vector4) const
{
return BasicVector4<Element>(
static_cast<Element>(_m[0] * vector4[0] + _m[4] * vector4[1] + _m[8] * vector4[2] + _m[12] * vector4[3]),
static_cast<Element>(_m[1] * vector4[0] + _m[5] * vector4[1] + _m[9] * vector4[2] + _m[13] * vector4[3]),
static_cast<Element>(_m[2] * vector4[0] + _m[6] * vector4[1] + _m[10] * vector4[2] + _m[14] * vector4[3]),
static_cast<Element>(_m[3] * vector4[0] + _m[7] * vector4[1] + _m[11] * vector4[2] + _m[15] * vector4[3])
);
Eigen::Matrix<T, 4, 1> eVec(&vector4.x());
auto result = _transform * eVec;
return BasicVector4<T>(result[0], result[1], result[2], result[3]);
}

inline Vector3 Matrix4::getEulerAnglesXYZ() const
Expand Down

0 comments on commit c137351

Please sign in to comment.