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
28struct 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
39float 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. */
55struct 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
92struct 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
108struct 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... */
118struct 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