Line data Source code
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 42 : struct vector matrix_eval(struct matrix mat, struct vector v)
29 : {
30 : struct vector result;
31 42 : 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 42 : 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 42 : 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 42 : return result;
35 : }
36 :
37 : //XXX: should probably pass by reference and we could
38 : //probably reuse this computation in matrix_invert
39 36 : float matrix_det(struct matrix mat)
40 : {
41 : float det;
42 108 : det = mat.m[0][0]*mat.m[1][1]*mat.m[2][2] +
43 72 : mat.m[0][1]*mat.m[1][2]*mat.m[2][0] +
44 72 : mat.m[0][2]*mat.m[1][0]*mat.m[2][1] -
45 72 : mat.m[0][0]*mat.m[1][2]*mat.m[2][1] -
46 36 : mat.m[0][1]*mat.m[1][0]*mat.m[2][2] -
47 36 : mat.m[0][2]*mat.m[1][1]*mat.m[2][0];
48 36 : 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 36 : 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 36 : float det = matrix_det(mat);
64 :
65 36 : if (det == 0) {
66 0 : dest_mat.invalid = true;
67 : } else {
68 36 : dest_mat.invalid = false;
69 : }
70 :
71 36 : det = 1/det;
72 :
73 144 : for (j = 0; j < 3; j++) {
74 432 : for (i = 0; i < 3; i++) {
75 : double p;
76 324 : int ai = a[i];
77 324 : int aj = a[j];
78 324 : int bi = b[i];
79 324 : int bj = b[j];
80 :
81 648 : p = mat.m[ai][aj] * mat.m[bi][bj] -
82 324 : mat.m[ai][bj] * mat.m[bi][aj];
83 324 : if (((i + j) & 1) != 0)
84 144 : p = -p;
85 :
86 324 : dest_mat.m[j][i] = det * p;
87 : }
88 : }
89 36 : return dest_mat;
90 : }
91 :
92 0 : struct matrix matrix_identity(void)
93 : {
94 : struct matrix i;
95 0 : i.m[0][0] = 1;
96 0 : i.m[0][1] = 0;
97 0 : i.m[0][2] = 0;
98 0 : i.m[1][0] = 0;
99 0 : i.m[1][1] = 1;
100 0 : i.m[1][2] = 0;
101 0 : i.m[2][0] = 0;
102 0 : i.m[2][1] = 0;
103 0 : i.m[2][2] = 1;
104 0 : i.invalid = false;
105 0 : return i;
106 : }
107 :
108 0 : struct matrix matrix_invalid(void)
109 : {
110 0 : struct matrix inv = matrix_identity();
111 0 : inv.invalid = true;
112 0 : return inv;
113 : }
114 :
115 :
116 : /* from pixman */
117 : /* MAT3per... */
118 50 : struct matrix matrix_multiply(struct matrix a, struct matrix b)
119 : {
120 : struct matrix result;
121 : int dx, dy;
122 : int o;
123 200 : for (dy = 0; dy < 3; dy++) {
124 600 : for (dx = 0; dx < 3; dx++) {
125 450 : double v = 0;
126 1800 : for (o = 0; o < 3; o++) {
127 1350 : v += a.m[dy][o] * b.m[o][dx];
128 : }
129 450 : result.m[dy][dx] = v;
130 : }
131 : }
132 50 : result.invalid = a.invalid || b.invalid;
133 50 : return result;
134 : }
135 :
136 :
|
