/*

 * 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.  Sun designates this

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

 * by Sun 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 Sun Microsystems, Inc., 4150 Network Circle, Santa Clara,

 * CA 95054 USA or visit www.sun.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 cms

//  Copyright (C) 1998-2006 Marti Maria

//

// 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 "lcms.h"

/*

Gamut check by default is a catching of 0xFFFF/0xFFFF/0xFFFF PCS values, used

internally by lcms to hold invalid values. Matrix LUT's, operates in a way that

unencodeable values are marked as this combination, if PCS is XYZ, this is a very

high value since encoding is a 1.15 fixed point, something like 1.9997, 1.9997, 1.9997

not a very common color after all. Lab PCS is not to be a problem, since L>100 are truely

undefined. There is a posibility than ICC comitee defines L>100 as a valid means

to use highlights, then it will be lost.

(1.10 - Actually ICC did it, so this should be checked for full ICC 4.0 support)

*/

BOOL _cmsEndPointsBySpace(icColorSpaceSignature Space, WORD **White, WORD **Black,

                            int *nOutputs)

{

       // Only most common spaces

       static WORD RGBblack[4]  = { 0, 0, 0 };

       static WORD RGBwhite[4]  = { 0xffff, 0xffff, 0xffff };

       static WORD CMYKblack[4] = { 0xffff, 0xffff, 0xffff, 0xffff };   // 400% of ink

       static WORD CMYKwhite[4] = { 0, 0, 0, 0 };

       static WORD LABblack[4]  = { 0, 0x8000, 0x8000 };

       static WORD LABwhite[4]  = { 0xFF00, 0x8000, 0x8000 };

       static WORD CMYblack[4]  = { 0xffff, 0xffff, 0xffff };

       static WORD CMYwhite[4]  = { 0, 0, 0 };

       static WORD Grayblack[4] = { 0 };

       static WORD GrayWhite[4] = { 0xffff };

       switch (Space) {

       case icSigGrayData: if (White)    *White = GrayWhite;

                           if (Black)    *Black = Grayblack;

                           if (nOutputs) *nOutputs = 1;

                           return TRUE;

       case icSigRgbData:  if (White)    *White = RGBwhite;

                           if (Black)    *Black = RGBblack;

                           if (nOutputs) *nOutputs = 3;

                           return TRUE;

       case icSigLabData:  if (White)    *White = LABwhite;

                           if (Black)    *Black = LABblack;

                           if (nOutputs) *nOutputs = 3;

                           return TRUE;

       case icSigCmykData: if (White)    *White = CMYKwhite;

                           if (Black)    *Black = CMYKblack;

                           if (nOutputs) *nOutputs = 4;

                           return TRUE;

       case icSigCmyData:  if (White)    *White = CMYwhite;

                           if (Black)    *Black = CMYblack;

                           if (nOutputs) *nOutputs = 3;

                           return TRUE;

       default:;

       }

  return FALSE;

}

WORD *_cmsWhiteBySpace(icColorSpaceSignature Space)

{

       WORD *White= NULL, *Black = NULL;

       int Dummy;

       static WORD Default[MAXCHANNELS];

       if (_cmsEndPointsBySpace(Space, &White, &Black, &Dummy))

              return White;

       return Default;

}

WORD Clamp_L(Fixed32 in)

{

       if (in == 0xFFFF) return 0xFFFFU;  // Marker

       if (in > 0xFF00) return 0xFF00U;  // L* = 100.0

       return (WORD) in;

}

#define ENCODE_AB(x) (WORD) (((x) + 128.0) * 256.0 + 0.5)

WORD Clamp_ab(Fixed32 in)

{

       if (in == 0xFFFF) return 0xFFFFU;             // Marker

       if (in < 0) return ENCODE_AB(-128.0);         // Max negative number

       if (in > 0xFFFF) return ENCODE_AB(+127.9961); // Max positive number

       return (WORD) in;

}

// Returns dE on two Lab values

double LCMSEXPORT cmsDeltaE(LPcmsCIELab Lab1, LPcmsCIELab Lab2)

