lolengine
/
lol
mirror of https://github.com/lolengine/lol
1
0
Fork 0
Code Issues 0 Releases 0 Wiki Activity
Browse Source

math: add inversion code for 2×2 and 3×3 matrices, and transposition 

code for all matrices.
legacy
Sam Hocevar sam 12 years ago
parent
ad0e3f8b6d
commit
eb51928415
7 changed files with 306 additions and 106 deletions
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. +23
    -18
      src/lol/math/vector.h
  2. +109
    -6
      src/math/vector.cpp
  3. +2
    -2
      test/Makefile.am
  4. +2
    -2
      test/benchmark/vector.cpp
  5. +167
    -0
      test/unit/matrix.cpp
  6. +1
    -77
      test/unit/vector.cpp
  7. +2
    -1
      win32/testsuite.vcxproj

+ 23
- 18
src/lol/math/vector.h View File

@@ -1354,9 +1354,6 @@ template <typename T> struct Mat2
inline Vec2<T>& operator[](size_t n) { return (&v0)[n]; }
inline Vec2<T> const& operator[](size_t n) const { return (&v0)[n]; }

T det() const;
Mat2<T> invert() const;

/* Helpers for transformation matrices */
static Mat2<T> rotate(T angle);

@@ -1436,9 +1433,6 @@ template <typename T> struct Mat3
inline Vec3<T>& operator[](size_t n) { return (&v0)[n]; }
inline Vec3<T> const& operator[](size_t n) const { return (&v0)[n]; }

T det() const;
Mat3<T> invert() const;

/* Helpers for transformation matrices */
static Mat3<T> rotate(T angle, T x, T y, T z);
static Mat3<T> rotate(T angle, Vec3<T> v);
@@ -1521,9 +1515,6 @@ template <typename T> struct Mat4
inline Vec4<T>& operator[](size_t n) { return (&v0)[n]; }
inline Vec4<T> const& operator[](size_t n) const { return (&v0)[n]; }

T det() const;
Mat4<T> invert() const;

/* Helpers for transformation matrices */
static Mat4<T> translate(T x, T y, T z);
static Mat4<T> translate(Vec3<T> v);
@@ -1531,16 +1522,18 @@ template <typename T> struct Mat4
static Mat4<T> rotate(T angle, Vec3<T> v);
static Mat4<T> rotate(Quat<T> q);

static inline Mat4<T> translate(Mat4<T> mat, Vec3<T> v)
static inline Mat4<T> translate(Mat4<T> const &mat, Vec3<T> v)
{
return translate(v) * mat;
}

static inline Mat4<T> rotate(Mat4<T> mat, T angle, Vec3<T> v)
static inline Mat4<T> rotate(Mat4<T> &mat, T angle, Vec3<T> v)
{
return rotate(angle, v) * mat;
}

static Mat3<T> normal(Mat4<T> const &mat);

/* Helpers for view matrices */
static Mat4<T> lookat(Vec3<T> eye, Vec3<T> center, Vec3<T> up);

@@ -1556,37 +1549,37 @@ template <typename T> struct Mat4
friend std::ostream &operator<<(std::ostream &stream, Mat4<U> const &m);
#endif

inline Mat4<T> operator +(Mat4<T> const m) const
inline Mat4<T> operator +(Mat4<T> const &m) const
{
return Mat4<T>(v0 + m[0], v1 + m[1], v2 + m[2], v3 + m[3]);
}

inline Mat4<T> operator +=(Mat4<T> const m)
inline Mat4<T> operator +=(Mat4<T> const &m)
{
return *this = *this + m;
}

inline Mat4<T> operator -(Mat4<T> const m) const
inline Mat4<T> operator -(Mat4<T> const &m) const
{
return Mat4<T>(v0 - m[0], v1 - m[1], v2 - m[2], v3 - m[3]);
}

inline Mat4<T> operator -=(Mat4<T> const m)
inline Mat4<T> operator -=(Mat4<T> const &m)
{
return *this = *this - m;
}

inline Mat4<T> operator *(Mat4<T> const m) const
inline Mat4<T> operator *(Mat4<T> const &m) const
{
return Mat4<T>(*this * m[0], *this * m[1], *this * m[2], *this * m[3]);
}

inline Mat4<T> operator *=(Mat4<T> const m)
inline Mat4<T> operator *=(Mat4<T> const &m)
{
return *this = *this * m;
}

inline Vec4<T> operator *(Vec4<T> const m) const
inline Vec4<T> operator *(Vec4<T> const &m) const
{
Vec4<T> ret;
for (int j = 0; j < 4; j++)
@@ -1602,6 +1595,18 @@ template <typename T> struct Mat4
Vec4<T> v0, v1, v2, v3;
};

