checkpoint
[anneal.git] / src / edu / berkeley / qfat / geom / Matrix.java
index 7d14a71..514f7d7 100644 (file)
@@ -1,41 +1,44 @@
 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 ]
-    //  [ m n o p ]   [ 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 Matrix(float scalex, float scaley, float scalez) {
-        a = scalex;
-        f = scaley;
-        k = scalez;
-        l = h = d = e = b = i = c = j = g = 0;            
-        m = n = o = 0;
-        p = 1;
-    }
-    public Matrix(Vec translate) {
-        d = translate.x; h = translate.y; l = translate.z;
-        a = f = k = 1;
-        b = c = e = g = i = j = 0;
-        m = n = o = 0;
-        p = 1;
-    }
+
+    /** the zero matrix */
+    public static final Matrix ZERO = new Matrix(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
+
+    /** the identity matrix */
+    public static final Matrix ONE  = new Matrix(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1);
+
+    /** a scaling matrix (uniform in all dimensions) */
+    public static Matrix scale(float scale) { return scale(scale, scale, scale); }
+
+    /** a scaling matrix */
+    public static Matrix scale(float scalex, float scaley, float scalez) {
+        return new Matrix(scalex, 0, 0, 0, 0, scaley, 0, 0, 0, 0, scalez, 0, 0, 0, 0, 1); }
+
     public Matrix(float a, float b, float c, float d, float e, float f, float g,
                   float h, float i, float j, float k, float l, float m, float n, float o, float p) {
         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 static Matrix translate(Vec translate) {
+        return new Matrix(1, 0, 0, translate.x,
+                          0, 1, 0, translate.y,
+                          0, 0, 1, translate.z,
+                          0, 0, 0, 1);
+    }
+
     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);
     }
@@ -56,28 +59,29 @@ public class Matrix {
         return ret;
     }
 
-    public Matrix(Vec axis, float angle) {
-        double q = Math.cos(angle);
-        double s = Math.sin(angle);
-        double t = 1.0 - q;
-        a = (float)(q + axis.x*axis.x*t);
-        f = (float)(q + axis.y*axis.y*t);
-        k = (float)(q + axis.z*axis.z*t);
-        double tmp1 = axis.x*axis.y*t;
-        double tmp2 = axis.z*s;
-        e = (float)(tmp1 + tmp2);
-        b = (float)(tmp1 - tmp2);
+    public static Matrix rotate(Vec axis, float angle) {
+        float q = (float)Math.cos(angle);
+        float s = (float)Math.sin(angle);
+        float t = 1 - q;
+        float a = (float)(q + axis.x*axis.x*t);
+        float f = (float)(q + axis.y*axis.y*t);
+        float k = (float)(q + axis.z*axis.z*t);
+        float tmp1 = axis.x*axis.y*t;
+        float tmp2 = axis.z*s;
+        float e = (float)(tmp1 + tmp2);
+        float b = (float)(tmp1 - tmp2);
         tmp1 = axis.x*axis.z*t;
         tmp2 = axis.y*s;
-        i = (float)(tmp1 - tmp2);
-        c = (float)(tmp1 + tmp2);
+        float i = (float)(tmp1 - tmp2);
+        float c = (float)(tmp1 + tmp2);
         tmp1 = axis.y*axis.z*t;
         tmp2 = axis.x*s;
-        j = (float)(tmp1 + tmp2);
-        g = (float)(tmp1 - tmp2);
-        d = h = l = 0;
-        m = n = o = 0;
-        p = 1;
+        float j = (float)(tmp1 + tmp2);
+        float g = (float)(tmp1 - tmp2);
+        return new Matrix(a, b, c, 0,
+                          e, f, g, 0,
+                          i, j, k, 0,
+                          0, 0, 0, 1);
     }
     public Point times(Point p) {
         // discards bottom row
@@ -95,4 +99,68 @@ public class Matrix {
     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());
+    }
 }