Cynopolis/Vector3D
1
0
Fork 0
Code Issues 3 Pull Requests 1 Actions Packages Projects Releases Wiki Activity

Working on polishing functions

Browse Source
This commit is contained in:
Quinn Henthorne
2024-12-11 17:30:47 -05:00
parent 9cbaeb2c27
commit 7d7d89c676
3 changed files with 328 additions and 175 deletions
Download Patch File Download Diff File Expand all files Collapse all files

397
Matrix.hpp
View File

@@ -15,6 +15,7 @@ public:
*/
Matrix(const std::array<float, rows * columns> &array);
Matrix(const Matrix<rows, columns> &other);
// TODO: Figure out how to do this
/**
* @brief Initialize a matrix directly with any number of arguments
@@ -28,8 +29,8 @@ public:
* @param result A buffer to store the result into
* @note there is no problem if result == this
*/
void Add(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const;
Matrix<rows, columns> &Add(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const;
/**
* @brief Element-wise subtract matrix
@@ -37,8 +38,8 @@ public:
* @param result A buffer to store the result into
* @note there is no problem if result == this
*/
void Sub(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const;
Matrix<rows, columns> &Sub(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const;
/**
* @brief Matrix multiply the two matrices
@@ -46,8 +47,8 @@ public:
* @param result A buffer to store the result into
*/
template <uint8_t other_columns>
void Mult(const Matrix<rows, columns> &other,
Matrix<columns, other_columns> &result) const;
Matrix<rows, columns> &Mult(const Matrix<columns, other_columns> &other,
Matrix<rows, other_columns> &result) const;
/**
* @brief Multiply the matrix by a scalar
@@ -55,37 +56,14 @@ public:
* @param result A buffer to store the result into
* @note there is no problem if result == this
*/
void Mult(float scalar, Matrix<rows, columns> &result) const;
/**
* @brief Invert this matrix
* @param result A buffer to store the result into
* @warning this is super slow! Only call it if you absolutely have to!!!
*/
void Invert(Matrix<rows, columns> &result) const;
/**
* @brief Transpose this matrix
* @param result A buffer to store the result into
*/
void Transpose(Matrix<columns, rows> &result) const;
Matrix<rows, columns> &Mult(float scalar,
Matrix<rows, columns> &result) const;
/**
* @brief Square this matrix
* @param result A buffer to store the result into
*/
void Square(Matrix<rows, columns> &result) const;
/**
* @return Get the determinant of the matrix
*/
float Det() const;
/**
* @brief Calculate the eigenvalues for a square matrix
* @param result a buffer to store the result into
*/
void EigenValues(Matrix<rows, 1> &eigenvalues) const;
Matrix<rows, columns> &Square(Matrix<rows, rows> &result) const;
/**
* @brief Element-wise multiply the two matrices
@@ -93,8 +71,8 @@ public:
* @param result A buffer to store the result into
* @note there is no problem if result == this
*/
void ElementMultiply(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const;
Matrix<rows, columns> &ElementMultiply(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const;
/**
* @brief Element-wise divide the two matrices
@@ -102,8 +80,59 @@ public:
* @param result A buffer to store the result into
* @note there is no problem if result == this
*/
void ElementDivide(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const;
Matrix<rows, columns> &ElementDivide(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const;
/**
* @return Get the determinant of the matrix
* @note for right now only 2x2 and 3x3 matrices are supported
*/
float Det() const;
/**
* @brief Invert this matrix
* @param result A buffer to store the result into
* @warning this is super slow! Only call it if you absolutely have to!!!
*/
Matrix<rows, columns> &Invert(Matrix<rows, columns> &result) const;
/**
* @brief Transpose this matrix
* @param result A buffer to store the result into
*/
Matrix<columns, rows> &Transpose(Matrix<columns, rows> &result) const;
/**
* @brief reduce the matrix so the sum of its elements equal 1
* @param result a buffer to store the result into
*/
Matrix<rows, columns> &Normalize(Matrix<rows, columns> &result) const;
/**
* @brief Get a row from the matrix
* @param row_index the row index to get
* @param row a buffer to write the row into
*/
Matrix<1, columns> &GetRow(uint8_t row_index, Matrix<1, columns> &row) const;
/**
* @brief Get a row from the matrix
* @param column_index the row index to get
* @param column a buffer to write the row into
*/
Matrix<rows, 1> &GetColumn(uint8_t column_index,
Matrix<rows, 1> &column) const;
/**
* @brief Get the number of rows in this matrix
*/
constexpr uint8_t GetRowSize() { return rows; }
/**
* @brief Get the number of columns in this matrix
*/
constexpr uint8_t GetColumnSize() { return columns; }
void ToString(std::string &stringBuffer) const;
/**
* @brief Get an element from the matrix
@@ -118,6 +147,10 @@ public:
* internal array
*/
std::array<float, columns> &operator[](uint8_t row_index) {
if (row_index > rows - 1) {
return this->matrix[0]; // TODO: We should throw something here instead of
// failing quietly.
}
return this->matrix[row_index];
}
@@ -127,34 +160,10 @@ public:
this->matrix[row_idx][column_idx] = other.Get(row_idx, column_idx);
}
}
// return a reference to ourselves so you can chain together these functions
return *this;
}
/**
* @brief Get a row from the matrix
* @param row_index the row index to get
* @param row a buffer to write the row into
*/
void GetRow(uint8_t row_index, Matrix<1, columns> &row) const;
/**
* @brief Get a row from the matrix
* @param column_index the row index to get
* @param column a buffer to write the row into
*/
void GetColumn(uint8_t column_index, Matrix<rows, 1> &column) const;
/**
* @brief Get the number of rows in this matrix
*/
constexpr uint8_t GetRowSize() { return rows; }
/**
* @brief Get the number of columns in this matrix
*/
constexpr uint8_t GetColumnSize() { return columns; }
void ToString(std::string &stringBuffer) const;
private:
/**
* @brief take the dot product of the two vectors
@@ -163,25 +172,21 @@ private:
static float dotProduct(const Matrix<1, vector_size> &vec1,
const Matrix<1, vector_size> &vec2);
template <uint8_t vector_size>
static float dotProduct(const Matrix<vector_size, 1> &vec1,
const Matrix<vector_size, 1> &vec2);
/**
* @brief Set all elements in this matrix to zero
*/
void zeroMatrix();
void matrixOfMinors(Matrix<rows, columns> &result) const;
Matrix<rows, columns> &matrixOfMinors(Matrix<rows, columns> &result) const;
void minorMatrix(Matrix<rows - 1, columns - 1> &result, uint8_t row_idx,
uint8_t column_idx) const;
Matrix<rows - 1, columns - 1> &
minorMatrix(Matrix<rows - 1, columns - 1> &result, uint8_t row_idx,
uint8_t column_idx) const;
void adjugate(Matrix<rows, columns> &result) const;
/**
* @brief reduce the matrix so the sum of its elements equal 1
* @param result a buffer to store the result into
*/
void normalize(Matrix<rows, columns> &result) const;
constexpr bool isSquare() { return rows == columns; }
Matrix<rows, columns> &adjugate(Matrix<rows, columns> &result) const;
void setMatrixToArray(const std::array<float, rows * columns> &array);
@@ -232,31 +237,46 @@ Matrix<rows, columns>::Matrix(const std::array<float, rows * columns> &array) {
// }
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::Add(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const {
Matrix<rows, columns>::Matrix(const Matrix<rows, columns> &other) {
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
this->matrix[row_idx][column_idx] = other.Get(row_idx, column_idx);
}
}
}
template <uint8_t rows, uint8_t columns>
Matrix<rows, columns> &
Matrix<rows, columns>::Add(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const {
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
result[row_idx][column_idx] =
this->Get(row_idx, column_idx) + other.Get(row_idx, column_idx);
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::Sub(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const {
Matrix<rows, columns> &
Matrix<rows, columns>::Sub(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const {
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
result[row_idx][column_idx] =
this->Get(row_idx, column_idx) - other.Get(row_idx, column_idx);
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
template <uint8_t other_columns>
void Matrix<rows, columns>::Mult(const Matrix<rows, columns> &other,
Matrix<columns, other_columns> &result) const {
Matrix<rows, columns> &
Matrix<rows, columns>::Mult(const Matrix<columns, other_columns> &other,
Matrix<rows, other_columns> &result) const {
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
// get our row
@@ -274,20 +294,25 @@ void Matrix<rows, columns>::Mult(const Matrix<rows, columns> &other,
Matrix<rows, columns>::dotProduct(this_row, other_column_t);
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::Mult(float scalar,
Matrix<rows, columns> &result) const {
Matrix<rows, columns> &
Matrix<rows, columns>::Mult(float scalar, Matrix<rows, columns> &result) const {
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
result[row_idx][column_idx] = this->Get(row_idx, column_idx) * scalar;
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
Matrix<rows, columns> &
Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
// since all matrix sizes have to be statically specified at compile time we
// can do this
static_assert(rows == columns,
@@ -295,10 +320,11 @@ void Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
// unfortunately we can't calculate this at compile time so we'll just reurn
// zeros
float determinant{this->Det()};
if (this->Det() < 0) {
// you can't invert a matrix with a negative determinant
result.zeroMatrix();
return;
return result;
}
// TODO: This algorithm is really inneficient because of the matrix of minors.
@@ -311,138 +337,139 @@ void Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
// now adjugate the matrix and save it in our output
minors.adjugate(result);
float determinant = this->Det();
// scale the result by 1/determinant and we have our answer
result.Mult(1 / determinant);
return result;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::Transpose(Matrix<columns, rows> &result) const {
Matrix<columns, rows> &
Matrix<rows, columns>::Transpose(Matrix<columns, rows> &result) const {
for (uint8_t column_idx{0}; column_idx < rows; column_idx++) {
for (uint8_t row_idx{0}; row_idx < columns; row_idx++) {
result[row_idx][column_idx] = this->Get(column_idx, row_idx);
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::Square(Matrix<rows, columns> &result) const {
static_assert(this->isSquare(), "You can't square an non-square matrix.");
Matrix<rows, columns> &
Matrix<rows, columns>::Square(Matrix<rows, rows> &result) const {
// TODO: Because template requirements are checked before static_assert, this
// never throws an error and fails at the Mult call instead.
static_assert(rows == columns, "You can't square a non-square matrix.");
this->Mult(this, result);
this->Mult(*this, result);
return result;
}
// explicitly define the determinant for a 3x3 matrix because it is definitely
// the fastest way to calculte a 2x2 matrix determinant
template <> float Matrix<2, 2>::Det() const {
return this->matrix[0][0] * this->matrix[1][1] -
this->matrix[0][1] * this->matrix[1][1];
}
// explicitly define the determinant for a 3x3 matrix because it will probably
// be faster than the jacobi method for nxn matrices
template <> float Matrix<3, 3>::Det() const {
float a{this->matrix[0][0]};
float b{this->matrix[0][1]};
float c{this->matrix[0][2]};
Matrix<2, 2> minors{};
this->minorMatrix(minors, 0, 0);
float det = a * minors.Det();
this->minorMatrix(minors, 0, 1);
det -= b * minors.Det();
this->minorMatrix(minors, 0, 2);
det += c * minors.Det();
return det;
}
template <uint8_t rows, uint8_t columns>
float Matrix<rows, columns>::Det() const {
static_assert(this->isSquare(),
"You can't take the determinant of a non-square matrix.");
Matrix<1, columns> eigenValues{};
this->EigenValues(eigenValues);
float determinant{1};
for (uint8_t i{0}; i < columns; i++) {
determinant *= eigenValues.Get(0, i);
}
return determinant;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::EigenValues(Matrix<rows, 1> &eigenvalues) const {
static_assert(rows == columns,
"Eigenvalues can only be computed for square matrices.");
// I got this code from:
// https://www.quora.com/What-is-the-C-code-for-finding-eigenvalues
Matrix<rows, 1> v{};
Matrix<rows, 1> Av{};
Matrix<rows, 1> z{};
"You can't take the determinant of a non-square matrix.");
// static_assert(
// false,
// "Right now this operation isn't supported for matrices bigger than
// 3x3");
// Matrix<1, columns> eigenValues{};
// this->EigenValues(eigenValues);
float d = 0.0;
float d_old = 0.0;
constexpr float convergence_value{1e-6};
constexpr uint16_t max_iterations{500};
// float determinant{1};
// for (uint8_t i{0}; i < columns; i++) {
// determinant *= eigenValues.Get(0, i);
// }
// Initialize v as a random vector
for (int i = 0; i < rows; i++) {
v[0][i] = rand() / RAND_MAX;
}
// run this loop until the eigenvalues converge or we give up
for (uint16_t k{0}; k < max_iterations; k++) {
/* Multiply A by v */
for (int i = 0; i < rows; i++) {
Av[0][i] = 0.0;
for (int j = 0; j < rows; j++) {
Av[0][i] += this->Get(0, i * rows + j) * v[0][j];
}
}
// Calculate the eigenvalue and update v
d_old = d;
d = dot_product(v, Av, rows);
for (int i = 0; i < rows; i++) {
z[0][i] = Av[0][i] - d * v[0][i];
}
z.normalize(z);
for (int i = 0; i < rows; i++) {
v[0][i] = z[0][i];
}
/* Check for convergence */
if (std::fabs(d - d_old) < convergence_value) {
eigenvalues[0][k] = d;
k++;
d = 0.0;
for (int i = 0; i < rows; i++) {
v[0][i] = rand() / RAND_MAX;
}
}
}
// return determinant;
return 0;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::ElementMultiply(
const Matrix<rows, columns> &other, Matrix<rows, columns> &result) const {
Matrix<rows, columns> &
Matrix<rows, columns>::ElementMultiply(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const {
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
result[row_idx][column_idx] =
this->Get(row_idx, column_idx) * other.Get(row_idx, column_idx);
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::ElementDivide(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const {
Matrix<rows, columns> &
Matrix<rows, columns>::ElementDivide(const Matrix<rows, columns> &other,
Matrix<rows, columns> &result) const {
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
result[row_idx][column_idx] =
this->Get(row_idx, column_idx) / other.Get(row_idx, column_idx);
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
float Matrix<rows, columns>::Get(uint8_t row_index,
uint8_t column_index) const {
if (row_index > rows - 1 || column_index > columns - 1) {
return 0; // TODO: We should throw something here instead of failing quietly
}
return this->matrix[row_index][column_index];
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::GetRow(uint8_t row_index,
Matrix<1, columns> &row) const {
Matrix<1, columns> &
Matrix<rows, columns>::GetRow(uint8_t row_index,
Matrix<1, columns> &row) const {
row = Matrix<1, columns>(this->matrix[row_index]);
return row;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::GetColumn(uint8_t column_index,
Matrix<rows, 1> &column) const {
Matrix<rows, 1> &
Matrix<rows, columns>::GetColumn(uint8_t column_index,
Matrix<rows, 1> &column) const {
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
column[row_idx][0] = this->Get(row_idx, column_index);
}
return column;
}
template <uint8_t rows, uint8_t columns>
@@ -471,6 +498,18 @@ float Matrix<rows, columns>::dotProduct(const Matrix<1, vector_size> &vec1,
return sum;
}
template <uint8_t rows, uint8_t columns>
template <uint8_t vector_size>
float Matrix<rows, columns>::dotProduct(const Matrix<vector_size, 1> &vec1,
const Matrix<vector_size, 1> &vec2) {
float sum{0};
for (uint8_t i{0}; i < vector_size; i++) {
sum += vec1.Get(i, 0) * vec2.Get(i, 0);
}
return sum;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::zeroMatrix() {
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
@@ -481,8 +520,8 @@ void Matrix<rows, columns>::zeroMatrix() {
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::matrixOfMinors(
Matrix<rows, columns> &result) const {
Matrix<rows, columns> &
Matrix<rows, columns>::matrixOfMinors(Matrix<rows, columns> &result) const {
Matrix<rows - 1, columns - 1> minorMatrix{};
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
@@ -491,12 +530,14 @@ void Matrix<rows, columns>::matrixOfMinors(
result[row_idx][column_idx] = minorMatrix.Det();
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::minorMatrix(Matrix<rows - 1, columns - 1> &result,
uint8_t row_idx,
uint8_t column_idx) const {
Matrix<rows - 1, columns - 1> &
Matrix<rows, columns>::minorMatrix(Matrix<rows - 1, columns - 1> &result,
uint8_t row_idx, uint8_t column_idx) const {
std::array<float, (rows - 1) * (columns - 1)> subArray{};
for (uint8_t row_iter{0}; row_iter < rows; row_iter++) {
@@ -511,10 +552,12 @@ void Matrix<rows, columns>::minorMatrix(Matrix<rows - 1, columns - 1> &result,
}
result = Matrix<rows - 1, columns - 1>{subArray};
return result;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::adjugate(Matrix<rows, columns> &result) const {
Matrix<rows, columns> &
Matrix<rows, columns>::adjugate(Matrix<rows, columns> &result) const {
for (uint8_t row_iter{0}; row_iter < rows; row_iter++) {
for (uint8_t column_iter{0}; column_iter < columns; column_iter++) {
float sign = ((row_iter + 1) % 2) == 0 ? -1 : 1;
@@ -522,26 +565,34 @@ void Matrix<rows, columns>::adjugate(Matrix<rows, columns> &result) const {
result[row_iter][column_iter] = this->Get(row_iter, column_iter) * sign;
}
}
return result;
}
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::normalize(Matrix<rows, columns> &result) const {
Matrix<rows, columns> &
Matrix<rows, columns>::Normalize(Matrix<rows, columns> &result) const {
float sum{0};
for (uint8_t column_idx{0}; column_idx < rows; column_idx++) {
for (uint8_t row_idx{0}; row_idx < columns; row_idx++) {
sum += this->Get(row_idx, column_idx);
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
float val{this->Get(row_idx, column_idx)};
sum += val * val;
}
}
if (sum == 0) {
// this wouldn't do anything anyways
result.zeroMatrix();
return;
return result;
}
for (uint8_t column_idx{0}; column_idx < rows; column_idx++) {
for (uint8_t row_idx{0}; row_idx < columns; row_idx++) {
sum = sqrt(sum);
for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
result[row_idx][column_idx] = this->Get(row_idx, column_idx) / sum;
}
}
return result;
}