/*

 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.

 *

 * This code is free software; you can redistribute it and/or modify it

 * under the terms of the GNU General Public License version 2 only, as

 * published by the Free Software Foundation.  Oracle designates this

 * particular file as subject to the "Classpath" exception as provided

 * by Oracle in the LICENSE file that accompanied this code.

 *

 * This code is distributed in the hope that it will be useful, but WITHOUT

 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License

 * version 2 for more details (a copy is included in the LICENSE file that

 * accompanied this code).

 *

 * You should have received a copy of the GNU General Public License version

 * 2 along with this work; if not, write to the Free Software Foundation,

 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.

 *

 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA

 * or visit www.oracle.com if you need additional information or have any

 * questions.

 */

// This file is available under and governed by the GNU General Public

// License version 2 only, as published by the Free Software Foundation.

// However, the following notice accompanied the original version of this

// file:

//

//---------------------------------------------------------------------------------

//

//  Little Color Management System

//  Copyright (c) 1998-2016 Marti Maria Saguer

//

// Permission is hereby granted, free of charge, to any person obtaining

// a copy of this software and associated documentation files (the "Software"),

// to deal in the Software without restriction, including without limitation

// the rights to use, copy, modify, merge, publish, distribute, sublicense,

// and/or sell copies of the Software, and to permit persons to whom the Software

// is furnished to do so, subject to the following conditions:

//

// The above copyright notice and this permission notice shall be included in

// all copies or substantial portions of the Software.

//

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,

// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO

// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND

// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE

// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION

// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION

// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

//

//---------------------------------------------------------------------------------

//

#include "lcms2_internal.h"

//      inter PCS conversions XYZ <-> CIE L* a* b*

/*

       CIE 15:2004 CIELab is defined as:

       L* = 116*f(Y/Yn) - 16                     0 <= L* <= 100

       a* = 500*[f(X/Xn) - f(Y/Yn)]

       b* = 200*[f(Y/Yn) - f(Z/Zn)]

       and

              f(t) = t^(1/3)                     1 >= t >  (24/116)^3

                     (841/108)*t + (16/116)      0 <= t <= (24/116)^3

       Reverse transform is:

       X = Xn*[a* / 500 + (L* + 16) / 116] ^ 3   if (X/Xn) > (24/116)

         = Xn*(a* / 500 + L* / 116) / 7.787      if (X/Xn) <= (24/116)

       PCS in Lab2 is encoded as:

              8 bit Lab PCS:

                     L*      0..100 into a 0..ff byte.

                     a*      t + 128 range is -128.0  +127.0

                     b*

             16 bit Lab PCS:

                     L*     0..100  into a 0..ff00 word.

                     a*     t + 128  range is  -128.0  +127.9961

                     b*

Interchange Space   Component     Actual Range        Encoded Range

CIE XYZ             X             0 -> 1.99997        0x0000 -> 0xffff

CIE XYZ             Y             0 -> 1.99997        0x0000 -> 0xffff

CIE XYZ             Z             0 -> 1.99997        0x0000 -> 0xffff

Version 2,3

-----------

CIELAB (16 bit)     L*            0 -> 100.0          0x0000 -> 0xff00

CIELAB (16 bit)     a*            -128.0 -> +127.996  0x0000 -> 0x8000 -> 0xffff

CIELAB (16 bit)     b*            -128.0 -> +127.996  0x0000 -> 0x8000 -> 0xffff

Version 4

---------

CIELAB (16 bit)     L*            0 -> 100.0          0x0000 -> 0xffff

CIELAB (16 bit)     a*            -128.0 -> +127      0x0000 -> 0x8080 -> 0xffff

CIELAB (16 bit)     b*            -128.0 -> +127      0x0000 -> 0x8080 -> 0xffff

*/

// Conversions

void CMSEXPORT cmsXYZ2xyY(cmsCIExyY* Dest, const cmsCIEXYZ* Source)

{

    cmsFloat64Number ISum;

    ISum = 1./(Source -> X + Source -> Y + Source -> Z);

    Dest -> x = (Source -> X) * ISum;

    Dest -> y = (Source -> Y) * ISum;

    Dest -> Y = Source -> Y;

}

