//
// Copyright 2015 The ANGLE Project Authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.
//
// Matrix:
// Utility class implementing various matrix operations.
// Supports matrices with minimum 2 and maximum 4 number of rows/columns.
//
// TODO: Check if we can merge Matrix.h in sample_util with this and replace it with this implementation.
// TODO: Rename this file to Matrix.h once we remove Matrix.h in sample_util.
#ifndef COMMON_MATRIX_UTILS_H_
#define COMMON_MATRIX_UTILS_H_
#include <vector>
#include "common/debug.h"
#include "common/mathutil.h"
namespace angle
{
template<typename T>
class Matrix
{
public:
Matrix(const std::vector<T> &elements, const unsigned int &numRows, const unsigned int &numCols)
: mElements(elements),
mRows(numRows),
mCols(numCols)
{
ASSERT(rows() >= 1 && rows() <= 4);
ASSERT(columns() >= 1 && columns() <= 4);
}
Matrix(const std::vector<T> &elements, const unsigned int &size)
: mElements(elements),
mRows(size),
mCols(size)
{
ASSERT(rows() >= 1 && rows() <= 4);
ASSERT(columns() >= 1 && columns() <= 4);
}
Matrix(const T *elements, const unsigned int &size)
: mRows(size),
mCols(size)
{
ASSERT(rows() >= 1 && rows() <= 4);
ASSERT(columns() >= 1 && columns() <= 4);
for (size_t i = 0; i < size * size; i++)
mElements.push_back(elements[i]);
}
const T &operator()(const unsigned int &rowIndex, const unsigned int &columnIndex) const
{
return mElements[rowIndex * columns() + columnIndex];
}
T &operator()(const unsigned int &rowIndex, const unsigned int &columnIndex)
{
return mElements[rowIndex * columns() + columnIndex];
}
const T &at(const unsigned int &rowIndex, const unsigned int &columnIndex) const
{
return operator()(rowIndex, columnIndex);
}
Matrix<T> operator*(const Matrix<T> &m)
{
ASSERT(columns() == m.rows());
unsigned int resultRows = rows();
unsigned int resultCols = m.columns();
Matrix<T> result(std::vector<T>(resultRows * resultCols), resultRows, resultCols);
for (unsigned int i = 0; i < resultRows; i++)
{
for (unsigned int j = 0; j < resultCols; j++)
{
T tmp = 0.0f;
for (unsigned int k = 0; k < columns(); k++)
tmp += at(i, k) * m(k, j);
result(i, j) = tmp;
}
}
return result;
}
unsigned int size() const
{
ASSERT(rows() == columns());
return rows();
}
unsigned int rows() const { return mRows; }
unsigned int columns() const { return mCols; }
std::vector<T> elements() const { return mElements; }
Matrix<T> compMult(const Matrix<T> &mat1) const
{
Matrix result(std::vector<T>(mElements.size()), size());
for (unsigned int i = 0; i < columns(); i++)
for (unsigned int j = 0; j < rows(); j++)
result(i, j) = at(i, j) * mat1(i, j);
return result;
}
Matrix<T> outerProduct(const Matrix<T> &mat1) const
{
unsigned int cols = mat1.columns();
Matrix result(std::vector<T>(rows() * cols), rows(), cols);
for (unsigned int i = 0; i < rows(); i++)
for (unsigned int j = 0; j < cols; j++)
result(i, j) = at(i, 0) * mat1(0, j);
return result;
}
Matrix<T> transpose() const
{
Matrix result(std::vector<T>(mElements.size()), columns(), rows());
for (unsigned int i = 0; i < columns(); i++)
for (unsigned int j = 0; j < rows(); j++)
result(i, j) = at(j, i);
return result;
}
T determinant() const
{
ASSERT(rows() == columns());
switch (size())
{
case 2:
return at(0, 0) * at(1, 1) - at(0, 1) * at(1, 0);
case 3:
return at(0, 0) * at(1, 1) * at(2, 2) +
at(0, 1) * at(1, 2) * at(2, 0) +
at(0, 2) * at(1, 0) * at(2, 1) -
at(0, 2) * at(1, 1) * at(2, 0) -
at(0, 1) * at(1, 0) * at(2, 2) -
at(0, 0) * at(1, 2) * at(2, 1);
case 4:
{
const float minorMatrices[4][3 * 3] =
{
{
at(1, 1), at(2, 1), at(3, 1),
at(1, 2), at(2, 2), at(3, 2),
at(1, 3), at(2, 3), at(3, 3),
},
{
at(1, 0), at(2, 0), at(3, 0),
at(1, 2), at(2, 2), at(3, 2),
at(1, 3), at(2, 3), at(3, 3),
},
{
at(1, 0), at(2, 0), at(3, 0),
at(1, 1), at(2, 1), at(3, 1),
at(1, 3), at(2, 3), at(3, 3),
},
{
at(1, 0), at(2, 0), at(3, 0),
at(1, 1), at(2, 1), at(3, 1),
at(1, 2), at(2, 2), at(3, 2),
}
};
return at(0, 0) * Matrix<T>(minorMatrices[0], 3).determinant() -
at(0, 1) * Matrix<T>(minorMatrices[1], 3).determinant() +
at(0, 2) * Matrix<T>(minorMatrices[2], 3).determinant() -
at(0, 3) * Matrix<T>(minorMatrices[3], 3).determinant();
}
default:
UNREACHABLE();
break;
}
return T();
}
Matrix<T> inverse() const
{
ASSERT(rows() == columns());
Matrix<T> cof(std::vector<T>(mElements.size()), rows(), columns());
switch (size())
{
case 2:
cof(0, 0) = at(1, 1);
cof(0, 1) = -at(1, 0);
cof(1, 0) = -at(0, 1);
cof(1, 1) = at(0, 0);
break;
case 3:
cof(0, 0) = at(1, 1) * at(2, 2) -
at(2, 1) * at(1, 2);
cof(0, 1) = -(at(1, 0) * at(2, 2) -
at(2, 0) * at(1, 2));
cof(0, 2) = at(1, 0) * at(2, 1) -
at(2, 0) * at(1, 1);
cof(1, 0) = -(at(0, 1) * at(2, 2) -
at(2, 1) * at(0, 2));
cof(1, 1) = at(0, 0) * at(2, 2) -
at(2, 0) * at(0, 2);
cof(1, 2) = -(at(0, 0) * at(2, 1) -
at(2, 0) * at(0, 1));
cof(2, 0) = at(0, 1) * at(1, 2) -
at(1, 1) * at(0, 2);
cof(2, 1) = -(at(0, 0) * at(1, 2) -
at(1, 0) * at(0, 2));
cof(2, 2) = at(0, 0) * at(1, 1) -
at(1, 0) * at(0, 1);
break;
case 4:
cof(0, 0) = at(1, 1) * at(2, 2) * at(3, 3) +
at(2, 1) * at(3, 2) * at(1, 3) +
at(3, 1) * at(1, 2) * at(2, 3) -
at(1, 1) * at(3, 2) * at(2, 3) -
at(2, 1) * at(1, 2) * at(3, 3) -
at(3, 1) * at(2, 2) * at(1, 3);
cof(0, 1) = -(at(1, 0) * at(2, 2) * at(3, 3) +
at(2, 0) * at(3, 2) * at(1, 3) +
at(3, 0) * at(1, 2) * at(2, 3) -
at(1, 0) * at(3, 2) * at(2, 3) -
at(2, 0) * at(1, 2) * at(3, 3) -
at(3, 0) * at(2, 2) * at(1, 3));
cof(0, 2) = at(1, 0) * at(2, 1) * at(3, 3) +
at(2, 0) * at(3, 1) * at(1, 3) +
at(3, 0) * at(1, 1) * at(2, 3) -
at(1, 0) * at(3, 1) * at(2, 3) -
at(2, 0) * at(1, 1) * at(3, 3) -
at(3, 0) * at(2, 1) * at(1, 3);
cof(0, 3) = -(at(1, 0) * at(2, 1) * at(3, 2) +
at(2, 0) * at(3, 1) * at(1, 2) +
at(3, 0) * at(1, 1) * at(2, 2) -
at(1, 0) * at(3, 1) * at(2, 2) -
at(2, 0) * at(1, 1) * at(3, 2) -
at(3, 0) * at(2, 1) * at(1, 2));
cof(1, 0) = -(at(0, 1) * at(2, 2) * at(3, 3) +
at(2, 1) * at(3, 2) * at(0, 3) +
at(3, 1) * at(0, 2) * at(2, 3) -
at(0, 1) * at(3, 2) * at(2, 3) -
at(2, 1) * at(0, 2) * at(3, 3) -
at(3, 1) * at(2, 2) * at(0, 3));
cof(1, 1) = at(0, 0) * at(2, 2) * at(3, 3) +
at(2, 0) * at(3, 2) * at(0, 3) +
at(3, 0) * at(0, 2) * at(2, 3) -
at(0, 0) * at(3, 2) * at(2, 3) -
at(2, 0) * at(0, 2) * at(3, 3) -
at(3, 0) * at(2, 2) * at(0, 3);
cof(1, 2) = -(at(0, 0) * at(2, 1) * at(3, 3) +
at(2, 0) * at(3, 1) * at(0, 3) +
at(3, 0) * at(0, 1) * at(2, 3) -
at(0, 0) * at(3, 1) * at(2, 3) -
at(2, 0) * at(0, 1) * at(3, 3) -
at(3, 0) * at(2, 1) * at(0, 3));
cof(1, 3) = at(0, 0) * at(2, 1) * at(3, 2) +
at(2, 0) * at(3, 1) * at(0, 2) +
at(3, 0) * at(0, 1) * at(2, 2) -
at(0, 0) * at(3, 1) * at(2, 2) -
at(2, 0) * at(0, 1) * at(3, 2) -
at(3, 0) * at(2, 1) * at(0, 2);
cof(2, 0) = at(0, 1) * at(1, 2) * at(3, 3) +
at(1, 1) * at(3, 2) * at(0, 3) +
at(3, 1) * at(0, 2) * at(1, 3) -
at(0, 1) * at(3, 2) * at(1, 3) -
at(1, 1) * at(0, 2) * at(3, 3) -
at(3, 1) * at(1, 2) * at(0, 3);
cof(2, 1) = -(at(0, 0) * at(1, 2) * at(3, 3) +
at(1, 0) * at(3, 2) * at(0, 3) +
at(3, 0) * at(0, 2) * at(1, 3) -
at(0, 0) * at(3, 2) * at(1, 3) -
at(1, 0) * at(0, 2) * at(3, 3) -
at(3, 0) * at(1, 2) * at(0, 3));
cof(2, 2) = at(0, 0) * at(1, 1) * at(3, 3) +
at(1, 0) * at(3, 1) * at(0, 3) +
at(3, 0) * at(0, 1) * at(1, 3) -
at(0, 0) * at(3, 1) * at(1, 3) -
at(1, 0) * at(0, 1) * at(3, 3) -
at(3, 0) * at(1, 1) * at(0, 3);
cof(2, 3) = -(at(0, 0) * at(1, 1) * at(3, 2) +
at(1, 0) * at(3, 1) * at(0, 2) +
at(3, 0) * at(0, 1) * at(1, 2) -
at(0, 0) * at(3, 1) * at(1, 2) -
at(1, 0) * at(0, 1) * at(3, 2) -
at(3, 0) * at(1, 1) * at(0, 2));
cof(3, 0) = -(at(0, 1) * at(1, 2) * at(2, 3) +
at(1, 1) * at(2, 2) * at(0, 3) +
at(2, 1) * at(0, 2) * at(1, 3) -
at(0, 1) * at(2, 2) * at(1, 3) -
at(1, 1) * at(0, 2) * at(2, 3) -
at(2, 1) * at(1, 2) * at(0, 3));
cof(3, 1) = at(0, 0) * at(1, 2) * at(2, 3) +
at(1, 0) * at(2, 2) * at(0, 3) +
at(2, 0) * at(0, 2) * at(1, 3) -
at(0, 0) * at(2, 2) * at(1, 3) -
at(1, 0) * at(0, 2) * at(2, 3) -
at(2, 0) * at(1, 2) * at(0, 3);
cof(3, 2) = -(at(0, 0) * at(1, 1) * at(2, 3) +
at(1, 0) * at(2, 1) * at(0, 3) +
at(2, 0) * at(0, 1) * at(1, 3) -
at(0, 0) * at(2, 1) * at(1, 3) -
at(1, 0) * at(0, 1) * at(2, 3) -
at(2, 0) * at(1, 1) * at(0, 3));
cof(3, 3) = at(0, 0) * at(1, 1) * at(2, 2) +
at(1, 0) * at(2, 1) * at(0, 2) +
at(2, 0) * at(0, 1) * at(1, 2) -
at(0, 0) * at(2, 1) * at(1, 2) -
at(1, 0) * at(0, 1) * at(2, 2) -
at(2, 0) * at(1, 1) * at(0, 2);
break;
default:
UNREACHABLE();
break;
}
// The inverse of A is the transpose of the cofactor matrix times the reciprocal of the determinant of A.
Matrix<T> adjugateMatrix(cof.transpose());
T det = determinant();
Matrix<T> result(std::vector<T>(mElements.size()), rows(), columns());
for (unsigned int i = 0; i < rows(); i++)
for (unsigned int j = 0; j < columns(); j++)
result(i, j) = det ? adjugateMatrix(i, j) / det : T();
return result;
}
void setToIdentity()
{
ASSERT(rows() == columns());
const auto one = T(1);
const auto zero = T(0);
for (auto &e : mElements)
e = zero;
for (unsigned int i = 0; i < rows(); ++i)
{
const auto pos = i * columns() + (i % columns());
mElements[pos] = one;
}
}
template <unsigned int Size>
static void setToIdentity(T(&matrix)[Size])
{
static_assert(gl::iSquareRoot<Size>() != 0, "Matrix is not square.");
const auto cols = gl::iSquareRoot<Size>();
const auto one = T(1);
const auto zero = T(0);
for (auto &e : matrix)
e = zero;
for (unsigned int i = 0; i < cols; ++i)
{
const auto pos = i * cols + (i % cols);
matrix[pos] = one;
}
}
private:
std::vector<T> mElements;
unsigned int mRows;
unsigned int mCols;
};
} // namespace angle
#endif // COMMON_MATRIX_UTILS_H_