/*

 * 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

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

 *

 * 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.

 *

 */

// 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-2007 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"

// Gamma handling.

LPGAMMATABLE LCMSEXPORT cmsAllocGamma(int nEntries);

void         LCMSEXPORT cmsFreeGamma(LPGAMMATABLE Gamma);

void         LCMSEXPORT cmsFreeGammaTriple(LPGAMMATABLE Gamma[3]);

LPGAMMATABLE LCMSEXPORT cmsBuildGamma(int nEntries, double Gamma);

LPGAMMATABLE LCMSEXPORT cmsDupGamma(LPGAMMATABLE Src);

LPGAMMATABLE LCMSEXPORT cmsReverseGamma(int nResultSamples, LPGAMMATABLE InGamma);

LPGAMMATABLE LCMSEXPORT cmsBuildParametricGamma(int nEntries, int Type, double Params[]);

LPGAMMATABLE LCMSEXPORT cmsJoinGamma(LPGAMMATABLE InGamma, LPGAMMATABLE OutGamma);

LPGAMMATABLE LCMSEXPORT cmsJoinGammaEx(LPGAMMATABLE InGamma, LPGAMMATABLE OutGamma, int nPoints);

LCMSBOOL         LCMSEXPORT cmsSmoothGamma(LPGAMMATABLE Tab, double lambda);

LCMSBOOL         cdecl _cmsSmoothEndpoints(LPWORD Table, int nPoints);

// Sampled curves

LPSAMPLEDCURVE cdecl cmsAllocSampledCurve(int nItems);

void           cdecl cmsFreeSampledCurve(LPSAMPLEDCURVE p);

void           cdecl cmsEndpointsOfSampledCurve(LPSAMPLEDCURVE p, double* Min, double* Max);

void           cdecl cmsClampSampledCurve(LPSAMPLEDCURVE p, double Min, double Max);

LCMSBOOL       cdecl cmsSmoothSampledCurve(LPSAMPLEDCURVE Tab, double SmoothingLambda);

void           cdecl cmsRescaleSampledCurve(LPSAMPLEDCURVE p, double Min, double Max, int nPoints);

LPSAMPLEDCURVE cdecl cmsJoinSampledCurves(LPSAMPLEDCURVE X, LPSAMPLEDCURVE Y, int nResultingPoints);

double LCMSEXPORT cmsEstimateGamma(LPGAMMATABLE t);

double LCMSEXPORT cmsEstimateGammaEx(LPWORD GammaTable, int nEntries, double Thereshold);

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

#define MAX_KNOTS   4096

typedef float vec[MAX_KNOTS+1];

// Ciclic-redundant-check for assuring table is a true representation of parametric curve

// The usual polynomial, which is used for AAL5, FDDI, and probably Ethernet

#define QUOTIENT 0x04c11db7

static

unsigned int Crc32(unsigned int result, LPVOID ptr, int len)

{

    int          i,j;

    BYTE         octet;

    LPBYTE       data = (LPBYTE) ptr;

    for (i=0; i < len; i++) {

        octet = *data++;

        for (j=0; j < 8; j++) {

            if (result & 0x80000000) {

                result = (result << 1) ^ QUOTIENT ^ (octet >> 7);

            }

            else

            {

                result = (result << 1) ^ (octet >> 7);

            }

            octet <<= 1;

        }

    }

    return result;

}

// Get CRC of gamma table

unsigned int _cmsCrc32OfGammaTable(LPGAMMATABLE Table)

{

    unsigned int crc = ~0U;

    crc = Crc32(crc, &Table -> Seed.Type,  sizeof(int));

    crc = Crc32(crc, Table ->Seed.Params,  sizeof(double)*10);

    crc = Crc32(crc, &Table ->nEntries,    sizeof(int));

    crc = Crc32(crc, Table ->GammaTable,   sizeof(WORD) * Table -> nEntries);

    return ~crc;

}

LPGAMMATABLE LCMSEXPORT cmsAllocGamma(int nEntries)

{

       LPGAMMATABLE p;

       size_t size;

       if (nEntries > 65530 || nEntries <= 0) {

                cmsSignalError(LCMS_ERRC_ABORTED, "Couldn't create gammatable of more than 65530 entries");

                return NULL;

       }

       size = sizeof(GAMMATABLE) + (sizeof(WORD) * (nEntries-1));

       p = (LPGAMMATABLE) _cmsMalloc(size);

       if (!p) return NULL;

       ZeroMemory(p, size);

       p -> Seed.Type     = 0;

       p -> nEntries = nEntries;

       return p;

}

void LCMSEXPORT cmsFreeGamma(LPGAMMATABLE Gamma)

{

       if (Gamma)  _cmsFree(Gamma);

}

void LCMSEXPORT cmsFreeGammaTriple(LPGAMMATABLE Gamma[3])

{

    cmsFreeGamma(Gamma[0]);

    cmsFreeGamma(Gamma[1]);

    cmsFreeGamma(Gamma[2]);

    Gamma[0] = Gamma[1] = Gamma[2] = NULL;

}

// Duplicate a gamma table