{

        double dL, da, db;

        if (Lab1 -> L < 0 ||

            Lab2 -> L < 0) return 65536.;

        if (Lab1 -> a < -200 || Lab1 -> a > 200) return 65536.;

        if (Lab1 -> b < -200 || Lab1 -> b > 200) return 65536.;

        if (Lab2 -> a < -200 || Lab2 -> a > 200) return 65536.;

        if (Lab2 -> b < -200 || Lab2 -> b > 200) return 65536.;

        if (Lab1 ->L == 0 && Lab2 ->L == 0) return 0;

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

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

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

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

}

// Square

static

double Sqr(double v)

{

    return v *  v;

}

// Return the CIE94 Delta E

double LCMSEXPORT cmsCIE94DeltaE(LPcmsCIELab Lab1, LPcmsCIELab Lab2)

{

    cmsCIELCh LCh1, LCh2;

    double dE, dL, dC, dh, dhsq;

    double c12, sc, sh;

    if (Lab1 ->L == 0 && Lab2 ->L == 0) return 0;

    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);

    sc = 1.0 + (0.048 * c12);

    sh = 1.0 + (0.014 * c12);

    return sqrt(Sqr(dL)  + Sqr(dC) / Sqr(sc) + Sqr(dh) / Sqr(sh));

}

// Auxiliary

static

double ComputeLBFD(LPcmsCIELab Lab)

{

  double yt;

  if (Lab->L > 7.996969)

        yt = (Sqr((Lab->L+16)/116)*((Lab->L+16)/116))*100;

  else

        yt = 100 * (Lab->L / 903.3);

  return (54.6 * (LOGE * (log(yt + 1.5))) - 9.6);

}

// bfd - gets BFD(1:1) difference between Lab1, Lab2

double LCMSEXPORT cmsBFDdeltaE(LPcmsCIELab Lab1, LPcmsCIELab Lab2)

{

    double lbfd1,lbfd2,AveC,Aveh,dE,deltaL,

        deltaC,deltah,dc,t,g,dh,rh,rc,rt,bfd;

    cmsCIELCh LCh1, LCh2;

    if (Lab1 ->L == 0 && Lab2 ->L == 0) return 0;

    lbfd1 = ComputeLBFD(Lab1);

    lbfd2 = ComputeLBFD(Lab2);

    deltaL = lbfd2 - lbfd1;

    cmsLab2LCh(&LCh1, Lab1);

    cmsLab2LCh(&LCh2, Lab2);

    deltaC = LCh2.C - LCh1.C;

    AveC = (LCh1.C+LCh2.C)/2;

    Aveh = (LCh1.h+LCh2.h)/2;

    dE = cmsDeltaE(Lab1, Lab2);

    if (Sqr(dE)>(Sqr(Lab2->L-Lab1->L)+Sqr(deltaC)))

        deltah = sqrt(Sqr(dE)-Sqr(Lab2->L-Lab1->L)-Sqr(deltaC));

    else

        deltah =0;

    dc   = 0.035 * AveC / (1 + 0.00365 * AveC)+0.521;

    g    = sqrt(Sqr(Sqr(AveC))/(Sqr(Sqr(AveC))+14000));

    t    = 0.627+(0.055*cos((Aveh-254)/(180/M_PI))-

        0.040*cos((2*Aveh-136)/(180/M_PI))+

        0.070*cos((3*Aveh-31)/(180/M_PI))+

        0.049*cos((4*Aveh+114)/(180/M_PI))-

        0.015*cos((5*Aveh-103)/(180/M_PI)));

    dh    = dc*(g*t+1-g);

    rh    = -0.260*cos((Aveh-308)/(180/M_PI))-

        0.379*cos((2*Aveh-160)/(180/M_PI))-

        0.636*cos((3*Aveh+254)/(180/M_PI))+

        0.226*cos((4*Aveh+140)/(180/M_PI))-

        0.194*cos((5*Aveh+280)/(180/M_PI));

    rc = sqrt((AveC*AveC*AveC*AveC*AveC*AveC)/((AveC*AveC*AveC*AveC*AveC*AveC)+70000000));

    rt = rh*rc;

    bfd = sqrt(Sqr(deltaL)+Sqr(deltaC/dc)+Sqr(deltah/dh)+(rt*(deltaC/dc)*(deltah/dh)));

    return bfd;

}

