001/*
002 * Renderer 9. The MIT License.
003 * Copyright (c) 2022 rlkraft@pnw.edu
004 * See LICENSE for details.
005*/
006
007package renderer.scene;
008
009/**
010   A {@code Matrix} object has four {@link Vector} objects.
011<p>
012   The four {@link Vector} objects represent the four column vectors
013   of the 4-by-4 matrix (as in a Linear Algebra course).
014<p>
015   In computer graphics, the points and vectors of 3-dimensional space
016   are represented using 4-dimensional homogeneous coordinates.
017   So each transformation of 3-dimensional space is represented by
018   a 4-by-4 (homogeneous) matrix.
019<p>
020   A 4-by-4 matrix represents a transformation of 3-dimensional space.
021   The most common transformations are translation, rotation, and
022   scaling. A 4-by-4 matrix can also represent a projection transformation.
023*/
024public final class Matrix
025{
026   public final Vector v1, v2, v3, v4; // these are column vectors
027
028   /**
029      Construct an arbitrary 4-by-4 {@code Matrix}
030      with the given column {@link Vector}s.
031      <p>
032      Notice that this is a private constructor. Other
033      objects should use the static facory methods to
034      create new {@code Matrix} objects.
035
036      @param v1  1st column {@link Vector} for the new {@code Matrix}
037      @param v2  2nd column {@link Vector} for the new {@code Matrix}
038      @param v3  3rd column {@link Vector} for the new {@code Matrix}
039      @param v4  4th column {@link Vector} for the new {@code Matrix}
040      @return a new {@code Matrix} object
041   */
042   private Matrix(final Vector v1, final Vector v2,
043                  final Vector v3, final Vector v4)
044   {
045      this.v1 = v1;  // Notice that we are not making
046      this.v2 = v2;  // copies of the column vectors,
047      this.v3 = v3;  // We are just making references
048      this.v4 = v4;  // to them.
049   }
050
051
052   /**
053      This is a static facory method.
054      <p>
055      Construct an arbitrary 4-by-4 {@code Matrix}
056      using the given column {@link Vector}s.
057
058      @param c1  1st column {@link Vector} for the new {@code Matrix}
059      @param c2  2nd column {@link Vector} for the new {@code Matrix}
060      @param c3  3rd column {@link Vector} for the new {@code Matrix}
061      @param c4  4th column {@link Vector} for the new {@code Matrix}
062      @return a new {@code Matrix} object
063   */
064   public static Matrix buildFromColumns(final Vector c1, final Vector c2,
065                                         final Vector c3, final Vector c4)
066   {
067      return new Matrix(c1, c2, c3, c4);
068   }
069
070
071   /**
072      This is a static facory method.
073      <p>
074      Construct an arbitrary 4-by-4 {@code Matrix}
075      using the given row {@link Vector}s.
076
077      @param r1  1st row {@link Vector} for the new {@code Matrix}
078      @param r2  2nd row {@link Vector} for the new {@code Matrix}
079      @param r3  3rd row {@link Vector} for the new {@code Matrix}
080      @param r4  4th row {@link Vector} for the new {@code Matrix}
081      @return a new {@code Matrix} object
082   */
083   public static Matrix buildFromRows(final Vector r1, final Vector r2,
084                                      final Vector r3, final Vector r4)
085   {
086      Vector c1 = new Vector(r1.x, r2.x, r3.x, r4.x);
087      Vector c2 = new Vector(r1.y, r2.y, r3.y, r4.y);
088      Vector c3 = new Vector(r1.z, r2.z, r3.z, r4.z);
089      Vector c4 = new Vector(r1.w, r2.w, r3.w, r4.w);
090      return new Matrix(c1, c2, c3, c4);
091   }
092
093
094   /**
095      This is a static facory method.
096      <p>
097      Construct an identity {@code Matrix}.
098
099      @return a new {@code Matrix} object containing an identity {@code Matrix}
100   */
101   public static Matrix identity()
102   {
103      return scale(1.0, 1.0, 1.0);
104   }
105
106
107   /**
108      This is a static facory method.
109      <p>
110      Construct a translation {@code Matrix} that translates by the
111      given amounts in the {@code x}, {@code y}, and {@code z} directions..
112
113      @param x  translation factor for the x-direction
114      @param y  translation factor for the y-direction
115      @param z  translation factor for the z-direction
116      @return a new {@code Matrix} object containing a translation {@code Matrix}
117   */
118   public static Matrix translate(final double x, final double y, final double z)
119   {
120      return new Matrix(new Vector(1.0, 0.0, 0.0, 0.0),
121                        new Vector(0.0, 1.0, 0.0, 0.0),
122                        new Vector(0.0, 0.0, 1.0, 0.0),
123                        new Vector(  x,   y,   z, 1.0));
124   }
125
126
127   /**
128      This is a static facory method.
129      <p>
130      Construct a diagonal {@code Matrix} with the given number
131      on the diagonal.
132      <p>
133      This is also a uniform scaling matrix.
134
135      @param d  the diagonal value for the new {@code Matrix}
136      @return a new {@code Matrix} object containing a scaling {@code Matrix}
137   */
138   public static Matrix scale(final double d)
139   {
140      return scale(d, d, d);
141   }
142
143
144   /**
145      This is a static facory method.
146      <p>
147      Construct a (diagonal) {@code Matrix} that scales in
148      the x, y, and z directions by the given factors.
149
150      @param x  scale factor for the x-direction
151      @param y  scale factor for the y-direction
152      @param z  scale factor for the z-direction
153      @return a new {@code Matrix} object containing a scaling {@code Matrix}
154   */
155   public static Matrix scale(final double x, final double y, final double z)
156   {
157      return new Matrix(new Vector(  x, 0.0, 0.0, 0.0),
158                        new Vector(0.0,   y, 0.0, 0.0),
159                        new Vector(0.0, 0.0,   z, 0.0),
160                        new Vector(0.0, 0.0, 0.0, 1.0));
161   }
162
163
164   /**
165      This is a static facory method.
166      <p>
167      Construct a rotation {@code Matrix} that rotates around
168      the x-axis by the angle {@code theta}.
169
170      @param theta  angle (in degrees) to rotate by around the x-axis
171      @return a new {@code Matrix} object containing a rotation {@code Matrix}
172   */
173   public static Matrix rotateX(final double theta)
174   {
175      return rotate(theta, 1,0,0);
176   }
177
178
179   /**
180      This is a static facory method.
181      <p>
182      Construct a rotation {@code Matrix} that rotates around
183      the y-axis by the angle {@code theta}.
184
185      @param theta  angle (in degrees) to rotate by around the y-axis
186      @return a new {@code Matrix} object containing a rotation {@code Matrix}
187   */
188   public static Matrix rotateY(final double theta)
189   {
190      return rotate(theta, 0,1,0);
191   }
192
193
194   /**
195      This is a static facory method.
196      <p>
197      Construct a rotation {@code Matrix} that rotates around
198      the z-axis by the angle {@code theta}.
199
200      @param theta  angle (in degrees) to rotate by around the z-axis
201      @return a new {@code Matrix} object containing a rotation {@code Matrix}
202   */
203   public static Matrix rotateZ(final double theta)
204   {
205      return rotate(theta, 0,0,1);
206   }
207
208
209   /**
210      This is a static facory method.
211      <p>
212      Construct a rotation {@code Matrix} that rotates around
213      the axis vector {@code (x,y,z)} by the angle {@code theta}.
214      <p>
215      See
216      <a href="https://www.opengl.org/sdk/docs/man2/xhtml/glRotate.xml" target="_top">
217               https://www.opengl.org/sdk/docs/man2/xhtml/glRotate.xml</a>
218
219      @param theta  angle (in degrees) to rotate by around the axis vector
220      @param x      x-component of the axis vector for the rotation
221      @param y      y-component of the axis vector for the rotation
222      @param z      z-component of the axis vector for the rotation
223      @return a new {@code Matrix} object containing a rotation {@code Matrix}
224   */
225   public static Matrix rotate(final double theta, final double x, final double y, final double z)
226   {
227      final double norm = Math.sqrt(x*x + y*y + z*z);
228      final double ux = x/norm;
229      final double uy = y/norm;
230      final double uz = z/norm;
231
232      final double c = Math.cos( (Math.PI/180.0)*theta );
233      final double s = Math.sin( (Math.PI/180.0)*theta );
234
235      return new Matrix(
236        new Vector(ux*ux*(1-c)+c,      uy*ux*(1-c)+(uz*s), uz*ux*(1-c)-(uy*s), 0.0),
237        new Vector(ux*uy*(1-c)-(uz*s), uy*uy*(1-c)+c,      uz*uy*(1-c)+(ux*s), 0.0),
238        new Vector(ux*uz*(1-c)+(uy*s), uy*uz*(1-c)-(ux*s), uz*uz*(1-c)+c,      0.0),
239        new Vector(0.0,                0.0,                0.0,                1.0));
240   }
241
242
243   /**
244      A scalar times this {@code Matrix} returns a new {@code Matrix}.
245
246      @param s  scalar value to multiply this {@code Matrix} by
247      @return a new {@code Matrix} object containing the scalar s times this {@code Matrix}
248   */
249   public Matrix times(final double s) // return s * this
250   {
251      return new Matrix(v1.times(s), v2.times(s), v3.times(s), v4.times(s));
252   }
253
254
255   /**
256      This {@code Matrix} times a {@link Vector} returns a new {@link Vector}.
257
258      @param v  {@link Vector} to be multiplied by this {@code Matrix}
259      @return new {@link Vector} object containing this {@code Matrix} times the {@link Vector} v
260   */
261   public Vector times(final Vector v) // return this * v
262   {
263      /*
264      return v1.times(v.x).plus(v2.times(v.y).plus(v3.times(v.z).plus(v4.times(v.w))));
265      */
266      /*
267      // Here is what this works out to be.
268      final Vector v1x = this.v1.times(v.x);
269      final Vector v2y = this.v2.times(v.y);
270      final Vector v3z = this.v3.times(v.z);
271      final Vector v4w = this.v4.times(v.w);
272      final Vector sum1 = v1x.plus(v2y);
273      final Vector sum2 = sum1.plus(v3z);
274      final Vector sum3 = sum2.plus(v4w);
275      return sum3;
276      */
277      // dot product of each row of this matrix with the vector v
278      final double x = (v1.x * v.x) + (v2.x * v.y) + (v3.x * v.z) + (v4.x * v.w);
279      final double y = (v1.y * v.x) + (v2.y * v.y) + (v3.y * v.z) + (v4.y * v.w);
280      final double z = (v1.z * v.x) + (v2.z * v.y) + (v3.z * v.z) + (v4.z * v.w);
281      final double w = (v1.w * v.x) + (v2.w * v.y) + (v3.w * v.z) + (v4.w * v.w);
282      return new Vector(x, y, z, w);
283   }
284
285
286   /**
287      This {@code Matrix} times {@code Matrix} {@code m} returns a new {@code Matrix}.
288
289      @param m  {@code Matrix} value to be multiplied on the right of this {@code Matrix}
290      @return new {@code Matrix} object containing this {@code Matrix} times {@code Matrix} {@code m}
291   */
292   public Matrix times(final Matrix m) // return this * m
293   {
294      return new Matrix(this.times(m.v1),   // 1st column vector of the result
295                        this.times(m.v2),   // 2nd column vector of the result
296                        this.times(m.v3),   // 3rd column vector of the result
297                        this.times(m.v4) ); // 4th column vector of the result
298   }
299
300
301   /**
302      This {@code Matrix} times a {@link Vertex} returns a new {@link Vertex}.
303
304      @param v  {@link Vertex} to be multiplied by this {@code Matrix}
305      @return new {@link Vertex} object containing this {@code Matrix} times the {@link Vertex} v
306   */
307   public Vertex times(final Vertex v) // return this * v
308   {
309      /*
310      final Vector v = v1.times(v.x).plus(v2.times(v.y).plus(v3.times(v.z).plus(v4.times(v.w))));
311      return new Vertex(v.x, v.y, v.z, v.w);
312      */
313      // dot product of each row of this matrix with the vertex v
314      final double x = (v1.x * v.x) + (v2.x * v.y)  + (v3.x * v.z) + (v4.x * v.w);
315      final double y = (v1.y * v.x) + (v2.y * v.y)  + (v3.y * v.z) + (v4.y * v.w);
316      final double z = (v1.z * v.x) + (v2.z * v.y)  + (v3.z * v.z) + (v4.z * v.w);
317      final double w = (v1.w * v.x) + (v2.w * v.y)  + (v3.w * v.z) + (v4.w * v.w);
318      return new Vertex(x, y, z, w);
319   }
320
321
322   /**
323      Assuming that the 3-by-3 "rotation part" of this 4-by-4
324      {@code Matrix} represents a pure rotation, return the
325      rotation's three Euler angles, in radians, in the
326      order {@code [x, y, z]} for rotations in the order
327      {@code R_z * R_y * R_x}.
328   <p>
329      A 3-by-3 matrix is a rotation matrix if its inverse is
330      equal to its transpose and its determinant is equal to 1.
331   <p>
332      See <a href="http://eecs.qmul.ac.uk/~gslabaugh/publications/euler.pdf" target="_top">
333                   http://eecs.qmul.ac.uk/~gslabaugh/publications/euler.pdf</a>
334
335      @return an array of 3 doubles which are this rotation's Euler angles in radians
336   */
337   public double[] rot2euler()
338   {
339      final double r_11 = v1.x,
340                   r_12 = v2.x,
341                   r_13 = v3.x,
342                   r_21 = v1.y,
343                   r_31 = v1.z,
344                   r_32 = v2.z,
345                   r_33 = v3.z;
346
347      final double r_x,
348                   r_y,
349                   r_z;
350
351      if (r_31 != 1.0 && r_31 != -1.0)
352      {
353         r_y = -Math.asin(r_31);
354         r_x = Math.atan2(r_32 / Math.cos(r_y),
355                          r_33 / Math.cos(r_y));
356         r_z = Math.atan2(r_21 / Math.cos(r_y),
357                          r_11 / Math.cos(r_y));
358      }
359      else
360      {
361         if (r_31 == -1.0)
362         {
363            r_y = Math.PI / 2.0;
364            r_x = Math.atan2(r_12, r_13);
365         }
366         else // r_31 == 1.0
367         {
368            r_y = -Math.PI / 2.0;
369            r_x = Math.atan2(-r_12, -r_13);
370         }
371         r_z = 0.0;
372      }
373
374      return new double[]{r_x, r_y, r_z};
375   }
376
377
378   /**
379      Assuming that this {@code Matrix} represents a 3D rotation,
380      return the rotation matrix formed by multiplying this matrix's
381      three Euler angle rotations in the order {@code R_z * R_y * R_x}.
382      <p>
383      This is mainly for debugging. If this matrix is really a pure
384      rotation, then this method will return a copy of this matrix.
385
386      @return the "eulerized" version of this {@code Matrix}
387   */
388   public Matrix eulerize()
389   {
390      double[] euler = this.rot2euler();
391      return Matrix.rotateZ(euler[2]*(180.0/Math.PI)).times(
392             Matrix.rotateY(euler[1]*(180.0/Math.PI)).times(
393             Matrix.rotateX(euler[0]*(180.0/Math.PI))));
394   }
395
396
397   /**
398      For debugging.
399
400      @return {@link String} representation of this {@code Matrix} object
401   */
402   @Override
403   public String toString()
404   {
405      String result = "";
406      final int p = 5;      // the precision for the following format string
407      final int w = p + 4;  // the width for the following format string
408      final String format = "% "+w+"."+p+"f  % "+w+"."+p+"f  % "+w+"."+p+"f  % "+w+"."+p+"f";
409      result += String.format("[[" + format + " ]\n",  v1.x, v2.x, v3.x, v4.x);
410      result += String.format(" [" + format + " ]\n",  v1.y, v2.y, v3.y, v4.y);
411      result += String.format(" [" + format + " ]\n",  v1.z, v2.z, v3.z, v4.z);
412      result += String.format(" [" + format + " ]]",   v1.w, v2.w, v3.w, v4.w);
413    //result += String.format("[[% .5f  % .5f  % .5f  % .5f ]\n",  v1.x, v2.x, v3.x, v4.x);
414    //result += String.format(" [% .5f  % .5f  % .5f  % .5f ]\n",  v1.y, v2.y, v3.y, v4.y);
415    //result += String.format(" [% .5f  % .5f  % .5f  % .5f ]\n",  v1.z, v2.z, v3.z, v4.z);
416    //result += String.format(" [% .5f  % .5f  % .5f  % .5f ]]",   v1.w, v2.w, v3.w, v4.w);
417      return result;
418   }
419}