Subversion Repositories shark

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
56 pj 1
/* $Id: m_matrix.c,v 1.1 2003-02-28 11:48:05 pj Exp $ */
2
 
3
/*
4
 * Mesa 3-D graphics library
5
 * Version:  4.1
6
 *
7
 * Copyright (C) 1999-2002  Brian Paul   All Rights Reserved.
8
 *
9
 * Permission is hereby granted, free of charge, to any person obtaining a
10
 * copy of this software and associated documentation files (the "Software"),
11
 * to deal in the Software without restriction, including without limitation
12
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
13
 * and/or sell copies of the Software, and to permit persons to whom the
14
 * Software is furnished to do so, subject to the following conditions:
15
 *
16
 * The above copyright notice and this permission notice shall be included
17
 * in all copies or substantial portions of the Software.
18
 *
19
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
20
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
22
 * BRIAN PAUL BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN
23
 * AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
24
 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
25
 */
26
 
27
 
28
/*
29
 * Matrix operations
30
 *
31
 * NOTES:
32
 * 1. 4x4 transformation matrices are stored in memory in column major order.
33
 * 2. Points/vertices are to be thought of as column vectors.
34
 * 3. Transformation of a point p by a matrix M is: p' = M * p
35
 */
36
 
37
#include "glheader.h"
38
#include "imports.h"
39
#include "macros.h"
40
#include "imports.h"
41
#include "mmath.h"
42
 
43
#include "m_matrix.h"
44
 
45
 
46
static const char *types[] = {
47
   "MATRIX_GENERAL",
48
   "MATRIX_IDENTITY",
49
   "MATRIX_3D_NO_ROT",
50
   "MATRIX_PERSPECTIVE",
51
   "MATRIX_2D",
52
   "MATRIX_2D_NO_ROT",
53
   "MATRIX_3D"
54
};
55
 
56
 
57
static GLfloat Identity[16] = {
58
   1.0, 0.0, 0.0, 0.0,
59
   0.0, 1.0, 0.0, 0.0,
60
   0.0, 0.0, 1.0, 0.0,
61
   0.0, 0.0, 0.0, 1.0
62
};
63
 
64
 
65
 
66
 
67
/*
68
 * This matmul was contributed by Thomas Malik
69
 *
70
 * Perform a 4x4 matrix multiplication  (product = a x b).
71
 * Input:  a, b - matrices to multiply
72
 * Output:  product - product of a and b
73
 * WARNING: (product != b) assumed
74
 * NOTE:    (product == a) allowed
75
 *
76
 * KW: 4*16 = 64 muls
77
 */
78
#define A(row,col)  a[(col<<2)+row]
79
#define B(row,col)  b[(col<<2)+row]
80
#define P(row,col)  product[(col<<2)+row]
81
 
82
static void matmul4( GLfloat *product, const GLfloat *a, const GLfloat *b )
83
{
84
   GLint i;
85
   for (i = 0; i < 4; i++) {
86
      const GLfloat ai0=A(i,0),  ai1=A(i,1),  ai2=A(i,2),  ai3=A(i,3);
87
      P(i,0) = ai0 * B(0,0) + ai1 * B(1,0) + ai2 * B(2,0) + ai3 * B(3,0);
88
      P(i,1) = ai0 * B(0,1) + ai1 * B(1,1) + ai2 * B(2,1) + ai3 * B(3,1);
89
      P(i,2) = ai0 * B(0,2) + ai1 * B(1,2) + ai2 * B(2,2) + ai3 * B(3,2);
90
      P(i,3) = ai0 * B(0,3) + ai1 * B(1,3) + ai2 * B(2,3) + ai3 * B(3,3);
91
   }
92
}
93
 
94
 
95
/* Multiply two matrices known to occupy only the top three rows, such
96
 * as typical model matrices, and ortho matrices.
97
 */
98
static void matmul34( GLfloat *product, const GLfloat *a, const GLfloat *b )
99
{
100
   GLint i;
101
   for (i = 0; i < 3; i++) {
102
      const GLfloat ai0=A(i,0),  ai1=A(i,1),  ai2=A(i,2),  ai3=A(i,3);
103
      P(i,0) = ai0 * B(0,0) + ai1 * B(1,0) + ai2 * B(2,0);
104
      P(i,1) = ai0 * B(0,1) + ai1 * B(1,1) + ai2 * B(2,1);
105
      P(i,2) = ai0 * B(0,2) + ai1 * B(1,2) + ai2 * B(2,2);
106
      P(i,3) = ai0 * B(0,3) + ai1 * B(1,3) + ai2 * B(2,3) + ai3;
107
   }
108
   P(3,0) = 0;
109
   P(3,1) = 0;
110
   P(3,2) = 0;
111
   P(3,3) = 1;
112
}
113
 
114
 
115
#undef A
116
#undef B
117
#undef P
118
 
119
 
120
/*
121
 * Multiply a matrix by an array of floats with known properties.
122
 */
123
static void matrix_multf( GLmatrix *mat, const GLfloat *m, GLuint flags )
124
{
125
   mat->flags |= (flags | MAT_DIRTY_TYPE | MAT_DIRTY_INVERSE);
126
 
127
   if (TEST_MAT_FLAGS(mat, MAT_FLAGS_3D))
128
      matmul34( mat->m, mat->m, m );
129
   else
130
      matmul4( mat->m, mat->m, m );
131
}
132
 
133
 
134
static void print_matrix_floats( const GLfloat m[16] )
135
{
136
   int i;
137
   for (i=0;i<4;i++) {
138
      _mesa_debug(NULL,"\t%f %f %f %f\n", m[i], m[4+i], m[8+i], m[12+i] );
139
   }
140
}
141
 
142
void
143
_math_matrix_print( const GLmatrix *m )
144
{
145
   _mesa_debug(NULL, "Matrix type: %s, flags: %x\n", types[m->type], m->flags);
146
   print_matrix_floats(m->m);
147
   _mesa_debug(NULL, "Inverse: \n");
148
   if (m->inv) {
149
      GLfloat prod[16];
150
      print_matrix_floats(m->inv);
151
      matmul4(prod, m->m, m->inv);
152
      _mesa_debug(NULL, "Mat * Inverse:\n");
153
      print_matrix_floats(prod);
154
   }
155
   else {
156
      _mesa_debug(NULL, "  - not available\n");
157
   }
158
}
159
 
160
 
161
 
162
 
163
#define SWAP_ROWS(a, b) { GLfloat *_tmp = a; (a)=(b); (b)=_tmp; }
164
#define MAT(m,r,c) (m)[(c)*4+(r)]
165
 