template<typename T> T determinant(Mat2<T> const &);
template<typename T> T determinant(Mat3<T> const &);
template<typename T> T determinant(Mat4<T> const &);

template<typename T> Mat2<T> transpose(Mat2<T> const &);
template<typename T> Mat3<T> transpose(Mat3<T> const &);
template<typename T> Mat4<T> transpose(Mat4<T> const &);

template<typename T> Mat2<T> inverse(Mat2<T> const &);
template<typename T> Mat3<T> inverse(Mat3<T> const &);
template<typename T> Mat4<T> inverse(Mat4<T> const &);

/*
* Arbitrarily-sized square matrices; for now this only supports
* naive inversion and is used for the Remez inversion method.


+ 109
- 6
src/math/vector.cpp View File

@@ -36,6 +36,12 @@ using namespace std;
namespace lol
{

static inline float det2(float a, float b,
float c, float d)
{
return a * d - b * c;
}

static inline float det3(float a, float b, float c,
float d, float e, float f,
float g, float h, float i)
@@ -45,7 +51,20 @@ static inline float det3(float a, float b, float c,
+ c * (d * h - g * e);
}

static inline float cofact3(mat4 const &mat, int i, int j)
static inline float cofact(mat2 const &mat, int i, int j)
{
return mat[(i + 1) & 1][(j + 1) & 1] * (((i + j) & 1) ? -1.0f : 1.0f);
}

static inline float cofact(mat3 const &mat, int i, int j)
{
return det2(mat[(i + 1) % 3][(j + 1) % 3],
mat[(i + 2) % 3][(j + 1) % 3],
mat[(i + 1) % 3][(j + 2) % 3],
mat[(i + 2) % 3][(j + 2) % 3]);
}

static inline float cofact(mat4 const &mat, int i, int j)
{
return det3(mat[(i + 1) & 3][(j + 1) & 3],
mat[(i + 2) & 3][(j + 1) & 3],
@@ -58,24 +77,91 @@ static inline float cofact3(mat4 const &mat, int i, int j)
mat[(i + 3) & 3][(j + 3) & 3]) * (((i + j) & 1) ? -1.0f : 1.0f);
}

template<> float mat4::det() const
template<> float determinant(mat2 const &mat)
{
return mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
}

template<> mat2 transpose(mat2 const &mat)
{
mat2 ret;
for (int j = 0; j < 2; j++)
for (int i = 0; i < 2; i++)
ret[j][i] = mat[i][j];
return ret;
}

template<> mat2 inverse(mat2 const &mat)
{
mat2 ret;
float d = determinant(mat);
if (d)
{
d = 1.0f / d;
for (int j = 0; j < 2; j++)
for (int i = 0; i < 2; i++)
ret[j][i] = cofact(mat, i, j) * d;
}
return ret;
}

template<> float determinant(mat3 const &mat)
{
return det3(mat[0][0], mat[0][1], mat[0][2],
mat[1][0], mat[1][1], mat[1][2],
mat[2][0], mat[2][1], mat[2][2]);
}

template<> mat3 transpose(mat3 const &mat)
{
mat3 ret;
for (int j = 0; j < 3; j++)
for (int i = 0; i < 3; i++)
ret[j][i] = mat[i][j];
return ret;
}

template<> mat3 inverse(mat3 const &mat)
{
mat3 ret;
float d = determinant(mat);
if (d)
{
d = 1.0f / d;
for (int j = 0; j < 3; j++)
for (int i = 0; i < 3; i++)
ret[j][i] = cofact(mat, i, j) * d;
}
return ret;
}

template<> float determinant(mat4 const &mat)
{
float ret = 0;
for (int n = 0; n < 4; n++)
ret += (*this)[n][0] * cofact3(*this, n, 0);
ret += mat[n][0] * cofact(mat, n, 0);
return ret;
}

template<> mat4 mat4::invert() const
template<> mat4 transpose(mat4 const &mat)
{
mat4 ret;
float d = det();
for (int j = 0; j < 4; j++)
for (int i = 0; i < 4; i++)
ret[j][i] = mat[i][j];
return ret;
}

template<> mat4 inverse(mat4 const &mat)
{
mat4 ret;
float d = determinant(mat);
if (d)
{
d = 1.0f / d;
for (int j = 0; j < 4; j++)
for (int i = 0; i < 4; i++)
ret[j][i] = cofact3(*this, i, j) * d;
ret[j][i] = cofact(mat, i, j) * d;
}
return ret;
}
@@ -120,6 +206,23 @@ template<> void quat::printf() const
Log::Debug("[ %6.6f %6.6f %6.6f %6.6f ]\n", x, y, z, w);
}

template<> void mat2::printf() const
{
mat2 const &p = *this;

Log::Debug("[ %6.6f %6.6f\n", p[0][0], p[1][0]);
Log::Debug(" %6.6f %6.6f ]\n", p[0][1], p[1][1]);
}

template<> void mat3::printf() const
{
mat3 const &p = *this;

Log::Debug("[ %6.6f %6.6f %6.6f\n", p[0][0], p[1][0], p[2][0]);
Log::Debug(" %6.6f %6.6f %6.6f\n", p[0][1], p[1][1], p[2][1]);
Log::Debug(" %6.6f %6.6f %6.6f ]\n", p[0][2], p[1][2], p[2][2]);
}

template<> void mat4::printf() const
{
mat4 const &p = *this;


+ 2
- 2
test/Makefile.am View File

@@ -22,8 +22,8 @@ noinst_PROGRAMS = quad benchsuite testsuite
TESTS = testsuite

testsuite_SOURCES = testsuite.cpp \
unit/vector.cpp unit/half.cpp unit/trig.cpp unit/build.cpp \
unit/real.cpp unit/image.cpp unit/quat.cpp unit/cmplx.cpp
unit/vector.cpp unit/matrix.cpp unit/half.cpp unit/trig.cpp \
unit/build.cpp unit/real.cpp unit/image.cpp unit/quat.cpp unit/cmplx.cpp
testsuite_CPPFLAGS = @LOL_CFLAGS@ @PIPI_CFLAGS@
testsuite_LDFLAGS = $(top_builddir)/src/liblol.a @LOL_LIBS@ @PIPI_LIBS@
testsuite_DEPENDENCIES = $(top_builddir)/src/liblol.a


+ 2
- 2
test/benchmark/vector.cpp View File

@@ -52,7 +52,7 @@ void bench_matrix(int mode)
/* Determinant */
timer.GetMs();
for (size_t i = 0; i < MATRIX_TABLE_SIZE; i++)
pf[i] = pm[i].det();
pf[i] = determinant(pm[i]);
result[1] += timer.GetMs();

