2 * Copyright (c) VMware, Inc.
4 * Permission is hereby granted, free of charge, to any person obtaining a
5 * copy of this software and associated documentation files (the "Software"),
6 * to deal in the Software without restriction, including without limitation
7 * on the rights to use, copy, modify, merge, publish, distribute, sub
8 * license, and/or sell copies of the Software, and to permit persons to whom
9 * the Software is furnished to do so, subject to the following conditions:
11 * The above copyright notice and this permission notice (including the next
12 * paragraph) shall be included in all copies or substantial portions of the
15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
18 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
20 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
21 * DEALINGS IN THE SOFTWARE.
28 #include "piglit-matrix.h"
31 #define DEG_TO_RAD(D) ((D) * M_PI / 180.0)
35 * Create a scaling matrix.
38 piglit_identity_matrix(float mat
[16])
63 * Create a scaling matrix.
66 piglit_scale_matrix(float mat
[16], float sx
, float sy
, float sz
)
91 * Create a translation matrix.
94 piglit_translation_matrix(float mat
[16], float tx
, float ty
, float tz
)
119 piglit_rotation_matrix(float mat
[16], float angle
, float x
, float y
, float z
)
121 /* Implementation borrowed from Mesa */
122 float xx
, yy
, zz
, xy
, yz
, zx
, xs
, ys
, zs
, one_c
, s
, c
;
124 bool optimized
= false;
126 s
= (float) sin(DEG_TO_RAD(angle
));
127 c
= (float) cos(DEG_TO_RAD(angle
));
129 piglit_identity_matrix(m
);
131 #define M(row,col) m[col*4+row]
137 /* rotate only around z-axis */
150 else if (z
== 0.0F
) {
152 /* rotate only around y-axis */
165 else if (y
== 0.0F
) {
168 /* rotate only around x-axis */
183 const float mag
= sqrtf(x
* x
+ y
* y
+ z
* z
);
186 /* no rotation, leave mat as-is */
196 * Arbitrary axis rotation matrix.
198 * This is composed of 5 matrices, Rz, Ry, T, Ry', Rz', multiplied
199 * like so: Rz * Ry * T * Ry' * Rz'. T is the final rotation
200 * (which is about the X-axis), and the two composite transforms
201 * Ry' * Rz' and Rz * Ry are (respectively) the rotations necessary
202 * from the arbitrary axis to the X-axis then back. They are
203 * all elementary rotations.
205 * Rz' is a rotation about the Z-axis, to bring the axis vector
206 * into the x-z plane. Then Ry' is applied, rotating about the
207 * Y-axis to bring the axis vector parallel with the X-axis. The
208 * rotation about the X-axis is then performed. Ry and Rz are
209 * simply the respective inverse transforms to bring the arbitrary
210 * axis back to its original orientation. The first transforms
211 * Rz' and Ry' are considered inverses, since the data from the
212 * arbitrary axis gives you info on how to get to it, not how
213 * to get away from it, and an inverse must be applied.
215 * The basic calculation used is to recognize that the arbitrary
216 * axis vector (x, y, z), since it is of unit length, actually
217 * represents the sines and cosines of the angles to rotate the
218 * X-axis to the same orientation, with theta being the angle about
219 * Z and phi the angle about Y (in the order described above)
222 * cos ( theta ) = x / sqrt ( 1 - z^2 )
223 * sin ( theta ) = y / sqrt ( 1 - z^2 )
225 * cos ( phi ) = sqrt ( 1 - z^2 )
228 * Note that cos ( phi ) can further be inserted to the above
231 * cos ( theta ) = x / cos ( phi )
232 * sin ( theta ) = y / sin ( phi )
234 * ...etc. Because of those relations and the standard trigonometric
235 * relations, it is pssible to reduce the transforms down to what
236 * is used below. It may be that any primary axis chosen will give the
237 * same results (modulo a sign convention) using thie method.
239 * Particularly nice is to notice that all divisions that might
240 * have caused trouble when parallel to certain planes or
241 * axis go away with care paid to reducing the expressions.
242 * After checking, it does perform correctly under all cases, since
243 * in all the cases of division where the denominator would have
244 * been zero, the numerator would have been zero as well, giving
245 * the expected result.
259 /* We already hold the identity-matrix so we can skip some statements */
260 M(0,0) = (one_c
* xx
) + c
;
261 M(0,1) = (one_c
* xy
) - zs
;
262 M(0,2) = (one_c
* zx
) + ys
;
265 M(1,0) = (one_c
* xy
) + zs
;
266 M(1,1) = (one_c
* yy
) + c
;
267 M(1,2) = (one_c
* yz
) - xs
;
270 M(2,0) = (one_c
* zx
) - ys
;
271 M(2,1) = (one_c
* yz
) + xs
;
272 M(2,2) = (one_c
* zz
) + c
;
287 piglit_ortho_matrix(float mat
[16],
288 float left
, float right
,
289 float bottom
, float top
,
290 float nearval
, float farval
)
292 #define M(row,col) mat[col*4+row]
293 M(0,0) = 2.0F
/ (right
-left
);
296 M(0,3) = -(right
+left
) / (right
-left
);
299 M(1,1) = 2.0F
/ (top
-bottom
);
301 M(1,3) = -(top
+bottom
) / (top
-bottom
);
305 M(2,2) = -2.0F
/ (farval
-nearval
);
306 M(2,3) = -(farval
+nearval
) / (farval
-nearval
);
317 piglit_frustum_matrix(float mat
[16],
318 float left
, float right
,
319 float bottom
, float top
,
320 float nearval
, float farval
)
322 float x
, y
, a
, b
, c
, d
;
324 x
= (2.0F
*nearval
) / (right
-left
);
325 y
= (2.0F
*nearval
) / (top
-bottom
);
326 a
= (right
+left
) / (right
-left
);
327 b
= (top
+bottom
) / (top
-bottom
);
328 c
= -(farval
+nearval
) / ( farval
-nearval
);
329 d
= -(2.0F
*farval
*nearval
) / (farval
-nearval
); /* error? */
331 #define M(row,col) mat[col*4+row]
332 M(0,0) = x
; M(0,1) = 0.0F
; M(0,2) = a
; M(0,3) = 0.0F
;
333 M(1,0) = 0.0F
; M(1,1) = y
; M(1,2) = b
; M(1,3) = 0.0F
;
334 M(2,0) = 0.0F
; M(2,1) = 0.0F
; M(2,2) = c
; M(2,3) = d
;
335 M(3,0) = 0.0F
; M(3,1) = 0.0F
; M(3,2) = -1.0F
; M(3,3) = 0.0F
;
341 piglit_matrix_mul_matrix(float product
[16],
342 const float a
[16], const float b
[16])
344 #define ELEM(MAT, ROW, COL) MAT[(COL) * 4 + (ROW)]
348 for (i
= 0; i
< 4; i
++) {
349 for (j
= 0; j
< 4; j
++) {
350 ELEM(tmp
, i
, j
) = (ELEM(a
, i
, 0) * ELEM(b
, 0, j
) +
351 ELEM(a
, i
, 1) * ELEM(b
, 1, j
) +
352 ELEM(a
, i
, 2) * ELEM(b
, 2, j
) +
353 ELEM(a
, i
, 3) * ELEM(b
, 3, j
));
357 for (i
= 0; i
< 16; i
++) {
365 * Compute "out = mat * in" where in and out are column vectors
366 * Typically used to transform homogeneous coordinates by a matrix.
369 piglit_matrix_mul_vector(float out
[4],
373 const float in0
= in
[0], in1
= in
[1], in2
= in
[2], in3
= in
[3];
374 #define M(row,col) mat[row + col*4]
375 out
[0] = M(0,0) * in0
+ M(0,1) * in1
+ M(0,2) * in2
+ M(0,3) * in3
;
376 out
[1] = M(1,0) * in0
+ M(1,1) * in1
+ M(1,2) * in2
+ M(1,3) * in3
;
377 out
[2] = M(2,0) * in0
+ M(2,1) * in1
+ M(2,2) * in2
+ M(2,3) * in3
;
378 out
[3] = M(3,0) * in0
+ M(3,1) * in1
+ M(3,2) * in2
+ M(3,3) * in3
;
384 * Transfrom NDC coordinate to window coordinate using a viewport.
387 piglit_ndc_to_window(float win
[3],
389 int vp_left
, int vp_bottom
, int vp_width
, int vp_height
)
391 float x
= ndc
[0] * 0.5 + 0.5;
392 float y
= ndc
[1] * 0.5 + 0.5;
393 float z
= ndc
[2] * 0.5 + 0.5;
394 win
[0] = vp_left
+ x
* vp_width
;
395 win
[1] = vp_bottom
+ y
* vp_height
;
401 * Transform an object coordinate to a window coordinate using a
402 * modelview matrix, projection matrix and viewport.
403 * \return true for success, false if coordinate is clipped away
406 piglit_project_to_window(float win
[3],
408 const float modelview
[16],
409 const float projection
[16],
410 int vp_left
, int vp_bottom
,
411 int vp_width
, int vp_height
)
413 float eye
[4], clip
[4], ndc
[4];
415 /* eye coord = modelview * object */
416 piglit_matrix_mul_vector(eye
, modelview
, obj
);
418 /* clip coord = projection * eye */
419 piglit_matrix_mul_vector(clip
, projection
, eye
);
421 /* view volume clipping */
422 if ( clip
[0] > clip
[3] ||
423 -clip
[0] > clip
[3] ||
425 -clip
[1] > clip
[3] ||
427 -clip
[2] > clip
[3]) {
432 /* ndc = clip / clip.w (divide by w) */
433 ndc
[0] = clip
[0] / clip
[3];
434 ndc
[1] = clip
[1] / clip
[3];
435 ndc
[2] = clip
[2] / clip
[3];
438 /* window = viewport_map(ndc) */
439 piglit_ndc_to_window(win
, ndc
, vp_left
, vp_bottom
, vp_width
, vp_height
);
446 piglit_print_matrix(const float mat
[16])
448 printf("%f %f %f %f\n", mat
[0], mat
[4], mat
[8], mat
[12]);
449 printf("%f %f %f %f\n", mat
[1], mat
[5], mat
[9], mat
[13]);
450 printf("%f %f %f %f\n", mat
[2], mat
[6], mat
[10], mat
[14]);
451 printf("%f %f %f %f\n", mat
[3], mat
[7], mat
[11], mat
[15]);
456 * Invert matrix using cramer's rule.
457 * This assumes that the matrix is non-singular (or non-near-singular).
460 piglit_matrix_inverse(float inv
[16], const float m
[16])
462 inv
[0] = m
[5] * m
[10] * m
[15] - m
[5] * m
[11] * m
[14] -
463 m
[9] * m
[6] * m
[15] + m
[9] * m
[7] * m
[14] +
464 m
[13] * m
[6] * m
[11] - m
[13] * m
[7] * m
[10];
466 inv
[4] = -m
[4] * m
[10] * m
[15] + m
[4] * m
[11] * m
[14] +
467 m
[8] * m
[6] * m
[15] - m
[8] * m
[7] * m
[14] -
468 m
[12] * m
[6] * m
[11] + m
[12] * m
[7] * m
[10];
470 inv
[8] = m
[4] * m
[9] * m
[15] - m
[4] * m
[11] * m
[13] -
471 m
[8] * m
[5] * m
[15] + m
[8] * m
[7] * m
[13] +
472 m
[12] * m
[5] * m
[11] - m
[12] * m
[7] * m
[9];
474 inv
[12] = -m
[4] * m
[9] * m
[14] + m
[4] * m
[10] * m
[13] +
475 m
[8] * m
[5] * m
[14] - m
[8] * m
[6] * m
[13] -
476 m
[12] * m
[5] * m
[10] + m
[12] * m
[6] * m
[9];
478 inv
[1] = -m
[1] * m
[10] * m
[15] + m
[1] * m
[11] * m
[14] +
479 m
[9] * m
[2] * m
[15] - m
[9] * m
[3] * m
[14] -
480 m
[13] * m
[2] * m
[11] + m
[13] * m
[3] * m
[10];
482 inv
[5] = m
[0] * m
[10] * m
[15] - m
[0] * m
[11] * m
[14] -
483 m
[8] * m
[2] * m
[15] + m
[8] * m
[3] * m
[14] +
484 m
[12] * m
[2] * m
[11] - m
[12] * m
[3] * m
[10];
486 inv
[9] = -m
[0] * m
[9] * m
[15] + m
[0] * m
[11] * m
[13] +
487 m
[8] * m
[1] * m
[15] - m
[8] * m
[3] * m
[13] -
488 m
[12] * m
[1] * m
[11] + m
[12] * m
[3] * m
[9];
490 inv
[13] = m
[0] * m
[9] * m
[14] - m
[0] * m
[10] * m
[13] -
491 m
[8] * m
[1] * m
[14] + m
[8] * m
[2] * m
[13] +
492 m
[12] * m
[1] * m
[10] - m
[12] * m
[2] * m
[9];
494 inv
[2] = m
[1] * m
[6] * m
[15] - m
[1] * m
[7] * m
[14] -
495 m
[5] * m
[2] * m
[15] + m
[5] * m
[3] * m
[14] +
496 m
[13] * m
[2] * m
[7] - m
[13] * m
[3] * m
[6];
498 inv
[6] = -m
[0] * m
[6] * m
[15] + m
[0] * m
[7] * m
[14] +
499 m
[4] * m
[2] * m
[15] - m
[4] * m
[3] * m
[14] -
500 m
[12] * m
[2] * m
[7] + m
[12] * m
[3] * m
[6];
502 inv
[10] = m
[0] * m
[5] * m
[15] - m
[0] * m
[7] * m
[13] -
503 m
[4] * m
[1] * m
[15] + m
[4] * m
[3] * m
[13] +
504 m
[12] * m
[1] * m
[7] - m
[12] * m
[3] * m
[5];
506 inv
[14] = -m
[0] * m
[5] * m
[14] + m
[0] * m
[6] * m
[13] +
507 m
[4] * m
[1] * m
[14] - m
[4] * m
[2] * m
[13] -
508 m
[12] * m
[1] * m
[6] + m
[12] * m
[2] * m
[5];
510 inv
[3] = -m
[1] * m
[6] * m
[11] + m
[1] * m
[7] * m
[10] +
511 m
[5] * m
[2] * m
[11] - m
[5] * m
[3] * m
[10] -
512 m
[9] * m
[2] * m
[7] + m
[9] * m
[3] * m
[6];
514 inv
[7] = m
[0] * m
[6] * m
[11] - m
[0] * m
[7] * m
[10] -
515 m
[4] * m
[2] * m
[11] + m
[4] * m
[3] * m
[10] +
516 m
[8] * m
[2] * m
[7] - m
[8] * m
[3] * m
[6];
518 inv
[11] = -m
[0] * m
[5] * m
[11] + m
[0] * m
[7] * m
[9] +
519 m
[4] * m
[1] * m
[11] - m
[4] * m
[3] * m
[9] -
520 m
[8] * m
[1] * m
[7] + m
[8] * m
[3] * m
[5];
522 inv
[15] = m
[0] * m
[5] * m
[10] - m
[0] * m
[6] * m
[9] -
523 m
[4] * m
[1] * m
[10] + m
[4] * m
[2] * m
[9] +
524 m
[8] * m
[1] * m
[6] - m
[8] * m
[2] * m
[5];
526 float det
= m
[0] * inv
[0] + m
[1] * inv
[4] + m
[2] * inv
[8] + m
[3] * inv
[12];
528 assert(fabsf(det
) > 1e-10f
);
530 for (int i
= 0; i
< 16; i
++)
531 inv
[i
] = inv
[i
] / det
;
536 piglit_matrix_transpose(float out
[16], const float m
[16])
538 for (int i
= 0; i
< 4; ++i
)
539 for (int j
= 0; j
< 4; ++j
)
540 out
[i
+4*j
] = m
[4*i
+j
];
546 * We could add a simple set of matrix stack functions.
547 * Could add scale/translate/rotate functions that accumulate onto
548 * the incoming matrix, similar to glScale, glTranslate, etc.