166
/*
167
 * Compute inverse of 4x4 transformation matrix.
168
 * Code contributed by Jacques Leroy jle@star.be
169
 * Return GL_TRUE for success, GL_FALSE for failure (singular matrix)
170
 */
171
static GLboolean invert_matrix_general( GLmatrix *mat )
172
{
173
   const GLfloat *m = mat->m;
174
   GLfloat *out = mat->inv;
175
   GLfloat wtmp[4][8];
176
   GLfloat m0, m1, m2, m3, s;
177
   GLfloat *r0, *r1, *r2, *r3;
178
 
179
   r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];
180
 
181
   r0[0] = MAT(m,0,0), r0[1] = MAT(m,0,1),
182
   r0[2] = MAT(m,0,2), r0[3] = MAT(m,0,3),
183
   r0[4] = 1.0, r0[5] = r0[6] = r0[7] = 0.0,
184
 
185
   r1[0] = MAT(m,1,0), r1[1] = MAT(m,1,1),
186
   r1[2] = MAT(m,1,2), r1[3] = MAT(m,1,3),
187
   r1[5] = 1.0, r1[4] = r1[6] = r1[7] = 0.0,
188
 
189
   r2[0] = MAT(m,2,0), r2[1] = MAT(m,2,1),
190
   r2[2] = MAT(m,2,2), r2[3] = MAT(m,2,3),
191
   r2[6] = 1.0, r2[4] = r2[5] = r2[7] = 0.0,
192
 
193
   r3[0] = MAT(m,3,0), r3[1] = MAT(m,3,1),
194
   r3[2] = MAT(m,3,2), r3[3] = MAT(m,3,3),
195
   r3[7] = 1.0, r3[4] = r3[5] = r3[6] = 0.0;
196
 
197
   /* choose pivot - or die */
198
   if (fabs(r3[0])>fabs(r2[0])) SWAP_ROWS(r3, r2);
199
   if (fabs(r2[0])>fabs(r1[0])) SWAP_ROWS(r2, r1);
200
   if (fabs(r1[0])>fabs(r0[0])) SWAP_ROWS(r1, r0);
201
   if (0.0 == r0[0])  return GL_FALSE;
202
 
203
   /* eliminate first variable     */
204
   m1 = r1[0]/r0[0]; m2 = r2[0]/r0[0]; m3 = r3[0]/r0[0];
205
   s = r0[1]; r1[1] -= m1 * s; r2[1] -= m2 * s; r3[1] -= m3 * s;
206
   s = r0[2]; r1[2] -= m1 * s; r2[2] -= m2 * s; r3[2] -= m3 * s;
207
   s = r0[3]; r1[3] -= m1 * s; r2[3] -= m2 * s; r3[3] -= m3 * s;
208
   s = r0[4];
209
   if (s != 0.0) { r1[4] -= m1 * s; r2[4] -= m2 * s; r3[4] -= m3 * s; }
210
   s = r0[5];
211
   if (s != 0.0) { r1[5] -= m1 * s; r2[5] -= m2 * s; r3[5] -= m3 * s; }
212
   s = r0[6];
213
   if (s != 0.0) { r1[6] -= m1 * s; r2[6] -= m2 * s; r3[6] -= m3 * s; }
214
   s = r0[7];
215
   if (s != 0.0) { r1[7] -= m1 * s; r2[7] -= m2 * s; r3[7] -= m3 * s; }
216
 
217
   /* choose pivot - or die */
218
   if (fabs(r3[1])>fabs(r2[1])) SWAP_ROWS(r3, r2);
219
   if (fabs(r2[1])>fabs(r1[1])) SWAP_ROWS(r2, r1);
220
   if (0.0 == r1[1])  return GL_FALSE;
221
 
222
   /* eliminate second variable */
223
   m2 = r2[1]/r1[1]; m3 = r3[1]/r1[1];
224
   r2[2] -= m2 * r1[2]; r3[2] -= m3 * r1[2];
225
   r2[3] -= m2 * r1[3]; r3[3] -= m3 * r1[3];
226
   s = r1[4]; if (0.0 != s) { r2[4] -= m2 * s; r3[4] -= m3 * s; }
227
   s = r1[5]; if (0.0 != s) { r2[5] -= m2 * s; r3[5] -= m3 * s; }
228
   s = r1[6]; if (0.0 != s) { r2[6] -= m2 * s; r3[6] -= m3 * s; }
229
   s = r1[7]; if (0.0 != s) { r2[7] -= m2 * s; r3[7] -= m3 * s; }
230
 
231
   /* choose pivot - or die */
232
   if (fabs(r3[2])>fabs(r2[2])) SWAP_ROWS(r3, r2);
233
   if (0.0 == r2[2])  return GL_FALSE;
234
 
235
   /* eliminate third variable */
236
   m3 = r3[2]/r2[2];
237
   r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
238
   r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6],
239
   r3[7] -= m3 * r2[7];
240
 
241
   /* last check */
242
   if (0.0 == r3[3]) return GL_FALSE;
243
 
244
   s = 1.0F/r3[3];             /* now back substitute row 3 */
245
   r3[4] *= s; r3[5] *= s; r3[6] *= s; r3[7] *= s;
246
 
247
   m2 = r2[3];                 /* now back substitute row 2 */
248
   s  = 1.0F/r2[2];
249
   r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
250
   r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
251
   m1 = r1[3];
252
   r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
253
   r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
254
   m0 = r0[3];
255
   r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
256
   r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;
257
 
258
   m1 = r1[2];                 /* now back substitute row 1 */
259
   s  = 1.0F/r1[1];
260
   r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
261
   r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
262
   m0 = r0[2];
263
   r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
264
   r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;
265
 
266
   m0 = r0[1];                 /* now back substitute row 0 */
267
   s  = 1.0F/r0[0];
268
   r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
269
   r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);
270
 
271
   MAT(out,0,0) = r0[4]; MAT(out,0,1) = r0[5],
272
   MAT(out,0,2) = r0[6]; MAT(out,0,3) = r0[7],
273
   MAT(out,1,0) = r1[4]; MAT(out,1,1) = r1[5],
274
   MAT(out,1,2) = r1[6]; MAT(out,1,3) = r1[7],
275
   MAT(out,2,0) = r2[4]; MAT(out,2,1) = r2[5],
276
   MAT(out,2,2) = r2[6]; MAT(out,2,3) = r2[7],
