LibGfx/Matrix: Add inverse() and friends

Matrix inversion comes in quite handy in 3D projections, so let's add
`Matrix<N,T>.inverse()`. To support matrix inversion, the following
methods are added:

* `Matrix.first_minor()`
  See: https://en.wikipedia.org/wiki/Minor_(linear_algebra)
* `Matrix.adjugate()`
  See: https://en.wikipedia.org/wiki/Adjugate_matrix
* `Matrix.determinant()`
  See: https://en.wikipedia.org/wiki/Determinant
* `Matrix.inverse()`
  See: https://en.wikipedia.org/wiki/Invertible_matrix
* `Matrix.operator/()`
  To support easy matrix division :-)

Code loosely based on an implementation listed here:
https://www.geeksforgeeks.org/adjoint-inverse-matrix/
This commit is contained in:
Jelle Raaijmakers 2021-05-23 21:43:29 +02:00 committed by Linus Groh
parent 06c835f857
commit 22d8778437
Notes: sideshowbarker 2024-07-18 17:28:40 +09:00

View file

@ -75,6 +75,70 @@ public:
return product;
}
constexpr Matrix operator/(T divisor) const
{
Matrix division;
for (size_t i = 0; i < N; ++i) {
for (size_t j = 0; j < N; ++j) {
division.m_elements[i][j] = m_elements[i][j] / divisor;
}
}
return division;
}
constexpr Matrix adjugate() const
{
if constexpr (N == 1)
return Matrix(1);
Matrix adjugate;
for (size_t i = 0; i < N; ++i) {
for (size_t j = 0; j < N; ++j) {
int sign = (i + j) % 2 == 0 ? 1 : -1;
adjugate.m_elements[j][i] = sign * first_minor(i, j);
}
}
return adjugate;
}
constexpr T determinant() const
{
if constexpr (N == 1) {
return m_elements[0][0];
} else {
T result = {};
int sign = 1;
for (size_t j = 0; j < N; ++j) {
result += sign * m_elements[0][j] * first_minor(0, j);
sign *= -1;
}
return result;
}
}
constexpr T first_minor(size_t skip_row, size_t skip_column) const
{
static_assert(N > 1);
VERIFY(skip_row < N);
VERIFY(skip_column < N);
Matrix<N - 1, T> first_minor;
constexpr auto new_size = N - 1;
size_t k = 0;
for (size_t i = 0; i < N; ++i) {
for (size_t j = 0; j < N; ++j) {
if (i == skip_row || j == skip_column)
continue;
first_minor.elements()[k / new_size][k % new_size] = m_elements[i][j];
++k;
}
}
return first_minor.determinant();
}
constexpr static Matrix identity()
{
Matrix result;
@ -89,6 +153,13 @@ public:
return result;
}
constexpr Matrix inverse() const
{
auto det = determinant();
VERIFY(det != 0);
return adjugate() / det;
}
constexpr Matrix transpose() const
{
Matrix result;