Home | History | Annotate | Download | only in src
      1 //---------------------------------------------------------------------------------
      2 //
      3 //  Little Color Management System
      4 //  Copyright (c) 1998-2012 Marti Maria Saguer
      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 //---------------------------------------------------------------------------------
     25 //
     26 
     27 #include "lcms2_internal.h"
     28 
     29 
     30 #define DSWAP(x, y)     {cmsFloat64Number tmp = (x); (x)=(y); (y)=tmp;}
     31 
     32 
     33 // Initiate a vector
     34 void CMSEXPORT _cmsVEC3init(cmsVEC3* r, cmsFloat64Number x, cmsFloat64Number y, cmsFloat64Number z)
     35 {
     36     r -> n[VX] = x;
     37     r -> n[VY] = y;
     38     r -> n[VZ] = z;
     39 }
     40 
     41 // Vector substraction
     42 void CMSEXPORT _cmsVEC3minus(cmsVEC3* r, const cmsVEC3* a, const cmsVEC3* b)
     43 {
     44   r -> n[VX] = a -> n[VX] - b -> n[VX];
     45   r -> n[VY] = a -> n[VY] - b -> n[VY];
     46   r -> n[VZ] = a -> n[VZ] - b -> n[VZ];
     47 }
     48 
     49 // Vector cross product
     50 void CMSEXPORT _cmsVEC3cross(cmsVEC3* r, const cmsVEC3* u, const cmsVEC3* v)
     51 {
     52     r ->n[VX] = u->n[VY] * v->n[VZ] - v->n[VY] * u->n[VZ];
     53     r ->n[VY] = u->n[VZ] * v->n[VX] - v->n[VZ] * u->n[VX];
     54     r ->n[VZ] = u->n[VX] * v->n[VY] - v->n[VX] * u->n[VY];
     55 }
     56 
     57 // Vector dot product
     58 cmsFloat64Number CMSEXPORT _cmsVEC3dot(const cmsVEC3* u, const cmsVEC3* v)
     59 {
     60     return u->n[VX] * v->n[VX] + u->n[VY] * v->n[VY] + u->n[VZ] * v->n[VZ];
     61 }
     62 
     63 // Euclidean length
     64 cmsFloat64Number CMSEXPORT _cmsVEC3length(const cmsVEC3* a)
     65 {
     66     return sqrt(a ->n[VX] * a ->n[VX] +
     67                 a ->n[VY] * a ->n[VY] +
     68                 a ->n[VZ] * a ->n[VZ]);
     69 }
     70 
     71 // Euclidean distance
     72 cmsFloat64Number CMSEXPORT _cmsVEC3distance(const cmsVEC3* a, const cmsVEC3* b)
     73 {
     74     cmsFloat64Number d1 = a ->n[VX] - b ->n[VX];
     75     cmsFloat64Number d2 = a ->n[VY] - b ->n[VY];
     76     cmsFloat64Number d3 = a ->n[VZ] - b ->n[VZ];
     77 
     78     return sqrt(d1*d1 + d2*d2 + d3*d3);
     79 }
     80 
     81 
     82 
     83 // 3x3 Identity
     84 void CMSEXPORT _cmsMAT3identity(cmsMAT3* a)
     85 {
     86     _cmsVEC3init(&a-> v[0], 1.0, 0.0, 0.0);
     87     _cmsVEC3init(&a-> v[1], 0.0, 1.0, 0.0);
     88     _cmsVEC3init(&a-> v[2], 0.0, 0.0, 1.0);
     89 }
     90 
     91 static
     92 cmsBool CloseEnough(cmsFloat64Number a, cmsFloat64Number b)
     93 {
     94     return fabs(b - a) < (1.0 / 65535.0);
     95 }
     96 
     97 
     98 cmsBool CMSEXPORT _cmsMAT3isIdentity(const cmsMAT3* a)
     99 {
    100     cmsMAT3 Identity;
    101     int i, j;
    102 
    103     _cmsMAT3identity(&Identity);
    104 
    105     for (i=0; i < 3; i++)
    106         for (j=0; j < 3; j++)
    107             if (!CloseEnough(a ->v[i].n[j], Identity.v[i].n[j])) return FALSE;
    108 
    109     return TRUE;
    110 }
    111 
    112 
    113 // Multiply two matrices
    114 void CMSEXPORT _cmsMAT3per(cmsMAT3* r, const cmsMAT3* a, const cmsMAT3* b)
    115 {
    116 #define ROWCOL(i, j) \
    117     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]
    118 
    119     _cmsVEC3init(&r-> v[0], ROWCOL(0,0), ROWCOL(0,1), ROWCOL(0,2));
    120     _cmsVEC3init(&r-> v[1], ROWCOL(1,0), ROWCOL(1,1), ROWCOL(1,2));
    121     _cmsVEC3init(&r-> v[2], ROWCOL(2,0), ROWCOL(2,1), ROWCOL(2,2));
    122 
    123 #undef ROWCOL //(i, j)
    124 }
    125 
    126 
    127 
    128 // Inverse of a matrix b = a^(-1)
    129 cmsBool  CMSEXPORT _cmsMAT3inverse(const cmsMAT3* a, cmsMAT3* b)
    130 {
    131    cmsFloat64Number det, c0, c1, c2;
    132 
    133    c0 =  a -> v[1].n[1]*a -> v[2].n[2] - a -> v[1].n[2]*a -> v[2].n[1];
    134    c1 = -a -> v[1].n[0]*a -> v[2].n[2] + a -> v[1].n[2]*a -> v[2].n[0];
    135    c2 =  a -> v[1].n[0]*a -> v[2].n[1] - a -> v[1].n[1]*a -> v[2].n[0];
    136 
    137    det = a -> v[0].n[0]*c0 + a -> v[0].n[1]*c1 + a -> v[0].n[2]*c2;
    138 
    139    if (fabs(det) < MATRIX_DET_TOLERANCE) return FALSE;  // singular matrix; can't invert
    140 
    141    b -> v[0].n[0] = c0/det;
    142    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;
    143    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;
    144    b -> v[1].n[0] = c1/det;
    145    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;
    146    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;
    147    b -> v[2].n[0] = c2/det;
    148    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;
    149    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;
    150 
    151    return TRUE;
    152 }
    153 
    154 
    155 // Solve a system in the form Ax = b
    156 cmsBool  CMSEXPORT _cmsMAT3solve(cmsVEC3* x, cmsMAT3* a, cmsVEC3* b)
    157 {
    158     cmsMAT3 m, a_1;
    159 
    160     memmove(&m, a, sizeof(cmsMAT3));
    161 
    162     if (!_cmsMAT3inverse(&m, &a_1)) return FALSE;  // Singular matrix
    163 
    164     _cmsMAT3eval(x, &a_1, b);
    165     return TRUE;
    166 }
    167 
    168 // Evaluate a vector across a matrix
    169 void CMSEXPORT _cmsMAT3eval(cmsVEC3* r, const cmsMAT3* a, const cmsVEC3* v)
    170 {
    171     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];
    172     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];
    173     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];
    174 }
    175 
    176