void CMSEXPORT cmsxyY2XYZ(cmsCIEXYZ* Dest, const cmsCIExyY* Source)

{

    Dest -> X = (Source -> x / Source -> y) * Source -> Y;

    Dest -> Y = Source -> Y;

    Dest -> Z = ((1 - Source -> x - Source -> y) / Source -> y) * Source -> Y;

}

/*

       The break point (24/116)^3 = (6/29)^3 is a very small amount of tristimulus

       primary (0.008856).  Generally, this only happens for

       nearly ideal blacks and for some orange / amber colors in transmission mode.

       For example, the Z value of the orange turn indicator lamp lens on an

       automobile will often be below this value.  But the Z does not

       contribute to the perceived color directly.

*/

static

cmsFloat64Number f(cmsFloat64Number t)

{

    const cmsFloat64Number Limit = (24.0/116.0) * (24.0/116.0) * (24.0/116.0);

    if (t <= Limit)

        return (841.0/108.0) * t + (16.0/116.0);

    else

        return pow(t, 1.0/3.0);

}

static

cmsFloat64Number f_1(cmsFloat64Number t)

{

    const cmsFloat64Number Limit = (24.0/116.0);

    if (t <= Limit) {

        return (108.0/841.0) * (t - (16.0/116.0));

    }

    return t * t * t;

}

// Standard XYZ to Lab. it can handle negative XZY numbers in some cases

void CMSEXPORT cmsXYZ2Lab(const cmsCIEXYZ* WhitePoint, cmsCIELab* Lab, const cmsCIEXYZ* xyz)

{

    cmsFloat64Number fx, fy, fz;

    if (WhitePoint == NULL)

        WhitePoint = cmsD50_XYZ();

    fx = f(xyz->X / WhitePoint->X);

    fy = f(xyz->Y / WhitePoint->Y);

    fz = f(xyz->Z / WhitePoint->Z);

    Lab->L = 116.0*fy - 16.0;

    Lab->a = 500.0*(fx - fy);

    Lab->b = 200.0*(fy - fz);

}

// Standard XYZ to Lab. It can return negative XYZ in some cases

void CMSEXPORT cmsLab2XYZ(const cmsCIEXYZ* WhitePoint, cmsCIEXYZ* xyz,  const cmsCIELab* Lab)

{

    cmsFloat64Number x, y, z;

    if (WhitePoint == NULL)

        WhitePoint = cmsD50_XYZ();

    y = (Lab-> L + 16.0) / 116.0;

    x = y + 0.002 * Lab -> a;

    z = y - 0.005 * Lab -> b;

    xyz -> X = f_1(x) * WhitePoint -> X;

    xyz -> Y = f_1(y) * WhitePoint -> Y;

    xyz -> Z = f_1(z) * WhitePoint -> Z;

}

static

cmsFloat64Number L2float2(cmsUInt16Number v)

{

    return (cmsFloat64Number) v / 652.800;

}

// the a/b part

static

cmsFloat64Number ab2float2(cmsUInt16Number v)

{

    return ((cmsFloat64Number) v / 256.0) - 128.0;

}

static

cmsUInt16Number L2Fix2(cmsFloat64Number L)

{

    return _cmsQuickSaturateWord(L *  652.8);

}

static

cmsUInt16Number ab2Fix2(cmsFloat64Number ab)

{

    return _cmsQuickSaturateWord((ab + 128.0) * 256.0);

}

static

cmsFloat64Number L2float4(cmsUInt16Number v)

{

    return (cmsFloat64Number) v / 655.35;

}

// the a/b part

static

cmsFloat64Number ab2float4(cmsUInt16Number v)

{

    return ((cmsFloat64Number) v / 257.0) - 128.0;

}

void CMSEXPORT cmsLabEncoded2FloatV2(cmsCIELab* Lab, const cmsUInt16Number wLab[3])

{

        Lab->L = L2float2(wLab[0]);

        Lab->a = ab2float2(wLab[1]);

        Lab->b = ab2float2(wLab[2]);

}