277
   MAT(out,3,0) = r3[4]; MAT(out,3,1) = r3[5],
278
   MAT(out,3,2) = r3[6]; MAT(out,3,3) = r3[7];
279
 
280
   return GL_TRUE;
281
}
282
#undef SWAP_ROWS
283
 
284
 
285
/* Adapted from graphics gems II.
286
 */
287
static GLboolean invert_matrix_3d_general( GLmatrix *mat )
288
{
289
   const GLfloat *in = mat->m;
290
   GLfloat *out = mat->inv;
291
   GLfloat pos, neg, t;
292
   GLfloat det;
293
 
294
   /* Calculate the determinant of upper left 3x3 submatrix and
295
    * determine if the matrix is singular.
296
    */
297
   pos = neg = 0.0;
298
   t =  MAT(in,0,0) * MAT(in,1,1) * MAT(in,2,2);
299
   if (t >= 0.0) pos += t; else neg += t;
300
 
301
   t =  MAT(in,1,0) * MAT(in,2,1) * MAT(in,0,2);
302
   if (t >= 0.0) pos += t; else neg += t;
303
 
304
   t =  MAT(in,2,0) * MAT(in,0,1) * MAT(in,1,2);
305
   if (t >= 0.0) pos += t; else neg += t;
306
 
307
   t = -MAT(in,2,0) * MAT(in,1,1) * MAT(in,0,2);
308
   if (t >= 0.0) pos += t; else neg += t;
309
 
310
   t = -MAT(in,1,0) * MAT(in,0,1) * MAT(in,2,2);
311
   if (t >= 0.0) pos += t; else neg += t;
312
 
313
   t = -MAT(in,0,0) * MAT(in,2,1) * MAT(in,1,2);
314
   if (t >= 0.0) pos += t; else neg += t;
315
 
316
   det = pos + neg;
317
 
318
   if (det*det < 1e-25)
319
      return GL_FALSE;
320
 
321
   det = 1.0F / det;
322
   MAT(out,0,0) = (  (MAT(in,1,1)*MAT(in,2,2) - MAT(in,2,1)*MAT(in,1,2) )*det);
323
   MAT(out,0,1) = (- (MAT(in,0,1)*MAT(in,2,2) - MAT(in,2,1)*MAT(in,0,2) )*det);
324
   MAT(out,0,2) = (  (MAT(in,0,1)*MAT(in,1,2) - MAT(in,1,1)*MAT(in,0,2) )*det);
325
   MAT(out,1,0) = (- (MAT(in,1,0)*MAT(in,2,2) - MAT(in,2,0)*MAT(in,1,2) )*det);
326
   MAT(out,1,1) = (  (MAT(in,0,0)*MAT(in,2,2) - MAT(in,2,0)*MAT(in,0,2) )*det);
327
   MAT(out,1,2) = (- (MAT(in,0,0)*MAT(in,1,2) - MAT(in,1,0)*MAT(in,0,2) )*det);
328
   MAT(out,2,0) = (  (MAT(in,1,0)*MAT(in,2,1) - MAT(in,2,0)*MAT(in,1,1) )*det);
329
   MAT(out,2,1) = (- (MAT(in,0,0)*MAT(in,2,1) - MAT(in,2,0)*MAT(in,0,1) )*det);
330
   MAT(out,2,2) = (  (MAT(in,0,0)*MAT(in,1,1) - MAT(in,1,0)*MAT(in,0,1) )*det);
331
 
332
   /* Do the translation part */
333
   MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0) +
334
                     MAT(in,1,3) * MAT(out,0,1) +
335
                     MAT(in,2,3) * MAT(out,0,2) );
336
   MAT(out,1,3) = - (MAT(in,0,3) * MAT(out,1,0) +
337
                     MAT(in,1,3) * MAT(out,1,1) +
338
                     MAT(in,2,3) * MAT(out,1,2) );
339
   MAT(out,2,3) = - (MAT(in,0,3) * MAT(out,2,0) +
340
                     MAT(in,1,3) * MAT(out,2,1) +
341
                     MAT(in,2,3) * MAT(out,2,2) );
342
 
343
   return GL_TRUE;
344
}
345
 
346
 
347
static GLboolean invert_matrix_3d( GLmatrix *mat )
348
{
349
   const GLfloat *in = mat->m;
350
   GLfloat *out = mat->inv;
351
 
352
   if (!TEST_MAT_FLAGS(mat, MAT_FLAGS_ANGLE_PRESERVING)) {
353
      return invert_matrix_3d_general( mat );
354
   }
355
 
356
   if (mat->flags & MAT_FLAG_UNIFORM_SCALE) {
357
      GLfloat scale = (MAT(in,0,0) * MAT(in,0,0) +
358
                       MAT(in,0,1) * MAT(in,0,1) +
359
                       MAT(in,0,2) * MAT(in,0,2));
360
 
361
      if (scale == 0.0)
362
         return GL_FALSE;
363
 
364
      scale = 1.0F / scale;
365
 
366
      /* Transpose and scale the 3 by 3 upper-left submatrix. */
367
      MAT(out,0,0) = scale * MAT(in,0,0);
368
      MAT(out,1,0) = scale * MAT(in,0,1);
369
      MAT(out,2,0) = scale * MAT(in,0,2);
370
      MAT(out,0,1) = scale * MAT(in,1,0);
371
      MAT(out,1,1) = scale * MAT(in,1,1);
372
      MAT(out,2,1) = scale * MAT(in,1,2);
373
      MAT(out,0,2) = scale * MAT(in,2,0);
374
      MAT(out,1,2) = scale * MAT(in,2,1);
375
      MAT(out,2,2) = scale * MAT(in,2,2);
376
   }
377
   else if (mat->flags & MAT_FLAG_ROTATION) {
378
      /* Transpose the 3 by 3 upper-left submatrix. */
379
      MAT(out,0,0) = MAT(in,0,0);
380
      MAT(out,1,0) = MAT(in,0,1);
381
      MAT(out,2,0) = MAT(in,0,2);
382
      MAT(out,0,1) = MAT(in,1,0);
383
      MAT(out,1,1) = MAT(in,1,1);
384
      MAT(out,2,1) = MAT(in,1,2);
385
      MAT(out,0,2) = MAT(in,2,0);
386
      MAT(out,1,2) = MAT(in,2,1);
387
      MAT(out,2,2) = MAT(in,2,2);
388
   }
389
   else {
390
      /* pure translation */
391
      MEMCPY( out, Identity, sizeof(Identity) );
392
      MAT(out,0,3) = - MAT(in,0,3);
393
      MAT(out,1,3) = - MAT(in,1,3);
394
      MAT(out,2,3) = - MAT(in,2,3);
395
      return GL_TRUE;
396
   }
397
 
398
   if (mat->flags & MAT_FLAG_TRANSLATION) {
399
      /* Do the translation part */
400
      MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0) +
401
                        MAT(in,1,3) * MAT(out,0,1) +
402
                        MAT(in,2,3) * MAT(out,0,2) );
