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 } else { 68 dest_mat.invalid = false; 69 } 70 71 det = 1/det; 72 73 for (j = 0; j < 3; j++) { 74 for (i = 0; i < 3; i++) { 75 double p; 76 int ai = a[i]; 77 int aj = a[j]; 78 int bi = b[i]; 79 int bj = b[j]; 80 81 p = mat.m[ai][aj] * mat.m[bi][bj] - 82 mat.m[ai][bj] * mat.m[bi][aj]; 83 if (((i + j) & 1) != 0) 84 p = -p; 85 86 dest_mat.m[j][i] = det * p; 87 } 88 } 89 return dest_mat; 90 } 91 92 struct matrix matrix_identity(void) 93 { 94 struct matrix i; 95 i.m[0][0] = 1; 96 i.m[0][1] = 0; 97 i.m[0][2] = 0; 98 i.m[1][0] = 0; 99 i.m[1][1] = 1; 100 i.m[1][2] = 0; 101 i.m[2][0] = 0; 102 i.m[2][1] = 0; 103 i.m[2][2] = 1; 104 i.invalid = false; 105 return i; 106 } 107 108 struct matrix matrix_invalid(void) 109 { 110 struct matrix inv = matrix_identity(); 111 inv.invalid = true; 112 return inv; 113 } 114 115 116 /* from pixman */ 117 /* MAT3per... */ 118 struct matrix matrix_multiply(struct matrix a, struct matrix b) 119 { 120 struct matrix result; 121 int dx, dy; 122 int o; 123 for (dy = 0; dy < 3; dy++) { 124 for (dx = 0; dx < 3; dx++) { 125 double v = 0; 126 for (o = 0; o < 3; o++) { 127 v += a.m[dy][o] * b.m[o][dx]; 128 } 129 result.m[dy][dx] = v; 130 } 131 } 132 result.invalid = a.invalid || b.invalid; 133 return result; 134 } 135 136 137
