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