m_matrix.c

winNT技术操作系统,国外开放的原代码和LIUX一样

   if (0.0 == r1[1])  return GL_FALSE;

   /* eliminate second variable */
   m2 = r2[1]/r1[1]; m3 = r3[1]/r1[1];
   r2[2] -= m2 * r1[2]; r3[2] -= m3 * r1[2];
   r2[3] -= m2 * r1[3]; r3[3] -= m3 * r1[3];
   s = r1[4]; if (0.0 != s) { r2[4] -= m2 * s; r3[4] -= m3 * s; }
   s = r1[5]; if (0.0 != s) { r2[5] -= m2 * s; r3[5] -= m3 * s; }
   s = r1[6]; if (0.0 != s) { r2[6] -= m2 * s; r3[6] -= m3 * s; }
   s = r1[7]; if (0.0 != s) { r2[7] -= m2 * s; r3[7] -= m3 * s; }

   /* choose pivot - or die */
   if (fabs(r3[2])>fabs(r2[2])) SWAP_ROWS(r3, r2);
   if (0.0 == r2[2])  return GL_FALSE;

   /* eliminate third variable */
   m3 = r3[2]/r2[2];
   r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
   r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6],
   r3[7] -= m3 * r2[7];

   /* last check */
   if (0.0 == r3[3]) return GL_FALSE;

   s = 1.0F/r3[3];             /* now back substitute row 3 */
   r3[4] *= s; r3[5] *= s; r3[6] *= s; r3[7] *= s;

   m2 = r2[3];                 /* now back substitute row 2 */
   s  = 1.0F/r2[2];
   r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
   r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
   m1 = r1[3];
   r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
   r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
   m0 = r0[3];
   r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
   r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;

   m1 = r1[2];                 /* now back substitute row 1 */
   s  = 1.0F/r1[1];
   r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
   r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
   m0 = r0[2];
   r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
   r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;

   m0 = r0[1];                 /* now back substitute row 0 */
   s  = 1.0F/r0[0];
   r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
   r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);

   MAT(out,0,0) = r0[4]; MAT(out,0,1) = r0[5],
   MAT(out,0,2) = r0[6]; MAT(out,0,3) = r0[7],
   MAT(out,1,0) = r1[4]; MAT(out,1,1) = r1[5],
   MAT(out,1,2) = r1[6]; MAT(out,1,3) = r1[7],
   MAT(out,2,0) = r2[4]; MAT(out,2,1) = r2[5],
   MAT(out,2,2) = r2[6]; MAT(out,2,3) = r2[7],
   MAT(out,3,0) = r3[4]; MAT(out,3,1) = r3[5],
   MAT(out,3,2) = r3[6]; MAT(out,3,3) = r3[7];

   return GL_TRUE;
}
#undef SWAP_ROWS

/**
 * Compute inverse of a general 3d transformation matrix.
 * 
 * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 * stored in the GLmatrix::inv attribute.
 * 
 * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 *
 * \author Adapted from graphics gems II.
 *
 * Calculates the inverse of the upper left by first calculating its
 * determinant and multiplying it to the symmetric adjust matrix of each
 * element. Finally deals with the translation part by transforming the
 * original translation vector using by the calculated submatrix inverse.
 */
