• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
matrix_eval(struct matrix mat,struct vector v)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
matrix_det(struct matrix mat)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. */
matrix_invert(struct matrix mat)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 
matrix_identity(void)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 
matrix_invalid(void)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... */
matrix_multiply(struct matrix a,struct matrix b)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