gfx/qcms/matrix.c

   1 /* vim: set ts=8 sw=8 noexpandtab: */
   2 //  qcms
   3 //  Copyright (C) 2009 Mozilla Foundation
   4 //  Copyright (C) 1998-2007 Marti Maria
   5 //
   6 // Permission is hereby granted, free of charge, to any person obtaining
   7 // a copy of this software and associated documentation files (the "Software"),
   8 // to deal in the Software without restriction, including without limitation
   9 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  10 // and/or sell copies of the Software, and to permit persons to whom the Software
  11 // is furnished to do so, subject to the following conditions:
  12 //
  13 // The above copyright notice and this permission notice shall be included in
  14 // all copies or substantial portions of the Software.
  15 //
  16 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
  18 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  20 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  21 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  22 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  23
  24 #include <stdlib.h>
  25 #include "qcmsint.h"
  26 #include "matrix.h"
  27
  28 struct vector matrix_eval(struct matrix mat, struct vector v)
  29 {
  30         struct vector result;
  31         result.v[0] = mat.m[0][0]*v.v[0] + mat.m[0][1]*v.v[1] + mat.m[0][2]*v.v[2];
  32         result.v[1] = mat.m[1][0]*v.v[0] + mat.m[1][1]*v.v[1] + mat.m[1][2]*v.v[2];
  33         result.v[2] = mat.m[2][0]*v.v[0] + mat.m[2][1]*v.v[1] + mat.m[2][2]*v.v[2];
  34         return result;
  35 }
  36
  37 //XXX: should probably pass by reference and we could
  38 //probably reuse this computation in matrix_invert
  39 float matrix_det(struct matrix mat)
  40 {
  41         float det;
  42         det = mat.m[0][0]*mat.m[1][1]*mat.m[2][2] +
  43                 mat.m[0][1]*mat.m[1][2]*mat.m[2][0] +
  44                 mat.m[0][2]*mat.m[1][0]*mat.m[2][1] -
  45                 mat.m[0][0]*mat.m[1][2]*mat.m[2][1] -
  46                 mat.m[0][1]*mat.m[1][0]*mat.m[2][2] -
  47                 mat.m[0][2]*mat.m[1][1]*mat.m[2][0];
  48         return det;
  49 }
  50
  51 /* from pixman and cairo and Mathematics for Game Programmers */
  52 /* lcms uses gauss-jordan elimination with partial pivoting which is
  53  * less efficient and not as numerically stable. See Mathematics for
  54  * Game Programmers. */
  55 struct matrix matrix_invert(struct matrix mat)
  56 {
  57         struct matrix dest_mat;
  58         int i,j;
  59         static int a[3] = { 2, 2, 1 };
  60         static int b[3] = { 1, 0, 0 };
  61
  62         /* inv  (A) = 1/det (A) * adj (A) */
  63         float det = matrix_det(mat);
  64
  65         if (det == 0) {
  66                 dest_mat.invalid = true;
  67                 return dest_mat;
  68         }
  69
  70         dest_mat.invalid = false;
  71
  72         det = 1/det;
  73
  74         for (j = 0; j < 3; j++) {
  75                 for (i = 0; i < 3; i++) {
  76                         double p;
  77                         int ai = a[i];
  78                         int aj = a[j];
  79                         int bi = b[i];
  80                         int bj = b[j];
  81
  82                         p = mat.m[ai][aj] * mat.m[bi][bj] -
  83                                 mat.m[ai][bj] * mat.m[bi][aj];
  84                         if (((i + j) & 1) != 0)
  85                                 p = -p;
  86
  87                         dest_mat.m[j][i] = det * p;
  88                 }
  89         }
  90         return dest_mat;
  91 }
  92
  93 struct matrix matrix_identity(void)
  94 {
  95         struct matrix i;
  96         i.m[0][0] = 1;
  97         i.m[0][1] = 0;
  98         i.m[0][2] = 0;
  99         i.m[1][0] = 0;
 100         i.m[1][1] = 1;
 101         i.m[1][2] = 0;
 102         i.m[2][0] = 0;
 103         i.m[2][1] = 0;
 104         i.m[2][2] = 1;
 105         i.invalid = false;
 106         return i;
 107 }
 108
 109 struct matrix matrix_invalid(void)
 110 {
 111         struct matrix inv = matrix_identity();
 112         inv.invalid = true;
 113         return inv;
 114 }
 115
 116
 117 /* from pixman */
 118 /* MAT3per... */
 119 struct matrix matrix_multiply(struct matrix a, struct matrix b)
 120 {
 121         struct matrix result;
 122         int dx, dy;
 123         int o;
 124         for (dy = 0; dy < 3; dy++) {
 125                 for (dx = 0; dx < 3; dx++) {
 126                         double v = 0;
 127                         for (o = 0; o < 3; o++) {
 128                                 v += a.m[dy][o] * b.m[o][dx];
 129                         }
 130                         result.m[dy][dx] = v;
 131                 }
 132         }
 133         result.invalid = a.invalid || b.invalid;
 134         return result;
 135 }
 136
 137