403
      MAT(out,1,3) = - (MAT(in,0,3) * MAT(out,1,0) +
404
                        MAT(in,1,3) * MAT(out,1,1) +
405
                        MAT(in,2,3) * MAT(out,1,2) );
406
      MAT(out,2,3) = - (MAT(in,0,3) * MAT(out,2,0) +
407
                        MAT(in,1,3) * MAT(out,2,1) +
408
                        MAT(in,2,3) * MAT(out,2,2) );
409
   }
410
   else {
411
      MAT(out,0,3) = MAT(out,1,3) = MAT(out,2,3) = 0.0;
412
   }
413
 
414
   return GL_TRUE;
415
}
416
 
417
 
418
 
419
static GLboolean invert_matrix_identity( GLmatrix *mat )
420
{
421
   MEMCPY( mat->inv, Identity, sizeof(Identity) );
422
   return GL_TRUE;
423
}
424
 
425
 
426
static GLboolean invert_matrix_3d_no_rot( GLmatrix *mat )
427
{
428
   const GLfloat *in = mat->m;
429
   GLfloat *out = mat->inv;
430
 
431
   if (MAT(in,0,0) == 0 || MAT(in,1,1) == 0 || MAT(in,2,2) == 0 )
432
      return GL_FALSE;
433
 
434
   MEMCPY( out, Identity, 16 * sizeof(GLfloat) );
435
   MAT(out,0,0) = 1.0F / MAT(in,0,0);
436
   MAT(out,1,1) = 1.0F / MAT(in,1,1);
437
   MAT(out,2,2) = 1.0F / MAT(in,2,2);
438
 
439
   if (mat->flags & MAT_FLAG_TRANSLATION) {
440
      MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0));
441
      MAT(out,1,3) = - (MAT(in,1,3) * MAT(out,1,1));
442
      MAT(out,2,3) = - (MAT(in,2,3) * MAT(out,2,2));
443
   }
444
 
445
   return GL_TRUE;
446
}
447
 
448
 
449
static GLboolean invert_matrix_2d_no_rot( GLmatrix *mat )
450
{
451
   const GLfloat *in = mat->m;
452
   GLfloat *out = mat->inv;
453
 
454
   if (MAT(in,0,0) == 0 || MAT(in,1,1) == 0)
455
      return GL_FALSE;
456
 
457
   MEMCPY( out, Identity, 16 * sizeof(GLfloat) );
458
   MAT(out,0,0) = 1.0F / MAT(in,0,0);
459
   MAT(out,1,1) = 1.0F / MAT(in,1,1);
460
 
461
   if (mat->flags & MAT_FLAG_TRANSLATION) {
462
      MAT(out,0,3) = - (MAT(in,0,3) * MAT(out,0,0));
463
      MAT(out,1,3) = - (MAT(in,1,3) * MAT(out,1,1));
464
   }
465
 
466
   return GL_TRUE;
467
}
468
 
469
 
470
#if 0
471
/* broken */
472
static GLboolean invert_matrix_perspective( GLmatrix *mat )
473
{
474
   const GLfloat *in = mat->m;
475
   GLfloat *out = mat->inv;
476
 
477
   if (MAT(in,2,3) == 0)
478
      return GL_FALSE;
479
 
480
   MEMCPY( out, Identity, 16 * sizeof(GLfloat) );
481
 
482
   MAT(out,0,0) = 1.0F / MAT(in,0,0);
483
   MAT(out,1,1) = 1.0F / MAT(in,1,1);
484
 
485
   MAT(out,0,3) = MAT(in,0,2);
486
   MAT(out,1,3) = MAT(in,1,2);
487
 
488
   MAT(out,2,2) = 0;
489
   MAT(out,2,3) = -1;
490
 
491
   MAT(out,3,2) = 1.0F / MAT(in,2,3);
492
   MAT(out,3,3) = MAT(in,2,2) * MAT(out,3,2);
493
 
494
   return GL_TRUE;
495
}
496
#endif
497
 
498
 
499
typedef GLboolean (*inv_mat_func)( GLmatrix *mat );
500
 
501
 
502
static inv_mat_func inv_mat_tab[7] = {
503
   invert_matrix_general,
504
   invert_matrix_identity,
505
   invert_matrix_3d_no_rot,
506
#if 0
507
   /* Don't use this function for now - it fails when the projection matrix
508
    * is premultiplied by a translation (ala Chromium's tilesort SPU).
509
    */
510
   invert_matrix_perspective,
511
#else
512
   invert_matrix_general,
513
#endif
514
   invert_matrix_3d,            /* lazy! */
515
   invert_matrix_2d_no_rot,
516
   invert_matrix_3d
517
};
518
 
519
 
520
static GLboolean matrix_invert( GLmatrix *mat )
521
{
522
   if (inv_mat_tab[mat->type](mat)) {
523
      mat->flags &= ~MAT_FLAG_SINGULAR;
524
      return GL_TRUE;
525
   } else {
526
      mat->flags |= MAT_FLAG_SINGULAR;
527
      MEMCPY( mat->inv, Identity, sizeof(Identity) );
528
      return GL_FALSE;
529
   }
530
}
531
 
532
 
533
 
534
 
535
 
536
 
537
/*
538
 * Generate a 4x4 transformation matrix from glRotate parameters, and
539
 * postmultiply the input matrix by it.
540
 * This function contributed by Erich Boleyn (erich@uruk.org).
541
 * Optimizatios contributed by Rudolf Opalla (rudi@khm.de).
542
 */
543
void
544
_math_matrix_rotate( GLmatrix *mat,
545
                     GLfloat angle, GLfloat x, GLfloat y, GLfloat z )
