1ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//---------------------------------------------------------------------------------
2ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//
3ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//  Little Color Management System
4ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//  Copyright (c) 1998-2012 Marti Maria Saguer
5ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//
6ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Permission is hereby granted, free of charge, to any person obtaining
7ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// a copy of this software and associated documentation files (the "Software"),
8ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// to deal in the Software without restriction, including without limitation
9ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// the rights to use, copy, modify, merge, publish, distribute, sublicense,
10ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// and/or sell copies of the Software, and to permit persons to whom the Software
11ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// is furnished to do so, subject to the following conditions:
12ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//
13ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// The above copyright notice and this permission notice shall be included in
14ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// all copies or substantial portions of the Software.
15ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//
16ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
18ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
20ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
21ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
22ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//
24ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//---------------------------------------------------------------------------------
25ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov//
26ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
27ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov#include "lcms2_internal.h"
28ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
29ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
30ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov#define DSWAP(x, y)     {cmsFloat64Number tmp = (x); (x)=(y); (y)=tmp;}
31ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
32ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
33ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Initiate a vector
34ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganovvoid CMSEXPORT _cmsVEC3init(cmsVEC3* r, cmsFloat64Number x, cmsFloat64Number y, cmsFloat64Number z)
35ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
36ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    r -> n[VX] = x;
37ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    r -> n[VY] = y;
38ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    r -> n[VZ] = z;
39ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
40ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
41ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Vector substraction
42ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganovvoid CMSEXPORT _cmsVEC3minus(cmsVEC3* r, const cmsVEC3* a, const cmsVEC3* b)
43ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
44ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov  r -> n[VX] = a -> n[VX] - b -> n[VX];
45ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov  r -> n[VY] = a -> n[VY] - b -> n[VY];
46ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov  r -> n[VZ] = a -> n[VZ] - b -> n[VZ];
47ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
48ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
49ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Vector cross product
50ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganovvoid CMSEXPORT _cmsVEC3cross(cmsVEC3* r, const cmsVEC3* u, const cmsVEC3* v)
51ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
52ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    r ->n[VX] = u->n[VY] * v->n[VZ] - v->n[VY] * u->n[VZ];
53ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    r ->n[VY] = u->n[VZ] * v->n[VX] - v->n[VZ] * u->n[VX];
54ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    r ->n[VZ] = u->n[VX] * v->n[VY] - v->n[VX] * u->n[VY];
55ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
56ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
57ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Vector dot product
58ee451cb395940862dad63c85adfe8f2fd55e864cSvet GanovcmsFloat64Number CMSEXPORT _cmsVEC3dot(const cmsVEC3* u, const cmsVEC3* v)
59ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
60ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    return u->n[VX] * v->n[VX] + u->n[VY] * v->n[VY] + u->n[VZ] * v->n[VZ];
61ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
62ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
63ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Euclidean length
64ee451cb395940862dad63c85adfe8f2fd55e864cSvet GanovcmsFloat64Number CMSEXPORT _cmsVEC3length(const cmsVEC3* a)
65ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
66ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    return sqrt(a ->n[VX] * a ->n[VX] +
67ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov                a ->n[VY] * a ->n[VY] +
68ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov                a ->n[VZ] * a ->n[VZ]);
69ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
70ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
71ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Euclidean distance
72ee451cb395940862dad63c85adfe8f2fd55e864cSvet GanovcmsFloat64Number CMSEXPORT _cmsVEC3distance(const cmsVEC3* a, const cmsVEC3* b)
73ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
74ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    cmsFloat64Number d1 = a ->n[VX] - b ->n[VX];
75ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    cmsFloat64Number d2 = a ->n[VY] - b ->n[VY];
76ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    cmsFloat64Number d3 = a ->n[VZ] - b ->n[VZ];
77ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
78ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    return sqrt(d1*d1 + d2*d2 + d3*d3);
79ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
80ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
81ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
82ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
83ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// 3x3 Identity
84ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganovvoid CMSEXPORT _cmsMAT3identity(cmsMAT3* a)
85ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
86ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    _cmsVEC3init(&a-> v[0], 1.0, 0.0, 0.0);
87ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    _cmsVEC3init(&a-> v[1], 0.0, 1.0, 0.0);
88ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    _cmsVEC3init(&a-> v[2], 0.0, 0.0, 1.0);
89ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
90ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
91ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganovstatic
92ee451cb395940862dad63c85adfe8f2fd55e864cSvet GanovcmsBool CloseEnough(cmsFloat64Number a, cmsFloat64Number b)
93ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
94ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    return fabs(b - a) < (1.0 / 65535.0);
95ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
96ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
97ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
98ee451cb395940862dad63c85adfe8f2fd55e864cSvet GanovcmsBool CMSEXPORT _cmsMAT3isIdentity(const cmsMAT3* a)
99ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
100ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    cmsMAT3 Identity;
101ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    int i, j;
102ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
103ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    _cmsMAT3identity(&Identity);
104ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
105ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    for (i=0; i < 3; i++)
106ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov        for (j=0; j < 3; j++)
107ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov            if (!CloseEnough(a ->v[i].n[j], Identity.v[i].n[j])) return FALSE;
108ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
109ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    return TRUE;
110ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
111ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
112ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
113ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Multiply two matrices
114ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganovvoid CMSEXPORT _cmsMAT3per(cmsMAT3* r, const cmsMAT3* a, const cmsMAT3* b)
115ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
116ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov#define ROWCOL(i, j) \
117ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    a->v[i].n[0]*b->v[0].n[j] + a->v[i].n[1]*b->v[1].n[j] + a->v[i].n[2]*b->v[2].n[j]
118ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
119ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    _cmsVEC3init(&r-> v[0], ROWCOL(0,0), ROWCOL(0,1), ROWCOL(0,2));
120ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    _cmsVEC3init(&r-> v[1], ROWCOL(1,0), ROWCOL(1,1), ROWCOL(1,2));
121ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    _cmsVEC3init(&r-> v[2], ROWCOL(2,0), ROWCOL(2,1), ROWCOL(2,2));
122ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
123ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov#undef ROWCOL //(i, j)
124ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
125ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
126ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
127ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
128ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Inverse of a matrix b = a^(-1)
129ee451cb395940862dad63c85adfe8f2fd55e864cSvet GanovcmsBool  CMSEXPORT _cmsMAT3inverse(const cmsMAT3* a, cmsMAT3* b)
130ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
131ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   cmsFloat64Number det, c0, c1, c2;
132ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
133ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   c0 =  a -> v[1].n[1]*a -> v[2].n[2] - a -> v[1].n[2]*a -> v[2].n[1];
134ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   c1 = -a -> v[1].n[0]*a -> v[2].n[2] + a -> v[1].n[2]*a -> v[2].n[0];
135ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   c2 =  a -> v[1].n[0]*a -> v[2].n[1] - a -> v[1].n[1]*a -> v[2].n[0];
136ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
137ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   det = a -> v[0].n[0]*c0 + a -> v[0].n[1]*c1 + a -> v[0].n[2]*c2;
138ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
139ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   if (fabs(det) < MATRIX_DET_TOLERANCE) return FALSE;  // singular matrix; can't invert
140ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
141ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   b -> v[0].n[0] = c0/det;
142ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   b -> v[0].n[1] = (a -> v[0].n[2]*a -> v[2].n[1] - a -> v[0].n[1]*a -> v[2].n[2])/det;
143ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   b -> v[0].n[2] = (a -> v[0].n[1]*a -> v[1].n[2] - a -> v[0].n[2]*a -> v[1].n[1])/det;
144ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   b -> v[1].n[0] = c1/det;
145ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   b -> v[1].n[1] = (a -> v[0].n[0]*a -> v[2].n[2] - a -> v[0].n[2]*a -> v[2].n[0])/det;
146ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   b -> v[1].n[2] = (a -> v[0].n[2]*a -> v[1].n[0] - a -> v[0].n[0]*a -> v[1].n[2])/det;
147ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   b -> v[2].n[0] = c2/det;
148ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   b -> v[2].n[1] = (a -> v[0].n[1]*a -> v[2].n[0] - a -> v[0].n[0]*a -> v[2].n[1])/det;
149ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   b -> v[2].n[2] = (a -> v[0].n[0]*a -> v[1].n[1] - a -> v[0].n[1]*a -> v[1].n[0])/det;
150ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
151ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov   return TRUE;
152ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
153ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
154ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
155ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Solve a system in the form Ax = b
156ee451cb395940862dad63c85adfe8f2fd55e864cSvet GanovcmsBool  CMSEXPORT _cmsMAT3solve(cmsVEC3* x, cmsMAT3* a, cmsVEC3* b)
157ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
158ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    cmsMAT3 m, a_1;
159ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
160ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    memmove(&m, a, sizeof(cmsMAT3));
161ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
162ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    if (!_cmsMAT3inverse(&m, &a_1)) return FALSE;  // Singular matrix
163ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
164ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    _cmsMAT3eval(x, &a_1, b);
165ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    return TRUE;
166ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
167ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
168ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov// Evaluate a vector across a matrix
169ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganovvoid CMSEXPORT _cmsMAT3eval(cmsVEC3* r, const cmsMAT3* a, const cmsVEC3* v)
170ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov{
171ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    r->n[VX] = a->v[0].n[VX]*v->n[VX] + a->v[0].n[VY]*v->n[VY] + a->v[0].n[VZ]*v->n[VZ];
172ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    r->n[VY] = a->v[1].n[VX]*v->n[VX] + a->v[1].n[VY]*v->n[VY] + a->v[1].n[VZ]*v->n[VZ];
173ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov    r->n[VZ] = a->v[2].n[VX]*v->n[VX] + a->v[2].n[VY]*v->n[VY] + a->v[2].n[VZ]*v->n[VZ];
174ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov}
175ee451cb395940862dad63c85adfe8f2fd55e864cSvet Ganov
176