checkpoint

darcs-hash:20071216010803-5007d-2514cc7d3257bc9a45a179cb1791b52d0f1562d7.gz

diff --git a/src/edu/berkeley/qfat/geom/Matrix.java b/src/edu/berkeley/qfat/geom/Matrix.java
index 514f7d7..a00eb4a 100644 (file)
--- a/src/edu/berkeley/qfat/geom/Matrix.java
+++ b/src/edu/berkeley/qfat/geom/Matrix.java
@@ -1,5 +1,7 @@
 package edu.berkeley.qfat.geom;
 
 package edu.berkeley.qfat.geom;
 

+// FEATURE: precompute/cache determinant?
+

 /** an affine matrix; immutable */
 public class Matrix {
 
 /** an affine matrix; immutable */
 public class Matrix {
 

@@ -19,6 +21,12 @@ public class Matrix {
     /** 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);
 
     /** 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);
 

+    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;
+    }
+

     /** a scaling matrix (uniform in all dimensions) */
     public static Matrix scale(float scale) { return scale(scale, scale, scale); }
 
     /** a scaling matrix (uniform in all dimensions) */
     public static Matrix scale(float scale) { return scale(scale, scale, scale); }
 

@@ -26,12 +34,7 @@ public class 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 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;
-    }
-
+    /** a translation matrix */

     public static Matrix translate(Vec translate) {
         return new Matrix(1, 0, 0, translate.x,
                           0, 1, 0, translate.y,
     public static Matrix translate(Vec translate) {
         return new Matrix(1, 0, 0, translate.x,
                           0, 1, 0, translate.y,

@@ -39,26 +42,7 @@ public class Matrix {
                           0, 0, 0, 1);
     }
 
                           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);
-    }
-    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;
-    }
-
+    /** a rotation matrix, <tt>angle</tt> degrees around <tt>axis</tt> */

     public static Matrix rotate(Vec axis, float angle) {
         float q = (float)Math.cos(angle);
         float s = (float)Math.sin(angle);
     public static Matrix rotate(Vec axis, float angle) {
         float q = (float)Math.cos(angle);
         float s = (float)Math.sin(angle);

@@ -83,24 +67,95 @@ public class Matrix {
                           i, j, k, 0,
                           0, 0, 0, 1);
     }
                           i, j, k, 0,
                           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);
+    }
+    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;
+    }
+
+    /** discards bottom row */

     public Point times(Point p) {
         // discards bottom row
         return new Point(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 times(Point p) {
         // discards bottom row
         return new Point(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);
     }

+
+    /** discards bottom row */

     public Vec times(Vec p) {
     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);
     }
         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; }

 
 

+    /** multiply by another matrix */
+    public Matrix times(Matrix z) {
+        float t00 = a;
+        float t01 = b;
+        float t02 = c;
+        float t03 = d;
+        float t10 = e;
+        float t11 = f;
+        float t12 = g;
+        float t13 = h;
+        float t20 = i;
+        float t21 = j;
+        float t22 = k;
+        float t23 = l;
+        float t30 = m;
+        float t31 = n;
+        float t32 = o;
+        float t33 = p;
+        float m00 = z.a;
+        float m01 = z.b;
+        float m02 = z.c;
+        float m03 = z.d;
+        float m10 = z.e;
+        float m11 = z.f;
+        float m12 = z.g;
+        float m13 = z.h;
+        float m20 = z.i;
+        float m21 = z.j;
+        float m22 = z.k;
+        float m23 = z.l;
+        float m30 = z.m;
+        float m31 = z.n;
+        float m32 = z.o;
+        float m33 = z.p;
+        return new Matrix(t00*m00 + t01*m10 + t02*m20 + t03*m30,
+                          t00*m01 + t01*m11 + t02*m21 + t03*m31,
+                          t00*m02 + t01*m12 + t02*m22 + t03*m32,
+                          t00*m03 + t01*m13 + t02*m23 + t03*m33,
+                          t10*m00 + t11*m10 + t12*m20 + t13*m30,
+                          t10*m01 + t11*m11 + t12*m21 + t13*m31,
+                          t10*m02 + t11*m12 + t12*m22 + t13*m32,
+                          t10*m03 + t11*m13 + t12*m23 + t13*m33,
+                          t20*m00 + t21*m10 + t22*m20 + t23*m30,
+                          t20*m01 + t21*m11 + t22*m21 + t23*m31,
+                          t20*m02 + t21*m12 + t22*m22 + t23*m32,
+                          t20*m03 + t21*m13 + t22*m23 + t23*m33,
+                          t30*m00 + t31*m10 + t32*m20 + t33*m30,
+                          t30*m01 + t31*m11 + t32*m21 + t33*m31,
+                          t30*m02 + t31*m12 + t32*m22 + t33*m32,
+                          t30*m03 + t31*m13 + t32*m23 + t33*m33);
+    }

 
 

+    /** compute the determinant of the matrix */

     public float determinant() {
         float m00 = a;
         float m01 = b;
     public float determinant() {
         float m00 = a;
         float m01 = b;

@@ -127,7 +182,10 @@ public class Matrix {
             m02 * m10 * m21 * m33-m00 * m12 * m21 * m33-m01 * m10 * m22 * m33+m00 * m11    * m22 * m33;
     }
 
             m02 * m10 * m21 * m33-m00 * m12 * m21 * m33-m01 * m10 * m22 * m33+m00 * m11    * m22 * m33;
     }
 

+    /** compute the inverse of the matrix, returns null if not invertible */

     public Matrix inverse() {
     public Matrix inverse() {

+        float determinant = determinant();
+        if (determinant==0) return null;

         float m00 = a;
         float m01 = b;
         float m02 = c;
         float m00 = a;
         float m01 = b;
         float m02 = c;

@@ -161,6 +219,7 @@ public class Matrix {
                        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)
                        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());
+            .times(1/determinant);

     }
     }

+

 }
 }

diff --git a/src/edu/berkeley/qfat/geom/Vec.java b/src/edu/berkeley/qfat/geom/Vec.java
index 0072911..9d1d30b 100644 (file)
--- a/src/edu/berkeley/qfat/geom/Vec.java
+++ b/src/edu/berkeley/qfat/geom/Vec.java
@@ -10,7 +10,6 @@ public final class Vec {
     public Vec cross(Vec v) { return new Vec(y*v.z-z*v.y, z*v.x-x*v.z, x*v.y-y*v.x); }
     public Vec plus(Vec v) { return new Vec(x+v.x, y+v.y, z+v.z); }
     public Vec norm() { return mag()==0 ? this : div(mag()); }
     public Vec cross(Vec v) { return new Vec(y*v.z-z*v.y, z*v.x-x*v.z, x*v.y-y*v.x); }
     public Vec plus(Vec v) { return new Vec(x+v.x, y+v.y, z+v.z); }
     public Vec norm() { return mag()==0 ? this : div(mag()); }

-    public Vec times(Matrix m) { return m.apply(this); }

     public float mag() { return (float)Math.sqrt(x*x+y*y+z*z); }
     public float dot(Vec v) { return x*v.x + y*v.y + z*v.z; }
     public Vec times(float mag) { return new Vec(x*mag, y*mag, z*mag); }
     public float mag() { return (float)Math.sqrt(x*x+y*y+z*z); }
     public float dot(Vec v) { return x*v.x + y*v.y + z*v.z; }
     public Vec times(float mag) { return new Vec(x*mag, y*mag, z*mag); }