src/edu/berkeley/qfat/geom/Matrix.java

   1 package edu.berkeley.qfat.geom;
   2
   3 /** an affine matrix; immutable */
   4 public class Matrix {
   5
   6     /**
   7      *  <pre>
   8      *  [ a b c d ]   [ x ]
   9      *  [ e f g h ]   [ y ]
  10      *  [ i j k l ]   [ z ]
  11      *  [ m n o p ]   [ 1 ]
  12      *  </pre>
  13      */
  14     public final float a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p;
  15
  16     /** the zero matrix */
  17     public static final Matrix ZERO = new Matrix(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
  18
  19     /** the identity matrix */
  20     public static final Matrix ONE  = new Matrix(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1);
  21
  22     /** a scaling matrix (uniform in all dimensions) */
  23     public static Matrix scale(float scale) { return scale(scale, scale, scale); }
  24
  25     /** a scaling matrix */
  26     public static Matrix scale(float scalex, float scaley, float scalez) {
  27         return new Matrix(scalex, 0, 0, 0, 0, scaley, 0, 0, 0, 0, scalez, 0, 0, 0, 0, 1); }
  28
  29     public Matrix(float a, float b, float c, float d, float e, float f, float g,
  30                   float h, float i, float j, float k, float l, float m, float n, float o, float p) {
  31         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;
  32         this.j = j; this.k = k; this.l = l; this.m = m; this.n = n; this.o = o; this.p = p;
  33     }
  34
  35     public static Matrix translate(Vec translate) {
  36         return new Matrix(1, 0, 0, translate.x,
  37                           0, 1, 0, translate.y,
  38                           0, 0, 1, translate.z,
  39                           0, 0, 0, 1);
  40     }
  41
  42     public Matrix plus(Matrix x) {
  43         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);
  44     }
  45     public Matrix minus(Matrix x) {
  46         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);
  47     }
  48     public Matrix times(float x) {
  49         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);
  50     }
  51
  52     /** computes (v^T)(this)(v) */
  53     public float preAndPostMultiply(Point point) {
  54         float ret =
  55             ((a*point.x + b*point.y + c*point.z + d) * point.x) +
  56             ((e*point.x + f*point.y + g*point.z + h) * point.y) +
  57             ((i*point.x + j*point.y + k*point.z + l) * point.z) +
  58             ((m*point.x + n*point.y + o*point.z + p) * 1);
  59         return ret;
  60     }
  61
  62     public static Matrix rotate(Vec axis, float angle) {
  63         float q = (float)Math.cos(angle);
  64         float s = (float)Math.sin(angle);
  65         float t = 1 - q;
  66         float a = (float)(q + axis.x*axis.x*t);
  67         float f = (float)(q + axis.y*axis.y*t);
  68         float k = (float)(q + axis.z*axis.z*t);
  69         float tmp1 = axis.x*axis.y*t;
  70         float tmp2 = axis.z*s;
  71         float e = (float)(tmp1 + tmp2);
  72         float b = (float)(tmp1 - tmp2);
  73         tmp1 = axis.x*axis.z*t;
  74         tmp2 = axis.y*s;
  75         float i = (float)(tmp1 - tmp2);
  76         float c = (float)(tmp1 + tmp2);
  77         tmp1 = axis.y*axis.z*t;
  78         tmp2 = axis.x*s;
  79         float j = (float)(tmp1 + tmp2);
  80         float g = (float)(tmp1 - tmp2);
  81         return new Matrix(a, b, c, 0,
  82                           e, f, g, 0,
  83                           i, j, k, 0,
  84                           0, 0, 0, 1);
  85     }
  86     public Point times(Point p) {
  87         // discards bottom row
  88         return new Point(a*p.x + b*p.y + c*p.z + d,
  89                          e*p.x + f*p.y + g*p.z + h,
  90                          i*p.x + j*p.y + k*p.z + l);
  91     }
  92     public Vec times(Vec p) {
  93         // discards bottom row
  94         return new Vec(a*p.x + b*p.y + c*p.z + d,
  95                        e*p.x + f*p.y + g*p.z + h,
  96                        i*p.x + j*p.y + k*p.z + l);
  97     }
  98     public Point apply(Point p) { return p; }
  99     public Vec apply(Vec v) { return v; }
 100     public Matrix invert() { return this; }
 101     public Matrix times(Matrix m) { return this; }
 102
 103
 104     public float determinant() {
 105         float m00 = a;
 106         float m01 = b;
 107         float m02 = c;
 108         float m03 = d;
 109         float m10 = e;
 110         float m11 = f;
 111         float m12 = g;
 112         float m13 = h;
 113         float m20 = i;
 114         float m21 = j;
 115         float m22 = k;
 116         float m23 = l;
 117         float m30 = m;
 118         float m31 = n;
 119         float m32 = o;
 120         float m33 = p;
 121         return
 122             m03 * m12 * m21 * m30-m02 * m13 * m21 * m30-m03 * m11 * m22 * m30+m01 * m13    * m22 * m30+
 123             m02 * m11 * m23 * m30-m01 * m12 * m23 * m30-m03 * m12 * m20 * m31+m02 * m13    * m20 * m31+
 124             m03 * m10 * m22 * m31-m00 * m13 * m22 * m31-m02 * m10 * m23 * m31+m00 * m12    * m23 * m31+
 125             m03 * m11 * m20 * m32-m01 * m13 * m20 * m32-m03 * m10 * m21 * m32+m00 * m13    * m21 * m32+
 126             m01 * m10 * m23 * m32-m00 * m11 * m23 * m32-m02 * m11 * m20 * m33+m01 * m12    * m20 * m33+
 127             m02 * m10 * m21 * m33-m00 * m12 * m21 * m33-m01 * m10 * m22 * m33+m00 * m11    * m22 * m33;
 128     }
 129
 130     public Matrix inverse() {
 131         float m00 = a;
 132         float m01 = b;
 133         float m02 = c;
 134         float m03 = d;
 135         float m10 = e;
 136         float m11 = f;
 137         float m12 = g;
 138         float m13 = h;
 139         float m20 = i;
 140         float m21 = j;
 141         float m22 = k;
 142         float m23 = l;
 143         float m30 = m;
 144         float m31 = n;
 145         float m32 = o;
 146         float m33 = p;
 147         return
 148             new Matrix(m12*m23*m31 - m13*m22*m31 + m13*m21*m32 - m11*m23*m32 - m12*m21*m33 + m11*m22*m33,
 149                        m03*m22*m31 - m02*m23*m31 - m03*m21*m32 + m01*m23*m32 + m02*m21*m33 - m01*m22*m33,
 150                        m02*m13*m31 - m03*m12*m31 + m03*m11*m32 - m01*m13*m32 - m02*m11*m33 + m01*m12*m33,
 151                        m03*m12*m21 - m02*m13*m21 - m03*m11*m22 + m01*m13*m22 + m02*m11*m23 - m01*m12*m23,
 152                        m13*m22*m30 - m12*m23*m30 - m13*m20*m32 + m10*m23*m32 + m12*m20*m33 - m10*m22*m33,
 153                        m02*m23*m30 - m03*m22*m30 + m03*m20*m32 - m00*m23*m32 - m02*m20*m33 + m00*m22*m33,
 154                        m03*m12*m30 - m02*m13*m30 - m03*m10*m32 + m00*m13*m32 + m02*m10*m33 - m00*m12*m33,
 155                        m02*m13*m20 - m03*m12*m20 + m03*m10*m22 - m00*m13*m22 - m02*m10*m23 + m00*m12*m23,
 156                        m11*m23*m30 - m13*m21*m30 + m13*m20*m31 - m10*m23*m31 - m11*m20*m33 + m10*m21*m33,
 157                        m03*m21*m30 - m01*m23*m30 - m03*m20*m31 + m00*m23*m31 + m01*m20*m33 - m00*m21*m33,
 158                        m01*m13*m30 - m03*m11*m30 + m03*m10*m31 - m00*m13*m31 - m01*m10*m33 + m00*m11*m33,
 159                        m03*m11*m20 - m01*m13*m20 - m03*m10*m21 + m00*m13*m21 + m01*m10*m23 - m00*m11*m23,
 160                        m12*m21*m30 - m11*m22*m30 - m12*m20*m31 + m10*m22*m31 + m11*m20*m32 - m10*m21*m32,
 161                        m01*m22*m30 - m02*m21*m30 + m02*m20*m31 - m00*m22*m31 - m01*m20*m32 + m00*m21*m32,
 162                        m02*m11*m30 - m01*m12*m30 - m02*m10*m31 + m00*m12*m31 + m01*m10*m32 - m00*m11*m32,
 163                        m01*m12*m20 - m02*m11*m20 + m02*m10*m21 - m00*m12*m21 - m01*m10*m22 + m00*m11*m22)
 164             .times(1/determinant());
 165     }
 166 }