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