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