aboutsummaryrefslogtreecommitdiffstats log msg author committer range
path: root/matrix.c
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
Diffstat (limited to 'matrix.c')
-rw-r--r--matrix.c57
1 files changed, 57 insertions, 0 deletions
 diff --git a/matrix.c b/matrix.cnew file mode 100644index 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;+}