546
{
547
   GLfloat xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c, s, c;
548
   GLfloat m[16];
549
   GLboolean optimized;
550
 
551
   s = (GLfloat) sin( angle * DEG2RAD );
552
   c = (GLfloat) cos( angle * DEG2RAD );
553
 
554
   MEMCPY(m, Identity, sizeof(GLfloat)*16);
555
   optimized = GL_FALSE;
556
 
557
#define M(row,col)  m[col*4+row]
558
 
559
   if (x == 0.0F) {
560
      if (y == 0.0F) {
561
         if (z != 0.0F) {
562
            optimized = GL_TRUE;
563
            /* rotate only around z-axis */
564
            M(0,0) = c;
565
            M(1,1) = c;
566
            if (z < 0.0F) {
567
               M(0,1) = s;
568
               M(1,0) = -s;
569
            }
570
            else {
571
               M(0,1) = -s;
572
               M(1,0) = s;
573
            }
574
         }
575
      }
576
      else if (z == 0.0F) {
577
         optimized = GL_TRUE;
578
         /* rotate only around y-axis */
579
         M(0,0) = c;
580
         M(2,2) = c;
581
         if (y < 0.0F) {
582
            M(0,2) = -s;
583
            M(2,0) = s;
584
         }
585
         else {
586
            M(0,2) = s;
587
            M(2,0) = -s;
588
         }
589
      }
590
   }
591
   else if (y == 0.0F) {
592
      if (z == 0.0F) {
593
         optimized = GL_TRUE;
594
         /* rotate only around x-axis */
595
         M(1,1) = c;
596
         M(2,2) = c;
597
         if (y < 0.0F) {
598
            M(1,2) = s;
599
            M(2,1) = -s;
600
         }
601
         else {
602
            M(1,2) = -s;
603
            M(2,1) = s;
604
         }
605
      }
606
   }
607
 
608
   if (!optimized) {
609
      const GLfloat mag = (GLfloat) GL_SQRT(x * x + y * y + z * z);
610
 
611
      if (mag <= 1.0e-4) {
612
         /* no rotation, leave mat as-is */
613
         return;
614
      }
615
 
616
      x /= mag;
617
      y /= mag;
618
      z /= mag;
619
 
620
 
621
      /*
622
       *     Arbitrary axis rotation matrix.
623
       *
624
       *  This is composed of 5 matrices, Rz, Ry, T, Ry', Rz', multiplied
625
       *  like so:  Rz * Ry * T * Ry' * Rz'.  T is the final rotation
626
       *  (which is about the X-axis), and the two composite transforms
627
       *  Ry' * Rz' and Rz * Ry are (respectively) the rotations necessary
628
       *  from the arbitrary axis to the X-axis then back.  They are
629
       *  all elementary rotations.
630
       *
631
       *  Rz' is a rotation about the Z-axis, to bring the axis vector
632
       *  into the x-z plane.  Then Ry' is applied, rotating about the
633
       *  Y-axis to bring the axis vector parallel with the X-axis.  The
634
       *  rotation about the X-axis is then performed.  Ry and Rz are
635
       *  simply the respective inverse transforms to bring the arbitrary
636
       *  axis back to it's original orientation.  The first transforms
637
       *  Rz' and Ry' are considered inverses, since the data from the
638
       *  arbitrary axis gives you info on how to get to it, not how
639
       *  to get away from it, and an inverse must be applied.
640
       *
641
       *  The basic calculation used is to recognize that the arbitrary
642
       *  axis vector (x, y, z), since it is of unit length, actually
643
       *  represents the sines and cosines of the angles to rotate the
644
       *  X-axis to the same orientation, with theta being the angle about
645
       *  Z and phi the angle about Y (in the order described above)
646
       *  as follows:
647
       *
648
       *  cos ( theta ) = x / sqrt ( 1 - z^2 )
649
       *  sin ( theta ) = y / sqrt ( 1 - z^2 )
650
       *
651
       *  cos ( phi ) = sqrt ( 1 - z^2 )
652
       *  sin ( phi ) = z
653
       *
654
       *  Note that cos ( phi ) can further be inserted to the above
655
       *  formulas:
656
       *
657
       *  cos ( theta ) = x / cos ( phi )
658
       *  sin ( theta ) = y / sin ( phi )
659
       *
660
       *  ...etc.  Because of those relations and the standard trigonometric
661
       *  relations, it is pssible to reduce the transforms down to what
662
       *  is used below.  It may be that any primary axis chosen will give the
663
       *  same results (modulo a sign convention) using thie method.
664
       *
665
       *  Particularly nice is to notice that all divisions that might
666
       *  have caused trouble when parallel to certain planes or
667
       *  axis go away with care paid to reducing the expressions.
668
       *  After checking, it does perform correctly under all cases, since
669
       *  in all the cases of division where the denominator would have
670
       *  been zero, the numerator would have been zero as well, giving
671
       *  the expected result.
672
       */
673
 
674
      xx = x * x;
675
      yy = y * y;
676
      zz = z * z;
677
      xy = x * y;
678
      yz = y * z;
679
      zx = z * x;
680
      xs = x * s;
681
      ys = y * s;
682
      zs = z * s;
683
      one_c = 1.0F - c;
684
 
685
      /* We already hold the identity-matrix so we can skip some statements */
686
      M(0,0) = (one_c * xx) + c;
687
      M(0,1) = (one_c * xy) - zs;
688
      M(0,2) = (one_c * zx) + ys;
689
/*    M(0,3) = 0.0F; */
690
 
691
      M(1,0) = (one_c * xy) + zs;
692
      M(1,1) = (one_c * yy) + c;
693
      M(1,2) = (one_c * yz) - xs;
694
/*    M(1,3) = 0.0F; */
695
 
696
      M(2,0) = (one_c * zx) - ys;
697
      M(2,1) = (one_c * yz) + xs;
698
      M(2,2) = (one_c * zz) + c;
699
/*    M(2,3) = 0.0F; */
700
 
701
/*
702
      M(3,0) = 0.0F;
703
      M(3,1) = 0.0F;
704
      M(3,2) = 0.0F;
705
      M(3,3) = 1.0F;
706
*/
707
   }
708
#undef M
709
 
710
   matrix_multf( mat, m, MAT_FLAG_ROTATION );
711
}
712
 
713
 
714
 
715
void
716
_math_matrix_frustum( GLmatrix *mat,
717
                      GLfloat left, GLfloat right,
718
                      GLfloat bottom, GLfloat top,
719
                      GLfloat nearval, GLfloat farval )
