libcruft-util/matrix.cpp

399 lines
10 KiB
C++
Raw Normal View History

/*
2015-04-13 18:05:28 +10:00
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
2015-04-13 18:05:28 +10:00
* http://www.apache.org/licenses/LICENSE-2.0
2014-08-18 22:14:31 +10:00
*
2015-04-13 18:05:28 +10:00
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
2015-10-30 23:40:13 +11:00
* Copyright 2011-2015 Danny Robson <danny@nerdcruft.net>
*/
2011-05-23 17:18:52 +10:00
#include "matrix.hpp"
#include "debug.hpp"
#include "iterator.hpp"
#include "point.hpp"
2011-05-23 17:18:52 +10:00
#include <cstring>
2014-12-15 13:44:33 +11:00
#include <cmath>
2011-05-23 17:18:52 +10:00
2011-08-15 20:11:42 +10:00
using namespace util;
2011-05-23 17:18:52 +10:00
///////////////////////////////////////////////////////////////////////////////
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix<S,T>
matrix<S,T>::transposed (void) const
2014-12-15 13:42:44 +11:00
{
2015-10-30 23:40:13 +11:00
matrix<S,T> m;
2015-10-30 23:40:13 +11:00
for (size_t i = 0; i < S; ++i)
for (size_t j = 0; j < S; ++j)
2014-12-15 13:42:44 +11:00
m.values[i][j] = values[j][i];
2014-12-15 13:42:44 +11:00
return m;
}
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix<S,T>&
matrix<S,T>::transpose (void)
2014-12-15 13:42:44 +11:00
{
2015-10-30 23:40:13 +11:00
for (size_t i = 0; i < S; ++i)
for (size_t j = i + 1; j < S; ++j)
2014-12-15 13:42:44 +11:00
std::swap (values[i][j], values[j][i]);
return *this;
}
2014-07-07 15:17:45 +10:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix<S,T>&
matrix<S,T>::invert (void)
{
return *this = inverse ();
2014-08-19 20:45:28 +10:00
}
//-----------------------------------------------------------------------------
//template <size_t S, typename T>
//matrix<S,T>&
//matrix<S,T>::invert_affine (void)
//{
// CHECK (is_affine ());
//
// // inv ([ M b ] == [ inv(M) -inv(M).b ]
// // [ 0 1 ]) [ 0 1 ]
//
// // Invert the 3x3 M
// T A = (values[1][1] * values[2][2] - values[1][2] * values[2][1]);
// T B = (values[1][2] * values[2][0] - values[1][0] * values[2][2]);
// T C = (values[1][0] * values[2][1] - values[1][1] * values[2][0]);
//
// T D = (values[0][2] * values[2][1] - values[0][1] * values[2][2]);
// T E = (values[0][0] * values[2][2] - values[0][2] * values[2][0]);
// T F = (values[2][0] * values[0][1] - values[0][0] * values[2][1]);
//
// T G = (values[0][1] * values[1][2] - values[0][2] * values[1][1]);
// T H = (values[0][2] * values[1][0] - values[0][0] * values[1][2]);
// T K = (values[0][0] * values[1][1] - values[0][1] * values[1][0]);
//
// T d = values[0][0] * A + values[0][1] * B + values[0][2] * C;
// CHECK_NEQ (d, 0.0);
//
// values[0][0] = A / d;
// values[0][1] = D / d;
// values[0][2] = G / d;
// values[1][0] = B / d;
// values[1][1] = E / d;
// values[1][2] = H / d;
// values[2][0] = C / d;
// values[2][1] = F / d;
// values[2][2] = K / d;
//
// // Multiply the b
// T b0 = - values[0][0] * values[0][3] - values[0][1] * values[1][3] - values[0][2] * values[2][3];
// T b1 = - values[1][0] * values[0][3] - values[1][1] * values[1][3] - values[1][2] * values[2][3];
// T b2 = - values[2][0] * values[0][3] - values[2][1] * values[1][3] - values[2][2] * values[2][3];
//
// values[0][3] = b0;
// values[1][3] = b1;
// values[2][3] = b2;
//
// return *this;
//}
2014-08-19 20:45:28 +10:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
T
util::matrix<S,T>::determinant (void) const
{
return util::determinant (*this);
2011-05-23 17:18:52 +10:00
}
2014-08-19 20:45:28 +10:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
util::matrix<S,T>
util::matrix<S,T>::inverse (void) const
{
return util::inverse (*this);
2014-08-19 20:45:28 +10:00
}
2014-07-07 15:17:45 +10:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix<S,T>
matrix<S,T>::operator* (const matrix<S,T> &rhs) const
{
2015-10-30 23:40:13 +11:00
matrix<S,T> m;
2011-05-23 17:18:52 +10:00
2015-10-30 23:40:13 +11:00
for (unsigned row = 0; row < S; ++row) {
for (unsigned col = 0; col < S; ++col) {
m.values[row][col] = T {0};
2015-10-30 23:40:13 +11:00
for (unsigned inner = 0; inner < S; ++inner)
m.values[row][col] += values[row][inner] * rhs.values[inner][col];
}
}
2011-05-23 17:18:52 +10:00
return m;
2011-05-23 17:18:52 +10:00
}
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix<S,T>&
matrix<S,T>::operator*=(const matrix<S,T> &rhs)
{
return *this = *this * rhs;
}
///////////////////////////////////////////////////////////////////////////////
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
vector<S,T>
matrix<S,T>::operator* (const vector<S,T> &rhs) const
{
vector<S,T> out;
for (size_t i = 0; i < S; ++i)
out[i] = dot (rhs, values[i]);
return out;
2014-08-19 20:46:00 +10:00
}
2015-01-13 18:37:53 +11:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
point<S,T>
matrix<S,T>::operator* (const point<S,T> &rhs) const
2015-01-13 18:37:53 +11:00
{
point<S,T> out;
for (size_t i = 0; i < S; ++i)
out[i] = dot (rhs, values[i]);
return out;
2015-01-13 18:37:53 +11:00
}
2014-07-07 15:17:45 +10:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
2011-05-23 17:18:52 +10:00
bool
matrix<S,T>::is_affine (void) const
{
for (size_t i = 0; i < S - 1; ++i)
if (!exactly_zero (values[S-1][i]))
return false;
return exactly_equal (values[S-1][S-1], T{1});
2011-05-23 17:18:52 +10:00
}
2014-12-15 13:44:33 +11:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix4<T>
matrix<S,T>::ortho (T left, T right,
T bottom, T top,
T near, T far)
2014-12-15 13:44:33 +11:00
{
CHECK_GT (far, near);
T tx = - (right + left) / (right - left);
T ty = - (top + bottom) / (top - bottom);
T tz = - (far + near) / (far - near);
T rl = 2 / (right - left);
T tb = 2 / (top - bottom);
T fn = 2 / (far - near);
return { {
{ rl, 0, 0, tx },
{ 0, tb, 0, ty },
{ 0, 0, fn, tz },
{ 0, 0, 0, 1 },
} };
}
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix4<T>
matrix<S,T>::ortho2D (T left, T right,
T bottom, T top)
2014-12-15 13:44:33 +11:00
{
return ortho (left, right, bottom, top, -1, 1);
}
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix4<T>
matrix<S,T>::perspective (T fov, T aspect, range<T> Z)
2014-12-15 13:44:33 +11:00
{
T f = std::tan (fov / 2);
T tx = 1 / (f * aspect);
T ty = 1 / f;
T z1 = (Z.max + Z.min) / (Z.min - Z.max);
T z2 = (2 * Z.max * Z.min) / (Z.min - Z.max);
2014-12-15 13:44:33 +11:00
return { {
{ tx, 0, 0, 0 },
{ 0, ty, 0, 0 },
{ 0, 0, z1, z2 },
{ 0, 0, -1, 0 }
} };
}
//-----------------------------------------------------------------------------
// Emulates gluLookAt
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix4<T>
matrix<S,T>::look_at (util::point<3,T> eye,
util::point<3,T> centre,
util::vector<3,T> up)
2014-12-15 13:44:33 +11:00
{
const auto f = normalised (centre - eye);
const auto s = normalised (cross (f, up));
2014-12-15 13:44:33 +11:00
const auto u = cross (s, f);
return { {
{ s.x, s.y, s.z, -dot (s, eye) },
{ u.x, u.y, u.z, -dot (u, eye) },
{ -f.x, -f.y, -f.z, dot (f, eye) },
{ 0, 0, 0, 1 },
} };
}
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix4<T>
matrix<S,T>::translation (util::vector<2,T> v)
{
return translation ({v.x, v.y, 0});
}
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix4<T>
matrix<S,T>::translation (util::vector<3,T> v)
{
return { {
{ 1.f, 0.f, 0.f, v.x },
{ 0.f, 1.f, 0.f, v.y },
{ 0.f, 0.f, 1.f, v.z },
{ 0.f, 0.f, 0.f, 1.f },
} };
}
2015-01-19 19:12:44 +11:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix4<T>
matrix<S,T>::scale (T mag)
2015-01-19 19:12:44 +11:00
{
return scale (vector<3,T> (mag));
}
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix4<T>
matrix<S,T>::scale (util::vector<3,T> v)
{
return { {
{ v.x, 0.f, 0.f, 0.f },
{ 0.f, v.y, 0.f, 0.f },
{ 0.f, 0.f, v.z, 0.f },
{ 0.f, 0.f, 0.f, 1.f }
} };
}
2014-12-30 01:32:02 +11:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
matrix4<T>
matrix<S,T>::rotation (T angle, util::vector<3,T> about)
2014-12-30 01:32:02 +11:00
{
CHECK (is_normalised (about));
2014-12-30 01:32:02 +11:00
T c = std::cos (angle);
T s = std::sin (angle);
T x = about.x,
y = about.y,
z = about.z;
return { {
{ x * x * (1 - c) + c,
x * y * (1 - c) - z * s,
x * z * (1 - c) + y * s,
0
},
{ y * x * (1 - c) + z * s,
y * y * (1 - c) + c,
y * z * (1 - c) - x * s,
0
},
{ z * x * (1 - c) - y * s,
z * y * (1 - c) + x * s,
z * z * (1 - c) + c,
0
},
{ 0, 0, 0, 1 }
} };
}
2014-07-07 15:17:45 +10:00
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template <size_t S, typename T>
const matrix<S,T>
matrix<S,T>::ZEROES { 0 };
2011-05-23 17:18:52 +10:00
2014-07-07 15:17:45 +10:00
//-----------------------------------------------------------------------------
template struct util::matrix<2,float>;
template struct util::matrix<2,double>;
template struct util::matrix<3,float>;
template struct util::matrix<3,double>;
template struct util::matrix<4,float>;
template struct util::matrix<4,double>;
2011-05-23 17:18:52 +10:00
///////////////////////////////////////////////////////////////////////////////
template <size_t S, typename T>
std::ostream&
util::operator<< (std::ostream &os, const matrix<S,T> &m)
{
os << "{ ";
for (size_t i = 0; i < S; ++i) {
os << "{ ";
std::copy (m[i], m[i]+S, util::infix_iterator<float> (os, ", "));
os << ((i == S - 1) ? " }" : " }, ");
}
return os << " }";
2011-05-23 17:18:52 +10:00
}
//-----------------------------------------------------------------------------
2015-10-30 23:40:13 +11:00
template std::ostream& util::operator<< (std::ostream&, const matrix<4,float>&);
template std::ostream& util::operator<< (std::ostream&, const matrix<4,double>&);