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

//  Little Color Management System

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

//

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

//

#include "lcms2_internal.h"

// Tone curves are powerful constructs that can contain curves specified in diverse ways.

// The curve is stored in segments, where each segment can be sampled or specified by parameters.

// a 16.bit simplification of the *whole* curve is kept for optimization purposes. For float operation,

// each segment is evaluated separately. Plug-ins may be used to define new parametric schemes,

// each plug-in may define up to MAX_TYPES_IN_LCMS_PLUGIN functions types. For defining a function,

// the plug-in should provide the type id, how many parameters each type has, and a pointer to

// a procedure that evaluates the function. In the case of reverse evaluation, the evaluator will

// be called with the type id as a negative value, and a sampled version of the reversed curve

// will be built.

// ----------------------------------------------------------------- Implementation

// Maxim number of nodes

#define MAX_NODES_IN_CURVE   4097

#define MINUS_INF            (-1E22F)

#define PLUS_INF             (+1E22F)

// The list of supported parametric curves

typedef struct _cmsParametricCurvesCollection_st {

    int nFunctions;                                     // Number of supported functions in this chunk

    int FunctionTypes[MAX_TYPES_IN_LCMS_PLUGIN];        // The identification types

    int ParameterCount[MAX_TYPES_IN_LCMS_PLUGIN];       // Number of parameters for each function

    cmsParametricCurveEvaluator    Evaluator;           // The evaluator

    struct _cmsParametricCurvesCollection_st* Next; // Next in list

} _cmsParametricCurvesCollection;

// This is the default (built-in) evaluator

static cmsFloat64Number DefaultEvalParametricFn(cmsInt32Number Type, const cmsFloat64Number Params[], cmsFloat64Number R);

// The built-in list

static _cmsParametricCurvesCollection DefaultCurves = {

    9,                                  // # of curve types

    { 1, 2, 3, 4, 5, 6, 7, 8, 108 },    // Parametric curve ID

    { 1, 3, 4, 5, 7, 4, 5, 5, 1 },      // Parameters by type

    DefaultEvalParametricFn,            // Evaluator

    NULL                                // Next in chain

};

// The linked list head

static _cmsParametricCurvesCollection* ParametricCurves = &DefaultCurves;

// As a way to install new parametric curves

cmsBool _cmsRegisterParametricCurvesPlugin(cmsPluginBase* Data)

{

    cmsPluginParametricCurves* Plugin = (cmsPluginParametricCurves*) Data;

    _cmsParametricCurvesCollection* fl;

    if (Data == NULL) {

          ParametricCurves =  &DefaultCurves;

          return TRUE;

    }

    fl = (_cmsParametricCurvesCollection*) _cmsPluginMalloc(sizeof(_cmsParametricCurvesCollection));

    if (fl == NULL) return FALSE;

    // Copy the parameters

    fl ->Evaluator  = Plugin ->Evaluator;

    fl ->nFunctions = Plugin ->nFunctions;

    // Make sure no mem overwrites

    if (fl ->nFunctions > MAX_TYPES_IN_LCMS_PLUGIN)

        fl ->nFunctions = MAX_TYPES_IN_LCMS_PLUGIN;

    // Copy the data

    memmove(fl->FunctionTypes,  Plugin ->FunctionTypes,   fl->nFunctions * sizeof(cmsUInt32Number));

    memmove(fl->ParameterCount, Plugin ->ParameterCount,  fl->nFunctions * sizeof(cmsUInt32Number));

    // Keep linked list

    fl ->Next = ParametricCurves;

    ParametricCurves = fl;

    // All is ok

    return TRUE;

}

// Search in type list, return position or -1 if not found

int IsInSet(int Type, _cmsParametricCurvesCollection* c)

    int i;

    for (i=0; i < c ->nFunctions; i++)

        if (abs(Type) == c ->FunctionTypes[i]) return i;

    return -1;

}

// Search for the collection which contains a specific type

static

_cmsParametricCurvesCollection *GetParametricCurveByType(int Type, int* index)

