/*
 * mx3.c: routines for manipulating 3x3 matrices
 * These matrices are commonly used to represent a 2-D projective mapping.
 */

#include <math.h>
#include "mx3.h"

/* mx3_adjoint: b = adjoint(a), returns determinant(a) */

double mx3d_adjoint(double a[3][3], double b[3][3])
{
    b[0][0] = DET2(a[1][1], a[1][2], a[2][1], a[2][2]);
    b[1][0] = DET2(a[1][2], a[1][0], a[2][2], a[2][0]);
    b[2][0] = DET2(a[1][0], a[1][1], a[2][0], a[2][1]);
    b[0][1] = DET2(a[2][1], a[2][2], a[0][1], a[0][2]);
    b[1][1] = DET2(a[2][2], a[2][0], a[0][2], a[0][0]);
    b[2][1] = DET2(a[2][0], a[2][1], a[0][0], a[0][1]);
    b[0][2] = DET2(a[0][1], a[0][2], a[1][1], a[1][2]);
    b[1][2] = DET2(a[0][2], a[0][0], a[1][2], a[1][0]);
    b[2][2] = DET2(a[0][0], a[0][1], a[1][0], a[1][1]);
    return a[0][0]*b[0][0] + a[0][1]*b[0][1] + a[0][2]*b[0][2];
}

/* mx3_mul: matrix multiply: c = a*b */

void mx3d_mul(double a[3][3], double b[3][3], double c[3][3])
{
    c[0][0] = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0];
    c[0][1] = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1];
    c[0][2] = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2];
    c[1][0] = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0];
    c[1][1] = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1];
    c[1][2] = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2];
    c[2][0] = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0];
    c[2][1] = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1];
    c[2][2] = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2];
}

/* mx3d_transform: transform point p by matrix a: q = p*a */

void mx3d_transform(double p[3], double a[3][3], double q[3])
{
    q[0] = p[0]*a[0][0] + p[1]*a[1][0] + p[2]*a[2][0];
    q[1] = p[0]*a[0][1] + p[1]*a[1][1] + p[2]*a[2][1];
    q[2] = p[0]*a[0][2] + p[1]*a[1][2] + p[2]*a[2][2];
}