static GLboolean invert_matrix_3d_general( GLmatrix *mat )
{
   const GLfloat *in = mat->m;
   GLfloat *out = mat->inv;
   GLfloat pos, neg, t;
   GLfloat det;

   /* Calculate the determinant of upper left 3x3 submatrix and
    * determine if the matrix is singular.
    */
   pos = neg = 0.0;
   t =  MAT(in,0,0) * MAT(in,1,1) * MAT(in,2,2);
   if (t >= 0.0) pos += t; else neg += t;

   t =  MAT(in,1,0) * MAT(in,2,1) * MAT(in,0,2);
   if (t >= 0.0) pos += t; else neg += t;

   t =  MAT(in,2,0) * MAT(in,0,1) * MAT(in,1,2);
   if (t >= 0.0) pos += t; else neg += t;

   t = -MAT(in,2,0) * MAT(in,1,1) * MAT(in,0,2);
   if (t >= 0.0) pos += t; else neg += t;

   t = -MAT(in,1,0) * MAT(in,0,1) * MAT(in,2,2);
   if (t >= 0.0) pos += t; else neg += t;

   t = -MAT(in,0,0) * MAT(in,2,1) * MAT(in,1,2);
   if (t >= 0.0) pos += t; else neg += t;

   det = pos + neg;

   if (det*det < 1e-25)
      return GL_FALSE;

   det = 1.0F / det;
   MAT(out,0,0) = (  (MAT(in,1,1)*MAT(in,2,2) - MAT(in,2,1)*MAT(in,1,2) )*det);
   MAT(out,0,1) = (- (MAT(in,0,1)*MAT(in,2,2) - MAT(in,2,1)*MAT(in,0,2) )*det);
   MAT(out,0,2) = (  (MAT(in,0,1)*MAT(in,1,2) - MAT(in,1,1)*MAT(in,0,2) )*det);
   MAT(out,1,0) = (- (MAT(in,1,0)*MAT(in,2,2) - MAT(in,2,0)*MAT(in,1,2) )*det);
   MAT(out,1,1) = (  (MAT(in,0,0)*MAT(in,2,2) - MAT(in,2,0)*MAT(in,0,2) )*det);
   MAT(out,1,2) = (- (MAT(in,0,0)*MAT(in,1,2) - MAT(in,1,0)*MAT(in,0,2) )*det);
   MAT(out,2,0) = (  (MAT(in,1,0)*MAT(in,2,1) - MAT(in,2,0)*MAT(in,1,1) )*det);
   MAT(out,2,1) = (- (MAT(in,0,0)*MAT(in,2,1) - MAT(in,2,0)*MAT(in,0,1) )*det);
   MAT(out,2,2) = (  (MAT(in,0,0)*MAT(in,1,1) - MAT(in,1,0)*MAT(in,0,1) )*det);

   /* Do the translation part */
   MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0) +
		     MAT(in,1,3) * MAT(out,0,1) +
		     MAT(in,2,3) * MAT(out,0,2) );
   MAT(out,1,3) = - (MAT(in,0,3) * MAT(out,1,0) +
		     MAT(in,1,3) * MAT(out,1,1) +
		     MAT(in,2,3) * MAT(out,1,2) );
   MAT(out,2,3) = - (MAT(in,0,3) * MAT(out,2,0) +
		     MAT(in,1,3) * MAT(out,2,1) +
		     MAT(in,2,3) * MAT(out,2,2) );

   return GL_TRUE;
}

/**
 * Compute inverse of a 3d transformation matrix.
 * 
 * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 * stored in the GLmatrix::inv attribute.
 * 
 * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 *
 * If the matrix is not an angle preserving matrix then calls
 * invert_matrix_3d_general for the actual calculation. Otherwise calculates
 * the inverse matrix analyzing and inverting each of the scaling, rotation and
 * translation parts.
 */