void CMSEXPORT cmsLabEncoded2Float(cmsCIELab* Lab, const cmsUInt16Number wLab[3])

{

        Lab->L = L2float4(wLab[0]);

        Lab->a = ab2float4(wLab[1]);

        Lab->b = ab2float4(wLab[2]);

}

static

cmsFloat64Number Clamp_L_doubleV2(cmsFloat64Number L)

{

    const cmsFloat64Number L_max = (cmsFloat64Number) (0xFFFF * 100.0) / 0xFF00;

    if (L < 0) L = 0;

    if (L > L_max) L = L_max;

    return L;

}

static

cmsFloat64Number Clamp_ab_doubleV2(cmsFloat64Number ab)

{

    if (ab < MIN_ENCODEABLE_ab2) ab = MIN_ENCODEABLE_ab2;

    if (ab > MAX_ENCODEABLE_ab2) ab = MAX_ENCODEABLE_ab2;

    return ab;

}

void CMSEXPORT cmsFloat2LabEncodedV2(cmsUInt16Number wLab[3], const cmsCIELab* fLab)

{

    cmsCIELab Lab;

    Lab.L = Clamp_L_doubleV2(fLab ->L);

    Lab.a = Clamp_ab_doubleV2(fLab ->a);

    Lab.b = Clamp_ab_doubleV2(fLab ->b);

    wLab[0] = L2Fix2(Lab.L);

    wLab[1] = ab2Fix2(Lab.a);

    wLab[2] = ab2Fix2(Lab.b);

}

static

cmsFloat64Number Clamp_L_doubleV4(cmsFloat64Number L)

{

    if (L < 0) L = 0;

    if (L > 100.0) L = 100.0;

    return L;

}

static

cmsFloat64Number Clamp_ab_doubleV4(cmsFloat64Number ab)

{

    if (ab < MIN_ENCODEABLE_ab4) ab = MIN_ENCODEABLE_ab4;

    if (ab > MAX_ENCODEABLE_ab4) ab = MAX_ENCODEABLE_ab4;

    return ab;

}

static

cmsUInt16Number L2Fix4(cmsFloat64Number L)

{

    return _cmsQuickSaturateWord(L *  655.35);

}

static

cmsUInt16Number ab2Fix4(cmsFloat64Number ab)

{

    return _cmsQuickSaturateWord((ab + 128.0) * 257.0);

}

void CMSEXPORT cmsFloat2LabEncoded(cmsUInt16Number wLab[3], const cmsCIELab* fLab)

{

    cmsCIELab Lab;

    Lab.L = Clamp_L_doubleV4(fLab ->L);

    Lab.a = Clamp_ab_doubleV4(fLab ->a);

    Lab.b = Clamp_ab_doubleV4(fLab ->b);

    wLab[0] = L2Fix4(Lab.L);

    wLab[1] = ab2Fix4(Lab.a);

    wLab[2] = ab2Fix4(Lab.b);

}

// Auxiliary: convert to Radians

static

cmsFloat64Number RADIANS(cmsFloat64Number deg)

{

    return (deg * M_PI) / 180.;

}

// Auxiliary: atan2 but operating in degrees and returning 0 if a==b==0

static

cmsFloat64Number atan2deg(cmsFloat64Number a, cmsFloat64Number b)

{

   cmsFloat64Number h;

   if (a == 0 && b == 0)

            h   = 0;

    else

            h = atan2(a, b);

    h *= (180. / M_PI);

    while (h > 360.)

        h -= 360.;

    while ( h < 0)

        h += 360.;

    return h;

}

// Auxiliary: Square

static

cmsFloat64Number Sqr(cmsFloat64Number v)

{

    return v *  v;

}

// From cylindrical coordinates. No check is performed, then negative values are allowed

void CMSEXPORT cmsLab2LCh(cmsCIELCh* LCh, const cmsCIELab* Lab)

{

    LCh -> L = Lab -> L;

    LCh -> C = pow(Sqr(Lab ->a) + Sqr(Lab ->b), 0.5);

    LCh -> h = atan2deg(Lab ->b, Lab ->a);

}