720
{
721
   GLfloat x, y, a, b, c, d;
722
   GLfloat m[16];
723
 
724
   x = (2.0F*nearval) / (right-left);
725
   y = (2.0F*nearval) / (top-bottom);
726
   a = (right+left) / (right-left);
727
   b = (top+bottom) / (top-bottom);
728
   c = -(farval+nearval) / ( farval-nearval);
729
   d = -(2.0F*farval*nearval) / (farval-nearval);  /* error? */
730
 
731
#define M(row,col)  m[col*4+row]
732
   M(0,0) = x;     M(0,1) = 0.0F;  M(0,2) = a;      M(0,3) = 0.0F;
733
   M(1,0) = 0.0F;  M(1,1) = y;     M(1,2) = b;      M(1,3) = 0.0F;
734
   M(2,0) = 0.0F;  M(2,1) = 0.0F;  M(2,2) = c;      M(2,3) = d;
735
   M(3,0) = 0.0F;  M(3,1) = 0.0F;  M(3,2) = -1.0F;  M(3,3) = 0.0F;
736
#undef M
737
 
738
   matrix_multf( mat, m, MAT_FLAG_PERSPECTIVE );
739
}
740
 
741
void
742
_math_matrix_ortho( GLmatrix *mat,
743
                    GLfloat left, GLfloat right,
744
                    GLfloat bottom, GLfloat top,
745
                    GLfloat nearval, GLfloat farval )
746
{
747
   GLfloat x, y, z;
748
   GLfloat tx, ty, tz;
749
   GLfloat m[16];
750
 
751
   x = 2.0F / (right-left);
752
   y = 2.0F / (top-bottom);
753
   z = -2.0F / (farval-nearval);
754
   tx = -(right+left) / (right-left);
755
   ty = -(top+bottom) / (top-bottom);
756
   tz = -(farval+nearval) / (farval-nearval);
757
 
758
#define M(row,col)  m[col*4+row]
759
   M(0,0) = x;     M(0,1) = 0.0F;  M(0,2) = 0.0F;  M(0,3) = tx;
760
   M(1,0) = 0.0F;  M(1,1) = y;     M(1,2) = 0.0F;  M(1,3) = ty;
761
   M(2,0) = 0.0F;  M(2,1) = 0.0F;  M(2,2) = z;     M(2,3) = tz;
762
   M(3,0) = 0.0F;  M(3,1) = 0.0F;  M(3,2) = 0.0F;  M(3,3) = 1.0F;
763
#undef M
764
 
765
   matrix_multf( mat, m, (MAT_FLAG_GENERAL_SCALE|MAT_FLAG_TRANSLATION));
766
}
767
 
768
 
769
#define ZERO(x) (1<<x)
770
#define ONE(x)  (1<<(x+16))
771
 
772
#define MASK_NO_TRX      (ZERO(12) | ZERO(13) | ZERO(14))
773
#define MASK_NO_2D_SCALE ( ONE(0)  | ONE(5))
774
 
775
#define MASK_IDENTITY    ( ONE(0)  | ZERO(4)  | ZERO(8)  | ZERO(12) |\
776
                          ZERO(1)  |  ONE(5)  | ZERO(9)  | ZERO(13) |\
777
                          ZERO(2)  | ZERO(6)  |  ONE(10) | ZERO(14) |\
778
                          ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
779
 
780
#define MASK_2D_NO_ROT   (           ZERO(4)  | ZERO(8)  |           \
781
                          ZERO(1)  |            ZERO(9)  |           \
782
                          ZERO(2)  | ZERO(6)  |  ONE(10) | ZERO(14) |\
783
                          ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
784
 
785
#define MASK_2D          (                      ZERO(8)  |           \
786
                                                ZERO(9)  |           \
787
                          ZERO(2)  | ZERO(6)  |  ONE(10) | ZERO(14) |\
788
                          ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
789
 
790
 
791
#define MASK_3D_NO_ROT   (           ZERO(4)  | ZERO(8)  |           \
792
                          ZERO(1)  |            ZERO(9)  |           \
793
                          ZERO(2)  | ZERO(6)  |                      \
794
                          ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
795
 
796
#define MASK_3D          (                                           \
797
                                                                     \
798
                                                                     \
799
                          ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
800
 
801
 
802
#define MASK_PERSPECTIVE (           ZERO(4)  |            ZERO(12) |\
803
                          ZERO(1)  |                       ZERO(13) |\
804
                          ZERO(2)  | ZERO(6)  |                      \
805
                          ZERO(3)  | ZERO(7)  |            ZERO(15) )
806
 
807
#define SQ(x) ((x)*(x))
808
 
809
/* Determine type and flags from scratch.  This is expensive enough to
810
 * only want to do it once.
811
 */
