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

@@ -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(int size);
-    Matrix(double values[], int width);
+    Matrix(double values[], int size);
 
-    double get(int x, int y) const { return this->data[this->width * x + y]; };
+    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->width * x + y] = v; };
+    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 */

source/matrix.cpp

@@ -7,6 +7,7 @@
  *
  */
 
+#include <stdio.h>
 #include <matrix.h>
 #include <tuples.h>
 #include <math_helper.h>
@@ -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);
+    Matrix ret = Matrix(this->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++)
         {
             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;
+}

tests/matrix_test.cpp

@@ -171,3 +171,220 @@ TEST(MatrixTest, Transposing_this_identity_matrix)
 
     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);
+}