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 (at) 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 (at) uruk.org). 791 * Optimizations contributed by Rudolf Opalla (rudi (at) 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, 16 * sizeof(GLfloat)); 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 } 1612
