X-Git-Url: http://git.megacz.com/?p=anneal.git;a=blobdiff_plain;f=src%2Fedu%2Fberkeley%2Fqfat%2Fgeom%2FMatrix.java;h=37137617786853e7b1eb8d66b490bd8ede5d0086;hp=19e1a9af38f172b09395a793cbf886f229617908;hb=64c6939d78acfa06e7bc380cb713e3800bf16be5;hpb=de98420094a537a70140842ff185e6d113ad270e

diff --git a/src/edu/berkeley/qfat/geom/Matrix.java b/src/edu/berkeley/qfat/geom/Matrix.java
index 19e1a9a..3713761 100644
--- a/src/edu/berkeley/qfat/geom/Matrix.java
+++ b/src/edu/berkeley/qfat/geom/Matrix.java
@@ -1,4 +1,5 @@
 package edu.berkeley.qfat.geom;
+import javax.media.opengl.*;
 
 // FEATURE: precompute/cache determinant?
 
@@ -16,10 +17,21 @@ public class Matrix {
     public final float a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p;
 
     /** 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);
+    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);
+    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 NEGATIVE_ONE  = new Matrix(-1,0,0,0,0,-1,0,0,0,0,-1,0,0,0,0,1);
+
+    public Matrix(double a, double b, double c, double d, double e, double f, double g,
+                  double h, double i, double j, double k, double l, double m, double n, double o, double p) {
+        this((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);
+    }
 
     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) {
@@ -68,6 +80,22 @@ public class Matrix {
                           0, 0, 0, 1);
     }
 
+    /** a reflection matrix across the plane passing through the origin with the specified normal */
+    public static Matrix reflect(Vec v) {
+        Vec reflectionPlaneNormal = v.norm();
+        float a = reflectionPlaneNormal.x;
+        float b = reflectionPlaneNormal.y;
+        float c = reflectionPlaneNormal.z;
+        return
+                new Matrix( 1-2*a*a,  -2*a*b,  -2*a*c, 0,
+                            -2*a*b,  1-2*b*b,  -2*b*c, 0,
+                            -2*a*c,   -2*b*c, 1-2*c*c, 0,
+                            0,       0,       0,       1);
+    }
+
+    /** returns the translational component of this matrix */
+    public Vec getTranslationalComponent() { return this.times(Point.ZERO).minus(Point.ZERO); }
+
     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);
     }
@@ -88,19 +116,16 @@ public class Matrix {
         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);
+        double x = a*p.x + b*p.y + c*p.z + d;
+        double y = e*p.x + f*p.y + g*p.z + h;
+        double z = i*p.x + j*p.y + k*p.z + l;
+        double q = m*p.x + n*p.y + o*p.z + this.p;
+        return new Point(x/q, y/q, z/q);
     }
 
-    /** discards bottom row */
     public Vec times(Vec p) {
-        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 this.times(Matrix.translate(p)).getTranslationalComponent();
     }
 
     /** multiply by another matrix */
