[PULL] Surface transformations v1

[Bill Spitzak]

I have no idea if you need this, but here is free code to do the matrix 
inversion. This is code I wrote and is public domain (although I 
certainly referred to a book to generate it, but I can't see that being 
a copyright problem).

The matrix struct uses aRC as the variable name for the entry in row R 
and column C. It is pretty much this:

struct Matrix {
   // Note transposition to OpenGL order!
   float a00, a10, a20, a30,
	a01, a11, a21, a31,
	a02, a12, a22, a32,
	a03, a13, a23, a33;
};

These are C++ methods on the matrix object, but you can change it to C 
by adding a "this" argument and adding it to all the axy references:

   float determinant(void) const {
   return
      a01*a23*a32*a10 -a01*a22*a33*a10 -a23*a31*a02*a10 +a22*a31*a03*a10
     -a00*a23*a32*a11 +a00*a22*a33*a11 +a23*a30*a02*a11 -a22*a30*a03*a11
     -a01*a23*a30*a12 +a00*a23*a31*a12 +a01*a22*a30*a13 -a00*a22*a31*a13
     -a33*a02*a11*a20 +a32*a03*a11*a20 +a01*a33*a12*a20 -a31*a03*a12*a20
     -a01*a32*a13*a20 +a31*a02*a13*a20 +a33*a02*a10*a21 -a32*a03*a10*a21
     -a00*a33*a12*a21 +a30*a03*a12*a21 +a00*a32*a13*a21 -a30*a02*a13*a21;
   }

// Usually code needs to figure out the determinant anyway, so this
// takes it as an argument. o must not be the same as this.
void invert(float det, Matrix& o) const {
   float d = 1.0f/det;
   o.a00 = (-a23*a32*a11 +a22*a33*a11 +a23*a31*a12 -a22*a31*a13 
-a33*a12*a21 +a32*a13*a21) * d;
   o.a01 = ( a01*a23*a32 -a01*a22*a33 -a23*a31*a02 +a22*a31*a03 
+a33*a02*a21 -a32*a03*a21) * d;
   o.a02 = (-a33*a02*a11 +a32*a03*a11 +a01*a33*a12 -a31*a03*a12 
-a01*a32*a13 +a31*a02*a13) * d;
   m.a03 = ( a23*a02*a11 -a22*a03*a11 -a01*a23*a12 +a01*a22*a13 
+a03*a12*a21 -a02*a13*a21) * d;

   o.a10 = ( a23*a32*a10 -a22*a33*a10 -a23*a30*a12 +a22*a30*a13 
+a33*a12*a20 -a32*a13*a20) * d;
   o.a11 = (-a00*a23*a32 +a00*a22*a33 +a23*a30*a02 -a22*a30*a03 
-a33*a02*a20 +a32*a03*a20) * d;
   o.a12 = ( a33*a02*a10 -a32*a03*a10 -a00*a33*a12 +a30*a03*a12 
+a00*a32*a13 -a30*a02*a13) * d;
   o.a13 = (-a23*a02*a10 +a22*a03*a10 +a00*a23*a12 -a00*a22*a13 
-a03*a12*a20 +a02*a13*a20) * d;

   o.a20 = (-a23*a31*a10 +a23*a30*a11 -a33*a11*a20 +a31*a13*a20 
+a33*a10*a21 -a30*a13*a21) * d;
   o.a21 = (-a01*a23*a30 +a00*a23*a31 +a01*a33*a20 -a31*a03*a20 
-a00*a33*a21 +a30*a03*a21) * d;
   o.a22 = (-a01*a33*a10 +a31*a03*a10 +a00*a33*a11 -a30*a03*a11 
+a01*a30*a13 -a00*a31*a13) * d;
   o.a23 = ( a01*a23*a10 -a00*a23*a11 +a03*a11*a20 -a01*a13*a20 
-a03*a10*a21 +a00*a13*a21) * d;

   o.a30 = ( a22*a31*a10 -a22*a30*a11 +a32*a11*a20 -a31*a12*a20 
-a32*a10*a21 +a30*a12*a21) * d;
   o.a31 = ( a01*a22*a30 -a00*a22*a31 -a01*a32*a20 +a31*a02*a20 
+a00*a32*a21 -a30*a02*a21) * d;
   o.a32 = ( a01*a32*a10 -a31*a02*a10 -a00*a32*a11 +a30*a02*a11 
-a01*a30*a12 +a00*a31*a12) * d;
   o.a33 = (-a01*a22*a10 +a00*a22*a11 -a02*a11*a20 +a01*a12*a20 
+a02*a10*a21 -a00*a12*a21) * d;
}

Pekka Paalanen wrote:

> This series also includes X/MIT-like licenced code for inverting a 4x4
> matrix by LU-decomposition. I could not find one in the web, so I
> implemented it from a book (*not* Numerical Recipes in C or C++, since
> they impose licence issues). The reference is in the source. Anyone
> working on their hobby graphics engine might be interested ;-)
> For testing the matrix inversion alone, there is a standalone unit test
> application.