aboutsummaryrefslogtreecommitdiffstats
path: root/matrix.c
diff options
context:
space:
mode:
Diffstat (limited to 'matrix.c')
-rw-r--r--matrix.c57
1 files changed, 57 insertions, 0 deletions
diff --git a/matrix.c b/matrix.c
new file mode 100644
index 0000000..3d86681
--- /dev/null
+++ b/matrix.c
@@ -0,0 +1,57 @@
+/*
+ * Matrix inversion for a 3x3 matrix
+ */
+
+#include "color.h"
+
+/*
+ * X = A^-1
+ *
+ * This uses the "method of adjugate" to calculate the inverse
+ * matrix. If the matrix is singular (meaning the columns are
+ * linearly equivalent "in a straight line") then the result will
+ * contain infinities.
+ */
+void color_mat_inverse(double X[3][3], const double A[3][3])
+{
+ double a[3][3]; /* Adjugate matrix */
+ double det; /* Determinant */
+ double invdet; /* Inverse determinant */
+ int i, j;
+
+ /*
+ * The "matrix of cofactors" is the determinate of each
+ * individual submatrix ignoring the row and column for
+ * which the value is computed, multiplied with (-1)^(i+j)
+ *
+ * The "adjugate" is the transpose of this matrix.
+ */
+
+ a[0][0] = A[1][1]*A[2][2] - A[1][2]*A[2][1]; /* \ - / */
+ a[0][1] = A[0][2]*A[2][1] - A[0][1]*A[2][2]; /* / - \ */
+ a[0][2] = A[0][1]*A[1][2] - A[0][2]*A[1][1]; /* \ - / */
+ a[1][0] = A[1][2]*A[2][0] - A[1][0]*A[2][2]; /* / - \ */
+ a[1][1] = A[0][0]*A[2][2] - A[0][2]*A[2][0]; /* \ - / */
+ a[1][2] = A[0][2]*A[1][0] - A[0][0]*A[1][2]; /* / - \ */
+ a[2][0] = A[1][0]*A[2][1] - A[1][1]*A[2][0]; /* \ - / */
+ a[2][1] = A[0][1]*A[2][0] - A[0][0]*A[2][1]; /* / - \ */
+ a[2][2] = A[0][0]*A[1][1] - A[0][1]*A[1][0]; /* \ - / */
+
+ /*
+ * The determinant is now simply the inner product of the
+ * first column of the adjugate with top row of the original
+ * matrix (or any other row/column)
+ */
+ det = A[0][0]*a[0][0] + A[1][0]*a[0][1] + A[2][0]*a[0][2];
+ invdet = 1.0/det;
+
+ /*
+ * The inverse is the adjugate divided with the determinant.
+ * If the matrix is singular (non-invertible), then
+ * det = 0 and invdet = infinity, so we get a matrix of
+ * infinities here.
+ */
+ for (i = 0; i < 3; i++)
+ for (j = 0; j < 3; j++)
+ X[i][j] = a[i][j] * invdet;
+}