static GLboolean invert_matrix_3d( GLmatrix *mat )
{
   const GLfloat *in = mat->m;
   GLfloat *out = mat->inv;

   if (!TEST_MAT_FLAGS(mat, MAT_FLAGS_ANGLE_PRESERVING)) {
      return invert_matrix_3d_general( mat );
   }

   if (mat->flags & MAT_FLAG_UNIFORM_SCALE) {
      GLfloat scale = (MAT(in,0,0) * MAT(in,0,0) +
                       MAT(in,0,1) * MAT(in,0,1) +
                       MAT(in,0,2) * MAT(in,0,2));

      if (scale == 0.0)
         return GL_FALSE;

      scale = 1.0F / scale;

      /* Transpose and scale the 3 by 3 upper-left submatrix. */
      MAT(out,0,0) = scale * MAT(in,0,0);
      MAT(out,1,0) = scale * MAT(in,0,1);
      MAT(out,2,0) = scale * MAT(in,0,2);
      MAT(out,0,1) = scale * MAT(in,1,0);
      MAT(out,1,1) = scale * MAT(in,1,1);
      MAT(out,2,1) = scale * MAT(in,1,2);
      MAT(out,0,2) = scale * MAT(in,2,0);
      MAT(out,1,2) = scale * MAT(in,2,1);
      MAT(out,2,2) = scale * MAT(in,2,2);
   }
   else if (mat->flags & MAT_FLAG_ROTATION) {
      /* Transpose the 3 by 3 upper-left submatrix. */
      MAT(out,0,0) = MAT(in,0,0);
      MAT(out,1,0) = MAT(in,0,1);
      MAT(out,2,0) = MAT(in,0,2);
      MAT(out,0,1) = MAT(in,1,0);
      MAT(out,1,1) = MAT(in,1,1);
      MAT(out,2,1) = MAT(in,1,2);
      MAT(out,0,2) = MAT(in,2,0);
      MAT(out,1,2) = MAT(in,2,1);
      MAT(out,2,2) = MAT(in,2,2);
   }
   else {
      /* pure translation */
      MEMCPY( out, Identity, sizeof(Identity) );
      MAT(out,0,3) = - MAT(in,0,3);
      MAT(out,1,3) = - MAT(in,1,3);
      MAT(out,2,3) = - MAT(in,2,3);
      return GL_TRUE;
   }

   if (mat->flags & MAT_FLAG_TRANSLATION) {
      /* Do the translation part */
      MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0) +
			MAT(in,1,3) * MAT(out,0,1) +
			MAT(in,2,3) * MAT(out,0,2) );
      MAT(out,1,3) = - (MAT(in,0,3) * MAT(out,1,0) +
			MAT(in,1,3) * MAT(out,1,1) +
			MAT(in,2,3) * MAT(out,1,2) );
      MAT(out,2,3) = - (MAT(in,0,3) * MAT(out,2,0) +
			MAT(in,1,3) * MAT(out,2,1) +
			MAT(in,2,3) * MAT(out,2,2) );
   }
   else {
      MAT(out,0,3) = MAT(out,1,3) = MAT(out,2,3) = 0.0;
   }

   return GL_TRUE;
}

/**
 * Compute inverse of an identity transformation matrix.
 * 
 * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 * stored in the GLmatrix::inv attribute.
 * 
 * \return always GL_TRUE.
 *
 * Simply copies Identity into GLmatrix::inv.
 */
static GLboolean invert_matrix_identity( GLmatrix *mat )
{
   MEMCPY( mat->inv, Identity, sizeof(Identity) );
   return GL_TRUE;
}

/**
 * Compute inverse of a no-rotation 3d transformation matrix.
 * 
 * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 * stored in the GLmatrix::inv attribute.
 * 
 * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 *
 * Calculates the 
 */
static GLboolean invert_matrix_3d_no_rot( GLmatrix *mat )
{
   const GLfloat *in = mat->m;
   GLfloat *out = mat->inv;

   if (MAT(in,0,0) == 0 || MAT(in,1,1) == 0 || MAT(in,2,2) == 0 )
      return GL_FALSE;

   MEMCPY( out, Identity, 16 * sizeof(GLfloat) );
   MAT(out,0,0) = 1.0F / MAT(in,0,0);
   MAT(out,1,1) = 1.0F / MAT(in,1,1);
   MAT(out,2,2) = 1.0F / MAT(in,2,2);

   if (mat->flags & MAT_FLAG_TRANSLATION) {
      MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0));
      MAT(out,1,3) = - (MAT(in,1,3) * MAT(out,1,1));
      MAT(out,2,3) = - (MAT(in,2,3) * MAT(out,2,2));
   }

   return GL_TRUE;
}