LPGAMMATABLE  LCMSEXPORT cmsDupGamma(LPGAMMATABLE In)

{

       LPGAMMATABLE Ptr;

       size_t size;

       Ptr = cmsAllocGamma(In -> nEntries);

       if (Ptr == NULL) return NULL;

       size = sizeof(GAMMATABLE) + (sizeof(WORD) * (In -> nEntries-1));

       CopyMemory(Ptr, In, size);

       return Ptr;

}

// Handle gamma using interpolation tables. The resulting curves can become

// very stange, but are pleasent to eye.

LPGAMMATABLE LCMSEXPORT cmsJoinGamma(LPGAMMATABLE InGamma,

                          LPGAMMATABLE OutGamma)

{

       register int i;

       L16PARAMS L16In, L16Out;

       LPWORD InPtr, OutPtr;

       LPGAMMATABLE p;

       p = cmsAllocGamma(256);

       if (!p) return NULL;

       cmsCalcL16Params(InGamma -> nEntries, &L16In);

       InPtr  = InGamma -> GammaTable;

       cmsCalcL16Params(OutGamma -> nEntries, &L16Out);

       OutPtr = OutGamma-> GammaTable;

       for (i=0; i < 256; i++)

       {

              WORD wValIn, wValOut;

              wValIn  = cmsLinearInterpLUT16(RGB_8_TO_16(i), InPtr, &L16In);

              wValOut = cmsReverseLinearInterpLUT16(wValIn, OutPtr, &L16Out);

              p -> GammaTable[i] = wValOut;

       }

       return p;

}

// New method, using smoothed parametric curves. This works FAR better.

// We want to get

//

//      y = f(g^-1(x))      ; f = ingamma, g = outgamma

//

// And this can be parametrized as

//

//      y = f(t)

//      x = g(t)

LPGAMMATABLE LCMSEXPORT cmsJoinGammaEx(LPGAMMATABLE InGamma,

                                       LPGAMMATABLE OutGamma, int nPoints)

{

    LPSAMPLEDCURVE x, y, r;

    LPGAMMATABLE res;

    x = cmsConvertGammaToSampledCurve(InGamma,  nPoints);

    y = cmsConvertGammaToSampledCurve(OutGamma, nPoints);

    r = cmsJoinSampledCurves(y, x, nPoints);

    // Does clean "hair"

    cmsSmoothSampledCurve(r, 0.001);

    cmsClampSampledCurve(r, 0.0, 65535.0);

    cmsFreeSampledCurve(x);

    cmsFreeSampledCurve(y);

    res = cmsConvertSampledCurveToGamma(r, 65535.0);

    cmsFreeSampledCurve(r);

    return res;

}

// Reverse a gamma table

LPGAMMATABLE LCMSEXPORT cmsReverseGamma(int nResultSamples, LPGAMMATABLE InGamma)

{

       register int i;

       L16PARAMS L16In;

       LPWORD InPtr;

       LPGAMMATABLE p;

       // Try to reverse it analytically whatever possible

       if (InGamma -> Seed.Type > 0 && InGamma -> Seed.Type <= 5 &&

            _cmsCrc32OfGammaTable(InGamma) == InGamma -> Seed.Crc32) {

                return cmsBuildParametricGamma(nResultSamples, -(InGamma -> Seed.Type), InGamma ->Seed.Params);

       }

       // Nope, reverse the table

       p = cmsAllocGamma(nResultSamples);

       if (!p) return NULL;

       cmsCalcL16Params(InGamma -> nEntries, &L16In);

       InPtr  = InGamma -> GammaTable;

       for (i=0; i < nResultSamples; i++)

       {

              WORD wValIn, wValOut;

              wValIn = _cmsQuantizeVal(i, nResultSamples);

              wValOut = cmsReverseLinearInterpLUT16(wValIn, InPtr, &L16In);

              p -> GammaTable[i] = wValOut;

       }

       return p;

}

// Parametric curves

//

// Parameters goes as: Gamma, a, b, c, d, e, f

// Type is the ICC type +1

// if type is negative, then the curve is analyticaly inverted

LPGAMMATABLE LCMSEXPORT cmsBuildParametricGamma(int nEntries, int Type, double Params[])

