checkpoint
[anneal.git] / src / edu / berkeley / qfat / geom / Matrix.java
index f617376..6712b28 100644 (file)
@@ -1,21 +1,23 @@
 package edu.berkeley.qfat.geom;
 
-/** affine matrix; immutable */
+/** an affine matrix; immutable */
 public class Matrix {
-    //
-    //  [ a b c d ]   [ x ]
-    //  [ e f g h ]   [ y ]
-    //  [ i j k l ]   [ z ]
-    //  [ 0 0 0 1 ]   [ 1 ]
-    //
+
+    /**
+     *  <pre>
+     *  [ a b c d ]   [ x ]
+     *  [ e f g h ]   [ y ]
+     *  [ i j k l ]   [ z ]
+     *  [ m n o p ]   [ 1 ]
+     *  </pre>
+     */
     public final float a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p;
-    public Matrix() { this(1); }
-    public Matrix(float scale) {
-        a = f = k = scale;
-        l = h = d = e = b = i = c = j = g = 0;            
-        m = n = o = 0;
-        p = 1;
-    }
+
+    public static final Matrix ZERO = new Matrix(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
+    public static final Matrix ONE  = new Matrix(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1);
+
+    public static Matrix scale(float scale) { return new Matrix(scale, 0, 0, 0, 0, scale, 0, 0, 0, 0, scale, 0, 0, 0, 0, 1); }
+
     public Matrix(float scalex, float scaley, float scalez) {
         a = scalex;
         f = scaley;
@@ -36,9 +38,26 @@ public class Matrix {
         this.a = a; this.b = b; this.c = c; this.d = d; this.e = e; this.f = f; this.g = g; this.h = h; this.i = i;
         this.j = j; this.k = k; this.l = l; this.m = m; this.n = n; this.o = o; this.p = p;
     }
+    public Matrix plus(Matrix x) {
+        return new Matrix(a+x.a, b+x.b, c+x.c, d+x.d, e+x.e, f+x.f, g+x.g, h+x.h, i+x.i, j+x.j, k+x.k, l+x.l, m+x.m, n+x.n, o+x.o, p+x.p);
+    }
+    public Matrix minus(Matrix x) {
+        return new Matrix(a-x.a, b-x.b, c-x.c, d-x.d, e-x.e, f-x.f, g-x.g, h-x.h, i-x.i, j-x.j, k-x.k, l-x.l, m-x.m, n-x.n, o-x.o, p-x.p);
+    }
     public Matrix times(float x) {
         return new Matrix(a*x, b*x, c*x, d*x, e*x, f*x, g*x, h*x, i*x, j*x, k*x, l*x, m*x, n*x, o*x, p*x);
     }
+
+    /** computes (v^T)(this)(v) */
+    public float preAndPostMultiply(Point point) {
+        float ret =
+            ((a*point.x + b*point.y + c*point.z + d) * point.x) +
+            ((e*point.x + f*point.y + g*point.z + h) * point.y) +
+            ((i*point.x + j*point.y + k*point.z + l) * point.z) + 
+            ((m*point.x + n*point.y + o*point.z + p) * 1);
+        return ret;
+    }
+
     public Matrix(Vec axis, float angle) {
         double q = Math.cos(angle);
         double s = Math.sin(angle);
@@ -68,8 +87,78 @@ public class Matrix {
                          e*p.x + f*p.y + g*p.z + h,
                          i*p.x + j*p.y + k*p.z + l);
     }
+    public Vec times(Vec p) {
+        // discards bottom row
+        return new Vec(a*p.x + b*p.y + c*p.z + d,
+                       e*p.x + f*p.y + g*p.z + h,
+                       i*p.x + j*p.y + k*p.z + l);
+    }
     public Point apply(Point p) { return p; }
     public Vec apply(Vec v) { return v; }
     public Matrix invert() { return this; }
     public Matrix times(Matrix m) { return this; }
+
+
+    public float determinant() {
+        float m00 = a;
+        float m01 = b;
+        float m02 = c;
+        float m03 = d;
+        float m10 = e;
+        float m11 = f;
+        float m12 = g;
+        float m13 = h;
+        float m20 = i;
+        float m21 = j;
+        float m22 = k;
+        float m23 = l;
+        float m30 = m;
+        float m31 = n;
+        float m32 = o;
+        float m33 = p;
+        return
+            m03 * m12 * m21 * m30-m02 * m13 * m21 * m30-m03 * m11 * m22 * m30+m01 * m13    * m22 * m30+
+            m02 * m11 * m23 * m30-m01 * m12 * m23 * m30-m03 * m12 * m20 * m31+m02 * m13    * m20 * m31+
+            m03 * m10 * m22 * m31-m00 * m13 * m22 * m31-m02 * m10 * m23 * m31+m00 * m12    * m23 * m31+
+            m03 * m11 * m20 * m32-m01 * m13 * m20 * m32-m03 * m10 * m21 * m32+m00 * m13    * m21 * m32+
+            m01 * m10 * m23 * m32-m00 * m11 * m23 * m32-m02 * m11 * m20 * m33+m01 * m12    * m20 * m33+
+            m02 * m10 * m21 * m33-m00 * m12 * m21 * m33-m01 * m10 * m22 * m33+m00 * m11    * m22 * m33;
+    }
+
+    public Matrix inverse() {
+        float m00 = a;
+        float m01 = b;
+        float m02 = c;
+        float m03 = d;
+        float m10 = e;
+        float m11 = f;
+        float m12 = g;
+        float m13 = h;
+        float m20 = i;
+        float m21 = j;
+        float m22 = k;
+        float m23 = l;
+        float m30 = m;
+        float m31 = n;
+        float m32 = o;
+        float m33 = p;
+        return
+            new Matrix(m12*m23*m31 - m13*m22*m31 + m13*m21*m32 - m11*m23*m32 - m12*m21*m33 + m11*m22*m33,
+                       m03*m22*m31 - m02*m23*m31 - m03*m21*m32 + m01*m23*m32 + m02*m21*m33 - m01*m22*m33,
+                       m02*m13*m31 - m03*m12*m31 + m03*m11*m32 - m01*m13*m32 - m02*m11*m33 + m01*m12*m33,
+                       m03*m12*m21 - m02*m13*m21 - m03*m11*m22 + m01*m13*m22 + m02*m11*m23 - m01*m12*m23,
+                       m13*m22*m30 - m12*m23*m30 - m13*m20*m32 + m10*m23*m32 + m12*m20*m33 - m10*m22*m33,
+                       m02*m23*m30 - m03*m22*m30 + m03*m20*m32 - m00*m23*m32 - m02*m20*m33 + m00*m22*m33,
+                       m03*m12*m30 - m02*m13*m30 - m03*m10*m32 + m00*m13*m32 + m02*m10*m33 - m00*m12*m33,
+                       m02*m13*m20 - m03*m12*m20 + m03*m10*m22 - m00*m13*m22 - m02*m10*m23 + m00*m12*m23,
+                       m11*m23*m30 - m13*m21*m30 + m13*m20*m31 - m10*m23*m31 - m11*m20*m33 + m10*m21*m33,
+                       m03*m21*m30 - m01*m23*m30 - m03*m20*m31 + m00*m23*m31 + m01*m20*m33 - m00*m21*m33,
+                       m01*m13*m30 - m03*m11*m30 + m03*m10*m31 - m00*m13*m31 - m01*m10*m33 + m00*m11*m33,
+                       m03*m11*m20 - m01*m13*m20 - m03*m10*m21 + m00*m13*m21 + m01*m10*m23 - m00*m11*m23,
+                       m12*m21*m30 - m11*m22*m30 - m12*m20*m31 + m10*m22*m31 + m11*m20*m32 - m10*m21*m32,
+                       m01*m22*m30 - m02*m21*m30 + m02*m20*m31 - m00*m22*m31 - m01*m20*m32 + m00*m21*m32,
+                       m02*m11*m30 - m01*m12*m30 - m02*m10*m31 + m00*m12*m31 + m01*m10*m32 - m00*m11*m32,
+                       m01*m12*m20 - m02*m11*m20 + m02*m10*m21 - m00*m12*m21 - m01*m10*m22 + m00*m11*m22)
+            .times(1/determinant());
+    }
 }