/* Multiply matrices */
@@ -70,7 +70,7 @@ void bench_matrix(int mode)
/* Invert matrix */
timer.GetMs();
for (size_t i = 0; i < MATRIX_TABLE_SIZE; i++)
pm[i] = pm[i].invert();
pm[i] = inverse(pm[i]);
result[4] += timer.GetMs();
}



+ 167
- 0
test/unit/matrix.cpp View File

@@ -0,0 +1,167 @@
//
// Lol Engine
//
// Copyright: (c) 2010-2012 Sam Hocevar <sam@hocevar.net>
// This program is free software; you can redistribute it and/or
// modify it under the terms of the Do What The Fuck You Want To
// Public License, Version 2, as published by Sam Hocevar. See
// http://sam.zoy.org/projects/COPYING.WTFPL for more details.
//

#if defined HAVE_CONFIG_H
# include "config.h"
#endif

#include "core.h"
#include "lol/unit.h"

namespace lol
{

LOLUNIT_FIXTURE(MatrixTest)
{
void SetUp()
{
id2 = mat2(1.0f);
tri2 = mat2(vec2(1.0f, 0.0f),
vec2(7.0f, 2.0f));
inv2 = mat2(vec2(4.0f, 3.0f),
vec2(3.0f, 2.0f));

id3 = mat3(1.0f);
tri3 = mat3(vec3(1.0f, 0.0f, 0.0f),
vec3(7.0f, 2.0f, 0.0f),
vec3(1.0f, 5.0f, 3.0f));
inv3 = mat3(vec3(2.0f, 3.0f, 5.0f),
vec3(3.0f, 2.0f, 3.0f),
vec3(9.0f, 5.0f, 7.0f));

id4 = mat4(1.0f);
tri4 = mat4(vec4(1.0f, 0.0f, 0.0f, 0.0f),
vec4(7.0f, 2.0f, 0.0f, 0.0f),
vec4(1.0f, 5.0f, 3.0f, 0.0f),
vec4(8.0f, 9.0f, 2.0f, 4.0f));
inv4 = mat4(vec4( 1.0f, 1.0f, 2.0f, -1.0f),
vec4(-2.0f, -1.0f, -2.0f, 2.0f),
vec4( 4.0f, 2.0f, 5.0f, -4.0f),
vec4( 5.0f, -3.0f, -7.0f, -6.0f));
}

void TearDown() {}

LOLUNIT_TEST(Determinant)
{
float d1, d2;

d1 = determinant(tri2);
LOLUNIT_ASSERT_EQUAL(d1, 2.0f);
d2 = determinant(inv2);
LOLUNIT_ASSERT_EQUAL(d2, -1.0f);

d1 = determinant(tri3);
LOLUNIT_ASSERT_EQUAL(d1, 6.0f);
d2 = determinant(inv3);
LOLUNIT_ASSERT_EQUAL(d2, 1.0f);

d1 = determinant(tri4);
LOLUNIT_ASSERT_EQUAL(d1, 24.0f);
d2 = determinant(inv4);
LOLUNIT_ASSERT_EQUAL(d2, -1.0f);
}

LOLUNIT_TEST(Multiplication)
{
mat4 m0 = id4;
mat4 m1 = id4;
mat4 m2 = m0 * m1;

LOLUNIT_ASSERT_EQUAL(m2[0][0], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][0], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][0], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][0], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][1], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][1], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][2], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][2], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][3], 1.0f);
}