// To cylindrical coordinates. No check is performed, then negative values are allowed

void CMSEXPORT cmsLCh2Lab(cmsCIELab* Lab, const cmsCIELCh* LCh)

{

    cmsFloat64Number h = (LCh -> h * M_PI) / 180.0;

    Lab -> L = LCh -> L;

    Lab -> a = LCh -> C * cos(h);

    Lab -> b = LCh -> C * sin(h);

}

// In XYZ All 3 components are encoded using 1.15 fixed point

static

cmsUInt16Number XYZ2Fix(cmsFloat64Number d)

{

    return _cmsQuickSaturateWord(d * 32768.0);

}

void CMSEXPORT cmsFloat2XYZEncoded(cmsUInt16Number XYZ[3], const cmsCIEXYZ* fXYZ)

{

    cmsCIEXYZ xyz;

    xyz.X = fXYZ -> X;

    xyz.Y = fXYZ -> Y;

    xyz.Z = fXYZ -> Z;

    // Clamp to encodeable values.

    if (xyz.Y <= 0) {

        xyz.X = 0;

        xyz.Y = 0;

        xyz.Z = 0;

    }

    if (xyz.X > MAX_ENCODEABLE_XYZ)

        xyz.X = MAX_ENCODEABLE_XYZ;

    if (xyz.X < 0)

        xyz.X = 0;

    if (xyz.Y > MAX_ENCODEABLE_XYZ)

        xyz.Y = MAX_ENCODEABLE_XYZ;

    if (xyz.Y < 0)

        xyz.Y = 0;

    if (xyz.Z > MAX_ENCODEABLE_XYZ)

        xyz.Z = MAX_ENCODEABLE_XYZ;

    if (xyz.Z < 0)

        xyz.Z = 0;

    XYZ[0] = XYZ2Fix(xyz.X);

    XYZ[1] = XYZ2Fix(xyz.Y);

    XYZ[2] = XYZ2Fix(xyz.Z);

}

//  To convert from Fixed 1.15 point to cmsFloat64Number

static

cmsFloat64Number XYZ2float(cmsUInt16Number v)

{

    cmsS15Fixed16Number fix32;

    // From 1.15 to 15.16

    fix32 = v << 1;

    // From fixed 15.16 to cmsFloat64Number

    return _cms15Fixed16toDouble(fix32);

}

void CMSEXPORT cmsXYZEncoded2Float(cmsCIEXYZ* fXYZ, const cmsUInt16Number XYZ[3])

{

    fXYZ -> X = XYZ2float(XYZ[0]);

    fXYZ -> Y = XYZ2float(XYZ[1]);

    fXYZ -> Z = XYZ2float(XYZ[2]);

}

// Returns dE on two Lab values

cmsFloat64Number CMSEXPORT cmsDeltaE(const cmsCIELab* Lab1, const cmsCIELab* Lab2)

{

    cmsFloat64Number dL, da, db;

    dL = fabs(Lab1 -> L - Lab2 -> L);

    da = fabs(Lab1 -> a - Lab2 -> a);

    db = fabs(Lab1 -> b - Lab2 -> b);

    return pow(Sqr(dL) + Sqr(da) + Sqr(db), 0.5);

}

// Return the CIE94 Delta E

cmsFloat64Number CMSEXPORT cmsCIE94DeltaE(const cmsCIELab* Lab1, const cmsCIELab* Lab2)

{

    cmsCIELCh LCh1, LCh2;

    cmsFloat64Number dE, dL, dC, dh, dhsq;

    cmsFloat64Number c12, sc, sh;

    dL = fabs(Lab1 ->L - Lab2 ->L);

    cmsLab2LCh(&LCh1, Lab1);

    cmsLab2LCh(&LCh2, Lab2);

    dC  = fabs(LCh1.C - LCh2.C);

    dE  = cmsDeltaE(Lab1, Lab2);

    dhsq = Sqr(dE) - Sqr(dL) - Sqr(dC);

    if (dhsq < 0)

        dh = 0;

    else

        dh = pow(dhsq, 0.5);

    c12 = sqrt(LCh1.C * LCh2.C);