optimizations to _math_matrix_rotate() (Rudolf Opalla)

1 files changed, 155 insertions, 97 deletions

diff --git a/src/mesa/math/m_matrix.c b/src/mesa/math/m_matrix.c
index 2a69336b4d..95b6ebed56 100644
--- a/src/mesa/math/m_matrix.c
+++ b/src/mesa/math/m_matrix.c
@@ -1,8 +1,8 @@
-/* $Id: m_matrix.c,v 1.12 2002/06/29 19:48:17 brianp Exp $ */
+/* $Id: m_matrix.c,v 1.13 2002/09/12 16:26:04 brianp Exp $ */
 
 /*
  * Mesa 3-D graphics library
- * Version:  4.0.2
+ * Version:  4.1
  *
  * Copyright (C) 1999-2002  Brian Paul   All Rights Reserved.
  *
@@ -539,123 +539,181 @@ static GLboolean matrix_invert( GLmatrix *mat )
 /*
  * Generate a 4x4 transformation matrix from glRotate parameters, and
  * postmultiply the input matrix by it.
+ * This function contributed by Erich Boleyn (erich@uruk.org).
+ * Optimizatios contributed by Rudolf Opalla (rudi@khm.de).
  */
 void
 _math_matrix_rotate( GLmatrix *mat,
 		     GLfloat angle, GLfloat x, GLfloat y, GLfloat z )
 {
-   /* This function contributed by Erich Boleyn (erich@uruk.org) */
-   GLfloat mag, s, c;
-   GLfloat xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c;
+   GLfloat xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c, s, c;
    GLfloat m[16];
+   GLboolean optimized;
 
    s = (GLfloat) sin( angle * DEG2RAD );
    c = (GLfloat) cos( angle * DEG2RAD );
 
-   mag = (GLfloat) GL_SQRT( x*x + y*y + z*z );
+   MEMCPY(m, Identity, sizeof(GLfloat)*16);
+   optimized = GL_FALSE;
 
-   if (mag <= 1.0e-4) {
-      /* generate an identity matrix and return */
-      MEMCPY(m, Identity, sizeof(GLfloat)*16);
-      return;
-   }
+#define M(row,col)  m[col*4+row]
 
-   x /= mag;
-   y /= mag;
-   z /= mag;
+   if (x == 0.0F) {
+      if (y == 0.0F) {
+         if (z != 0.0F) {
+            optimized = GL_TRUE;
+            /* rotate only around z-axis */
+            M(0,0) = c;
+            M(1,1) = c;
+            if (z < 0.0F) {
+               M(0,1) = s;
+               M(1,0) = -s;
+            }
+            else {
+               M(0,1) = -s;
+               M(1,0) = s;
+            }
+         }
+      }
+      else if (z == 0.0F) {
+         optimized = GL_TRUE;
+         /* rotate only around y-axis */
+         M(0,0) = c;
+         M(2,2) = c;
+         if (y < 0.0F) {
+            M(0,2) = -s;
+            M(2,0) = s;
+         }
+         else {
+            M(0,2) = s;
+            M(2,0) = -s;
+         }
+      }
+   }
+   else if (y == 0.0F) {
+      if (z == 0.0F) {
+         optimized = GL_TRUE;
+         /* rotate only around x-axis */
+         M(1,1) = c;
+         M(2,2) = c;
+         if (y < 0.0F) {
+            M(1,2) = s;
+            M(2,1) = -s;
+         }
+         else {
+            M(1,2) = -s;
+            M(2,1) = s;
+         }
+      }
+   }
 
-#define M(row,col)  m[col*4+row]
+   if (!optimized) {
+      const GLfloat mag = (GLfloat) GL_SQRT(x * x + y * y + z * z);
 
-   /*
-    *     Arbitrary axis rotation matrix.
-    *
-    *  This is composed of 5 matrices, Rz, Ry, T, Ry', Rz', multiplied
-    *  like so:  Rz * Ry * T * Ry' * Rz'.  T is the final rotation
-    *  (which is about the X-axis), and the two composite transforms
-    *  Ry' * Rz' and Rz * Ry are (respectively) the rotations necessary
-    *  from the arbitrary axis to the X-axis then back.  They are
-    *  all elementary rotations.
-    *
-    *  Rz' is a rotation about the Z-axis, to bring the axis vector
-    *  into the x-z plane.  Then Ry' is applied, rotating about the
-    *  Y-axis to bring the axis vector parallel with the X-axis.  The
-    *  rotation about the X-axis is then performed.  Ry and Rz are
-    *  simply the respective inverse transforms to bring the arbitrary
-    *  axis back to it's original orientation.  The first transforms
-    *  Rz' and Ry' are considered inverses, since the data from the
-    *  arbitrary axis gives you info on how to get to it, not how
-    *  to get away from it, and an inverse must be applied.
-    *
-    *  The basic calculation used is to recognize that the arbitrary
-    *  axis vector (x, y, z), since it is of unit length, actually
-    *  represents the sines and cosines of the angles to rotate the
-    *  X-axis to the same orientation, with theta being the angle about
-    *  Z and phi the angle about Y (in the order described above)
-    *  as follows:
-    *
-    *  cos ( theta ) = x / sqrt ( 1 - z^2 )
-    *  sin ( theta ) = y / sqrt ( 1 - z^2 )
-    *
-    *  cos ( phi ) = sqrt ( 1 - z^2 )
-    *  sin ( phi ) = z
-    *
-    *  Note that cos ( phi ) can further be inserted to the above
-    *  formulas:
-    *
-    *  cos ( theta ) = x / cos ( phi )
-    *  sin ( theta ) = y / sin ( phi )
-    *
-    *  ...etc.  Because of those relations and the standard trigonometric
-    *  relations, it is pssible to reduce the transforms down to what
-    *  is used below.  It may be that any primary axis chosen will give the
-    *  same results (modulo a sign convention) using thie method.
-    *
-    *  Particularly nice is to notice that all divisions that might
-    *  have caused trouble when parallel to certain planes or
-    *  axis go away with care paid to reducing the expressions.
-    *  After checking, it does perform correctly under all cases, since
-    *  in all the cases of division where the denominator would have
-    *  been zero, the numerator would have been zero as well, giving
-    *  the expected result.
-    */
+      if (mag <= 1.0e-4) {
+         /* no rotation, leave mat as-is */
+         return;
+      }
 
-   xx = x * x;
-   yy = y * y;
-   zz = z * z;
-   xy = x * y;
-   yz = y * z;
-   zx = z * x;
-   xs = x * s;
-   ys = y * s;
-   zs = z * s;
-   one_c = 1.0F - c;
-
-   M(0,0) = (one_c * xx) + c;
-   M(0,1) = (one_c * xy) - zs;
-   M(0,2) = (one_c * zx) + ys;
-   M(0,3) = 0.0F;
-
-   M(1,0) = (one_c * xy) + zs;
-   M(1,1) = (one_c * yy) + c;
-   M(1,2) = (one_c * yz) - xs;
-   M(1,3) = 0.0F;
-
-   M(2,0) = (one_c * zx) - ys;
-   M(2,1) = (one_c * yz) + xs;
-   M(2,2) = (one_c * zz) + c;
-   M(2,3) = 0.0F;
-
-   M(3,0) = 0.0F;
-   M(3,1) = 0.0F;
-   M(3,2) = 0.0F;
-   M(3,3) = 1.0F;
+      x /= mag;
+      y /= mag;
+      z /= mag;
+
+
+      /*
+       *     Arbitrary axis rotation matrix.
+       *
+       *  This is composed of 5 matrices, Rz, Ry, T, Ry', Rz', multiplied
+       *  like so:  Rz * Ry * T * Ry' * Rz'.  T is the final rotation
+       *  (which is about the X-axis), and the two composite transforms
+       *  Ry' * Rz' and Rz * Ry are (respectively) the rotations necessary
+       *  from the arbitrary axis to the X-axis then back.  They are
+       *  all elementary rotations.
+       *
+       *  Rz' is a rotation about the Z-axis, to bring the axis vector
+       *  into the x-z plane.  Then Ry' is applied, rotating about the
+       *  Y-axis to bring the axis vector parallel with the X-axis.  The
+       *  rotation about the X-axis is then performed.  Ry and Rz are
+       *  simply the respective inverse transforms to bring the arbitrary
+       *  axis back to it's original orientation.  The first transforms
+       *  Rz' and Ry' are considered inverses, since the data from the
+       *  arbitrary axis gives you info on how to get to it, not how
+       *  to get away from it, and an inverse must be applied.
+       *
+       *  The basic calculation used is to recognize that the arbitrary
+       *  axis vector (x, y, z), since it is of unit length, actually
+       *  represents the sines and cosines of the angles to rotate the
+       *  X-axis to the same orientation, with theta being the angle about
+       *  Z and phi the angle about Y (in the order described above)
+       *  as follows:
+       *
+       *  cos ( theta ) = x / sqrt ( 1 - z^2 )
+       *  sin ( theta ) = y / sqrt ( 1 - z^2 )
+       *
+       *  cos ( phi ) = sqrt ( 1 - z^2 )
+       *  sin ( phi ) = z
+       *
+       *  Note that cos ( phi ) can further be inserted to the above
+       *  formulas:
+       *
+       *  cos ( theta ) = x / cos ( phi )
+       *  sin ( theta ) = y / sin ( phi )
+       *
+       *  ...etc.  Because of those relations and the standard trigonometric
+       *  relations, it is pssible to reduce the transforms down to what
+       *  is used below.  It may be that any primary axis chosen will give the
+       *  same results (modulo a sign convention) using thie method.
+       *
+       *  Particularly nice is to notice that all divisions that might
+       *  have caused trouble when parallel to certain planes or
+       *  axis go away with care paid to reducing the expressions.
+       *  After checking, it does perform correctly under all cases, since
+       *  in all the cases of division where the denominator would have
+       *  been zero, the numerator would have been zero as well, giving
+       *  the expected result.
+       */
+
+      xx = x * x;
+      yy = y * y;
+      zz = z * z;
+      xy = x * y;
+      yz = y * z;
+      zx = z * x;
+      xs = x * s;
+      ys = y * s;
+      zs = z * s;
+      one_c = 1.0F - c;
+
+      /* We already hold the identity-matrix so we can skip some statements */
+      M(0,0) = (one_c * xx) + c;
+      M(0,1) = (one_c * xy) - zs;
+      M(0,2) = (one_c * zx) + ys;
+/*    M(0,3) = 0.0F; */
+
+      M(1,0) = (one_c * xy) + zs;
+      M(1,1) = (one_c * yy) + c;
+      M(1,2) = (one_c * yz) - xs;
+/*    M(1,3) = 0.0F; */
+
+      M(2,0) = (one_c * zx) - ys;
+      M(2,1) = (one_c * yz) + xs;
+      M(2,2) = (one_c * zz) + c;
+/*    M(2,3) = 0.0F; */
 
+/*
+      M(3,0) = 0.0F;
+      M(3,1) = 0.0F;
+      M(3,2) = 0.0F;
+      M(3,3) = 1.0F;
+*/
+   }
 #undef M
 
    matrix_multf( mat, m, MAT_FLAG_ROTATION );
 }
 
 
+
 void
 _math_matrix_frustum( GLmatrix *mat,
 		      GLfloat left, GLfloat right,