{

    _cmsParametricCurvesCollection* c;

    int Position;

    for (c = ParametricCurves; c != NULL; c = c ->Next) {

        Position = IsInSet(Type, c);

        if (Position != -1) {

            if (index != NULL)

                *index = Position;

            return c;

    return NULL;

// Low level allocate, which takes care of memory details. nEntries may be zero, and in this case

// no optimation curve is computed. nSegments may also be zero in the inverse case, where only the

// optimization curve is given. Both features simultaneously is an error

static

cmsToneCurve* AllocateToneCurveStruct(cmsContext ContextID, cmsInt32Number nEntries,

                                      cmsInt32Number nSegments, const cmsCurveSegment* Segments,

                                      const cmsUInt16Number* Values)

{

    cmsToneCurve* p;

    int i;

    // We allow huge tables, which are then restricted for smoothing operations

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

        cmsSignalError(ContextID, cmsERROR_RANGE, "Couldn't create tone curve of more than 65530 entries");

        return NULL;

    }

    if (nEntries <= 0 && nSegments <= 0) {

        cmsSignalError(ContextID, cmsERROR_RANGE, "Couldn't create tone curve with zero segments and no table");

        return NULL;

    }

    // Allocate all required pointers, etc.

    p = (cmsToneCurve*) _cmsMallocZero(ContextID, sizeof(cmsToneCurve));

    if (!p) return NULL;

    // In this case, there are no segments

    if (nSegments <= 0) {

        p ->Segments = NULL;

        p ->Evals = NULL;

    }

    else {

        p ->Segments = (cmsCurveSegment*) _cmsCalloc(ContextID, nSegments, sizeof(cmsCurveSegment));

        if (p ->Segments == NULL) goto Error;

        p ->Evals    = (cmsParametricCurveEvaluator*) _cmsCalloc(ContextID, nSegments, sizeof(cmsParametricCurveEvaluator));

        if (p ->Evals == NULL) goto Error;

    }

    p -> nSegments = nSegments;

    // This 16-bit table contains a limited precision representation of the whole curve and is kept for

    // increasing xput on certain operations.

    if (nEntries <= 0) {

        p ->Table16 = NULL;

    }

    else {

       p ->Table16 = (cmsUInt16Number*)  _cmsCalloc(ContextID, nEntries, sizeof(cmsUInt16Number));

       if (p ->Table16 == NULL) goto Error;

    }

    p -> nEntries  = nEntries;

    // Initialize members if requested

    if (Values != NULL && (nEntries > 0)) {

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

            p ->Table16[i] = Values[i];

    }

    // Initialize the segments stuff. The evaluator for each segment is located and a pointer to it

    // is placed in advance to maximize performance.

    if (Segments != NULL && (nSegments > 0)) {

        _cmsParametricCurvesCollection *c;

        p ->SegInterp = (cmsInterpParams**) _cmsCalloc(ContextID, nSegments, sizeof(cmsInterpParams*));

        if (p ->SegInterp == NULL) goto Error;

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

            // Type 0 is a special marker for table-based curves

            if (Segments[i].Type == 0)

                p ->SegInterp[i] = _cmsComputeInterpParams(ContextID, Segments[i].nGridPoints, 1, 1, NULL, CMS_LERP_FLAGS_FLOAT);

            memmove(&p ->Segments[i], &Segments[i], sizeof(cmsCurveSegment));

            if (Segments[i].Type == 0 && Segments[i].SampledPoints != NULL)

                p ->Segments[i].SampledPoints = (cmsFloat32Number*) _cmsDupMem(ContextID, Segments[i].SampledPoints, sizeof(cmsFloat32Number) * Segments[i].nGridPoints);

            else

                p ->Segments[i].SampledPoints = NULL;

            c = GetParametricCurveByType(Segments[i].Type, NULL);

            if (c != NULL)

                    p ->Evals[i] = c ->Evaluator;

        }

    }

    p ->InterpParams = _cmsComputeInterpParams(ContextID, p ->nEntries, 1, 1, p->Table16, CMS_LERP_FLAGS_16BITS);

    return p;

Error:

    if (p -> Segments) _cmsFree(ContextID, p ->Segments);

    if (p -> Evals) _cmsFree(ContextID, p -> Evals);

    if (p ->Table16) _cmsFree(ContextID, p ->Table16);

    _cmsFree(ContextID, p);

    return NULL;

// Parametric Fn using floating point

static

cmsFloat64Number DefaultEvalParametricFn(cmsInt32Number Type, const cmsFloat64Number Params[], cmsFloat64Number R)

    cmsFloat64Number e, Val, disc;

    switch (Type) {

    // X = Y ^ Gamma

    case 1:

        if (R < 0)

            Val = 0;

        else

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

        break;

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

    case -1:

        if (R < 0)

            Val = 0;

        else

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

        break;

    // CIE 122-1966

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

    // Y = 0               | else

    case 2:

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

        if (R >= disc ) {

            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:

         if (R < 0)

             Val = 0;

         else

             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:

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

        if (disc < 0)

            disc = 0;

        if (R >= disc) {

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

            if (e > 0)

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

            else

                Val = 0;

        }

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

            if (e > 0)

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

            else

                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:

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

        if (e < 0)

            disc = 0;

        else

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

        if (R >= disc) {

            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             | X < d

    case 5:

        if (R >= Params[4]) {

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

            if (e > 0)

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

            else

                Val = 0;

        }

        else

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

        break;

    // Reversed type 5

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

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

    case -5:

        disc = Params[3] * Params[4] + Params[6];

        if (R >= disc) {

            e = R - Params[5];

            if (e < 0)

                Val = 0;

            else

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

        }

        else {

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

        }

        break;

    // Types 6,7,8 comes from segmented curves as described in ICCSpecRevision_02_11_06_Float.pdf

    // Type 6 is basically identical to type 5 without d

    // Y = (a * X + b) ^ Gamma + c

    case 6:

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

        if (e < 0)

            Val = 0;

        else

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

        break;

    // ((Y - c) ^1/Gamma - b) / a

    case -6:

        e = R - Params[3];

        if (e < 0)

            Val = 0;

        else

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

        break;

    // Y = a * log (b * X^Gamma + c) + d

    case 7:

       e = Params[2] * pow(R, Params[0]) + Params[3];

       if (e <= 0)

           Val = 0;

       else

           Val = Params[1]*log10(e) + Params[4];

       break;

    // (Y - d) / a = log(b * X ^Gamma + c)

    // pow(10, (Y-d) / a) = b * X ^Gamma + c

    // pow((pow(10, (Y-d) / a) - c) / b, 1/g) = X

    case -7:

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

       break;

   //Y = a * b^(c*X+d) + e

   case 8:

       Val = (Params[0] * pow(Params[1], Params[2] * R + Params[3]) + Params[4]);

       break;

   // Y = (log((y-e) / a) / log(b) - d ) / c

   // a=0, b=1, c=2, d=3, e=4,

   case -8:

       disc = R - Params[4];

       if (disc < 0) Val = 0;

       else