Megapixels/src/matrix.h

static float XYZD50_to_D65[] = { 0.9555766f,  -0.0230393f, 0.0631636f,
                                -0.0282895f, 1.0099416f,  0.0210077f,
                                 0.0122982f,  -0.0204830f, 1.3299098f };

static float XYZD65_to_sRGB[] = {  3.2406f,  -1.5372f,   -0.4986f,
                                  -0.9689f,   1.8758f,    0.0415f,
                                   0.0557f,  -0.2040f,    1.0570f };

void multiply_matrices(float a[9], float b[9], float out[9]);

void
invert_matrix(const float in[9], float out[9])
{
        float det = in[0] * (in[4] * in[8] - in[5] * in[7]) -
                    in[1] * (in[3] * in[8] - in[5] * in[6]) +
                    in[2] * (in[3] * in[7] - in[4] * in[7]);
        out[0] = (in[4] * in[8] - in[7] * in[5]) / det;
        out[1] = (in[7] * in[2] - in[1] * in[8]) / det;
        out[2] = (in[1] * in[5] - in[4] * in[2]) / det;
        out[3] = (in[6] * in[5] - in[3] * in[8]) / det;
        out[4] = (in[0] * in[8] - in[6] * in[5]) / det;
        out[5] = (in[3] * in[2] - in[0] * in[5]) / det;
        out[6] = (in[3] * in[7] - in[6] * in[4]) / det;
        out[7] = (in[6] * in[1] - in[0] * in[7]) / det;
        out[8] = (in[0] * in[4] - in[3] * in[1]) / det;
}

void
transpose_matrix(const float in[9], float out[9])
{
        for (int i = 0; i < 3; ++i)
                for (int j = 0; j < 3; ++j)
                        out[i + j * 3] = in[j + i * 3];
}