/* * 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; }