1 package edu.berkeley.qfat.geom;
3 // FEATURE: precompute/cache determinant?
5 /** an affine matrix; immutable */
16 public final float a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p;
18 /** the zero matrix */
19 public static final Matrix ZERO = new Matrix(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
21 /** the identity matrix */
22 public static final Matrix ONE = new Matrix(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1);
24 public Matrix(float a, float b, float c, float d, float e, float f, float g,
25 float h, float i, float j, float k, float l, float m, float n, float o, float p) {
26 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;
27 this.j = j; this.k = k; this.l = l; this.m = m; this.n = n; this.o = o; this.p = p;
30 /** a scaling matrix (uniform in all dimensions) */
31 public static Matrix scale(float scale) { return scale(scale, scale, scale); }
33 /** a scaling matrix */
34 public static Matrix scale(float scalex, float scaley, float scalez) {
35 return new Matrix(scalex, 0, 0, 0, 0, scaley, 0, 0, 0, 0, scalez, 0, 0, 0, 0, 1); }
37 /** a translation matrix */
38 public static Matrix translate(Vec translate) {
39 return new Matrix(1, 0, 0, translate.x,
45 /** a rotation matrix, <tt>angle</tt> degrees around <tt>axis</tt> */
46 public static Matrix rotate(Vec axis, float angle) {
47 float q = (float)Math.cos(angle);
48 float s = (float)Math.sin(angle);
50 float a = (float)(q + axis.x*axis.x*t);
51 float f = (float)(q + axis.y*axis.y*t);
52 float k = (float)(q + axis.z*axis.z*t);
53 float tmp1 = axis.x*axis.y*t;
54 float tmp2 = axis.z*s;
55 float e = (float)(tmp1 + tmp2);
56 float b = (float)(tmp1 - tmp2);
57 tmp1 = axis.x*axis.z*t;
59 float i = (float)(tmp1 - tmp2);
60 float c = (float)(tmp1 + tmp2);
61 tmp1 = axis.y*axis.z*t;
63 float j = (float)(tmp1 + tmp2);
64 float g = (float)(tmp1 - tmp2);
65 return new Matrix(a, b, c, 0,
71 public Matrix plus(Matrix x) {
72 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);
74 public Matrix minus(Matrix x) {
75 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);
77 public Matrix times(float x) {
78 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);
81 /** computes (v^T)(this)(v) */
82 public float preAndPostMultiply(Point point) {
84 ((a*point.x + b*point.y + c*point.z + d) * point.x) +
85 ((e*point.x + f*point.y + g*point.z + h) * point.y) +
86 ((i*point.x + j*point.y + k*point.z + l) * point.z) +
87 ((m*point.x + n*point.y + o*point.z + p) * 1);
91 /** discards bottom row */
92 public Point times(Point p) {
93 // discards bottom row
94 return new Point(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);
99 /** discards bottom row */
100 public Vec times(Vec p) {
101 return new Vec(a*p.x + b*p.y + c*p.z + d,
102 e*p.x + f*p.y + g*p.z + h,
103 i*p.x + j*p.y + k*p.z + l);
106 /** multiply by another matrix */
107 public Matrix times(Matrix z) {
140 return new Matrix(t00*m00 + t01*m10 + t02*m20 + t03*m30,
141 t00*m01 + t01*m11 + t02*m21 + t03*m31,
142 t00*m02 + t01*m12 + t02*m22 + t03*m32,
143 t00*m03 + t01*m13 + t02*m23 + t03*m33,
144 t10*m00 + t11*m10 + t12*m20 + t13*m30,
145 t10*m01 + t11*m11 + t12*m21 + t13*m31,
146 t10*m02 + t11*m12 + t12*m22 + t13*m32,
147 t10*m03 + t11*m13 + t12*m23 + t13*m33,
148 t20*m00 + t21*m10 + t22*m20 + t23*m30,
149 t20*m01 + t21*m11 + t22*m21 + t23*m31,
150 t20*m02 + t21*m12 + t22*m22 + t23*m32,
151 t20*m03 + t21*m13 + t22*m23 + t23*m33,
152 t30*m00 + t31*m10 + t32*m20 + t33*m30,
153 t30*m01 + t31*m11 + t32*m21 + t33*m31,
154 t30*m02 + t31*m12 + t32*m22 + t33*m32,
155 t30*m03 + t31*m13 + t32*m23 + t33*m33);
158 /** compute the determinant of the matrix */
159 public float determinant() {
177 m03 * m12 * m21 * m30-m02 * m13 * m21 * m30-m03 * m11 * m22 * m30+m01 * m13 * m22 * m30+
178 m02 * m11 * m23 * m30-m01 * m12 * m23 * m30-m03 * m12 * m20 * m31+m02 * m13 * m20 * m31+
179 m03 * m10 * m22 * m31-m00 * m13 * m22 * m31-m02 * m10 * m23 * m31+m00 * m12 * m23 * m31+
180 m03 * m11 * m20 * m32-m01 * m13 * m20 * m32-m03 * m10 * m21 * m32+m00 * m13 * m21 * m32+
181 m01 * m10 * m23 * m32-m00 * m11 * m23 * m32-m02 * m11 * m20 * m33+m01 * m12 * m20 * m33+
182 m02 * m10 * m21 * m33-m00 * m12 * m21 * m33-m01 * m10 * m22 * m33+m00 * m11 * m22 * m33;
185 /** compute the inverse of the matrix, returns null if not invertible */
186 public Matrix inverse() {
187 float determinant = determinant();
188 if (determinant==0) return null;
206 new Matrix(m12*m23*m31 - m13*m22*m31 + m13*m21*m32 - m11*m23*m32 - m12*m21*m33 + m11*m22*m33,
207 m03*m22*m31 - m02*m23*m31 - m03*m21*m32 + m01*m23*m32 + m02*m21*m33 - m01*m22*m33,
208 m02*m13*m31 - m03*m12*m31 + m03*m11*m32 - m01*m13*m32 - m02*m11*m33 + m01*m12*m33,
209 m03*m12*m21 - m02*m13*m21 - m03*m11*m22 + m01*m13*m22 + m02*m11*m23 - m01*m12*m23,
210 m13*m22*m30 - m12*m23*m30 - m13*m20*m32 + m10*m23*m32 + m12*m20*m33 - m10*m22*m33,
211 m02*m23*m30 - m03*m22*m30 + m03*m20*m32 - m00*m23*m32 - m02*m20*m33 + m00*m22*m33,
212 m03*m12*m30 - m02*m13*m30 - m03*m10*m32 + m00*m13*m32 + m02*m10*m33 - m00*m12*m33,
213 m02*m13*m20 - m03*m12*m20 + m03*m10*m22 - m00*m13*m22 - m02*m10*m23 + m00*m12*m23,
214 m11*m23*m30 - m13*m21*m30 + m13*m20*m31 - m10*m23*m31 - m11*m20*m33 + m10*m21*m33,
215 m03*m21*m30 - m01*m23*m30 - m03*m20*m31 + m00*m23*m31 + m01*m20*m33 - m00*m21*m33,
216 m01*m13*m30 - m03*m11*m30 + m03*m10*m31 - m00*m13*m31 - m01*m10*m33 + m00*m11*m33,
217 m03*m11*m20 - m01*m13*m20 - m03*m10*m21 + m00*m13*m21 + m01*m10*m23 - m00*m11*m23,
218 m12*m21*m30 - m11*m22*m30 - m12*m20*m31 + m10*m22*m31 + m11*m20*m32 - m10*m21*m32,
219 m01*m22*m30 - m02*m21*m30 + m02*m20*m31 - m00*m22*m31 - m01*m20*m32 + m00*m21*m32,
220 m02*m11*m30 - m01*m12*m30 - m02*m10*m31 + m00*m12*m31 + m01*m10*m32 - m00*m11*m32,
221 m01*m12*m20 - m02*m11*m20 + m02*m10*m21 - m00*m12*m21 - m01*m10*m22 + m00*m11*m22)
222 .times(1/determinant);