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