cogl/cogl-matrix.c

   1 /*
   2  * Cogl
   3  *
   4  * An object oriented GL/GLES Abstraction/Utility Layer
   5  *
   6  * Copyright (C) 2009,2010,2011 Intel Corporation.
   7  *
   8  * This library is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU Lesser General Public
  10  * License as published by the Free Software Foundation; either
  11  * version 2 of the License, or (at your option) any later version.
  12  *
  13  * This library is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16  * Lesser General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU Lesser General Public
  19  * License along with this library. If not, see <http://www.gnu.org/licenses/>.
  20  *
  21  * Authors:
  22  *   Robert Bragg <robert@linux.intel.com>
  23  */
  24 /*
  25  * Copyright (C) 1999-2005  Brian Paul   All Rights Reserved.
  26  *
  27  * Permission is hereby granted, free of charge, to any person obtaining a
  28  * copy of this software and associated documentation files (the "Software"),
  29  * to deal in the Software without restriction, including without limitation
  30  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  31  * and/or sell copies of the Software, and to permit persons to whom the
  32  * Software is furnished to do so, subject to the following conditions:
  33  *
  34  * The above copyright notice and this permission notice shall be included
  35  * in all copies or substantial portions of the Software.
  36  *
  37  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  38  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  39  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  40  * BRIAN PAUL BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN
  41  * AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  42  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  43  */
  44
  45 /*
  46  * Note: a lot of this code is based on code that was taken from Mesa.
  47  *
  48  * Changes compared to the original code from Mesa:
  49  *
  50  * - instead of allocating matrix->m and matrix->inv using malloc, our
  51  *   public CoglMatrix typedef is large enough to directly contain the
  52  *   matrix, its inverse, a type and a set of flags.
  53  * - instead of having a _cogl_matrix_analyse which updates the type,
  54  *   flags and inverse, we have _cogl_matrix_update_inverse which
  55  *   essentially does the same thing (internally making use of
  56  *   _cogl_matrix_update_type_and_flags()) but with additional guards in
  57  *   place to bail out when the inverse matrix is still valid.
  58  * - when initializing a matrix with the identity matrix we don't
  59  *   immediately initialize the inverse matrix; rather we just set the
  60  *   dirty flag for the inverse (since it's likely the user won't request
  61  *   the inverse of the identity matrix)
  62  */
  63
  64 #ifdef HAVE_CONFIG_H
  65 #include "config.h"
  66 #endif
  67
  68 #include <cogl-util.h>
  69 #include <cogl-debug.h>
  70 #include <cogl-quaternion.h>
  71 #include <cogl-quaternion-private.h>
  72 #include <cogl-matrix.h>
  73 #include <cogl-matrix-private.h>
  74 #include <cogl-quaternion-private.h>
  75
  76 #include <glib.h>
  77 #include <math.h>
  78 #include <string.h>
  79
  80 #ifdef _COGL_SUPPORTS_GTYPE_INTEGRATION
  81 #include <cogl-gtype-private.h>
  82 COGL_GTYPE_DEFINE_BOXED ("Matrix", matrix,
  83                          cogl_matrix_copy,
  84                          cogl_matrix_free);
  85 #endif
  86
  87 /*
  88  * Symbolic names to some of the entries in the matrix
  89  *
  90  * These are handy for the viewport mapping, which is expressed as a matrix.
  91  */
  92 #define MAT_SX 0
  93 #define MAT_SY 5
  94 #define MAT_SZ 10
  95 #define MAT_TX 12
  96 #define MAT_TY 13
  97 #define MAT_TZ 14
  98
  99 /*
 100  * These identify different kinds of 4x4 transformation matrices and we use
 101  * this information to find fast-paths when available.
 102  */
 103 enum CoglMatrixType {
 104    COGL_MATRIX_TYPE_GENERAL,    /**< general 4x4 matrix */
 105    COGL_MATRIX_TYPE_IDENTITY,   /**< identity matrix */
 106    COGL_MATRIX_TYPE_3D_NO_ROT,  /**< orthogonal projection and others... */
 107    COGL_MATRIX_TYPE_PERSPECTIVE,        /**< perspective projection matrix */
 108    COGL_MATRIX_TYPE_2D,         /**< 2-D transformation */
 109    COGL_MATRIX_TYPE_2D_NO_ROT,  /**< 2-D scale & translate only */
 110    COGL_MATRIX_TYPE_3D,         /**< 3-D transformation */
 111    COGL_MATRIX_N_TYPES
 112 } ;
 113
 114 #define DEG2RAD (G_PI/180.0)
 115
 116 /* Dot product of two 2-element vectors */
 117 #define DOT2(A,B)  ( (A)[0]*(B)[0] + (A)[1]*(B)[1] )
 118
 119 /* Dot product of two 3-element vectors */
 120 #define DOT3(A,B)  ( (A)[0]*(B)[0] + (A)[1]*(B)[1] + (A)[2]*(B)[2] )
 121
 122 #define CROSS3(N, U, V) \
 123 do { \
 124     (N)[0] = (U)[1]*(V)[2] - (U)[2]*(V)[1]; \
 125     (N)[1] = (U)[2]*(V)[0] - (U)[0]*(V)[2]; \
 126     (N)[2] = (U)[0]*(V)[1] - (U)[1]*(V)[0]; \
 127 } while (0)
 128
 129 #define SUB_3V(DST, SRCA, SRCB) \
 130 do { \
 131     (DST)[0] = (SRCA)[0] - (SRCB)[0]; \
 132     (DST)[1] = (SRCA)[1] - (SRCB)[1]; \
 133     (DST)[2] = (SRCA)[2] - (SRCB)[2]; \
 134 } while (0)
 135
 136 #define LEN_SQUARED_3FV( V ) ((V)[0]*(V)[0]+(V)[1]*(V)[1]+(V)[2]*(V)[2])
 137
 138 /*
 139  * \defgroup MatFlags MAT_FLAG_XXX-flags
 140  *
 141  * Bitmasks to indicate different kinds of 4x4 matrices in CoglMatrix::flags
 142  */
 143 #define MAT_FLAG_IDENTITY       0     /*< is an identity matrix flag.
 144                                        *   (Not actually used - the identity
 145                                        *   matrix is identified by the absense
 146                                        *   of all other flags.)
 147                                        */
 148 #define MAT_FLAG_GENERAL        0x1   /*< is a general matrix flag */
 149 #define MAT_FLAG_ROTATION       0x2   /*< is a rotation matrix flag */
 150 #define MAT_FLAG_TRANSLATION    0x4   /*< is a translation matrix flag */
 151 #define MAT_FLAG_UNIFORM_SCALE  0x8   /*< is an uniform scaling matrix flag */
 152 #define MAT_FLAG_GENERAL_SCALE  0x10  /*< is a general scaling matrix flag */
 153 #define MAT_FLAG_GENERAL_3D     0x20  /*< general 3D matrix flag */
 154 #define MAT_FLAG_PERSPECTIVE    0x40  /*< is a perspective proj matrix flag */
 155 #define MAT_FLAG_SINGULAR       0x80  /*< is a singular matrix flag */
 156 #define MAT_DIRTY_TYPE          0x100  /*< matrix type is dirty */
 157 #define MAT_DIRTY_FLAGS         0x200  /*< matrix flags are dirty */
 158 #define MAT_DIRTY_INVERSE       0x400  /*< matrix inverse is dirty */
 159
 160 /* angle preserving matrix flags mask */
 161 #define MAT_FLAGS_ANGLE_PRESERVING (MAT_FLAG_ROTATION | \
 162                                     MAT_FLAG_TRANSLATION | \
 163                                     MAT_FLAG_UNIFORM_SCALE)
 164
 165 /* geometry related matrix flags mask */
 166 #define MAT_FLAGS_GEOMETRY (MAT_FLAG_GENERAL | \
 167                             MAT_FLAG_ROTATION | \
 168                             MAT_FLAG_TRANSLATION | \
 169                             MAT_FLAG_UNIFORM_SCALE | \
 170                             MAT_FLAG_GENERAL_SCALE | \
 171                             MAT_FLAG_GENERAL_3D | \
 172                             MAT_FLAG_PERSPECTIVE | \
 173                             MAT_FLAG_SINGULAR)
 174
 175 /* length preserving matrix flags mask */
 176 #define MAT_FLAGS_LENGTH_PRESERVING (MAT_FLAG_ROTATION | \
 177                                      MAT_FLAG_TRANSLATION)
 178
 179
 180 /* 3D (non-perspective) matrix flags mask */
 181 #define MAT_FLAGS_3D (MAT_FLAG_ROTATION | \
 182                       MAT_FLAG_TRANSLATION | \
 183                       MAT_FLAG_UNIFORM_SCALE | \
 184                       MAT_FLAG_GENERAL_SCALE | \
 185                       MAT_FLAG_GENERAL_3D)
 186
 187 /* dirty matrix flags mask */
 188 #define MAT_DIRTY_ALL      (MAT_DIRTY_TYPE | \
 189                             MAT_DIRTY_FLAGS | \
 190                             MAT_DIRTY_INVERSE)
 191
 192
 193 /*
 194  * Test geometry related matrix flags.
 195  *
 196  * @mat a pointer to a CoglMatrix structure.
 197  * @a flags mask.
 198  *
 199  * Returns: non-zero if all geometry related matrix flags are contained within
 200  * the mask, or zero otherwise.
 201  */
 202 #define TEST_MAT_FLAGS(mat, a)  \
 203     ((MAT_FLAGS_GEOMETRY & (~(a)) & ((mat)->flags) ) == 0)
 204
 205
 206
 207 /*
 208  * Names of the corresponding CoglMatrixType values.
 209  */
 210 static const char *types[] = {
 211    "COGL_MATRIX_TYPE_GENERAL",
 212    "COGL_MATRIX_TYPE_IDENTITY",
 213    "COGL_MATRIX_TYPE_3D_NO_ROT",
 214    "COGL_MATRIX_TYPE_PERSPECTIVE",
 215    "COGL_MATRIX_TYPE_2D",
 216    "COGL_MATRIX_TYPE_2D_NO_ROT",
 217    "COGL_MATRIX_TYPE_3D"
 218 };
 219
 220
 221 /*
 222  * Identity matrix.
 223  */
 224 static float identity[16] = {
 225    1.0, 0.0, 0.0, 0.0,
 226    0.0, 1.0, 0.0, 0.0,
 227    0.0, 0.0, 1.0, 0.0,
 228    0.0, 0.0, 0.0, 1.0
 229 };
 230
 231
 232 #define A(row,col)  a[(col<<2)+row]
 233 #define B(row,col)  b[(col<<2)+row]
 234 #define R(row,col)  result[(col<<2)+row]
 235
 236 /*
 237  * Perform a full 4x4 matrix multiplication.
 238  *
 239  * <note>It's assumed that @result != @b. @product == @a is allowed.</note>
 240  *
 241  * <note>KW: 4*16 = 64 multiplications</note>
 242  */
 243 static void
 244 matrix_multiply4x4 (float *result, const float *a, const float *b)
 245 {
 246   int i;
 247   for (i = 0; i < 4; i++)
 248     {
 249       const float ai0 = A(i,0),  ai1=A(i,1),  ai2=A(i,2),  ai3=A(i,3);
 250       R(i,0) = ai0 * B(0,0) + ai1 * B(1,0) + ai2 * B(2,0) + ai3 * B(3,0);
 251       R(i,1) = ai0 * B(0,1) + ai1 * B(1,1) + ai2 * B(2,1) + ai3 * B(3,1);
 252       R(i,2) = ai0 * B(0,2) + ai1 * B(1,2) + ai2 * B(2,2) + ai3 * B(3,2);
 253       R(i,3) = ai0 * B(0,3) + ai1 * B(1,3) + ai2 * B(2,3) + ai3 * B(3,3);
 254     }
 255 }
 256
 257 /*
 258  * Multiply two matrices known to occupy only the top three rows, such
 259  * as typical model matrices, and orthogonal matrices.
 260  *
 261  * @a matrix.
 262  * @b matrix.
 263  * @product will receive the product of \p a and \p b.
 264  */
 265 static void
 266 matrix_multiply3x4 (float *result, const float *a, const float *b)
 267 {
 268   int i;
 269   for (i = 0; i < 3; i++)
 270     {
 271       const float ai0 = A(i,0), ai1 = A(i,1), ai2 = A(i,2), ai3 = A(i,3);
 272       R(i,0) = ai0 * B(0,0) + ai1 * B(1,0) + ai2 * B(2,0);
 273       R(i,1) = ai0 * B(0,1) + ai1 * B(1,1) + ai2 * B(2,1);
 274       R(i,2) = ai0 * B(0,2) + ai1 * B(1,2) + ai2 * B(2,2);
 275       R(i,3) = ai0 * B(0,3) + ai1 * B(1,3) + ai2 * B(2,3) + ai3;
 276     }
 277   R(3,0) = 0;
 278   R(3,1) = 0;
 279   R(3,2) = 0;
 280   R(3,3) = 1;
 281 }
 282
 283 #undef A
 284 #undef B
 285 #undef R
 286
 287 /*
 288  * Multiply a matrix by an array of floats with known properties.
 289  *
 290  * @mat pointer to a CoglMatrix structure containing the left multiplication
 291  * matrix, and that will receive the product result.
 292  * @m right multiplication matrix array.
 293  * @flags flags of the matrix \p m.
 294  *
 295  * Joins both flags and marks the type and inverse as dirty.  Calls
 296  * matrix_multiply3x4() if both matrices are 3D, or matrix_multiply4x4()
 297  * otherwise.
 298  */
 299 static void
 300 matrix_multiply_array_with_flags (CoglMatrix *result,
 301                                   const float *array,
 302                                   unsigned int flags)
 303 {
 304   result->flags |= (flags | MAT_DIRTY_TYPE | MAT_DIRTY_INVERSE);
 305
 306   if (TEST_MAT_FLAGS (result, MAT_FLAGS_3D))
 307     matrix_multiply3x4 ((float *)result, (float *)result, array);
 308   else
 309     matrix_multiply4x4 ((float *)result, (float *)result, array);
 310 }
 311
 312 /* Joins both flags and marks the type and inverse as dirty.  Calls
 313  * matrix_multiply3x4() if both matrices are 3D, or matrix_multiply4x4()
 314  * otherwise.
 315  */
 316 static void
 317 _cogl_matrix_multiply (CoglMatrix *result,
 318                        const CoglMatrix *a,
 319                        const CoglMatrix *b)
 320 {
 321   result->flags = (a->flags |
 322                    b->flags |
 323                    MAT_DIRTY_TYPE |
 324                    MAT_DIRTY_INVERSE);
 325
 326   if (TEST_MAT_FLAGS(result, MAT_FLAGS_3D))
 327     matrix_multiply3x4 ((float *)result, (float *)a, (float *)b);
 328   else
 329     matrix_multiply4x4 ((float *)result, (float *)a, (float *)b);
 330 }
 331
 332 void
 333 cogl_matrix_multiply (CoglMatrix *result,
 334                       const CoglMatrix *a,
 335                       const CoglMatrix *b)
 336 {
 337   _cogl_matrix_multiply (result, a, b);
 338   _COGL_MATRIX_DEBUG_PRINT (result);
 339 }
 340
 341 #if 0
 342 /* Marks the matrix flags with general flag, and type and inverse dirty flags.
 343  * Calls matrix_multiply4x4() for the multiplication.
 344  */
 345 static void
 346 _cogl_matrix_multiply_array (CoglMatrix *result, const float *array)
 347 {
 348   result->flags |= (MAT_FLAG_GENERAL |
 349                   MAT_DIRTY_TYPE |
 350                   MAT_DIRTY_INVERSE |
 351                   MAT_DIRTY_FLAGS);
 352
 353   matrix_multiply4x4 ((float *)result, (float *)result, (float *)array);
 354 }
 355 #endif
 356
 357 /*
 358  * Print a matrix array.
 359  *
 360  * Called by _cogl_matrix_print() to print a matrix or its inverse.
 361  */
 362 static void
 363 print_matrix_floats (const float m[16])
 364 {
 365   int i;
 366   for (i = 0;i < 4; i++)
 367     g_print ("\t%f %f %f %f\n", m[i], m[4+i], m[8+i], m[12+i] );
 368 }
 369
 370 /*
 371  * Dumps the contents of a CoglMatrix structure.
 372  */
 373 void
 374 _cogl_matrix_print (const CoglMatrix *matrix)
 375 {
 376   if (!(matrix->flags & MAT_DIRTY_TYPE))
 377     {
 378       _COGL_RETURN_IF_FAIL (matrix->type < COGL_MATRIX_N_TYPES);
 379       g_print ("Matrix type: %s, flags: %x\n",
 380                types[matrix->type], (int)matrix->flags);
 381     }
 382   else
 383     g_print ("Matrix type: DIRTY, flags: %x\n", (int)matrix->flags);
 384
 385   print_matrix_floats ((float *)matrix);
 386   g_print ("Inverse: \n");
 387   if (!(matrix->flags & MAT_DIRTY_INVERSE))
 388     {
 389       float prod[16];
 390       print_matrix_floats (matrix->inv);
 391       matrix_multiply4x4 (prod, (float *)matrix, matrix->inv);
 392       g_print ("Mat * Inverse:\n");
 393       print_matrix_floats (prod);
 394     }
 395   else
 396     g_print ("  - not available\n");
 397 }
 398
 399 /*
 400  * References an element of 4x4 matrix.
 401  *
 402  * @m matrix array.
 403  * @c column of the desired element.
 404  * @r row of the desired element.
 405  *
 406  * Returns: value of the desired element.
 407  *
 408  * Calculate the linear storage index of the element and references it.
 409  */
 410 #define MAT(m,r,c) (m)[(c)*4+(r)]
 411
 412 /*
 413  * Swaps the values of two floating pointer variables.
 414  *
 415  * Used by invert_matrix_general() to swap the row pointers.
 416  */
 417 #define SWAP_ROWS(a, b) { float *_tmp = a; (a)=(b); (b)=_tmp; }
 418
 419 /*
 420  * Compute inverse of 4x4 transformation matrix.
 421  *
 422  * @mat pointer to a CoglMatrix structure. The matrix inverse will be
 423  * stored in the CoglMatrix::inv attribute.
 424  *
 425  * Returns: %TRUE for success, %FALSE for failure (\p singular matrix).
 426  *
 427  * \author
 428  * Code contributed by Jacques Leroy jle@star.be
 429  *
 430  * Calculates the inverse matrix by performing the gaussian matrix reduction
 431  * with partial pivoting followed by back/substitution with the loops manually
 432  * unrolled.
 433  */
 434 static gboolean
 435 invert_matrix_general (CoglMatrix *matrix)
 436 {
 437   const float *m = (float *)matrix;
 438   float *out = matrix->inv;
 439   float wtmp[4][8];
 440   float m0, m1, m2, m3, s;
 441   float *r0, *r1, *r2, *r3;
 442
 443   r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];
 444
 445   r0[0] = MAT (m, 0, 0), r0[1] = MAT (m, 0, 1),
 446     r0[2] = MAT (m, 0, 2), r0[3] = MAT (m, 0, 3),
 447     r0[4] = 1.0, r0[5] = r0[6] = r0[7] = 0.0,
 448
 449     r1[0] = MAT (m, 1, 0), r1[1] = MAT (m, 1, 1),
 450     r1[2] = MAT (m, 1, 2), r1[3] = MAT (m, 1, 3),
 451     r1[5] = 1.0, r1[4] = r1[6] = r1[7] = 0.0,
 452
 453     r2[0] = MAT (m, 2, 0), r2[1] = MAT (m, 2, 1),
 454     r2[2] = MAT (m, 2, 2), r2[3] = MAT (m, 2, 3),
 455     r2[6] = 1.0, r2[4] = r2[5] = r2[7] = 0.0,
 456
 457     r3[0] = MAT (m, 3, 0), r3[1] = MAT (m, 3, 1),
 458     r3[2] = MAT (m, 3, 2), r3[3] = MAT (m, 3, 3),
 459     r3[7] = 1.0, r3[4] = r3[5] = r3[6] = 0.0;
 460
 461   /* choose pivot - or die */
 462   if (fabsf (r3[0]) > fabsf (r2[0]))
 463     SWAP_ROWS (r3, r2);
 464   if (fabsf (r2[0]) > fabsf (r1[0]))
 465     SWAP_ROWS (r2, r1);
 466   if (fabsf (r1[0]) > fabsf (r0[0]))
 467     SWAP_ROWS (r1, r0);
 468   if (0.0 == r0[0])
 469     return FALSE;
 470
 471   /* eliminate first variable     */
 472   m1 = r1[0]/r0[0]; m2 = r2[0]/r0[0]; m3 = r3[0]/r0[0];
 473   s = r0[1]; r1[1] -= m1 * s; r2[1] -= m2 * s; r3[1] -= m3 * s;
 474   s = r0[2]; r1[2] -= m1 * s; r2[2] -= m2 * s; r3[2] -= m3 * s;
 475   s = r0[3]; r1[3] -= m1 * s; r2[3] -= m2 * s; r3[3] -= m3 * s;
 476   s = r0[4];
 477   if (s != 0.0) { r1[4] -= m1 * s; r2[4] -= m2 * s; r3[4] -= m3 * s; }
 478   s = r0[5];
 479   if (s != 0.0) { r1[5] -= m1 * s; r2[5] -= m2 * s; r3[5] -= m3 * s; }
 480   s = r0[6];
 481   if (s != 0.0) { r1[6] -= m1 * s; r2[6] -= m2 * s; r3[6] -= m3 * s; }
 482   s = r0[7];
 483   if (s != 0.0) { r1[7] -= m1 * s; r2[7] -= m2 * s; r3[7] -= m3 * s; }
 484
 485   /* choose pivot - or die */
 486   if (fabsf (r3[1]) > fabsf (r2[1]))
 487     SWAP_ROWS (r3, r2);
 488   if (fabsf (r2[1]) > fabsf (r1[1]))
 489     SWAP_ROWS (r2, r1);
 490   if (0.0 == r1[1])
 491     return FALSE;
 492
 493   /* eliminate second variable */
 494   m2 = r2[1] / r1[1]; m3 = r3[1] / r1[1];
 495   r2[2] -= m2 * r1[2]; r3[2] -= m3 * r1[2];
 496   r2[3] -= m2 * r1[3]; r3[3] -= m3 * r1[3];
 497   s = r1[4]; if (0.0 != s) { r2[4] -= m2 * s; r3[4] -= m3 * s; }
 498   s = r1[5]; if (0.0 != s) { r2[5] -= m2 * s; r3[5] -= m3 * s; }
 499   s = r1[6]; if (0.0 != s) { r2[6] -= m2 * s; r3[6] -= m3 * s; }
 500   s = r1[7]; if (0.0 != s) { r2[7] -= m2 * s; r3[7] -= m3 * s; }
 501
 502   /* choose pivot - or die */
 503   if (fabsf (r3[2]) > fabsf (r2[2]))
 504     SWAP_ROWS (r3, r2);
 505   if (0.0 == r2[2])
 506     return FALSE;
 507
 508   /* eliminate third variable */
 509   m3 = r3[2] / r2[2];
 510   r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
 511     r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6],
 512     r3[7] -= m3 * r2[7];
 513
 514   /* last check */
 515   if (0.0 == r3[3])
 516     return FALSE;
 517
 518   s = 1.0f / r3[3];             /* now back substitute row 3 */
 519   r3[4] *= s; r3[5] *= s; r3[6] *= s; r3[7] *= s;
 520
 521   m2 = r2[3];                 /* now back substitute row 2 */
 522   s  = 1.0f / r2[2];
 523   r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
 524     r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
 525   m1 = r1[3];
 526   r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
 527     r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
 528   m0 = r0[3];
 529   r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
 530     r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;
 531
 532   m1 = r1[2];                 /* now back substitute row 1 */
 533   s  = 1.0f / r1[1];
 534   r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
 535     r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
 536   m0 = r0[2];
 537   r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
 538     r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;
 539
 540   m0 = r0[1];                 /* now back substitute row 0 */
 541   s  = 1.0f / r0[0];
 542   r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
 543     r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);
 544
 545   MAT (out, 0, 0) = r0[4]; MAT (out, 0, 1) = r0[5],
 546     MAT (out, 0, 2) = r0[6]; MAT (out, 0, 3) = r0[7],
 547     MAT (out, 1, 0) = r1[4]; MAT (out, 1, 1) = r1[5],
 548     MAT (out, 1, 2) = r1[6]; MAT (out, 1, 3) = r1[7],
 549     MAT (out, 2, 0) = r2[4]; MAT (out, 2, 1) = r2[5],
 550     MAT (out, 2, 2) = r2[6]; MAT (out, 2, 3) = r2[7],
 551     MAT (out, 3, 0) = r3[4]; MAT (out, 3, 1) = r3[5],
 552     MAT (out, 3, 2) = r3[6]; MAT (out, 3, 3) = r3[7];
 553
 554   return TRUE;
 555 }
 556 #undef SWAP_ROWS
 557
 558 /*
 559  * Compute inverse of a general 3d transformation matrix.
 560  *
 561  * @mat pointer to a CoglMatrix structure. The matrix inverse will be
 562  * stored in the CoglMatrix::inv attribute.
 563  *
 564  * Returns: %TRUE for success, %FALSE for failure (\p singular matrix).
 565  *
 566  * \author Adapted from graphics gems II.
 567  *
 568  * Calculates the inverse of the upper left by first calculating its
 569  * determinant and multiplying it to the symmetric adjust matrix of each
 570  * element. Finally deals with the translation part by transforming the
 571  * original translation vector using by the calculated submatrix inverse.
 572  */
 573 static gboolean
 574 invert_matrix_3d_general (CoglMatrix *matrix)
 575 {
 576   const float *in = (float *)matrix;
 577   float *out = matrix->inv;
 578   float pos, neg, t;
 579   float det;
 580
 581   /* Calculate the determinant of upper left 3x3 submatrix and
 582    * determine if the matrix is singular.
 583    */
 584   pos = neg = 0.0;
 585   t =  MAT (in,0,0) * MAT (in,1,1) * MAT (in,2,2);
 586   if (t >= 0.0) pos += t; else neg += t;
 587
 588   t =  MAT (in,1,0) * MAT (in,2,1) * MAT (in,0,2);
 589   if (t >= 0.0) pos += t; else neg += t;
 590
 591   t =  MAT (in,2,0) * MAT (in,0,1) * MAT (in,1,2);
 592   if (t >= 0.0) pos += t; else neg += t;
 593
 594   t = -MAT (in,2,0) * MAT (in,1,1) * MAT (in,0,2);
 595   if (t >= 0.0) pos += t; else neg += t;
 596
 597   t = -MAT (in,1,0) * MAT (in,0,1) * MAT (in,2,2);
 598   if (t >= 0.0) pos += t; else neg += t;
 599
 600   t = -MAT (in,0,0) * MAT (in,2,1) * MAT (in,1,2);
 601   if (t >= 0.0) pos += t; else neg += t;
 602
 603   det = pos + neg;
 604
 605   if (det*det < 1e-25)
 606     return FALSE;
 607
 608   det = 1.0f / det;
 609   MAT (out,0,0) =
 610     (  (MAT (in, 1, 1)*MAT (in, 2, 2) - MAT (in, 2, 1)*MAT (in, 1, 2) )*det);
 611   MAT (out,0,1) =
 612     (- (MAT (in, 0, 1)*MAT (in, 2, 2) - MAT (in, 2, 1)*MAT (in, 0, 2) )*det);
 613   MAT (out,0,2) =
 614     (  (MAT (in, 0, 1)*MAT (in, 1, 2) - MAT (in, 1, 1)*MAT (in, 0, 2) )*det);
 615   MAT (out,1,0) =
 616     (- (MAT (in,1,0)*MAT (in,2,2) - MAT (in,2,0)*MAT (in,1,2) )*det);
 617   MAT (out,1,1) =
 618     (  (MAT (in,0,0)*MAT (in,2,2) - MAT (in,2,0)*MAT (in,0,2) )*det);
 619   MAT (out,1,2) =
 620     (- (MAT (in,0,0)*MAT (in,1,2) - MAT (in,1,0)*MAT (in,0,2) )*det);
 621   MAT (out,2,0) =
 622     (  (MAT (in,1,0)*MAT (in,2,1) - MAT (in,2,0)*MAT (in,1,1) )*det);
 623   MAT (out,2,1) =
 624     (- (MAT (in,0,0)*MAT (in,2,1) - MAT (in,2,0)*MAT (in,0,1) )*det);
 625   MAT (out,2,2) =
 626     (  (MAT (in,0,0)*MAT (in,1,1) - MAT (in,1,0)*MAT (in,0,1) )*det);
 627
 628   /* Do the translation part */
 629   MAT (out,0,3) = - (MAT (in, 0, 3) * MAT (out, 0, 0) +
 630                     MAT (in, 1, 3) * MAT (out, 0, 1) +
 631                     MAT (in, 2, 3) * MAT (out, 0, 2) );
 632   MAT (out,1,3) = - (MAT (in, 0, 3) * MAT (out, 1, 0) +
 633                     MAT (in, 1, 3) * MAT (out, 1, 1) +
 634                     MAT (in, 2, 3) * MAT (out, 1, 2) );
 635   MAT (out,2,3) = - (MAT (in, 0, 3) * MAT (out, 2 ,0) +
 636                     MAT (in, 1, 3) * MAT (out, 2, 1) +
 637                     MAT (in, 2, 3) * MAT (out, 2, 2) );
 638
 639   return TRUE;
 640 }
 641
 642 /*
 643  * Compute inverse of a 3d transformation matrix.
 644  *
 645  * @mat pointer to a CoglMatrix structure. The matrix inverse will be
 646  * stored in the CoglMatrix::inv attribute.
 647  *
 648  * Returns: %TRUE for success, %FALSE for failure (\p singular matrix).
 649  *
 650  * If the matrix is not an angle preserving matrix then calls
 651  * invert_matrix_3d_general for the actual calculation. Otherwise calculates
 652  * the inverse matrix analyzing and inverting each of the scaling, rotation and
 653  * translation parts.
 654  */
 655 static gboolean
 656 invert_matrix_3d (CoglMatrix *matrix)
 657 {
 658   const float *in = (float *)matrix;
 659   float *out = matrix->inv;
 660
 661   if (!TEST_MAT_FLAGS(matrix, MAT_FLAGS_ANGLE_PRESERVING))
 662     return invert_matrix_3d_general (matrix);
 663
 664   if (matrix->flags & MAT_FLAG_UNIFORM_SCALE)
 665     {
 666       float scale = (MAT (in, 0, 0) * MAT (in, 0, 0) +
 667                      MAT (in, 0, 1) * MAT (in, 0, 1) +
 668                      MAT (in, 0, 2) * MAT (in, 0, 2));
 669
 670       if (scale == 0.0)
 671         return FALSE;
 672
 673       scale = 1.0f / scale;
 674
 675       /* Transpose and scale the 3 by 3 upper-left submatrix. */
 676       MAT (out, 0, 0) = scale * MAT (in, 0, 0);
 677       MAT (out, 1, 0) = scale * MAT (in, 0, 1);
 678       MAT (out, 2, 0) = scale * MAT (in, 0, 2);
 679       MAT (out, 0, 1) = scale * MAT (in, 1, 0);
 680       MAT (out, 1, 1) = scale * MAT (in, 1, 1);
 681       MAT (out, 2, 1) = scale * MAT (in, 1, 2);
 682       MAT (out, 0, 2) = scale * MAT (in, 2, 0);
 683       MAT (out, 1, 2) = scale * MAT (in, 2, 1);
 684       MAT (out, 2, 2) = scale * MAT (in, 2, 2);
 685     }
 686   else if (matrix->flags & MAT_FLAG_ROTATION)
 687     {
 688       /* Transpose the 3 by 3 upper-left submatrix. */
 689       MAT (out, 0, 0) = MAT (in, 0, 0);
 690       MAT (out, 1, 0) = MAT (in, 0, 1);
 691       MAT (out, 2, 0) = MAT (in, 0, 2);
 692       MAT (out, 0, 1) = MAT (in, 1, 0);
 693       MAT (out, 1, 1) = MAT (in, 1, 1);
 694       MAT (out, 2, 1) = MAT (in, 1, 2);
 695       MAT (out, 0, 2) = MAT (in, 2, 0);
 696       MAT (out, 1, 2) = MAT (in, 2, 1);
 697       MAT (out, 2, 2) = MAT (in, 2, 2);
 698     }
 699   else
 700     {
 701       /* pure translation */
 702       memcpy (out, identity, 16 * sizeof (float));
 703       MAT (out, 0, 3) = - MAT (in, 0, 3);
 704       MAT (out, 1, 3) = - MAT (in, 1, 3);
 705       MAT (out, 2, 3) = - MAT (in, 2, 3);
 706       return TRUE;
 707     }
 708
 709   if (matrix->flags & MAT_FLAG_TRANSLATION)
 710     {
 711       /* Do the translation part */
 712       MAT (out,0,3) = - (MAT (in, 0, 3) * MAT (out, 0, 0) +
 713                         MAT (in, 1, 3) * MAT (out, 0, 1) +
 714                         MAT (in, 2, 3) * MAT (out, 0, 2) );
 715       MAT (out,1,3) = - (MAT (in, 0, 3) * MAT (out, 1, 0) +
 716                         MAT (in, 1, 3) * MAT (out, 1, 1) +
 717                         MAT (in, 2, 3) * MAT (out, 1, 2) );
 718       MAT (out,2,3) = - (MAT (in, 0, 3) * MAT (out, 2, 0) +
 719                         MAT (in, 1, 3) * MAT (out, 2, 1) +
 720                         MAT (in, 2, 3) * MAT (out, 2, 2) );
 721     }
 722   else
 723     MAT (out, 0, 3) = MAT (out, 1, 3) = MAT (out, 2, 3) = 0.0;
 724
 725   return TRUE;
 726 }
 727
 728 /*
 729  * Compute inverse of an identity transformation matrix.
 730  *
 731  * @mat pointer to a CoglMatrix structure. The matrix inverse will be
 732  * stored in the CoglMatrix::inv attribute.
 733  *
 734  * Returns: always %TRUE.
 735  *
 736  * Simply copies identity into CoglMatrix::inv.
 737  */
 738 static gboolean
 739 invert_matrix_identity (CoglMatrix *matrix)
 740 {
 741   memcpy (matrix->inv, identity, 16 * sizeof (float));
 742   return TRUE;
 743 }
 744
 745 /*
 746  * Compute inverse of a no-rotation 3d transformation matrix.
 747  *
 748  * @mat pointer to a CoglMatrix structure. The matrix inverse will be
 749  * stored in the CoglMatrix::inv attribute.
 750  *
 751  * Returns: %TRUE for success, %FALSE for failure (\p singular matrix).
 752  *
 753  * Calculates the
 754  */
 755 static gboolean
 756 invert_matrix_3d_no_rotation (CoglMatrix *matrix)
 757 {
 758   const float *in = (float *)matrix;
 759   float *out = matrix->inv;
 760
 761   if (MAT (in,0,0) == 0 || MAT (in,1,1) == 0 || MAT (in,2,2) == 0)
 762     return FALSE;
 763
 764   memcpy (out, identity, 16 * sizeof (float));
 765   MAT (out,0,0) = 1.0f / MAT (in,0,0);
 766   MAT (out,1,1) = 1.0f / MAT (in,1,1);
 767   MAT (out,2,2) = 1.0f / MAT (in,2,2);
 768
 769   if (matrix->flags & MAT_FLAG_TRANSLATION)
 770     {
 771       MAT (out,0,3) = - (MAT (in,0,3) * MAT (out,0,0));
 772       MAT (out,1,3) = - (MAT (in,1,3) * MAT (out,1,1));
 773       MAT (out,2,3) = - (MAT (in,2,3) * MAT (out,2,2));
 774     }
 775
 776   return TRUE;
 777 }
 778
 779 /*
 780  * Compute inverse of a no-rotation 2d transformation matrix.
 781  *
 782  * @mat pointer to a CoglMatrix structure. The matrix inverse will be
 783  * stored in the CoglMatrix::inv attribute.
 784  *
 785  * Returns: %TRUE for success, %FALSE for failure (\p singular matrix).
 786  *
 787  * Calculates the inverse matrix by applying the inverse scaling and
 788  * translation to the identity matrix.
 789  */
 790 static gboolean
 791 invert_matrix_2d_no_rotation (CoglMatrix *matrix)
 792 {
 793   const float *in = (float *)matrix;
 794   float *out = matrix->inv;
 795
 796   if (MAT (in, 0, 0) == 0 || MAT (in, 1, 1) == 0)
 797     return FALSE;
 798
 799   memcpy (out, identity, 16 * sizeof (float));
 800   MAT (out, 0, 0) = 1.0f / MAT (in, 0, 0);
 801   MAT (out, 1, 1) = 1.0f / MAT (in, 1, 1);
 802
 803   if (matrix->flags & MAT_FLAG_TRANSLATION)
 804     {
 805       MAT (out, 0, 3) = - (MAT (in, 0, 3) * MAT (out, 0, 0));
 806       MAT (out, 1, 3) = - (MAT (in, 1, 3) * MAT (out, 1, 1));
 807     }
 808
 809   return TRUE;
 810 }
 811
 812 #if 0
 813 /* broken */
 814 static gboolean
 815 invert_matrix_perspective (CoglMatrix *matrix)
 816 {
 817   const float *in = matrix;
 818   float *out = matrix->inv;
 819
 820   if (MAT (in,2,3) == 0)
 821     return FALSE;
 822
 823   memcpy( out, identity, 16 * sizeof(float) );
 824
 825   MAT (out, 0, 0) = 1.0f / MAT (in, 0, 0);
 826   MAT (out, 1, 1) = 1.0f / MAT (in, 1, 1);
 827
 828   MAT (out, 0, 3) = MAT (in, 0, 2);
 829   MAT (out, 1, 3) = MAT (in, 1, 2);
 830
 831   MAT (out,2,2) = 0;
 832   MAT (out,2,3) = -1;
 833
 834   MAT (out,3,2) = 1.0f / MAT (in,2,3);
 835   MAT (out,3,3) = MAT (in,2,2) * MAT (out,3,2);
 836
 837   return TRUE;
 838 }
 839 #endif
 840
 841 /*
 842  * Matrix inversion function pointer type.
 843  */
 844 typedef gboolean (*inv_mat_func)(CoglMatrix *matrix);
 845
 846 /*
 847  * Table of the matrix inversion functions according to the matrix type.
 848  */
 849 static inv_mat_func inv_mat_tab[7] = {
 850     invert_matrix_general,
 851     invert_matrix_identity,
 852     invert_matrix_3d_no_rotation,
 853 #if 0
 854     /* Don't use this function for now - it fails when the projection matrix
 855      * is premultiplied by a translation (ala Chromium's tilesort SPU).
 856      */
 857     invert_matrix_perspective,
 858 #else
 859     invert_matrix_general,
 860 #endif
 861     invert_matrix_3d,           /* lazy! */
 862     invert_matrix_2d_no_rotation,
 863     invert_matrix_3d
 864 };
 865
 866 #define ZERO(x) (1<<x)
 867 #define ONE(x)  (1<<(x+16))
 868
 869 #define MASK_NO_TRX      (ZERO(12) | ZERO(13) | ZERO(14))
 870 #define MASK_NO_2D_SCALE ( ONE(0)  | ONE(5))
 871
 872 #define MASK_IDENTITY    ( ONE(0)  | ZERO(4)  | ZERO(8)  | ZERO(12) |\
 873                           ZERO(1)  |  ONE(5)  | ZERO(9)  | ZERO(13) |\
 874                           ZERO(2)  | ZERO(6)  |  ONE(10) | ZERO(14) |\
 875                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
 876
 877 #define MASK_2D_NO_ROT   (           ZERO(4)  | ZERO(8)  |           \
 878                           ZERO(1)  |            ZERO(9)  |           \
 879                           ZERO(2)  | ZERO(6)  |  ONE(10) | ZERO(14) |\
 880                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
 881
 882 #define MASK_2D          (                      ZERO(8)  |           \
 883                           ZERO(9)  |           \
 884                           ZERO(2)  | ZERO(6)  |  ONE(10) | ZERO(14) |\
 885                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
 886
 887
 888 #define MASK_3D_NO_ROT   (           ZERO(4)  | ZERO(8)  |           \
 889                           ZERO(1)  |            ZERO(9)  |           \
 890                           ZERO(2)  | ZERO(6)  |                      \
 891                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
 892
 893 #define MASK_3D          (                                           \
 894                           \
 895                           \
 896                           ZERO(3)  | ZERO(7)  | ZERO(11) |  ONE(15) )
 897
 898
 899 #define MASK_PERSPECTIVE (           ZERO(4)  |            ZERO(12) |\
 900                           ZERO(1)  |                       ZERO(13) |\
 901                           ZERO(2)  | ZERO(6)  |                      \
 902                           ZERO(3)  | ZERO(7)  |            ZERO(15) )
 903
 904 #define SQ(x) ((x)*(x))
 905
 906 /*
 907  * Determine type and flags from scratch.
 908  *
 909  * This is expensive enough to only want to do it once.
 910  */
 911 static void
 912 analyse_from_scratch (CoglMatrix *matrix)
 913 {
 914   const float *m = (float *)matrix;
 915   unsigned int mask = 0;
 916   unsigned int i;
 917
 918   for (i = 0 ; i < 16 ; i++)
 919     {
 920       if (m[i] == 0.0) mask |= (1<<i);
 921     }
 922
 923   if (m[0] == 1.0f) mask |= (1<<16);
 924   if (m[5] == 1.0f) mask |= (1<<21);
 925   if (m[10] == 1.0f) mask |= (1<<26);
 926   if (m[15] == 1.0f) mask |= (1<<31);
 927
 928   matrix->flags &= ~MAT_FLAGS_GEOMETRY;
 929
 930   /* Check for translation - no-one really cares
 931   */
 932   if ((mask & MASK_NO_TRX) != MASK_NO_TRX)
 933     matrix->flags |= MAT_FLAG_TRANSLATION;
 934
 935   /* Do the real work
 936   */
 937   if (mask == (unsigned int) MASK_IDENTITY)
 938     matrix->type = COGL_MATRIX_TYPE_IDENTITY;
 939   else if ((mask & MASK_2D_NO_ROT) == (unsigned int) MASK_2D_NO_ROT)
 940     {
 941       matrix->type = COGL_MATRIX_TYPE_2D_NO_ROT;
 942
 943       if ((mask & MASK_NO_2D_SCALE) != MASK_NO_2D_SCALE)
 944         matrix->flags |= MAT_FLAG_GENERAL_SCALE;
 945     }
 946   else if ((mask & MASK_2D) == (unsigned int) MASK_2D)
 947     {
 948       float mm = DOT2 (m, m);
 949       float m4m4 = DOT2 (m+4,m+4);
 950       float mm4 = DOT2 (m,m+4);
 951
 952       matrix->type = COGL_MATRIX_TYPE_2D;
 953
 954       /* Check for scale */
 955       if (SQ (mm-1) > SQ (1e-6) ||
 956           SQ (m4m4-1) > SQ (1e-6))
 957         matrix->flags |= MAT_FLAG_GENERAL_SCALE;
 958
 959       /* Check for rotation */
 960       if (SQ (mm4) > SQ (1e-6))
 961         matrix->flags |= MAT_FLAG_GENERAL_3D;
 962       else
 963         matrix->flags |= MAT_FLAG_ROTATION;
 964
 965     }
 966   else if ((mask & MASK_3D_NO_ROT) == (unsigned int) MASK_3D_NO_ROT)
 967     {
 968       matrix->type = COGL_MATRIX_TYPE_3D_NO_ROT;
 969
 970       /* Check for scale */
 971       if (SQ (m[0]-m[5]) < SQ (1e-6) &&
 972           SQ (m[0]-m[10]) < SQ (1e-6))
 973         {
 974           if (SQ (m[0]-1.0) > SQ (1e-6))
 975             matrix->flags |= MAT_FLAG_UNIFORM_SCALE;
 976         }
 977       else
 978         matrix->flags |= MAT_FLAG_GENERAL_SCALE;
 979     }
 980   else if ((mask & MASK_3D) == (unsigned int) MASK_3D)
 981     {
 982       float c1 = DOT3 (m,m);
 983       float c2 = DOT3 (m+4,m+4);
 984       float c3 = DOT3 (m+8,m+8);
 985       float d1 = DOT3 (m, m+4);
 986       float cp[3];
 987
 988       matrix->type = COGL_MATRIX_TYPE_3D;
 989
 990       /* Check for scale */
 991       if (SQ (c1-c2) < SQ (1e-6) && SQ (c1-c3) < SQ (1e-6))
 992         {
 993           if (SQ (c1-1.0) > SQ (1e-6))
 994             matrix->flags |= MAT_FLAG_UNIFORM_SCALE;
 995           /* else no scale at all */
 996         }
 997       else
 998         matrix->flags |= MAT_FLAG_GENERAL_SCALE;
 999
1000       /* Check for rotation */
1001       if (SQ (d1) < SQ (1e-6))
1002         {
1003           CROSS3 ( cp, m, m+4);
1004           SUB_3V ( cp, cp, (m+8));
1005           if (LEN_SQUARED_3FV(cp) < SQ(1e-6))
1006             matrix->flags |= MAT_FLAG_ROTATION;
1007           else
1008             matrix->flags |= MAT_FLAG_GENERAL_3D;
1009         }
1010       else
1011         matrix->flags |= MAT_FLAG_GENERAL_3D; /* shear, etc */
1012     }
1013   else if ((mask & MASK_PERSPECTIVE) == MASK_PERSPECTIVE && m[11]==-1.0f)
1014     {
1015       matrix->type = COGL_MATRIX_TYPE_PERSPECTIVE;
1016       matrix->flags |= MAT_FLAG_GENERAL;
1017     }
1018   else
1019     {
1020       matrix->type = COGL_MATRIX_TYPE_GENERAL;
1021       matrix->flags |= MAT_FLAG_GENERAL;
1022     }
1023 }
1024
1025 /*
1026  * Analyze a matrix given that its flags are accurate.
1027  *
1028  * This is the more common operation, hopefully.
1029  */
1030 static void
1031 analyse_from_flags (CoglMatrix *matrix)
1032 {
1033   const float *m = (float *)matrix;
1034
1035   if (TEST_MAT_FLAGS(matrix, 0))
1036     matrix->type = COGL_MATRIX_TYPE_IDENTITY;
1037   else if (TEST_MAT_FLAGS(matrix, (MAT_FLAG_TRANSLATION |
1038                                    MAT_FLAG_UNIFORM_SCALE |
1039                                    MAT_FLAG_GENERAL_SCALE)))
1040     {
1041       if ( m[10] == 1.0f && m[14] == 0.0f )
1042         matrix->type = COGL_MATRIX_TYPE_2D_NO_ROT;
1043       else
1044         matrix->type = COGL_MATRIX_TYPE_3D_NO_ROT;
1045     }
1046   else if (TEST_MAT_FLAGS (matrix, MAT_FLAGS_3D))
1047     {
1048       if (                               m[ 8]==0.0f
1049           &&                             m[ 9]==0.0f
1050           && m[2]==0.0f && m[6]==0.0f && m[10]==1.0f && m[14]==0.0f)
1051         {
1052           matrix->type = COGL_MATRIX_TYPE_2D;
1053         }
1054       else
1055         matrix->type = COGL_MATRIX_TYPE_3D;
1056     }
1057   else if (                 m[4]==0.0f                 && m[12]==0.0f
1058            && m[1]==0.0f                               && m[13]==0.0f
1059            && m[2]==0.0f && m[6]==0.0f
1060            && m[3]==0.0f && m[7]==0.0f && m[11]==-1.0f && m[15]==0.0f)
1061     {
1062       matrix->type = COGL_MATRIX_TYPE_PERSPECTIVE;
1063     }
1064   else
1065     matrix->type = COGL_MATRIX_TYPE_GENERAL;
1066 }
1067
1068 /*
1069  * Analyze and update the type and flags of a matrix.
1070  *
1071  * If the matrix type is dirty then calls either analyse_from_scratch() or
1072  * analyse_from_flags() to determine its type, according to whether the flags
1073  * are dirty or not, respectively. If the matrix has an inverse and it's dirty
1074  * then calls matrix_invert(). Finally clears the dirty flags.
1075  */
1076 static void
1077 _cogl_matrix_update_type_and_flags (CoglMatrix *matrix)
1078 {
1079   if (matrix->flags & MAT_DIRTY_TYPE)
1080     {
1081       if (matrix->flags & MAT_DIRTY_FLAGS)
1082         analyse_from_scratch (matrix);
1083       else
1084         analyse_from_flags (matrix);
1085     }
1086
1087   matrix->flags &= ~(MAT_DIRTY_FLAGS | MAT_DIRTY_TYPE);
1088 }
1089
1090 /*
1091  * Compute inverse of a transformation matrix.
1092  *
1093  * @mat pointer to a CoglMatrix structure. The matrix inverse will be
1094  * stored in the CoglMatrix::inv attribute.
1095  *
1096  * Returns: %TRUE for success, %FALSE for failure (\p singular matrix).
1097  *
1098  * Calls the matrix inversion function in inv_mat_tab corresponding to the
1099  * given matrix type.  In case of failure, updates the MAT_FLAG_SINGULAR flag,
1100  * and copies the identity matrix into CoglMatrix::inv.
1101  */
1102 static gboolean
1103 _cogl_matrix_update_inverse (CoglMatrix *matrix)
1104 {
1105   if (matrix->flags & MAT_DIRTY_FLAGS ||
1106       matrix->flags & MAT_DIRTY_INVERSE)
1107     {
1108       _cogl_matrix_update_type_and_flags (matrix);
1109
1110       if (inv_mat_tab[matrix->type](matrix))
1111         matrix->flags &= ~MAT_FLAG_SINGULAR;
1112       else
1113         {
1114           matrix->flags |= MAT_FLAG_SINGULAR;
1115           memcpy (matrix->inv, identity, 16 * sizeof (float));
1116         }
1117
1118       matrix->flags &= ~MAT_DIRTY_INVERSE;
1119     }
1120
1121   if (matrix->flags & MAT_FLAG_SINGULAR)
1122     return FALSE;
1123   else
1124     return TRUE;
1125 }
1126
1127 gboolean
1128 cogl_matrix_get_inverse (const CoglMatrix *matrix, CoglMatrix *inverse)
1129 {
1130   if (_cogl_matrix_update_inverse ((CoglMatrix *)matrix))
1131     {
1132       cogl_matrix_init_from_array (inverse, matrix->inv);
1133       return TRUE;
1134     }
1135   else
1136     {
1137       cogl_matrix_init_identity (inverse);
1138       return FALSE;
1139     }
1140 }
1141
1142 /*
1143  * Generate a 4x4 transformation matrix from glRotate parameters, and
1144  * post-multiply the input matrix by it.
1145  *
1146  * \author
1147  * This function was contributed by Erich Boleyn (erich@uruk.org).
1148  * Optimizations contributed by Rudolf Opalla (rudi@khm.de).
1149  */
1150 static void
1151 _cogl_matrix_rotate (CoglMatrix *matrix,
1152                      float angle,
1153                      float x,
1154                      float y,
1155                      float z)
1156 {
1157   float xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c, s, c;
1158   float m[16];
1159   gboolean optimized;
1160
1161   s = sinf (angle * DEG2RAD);
1162   c = cosf (angle * DEG2RAD);
1163
1164   memcpy (m, identity, 16 * sizeof (float));
1165   optimized = FALSE;
1166
1167 #define M(row,col)  m[col*4+row]
1168
1169   if (x == 0.0f)
1170     {
1171       if (y == 0.0f)
1172         {
1173           if (z != 0.0f)
1174             {
1175               optimized = TRUE;
1176               /* rotate only around z-axis */
1177               M (0,0) = c;
1178               M (1,1) = c;
1179               if (z < 0.0f)
1180                 {
1181                   M (0,1) = s;
1182                   M (1,0) = -s;
1183                 }
1184               else
1185                 {
1186                   M (0,1) = -s;
1187                   M (1,0) = s;
1188                 }
1189             }
1190         }
1191       else if (z == 0.0f)
1192         {
1193           optimized = TRUE;
1194           /* rotate only around y-axis */
1195           M (0,0) = c;
1196           M (2,2) = c;
1197           if (y < 0.0f)
1198             {
1199               M (0,2) = -s;
1200               M (2,0) = s;
1201             }
1202           else
1203             {
1204               M (0,2) = s;
1205               M (2,0) = -s;
1206             }
1207         }
1208     }
1209   else if (y == 0.0f)
1210     {
1211       if (z == 0.0f)
1212         {
1213           optimized = TRUE;
1214           /* rotate only around x-axis */
1215           M (1,1) = c;
1216           M (2,2) = c;
1217           if (x < 0.0f)
1218             {
1219               M (1,2) = s;
1220               M (2,1) = -s;
1221             }
1222           else
1223             {
1224               M (1,2) = -s;
1225               M (2,1) = s;
1226             }
1227         }
1228     }
1229
1230   if (!optimized)
1231     {
1232       const float mag = sqrtf (x * x + y * y + z * z);
1233
1234       if (mag <= 1.0e-4)
1235         {
1236           /* no rotation, leave mat as-is */
1237           return;
1238         }
1239
1240       x /= mag;
1241       y /= mag;
1242       z /= mag;
1243
1244
1245       /*
1246        *     Arbitrary axis rotation matrix.
1247        *
1248        *  This is composed of 5 matrices, Rz, Ry, T, Ry', Rz', multiplied
1249        *  like so:  Rz * Ry * T * Ry' * Rz'.  T is the final rotation
1250        *  (which is about the X-axis), and the two composite transforms
1251        *  Ry' * Rz' and Rz * Ry are (respectively) the rotations necessary
1252        *  from the arbitrary axis to the X-axis then back.  They are
1253        *  all elementary rotations.
1254        *
1255        *  Rz' is a rotation about the Z-axis, to bring the axis vector
1256        *  into the x-z plane.  Then Ry' is applied, rotating about the
1257        *  Y-axis to bring the axis vector parallel with the X-axis.  The
1258        *  rotation about the X-axis is then performed.  Ry and Rz are
1259        *  simply the respective inverse transforms to bring the arbitrary
1260        *  axis back to it's original orientation.  The first transforms
1261        *  Rz' and Ry' are considered inverses, since the data from the
1262        *  arbitrary axis gives you info on how to get to it, not how
1263        *  to get away from it, and an inverse must be applied.
1264        *
1265        *  The basic calculation used is to recognize that the arbitrary
1266        *  axis vector (x, y, z), since it is of unit length, actually
1267        *  represents the sines and cosines of the angles to rotate the
1268        *  X-axis to the same orientation, with theta being the angle about
1269        *  Z and phi the angle about Y (in the order described above)
1270        *  as follows:
1271        *
1272        *  cos ( theta ) = x / sqrt ( 1 - z^2 )
1273        *  sin ( theta ) = y / sqrt ( 1 - z^2 )
1274        *
1275        *  cos ( phi ) = sqrt ( 1 - z^2 )
1276        *  sin ( phi ) = z
1277        *
1278        *  Note that cos ( phi ) can further be inserted to the above
1279        *  formulas:
1280        *
1281        *  cos ( theta ) = x / cos ( phi )
1282        *  sin ( theta ) = y / sin ( phi )
1283        *
1284        *  ...etc.  Because of those relations and the standard trigonometric
1285        *  relations, it is pssible to reduce the transforms down to what
1286        *  is used below.  It may be that any primary axis chosen will give the
1287        *  same results (modulo a sign convention) using thie method.
1288        *
1289        *  Particularly nice is to notice that all divisions that might
1290        *  have caused trouble when parallel to certain planes or
1291        *  axis go away with care paid to reducing the expressions.
1292        *  After checking, it does perform correctly under all cases, since
1293        *  in all the cases of division where the denominator would have
1294        *  been zero, the numerator would have been zero as well, giving
1295        *  the expected result.
1296        */
1297
1298       xx = x * x;
1299       yy = y * y;
1300       zz = z * z;
1301       xy = x * y;
1302       yz = y * z;
1303       zx = z * x;
1304       xs = x * s;
1305       ys = y * s;
1306       zs = z * s;
1307       one_c = 1.0f - c;
1308
1309       /* We already hold the identity-matrix so we can skip some statements */
1310       M (0,0) = (one_c * xx) + c;
1311       M (0,1) = (one_c * xy) - zs;
1312       M (0,2) = (one_c * zx) + ys;
1313       /*    M (0,3) = 0.0f; */
1314
1315       M (1,0) = (one_c * xy) + zs;
1316       M (1,1) = (one_c * yy) + c;
1317       M (1,2) = (one_c * yz) - xs;
1318       /*    M (1,3) = 0.0f; */
1319
1320       M (2,0) = (one_c * zx) - ys;
1321       M (2,1) = (one_c * yz) + xs;
1322       M (2,2) = (one_c * zz) + c;
1323       /*    M (2,3) = 0.0f; */
1324
1325       /*
1326          M (3,0) = 0.0f;
1327          M (3,1) = 0.0f;
1328          M (3,2) = 0.0f;
1329          M (3,3) = 1.0f;
1330          */
1331     }
1332 #undef M
1333
1334   matrix_multiply_array_with_flags (matrix, m, MAT_FLAG_ROTATION);
1335 }
1336
1337 void
1338 cogl_matrix_rotate (CoglMatrix *matrix,
1339                     float angle,
1340                     float x,
1341                     float y,
1342                     float z)
1343 {
1344   _cogl_matrix_rotate (matrix, angle, x, y, z);
1345   _COGL_MATRIX_DEBUG_PRINT (matrix);
1346 }
1347
1348 /*
1349  * Apply a perspective projection matrix.
1350  *
1351  * Creates the projection matrix and multiplies it with matrix, marking the
1352  * MAT_FLAG_PERSPECTIVE flag.
1353  */
1354 static void
1355 _cogl_matrix_frustum (CoglMatrix *matrix,
1356                       float left,
1357                       float right,
1358                       float bottom,
1359                       float top,
1360                       float nearval,
1361                       float farval)
1362 {
1363   float x, y, a, b, c, d;
1364   float m[16];
1365
1366   x = (2.0f * nearval) / (right - left);
1367   y = (2.0f * nearval) / (top - bottom);
1368   a = (right + left) / (right - left);
1369   b = (top + bottom) / (top - bottom);
1370   c = -(farval + nearval) / ( farval - nearval);
1371   d = -(2.0f * farval * nearval) / (farval - nearval);  /* error? */
1372
1373 #define M(row,col)  m[col*4+row]
1374   M (0,0) = x;     M (0,1) = 0.0f;  M (0,2) = a;      M (0,3) = 0.0f;
1375   M (1,0) = 0.0f;  M (1,1) = y;     M (1,2) = b;      M (1,3) = 0.0f;
1376   M (2,0) = 0.0f;  M (2,1) = 0.0f;  M (2,2) = c;      M (2,3) = d;
1377   M (3,0) = 0.0f;  M (3,1) = 0.0f;  M (3,2) = -1.0f;  M (3,3) = 0.0f;
1378 #undef M
1379
1380   matrix_multiply_array_with_flags (matrix, m, MAT_FLAG_PERSPECTIVE);
1381 }
1382
1383 void
1384 cogl_matrix_frustum (CoglMatrix *matrix,
1385                      float       left,
1386                      float       right,
1387                      float       bottom,
1388                      float       top,
1389                      float       z_near,
1390                      float       z_far)
1391 {
1392   _cogl_matrix_frustum (matrix, left, right, bottom, top, z_near, z_far);
1393   _COGL_MATRIX_DEBUG_PRINT (matrix);
1394 }
1395
1396 void
1397 cogl_matrix_perspective (CoglMatrix *matrix,
1398                          float       fov_y,
1399                          float       aspect,
1400                          float       z_near,
1401                          float       z_far)
1402 {
1403   float ymax = z_near * tan (fov_y * G_PI / 360.0);
1404
1405   cogl_matrix_frustum (matrix,
1406                        -ymax * aspect,  /* left */
1407                        ymax * aspect,   /* right */
1408                        -ymax,           /* bottom */
1409                        ymax,            /* top */
1410                        z_near,
1411                        z_far);
1412   _COGL_MATRIX_DEBUG_PRINT (matrix);
1413 }
1414
1415 /*
1416  * Apply an orthographic projection matrix.
1417  *
1418  * Creates the projection matrix and multiplies it with matrix, marking the
1419  * MAT_FLAG_GENERAL_SCALE and MAT_FLAG_TRANSLATION flags.
1420  */
1421 static void
1422 _cogl_matrix_orthographic (CoglMatrix *matrix,
1423                            float x_1,
1424                            float y_1,
1425                            float x_2,
1426                            float y_2,
1427                            float nearval,
1428                            float farval)
1429 {
1430   float m[16];
1431
1432 #define M(row, col)  m[col * 4 + row]
1433   M (0,0) = 2.0f / (x_2 - x_1);
1434   M (0,1) = 0.0f;
1435   M (0,2) = 0.0f;
1436   M (0,3) = -(x_2 + x_1) / (x_2 - x_1);
1437
1438   M (1,0) = 0.0f;
1439   M (1,1) = 2.0f / (y_1 - y_2);
1440   M (1,2) = 0.0f;
1441   M (1,3) = -(y_1 + y_2) / (y_1 - y_2);
1442
1443   M (2,0) = 0.0f;
1444   M (2,1) = 0.0f;
1445   M (2,2) = -2.0f / (farval - nearval);
1446   M (2,3) = -(farval + nearval) / (farval - nearval);
1447
1448   M (3,0) = 0.0f;
1449   M (3,1) = 0.0f;
1450   M (3,2) = 0.0f;
1451   M (3,3) = 1.0f;
1452 #undef M
1453
1454   matrix_multiply_array_with_flags (matrix, m,
1455                                     (MAT_FLAG_GENERAL_SCALE |
1456                                      MAT_FLAG_TRANSLATION));
1457 }
1458
1459 void
1460 cogl_matrix_ortho (CoglMatrix *matrix,
1461                    float left,
1462                    float right,
1463                    float bottom,
1464                    float top,
1465                    float near,
1466                    float far)
1467 {
1468   _cogl_matrix_orthographic (matrix, left, top, right, bottom, near, far);
1469   _COGL_MATRIX_DEBUG_PRINT (matrix);
1470 }
1471
1472 void
1473 cogl_matrix_orthographic (CoglMatrix *matrix,
1474                           float x_1,
1475                           float y_1,
1476                           float x_2,
1477                           float y_2,
1478                           float near,
1479                           float far)
1480 {
1481   _cogl_matrix_orthographic (matrix, x_1, y_1, x_2, y_2, near, far);
1482   _COGL_MATRIX_DEBUG_PRINT (matrix);
1483 }
1484
1485 /*
1486  * Multiply a matrix with a general scaling matrix.
1487  *
1488  * Multiplies in-place the elements of matrix by the scale factors. Checks if
1489  * the scales factors are roughly the same, marking the MAT_FLAG_UNIFORM_SCALE
1490  * flag, or MAT_FLAG_GENERAL_SCALE. Marks the MAT_DIRTY_TYPE and
1491  * MAT_DIRTY_INVERSE dirty flags.
1492  */
1493 static void
1494 _cogl_matrix_scale (CoglMatrix *matrix, float x, float y, float z)
1495 {
1496   float *m = (float *)matrix;
1497   m[0] *= x;   m[4] *= y;   m[8]  *= z;
1498   m[1] *= x;   m[5] *= y;   m[9]  *= z;
1499   m[2] *= x;   m[6] *= y;   m[10] *= z;
1500   m[3] *= x;   m[7] *= y;   m[11] *= z;
1501
1502   if (fabsf (x - y) < 1e-8 && fabsf (x - z) < 1e-8)
1503     matrix->flags |= MAT_FLAG_UNIFORM_SCALE;
1504   else
1505     matrix->flags |= MAT_FLAG_GENERAL_SCALE;
1506
1507   matrix->flags |= (MAT_DIRTY_TYPE | MAT_DIRTY_INVERSE);
1508 }
1509
1510 void
1511 cogl_matrix_scale (CoglMatrix *matrix,
1512                    float sx,
1513                    float sy,
1514                    float sz)
1515 {
1516   _cogl_matrix_scale (matrix, sx, sy, sz);
1517   _COGL_MATRIX_DEBUG_PRINT (matrix);
1518 }
1519
1520 /*
1521  * Multiply a matrix with a translation matrix.
1522  *
1523  * Adds the translation coordinates to the elements of matrix in-place.  Marks
1524  * the MAT_FLAG_TRANSLATION flag, and the MAT_DIRTY_TYPE and MAT_DIRTY_INVERSE
1525  * dirty flags.
1526  */
1527 static void
1528 _cogl_matrix_translate (CoglMatrix *matrix, float x, float y, float z)
1529 {
1530   float *m = (float *)matrix;
1531   m[12] = m[0] * x + m[4] * y + m[8]  * z + m[12];
1532   m[13] = m[1] * x + m[5] * y + m[9]  * z + m[13];
1533   m[14] = m[2] * x + m[6] * y + m[10] * z + m[14];
1534   m[15] = m[3] * x + m[7] * y + m[11] * z + m[15];
1535
1536   matrix->flags |= (MAT_FLAG_TRANSLATION |
1537                     MAT_DIRTY_TYPE |
1538                     MAT_DIRTY_INVERSE);
1539 }
1540
1541 void
1542 cogl_matrix_translate (CoglMatrix *matrix,
1543                        float x,
1544                        float y,
1545                        float z)
1546 {
1547   _cogl_matrix_translate (matrix, x, y, z);
1548   _COGL_MATRIX_DEBUG_PRINT (matrix);
1549 }
1550
1551 #if 0
1552 /*
1553  * Set matrix to do viewport and depthrange mapping.
1554  * Transforms Normalized Device Coords to window/Z values.
1555  */
1556 static void
1557 _cogl_matrix_viewport (CoglMatrix *matrix,
1558                        float x, float y,
1559                        float width, float height,
1560                        float zNear, float zFar, float depthMax)
1561 {
1562   float *m = (float *)matrix;
1563   m[MAT_SX] = width / 2.0f;
1564   m[MAT_TX] = m[MAT_SX] + x;
1565   m[MAT_SY] = height / 2.0f;
1566   m[MAT_TY] = m[MAT_SY] + y;
1567   m[MAT_SZ] = depthMax * ((zFar - zNear) / 2.0f);
1568   m[MAT_TZ] = depthMax * ((zFar - zNear) / 2.0f + zNear);
1569   matrix->flags = MAT_FLAG_GENERAL_SCALE | MAT_FLAG_TRANSLATION;
1570   matrix->type = COGL_MATRIX_TYPE_3D_NO_ROT;
1571 }
1572 #endif
1573
1574 /*
1575  * Set a matrix to the identity matrix.
1576  *
1577  * @mat matrix.
1578  *
1579  * Copies ::identity into \p CoglMatrix::m, and into CoglMatrix::inv if
1580  * not NULL. Sets the matrix type to identity, resets the flags. It
1581  * doesn't initialize the inverse matrix, it just marks it dirty.
1582  */
1583 static void
1584 _cogl_matrix_init_identity (CoglMatrix *matrix)
1585 {
1586   memcpy (matrix, identity, 16 * sizeof (float));
1587
1588   matrix->type = COGL_MATRIX_TYPE_IDENTITY;
1589   matrix->flags = MAT_DIRTY_INVERSE;
1590 }
1591
1592 void
1593 cogl_matrix_init_identity (CoglMatrix *matrix)
1594 {
1595   _cogl_matrix_init_identity (matrix);
1596   _COGL_MATRIX_DEBUG_PRINT (matrix);
1597 }
1598
1599 #if 0
1600 /*
1601  * Test if the given matrix preserves vector lengths.
1602  */
1603 static gboolean
1604 _cogl_matrix_is_length_preserving (const CoglMatrix *m)
1605 {
1606   return TEST_MAT_FLAGS (m, MAT_FLAGS_LENGTH_PRESERVING);
1607 }
1608
1609 /*
1610  * Test if the given matrix does any rotation.
1611  * (or perhaps if the upper-left 3x3 is non-identity)
1612  */
1613 static gboolean
1614 _cogl_matrix_has_rotation (const CoglMatrix *matrix)
1615 {
1616   if (matrix->flags & (MAT_FLAG_GENERAL |
1617                        MAT_FLAG_ROTATION |
1618                        MAT_FLAG_GENERAL_3D |
1619                        MAT_FLAG_PERSPECTIVE))
1620     return TRUE;
1621   else
1622     return FALSE;
1623 }
1624
1625 static gboolean
1626 _cogl_matrix_is_general_scale (const CoglMatrix *matrix)
1627 {
1628   return (matrix->flags & MAT_FLAG_GENERAL_SCALE) ? TRUE : FALSE;
1629 }
1630
1631 static gboolean
1632 _cogl_matrix_is_dirty (const CoglMatrix *matrix)
1633 {
1634   return (matrix->flags & MAT_DIRTY_ALL) ? TRUE : FALSE;
1635 }
1636 #endif
1637
1638 /*
1639  * Loads a matrix array into CoglMatrix.
1640  *
1641  * @m matrix array.
1642  * @mat matrix.
1643  *
1644  * Copies \p m into CoglMatrix::m and marks the MAT_FLAG_GENERAL and
1645  * MAT_DIRTY_ALL
1646  * flags.
1647  */
1648 static void
1649 _cogl_matrix_init_from_array (CoglMatrix *matrix, const float *array)
1650 {
1651   memcpy (matrix, array, 16 * sizeof (float));
1652   matrix->flags = (MAT_FLAG_GENERAL | MAT_DIRTY_ALL);
1653 }
1654
1655 void
1656 cogl_matrix_init_from_array (CoglMatrix *matrix, const float *array)
1657 {
1658   _cogl_matrix_init_from_array (matrix, array);
1659   _COGL_MATRIX_DEBUG_PRINT (matrix);
1660 }
1661
1662 static void
1663 _cogl_matrix_init_from_quaternion (CoglMatrix *matrix,
1664                                    CoglQuaternion *quaternion)
1665 {
1666   float qnorm = _COGL_QUATERNION_NORM (quaternion);
1667   float s = (qnorm > 0.0f) ? (2.0f / qnorm) : 0.0f;
1668   float xs = quaternion->x * s;
1669   float ys = quaternion->y * s;
1670   float zs = quaternion->z * s;
1671   float wx = quaternion->w * xs;
1672   float wy = quaternion->w * ys;
1673   float wz = quaternion->w * zs;
1674   float xx = quaternion->x * xs;
1675   float xy = quaternion->x * ys;
1676   float xz = quaternion->x * zs;
1677   float yy = quaternion->y * ys;
1678   float yz = quaternion->y * zs;
1679   float zz = quaternion->z * zs;
1680
1681   matrix->xx = 1.0f - (yy + zz);
1682   matrix->yx = xy + wz;
1683   matrix->zx = xz - wy;
1684   matrix->xy = xy - wz;
1685   matrix->yy = 1.0f - (xx + zz);
1686   matrix->zy = yz + wx;
1687   matrix->xz = xz + wy;
1688   matrix->yz = yz - wx;
1689   matrix->zz = 1.0f - (xx + yy);
1690   matrix->xw = matrix->yw = matrix->zw = 0.0f;
1691   matrix->wx = matrix->wy = matrix->wz = 0.0f;
1692   matrix->ww = 1.0f;
1693
1694   matrix->flags = (MAT_FLAG_GENERAL | MAT_DIRTY_ALL);
1695 }
1696
1697 void
1698 cogl_matrix_init_from_quaternion (CoglMatrix *matrix,
1699                                   CoglQuaternion *quaternion)
1700 {
1701   _cogl_matrix_init_from_quaternion (matrix, quaternion);
1702 }
1703
1704 /*
1705  * Transpose a float matrix.
1706  */
1707 static void
1708 _cogl_matrix_util_transposef (float to[16], const float from[16])
1709 {
1710   to[0] = from[0];
1711   to[1] = from[4];
1712   to[2] = from[8];
1713   to[3] = from[12];
1714   to[4] = from[1];
1715   to[5] = from[5];
1716   to[6] = from[9];
1717   to[7] = from[13];
1718   to[8] = from[2];
1719   to[9] = from[6];
1720   to[10] = from[10];
1721   to[11] = from[14];
1722   to[12] = from[3];
1723   to[13] = from[7];
1724   to[14] = from[11];
1725   to[15] = from[15];
1726 }
1727
1728 void
1729 cogl_matrix_view_2d_in_frustum (CoglMatrix *matrix,
1730                                 float left,
1731                                 float right,
1732                                 float bottom,
1733                                 float top,
1734                                 float z_near,
1735                                 float z_2d,
1736                                 float width_2d,
1737                                 float height_2d)
1738 {
1739   float left_2d_plane = left / z_near * z_2d;
1740   float right_2d_plane = right / z_near * z_2d;
1741   float bottom_2d_plane = bottom / z_near * z_2d;
1742   float top_2d_plane = top / z_near * z_2d;
1743
1744   float width_2d_start = right_2d_plane - left_2d_plane;
1745   float height_2d_start = top_2d_plane - bottom_2d_plane;
1746
1747   /* Factors to scale from framebuffer geometry to frustum
1748    * cross-section geometry. */
1749   float width_scale = width_2d_start / width_2d;
1750   float height_scale = height_2d_start / height_2d;
1751
1752   cogl_matrix_translate (matrix,
1753                          left_2d_plane, top_2d_plane, -z_2d);
1754
1755   cogl_matrix_scale (matrix, width_scale, -height_scale, width_scale);
1756 }
1757
1758 /* Assuming a symmetric perspective matrix is being used for your
1759  * projective transform this convenience function lets you compose a
1760  * view transform such that geometry on the z=0 plane will map to
1761  * screen coordinates with a top left origin of (0,0) and with the
1762  * given width and height.
1763  */
1764 void
1765 cogl_matrix_view_2d_in_perspective (CoglMatrix *matrix,
1766                                     float fov_y,
1767                                     float aspect,
1768                                     float z_near,
1769                                     float z_2d,
1770                                     float width_2d,
1771                                     float height_2d)
1772 {
1773   float top = z_near * tan (fov_y * G_PI / 360.0);
1774   cogl_matrix_view_2d_in_frustum (matrix,
1775                                   -top * aspect,
1776                                   top * aspect,
1777                                   -top,
1778                                   top,
1779                                   z_near,
1780                                   z_2d,
1781                                   width_2d,
1782                                   height_2d);
1783 }
1784
1785 gboolean
1786 cogl_matrix_equal (gconstpointer v1, gconstpointer v2)
1787 {
1788   const CoglMatrix *a = v1;
1789   const CoglMatrix *b = v2;
1790
1791   _COGL_RETURN_VAL_IF_FAIL (v1 != NULL, FALSE);
1792   _COGL_RETURN_VAL_IF_FAIL (v2 != NULL, FALSE);
1793
1794   /* We want to avoid having a fuzzy _equal() function (e.g. that uses
1795    * an arbitrary epsilon value) since this function noteably conforms
1796    * to the prototype suitable for use with g_hash_table_new() and a
1797    * fuzzy hash function isn't really appropriate for comparing hash
1798    * table keys since it's possible that you could end up fetching
1799    * different values if you end up with multiple similar keys in use
1800    * at the same time. If you consider that fuzzyness allows cases
1801    * such as A == B == C but A != C then you could also end up loosing
1802    * values in a hash table.
1803    *
1804    * We do at least use the == operator to compare values though so
1805    * that -0 is considered equal to 0.
1806    */
1807
1808   /* XXX: We don't compare the flags, inverse matrix or padding */
1809   if (a->xx == b->xx &&
1810       a->xy == b->xy &&
1811       a->xz == b->xz &&
1812       a->xw == b->xw &&
1813       a->yx == b->yx &&
1814       a->yy == b->yy &&
1815       a->yz == b->yz &&
1816       a->yw == b->yw &&
1817       a->zx == b->zx &&
1818       a->zy == b->zy &&
1819       a->zz == b->zz &&
1820       a->zw == b->zw &&
1821       a->wx == b->wx &&
1822       a->wy == b->wy &&
1823       a->wz == b->wz &&
1824       a->ww == b->ww)
1825     return TRUE;
1826   else
1827     return FALSE;
1828 }
1829
1830 CoglMatrix *
1831 cogl_matrix_copy (const CoglMatrix *matrix)
1832 {
1833   if (G_LIKELY (matrix))
1834     return g_slice_dup (CoglMatrix, matrix);
1835
1836   return NULL;
1837 }
1838
1839 void
1840 cogl_matrix_free (CoglMatrix *matrix)
1841 {
1842   g_slice_free (CoglMatrix, matrix);
1843 }
1844
1845 const float *
1846 cogl_matrix_get_array (const CoglMatrix *matrix)
1847 {
1848   return (float *)matrix;
1849 }
1850
1851 void
1852 cogl_matrix_transform_point (const CoglMatrix *matrix,
1853                              float *x,
1854                              float *y,
1855                              float *z,
1856                              float *w)
1857 {
1858   float _x = *x, _y = *y, _z = *z, _w = *w;
1859
1860   *x = matrix->xx * _x + matrix->xy * _y + matrix->xz * _z + matrix->xw * _w;
1861   *y = matrix->yx * _x + matrix->yy * _y + matrix->yz * _z + matrix->yw * _w;
1862   *z = matrix->zx * _x + matrix->zy * _y + matrix->zz * _z + matrix->zw * _w;
1863   *w = matrix->wx * _x + matrix->wy * _y + matrix->wz * _z + matrix->ww * _w;
1864 }
1865
1866 typedef struct _Point2f
1867 {
1868   float x;
1869   float y;
1870 } Point2f;
1871
1872 typedef struct _Point3f
1873 {
1874   float x;
1875   float y;
1876   float z;
1877 } Point3f;
1878
1879 typedef struct _Point4f
1880 {
1881   float x;
1882   float y;
1883   float z;
1884   float w;
1885 } Point4f;
1886
1887 static void
1888 _cogl_matrix_transform_points_f2 (const CoglMatrix *matrix,
1889                                   size_t stride_in,
1890                                   const void *points_in,
1891                                   size_t stride_out,
1892                                   void *points_out,
1893                                   int n_points)
1894 {
1895   int i;
1896
1897   for (i = 0; i < n_points; i++)
1898     {
1899       Point2f p = *(Point2f *)((guint8 *)points_in + i * stride_in);
1900       Point3f *o = (Point3f *)((guint8 *)points_out + i * stride_out);
1901
1902       o->x = matrix->xx * p.x + matrix->xy * p.y + matrix->xw;
1903       o->y = matrix->yx * p.x + matrix->yy * p.y + matrix->yw;
1904       o->z = matrix->zx * p.x + matrix->zy * p.y + matrix->zw;
1905     }
1906 }
1907
1908 static void
1909 _cogl_matrix_project_points_f2 (const CoglMatrix *matrix,
1910                                 size_t stride_in,
1911                                 const void *points_in,
1912                                 size_t stride_out,
1913                                 void *points_out,
1914                                 int n_points)
1915 {
1916   int i;
1917
1918   for (i = 0; i < n_points; i++)
1919     {
1920       Point2f p = *(Point2f *)((guint8 *)points_in + i * stride_in);
1921       Point4f *o = (Point4f *)((guint8 *)points_out + i * stride_out);
1922
1923       o->x = matrix->xx * p.x + matrix->xy * p.y + matrix->xw;
1924       o->y = matrix->yx * p.x + matrix->yy * p.y + matrix->yw;
1925       o->z = matrix->zx * p.x + matrix->zy * p.y + matrix->zw;
1926       o->w = matrix->wx * p.x + matrix->wy * p.y + matrix->ww;
1927     }
1928 }
1929
1930 static void
1931 _cogl_matrix_transform_points_f3 (const CoglMatrix *matrix,
1932                                   size_t stride_in,
1933                                   const void *points_in,
1934                                   size_t stride_out,
1935                                   void *points_out,
1936                                   int n_points)
1937 {
1938   int i;
1939
1940   for (i = 0; i < n_points; i++)
1941     {
1942       Point3f p = *(Point3f *)((guint8 *)points_in + i * stride_in);
1943       Point3f *o = (Point3f *)((guint8 *)points_out + i * stride_out);
1944
1945       o->x = matrix->xx * p.x + matrix->xy * p.y +
1946              matrix->xz * p.z + matrix->xw;
1947       o->y = matrix->yx * p.x + matrix->yy * p.y +
1948              matrix->yz * p.z + matrix->yw;
1949       o->z = matrix->zx * p.x + matrix->zy * p.y +
1950              matrix->zz * p.z + matrix->zw;
1951     }
1952 }
1953
1954 static void
1955 _cogl_matrix_project_points_f3 (const CoglMatrix *matrix,
1956                                 size_t stride_in,
1957                                 const void *points_in,
1958                                 size_t stride_out,
1959                                 void *points_out,
1960                                 int n_points)
1961 {
1962   int i;
1963
1964   for (i = 0; i < n_points; i++)
1965     {
1966       Point3f p = *(Point3f *)((guint8 *)points_in + i * stride_in);
1967       Point4f *o = (Point4f *)((guint8 *)points_out + i * stride_out);
1968
1969       o->x = matrix->xx * p.x + matrix->xy * p.y +
1970              matrix->xz * p.z + matrix->xw;
1971       o->y = matrix->yx * p.x + matrix->yy * p.y +
1972              matrix->yz * p.z + matrix->yw;
1973       o->z = matrix->zx * p.x + matrix->zy * p.y +
1974              matrix->zz * p.z + matrix->zw;
1975       o->w = matrix->wx * p.x + matrix->wy * p.y +
1976              matrix->wz * p.z + matrix->ww;
1977     }
1978 }
1979
1980 static void
1981 _cogl_matrix_project_points_f4 (const CoglMatrix *matrix,
1982                                 size_t stride_in,
1983                                 const void *points_in,
1984                                 size_t stride_out,
1985                                 void *points_out,
1986                                 int n_points)
1987 {
1988   int i;
1989
1990   for (i = 0; i < n_points; i++)
1991     {
1992       Point4f p = *(Point4f *)((guint8 *)points_in + i * stride_in);
1993       Point4f *o = (Point4f *)((guint8 *)points_out + i * stride_out);
1994
1995       o->x = matrix->xx * p.x + matrix->xy * p.y +
1996              matrix->xz * p.z + matrix->xw * p.w;
1997       o->y = matrix->yx * p.x + matrix->yy * p.y +
1998              matrix->yz * p.z + matrix->yw * p.w;
1999       o->z = matrix->zx * p.x + matrix->zy * p.y +
2000              matrix->zz * p.z + matrix->zw * p.w;
2001       o->w = matrix->wx * p.x + matrix->wy * p.y +
2002              matrix->wz * p.z + matrix->ww * p.w;
2003     }
2004 }
2005
2006 void
2007 cogl_matrix_transform_points (const CoglMatrix *matrix,
2008                               int n_components,
2009                               size_t stride_in,
2010                               const void *points_in,
2011                               size_t stride_out,
2012                               void *points_out,
2013                               int n_points)
2014 {
2015   /* The results of transforming always have three components... */
2016   _COGL_RETURN_IF_FAIL (stride_out >= sizeof (Point3f));
2017
2018   if (n_components == 2)
2019     _cogl_matrix_transform_points_f2 (matrix,
2020                                       stride_in, points_in,
2021                                       stride_out, points_out,
2022                                       n_points);
2023   else
2024     {
2025       _COGL_RETURN_IF_FAIL (n_components == 3);
2026
2027       _cogl_matrix_transform_points_f3 (matrix,
2028                                         stride_in, points_in,
2029                                         stride_out, points_out,
2030                                         n_points);
2031     }
2032 }
2033
2034 void
2035 cogl_matrix_project_points (const CoglMatrix *matrix,
2036                             int n_components,
2037                             size_t stride_in,
2038                             const void *points_in,
2039                             size_t stride_out,
2040                             void *points_out,
2041                             int n_points)
2042 {
2043   if (n_components == 2)
2044     _cogl_matrix_project_points_f2 (matrix,
2045                                     stride_in, points_in,
2046                                     stride_out, points_out,
2047                                     n_points);
2048   else if (n_components == 3)
2049     _cogl_matrix_project_points_f3 (matrix,
2050                                     stride_in, points_in,
2051                                     stride_out, points_out,
2052                                     n_points);
2053   else
2054     {
2055       _COGL_RETURN_IF_FAIL (n_components == 4);
2056
2057       _cogl_matrix_project_points_f4 (matrix,
2058                                       stride_in, points_in,
2059                                       stride_out, points_out,
2060                                       n_points);
2061     }
2062 }
2063
2064 gboolean
2065 cogl_matrix_is_identity (const CoglMatrix *matrix)
2066 {
2067   if (!(matrix->flags & MAT_DIRTY_TYPE) &&
2068       matrix->type == COGL_MATRIX_TYPE_IDENTITY)
2069     return TRUE;
2070   else
2071     return memcmp (matrix, identity, sizeof (float) * 16) == 0;
2072 }
2073
2074 void
2075 cogl_matrix_look_at (CoglMatrix *matrix,
2076                      float eye_position_x,
2077                      float eye_position_y,
2078                      float eye_position_z,
2079                      float object_x,
2080                      float object_y,
2081                      float object_z,
2082                      float world_up_x,
2083                      float world_up_y,
2084                      float world_up_z)
2085 {
2086   CoglMatrix tmp;
2087   float forward[3];
2088   float side[3];
2089   float up[3];
2090
2091   /* Get a unit viewing direction vector */
2092   cogl_vector3_init (forward,
2093                      object_x - eye_position_x,
2094                      object_y - eye_position_y,
2095                      object_z - eye_position_z);
2096   cogl_vector3_normalize (forward);
2097
2098   cogl_vector3_init (up, world_up_x, world_up_y, world_up_z);
2099
2100   /* Take the sideways direction as being perpendicular to the viewing
2101    * direction and the word up vector. */
2102   cogl_vector3_cross_product (side, forward, up);
2103   cogl_vector3_normalize (side);
2104
2105   /* Now we have unit sideways and forward-direction vectors calculate
2106    * a new mutually perpendicular up vector. */
2107   cogl_vector3_cross_product (up, side, forward);
2108
2109   tmp.xx = side[0];
2110   tmp.yx = side[1];
2111   tmp.zx = side[2];
2112   tmp.wx = 0;
2113
2114   tmp.xy = up[0];
2115   tmp.yy = up[1];
2116   tmp.zy = up[2];
2117   tmp.wy = 0;
2118
2119   tmp.xz = -forward[0];
2120   tmp.yz = -forward[1];
2121   tmp.zz = -forward[2];
2122   tmp.wz = 0;
2123
2124   tmp.xw = 0;
2125   tmp.yw = 0;
2126   tmp.zw = 0;
2127   tmp.ww = 1;
2128
2129   tmp.flags = (MAT_FLAG_GENERAL_3D | MAT_DIRTY_TYPE | MAT_DIRTY_INVERSE);
2130
2131   cogl_matrix_translate (&tmp, -eye_position_x, -eye_position_y, -eye_position_z);
2132
2133   cogl_matrix_multiply (matrix, matrix, &tmp);
2134 }
2135
2136 void
2137 cogl_matrix_transpose (CoglMatrix *matrix)
2138 {
2139   float new_values[16];
2140
2141   /* We don't need to do anything if the matrix is the identity matrix */
2142   if (!(matrix->flags & MAT_DIRTY_TYPE) &&
2143       matrix->type == COGL_MATRIX_TYPE_IDENTITY)
2144     return;
2145
2146   _cogl_matrix_util_transposef (new_values, cogl_matrix_get_array (matrix));
2147
2148   cogl_matrix_init_from_array (matrix, new_values);
2149 }