812
static void analyse_from_scratch( GLmatrix *mat )
813
{
814
   const GLfloat *m = mat->m;
815
   GLuint mask = 0;
816
   GLuint i;
817
 
818
   for (i = 0 ; i < 16 ; i++) {
819
      if (m[i] == 0.0) mask |= (1<<i);
820
   }
821
 
822
   if (m[0] == 1.0F) mask |= (1<<16);
823
   if (m[5] == 1.0F) mask |= (1<<21);
824
   if (m[10] == 1.0F) mask |= (1<<26);
825
   if (m[15] == 1.0F) mask |= (1<<31);
826
 
827
   mat->flags &= ~MAT_FLAGS_GEOMETRY;
828
 
829
   /* Check for translation - no-one really cares
830
    */
831
   if ((mask & MASK_NO_TRX) != MASK_NO_TRX)
832
      mat->flags |= MAT_FLAG_TRANSLATION;
833
 
834
   /* Do the real work
835
    */
836
   if (mask == (GLuint) MASK_IDENTITY) {
837
      mat->type = MATRIX_IDENTITY;
838
   }
839
   else if ((mask & MASK_2D_NO_ROT) == (GLuint) MASK_2D_NO_ROT) {
840
      mat->type = MATRIX_2D_NO_ROT;
841
 
842
      if ((mask & MASK_NO_2D_SCALE) != MASK_NO_2D_SCALE)
843
         mat->flags = MAT_FLAG_GENERAL_SCALE;
844
   }
845
   else if ((mask & MASK_2D) == (GLuint) MASK_2D) {
846
      GLfloat mm = DOT2(m, m);
847
      GLfloat m4m4 = DOT2(m+4,m+4);
848
      GLfloat mm4 = DOT2(m,m+4);
849
 
850
      mat->type = MATRIX_2D;
851
 
852
      /* Check for scale */
853
      if (SQ(mm-1) > SQ(1e-6) ||
854
          SQ(m4m4-1) > SQ(1e-6))
855
         mat->flags |= MAT_FLAG_GENERAL_SCALE;
856
 
857
      /* Check for rotation */
858
      if (SQ(mm4) > SQ(1e-6))
859
         mat->flags |= MAT_FLAG_GENERAL_3D;
860
      else
861
         mat->flags |= MAT_FLAG_ROTATION;
862
 
863
   }
864
   else if ((mask & MASK_3D_NO_ROT) == (GLuint) MASK_3D_NO_ROT) {
865
      mat->type = MATRIX_3D_NO_ROT;
866
 
867
      /* Check for scale */
868
      if (SQ(m[0]-m[5]) < SQ(1e-6) &&
869
          SQ(m[0]-m[10]) < SQ(1e-6)) {
870
         if (SQ(m[0]-1.0) > SQ(1e-6)) {
871
            mat->flags |= MAT_FLAG_UNIFORM_SCALE;
872
         }
873
      }
874
      else {
875
         mat->flags |= MAT_FLAG_GENERAL_SCALE;
876
      }
877
   }
878
   else if ((mask & MASK_3D) == (GLuint) MASK_3D) {
879
      GLfloat c1 = DOT3(m,m);
880
      GLfloat c2 = DOT3(m+4,m+4);
881
      GLfloat c3 = DOT3(m+8,m+8);
882
      GLfloat d1 = DOT3(m, m+4);
883
      GLfloat cp[3];
884
 
885
      mat->type = MATRIX_3D;
886
 
887
      /* Check for scale */
888
      if (SQ(c1-c2) < SQ(1e-6) && SQ(c1-c3) < SQ(1e-6)) {
889
         if (SQ(c1-1.0) > SQ(1e-6))
890
            mat->flags |= MAT_FLAG_UNIFORM_SCALE;
891
         /* else no scale at all */
892
      }
893
      else {
894
         mat->flags |= MAT_FLAG_GENERAL_SCALE;
895
      }
896
 
897
      /* Check for rotation */
898
      if (SQ(d1) < SQ(1e-6)) {
899
         CROSS3( cp, m, m+4 );
900
         SUB_3V( cp, cp, (m+8) );
901
         if (LEN_SQUARED_3FV(cp) < SQ(1e-6))
902
            mat->flags |= MAT_FLAG_ROTATION;
903
         else
904
            mat->flags |= MAT_FLAG_GENERAL_3D;
905
      }
906
      else {
907
         mat->flags |= MAT_FLAG_GENERAL_3D; /* shear, etc */
908
      }
909
   }
910
   else if ((mask & MASK_PERSPECTIVE) == MASK_PERSPECTIVE && m[11]==-1.0F) {
911
      mat->type = MATRIX_PERSPECTIVE;
912
      mat->flags |= MAT_FLAG_GENERAL;
913
   }
914
   else {
915
      mat->type = MATRIX_GENERAL;
916
      mat->flags |= MAT_FLAG_GENERAL;
917
   }
918
}
919
 
920
 
921
/* Analyse a matrix given that its flags are accurate - this is the
922
 * more common operation, hopefully.
923
 */
924
static void analyse_from_flags( GLmatrix *mat )
925
{
926
   const GLfloat *m = mat->m;
927
 
928
   if (TEST_MAT_FLAGS(mat, 0)) {
929
      mat->type = MATRIX_IDENTITY;
930
   }
931
   else if (TEST_MAT_FLAGS(mat, (MAT_FLAG_TRANSLATION |
932
                                 MAT_FLAG_UNIFORM_SCALE |
933
                                 MAT_FLAG_GENERAL_SCALE))) {
934
      if ( m[10]==1.0F && m[14]==0.0F ) {
935
         mat->type = MATRIX_2D_NO_ROT;
936
      }
937
      else {
938
         mat->type = MATRIX_3D_NO_ROT;
939
      }
940
   }
941
   else if (TEST_MAT_FLAGS(mat, MAT_FLAGS_3D)) {
942
      if (                                 m[ 8]==0.0F
943
            &&                             m[ 9]==0.0F
944
            && m[2]==0.0F && m[6]==0.0F && m[10]==1.0F && m[14]==0.0F) {
945
         mat->type = MATRIX_2D;
946
      }
947
      else {
948
         mat->type = MATRIX_3D;
949
      }
950
   }
951
   else if (                 m[4]==0.0F                 && m[12]==0.0F
952
            && m[1]==0.0F                               && m[13]==0.0F
953
            && m[2]==0.0F && m[6]==0.0F
954
            && m[3]==0.0F && m[7]==0.0F && m[11]==-1.0F && m[15]==0.0F) {
955
      mat->type = MATRIX_PERSPECTIVE;
956
   }
957
   else {
958
      mat->type = MATRIX_GENERAL;
959
   }
960
}
961
 
962
 
963
void
964
_math_matrix_analyse( GLmatrix *mat )
965
{
966
   if (mat->flags & MAT_DIRTY_TYPE) {
967
      if (mat->flags & MAT_DIRTY_FLAGS)
968
         analyse_from_scratch( mat );
969
      else
970
         analyse_from_flags( mat );
971
   }
972
 
973
   if (mat->inv && (mat->flags & MAT_DIRTY_INVERSE)) {
974
      matrix_invert( mat );
975
   }
976
 
977
   mat->flags &= ~(MAT_DIRTY_FLAGS|
978
                   MAT_DIRTY_TYPE|
979
                   MAT_DIRTY_INVERSE);
980
}
981
 
982
 
983
void
984
_math_matrix_copy( GLmatrix *to, const GLmatrix *from )
985
{
986
   MEMCPY( to->m, from->m, sizeof(Identity) );
987
   to->flags = from->flags;
988
   to->type = from->type;
989
 
990
   if (to->inv != 0) {
991
      if (from->inv == 0) {
992
         matrix_invert( to );
993
      }
994
      else {
995
         MEMCPY(to->inv, from->inv, sizeof(GLfloat)*16);
996
      }
997
   }
998
}
999
 
1000
 
1001
void
1002
_math_matrix_scale( GLmatrix *mat, GLfloat x, GLfloat y, GLfloat z )
1003
{
1004
   GLfloat *m = mat->m;
1005
   m[0] *= x;   m[4] *= y;   m[8]  *= z;
1006
   m[1] *= x;   m[5] *= y;   m[9]  *= z;
1007
   m[2] *= x;   m[6] *= y;   m[10] *= z;
1008
   m[3] *= x;   m[7] *= y;   m[11] *= z;
1009
 
1010
   if (fabs(x - y) < 1e-8 && fabs(x - z) < 1e-8)
1011
      mat->flags |= MAT_FLAG_UNIFORM_SCALE;
1012
   else
1013
      mat->flags |= MAT_FLAG_GENERAL_SCALE;
1014
 
1015
   mat->flags |= (MAT_DIRTY_TYPE |
1016
                  MAT_DIRTY_INVERSE);
1017
}
1018
 
