1 /*
   2  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
   3  *
   4  * This code is free software; you can redistribute it and/or modify it
   5  * under the terms of the GNU General Public License version 2 only, as
   6  * published by the Free Software Foundation.  Oracle designates this
   7  * particular file as subject to the "Classpath" exception as provided
   8  * by Oracle in the LICENSE file that accompanied this code.
   9  *
  10  * This code is distributed in the hope that it will be useful, but WITHOUT
  11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13  * version 2 for more details (a copy is included in the LICENSE file that
  14  * accompanied this code).
  15  *
  16  * You should have received a copy of the GNU General Public License version
  17  * 2 along with this work; if not, write to the Free Software Foundation,
  18  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  19  *
  20  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
  21  * or visit www.oracle.com if you need additional information or have any
  22  * questions.
  23  */
  24 
  25 // This file is available under and governed by the GNU General Public
  26 // License version 2 only, as published by the Free Software Foundation.
  27 // However, the following notice accompanied the original version of this
  28 // file:
  29 //
  30 //---------------------------------------------------------------------------------
  31 //
  32 //  Little Color Management System
  33 //  Copyright (c) 1998-2017 Marti Maria Saguer
  34 //
  35 // Permission is hereby granted, free of charge, to any person obtaining
  36 // a copy of this software and associated documentation files (the "Software"),
  37 // to deal in the Software without restriction, including without limitation
  38 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  39 // and/or sell copies of the Software, and to permit persons to whom the Software
  40 // is furnished to do so, subject to the following conditions:
  41 //
  42 // The above copyright notice and this permission notice shall be included in
  43 // all copies or substantial portions of the Software.
  44 //
  45 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  46 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
  47 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  48 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  49 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  50 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  51 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  52 //
  53 //---------------------------------------------------------------------------------
  54 //
  55 
  56 #include "lcms2_internal.h"
  57 
  58 
  59 #define DSWAP(x, y)     {cmsFloat64Number tmp = (x); (x)=(y); (y)=tmp;}
  60 
  61 
  62 // Initiate a vector
  63 void CMSEXPORT _cmsVEC3init(cmsVEC3* r, cmsFloat64Number x, cmsFloat64Number y, cmsFloat64Number z)
  64 {
  65     r -> n[VX] = x;
  66     r -> n[VY] = y;
  67     r -> n[VZ] = z;
  68 }
  69 
  70 // Vector subtraction
  71 void CMSEXPORT _cmsVEC3minus(cmsVEC3* r, const cmsVEC3* a, const cmsVEC3* b)
  72 {
  73   r -> n[VX] = a -> n[VX] - b -> n[VX];
  74   r -> n[VY] = a -> n[VY] - b -> n[VY];
  75   r -> n[VZ] = a -> n[VZ] - b -> n[VZ];
  76 }
  77 
  78 // Vector cross product
  79 void CMSEXPORT _cmsVEC3cross(cmsVEC3* r, const cmsVEC3* u, const cmsVEC3* v)
  80 {
  81     r ->n[VX] = u->n[VY] * v->n[VZ] - v->n[VY] * u->n[VZ];
  82     r ->n[VY] = u->n[VZ] * v->n[VX] - v->n[VZ] * u->n[VX];
  83     r ->n[VZ] = u->n[VX] * v->n[VY] - v->n[VX] * u->n[VY];
  84 }
  85 
  86 // Vector dot product
  87 cmsFloat64Number CMSEXPORT _cmsVEC3dot(const cmsVEC3* u, const cmsVEC3* v)
  88 {
  89     return u->n[VX] * v->n[VX] + u->n[VY] * v->n[VY] + u->n[VZ] * v->n[VZ];
  90 }
  91 
  92 // Euclidean length
  93 cmsFloat64Number CMSEXPORT _cmsVEC3length(const cmsVEC3* a)
  94 {
  95     return sqrt(a ->n[VX] * a ->n[VX] +
  96                 a ->n[VY] * a ->n[VY] +
  97                 a ->n[VZ] * a ->n[VZ]);
  98 }
  99 
 100 // Euclidean distance
 101 cmsFloat64Number CMSEXPORT _cmsVEC3distance(const cmsVEC3* a, const cmsVEC3* b)
 102 {
 103     cmsFloat64Number d1 = a ->n[VX] - b ->n[VX];
 104     cmsFloat64Number d2 = a ->n[VY] - b ->n[VY];
 105     cmsFloat64Number d3 = a ->n[VZ] - b ->n[VZ];
 106 
 107     return sqrt(d1*d1 + d2*d2 + d3*d3);
 108 }
 109 
 110 
 111 
 112 // 3x3 Identity
 113 void CMSEXPORT _cmsMAT3identity(cmsMAT3* a)
 114 {
 115     _cmsVEC3init(&a-> v[0], 1.0, 0.0, 0.0);
 116     _cmsVEC3init(&a-> v[1], 0.0, 1.0, 0.0);
 117     _cmsVEC3init(&a-> v[2], 0.0, 0.0, 1.0);
 118 }
 119 
 120 static
 121 cmsBool CloseEnough(cmsFloat64Number a, cmsFloat64Number b)
 122 {
 123     return fabs(b - a) < (1.0 / 65535.0);
 124 }
 125 
 126 
 127 cmsBool CMSEXPORT _cmsMAT3isIdentity(const cmsMAT3* a)
 128 {
 129     cmsMAT3 Identity;
 130     int i, j;
 131 
 132     _cmsMAT3identity(&Identity);
 133 
 134     for (i=0; i < 3; i++)
 135         for (j=0; j < 3; j++)
 136             if (!CloseEnough(a ->v[i].n[j], Identity.v[i].n[j])) return FALSE;
 137 
 138     return TRUE;
 139 }
 140 
 141 
 142 // Multiply two matrices
 143 void CMSEXPORT _cmsMAT3per(cmsMAT3* r, const cmsMAT3* a, const cmsMAT3* b)
 144 {
 145 #define ROWCOL(i, j) \
 146     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]
 147 
 148     _cmsVEC3init(&r-> v[0], ROWCOL(0,0), ROWCOL(0,1), ROWCOL(0,2));
 149     _cmsVEC3init(&r-> v[1], ROWCOL(1,0), ROWCOL(1,1), ROWCOL(1,2));
 150     _cmsVEC3init(&r-> v[2], ROWCOL(2,0), ROWCOL(2,1), ROWCOL(2,2));
 151 
 152 #undef ROWCOL //(i, j)
 153 }
 154 
 155 
 156 
 157 // Inverse of a matrix b = a^(-1)
 158 cmsBool  CMSEXPORT _cmsMAT3inverse(const cmsMAT3* a, cmsMAT3* b)
 159 {
 160    cmsFloat64Number det, c0, c1, c2;
 161 
 162    c0 =  a -> v[1].n[1]*a -> v[2].n[2] - a -> v[1].n[2]*a -> v[2].n[1];
 163    c1 = -a -> v[1].n[0]*a -> v[2].n[2] + a -> v[1].n[2]*a -> v[2].n[0];
 164    c2 =  a -> v[1].n[0]*a -> v[2].n[1] - a -> v[1].n[1]*a -> v[2].n[0];
 165 
 166    det = a -> v[0].n[0]*c0 + a -> v[0].n[1]*c1 + a -> v[0].n[2]*c2;
 167 
 168    if (fabs(det) < MATRIX_DET_TOLERANCE) return FALSE;  // singular matrix; can't invert
 169 
 170    b -> v[0].n[0] = c0/det;
 171    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;
 172    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;
 173    b -> v[1].n[0] = c1/det;
 174    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;
 175    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;
 176    b -> v[2].n[0] = c2/det;
 177    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;
 178    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;
 179 
 180    return TRUE;
 181 }
 182 
 183 
 184 // Solve a system in the form Ax = b
 185 cmsBool  CMSEXPORT _cmsMAT3solve(cmsVEC3* x, cmsMAT3* a, cmsVEC3* b)
 186 {
 187     cmsMAT3 m, a_1;
 188 
 189     memmove(&m, a, sizeof(cmsMAT3));
 190 
 191     if (!_cmsMAT3inverse(&m, &a_1)) return FALSE;  // Singular matrix
 192 
 193     _cmsMAT3eval(x, &a_1, b);
 194     return TRUE;
 195 }
 196 
 197 // Evaluate a vector across a matrix
 198 void CMSEXPORT _cmsMAT3eval(cmsVEC3* r, const cmsMAT3* a, const cmsVEC3* v)
 199 {
 200     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];
 201     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];
 202     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];
 203 }
 204 
 205