//  cmc - CMC(1:1) difference between Lab1, Lab2

double LCMSEXPORT cmsCMCdeltaE(LPcmsCIELab Lab1, LPcmsCIELab Lab2)

{

  double dE,dL,dC,dh,sl,sc,sh,t,f,cmc;

  cmsCIELCh LCh1, LCh2;

  if (Lab1 ->L == 0 && Lab2 ->L == 0) return 0;

  cmsLab2LCh(&LCh1, Lab1);

  cmsLab2LCh(&LCh2, Lab2);

  dL = Lab2->L-Lab1->L;

  dC = LCh2.C-LCh1.C;

  dE = cmsDeltaE(Lab1, Lab2);

  if (Sqr(dE)>(Sqr(dL)+Sqr(dC)))

            dh = sqrt(Sqr(dE)-Sqr(dL)-Sqr(dC));

  else

            dh =0;

  if ((LCh1.h > 164) && (LCh1.h<345))

      t = 0.56 + fabs(0.2 * cos(((LCh1.h + 168)/(180/M_PI))));

  else

      t = 0.36 + fabs(0.4 * cos(((LCh1.h + 35 )/(180/M_PI))));

   sc  = 0.0638   * LCh1.C / (1 + 0.0131  * LCh1.C) + 0.638;

   sl  = 0.040975 * Lab1->L /(1 + 0.01765 * Lab1->L);

   if (Lab1->L<16)

         sl = 0.511;

   f   = sqrt((LCh1.C * LCh1.C * LCh1.C * LCh1.C)/((LCh1.C * LCh1.C * LCh1.C * LCh1.C)+1900));

   sh  = sc*(t*f+1-f);

   cmc = sqrt(Sqr(dL/sl)+Sqr(dC/sc)+Sqr(dh/sh));

   return cmc;

}

static

double atan2deg(double b, double a)

{

   double 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;

}

static

double RADIANES(double deg)

{

    return (deg * M_PI) / 180.;

}

// dE2000 The weightings KL, KC and KH can be modified to reflect the relative

// importance of lightness, chroma and hue in different industrial applications

double LCMSEXPORT cmsCIE2000DeltaE(LPcmsCIELab Lab1, LPcmsCIELab Lab2,

                                  double Kl, double Kc, double Kh)

{

    double L1  = Lab1->L;

    double a1  = Lab1->a;

    double b1  = Lab1->b;

    double C   = sqrt( Sqr(a1) + Sqr(b1) );

    double Ls = Lab2 ->L;

    double as = Lab2 ->a;

    double bs = Lab2 ->b;

    double Cs = sqrt( Sqr(as) + Sqr(bs) );

    double G = 0.5 * ( 1 - sqrt(pow((C + Cs) / 2 , 7.0) / (pow((C + Cs) / 2, 7.0) + pow(25.0, 7.0) ) ));

    double a_p = (1 + G ) * a1;

    double b_p = b1;

    double C_p = sqrt( Sqr(a_p) + Sqr(b_p));

    double h_p = atan2deg(a_p, b_p);

    double a_ps = (1 + G) * as;

    double b_ps = bs;

    double C_ps = sqrt(Sqr(a_ps) + Sqr(b_ps));

    double h_ps = atan2deg(a_ps, b_ps);

    double meanC_p =(C_p + C_ps) / 2;

    double meanh_p = fabs(h_ps-h_p) <= 180 ? (h_ps + h_p)/2 : (h_ps+h_p-360)/2;

    double delta_h = fabs(h_p - h_ps) <= 180 ? fabs(h_p - h_ps) : 360 - fabs(h_p - h_ps);

    double delta_L = fabs(L1 - Ls);

    double delta_C = fabs(C_p - C_ps);

    double delta_H =2 * sqrt(C_ps*C_p) * sin(RADIANES(delta_h) / 2);

    double T = 1 - 0.17 * cos(RADIANES(meanh_p-30))

                 + 0.24 * cos(RADIANES(2*meanh_p))

                 + 0.32 * cos(RADIANES(3*meanh_p + 6))

                 - 0.2  * cos(RADIANES(4*meanh_p - 63));

    double Sl = 1 + (0.015 * Sqr((Ls + L1) /2- 50) )/ sqrt(20 + Sqr( (Ls+L1)/2 - 50) );

    double Sc = 1 + 0.045 * (C_p + C_ps)/2;

    double Sh = 1 + 0.015 * ((C_ps + C_p)/2) * T;

    double delta_ro = 30 * exp( -Sqr(((meanh_p - 275 ) / 25)));

    double Rc = 2 * sqrt(( pow(meanC_p, 7.0) )/( pow(meanC_p, 7.0) + pow(25.0, 7.0)));

    double Rt = -sin(2 * RADIANES(delta_ro)) * Rc;

    double deltaE00 = sqrt( Sqr(delta_L /(Sl * Kl)) +

                            Sqr(delta_C/(Sc * Kc))  +

                            Sqr(delta_H/(Sh * Kh))  +

                            Rt*(delta_C/(Sc * Kc)) * (delta_H / (Sh * Kh)));

    return deltaE00;

}