1019
 
1020
void
1021
_math_matrix_translate( GLmatrix *mat, GLfloat x, GLfloat y, GLfloat z )
1022
{
1023
   GLfloat *m = mat->m;
1024
   m[12] = m[0] * x + m[4] * y + m[8]  * z + m[12];
1025
   m[13] = m[1] * x + m[5] * y + m[9]  * z + m[13];
1026
   m[14] = m[2] * x + m[6] * y + m[10] * z + m[14];
1027
   m[15] = m[3] * x + m[7] * y + m[11] * z + m[15];
1028
 
1029
   mat->flags |= (MAT_FLAG_TRANSLATION |
1030
                  MAT_DIRTY_TYPE |
1031
                  MAT_DIRTY_INVERSE);
1032
}
1033
 
1034
 
1035
void
1036
_math_matrix_loadf( GLmatrix *mat, const GLfloat *m )
1037
{
1038
   MEMCPY( mat->m, m, 16*sizeof(GLfloat) );
1039
   mat->flags = (MAT_FLAG_GENERAL | MAT_DIRTY);
1040
}
1041
 
1042
void
1043
_math_matrix_ctr( GLmatrix *m )
1044
{
1045
   m->m = (GLfloat *) ALIGN_MALLOC( 16 * sizeof(GLfloat), 16 );
1046
   if (m->m)
1047
      MEMCPY( m->m, Identity, sizeof(Identity) );
1048
   m->inv = NULL;
1049
   m->type = MATRIX_IDENTITY;
1050
   m->flags = 0;
1051
}
1052
 
1053
void
1054
_math_matrix_dtr( GLmatrix *m )
1055
{
1056
   if (m->m) {
1057
      ALIGN_FREE( m->m );
1058
      m->m = NULL;
1059
   }
1060
   if (m->inv) {
1061
      ALIGN_FREE( m->inv );
1062
      m->inv = NULL;
1063
   }
1064
}
1065
 
1066
 
1067
void
1068
_math_matrix_alloc_inv( GLmatrix *m )
1069
{
1070
   if (!m->inv) {
1071
      m->inv = (GLfloat *) ALIGN_MALLOC( 16 * sizeof(GLfloat), 16 );
1072
      if (m->inv)
1073
         MEMCPY( m->inv, Identity, 16 * sizeof(GLfloat) );
1074
   }
1075
}
1076
 
1077
 
1078
void
1079
_math_matrix_mul_matrix( GLmatrix *dest, const GLmatrix *a, const GLmatrix *b )
1080
{
1081
   dest->flags = (a->flags |
1082
                  b->flags |
1083
                  MAT_DIRTY_TYPE |
1084
                  MAT_DIRTY_INVERSE);
1085
 
1086
   if (TEST_MAT_FLAGS(dest, MAT_FLAGS_3D))
1087
      matmul34( dest->m, a->m, b->m );
1088
   else
1089
      matmul4( dest->m, a->m, b->m );
1090
}
1091
 
1092
 
1093
void
1094
_math_matrix_mul_floats( GLmatrix *dest, const GLfloat *m )
1095
{
1096
   dest->flags |= (MAT_FLAG_GENERAL |
1097
                   MAT_DIRTY_TYPE |
1098
                   MAT_DIRTY_INVERSE);
1099
 
1100
   matmul4( dest->m, dest->m, m );
1101
}
1102
 
1103
void
1104
_math_matrix_set_identity( GLmatrix *mat )
1105
{
1106
   MEMCPY( mat->m, Identity, 16*sizeof(GLfloat) );
1107
 
1108
   if (mat->inv)
1109
      MEMCPY( mat->inv, Identity, 16*sizeof(GLfloat) );
1110
 
1111
   mat->type = MATRIX_IDENTITY;
1112
   mat->flags &= ~(MAT_DIRTY_FLAGS|
1113
                   MAT_DIRTY_TYPE|
1114
                   MAT_DIRTY_INVERSE);
1115
}
1116
 
1117
 
1118
 
1119
void
1120
_math_transposef( GLfloat to[16], const GLfloat from[16] )
1121
{
1122
   to[0] = from[0];
1123
   to[1] = from[4];
1124
   to[2] = from[8];
1125
   to[3] = from[12];
1126
   to[4] = from[1];
1127
   to[5] = from[5];
1128
   to[6] = from[9];
1129
   to[7] = from[13];
1130
   to[8] = from[2];
1131
   to[9] = from[6];
1132
   to[10] = from[10];
1133
   to[11] = from[14];
1134
   to[12] = from[3];
1135
   to[13] = from[7];
1136
   to[14] = from[11];
1137
   to[15] = from[15];
1138
}
1139
 
1140
 
1141
void
1142
_math_transposed( GLdouble to[16], const GLdouble from[16] )
1143
{
1144
   to[0] = from[0];
1145
   to[1] = from[4];
1146
   to[2] = from[8];
1147
   to[3] = from[12];
1148
   to[4] = from[1];
1149
   to[5] = from[5];
1150
   to[6] = from[9];
1151
   to[7] = from[13];
1152
   to[8] = from[2];
1153
   to[9] = from[6];
1154
   to[10] = from[10];
1155
   to[11] = from[14];
1156
   to[12] = from[3];
1157
   to[13] = from[7];
1158
   to[14] = from[11];
1159
   to[15] = from[15];
1160
}
1161
 
1162
void
1163
_math_transposefd( GLfloat to[16], const GLdouble from[16] )
1164
{
1165
   to[0] = (GLfloat) from[0];
1166
   to[1] = (GLfloat) from[4];
1167
   to[2] = (GLfloat) from[8];
1168
   to[3] = (GLfloat) from[12];
1169
   to[4] = (GLfloat) from[1];
1170
   to[5] = (GLfloat) from[5];
1171
   to[6] = (GLfloat) from[9];
1172
   to[7] = (GLfloat) from[13];
1173
   to[8] = (GLfloat) from[2];
1174
   to[9] = (GLfloat) from[6];
1175
   to[10] = (GLfloat) from[10];
1176
   to[11] = (GLfloat) from[14];
1177
   to[12] = (GLfloat) from[3];
1178
   to[13] = (GLfloat) from[7];
1179
   to[14] = (GLfloat) from[11];
1180
   to[15] = (GLfloat) from[15];
1181
}