@@ -235,22 +260,42 @@ public class Matrix {
         if (!(oo instanceof Matrix)) return false;
         Matrix z = (Matrix)oo;
         return
-            a==z.a && 
-            b==z.b && 
-            c==z.c && 
-            d==z.d && 
-            e==z.e && 
-            f==z.f && 
-            g==z.g && 
-            h==z.h && 
-            i==z.i && 
-            j==z.j && 
-            k==z.k && 
-            l==z.l && 
-            m==z.m && 
-            n==z.n && 
-            o==z.o && 
-            p==z.p;
+            (a==z.a) && 
+            (b==z.b) && 
+            (c==z.c) && 
+            (d==z.d) && 
+            (e==z.e) && 
+            (f==z.f) && 
+            (g==z.g) && 
+            (h==z.h) && 
+            (i==z.i) && 
+            (j==z.j) && 
+            (k==z.k) && 
+            (l==z.l) && 
+            (m==z.m) && 
+            (n==z.n) && 
+            (o==z.o) && 
+            (p==z.p);
+    }
+
+    public boolean equalsModuloEpsilon(Matrix z, float epsilon) {
+        return
+            Math.abs(a-z.a)<epsilon && 
+            Math.abs(b-z.b)<epsilon && 
+            Math.abs(c-z.c)<epsilon && 
+            Math.abs(d-z.d)<epsilon && 
+            Math.abs(e-z.e)<epsilon && 
+            Math.abs(f-z.f)<epsilon && 
+            Math.abs(g-z.g)<epsilon && 
+            Math.abs(h-z.h)<epsilon && 
+            Math.abs(i-z.i)<epsilon && 
+            Math.abs(j-z.j)<epsilon && 
+            Math.abs(k-z.k)<epsilon && 
+            Math.abs(l-z.l)<epsilon && 
+            Math.abs(m-z.m)<epsilon && 
+            Math.abs(n-z.n)<epsilon && 
+            Math.abs(o-z.o)<epsilon && 
+            Math.abs(p-z.p)<epsilon;
     }
 
     public int hashCode() {
@@ -272,4 +317,72 @@ public class Matrix {
             Float.floatToIntBits(o) ^
             Float.floatToIntBits(p);
     }
+
+    /** retrieves the current objectspace-to-windowspace projection matrix from the GL context */
+    public static Matrix getProjectionMatrix(GL gl) {
+        int view[] = new int[4];
+        double mvmatrix[] = new double[16];
+        double projmatrix[] = new double[16];
+        gl.glGetIntegerv(GL.GL_VIEWPORT, view, 0);
+        gl.glGetDoublev(GL.GL_MODELVIEW_MATRIX, mvmatrix, 0);
+        gl.glGetDoublev(GL.GL_PROJECTION_MATRIX, projmatrix, 0);
+        Matrix m = new Matrix(mvmatrix[0],
+                              mvmatrix[4],
+                              mvmatrix[8],
+                              mvmatrix[12],
+                              mvmatrix[1],
+                              mvmatrix[5],
+                              mvmatrix[9],
+                              mvmatrix[13],
+                              mvmatrix[2],
+                              mvmatrix[6],
+                              mvmatrix[10],
+                              mvmatrix[14],
+                              mvmatrix[3],
+                              mvmatrix[7],
+                              mvmatrix[11],
+                              mvmatrix[15]);
+        Matrix p = new Matrix(projmatrix[0],
+                              projmatrix[4],
+                              projmatrix[8],
+                              projmatrix[12],
+                              projmatrix[1],
+                              projmatrix[5],
+                              projmatrix[9],
+                              projmatrix[13],
+                              projmatrix[2],
+                              projmatrix[6],
+                              projmatrix[10],
+                              projmatrix[14],
+                              projmatrix[3],
+                              projmatrix[7],
+                              projmatrix[11],
+                              projmatrix[15]);
+        Matrix z =
+            new Matrix(0.5*view[2], 0,           0,   view[0]+view[2]*0.5,
+                       0,           0.5*view[3], 0,   view[1]+view[3]*0.5,
+                       0,           0,           0.5, 0.5,
+                       0,           0,           0,   1);
+        return z.times(p).times(m);
+    }
+
+    /** returns the constraint-conjunction "(forall v)Mv=v" */
+    public AffineConstraint getAffineConstraint(float epsilon) {
+        AffineConstraint c1 = getConstraint(a-1, b,   c,   d  , epsilon);
+        AffineConstraint c2 = getConstraint(e,   f-1, g,   h  , epsilon);
+        AffineConstraint c3 = getConstraint(i,   j,   k-1, l  , epsilon);
+        // FIXME: bottom row constraint?
+        return c1.intersect(c2.intersect(c3, epsilon), epsilon);
+    }
+
+    /** the AffineConstraint representing ax+by+cz+d=0 */
+    private static AffineConstraint getConstraint(float a, float b, float c, float d, float epsilon) {
+        a = Math.abs(a) <= epsilon ? 0 : a;
+        b = Math.abs(b) <= epsilon ? 0 : b;
+        c = Math.abs(c) <= epsilon ? 0 : c;
+        d = Math.abs(d) <= epsilon ? 0 : d;
+        if (a!=0 || b!=0 || c!=0) return new Plane(a, b, c, d);
+        if (d==0) return new AffineConstraint.All();
+        return new AffineConstraint.Nothing();
+    }
 }