/**
 * Compute inverse of a no-rotation 2d transformation matrix.
 * 
 * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 * stored in the GLmatrix::inv attribute.
 * 
 * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 *
 * Calculates the inverse matrix by applying the inverse scaling and
 * translation to the identity matrix.
 */
static GLboolean invert_matrix_2d_no_rot( GLmatrix *mat )
{
   const GLfloat *in = mat->m;
   GLfloat *out = mat->inv;

   if (MAT(in,0,0) == 0 || MAT(in,1,1) == 0)
      return GL_FALSE;

   MEMCPY( out, Identity, 16 * sizeof(GLfloat) );
   MAT(out,0,0) = 1.0F / MAT(in,0,0);
   MAT(out,1,1) = 1.0F / MAT(in,1,1);

   if (mat->flags & MAT_FLAG_TRANSLATION) {
      MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0));
      MAT(out,1,3) = - (MAT(in,1,3) * MAT(out,1,1));
   }

   return GL_TRUE;
}

#if 0
/* broken */
static GLboolean invert_matrix_perspective( GLmatrix *mat )
{
   const GLfloat *in = mat->m;
   GLfloat *out = mat->inv;

   if (MAT(in,2,3) == 0)
      return GL_FALSE;

   MEMCPY( out, Identity, 16 * sizeof(GLfloat) );

   MAT(out,0,0) = 1.0F / MAT(in,0,0);
   MAT(out,1,1) = 1.0F / MAT(in,1,1);

   MAT(out,0,3) = MAT(in,0,2);
   MAT(out,1,3) = MAT(in,1,2);

   MAT(out,2,2) = 0;
   MAT(out,2,3) = -1;

   MAT(out,3,2) = 1.0F / MAT(in,2,3);
   MAT(out,3,3) = MAT(in,2,2) * MAT(out,3,2);

   return GL_TRUE;
}
#endif

/**
 * Matrix inversion function pointer type.
 */
typedef GLboolean (*inv_mat_func)( GLmatrix *mat );

/**
 * Table of the matrix inversion functions according to the matrix type.
 */
static inv_mat_func inv_mat_tab[7] = {
   invert_matrix_general,
   invert_matrix_identity,
   invert_matrix_3d_no_rot,
#if 0
   /* Don't use this function for now - it fails when the projection matrix
    * is premultiplied by a translation (ala Chromium's tilesort SPU).
    */
   invert_matrix_perspective,
#else
   invert_matrix_general,
#endif
   invert_matrix_3d,		/* lazy! */
   invert_matrix_2d_no_rot,
   invert_matrix_3d
};

/**
 * Compute inverse of a transformation matrix.
 * 
 * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 * stored in the GLmatrix::inv attribute.
 * 
 * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 *
 * Calls the matrix inversion function in inv_mat_tab corresponding to the
 * given matrix type.  In case of failure, updates the MAT_FLAG_SINGULAR flag,
 * and copies the identity matrix into GLmatrix::inv.
 */
static GLboolean matrix_invert( GLmatrix *mat )
{
   if (inv_mat_tab[mat->type](mat)) {
      mat->flags &= ~MAT_FLAG_SINGULAR;
      return GL_TRUE;
   } else {
      mat->flags |= MAT_FLAG_SINGULAR;
      MEMCPY( mat->inv, Identity, sizeof(Identity) );
      return GL_FALSE;
   }
}

/*@}*/


/**********************************************************************/
/** \name Matrix generation */
/*@{*/

/**
 * Generate a 4x4 transformation matrix from glRotate parameters, and
 * post-multiply the input matrix by it.
 *
 * \author
 * This function was contributed by Erich Boleyn ([e-mail]).
 * Optimizations contributed by Rudolf Opalla ([e-mail]).
 */
void
_math_matrix_rotate( GLmatrix *mat,
		     GLfloat angle, GLfloat x, GLfloat y, GLfloat z )
{
   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 );

   MEMCPY(m, Identity, sizeof(GLfloat)*16);
   optimized = GL_FALSE;