// Carefully,  clamp on CIELab space.

void LCMSEXPORT cmsClampLab(LPcmsCIELab Lab, double amax, double amin,

                                   double bmax, double bmin)

{

            // Whole Luma surface to zero

        if (Lab -> L < 0) {

                Lab-> L = Lab->a = Lab-> b = 0.0;

                return;

            }

            // Clamp white, DISCARD HIGHLIGHTS. This is done

            // in such way because icc spec doesn't allow the

            // use of L>100 as a highlight means.

            if (Lab->L > 100)

                        Lab -> L = 100;

            // Check out gamut prism, on a, b faces

            if (Lab -> a < amin || Lab->a > amax||

                Lab -> b < bmin || Lab->b > bmax) {

                 cmsCIELCh LCh;

                 double h, slope;

                 // Falls outside a, b limits. Transports to LCh space,

                 // and then do the clipping

                 if (Lab -> a == 0.0) { // Is hue exactly 90?

                        // atan will not work, so clamp here

                        Lab -> b = Lab->b < 0 ? bmin : bmax;

                        return;

                 }

                 cmsLab2LCh(&LCh, Lab);

                 slope = Lab -> b / Lab -> a;

                 h = LCh.h;

                 // There are 4 zones

                 if ((h >= 0. && h < 45.) ||

                     (h >= 315 && h <= 360.)) {

                     // clip by amax

                     Lab -> a = amax;

                     Lab -> b = amax * slope;

                 }

                 else

                 if (h >= 45. && h < 135)

                 {

                        // clip by bmax

                        Lab -> b = bmax;

                        Lab -> a = bmax / slope;

                 }

                 else

                 if (h >= 135 && h < 225) {

                        // clip by amin

                        Lab -> a = amin;

                        Lab -> b = amin * slope;

                 }

                 else

                 if (h >= 225 && h < 315) {

                        // clip by bmin

                        Lab -> b = bmin;

                        Lab -> a = bmin / slope;

                 }

                 else

                        cmsSignalError(LCMS_ERRC_ABORTED, "Invalid angle");

        }

}

// Several utilities -------------------------------------------------------

// Translate from our colorspace to ICC representation

icColorSpaceSignature LCMSEXPORT _cmsICCcolorSpace(int OurNotation)

{

       switch (OurNotation) {

       case 1:

       case PT_GRAY: return  icSigGrayData;

       case 2:

       case PT_RGB:  return  icSigRgbData;

       case PT_CMY:  return  icSigCmyData;

       case PT_CMYK: return  icSigCmykData;

       case PT_YCbCr:return  icSigYCbCrData;

       case PT_YUV:  return  icSigLuvData;

       case PT_XYZ:  return  icSigXYZData;