src/mesa/math/m_matrix.c

   1 /*
   2  * Mesa 3-D graphics library
   3  * Version:  6.3
   4  *
   5  * Copyright (C) 1999-2005  Brian Paul   All Rights Reserved.
   6  *
   7  * Permission is hereby granted, free of charge, to any person obtaining a
   8  * copy of this software and associated documentation files (the "Software"),
   9  * to deal in the Software without restriction, including without limitation
  10  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  11  * and/or sell copies of the Software, and to permit persons to whom the
  12  * Software is furnished to do so, subject to the following conditions:
  13  *
  14  * The above copyright notice and this permission notice shall be included
  15  * in all copies or substantial portions of the Software.
  16  *
  17  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  18  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  20  * BRIAN PAUL BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN
  21  * AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  22  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  23  */
  24
  25
  26 /**
  27  * \file m_matrix.c
  28  * Matrix operations.
  29  *
  30  * \note
  31  * -# 4x4 transformation matrices are stored in memory in column major order.
  32  * -# Points/vertices are to be thought of as column vectors.
  33  * -# Transformation of a point p by a matrix M is: p' = M * p
  34  */
  35
  36
  37 #include "main/glheader.h"
  38 #include "main/imports.h"
  39 #include "main/macros.h"
  40
  41 #include "m_matrix.h"
  42
  43
  44 /**
  45  * \defgroup MatFlags MAT_FLAG_XXX-flags
  46  *
  47  * Bitmasks to indicate different kinds of 4x4 matrices in GLmatrix::flags
  48  */
  49 /*@{*/
  50 #define MAT_FLAG_IDENTITY       0     /**< is an identity matrix flag.
  51                                        *   (Not actually used - the identity
  52                                        *   matrix is identified by the absense
  53                                        *   of all other flags.)
  54                                        */
  55 #define MAT_FLAG_GENERAL        0x1   /**< is a general matrix flag */
  56 #define MAT_FLAG_ROTATION       0x2   /**< is a rotation matrix flag */
  57 #define MAT_FLAG_TRANSLATION    0x4   /**< is a translation matrix flag */
  58 #define MAT_FLAG_UNIFORM_SCALE  0x8   /**< is an uniform scaling matrix flag */
  59 #define MAT_FLAG_GENERAL_SCALE  0x10  /**< is a general scaling matrix flag */
  60 #define MAT_FLAG_GENERAL_3D     0x20  /**< general 3D matrix flag */
  61 #define MAT_FLAG_PERSPECTIVE    0x40  /**< is a perspective proj matrix flag */
  62 #define MAT_FLAG_SINGULAR       0x80  /**< is a singular matrix flag */
  63 #define MAT_DIRTY_TYPE          0x100  /**< matrix type is dirty */
  64 #define MAT_DIRTY_FLAGS         0x200  /**< matrix flags are dirty */
  65 #define MAT_DIRTY_INVERSE       0x400  /**< matrix inverse is dirty */
  66
  67 /** angle preserving matrix flags mask */
  68 #define MAT_FLAGS_ANGLE_PRESERVING (MAT_FLAG_ROTATION | \
  69                                     MAT_FLAG_TRANSLATION | \
  70                                     MAT_FLAG_UNIFORM_SCALE)
  71
  72 /** geometry related matrix flags mask */
  73 #define MAT_FLAGS_GEOMETRY (MAT_FLAG_GENERAL | \
  74                             MAT_FLAG_ROTATION | \
  75                             MAT_FLAG_TRANSLATION | \
  76                             MAT_FLAG_UNIFORM_SCALE | \
  77                             MAT_FLAG_GENERAL_SCALE | \
  78                             MAT_FLAG_GENERAL_3D | \
  79                             MAT_FLAG_PERSPECTIVE | \
  80                             MAT_FLAG_SINGULAR)
  81
  82 /** length preserving matrix flags mask */
  83 #define MAT_FLAGS_LENGTH_PRESERVING (MAT_FLAG_ROTATION | \
  84                                      MAT_FLAG_TRANSLATION)
  85
  86
  87 /** 3D (non-perspective) matrix flags mask */
  88 #define MAT_FLAGS_3D (MAT_FLAG_ROTATION | \
  89                       MAT_FLAG_TRANSLATION | \
  90                       MAT_FLAG_UNIFORM_SCALE | \
  91                       MAT_FLAG_GENERAL_SCALE | \
  92                       MAT_FLAG_GENERAL_3D)
  93
  94 /** dirty matrix flags mask */
  95 #define MAT_DIRTY          (MAT_DIRTY_TYPE | \
  96                             MAT_DIRTY_FLAGS | \
  97                             MAT_DIRTY_INVERSE)
  98
  99 /*@}*/
 100
 101
 102 /**
 103  * Test geometry related matrix flags.
 104  *
 105  * \param mat a pointer to a GLmatrix structure.
 106  * \param a flags mask.
 107  *
 108  * \returns non-zero if all geometry related matrix flags are contained within
 109  * the mask, or zero otherwise.
 110  */
 111 #define TEST_MAT_FLAGS(mat, a)  \
 112     ((MAT_FLAGS_GEOMETRY & (~(a)) & ((mat)->flags) ) == 0)
 113
 114
 115
 116 /**
 117  * Names of the corresponding GLmatrixtype values.
 118  */
 119 static const char *types[] = {
 120    "MATRIX_GENERAL",
 121    "MATRIX_IDENTITY",
 122    "MATRIX_3D_NO_ROT",
 123    "MATRIX_PERSPECTIVE",
 124    "MATRIX_2D",
 125    "MATRIX_2D_NO_ROT",
 126    "MATRIX_3D"
 127 };
 128
 129
 130 /**
 131  * Identity matrix.
 132  */
 133 static GLfloat Identity[16] = {
 134    1.0, 0.0, 0.0, 0.0,
 135    0.0, 1.0, 0.0, 0.0,
 136    0.0, 0.0, 1.0, 0.0,
 137    0.0, 0.0, 0.0, 1.0
 138 };
 139
 140
 141
 142 /**********************************************************************/
 143 /** \name Matrix multiplication */
 144 /*@{*/
 145
 146 #define A(row,col)  a[(col<<2)+row]
 147 #define B(row,col)  b[(col<<2)+row]
 148 #define P(row,col)  product[(col<<2)+row]
 149
 150 /**
 151  * Perform a full 4x4 matrix multiplication.
 152  *
 153  * \param a matrix.
 154  * \param b matrix.
 155  * \param product will receive the product of \p a and \p b.
 156  *
 157  * \warning Is assumed that \p product != \p b. \p product == \p a is allowed.
 158  *
 159  * \note KW: 4*16 = 64 multiplications
 160  *
 161  * \author This \c matmul was contributed by Thomas Malik
 162  */
 163 static void matmul4( GLfloat *product, const GLfloat *a, const GLfloat *b )
 164 {
 165    GLint i;
 166    for (i = 0; i < 4; i++) {
 167       const GLfloat ai0=A(i,0),  ai1=A(i,1),  ai2=A(i,2),  ai3=A(i,3);
 168       P(i,0) = ai0 * B(0,0) + ai1 * B(1,0) + ai2 * B(2,0) + ai3 * B(3,0);
 169       P(i,1) = ai0 * B(0,1) + ai1 * B(1,1) + ai2 * B(2,1) + ai3 * B(3,1);
 170       P(i,2) = ai0 * B(0,2) + ai1 * B(1,2) + ai2 * B(2,2) + ai3 * B(3,2);
 171       P(i,3) = ai0 * B(0,3) + ai1 * B(1,3) + ai2 * B(2,3) + ai3 * B(3,3);
 172    }
 173 }
 174
 175 /**
 176  * Multiply two matrices known to occupy only the top three rows, such
 177  * as typical model matrices, and orthogonal matrices.
 178  *
 179  * \param a matrix.
 180  * \param b matrix.
 181  * \param product will receive the product of \p a and \p b.
 182  */
 183 static void matmul34( GLfloat *product, const GLfloat *a, const GLfloat *b )
 184 {
 185    GLint i;
 186    for (i = 0; i < 3; i++) {
 187       const GLfloat ai0=A(i,0),  ai1=A(i,1),  ai2=A(i,2),  ai3=A(i,3);
 188       P(i,0) = ai0 * B(0,0) + ai1 * B(1,0) + ai2 * B(2,0);
 189       P(i,1) = ai0 * B(0,1) + ai1 * B(1,1) + ai2 * B(2,1);
 190       P(i,2) = ai0 * B(0,2) + ai1 * B(1,2) + ai2 * B(2,2);
 191       P(i,3) = ai0 * B(0,3) + ai1 * B(1,3) + ai2 * B(2,3) + ai3;
 192    }
 193    P(3,0) = 0;
 194    P(3,1) = 0;
 195    P(3,2) = 0;
 196    P(3,3) = 1;
 197 }
 198
 199 #undef A
 200 #undef B
 201 #undef P
 202
 203 /**
 204  * Multiply a matrix by an array of floats with known properties.
 205  *
 206  * \param mat pointer to a GLmatrix structure containing the left multiplication
 207  * matrix, and that will receive the product result.
 208  * \param m right multiplication matrix array.
 209  * \param flags flags of the matrix \p m.
 210  *
 211  * Joins both flags and marks the type and inverse as dirty.  Calls matmul34()
 212  * if both matrices are 3D, or matmul4() otherwise.
 213  */
 214 static void matrix_multf( GLmatrix *mat, const GLfloat *m, GLuint flags )
 215 {
 216    mat->flags |= (flags | MAT_DIRTY_TYPE | MAT_DIRTY_INVERSE);
 217
 218    if (TEST_MAT_FLAGS(mat, MAT_FLAGS_3D))
 219       matmul34( mat->m, mat->m, m );
 220    else
 221       matmul4( mat->m, mat->m, m );
 222 }
 223
 224 /**
 225  * Matrix multiplication.
 226  *
 227  * \param dest destination matrix.
 228  * \param a left matrix.
 229  * \param b right matrix.
 230  *
 231  * Joins both flags and marks the type and inverse as dirty.  Calls matmul34()
 232  * if both matrices are 3D, or matmul4() otherwise.
 233  */
 234 void
 235 _math_matrix_mul_matrix( GLmatrix *dest, const GLmatrix *a, const GLmatrix *b )
 236 {
 237    dest->flags = (a->flags |
 238                   b->flags |
 239                   MAT_DIRTY_TYPE |
 240                   MAT_DIRTY_INVERSE);
 241
 242    if (TEST_MAT_FLAGS(dest, MAT_FLAGS_3D))
 243       matmul34( dest->m, a->m, b->m );
 244    else
 245       matmul4( dest->m, a->m, b->m );
 246 }
 247
 248 /**
 249  * Matrix multiplication.
 250  *
 251  * \param dest left and destination matrix.
 252  * \param m right matrix array.
 253  *
 254  * Marks the matrix flags with general flag, and type and inverse dirty flags.
 255  * Calls matmul4() for the multiplication.
 256  */
 257 void
 258 _math_matrix_mul_floats( GLmatrix *dest, const GLfloat *m )
 259 {
 260    dest->flags |= (MAT_FLAG_GENERAL |
 261                    MAT_DIRTY_TYPE |
 262                    MAT_DIRTY_INVERSE |
 263                    MAT_DIRTY_FLAGS);
 264
 265    matmul4( dest->m, dest->m, m );
 266 }
 267
 268 /*@}*/
 269
 270
 271 /**********************************************************************/
 272 /** \name Matrix output */
 273 /*@{*/
 274
 275 /**
 276  * Print a matrix array.
 277  *
 278  * \param m matrix array.
 279  *
 280  * Called by _math_matrix_print() to print a matrix or its inverse.
 281  */
 282 static void print_matrix_floats( const GLfloat m[16] )
 283 {
 284    int i;
 285    for (i=0;i<4;i++) {
 286       _mesa_debug(NULL,"\t%f %f %f %f\n", m[i], m[4+i], m[8+i], m[12+i] );
 287    }
 288 }
 289
 290 /**
 291  * Dumps the contents of a GLmatrix structure.
 292  *
 293  * \param m pointer to the GLmatrix structure.
 294  */
 295 void
 296 _math_matrix_print( const GLmatrix *m )
 297 {
 298    GLfloat prod[16];
 299
 300    _mesa_debug(NULL, "Matrix type: %s, flags: %x\n", types[m->type], m->flags);
 301    print_matrix_floats(m->m);
 302    _mesa_debug(NULL, "Inverse: \n");
 303    print_matrix_floats(m->inv);
 304    matmul4(prod, m->m, m->inv);
 305    _mesa_debug(NULL, "Mat * Inverse:\n");
 306    print_matrix_floats(prod);
 307 }
 308
 309 /*@}*/
 310
 311
 312 /**
 313  * References an element of 4x4 matrix.
 314  *
 315  * \param m matrix array.
 316  * \param c column of the desired element.
 317  * \param r row of the desired element.
 318  *
 319  * \return value of the desired element.
 320  *
 321  * Calculate the linear storage index of the element and references it.
 322  */
 323 #define MAT(m,r,c) (m)[(c)*4+(r)]
 324
 325
 326 /**********************************************************************/
 327 /** \name Matrix inversion */
 328 /*@{*/
 329
 330 /**
 331  * Swaps the values of two floating point variables.
 332  *
 333  * Used by invert_matrix_general() to swap the row pointers.
 334  */
 335 #define SWAP_ROWS(a, b) { GLfloat *_tmp = a; (a)=(b); (b)=_tmp; }
 336
 337 /**
 338  * Compute inverse of 4x4 transformation matrix.
 339  *
 340  * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 341  * stored in the GLmatrix::inv attribute.
 342  *
 343  * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 344  *
 345  * \author
 346  * Code contributed by Jacques Leroy jle@star.be
 347  *
 348  * Calculates the inverse matrix by performing the gaussian matrix reduction
 349  * with partial pivoting followed by back/substitution with the loops manually
 350  * unrolled.
 351  */
 352 static GLboolean invert_matrix_general( GLmatrix *mat )
 353 {
 354    const GLfloat *m = mat->m;
 355    GLfloat *out = mat->inv;
 356    GLfloat wtmp[4][8];
 357    GLfloat m0, m1, m2, m3, s;
 358    GLfloat *r0, *r1, *r2, *r3;
 359
 360    r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];
 361
 362    r0[0] = MAT(m,0,0), r0[1] = MAT(m,0,1),
 363    r0[2] = MAT(m,0,2), r0[3] = MAT(m,0,3),
 364    r0[4] = 1.0, r0[5] = r0[6] = r0[7] = 0.0,
 365
 366    r1[0] = MAT(m,1,0), r1[1] = MAT(m,1,1),
 367    r1[2] = MAT(m,1,2), r1[3] = MAT(m,1,3),
 368    r1[5] = 1.0, r1[4] = r1[6] = r1[7] = 0.0,
 369
 370    r2[0] = MAT(m,2,0), r2[1] = MAT(m,2,1),
 371    r2[2] = MAT(m,2,2), r2[3] = MAT(m,2,3),
 372    r2[6] = 1.0, r2[4] = r2[5] = r2[7] = 0.0,
 373
 374    r3[0] = MAT(m,3,0), r3[1] = MAT(m,3,1),
 375    r3[2] = MAT(m,3,2), r3[3] = MAT(m,3,3),
 376    r3[7] = 1.0, r3[4] = r3[5] = r3[6] = 0.0;
 377
 378    /* choose pivot - or die */
 379    if (FABSF(r3[0])>FABSF(r2[0])) SWAP_ROWS(r3, r2);
 380    if (FABSF(r2[0])>FABSF(r1[0])) SWAP_ROWS(r2, r1);
 381    if (FABSF(r1[0])>FABSF(r0[0])) SWAP_ROWS(r1, r0);
 382    if (0.0 == r0[0])  return GL_FALSE;
 383
 384    /* eliminate first variable     */
 385    m1 = r1[0]/r0[0]; m2 = r2[0]/r0[0]; m3 = r3[0]/r0[0];
 386    s = r0[1]; r1[1] -= m1 * s; r2[1] -= m2 * s; r3[1] -= m3 * s;
 387    s = r0[2]; r1[2] -= m1 * s; r2[2] -= m2 * s; r3[2] -= m3 * s;
 388    s = r0[3]; r1[3] -= m1 * s; r2[3] -= m2 * s; r3[3] -= m3 * s;
 389    s = r0[4];
 390    if (s != 0.0) { r1[4] -= m1 * s; r2[4] -= m2 * s; r3[4] -= m3 * s; }
 391    s = r0[5];
 392    if (s != 0.0) { r1[5] -= m1 * s; r2[5] -= m2 * s; r3[5] -= m3 * s; }
 393    s = r0[6];
 394    if (s != 0.0) { r1[6] -= m1 * s; r2[6] -= m2 * s; r3[6] -= m3 * s; }
 395    s = r0[7];
 396    if (s != 0.0) { r1[7] -= m1 * s; r2[7] -= m2 * s; r3[7] -= m3 * s; }
 397
 398    /* choose pivot - or die */
 399    if (FABSF(r3[1])>FABSF(r2[1])) SWAP_ROWS(r3, r2);
 400    if (FABSF(r2[1])>FABSF(r1[1])) SWAP_ROWS(r2, r1);
 401    if (0.0 == r1[1])  return GL_FALSE;
 402
 403    /* eliminate second variable */
 404    m2 = r2[1]/r1[1]; m3 = r3[1]/r1[1];
 405    r2[2] -= m2 * r1[2]; r3[2] -= m3 * r1[2];
 406    r2[3] -= m2 * r1[3]; r3[3] -= m3 * r1[3];
 407    s = r1[4]; if (0.0 != s) { r2[4] -= m2 * s; r3[4] -= m3 * s; }
 408    s = r1[5]; if (0.0 != s) { r2[5] -= m2 * s; r3[5] -= m3 * s; }
 409    s = r1[6]; if (0.0 != s) { r2[6] -= m2 * s; r3[6] -= m3 * s; }
 410    s = r1[7]; if (0.0 != s) { r2[7] -= m2 * s; r3[7] -= m3 * s; }
 411
 412    /* choose pivot - or die */
 413    if (FABSF(r3[2])>FABSF(r2[2])) SWAP_ROWS(r3, r2);
 414    if (0.0 == r2[2])  return GL_FALSE;
 415
 416    /* eliminate third variable */
 417    m3 = r3[2]/r2[2];
 418    r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
 419    r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6],
 420    r3[7] -= m3 * r2[7];
 421
 422    /* last check */
 423    if (0.0 == r3[3]) return GL_FALSE;
 424
 425    s = 1.0F/r3[3];             /* now back substitute row 3 */
 426    r3[4] *= s; r3[5] *= s; r3[6] *= s; r3[7] *= s;
 427
 428    m2 = r2[3];                 /* now back substitute row 2 */
 429    s  = 1.0F/r2[2];
 430    r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
 431    r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
 432    m1 = r1[3];
 433    r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
 434    r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
 435    m0 = r0[3];
 436    r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
 437    r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;
 438
 439    m1 = r1[2];                 /* now back substitute row 1 */
 440    s  = 1.0F/r1[1];
 441    r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
 442    r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
 443    m0 = r0[2];
 444    r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
 445    r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;
 446
 447    m0 = r0[1];                 /* now back substitute row 0 */
 448    s  = 1.0F/r0[0];
 449    r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
 450    r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);
 451
 452    MAT(out,0,0) = r0[4]; MAT(out,0,1) = r0[5],
 453    MAT(out,0,2) = r0[6]; MAT(out,0,3) = r0[7],
 454    MAT(out,1,0) = r1[4]; MAT(out,1,1) = r1[5],
 455    MAT(out,1,2) = r1[6]; MAT(out,1,3) = r1[7],
 456    MAT(out,2,0) = r2[4]; MAT(out,2,1) = r2[5],
 457    MAT(out,2,2) = r2[6]; MAT(out,2,3) = r2[7],
 458    MAT(out,3,0) = r3[4]; MAT(out,3,1) = r3[5],
 459    MAT(out,3,2) = r3[6]; MAT(out,3,3) = r3[7];
 460
 461    return GL_TRUE;
 462 }
 463 #undef SWAP_ROWS
 464
 465 /**
 466  * Compute inverse of a general 3d transformation matrix.
 467  *
 468  * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 469  * stored in the GLmatrix::inv attribute.
 470  *
 471  * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 472  *
 473  * \author Adapted from graphics gems II.
 474  *
 475  * Calculates the inverse of the upper left by first calculating its
 476  * determinant and multiplying it to the symmetric adjust matrix of each
 477  * element. Finally deals with the translation part by transforming the
 478  * original translation vector using by the calculated submatrix inverse.
 479  */
 480 static GLboolean invert_matrix_3d_general( GLmatrix *mat )
 481 {
 482    const GLfloat *in = mat->m;
 483    GLfloat *out = mat->inv;
 484    GLfloat pos, neg, t;
 485    GLfloat det;
 486
 487    /* Calculate the determinant of upper left 3x3 submatrix and
 488     * determine if the matrix is singular.
 489     */
 490    pos = neg = 0.0;
 491    t =  MAT(in,0,0) * MAT(in,1,1) * MAT(in,2,2);
 492    if (t >= 0.0) pos += t; else neg += t;
 493
 494    t =  MAT(in,1,0) * MAT(in,2,1) * MAT(in,0,2);
 495    if (t >= 0.0) pos += t; else neg += t;
 496
 497    t =  MAT(in,2,0) * MAT(in,0,1) * MAT(in,1,2);
 498    if (t >= 0.0) pos += t; else neg += t;
 499
 500    t = -MAT(in,2,0) * MAT(in,1,1) * MAT(in,0,2);
 501    if (t >= 0.0) pos += t; else neg += t;
 502
 503    t = -MAT(in,1,0) * MAT(in,0,1) * MAT(in,2,2);
 504    if (t >= 0.0) pos += t; else neg += t;
 505
 506    t = -MAT(in,0,0) * MAT(in,2,1) * MAT(in,1,2);
 507    if (t >= 0.0) pos += t; else neg += t;
 508
 509    det = pos + neg;
 510
 511    if (FABSF(det) < 1e-25)
 512       return GL_FALSE;
 513
 514    det = 1.0F / det;
 515    MAT(out,0,0) = (  (MAT(in,1,1)*MAT(in,2,2) - MAT(in,2,1)*MAT(in,1,2) )*det);
 516    MAT(out,0,1) = (- (MAT(in,0,1)*MAT(in,2,2) - MAT(in,2,1)*MAT(in,0,2) )*det);
 517    MAT(out,0,2) = (  (MAT(in,0,1)*MAT(in,1,2) - MAT(in,1,1)*MAT(in,0,2) )*det);
 518    MAT(out,1,0) = (- (MAT(in,1,0)*MAT(in,2,2) - MAT(in,2,0)*MAT(in,1,2) )*det);
 519    MAT(out,1,1) = (  (MAT(in,0,0)*MAT(in,2,2) - MAT(in,2,0)*MAT(in,0,2) )*det);
 520    MAT(out,1,2) = (- (MAT(in,0,0)*MAT(in,1,2) - MAT(in,1,0)*MAT(in,0,2) )*det);
 521    MAT(out,2,0) = (  (MAT(in,1,0)*MAT(in,2,1) - MAT(in,2,0)*MAT(in,1,1) )*det);
 522    MAT(out,2,1) = (- (MAT(in,0,0)*MAT(in,2,1) - MAT(in,2,0)*MAT(in,0,1) )*det);
 523    MAT(out,2,2) = (  (MAT(in,0,0)*MAT(in,1,1) - MAT(in,1,0)*MAT(in,0,1) )*det);
 524
 525    /* Do the translation part */
 526    MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0) +
 527                      MAT(in,1,3) * MAT(out,0,1) +
 528                      MAT(in,2,3) * MAT(out,0,2) );
 529    MAT(out,1,3) = - (MAT(in,0,3) * MAT(out,1,0) +
 530                      MAT(in,1,3) * MAT(out,1,1) +
 531                      MAT(in,2,3) * MAT(out,1,2) );
 532    MAT(out,2,3) = - (MAT(in,0,3) * MAT(out,2,0) +
 533                      MAT(in,1,3) * MAT(out,2,1) +
 534                      MAT(in,2,3) * MAT(out,2,2) );
 535
 536    return GL_TRUE;
 537 }
 538
 539 /**
 540  * Compute inverse of a 3d transformation matrix.
 541  *
 542  * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 543  * stored in the GLmatrix::inv attribute.
 544  *
 545  * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 546  *
 547  * If the matrix is not an angle preserving matrix then calls
 548  * invert_matrix_3d_general for the actual calculation. Otherwise calculates
 549  * the inverse matrix analyzing and inverting each of the scaling, rotation and
 550  * translation parts.
 551  */
 552 static GLboolean invert_matrix_3d( GLmatrix *mat )
 553 {
 554    const GLfloat *in = mat->m;
 555    GLfloat *out = mat->inv;
 556
 557    if (!TEST_MAT_FLAGS(mat, MAT_FLAGS_ANGLE_PRESERVING)) {
 558       return invert_matrix_3d_general( mat );
 559    }
 560
 561    if (mat->flags & MAT_FLAG_UNIFORM_SCALE) {
 562       GLfloat scale = (MAT(in,0,0) * MAT(in,0,0) +
 563                        MAT(in,0,1) * MAT(in,0,1) +
 564                        MAT(in,0,2) * MAT(in,0,2));
 565
 566       if (scale == 0.0)
 567          return GL_FALSE;
 568
 569       scale = 1.0F / scale;
 570
 571       /* Transpose and scale the 3 by 3 upper-left submatrix. */
 572       MAT(out,0,0) = scale * MAT(in,0,0);
 573       MAT(out,1,0) = scale * MAT(in,0,1);
 574       MAT(out,2,0) = scale * MAT(in,0,2);
 575       MAT(out,0,1) = scale * MAT(in,1,0);
 576       MAT(out,1,1) = scale * MAT(in,1,1);
 577       MAT(out,2,1) = scale * MAT(in,1,2);
 578       MAT(out,0,2) = scale * MAT(in,2,0);
 579       MAT(out,1,2) = scale * MAT(in,2,1);
 580       MAT(out,2,2) = scale * MAT(in,2,2);
 581    }
 582    else if (mat->flags & MAT_FLAG_ROTATION) {
 583       /* Transpose the 3 by 3 upper-left submatrix. */
 584       MAT(out,0,0) = MAT(in,0,0);
 585       MAT(out,1,0) = MAT(in,0,1);
 586       MAT(out,2,0) = MAT(in,0,2);
 587       MAT(out,0,1) = MAT(in,1,0);
 588       MAT(out,1,1) = MAT(in,1,1);
 589       MAT(out,2,1) = MAT(in,1,2);
 590       MAT(out,0,2) = MAT(in,2,0);
 591       MAT(out,1,2) = MAT(in,2,1);
 592       MAT(out,2,2) = MAT(in,2,2);
 593    }
 594    else {
 595       /* pure translation */
 596       memcpy( out, Identity, sizeof(Identity) );
 597       MAT(out,0,3) = - MAT(in,0,3);
 598       MAT(out,1,3) = - MAT(in,1,3);
 599       MAT(out,2,3) = - MAT(in,2,3);
 600       return GL_TRUE;
 601    }
 602
 603    if (mat->flags & MAT_FLAG_TRANSLATION) {
 604       /* Do the translation part */
 605       MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0) +
 606                         MAT(in,1,3) * MAT(out,0,1) +
 607                         MAT(in,2,3) * MAT(out,0,2) );
 608       MAT(out,1,3) = - (MAT(in,0,3) * MAT(out,1,0) +
 609                         MAT(in,1,3) * MAT(out,1,1) +
 610                         MAT(in,2,3) * MAT(out,1,2) );
 611       MAT(out,2,3) = - (MAT(in,0,3) * MAT(out,2,0) +
 612                         MAT(in,1,3) * MAT(out,2,1) +
 613                         MAT(in,2,3) * MAT(out,2,2) );
 614    }
 615    else {
 616       MAT(out,0,3) = MAT(out,1,3) = MAT(out,2,3) = 0.0;
 617    }
 618
 619    return GL_TRUE;
 620 }
 621
 622 /**
 623  * Compute inverse of an identity transformation matrix.
 624  *
 625  * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 626  * stored in the GLmatrix::inv attribute.
 627  *
 628  * \return always GL_TRUE.
 629  *
 630  * Simply copies Identity into GLmatrix::inv.
 631  */
 632 static GLboolean invert_matrix_identity( GLmatrix *mat )
 633 {
 634    memcpy( mat->inv, Identity, sizeof(Identity) );
 635    return GL_TRUE;
 636 }
 637
 638 /**
 639  * Compute inverse of a no-rotation 3d transformation matrix.
 640  *
 641  * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 642  * stored in the GLmatrix::inv attribute.
 643  *
 644  * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 645  *
 646  * Calculates the
 647  */
 648 static GLboolean invert_matrix_3d_no_rot( GLmatrix *mat )
 649 {
 650    const GLfloat *in = mat->m;
 651    GLfloat *out = mat->inv;
 652
 653    if (MAT(in,0,0) == 0 || MAT(in,1,1) == 0 || MAT(in,2,2) == 0 )
 654       return GL_FALSE;
 655
 656    memcpy( out, Identity, 16 * sizeof(GLfloat) );
 657    MAT(out,0,0) = 1.0F / MAT(in,0,0);
 658    MAT(out,1,1) = 1.0F / MAT(in,1,1);
 659    MAT(out,2,2) = 1.0F / MAT(in,2,2);
 660
 661    if (mat->flags & MAT_FLAG_TRANSLATION) {
 662       MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0));
 663       MAT(out,1,3) = - (MAT(in,1,3) * MAT(out,1,1));
 664       MAT(out,2,3) = - (MAT(in,2,3) * MAT(out,2,2));
 665    }
 666
 667    return GL_TRUE;
 668 }
 669
 670 /**
 671  * Compute inverse of a no-rotation 2d transformation matrix.
 672  *
 673  * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 674  * stored in the GLmatrix::inv attribute.
 675  *
 676  * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 677  *
 678  * Calculates the inverse matrix by applying the inverse scaling and
 679  * translation to the identity matrix.
 680  */
 681 static GLboolean invert_matrix_2d_no_rot( GLmatrix *mat )
 682 {
 683    const GLfloat *in = mat->m;
 684    GLfloat *out = mat->inv;
 685
 686    if (MAT(in,0,0) == 0 || MAT(in,1,1) == 0)
 687       return GL_FALSE;
 688
 689    memcpy( out, Identity, 16 * sizeof(GLfloat) );
 690    MAT(out,0,0) = 1.0F / MAT(in,0,0);
 691    MAT(out,1,1) = 1.0F / MAT(in,1,1);
 692
 693    if (mat->flags & MAT_FLAG_TRANSLATION) {
 694       MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0));
 695       MAT(out,1,3) = - (MAT(in,1,3) * MAT(out,1,1));
 696    }
 697
 698    return GL_TRUE;
 699 }
 700
 701 #if 0
 702 /* broken */
 703 static GLboolean invert_matrix_perspective( GLmatrix *mat )
 704 {
 705    const GLfloat *in = mat->m;
 706    GLfloat *out = mat->inv;
 707
 708    if (MAT(in,2,3) == 0)
 709       return GL_FALSE;
 710
 711    memcpy( out, Identity, 16 * sizeof(GLfloat) );
 712
 713    MAT(out,0,0) = 1.0F / MAT(in,0,0);
 714    MAT(out,1,1) = 1.0F / MAT(in,1,1);
 715
 716    MAT(out,0,3) = MAT(in,0,2);
 717    MAT(out,1,3) = MAT(in,1,2);
 718
 719    MAT(out,2,2) = 0;
 720    MAT(out,2,3) = -1;
 721
 722    MAT(out,3,2) = 1.0F / MAT(in,2,3);
 723    MAT(out,3,3) = MAT(in,2,2) * MAT(out,3,2);
 724
 725    return GL_TRUE;
 726 }
 727 #endif
 728
 729 /**
 730  * Matrix inversion function pointer type.
 731  */
 732 typedef GLboolean (*inv_mat_func)( GLmatrix *mat );
 733
 734 /**
 735  * Table of the matrix inversion functions according to the matrix type.
 736  */
 737 static inv_mat_func inv_mat_tab[7] = {
 738    invert_matrix_general,
 739    invert_matrix_identity,
 740    invert_matrix_3d_no_rot,
 741 #if 0
 742    /* Don't use this function for now - it fails when the projection matrix
 743     * is premultiplied by a translation (ala Chromium's tilesort SPU).
 744     */
 745    invert_matrix_perspective,
 746 #else
 747    invert_matrix_general,
 748 #endif
 749    invert_matrix_3d,            /* lazy! */
 750    invert_matrix_2d_no_rot,
 751    invert_matrix_3d
 752 };
 753
 754 /**
 755  * Compute inverse of a transformation matrix.
 756  *
 757  * \param mat pointer to a GLmatrix structure. The matrix inverse will be
 758  * stored in the GLmatrix::inv attribute.
 759  *
 760  * \return GL_TRUE for success, GL_FALSE for failure (\p singular matrix).
 761  *
 762  * Calls the matrix inversion function in inv_mat_tab corresponding to the
 763  * given matrix type.  In case of failure, updates the MAT_FLAG_SINGULAR flag,
 764  * and copies the identity matrix into GLmatrix::inv.
 765  */
 766 static GLboolean matrix_invert( GLmatrix *mat )
 767 {
 768    if (inv_mat_tab[mat->type](mat)) {
 769       mat->flags &= ~MAT_FLAG_SINGULAR;
 770       return GL_TRUE;
 771    } else {
 772       mat->flags |= MAT_FLAG_SINGULAR;
 773       memcpy( mat->inv, Identity, sizeof(Identity) );
 774       return GL_FALSE;
 775    }
 776 }
 777
 778 /*@}*/
 779
 780
 781 /**********************************************************************/
 782 /** \name Matrix generation */
 783 /*@{*/
 784
 785 /**
 786  * Generate a 4x4 transformation matrix from glRotate parameters, and
 787  * post-multiply the input matrix by it.
 788  *
 789  * \author
 790  * This function was contributed by Erich Boleyn (erich@uruk.org).
 791  * Optimizations contributed by Rudolf Opalla (rudi@khm.de).
 792  */
 793 void
 794 _math_matrix_rotate( GLmatrix *mat,
 795                      GLfloat angle, GLfloat x, GLfloat y, GLfloat z )
 796 {
 797    GLfloat xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c, s, c;
 798    GLfloat m[16];
 799    GLboolean optimized;
 800
 801    s = (GLfloat) sin( angle * DEG2RAD );
 802    c = (GLfloat) cos( angle * DEG2RAD );
 803
 804    memcpy(m, Identity, sizeof(GLfloat)*16);
 805    optimized = GL_FALSE;
 806
 807 #define M(row,col)  m[col*4+row]
 808
 809    if (x == 0.0F) {
 810       if (y == 0.0F) {
 811          if (z != 0.0F) {
 812             optimized = GL_TRUE;
 813             /* rotate only around z-axis */
 814             M(0,0) = c;
 815             M(1,1) = c;
 816             if (z < 0.0F) {
 817                M(0,1) = s;
 818                M(1,0) = -s;
 819             }
 820             else {
 821                M(0,1) = -s;
 822                M(1,0) = s;
 823             }
 824          }
 825       }
 826       else if (z == 0.0F) {
 827          optimized = GL_TRUE;
 828          /* rotate only around y-axis */
 829          M(0,0) = c;
 830          M(2,2) = c;
 831          if (y < 0.0F) {
 832             M(0,2) = -s;
 833             M(2,0) = s;
 834          }
 835          else {
 836             M(0,2) = s;
 837             M(2,0) = -s;
 838          }
 839       }
 840    }
 841    else if (y == 0.0F) {
 842       if (z == 0.0F) {
 843          optimized = GL_TRUE;
 844          /* rotate only around x-axis */
 845          M(1,1) = c;
 846          M(2,2) = c;
 847          if (x < 0.0F) {
 848             M(1,2) = s;
 849             M(2,1) = -s;
 850          }
 851          else {
 852             M(1,2) = -s;
 853             M(2,1) = s;
 854          }
 855       }
 856    }
 857
 858    if (!optimized) {
 859       const GLfloat mag = SQRTF(x * x + y * y + z * z);
 860
 861       if (mag <= 1.0e-4) {
 862          /* no rotation, leave mat as-is */
 863          return;
 864       }
 865
 866       x /= mag;
 867       y /= mag;
 868       z /= mag;
 869
 870
 871       /*
 872        *     Arbitrary axis rotation matrix.
 873        *
 874        *  This is composed of 5 matrices, Rz, Ry, T, Ry', Rz', multiplied
 875        *  like so:  Rz * Ry * T * Ry' * Rz'.  T is the final rotation
 876        *  (which is about the X-axis), and the two composite transforms
 877        *  Ry' * Rz' and Rz * Ry are (respectively) the rotations necessary
 878        *  from the arbitrary axis to the X-axis then back.  They are
 879        *  all elementary rotations.
 880        *
 881        *  Rz' is a rotation about the Z-axis, to bring the axis vector
 882        *  into the x-z plane.  Then Ry' is applied, rotating about the
 883        *  Y-axis to bring the axis vector parallel with the X-axis.  The
 884        *  rotation about the X-axis is then performed.  Ry and Rz are
 885        *  simply the respective inverse transforms to bring the arbitrary
 886        *  axis back to its original orientation.  The first transforms
 887        *  Rz' and Ry' are considered inverses, since the data from the
 888        *  arbitrary axis gives you info on how to get to it, not how
 889        *  to get away from it, and an inverse must be applied.
 890        *
 891        *  The basic calculation used is to recognize that the arbitrary
 892        *  axis vector (x, y, z), since it is of unit length, actually
 893        *  represents the sines and cosines of the angles to rotate the
 894        *  X-axis to the same orientation, with theta being the angle about
 895        *  Z and phi the angle about Y (in the order described above)
 896        *  as follows:
 897        *
 898        *  cos ( theta ) = x / sqrt ( 1 - z^2 )
 899        *  sin ( theta ) = y / sqrt ( 1 - z^2 )
 900        *
 901        *  cos ( phi ) = sqrt ( 1 - z^2 )
 902        *  sin ( phi ) = z
 903        *
 904        *  Note that cos ( phi ) can further be inserted to the above
 905        *  formulas:
 906        *
 907        *  cos ( theta ) = x / cos ( phi )
 908        *  sin ( theta ) = y / sin ( phi )
 909        *
 910        *  ...etc.  Because of those relations and the standard trigonometric
 911        *  relations, it is pssible to reduce the transforms down to what
 912        *  is used below.  It may be that any primary axis chosen will give the
 913        *  same results (modulo a sign convention) using thie method.
 914        *
 915        *  Particularly nice is to notice that all divisions that might
 916        *  have caused trouble when parallel to certain planes or
 917        *  axis go away with care paid to reducing the expressions.
 918        *  After checking, it does perform correctly under all cases, since
 919        *  in all the cases of division where the denominator would have
 920        *  been zero, the numerator would have been zero as well, giving
 921        *  the expected result.
 922        */
 923
 924       xx = x * x;
 925       yy = y * y;
 926       zz = z * z;
 927       xy = x * y;
 928       yz = y * z;
 929       zx = z * x;
 930       xs = x * s;
 931       ys = y * s;
 932       zs = z * s;
 933       one_c = 1.0F - c;
 934
 935       /* We already hold the identity-matrix so we can skip some statements */
 936       M(0,0) = (one_c * xx) + c;
 937       M(0,1) = (one_c * xy) - zs;
 938       M(0,2) = (one_c * zx) + ys;
 939 /*    M(0,3) = 0.0F; */
 940
 941       M(1,0) = (one_c * xy) + zs;
 942       M(1,1) = (one_c * yy) + c;
 943       M(1,2) = (one_c * yz) - xs;
 944 /*    M(1,3) = 0.0F; */
 945
 946       M(2,0) = (one_c * zx) - ys;
 947       M(2,1) = (one_c * yz) + xs;
 948       M(2,2) = (one_c * zz) + c;
 949 /*    M(2,3) = 0.0F; */
 950
 951 /*
 952       M(3,0) = 0.0F;
 953       M(3,1) = 0.0F;
 954       M(3,2) = 0.0F;
 955       M(3,3) = 1.0F;
 956 */
 957    }
 958 #undef M
 959
 960    matrix_multf( mat, m, MAT_FLAG_ROTATION );
 961 }
 962
 963 /**
 964  * Apply a perspective projection matrix.
 965  *
 966  * \param mat matrix to apply the projection.
 967  * \param left left clipping plane coordinate.
 968  * \param right right clipping plane coordinate.
 969  * \param bottom bottom clipping plane coordinate.
 970  * \param top top clipping plane coordinate.
 971  * \param nearval distance to the near clipping plane.
 972  * \param farval distance to the far clipping plane.
 973  *
 974  * Creates the projection matrix and multiplies it with \p mat, marking the
 975  * MAT_FLAG_PERSPECTIVE flag.
 976  */
 977 void
 978 _math_matrix_frustum( GLmatrix *mat,
 979                       GLfloat left, GLfloat right,
 980                       GLfloat bottom, GLfloat top,
 981                       GLfloat nearval, GLfloat farval )
 982 {
 983    GLfloat x, y, a, b, c, d;
 984    GLfloat m[16];
 985
 986    x = (2.0F*nearval) / (right-left);
 987    y = (2.0F*nearval) / (top-bottom);
 988    a = (right+left) / (right-left);
 989    b = (top+bottom) / (top-bottom);
 990    c = -(farval+nearval) / ( farval-nearval);
 991    d = -(2.0F*farval*nearval) / (farval-nearval);  /* error? */
 992
 993 #define M(row,col)  m[col*4+row]
 994    M(0,0) = x;     M(0,1) = 0.0F;  M(0,2) = a;      M(0,3) = 0.0F;
 995    M(1,0) = 0.0F;  M(1,1) = y;     M(1,2) = b;      M(1,3) = 0.0F;
 996    M(2,0) = 0.0F;  M(2,1) = 0.0F;  M(2,2) = c;      M(2,3) = d;
 997    M(3,0) = 0.0F;  M(3,1) = 0.0F;  M(3,2) = -1.0F;  M(3,3) = 0.0F;
 998 #undef M
 999
1000    matrix_multf( mat, m, MAT_FLAG_PERSPECTIVE );
1001 }
1002
1003 /**
1004  * Apply an orthographic projection matrix.
1005  *
1006  * \param mat matrix to apply the projection.
1007  * \param left left clipping plane coordinate.
1008  * \param right right clipping plane coordinate.
1009  * \param bottom bottom clipping plane coordinate.
1010  * \param top top clipping plane coordinate.
1011  * \param nearval distance to the near clipping plane.
1012  * \param farval distance to the far clipping plane.
1013  *
1014  * Creates the projection matrix and multiplies it with \p mat, marking the
1015  * MAT_FLAG_GENERAL_SCALE and MAT_FLAG_TRANSLATION flags.
1016  */
1017 void
1018 _math_matrix_ortho( GLmatrix *mat,
1019                     GLfloat left, GLfloat right,
1020                     GLfloat bottom, GLfloat top,
1021                     GLfloat nearval, GLfloat farval )
1022 {
1023    GLfloat m[16];
1024
1025 #define M(row,col)  m[col*4+row]
1026    M(0,0) = 2.0F / (right-left);
1027    M(0,1) = 0.0F;
1028    M(0,2) = 0.0F;
1029    M(0,3) = -(right+left) / (right-left);
1030
1031    M(1,0) = 0.0F;
1032    M(1,1) = 2.0F / (top-bottom);
1033    M(1,2) = 0.0F;
1034    M(1,3) = -(top+bottom) / (top-bottom);
1035
1036    M(2,0) = 0.0F;
1037    M(2,1) = 0.0F;
1038    M(2,2) = -2.0F / (farval-nearval);
1039    M(2,3) = -(farval+nearval) / (farval-nearval);
1040
1041    M(3,0) = 0.0F;
1042    M(3,1) = 0.0F;
1043    M(3,2) = 0.0F;
1044    M(3,3) = 1.0F;
1045 #undef M
1046
1047    matrix_multf( mat, m, (MAT_FLAG_GENERAL_SCALE|MAT_FLAG_TRANSLATION));
1048 }
1049
1050 /**
1051  * Multiply a matrix with a general scaling matrix.
1052  *
1053  * \param mat matrix.
1054  * \param x x axis scale factor.
1055  * \param y y axis scale factor.
1056  * \param z z axis scale factor.
1057  *
1058  * Multiplies in-place the elements of \p mat by the scale factors. Checks if
1059  * the scales factors are roughly the same, marking the MAT_FLAG_UNIFORM_SCALE
1060  * flag, or MAT_FLAG_GENERAL_SCALE. Marks the MAT_DIRTY_TYPE and
1061  * MAT_DIRTY_INVERSE dirty flags.
1062  */
1063 void
1064 _math_matrix_scale( GLmatrix *mat, GLfloat x, GLfloat y, GLfloat z )
1065 {
1066    GLfloat *m = mat->m;
1067    m[0] *= x;   m[4] *= y;   m[8]  *= z;
1068    m[1] *= x;   m[5] *= y;   m[9]  *= z;
1069    m[2] *= x;   m[6] *= y;   m[10] *= z;
1070    m[3] *= x;   m[7] *= y;   m[11] *= z;
1071
1072    if (FABSF(x - y) < 1e-8 && FABSF(x - z) < 1e-8)
1073       mat->flags |= MAT_FLAG_UNIFORM_SCALE;
1074    else
1075       mat->flags |= MAT_FLAG_GENERAL_SCALE;
1076
1077    mat->flags |= (MAT_DIRTY_TYPE |
1078                   MAT_DIRTY_INVERSE);
1079 }
1080
1081 /**
1082  * Multiply a matrix with a translation matrix.
1083  *
1084  * \param mat matrix.
1085  * \param x translation vector x coordinate.
1086  * \param y translation vector y coordinate.
1087  * \param z translation vector z coordinate.
1088  *
1089  * Adds the translation coordinates to the elements of \p mat in-place.  Marks
1090  * the MAT_FLAG_TRANSLATION flag, and the MAT_DIRTY_TYPE and MAT_DIRTY_INVERSE
1091  * dirty flags.
1092  */
1093 void
1094 _math_matrix_translate( GLmatrix *mat, GLfloat x, GLfloat y, GLfloat z )
1095 {
1096    GLfloat *m = mat->m;
1097    m[12] = m[0] * x + m[4] * y + m[8]  * z + m[12];
1098    m[13] = m[1] * x + m[5] * y + m[9]  * z + m[13];
1099    m[14] = m[2] * x + m[6] * y + m[10] * z + m[14];
1100    m[15] = m[3] * x + m[7] * y + m[11] * z + m[15];
1101
1102    mat->flags |= (MAT_FLAG_TRANSLATION |
1103                   MAT_DIRTY_TYPE |
1104                   MAT_DIRTY_INVERSE);
1105 }
1106
1107
1108 /**
1109  * Set matrix to do viewport and depthrange mapping.
1110  * Transforms Normalized Device Coords to window/Z values.
1111  */
1112 void
1113 _math_matrix_viewport(GLmatrix *m, GLint x, GLint y, GLint width, GLint height,
1114                       GLfloat zNear, GLfloat zFar, GLfloat depthMax)
1115 {
1116    m->m[MAT_SX] = (GLfloat) width / 2.0F;
1117    m->m[MAT_TX] = m->m[MAT_SX] + x;
1118    m->m[MAT_SY] = (GLfloat) height / 2.0F;
1119    m->m[MAT_TY] = m->m[MAT_SY] + y;
1120    m->m[MAT_SZ] = depthMax * ((zFar - zNear) / 2.0F);
1121    m->m[MAT_TZ] = depthMax * ((zFar - zNear) / 2.0F + zNear);
1122    m->flags = MAT_FLAG_GENERAL_SCALE | MAT_FLAG_TRANSLATION;
1123    m->type = MATRIX_3D_NO_ROT;
1124 }
1125
1126
1127 /**
1128  * Set a matrix to the identity matrix.
1129  *
1130  * \param mat matrix.
1131  *
1132  * Copies ::Identity into \p GLmatrix::m, and into GLmatrix::inv if not NULL.
1133  * Sets the matrix type to identity, and clear the dirty flags.
1134  */
1135 void
1136 _math_matrix_set_identity( GLmatrix *mat )
1137 {
1138    memcpy( mat->m, Identity, 16*sizeof(GLfloat) );
1139    memcpy( mat->inv, Identity, 16*sizeof(GLfloat) );
1140
1141    mat->type = MATRIX_IDENTITY;
1142    mat->flags &= ~(MAT_DIRTY_FLAGS|
1143                    MAT_DIRTY_TYPE|
1144                    MAT_DIRTY_INVERSE);
1145 }
1146
1147 /*@}*/
1148
1149
1150 /**********************************************************************/
1151 /** \name Matrix analysis */
1152 /*@{*/
1153
1154 #define ZERO(x) (1<<x)
1155 #define ONE(x)  (1<<(x+16))
1156
1157 #define MASK_NO_TRX      (ZERO(12) | ZERO(13) | ZERO(14))
1158 #define MASK_NO_2D_SCALE ( ONE(0)  | ONE(5))
1159
1160 #define MASK_IDENTITY    ( ONE(0)  | ZERO(4)  | ZERO(8)  | ZERO(12) |\
1161                           ZERO(1)  |  ONE(5)  | ZERO(9)  | ZERO(13) |\
1162                           ZERO(2)  | ZERO(6)  |  ONE(10) | ZERO(14) |\
1163                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
1164
1165 #define MASK_2D_NO_ROT   (           ZERO(4)  | ZERO(8)  |           \
1166                           ZERO(1)  |            ZERO(9)  |           \
1167                           ZERO(2)  | ZERO(6)  |  ONE(10) | ZERO(14) |\
1168                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
1169
1170 #define MASK_2D          (                      ZERO(8)  |           \
1171                                                 ZERO(9)  |           \
1172                           ZERO(2)  | ZERO(6)  |  ONE(10) | ZERO(14) |\
1173                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
1174
1175
1176 #define MASK_3D_NO_ROT   (           ZERO(4)  | ZERO(8)  |           \
1177                           ZERO(1)  |            ZERO(9)  |           \
1178                           ZERO(2)  | ZERO(6)  |                      \
1179                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
1180
1181 #define MASK_3D          (                                           \
1182                                                                      \
1183                                                                      \
1184                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
1185
1186
1187 #define MASK_PERSPECTIVE (           ZERO(4)  |            ZERO(12) |\
1188                           ZERO(1)  |                       ZERO(13) |\
1189                           ZERO(2)  | ZERO(6)  |                      \
1190                           ZERO(3)  | ZERO(7)  |            ZERO(15) )
1191
1192 #define SQ(x) ((x)*(x))
1193
1194 /**
1195  * Determine type and flags from scratch.
1196  *
1197  * \param mat matrix.
1198  *
1199  * This is expensive enough to only want to do it once.
1200  */
1201 static void analyse_from_scratch( GLmatrix *mat )
1202 {
1203    const GLfloat *m = mat->m;
1204    GLuint mask = 0;
1205    GLuint i;
1206
1207    for (i = 0 ; i < 16 ; i++) {
1208       if (m[i] == 0.0) mask |= (1<<i);
1209    }
1210
1211    if (m[0] == 1.0F) mask |= (1<<16);
1212    if (m[5] == 1.0F) mask |= (1<<21);
1213    if (m[10] == 1.0F) mask |= (1<<26);
1214    if (m[15] == 1.0F) mask |= (1<<31);
1215
1216    mat->flags &= ~MAT_FLAGS_GEOMETRY;
1217
1218    /* Check for translation - no-one really cares
1219     */
1220    if ((mask & MASK_NO_TRX) != MASK_NO_TRX)
1221       mat->flags |= MAT_FLAG_TRANSLATION;
1222
1223    /* Do the real work
1224     */
1225    if (mask == (GLuint) MASK_IDENTITY) {
1226       mat->type = MATRIX_IDENTITY;
1227    }
1228    else if ((mask & MASK_2D_NO_ROT) == (GLuint) MASK_2D_NO_ROT) {
1229       mat->type = MATRIX_2D_NO_ROT;
1230
1231       if ((mask & MASK_NO_2D_SCALE) != MASK_NO_2D_SCALE)
1232          mat->flags |= MAT_FLAG_GENERAL_SCALE;
1233    }
1234    else if ((mask & MASK_2D) == (GLuint) MASK_2D) {
1235       GLfloat mm = DOT2(m, m);
1236       GLfloat m4m4 = DOT2(m+4,m+4);
1237       GLfloat mm4 = DOT2(m,m+4);
1238
1239       mat->type = MATRIX_2D;
1240
1241       /* Check for scale */
1242       if (SQ(mm-1) > SQ(1e-6) ||
1243           SQ(m4m4-1) > SQ(1e-6))
1244          mat->flags |= MAT_FLAG_GENERAL_SCALE;
1245
1246       /* Check for rotation */
1247       if (SQ(mm4) > SQ(1e-6))
1248          mat->flags |= MAT_FLAG_GENERAL_3D;
1249       else
1250          mat->flags |= MAT_FLAG_ROTATION;
1251
1252    }
1253    else if ((mask & MASK_3D_NO_ROT) == (GLuint) MASK_3D_NO_ROT) {
1254       mat->type = MATRIX_3D_NO_ROT;
1255
1256       /* Check for scale */
1257       if (SQ(m[0]-m[5]) < SQ(1e-6) &&
1258           SQ(m[0]-m[10]) < SQ(1e-6)) {
1259          if (SQ(m[0]-1.0) > SQ(1e-6)) {
1260             mat->flags |= MAT_FLAG_UNIFORM_SCALE;
1261          }
1262       }
1263       else {
1264          mat->flags |= MAT_FLAG_GENERAL_SCALE;
1265       }
1266    }
1267    else if ((mask & MASK_3D) == (GLuint) MASK_3D) {
1268       GLfloat c1 = DOT3(m,m);
1269       GLfloat c2 = DOT3(m+4,m+4);
1270       GLfloat c3 = DOT3(m+8,m+8);
1271       GLfloat d1 = DOT3(m, m+4);
1272       GLfloat cp[3];
1273
1274       mat->type = MATRIX_3D;
1275
1276       /* Check for scale */
1277       if (SQ(c1-c2) < SQ(1e-6) && SQ(c1-c3) < SQ(1e-6)) {
1278          if (SQ(c1-1.0) > SQ(1e-6))
1279             mat->flags |= MAT_FLAG_UNIFORM_SCALE;
1280          /* else no scale at all */
1281       }
1282       else {
1283          mat->flags |= MAT_FLAG_GENERAL_SCALE;
1284       }
1285
1286       /* Check for rotation */
1287       if (SQ(d1) < SQ(1e-6)) {
1288          CROSS3( cp, m, m+4 );
1289          SUB_3V( cp, cp, (m+8) );
1290          if (LEN_SQUARED_3FV(cp) < SQ(1e-6))
1291             mat->flags |= MAT_FLAG_ROTATION;
1292          else
1293             mat->flags |= MAT_FLAG_GENERAL_3D;
1294       }
1295       else {
1296          mat->flags |= MAT_FLAG_GENERAL_3D; /* shear, etc */
1297       }
1298    }
1299    else if ((mask & MASK_PERSPECTIVE) == MASK_PERSPECTIVE && m[11]==-1.0F) {
1300       mat->type = MATRIX_PERSPECTIVE;
1301       mat->flags |= MAT_FLAG_GENERAL;
1302    }
1303    else {
1304       mat->type = MATRIX_GENERAL;
1305       mat->flags |= MAT_FLAG_GENERAL;
1306    }
1307 }
1308
1309 /**
1310  * Analyze a matrix given that its flags are accurate.
1311  *
1312  * This is the more common operation, hopefully.
1313  */
1314 static void analyse_from_flags( GLmatrix *mat )
1315 {
1316    const GLfloat *m = mat->m;
1317
1318    if (TEST_MAT_FLAGS(mat, 0)) {
1319       mat->type = MATRIX_IDENTITY;
1320    }
1321    else if (TEST_MAT_FLAGS(mat, (MAT_FLAG_TRANSLATION |
1322                                  MAT_FLAG_UNIFORM_SCALE |
1323                                  MAT_FLAG_GENERAL_SCALE))) {
1324       if ( m[10]==1.0F && m[14]==0.0F ) {
1325          mat->type = MATRIX_2D_NO_ROT;
1326       }
1327       else {
1328          mat->type = MATRIX_3D_NO_ROT;
1329       }
1330    }
1331    else if (TEST_MAT_FLAGS(mat, MAT_FLAGS_3D)) {
1332       if (                                 m[ 8]==0.0F
1333             &&                             m[ 9]==0.0F
1334             && m[2]==0.0F && m[6]==0.0F && m[10]==1.0F && m[14]==0.0F) {
1335          mat->type = MATRIX_2D;
1336       }
1337       else {
1338          mat->type = MATRIX_3D;
1339       }
1340    }
1341    else if (                 m[4]==0.0F                 && m[12]==0.0F
1342             && m[1]==0.0F                               && m[13]==0.0F
1343             && m[2]==0.0F && m[6]==0.0F
1344             && m[3]==0.0F && m[7]==0.0F && m[11]==-1.0F && m[15]==0.0F) {
1345       mat->type = MATRIX_PERSPECTIVE;
1346    }
1347    else {
1348       mat->type = MATRIX_GENERAL;
1349    }
1350 }
1351
1352 /**
1353  * Analyze and update a matrix.
1354  *
1355  * \param mat matrix.
1356  *
1357  * If the matrix type is dirty then calls either analyse_from_scratch() or
1358  * analyse_from_flags() to determine its type, according to whether the flags
1359  * are dirty or not, respectively. If the matrix has an inverse and it's dirty
1360  * then calls matrix_invert(). Finally clears the dirty flags.
1361  */
1362 void
1363 _math_matrix_analyse( GLmatrix *mat )
1364 {
1365    if (mat->flags & MAT_DIRTY_TYPE) {
1366       if (mat->flags & MAT_DIRTY_FLAGS)
1367          analyse_from_scratch( mat );
1368       else
1369          analyse_from_flags( mat );
1370    }
1371
1372    if (mat->inv && (mat->flags & MAT_DIRTY_INVERSE)) {
1373       matrix_invert( mat );
1374       mat->flags &= ~MAT_DIRTY_INVERSE;
1375    }
1376
1377    mat->flags &= ~(MAT_DIRTY_FLAGS | MAT_DIRTY_TYPE);
1378 }
1379
1380 /*@}*/
1381
1382
1383 /**
1384  * Test if the given matrix preserves vector lengths.
1385  */
1386 GLboolean
1387 _math_matrix_is_length_preserving( const GLmatrix *m )
1388 {
1389    return TEST_MAT_FLAGS( m, MAT_FLAGS_LENGTH_PRESERVING);
1390 }
1391
1392
1393 /**
1394  * Test if the given matrix does any rotation.
1395  * (or perhaps if the upper-left 3x3 is non-identity)
1396  */
1397 GLboolean
1398 _math_matrix_has_rotation( const GLmatrix *m )
1399 {
1400    if (m->flags & (MAT_FLAG_GENERAL |
1401                    MAT_FLAG_ROTATION |
1402                    MAT_FLAG_GENERAL_3D |
1403                    MAT_FLAG_PERSPECTIVE))
1404       return GL_TRUE;
1405    else
1406       return GL_FALSE;
1407 }
1408
1409
1410 GLboolean
1411 _math_matrix_is_general_scale( const GLmatrix *m )
1412 {
1413    return (m->flags & MAT_FLAG_GENERAL_SCALE) ? GL_TRUE : GL_FALSE;
1414 }
1415
1416
1417 GLboolean
1418 _math_matrix_is_dirty( const GLmatrix *m )
1419 {
1420    return (m->flags & MAT_DIRTY) ? GL_TRUE : GL_FALSE;
1421 }
1422
1423
1424 /**********************************************************************/
1425 /** \name Matrix setup */
1426 /*@{*/
1427
1428 /**
1429  * Copy a matrix.
1430  *
1431  * \param to destination matrix.
1432  * \param from source matrix.
1433  *
1434  * Copies all fields in GLmatrix, creating an inverse array if necessary.
1435  */
1436 void
1437 _math_matrix_copy( GLmatrix *to, const GLmatrix *from )
1438 {
1439    memcpy( to->m, from->m, sizeof(Identity) );
1440    memcpy(to->inv, from->inv, sizeof(*from->inv));
1441    to->flags = from->flags;
1442    to->type = from->type;
1443 }
1444
1445 /**
1446  * Loads a matrix array into GLmatrix.
1447  *
1448  * \param m matrix array.
1449  * \param mat matrix.
1450  *
1451  * Copies \p m into GLmatrix::m and marks the MAT_FLAG_GENERAL and MAT_DIRTY
1452  * flags.
1453  */
1454 void
1455 _math_matrix_loadf( GLmatrix *mat, const GLfloat *m )
1456 {
1457    memcpy( mat->m, m, 16*sizeof(GLfloat) );
1458    mat->flags = (MAT_FLAG_GENERAL | MAT_DIRTY);
1459 }
1460
1461 /**
1462  * Matrix constructor.
1463  *
1464  * \param m matrix.
1465  *
1466  * Initialize the GLmatrix fields.
1467  */
1468 void
1469 _math_matrix_ctr( GLmatrix *m )
1470 {
1471    m->m = (GLfloat *) _mesa_align_malloc( 16 * sizeof(GLfloat), 16 );
1472    if (m->m)
1473       memcpy( m->m, Identity, sizeof(Identity) );
1474    m->inv = (GLfloat *) _mesa_align_malloc( 16 * sizeof(GLfloat), 16 );
1475    if (m->inv)
1476       memcpy( m->inv, Identity, sizeof(Identity) );
1477    m->type = MATRIX_IDENTITY;
1478    m->flags = 0;
1479 }
1480
1481 /**
1482  * Matrix destructor.
1483  *
1484  * \param m matrix.
1485  *
1486  * Frees the data in a GLmatrix.
1487  */
1488 void
1489 _math_matrix_dtr( GLmatrix *m )
1490 {
1491    if (m->m) {
1492       _mesa_align_free( m->m );
1493       m->m = NULL;
1494    }
1495    if (m->inv) {
1496       _mesa_align_free( m->inv );
1497       m->inv = NULL;
1498    }
1499 }
1500
1501 /*@}*/
1502
1503
1504 /**********************************************************************/
1505 /** \name Matrix transpose */
1506 /*@{*/
1507
1508 /**
1509  * Transpose a GLfloat matrix.
1510  *
1511  * \param to destination array.
1512  * \param from source array.
1513  */
1514 void
1515 _math_transposef( GLfloat to[16], const GLfloat from[16] )
1516 {
1517    to[0] = from[0];
1518    to[1] = from[4];
1519    to[2] = from[8];
1520    to[3] = from[12];
1521    to[4] = from[1];
1522    to[5] = from[5];
1523    to[6] = from[9];
1524    to[7] = from[13];
1525    to[8] = from[2];
1526    to[9] = from[6];
1527    to[10] = from[10];
1528    to[11] = from[14];
1529    to[12] = from[3];
1530    to[13] = from[7];
1531    to[14] = from[11];
1532    to[15] = from[15];
1533 }
1534
1535 /**
1536  * Transpose a GLdouble matrix.
1537  *
1538  * \param to destination array.
1539  * \param from source array.
1540  */
1541 void
1542 _math_transposed( GLdouble to[16], const GLdouble from[16] )
1543 {
1544    to[0] = from[0];
1545    to[1] = from[4];
1546    to[2] = from[8];
1547    to[3] = from[12];
1548    to[4] = from[1];
1549    to[5] = from[5];
1550    to[6] = from[9];
1551    to[7] = from[13];
1552    to[8] = from[2];
1553    to[9] = from[6];
1554    to[10] = from[10];
1555    to[11] = from[14];
1556    to[12] = from[3];
1557    to[13] = from[7];
1558    to[14] = from[11];
1559    to[15] = from[15];
1560 }
1561
1562 /**
1563  * Transpose a GLdouble matrix and convert to GLfloat.
1564  *
1565  * \param to destination array.
1566  * \param from source array.
1567  */
1568 void
1569 _math_transposefd( GLfloat to[16], const GLdouble from[16] )
1570 {
1571    to[0] = (GLfloat) from[0];
1572    to[1] = (GLfloat) from[4];
1573    to[2] = (GLfloat) from[8];
1574    to[3] = (GLfloat) from[12];
1575    to[4] = (GLfloat) from[1];
1576    to[5] = (GLfloat) from[5];
1577    to[6] = (GLfloat) from[9];
1578    to[7] = (GLfloat) from[13];
1579    to[8] = (GLfloat) from[2];
1580    to[9] = (GLfloat) from[6];
1581    to[10] = (GLfloat) from[10];
1582    to[11] = (GLfloat) from[14];
1583    to[12] = (GLfloat) from[3];
1584    to[13] = (GLfloat) from[7];
1585    to[14] = (GLfloat) from[11];
1586    to[15] = (GLfloat) from[15];
1587 }
1588
1589 /*@}*/
1590
1591
1592 /**
1593  * Transform a 4-element row vector (1x4 matrix) by a 4x4 matrix.  This
1594  * function is used for transforming clipping plane equations and spotlight
1595  * directions.
1596  * Mathematically,  u = v * m.
1597  * Input:  v - input vector
1598  *         m - transformation matrix
1599  * Output:  u - transformed vector
1600  */
1601 void
1602 _mesa_transform_vector( GLfloat u[4], const GLfloat v[4], const GLfloat m[16] )
1603 {
1604    const GLfloat v0 = v[0], v1 = v[1], v2 = v[2], v3 = v[3];
1605 #define M(row,col)  m[row + col*4]
1606    u[0] = v0 * M(0,0) + v1 * M(1,0) + v2 * M(2,0) + v3 * M(3,0);
1607    u[1] = v0 * M(0,1) + v1 * M(1,1) + v2 * M(2,1) + v3 * M(3,1);
1608    u[2] = v0 * M(0,2) + v1 * M(1,2) + v2 * M(2,2) + v3 * M(3,2);
1609    u[3] = v0 * M(0,3) + v1 * M(1,3) + v2 * M(2,3) + v3 * M(3,3);
1610 #undef M
1611 }