LOLUNIT_TEST(Inverse2x2)
{
mat2 m0 = inv2;
mat2 m1 = inverse(m0);

mat2 m2 = m0 * m1;

LOLUNIT_ASSERT_EQUAL(m2[0][0], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][0], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][1], 1.0f);
}

LOLUNIT_TEST(Inverse3x3)
{
mat3 m0 = inv3;
mat3 m1 = inverse(m0);

mat3 m2 = m0 * m1;

LOLUNIT_ASSERT_EQUAL(m2[0][0], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][0], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][0], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][1], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][1], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][2], 1.0f);
}

LOLUNIT_TEST(Inverse4x4)
{
mat4 m0 = inv4;
mat4 m1 = inverse(m0);

mat4 m2 = m0 * m1;

LOLUNIT_ASSERT_EQUAL(m2[0][0], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][0], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][0], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][0], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][1], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][1], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][2], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][2], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][3], 1.0f);
}

mat2 tri2, id2, inv2;
mat3 tri3, id3, inv3;
mat4 tri4, id4, inv4;
};

} /* namespace lol */


+ 1
- 77
test/unit/vector.cpp View File

@@ -20,18 +20,7 @@ namespace lol

LOLUNIT_FIXTURE(VectorTest)
{
void SetUp()
{
identity = mat4(1.0f);
triangular = mat4(vec4(1.0f, 0.0f, 0.0f, 0.0f),
vec4(7.0f, 2.0f, 0.0f, 0.0f),
vec4(1.0f, 5.0f, 3.0f, 0.0f),
vec4(8.0f, 9.0f, 2.0f, 4.0f));
invertible = mat4(vec4( 1.0f, 1.0f, 2.0f, -1.0f),
vec4(-2.0f, -1.0f, -2.0f, 2.0f),
vec4( 4.0f, 2.0f, 5.0f, -4.0f),
vec4( 5.0f, -3.0f, -7.0f, -6.0f));
}
void SetUp() {}

void TearDown() {}

@@ -135,71 +124,6 @@ LOLUNIT_FIXTURE(VectorTest)
LOLUNIT_ASSERT_EQUAL(c.w, 0.0f);
LOLUNIT_ASSERT_EQUAL(a3, a1);
}

LOLUNIT_TEST(MatrixDeterminant)
{
float d1 = triangular.det();
LOLUNIT_ASSERT_EQUAL(d1, 24.0f);
float d2 = invertible.det();
LOLUNIT_ASSERT_EQUAL(d2, -1.0f);
}

LOLUNIT_TEST(MatrixMultiplication)
{
mat4 m0 = identity;
mat4 m1 = identity;
mat4 m2 = m0 * m1;

LOLUNIT_ASSERT_EQUAL(m2[0][0], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][0], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][0], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][0], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][1], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][1], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][2], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][2], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][3], 1.0f);
}

LOLUNIT_TEST(MatrixInverse)
{
mat4 m0 = invertible;
mat4 m1 = m0.invert();

mat4 m2 = m0 * m1;

LOLUNIT_ASSERT_EQUAL(m2[0][0], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][0], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][0], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][0], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][1], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][1], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][1], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][2], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][2], 1.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][2], 0.0f);

LOLUNIT_ASSERT_EQUAL(m2[0][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[1][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[2][3], 0.0f);
LOLUNIT_ASSERT_EQUAL(m2[3][3], 1.0f);
}

mat4 triangular, identity, invertible;
};

} /* namespace lol */


+ 2
- 1
win32/testsuite.vcxproj View File

@@ -32,6 +32,7 @@
<ClCompile Include="..\test\unit\cmplx.cpp" />
<ClCompile Include="..\test\unit\half.cpp" />
<ClCompile Include="..\test\unit\image.cpp" />
<ClCompile Include="..\test\unit\matrix.cpp" />
<ClCompile Include="..\test\unit\quat.cpp" />
<ClCompile Include="..\test\unit\real.cpp" />
<ClCompile Include="..\test\unit\trig.cpp" />
@@ -60,4 +61,4 @@
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
<ImportGroup Label="ExtensionTargets">
</ImportGroup>
</Project>
</Project>

Write Preview
Loading…
Cancel
Save