From 0621ca86e167cfe43041c1261fb41fb0a8e6d48d Mon Sep 17 00:00:00 2001 From: Godzil Date: Fri, 14 Feb 2020 19:19:36 +0000 Subject: [PATCH] Change width to size, as it is more correct. Add calculation of determinant, submatric, minorm cofactor and inverse of a matrix. --- source/include/matrix.h | 37 ++++--- source/matrix.cpp | 104 +++++++++++++++---- tests/matrix_test.cpp | 217 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 325 insertions(+), 33 deletions(-) diff --git a/source/include/matrix.h b/source/include/matrix.h index b83a30e..1589706 100644 --- a/source/include/matrix.h +++ b/source/include/matrix.h @@ -13,24 +13,30 @@ class Matrix { -private: +protected: /* 4x4 is the default */ double data[4*4]; - int width; + int size; public: - Matrix(int width); - Matrix(double values[], int width); + Matrix(int size); + Matrix(double values[], int size); - double get(int x, int y) const { return this->data[this->width * x + y]; }; - void set(int x, int y, double v) { this->data[this->width * x + y] = v; }; + double get(int x, int y) const { return this->data[this->size * x + y]; }; + void set(int x, int y, double v) { this->data[this->size * x + y] = v; }; Matrix identity(); Matrix transpose(); - + double determinant(); + Matrix submatrix(int row, int column); + Matrix inverse(); + double minor(int row, int column) { return this->submatrix(row, column).determinant(); } + double cofactor(int row, int column) { return (((column+row)&1)?-1:1) * this->minor(row, column); } bool operator==(const Matrix &b) const; bool operator!=(const Matrix &b) const; + bool isInvertible() { return this->determinant() != 0; } + Matrix operator*(const Matrix &b) const; Tuple operator*(const Tuple &b) const; }; @@ -42,14 +48,6 @@ public: Matrix4(double values[]) : Matrix(values, 4) { }; }; - -class Matrix2 : public Matrix -{ -public: - Matrix2() : Matrix(2) { }; - Matrix2(double values[]) : Matrix(values, 2) { }; -}; - class Matrix3 : public Matrix { public: @@ -57,4 +55,13 @@ public: Matrix3(double values[]) : Matrix(values, 3) { }; }; +class Matrix2 : public Matrix +{ +private: + using Matrix::data; +public: + Matrix2() : Matrix(2) { }; + Matrix2(double values[]) : Matrix(values, 2) { }; +}; + #endif /* DORAYME_MATRIX_H */ diff --git a/source/matrix.cpp b/source/matrix.cpp index 930191d..446d6ac 100644 --- a/source/matrix.cpp +++ b/source/matrix.cpp @@ -7,6 +7,7 @@ * */ +#include #include #include #include @@ -15,7 +16,7 @@ Matrix::Matrix(int width) { int i; - this->width = width; + this->size = width; for(i = 0; i < width*width; i++) { @@ -27,13 +28,13 @@ Matrix::Matrix(double values[], int width) { int x, y; - this->width = width; + this->size = width; - for(y = 0; y < this->width; y++) + for(y = 0; y < this->size; y++) { - for (x = 0 ; x < this->width ; x++) + for (x = 0 ; x < this->size ; x++) { - this->data[this->width * x + y] = values[this->width * x + y]; + this->data[this->size * x + y] = values[this->size * x + y]; } } }; @@ -41,13 +42,13 @@ Matrix::Matrix(double values[], int width) bool Matrix::operator==(const Matrix &b) const { int i; - if (this->width != b.width) + if (this->size != b.size) { /* If they are not the same size don't even bother */ return false; } - for(i = 0; i < this->width*this->width; i++) + for(i = 0; i < this->size*this->size; i++) { if (!double_equal(this->data[i], b.data[i])) { @@ -61,13 +62,13 @@ bool Matrix::operator==(const Matrix &b) const bool Matrix::operator!=(const Matrix &b) const { int i; - if (this->width != b.width) + if (this->size != b.size) { /* If they are not the same size don't even bother */ return true; } - for(i = 0; i < this->width*this->width; i++) + for(i = 0; i < this->size*this->size; i++) { if (!double_equal(this->data[i], b.data[i])) { @@ -80,16 +81,16 @@ bool Matrix::operator!=(const Matrix &b) const Matrix Matrix::operator*(const Matrix &b) const { int x, y, k; - Matrix ret = Matrix(this->width); + Matrix ret = Matrix(this->size); - if (this->width == b.width) + if (this->size == b.size) { - for (y = 0 ; y < this->width ; y++) + for (y = 0 ; y < this->size ; y++) { - for (x = 0 ; x < this->width ; x++) + for (x = 0 ; x < this->size ; x++) { double v = 0; - for (k = 0 ; k < this->width ; k++) + for (k = 0 ; k < this->size ; k++) { v += this->get(x, k) * b.get(k, y); } @@ -111,7 +112,7 @@ Tuple Matrix::operator*(const Tuple &b) const Matrix Matrix::identity() { int i; - for(i = 0; i < this->width; i++) + for(i = 0; i < this->size; i++) { this->set(i, i, 1); } @@ -121,13 +122,80 @@ Matrix Matrix::identity() Matrix Matrix::transpose() { int x, y; - Matrix ret = Matrix(this->width); - for (y = 0 ; y < this->width ; y++) + Matrix ret = Matrix(this->size); + for (y = 0 ; y < this->size ; y++) { - for (x = 0 ; x < this->width ; x++) + for (x = 0 ; x < this->size ; x++) { ret.set(y, x, this->get(x, y)); } } + return ret; +} + +Matrix Matrix::submatrix(int row, int column) +{ + int i, j; + int x = 0, y = 0; + Matrix ret = Matrix(this->size - 1); + for (i = 0 ; i < this->size ; i++) + { + if (i == row) + { + /* Skip that row */ + continue; + } + y = 0; + for (j = 0 ; j < this->size ; j++) + { + if (j == column) + { + /* skip that column */ + continue; + } + ret.set(x, y, this->get(i, j)); + y++; + } + x++; + } + return ret; +} + +double Matrix::determinant() +{ + double det = 0; + if (this->size == 2) + { + det = this->data[0] * this->data[3] - this->data[1] * this->data[2]; + } + else + { + int col; + for(col = 0; col < this->size; col++) + { + det = det + this->get(0, col) * this->cofactor(0, col); + } + } + return det; +} + +Matrix Matrix::inverse() +{ + Matrix ret = Matrix(this->size); + + if (this->isInvertible()) + { + int row, col; + double c; + for (row = 0; row < this->size; row++) + { + for (col = 0; col < this->size; col++) + { + c = this->cofactor(row, col); + ret.set(col, row, c / this->determinant()); + } + } + } + return ret; } \ No newline at end of file diff --git a/tests/matrix_test.cpp b/tests/matrix_test.cpp index cee15e6..26a6f05 100644 --- a/tests/matrix_test.cpp +++ b/tests/matrix_test.cpp @@ -170,4 +170,221 @@ TEST(MatrixTest, Transposing_this_identity_matrix) Matrix ident = Matrix4().identity(); ASSERT_EQ(ident.transpose(), ident); +} + +TEST(MatrixTest, Calculating_the_determinant_of_a_2x2_matrix) +{ + double valuesA[] = { 1, 5, + -3, 2 }; + Matrix2 A = Matrix2(valuesA); + + ASSERT_EQ(A.determinant(), 17); +} + +TEST(MatrixTest, A_submatrix_of_a_3x3_matrix_is_a_2x2_matrix) +{ + double valuesA[] = { 1, 5, 0, + -3, 2, 7, + 0, 6, -3 }; + double results[] = { -3, 2, + 0, 6 }; + Matrix3 A = Matrix3(valuesA); + + ASSERT_EQ(A.submatrix(0, 2), Matrix2(results)); +} + +TEST(MatrixTest, A_submatrix_of_a_4x4_matrix_is_a_3x3_matrix) +{ + double valuesA[] = { -6, 1, 1, 6, + -8, 5, 8, 6, + -1, 0, 8, 2, + -7, 1, -1, 1 }; + double results[] = { -6, 1, 6, + -8, 8, 6, + -7,-1, 1 }; + Matrix4 A = Matrix4(valuesA); + + ASSERT_EQ(A.submatrix(2, 1), Matrix3(results)); +} + +TEST(MatrixTest, Calculate_a_minor_of_a_3x3_matrix) +{ + double valuesA[] = { 3, 5, 0, + 2, -1, -7, + 6, -1, 5 }; + + Matrix3 A = Matrix3(valuesA); + + Matrix B = A.submatrix(1, 0); + + ASSERT_EQ(B.determinant(), 25); + ASSERT_EQ(A.minor(1, 0), 25); +} + +TEST(MatrixTest, Calculating_a_cofactor_of_a_3x3_matrix) +{ + double valuesA[] = { 3, 5, 0, + 2, -1, -7, + 6, -1, 5 }; + + Matrix3 A = Matrix3(valuesA); + + ASSERT_EQ(A.minor(0, 0), -12); + ASSERT_EQ(A.cofactor(0, 0), -12); + ASSERT_EQ(A.minor(1, 0), 25); + ASSERT_EQ(A.cofactor(1, 0), -25); +} + +TEST(MatrixTest, Calculating_the_determinant_of_a_3x3_matrix) +{ + double valuesA[] = { 1, 2, 6, + -5, 8, -4, + 2, 6, 4 }; + + Matrix A = Matrix3(valuesA); + + ASSERT_EQ(A.cofactor(0, 0), 56); + ASSERT_EQ(A.cofactor(0, 1), 12); + ASSERT_EQ(A.minor(0, 2), -46); + ASSERT_EQ(A.determinant(), -196); +} + +TEST(MatrixTest, Calculating_the_determinant_of_a_4x4_matrix) +{ + double valuesA[] = { -2, -8, 3, 5, + -3, 1, 7, 3, + 1, 2, -9, 6, + -6, 7, 7, -9 }; + + Matrix A = Matrix4(valuesA); + + ASSERT_EQ(A.cofactor(0, 0), 690); + ASSERT_EQ(A.cofactor(0, 1), 447); + ASSERT_EQ(A.cofactor(0, 2), 210); + ASSERT_EQ(A.minor(0, 3), -51); + ASSERT_EQ(A.determinant(), -4071); +} + +TEST(MatrixTest, Testing_an_invertible_matrix_for_invertibility) +{ + double valuesA[] = { 6, 4, 4, 4, + 5, 5, 7, 6, + 4, -9, 3, -7, + 9, 1, 7, -6 }; + Matrix A = Matrix4(valuesA); + + ASSERT_EQ(A.determinant(), -2120); + ASSERT_TRUE(A.isInvertible()); +} + +TEST(MatrixTest, Testing_an_noninvertible_matrix_for_invertibility) +{ + double valuesA[] = { -4, 2, -2, -3, + 9, 6, 2, 6, + 0, -5, 1, -5, + 0, 0, 0, 0 }; + Matrix A = Matrix4(valuesA); + + ASSERT_EQ(A.determinant(), 0); + ASSERT_FALSE(A.isInvertible()); +} + +TEST(MatrixTest, Calculating_the_inverse_of_a_matrix) +{ + double valuesA[] = { -5, 2, 6, -8, + 1, -5, 1, 8, + 7, 7, -6, -7, + 1, -3, 7, 4 }; + + double results[] = { 0.21805, 0.45113, 0.24060, -0.04511, + -0.80827, -1.45677, -0.44361, 0.52068, + -0.07895, -0.22368, -0.05263, 0.19737, + -0.52256, -0.81391, -0.30075, 0.30639 }; + + Matrix A = Matrix4(valuesA); + Matrix B = A.inverse(); + + ASSERT_EQ(A.determinant(), 532); + ASSERT_EQ(A.cofactor(2, 3), -160); + ASSERT_NEAR(B.get(3, 2), -160./532., DBL_EPSILON); + ASSERT_EQ(A.cofactor(3, 2), 105); + ASSERT_NEAR(B.get(2, 3), 105./532., DBL_EPSILON); + + /* Temporary lower the precision */ + set_equal_precision(0.00001); + + ASSERT_EQ(B, Matrix4(results)); + + /* Revert to default */ + set_equal_precision(FLT_EPSILON); +} + +TEST(MatrixTest, Calculating_the_inverse_of_another_matrix) +{ + double valuesA[] = { 8, -5, 9, 2, + 7, 5, 6, 1, + -6, 0, 9, 6, + -3, 0, -9, -4 }; + + double results[] = { -0.15385, -0.15385, -0.28205, -0.53846, + -0.07692, 0.12308, 0.02564, 0.03077, + 0.35897, 0.35897, 0.43590, 0.92308, + -0.69231, -0.69231, -0.76923, -1.92308 }; + + Matrix A = Matrix4(valuesA); + Matrix B = A.inverse(); + + + /* Temporary lower the precision */ + set_equal_precision(0.00001); + + ASSERT_EQ(B, Matrix4(results)); + + /* Revert to default */ + set_equal_precision(FLT_EPSILON); +} + +TEST(MatrixTest, Calculating_the_inverse_of_third_matrix) +{ + double valuesA[] = { 9, 3, 0, 9, + -5, -2, -6, -3, + -4, 9, 6, 4, + -7, 6, 6, 2 }; + + double results[] = { -0.04074, -0.07778, 0.14444, -0.22222, + -0.07778, 0.03333, 0.36667, -0.33333, + -0.02901, -0.14630, -0.10926, 0.12963, + 0.17778, 0.06667, -0.26667, 0.33333 }; + + Matrix A = Matrix4(valuesA); + Matrix B = A.inverse(); + + + /* Temporary lower the precision */ + set_equal_precision(0.00001); + + ASSERT_EQ(B, Matrix4(results)); + + /* Revert to default */ + set_equal_precision(FLT_EPSILON); +} + +TEST(MatrixTest, Multiplying_a_product_by_its_inverse) +{ + double valuesA[] = { 3, -9, 7, 3, + 3, -8, 2, -9, + -4, 4, 4, 1, + -6, 5, -1, 1 }; + + double valuesB[] = { 8, 2, 2, 2, + 3, -1, 7, 0, + 7, 0, 5, 4, + 6, -2, 0, 5 }; + + Matrix A = Matrix4(valuesA); + Matrix B = Matrix4(valuesB); + + Matrix C = A * B; + + ASSERT_EQ(C * B.inverse(), A); } \ No newline at end of file