{

        LPGAMMATABLE Table;

        double R, Val, dval, e;

        int i;

        int ParamsByType[] = { 0, 1, 3, 4, 5, 7 };

        Table = cmsAllocGamma(nEntries);

        if (NULL == Table) return NULL;

        Table -> Seed.Type = Type;

        CopyMemory(Table ->Seed.Params, Params, ParamsByType[abs(Type)] * sizeof(double));

        for (i=0; i < nEntries; i++) {

                R   = (double) i / (nEntries-1);

                switch (Type) {

                // X = Y ^ Gamma

                case 1:

                      Val = pow(R, Params[0]);

                      break;

                // Type 1 Reversed: X = Y ^1/gamma

                case -1:

                      Val = pow(R, 1/Params[0]);

                      break;

                // CIE 122-1966

                // Y = (aX + b)^Gamma  | X >= -b/a

                // Y = 0               | else

                case 2:

                    if (R >= -Params[2] / Params[1]) {

                              e = Params[1]*R + Params[2];

                              if (e > 0)

                                Val = pow(e, Params[0]);

                              else

                                Val = 0;

                    }

                    else

                              Val = 0;

                      break;

                // Type 2 Reversed

                // X = (Y ^1/g  - b) / a

                case -2:

                    Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1];

                    if (Val < 0)

                            Val = 0;

                    break;

                // IEC 61966-3

                // Y = (aX + b)^Gamma | X <= -b/a

                // Y = c              | else

                case 3:

                    if (R >= -Params[2] / Params[1]) {

                      e = Params[1]*R + Params[2];

                      Val = pow(e, Params[0]) + Params[3];

                    }

                    else

                      Val = Params[3];

                    break;

                // Type 3 reversed

                // X=((Y-c)^1/g - b)/a      | (Y>=c)

                // X=-b/a                   | (Y<c)

                case -3:

                    if (R >= Params[3])  {

                        e = R - Params[3];

                        Val = (pow(e, 1/Params[0]) - Params[2]) / Params[1];

                        if (Val < 0) Val = 0;

                    }

                    else {

                        Val = -Params[2] / Params[1];

                    }

                    break;

                // IEC 61966-2.1 (sRGB)

                // Y = (aX + b)^Gamma | X >= d

                // Y = cX             | X < d

                case 4:

                    if (R >= Params[4]) {

                              e = Params[1]*R + Params[2];

                              if (e > 0)

                                Val = pow(e, Params[0]);

                              else

                                Val = 0;

                    }

                      else

                              Val = R * Params[3];

                      break;

                // Type 4 reversed

                // X=((Y^1/g-b)/a)    | Y >= (ad+b)^g

                // X=Y/c              | Y< (ad+b)^g

                case -4:

                    if (R >= pow(Params[1] * Params[4] + Params[2], Params[0])) {

                        Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1];

                    }

                    else {

                        Val = R / Params[3];

                    }

                    break;

                // Y = (aX + b)^Gamma + e | X <= d

                // Y = cX + f             | else

                case 5:

                    if (R >= Params[4]) {

                        e = Params[1]*R + Params[2];

                        Val = pow(e, Params[0]) + Params[5];

                    }

                    else

                        Val = R*Params[3] + Params[6];

                    break;

                // Reversed type 5

                // X=((Y-e)1/g-b)/a   | Y >=(ad+b)^g+e)

                // X=(Y-f)/c          | else

                case -5:

                if (R >= pow(Params[1] * Params[4], Params[0]) + Params[5]) {

                    Val = pow(R - Params[5], 1/Params[0]) - Params[2] / Params[1];

                }

                else {

                    Val = (R - Params[6]) / Params[3];

                }

                break;

                default:

                        cmsSignalError(LCMS_ERRC_ABORTED, "Unsupported parametric curve type=%d", abs(Type)-1);

                        cmsFreeGamma(Table);

                        return NULL;

                }

        // Saturate

        dval = Val * 65535.0 + .5;

        if (dval > 65535.) dval = 65535.0;

        if (dval < 0) dval = 0;

        Table->GammaTable[i] = (WORD) floor(dval);

        }

        Table -> Seed.Crc32 = _cmsCrc32OfGammaTable(Table);

        return Table;

}

// Build a gamma table based on gamma constant

LPGAMMATABLE LCMSEXPORT cmsBuildGamma(int nEntries, double Gamma)

{

    return cmsBuildParametricGamma(nEntries, 1, &Gamma);

}

// From: Eilers, P.H.C. (1994) Smoothing and interpolation with finite

// differences. in: Graphic Gems IV, Heckbert, P.S. (ed.), Academic press.

//

// Smoothing and interpolation with second differences.

//

//   Input:  weights (w), data (y): vector from 1 to m.

//   Input:  smoothing parameter (lambda), length (m).

//   Output: smoothed vector (z): vector from 1 to m.

static

void smooth2(vec w, vec y, vec z, float lambda, int m)

{

  int i, i1, i2;

  vec c, d, e;

  d[1] = w[1] + lambda;

  c[1] = -2 * lambda / d[1];

  e[1] = lambda /d[1];

  z[1] = w[1] * y[1];

  d[2] = w[2] + 5 * lambda - d[1] * c[1] *  c[1];

  c[2] = (-4 * lambda - d[1] * c[1] * e[1]) / d[2];

  e[2] = lambda / d[2];

  z[2] = w[2] * y[2] - c[1] * z[1];

  for (i = 3; i < m - 1; i++) {

    i1 = i - 1; i2 = i - 2;

    d[i]= w[i] + 6 * lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];

    c[i] = (-4 * lambda -d[i1] * c[i1] * e[i1])/ d[i];

    e[i] = lambda / d[i];

    z[i] = w[i] * y[i] - c[i1] * z[i1] - e[i2] * z[i2];

  }

  i1 = m - 2; i2 = m - 3;

  d[m - 1] = w[m - 1] + 5 * lambda -c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];

  c[m - 1] = (-2 * lambda - d[i1] * c[i1] * e[i1]) / d[m - 1];