1//---------------------------------------------------------------------------------
2//
3//  Little Color Management System
4//  Copyright (c) 1998-2013 Marti Maria Saguer
5//
6// Permission is hereby granted, free of charge, to any person obtaining
7// a copy of this software and associated documentation files (the "Software"),
8// to deal in the Software without restriction, including without limitation
9// the rights to use, copy, modify, merge, publish, distribute, sublicense,
10// and/or sell copies of the Software, and to permit persons to whom the Software
11// is furnished to do so, subject to the following conditions:
12//
13// The above copyright notice and this permission notice shall be included in
14// all copies or substantial portions of the Software.
15//
16// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
18// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
20// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
21// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
22// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23//
24//---------------------------------------------------------------------------------
25//
26
27#include "lcms2_internal.h"
28
29// Tone curves are powerful constructs that can contain curves specified in diverse ways.
30// The curve is stored in segments, where each segment can be sampled or specified by parameters.
31// a 16.bit simplification of the *whole* curve is kept for optimization purposes. For float operation,
32// each segment is evaluated separately. Plug-ins may be used to define new parametric schemes,
33// each plug-in may define up to MAX_TYPES_IN_LCMS_PLUGIN functions types. For defining a function,
34// the plug-in should provide the type id, how many parameters each type has, and a pointer to
35// a procedure that evaluates the function. In the case of reverse evaluation, the evaluator will
36// be called with the type id as a negative value, and a sampled version of the reversed curve
37// will be built.
38
39// ----------------------------------------------------------------- Implementation
40// Maxim number of nodes
41#define MAX_NODES_IN_CURVE   4097
42#define MINUS_INF            (-1E22F)
43#define PLUS_INF             (+1E22F)
44
45// The list of supported parametric curves
46typedef struct _cmsParametricCurvesCollection_st {
47
48    int nFunctions;                                     // Number of supported functions in this chunk
49    int FunctionTypes[MAX_TYPES_IN_LCMS_PLUGIN];        // The identification types
50    int ParameterCount[MAX_TYPES_IN_LCMS_PLUGIN];       // Number of parameters for each function
51    cmsParametricCurveEvaluator    Evaluator;           // The evaluator
52
53    struct _cmsParametricCurvesCollection_st* Next; // Next in list
54
55} _cmsParametricCurvesCollection;
56
57// This is the default (built-in) evaluator
58static cmsFloat64Number DefaultEvalParametricFn(cmsInt32Number Type, const cmsFloat64Number Params[], cmsFloat64Number R);
59
60// The built-in list
61static _cmsParametricCurvesCollection DefaultCurves = {
62    9,                                  // # of curve types
63    { 1, 2, 3, 4, 5, 6, 7, 8, 108 },    // Parametric curve ID
64    { 1, 3, 4, 5, 7, 4, 5, 5, 1 },      // Parameters by type
65    DefaultEvalParametricFn,            // Evaluator
66    NULL                                // Next in chain
67};
68
69// Duplicates the zone of memory used by the plug-in in the new context
70static
71void DupPluginCurvesList(struct _cmsContext_struct* ctx,
72                                               const struct _cmsContext_struct* src)
73{
74   _cmsCurvesPluginChunkType newHead = { NULL };
75   _cmsParametricCurvesCollection*  entry;
76   _cmsParametricCurvesCollection*  Anterior = NULL;
77   _cmsCurvesPluginChunkType* head = (_cmsCurvesPluginChunkType*) src->chunks[CurvesPlugin];
78
79    _cmsAssert(head != NULL);
80
81    // Walk the list copying all nodes
82   for (entry = head->ParametricCurves;
83        entry != NULL;
84        entry = entry ->Next) {
85
86            _cmsParametricCurvesCollection *newEntry = ( _cmsParametricCurvesCollection *) _cmsSubAllocDup(ctx ->MemPool, entry, sizeof(_cmsParametricCurvesCollection));
87
88            if (newEntry == NULL)
89                return;
90
91            // We want to keep the linked list order, so this is a little bit tricky
92            newEntry -> Next = NULL;
93            if (Anterior)
94                Anterior -> Next = newEntry;
95
96            Anterior = newEntry;
97
98            if (newHead.ParametricCurves == NULL)
99                newHead.ParametricCurves = newEntry;
100    }
101
102  ctx ->chunks[CurvesPlugin] = _cmsSubAllocDup(ctx->MemPool, &newHead, sizeof(_cmsCurvesPluginChunkType));
103}
104
105// The allocator have to follow the chain
106void _cmsAllocCurvesPluginChunk(struct _cmsContext_struct* ctx,
107                                const struct _cmsContext_struct* src)
108{
109    _cmsAssert(ctx != NULL);
110
111    if (src != NULL) {
112
113        // Copy all linked list
114       DupPluginCurvesList(ctx, src);
115    }
116    else {
117        static _cmsCurvesPluginChunkType CurvesPluginChunk = { NULL };
118        ctx ->chunks[CurvesPlugin] = _cmsSubAllocDup(ctx ->MemPool, &CurvesPluginChunk, sizeof(_cmsCurvesPluginChunkType));
119    }
120}
121
122
123// The linked list head
124_cmsCurvesPluginChunkType _cmsCurvesPluginChunk = { NULL };
125
126// As a way to install new parametric curves
127cmsBool _cmsRegisterParametricCurvesPlugin(cmsContext ContextID, cmsPluginBase* Data)
128{
129    _cmsCurvesPluginChunkType* ctx = ( _cmsCurvesPluginChunkType*) _cmsContextGetClientChunk(ContextID, CurvesPlugin);
130    cmsPluginParametricCurves* Plugin = (cmsPluginParametricCurves*) Data;
131    _cmsParametricCurvesCollection* fl;
132
133    if (Data == NULL) {
134
135          ctx -> ParametricCurves =  NULL;
136          return TRUE;
137    }
138
139    fl = (_cmsParametricCurvesCollection*) _cmsPluginMalloc(ContextID, sizeof(_cmsParametricCurvesCollection));
140    if (fl == NULL) return FALSE;
141
142    // Copy the parameters
143    fl ->Evaluator  = Plugin ->Evaluator;
144    fl ->nFunctions = Plugin ->nFunctions;
145
146    // Make sure no mem overwrites
147    if (fl ->nFunctions > MAX_TYPES_IN_LCMS_PLUGIN)
148        fl ->nFunctions = MAX_TYPES_IN_LCMS_PLUGIN;
149
150    // Copy the data
151    memmove(fl->FunctionTypes,  Plugin ->FunctionTypes,   fl->nFunctions * sizeof(cmsUInt32Number));
152    memmove(fl->ParameterCount, Plugin ->ParameterCount,  fl->nFunctions * sizeof(cmsUInt32Number));
153
154    // Keep linked list
155    fl ->Next = ctx->ParametricCurves;
156    ctx->ParametricCurves = fl;
157
158    // All is ok
159    return TRUE;
160}
161
162
163// Search in type list, return position or -1 if not found
164static
165int IsInSet(int Type, _cmsParametricCurvesCollection* c)
166{
167    int i;
168
169    for (i=0; i < c ->nFunctions; i++)
170        if (abs(Type) == c ->FunctionTypes[i]) return i;
171
172    return -1;
173}
174
175
176// Search for the collection which contains a specific type
177static
178_cmsParametricCurvesCollection *GetParametricCurveByType(cmsContext ContextID, int Type, int* index)
179{
180    _cmsParametricCurvesCollection* c;
181    int Position;
182    _cmsCurvesPluginChunkType* ctx = ( _cmsCurvesPluginChunkType*) _cmsContextGetClientChunk(ContextID, CurvesPlugin);
183
184    for (c = ctx->ParametricCurves; c != NULL; c = c ->Next) {
185
186        Position = IsInSet(Type, c);
187
188        if (Position != -1) {
189            if (index != NULL)
190                *index = Position;
191            return c;
192        }
193    }
194    // If none found, revert for defaults
195    for (c = &DefaultCurves; c != NULL; c = c ->Next) {
196
197        Position = IsInSet(Type, c);
198
199        if (Position != -1) {
200            if (index != NULL)
201                *index = Position;
202            return c;
203        }
204    }
205
206    return NULL;
207}
208
209// Low level allocate, which takes care of memory details. nEntries may be zero, and in this case
210// no optimation curve is computed. nSegments may also be zero in the inverse case, where only the
211// optimization curve is given. Both features simultaneously is an error
212static
213cmsToneCurve* AllocateToneCurveStruct(cmsContext ContextID, cmsInt32Number nEntries,
214                                      cmsInt32Number nSegments, const cmsCurveSegment* Segments,
215                                      const cmsUInt16Number* Values)
216{
217    cmsToneCurve* p;
218    int i;
219
220    // We allow huge tables, which are then restricted for smoothing operations
221    if (nEntries > 65530 || nEntries < 0) {
222        cmsSignalError(ContextID, cmsERROR_RANGE, "Couldn't create tone curve of more than 65530 entries");
223        return NULL;
224    }
225
226    if (nEntries <= 0 && nSegments <= 0) {
227        cmsSignalError(ContextID, cmsERROR_RANGE, "Couldn't create tone curve with zero segments and no table");
228        return NULL;
229    }
230
231    // Allocate all required pointers, etc.
232    p = (cmsToneCurve*) _cmsMallocZero(ContextID, sizeof(cmsToneCurve));
233    if (!p) return NULL;
234
235    // In this case, there are no segments
236    if (nSegments <= 0) {
237        p ->Segments = NULL;
238        p ->Evals = NULL;
239    }
240    else {
241        p ->Segments = (cmsCurveSegment*) _cmsCalloc(ContextID, nSegments, sizeof(cmsCurveSegment));
242        if (p ->Segments == NULL) goto Error;
243
244        p ->Evals    = (cmsParametricCurveEvaluator*) _cmsCalloc(ContextID, nSegments, sizeof(cmsParametricCurveEvaluator));
245        if (p ->Evals == NULL) goto Error;
246    }
247
248    p -> nSegments = nSegments;
249
250    // This 16-bit table contains a limited precision representation of the whole curve and is kept for
251    // increasing xput on certain operations.
252    if (nEntries <= 0) {
253        p ->Table16 = NULL;
254    }
255    else {
256       p ->Table16 = (cmsUInt16Number*)  _cmsCalloc(ContextID, nEntries, sizeof(cmsUInt16Number));
257       if (p ->Table16 == NULL) goto Error;
258    }
259
260    p -> nEntries  = nEntries;
261
262    // Initialize members if requested
263    if (Values != NULL && (nEntries > 0)) {
264
265        for (i=0; i < nEntries; i++)
266            p ->Table16[i] = Values[i];
267    }
268
269    // Initialize the segments stuff. The evaluator for each segment is located and a pointer to it
270    // is placed in advance to maximize performance.
271    if (Segments != NULL && (nSegments > 0)) {
272
273        _cmsParametricCurvesCollection *c;
274
275        p ->SegInterp = (cmsInterpParams**) _cmsCalloc(ContextID, nSegments, sizeof(cmsInterpParams*));
276        if (p ->SegInterp == NULL) goto Error;
277
278        for (i=0; i< nSegments; i++) {
279
280            // Type 0 is a special marker for table-based curves
281            if (Segments[i].Type == 0)
282                p ->SegInterp[i] = _cmsComputeInterpParams(ContextID, Segments[i].nGridPoints, 1, 1, NULL, CMS_LERP_FLAGS_FLOAT);
283
284            memmove(&p ->Segments[i], &Segments[i], sizeof(cmsCurveSegment));
285
286            if (Segments[i].Type == 0 && Segments[i].SampledPoints != NULL)
287                p ->Segments[i].SampledPoints = (cmsFloat32Number*) _cmsDupMem(ContextID, Segments[i].SampledPoints, sizeof(cmsFloat32Number) * Segments[i].nGridPoints);
288            else
289                p ->Segments[i].SampledPoints = NULL;
290
291
292            c = GetParametricCurveByType(ContextID, Segments[i].Type, NULL);
293            if (c != NULL)
294                    p ->Evals[i] = c ->Evaluator;
295        }
296    }
297
298    p ->InterpParams = _cmsComputeInterpParams(ContextID, p ->nEntries, 1, 1, p->Table16, CMS_LERP_FLAGS_16BITS);
299    if (p->InterpParams != NULL)
300        return p;
301
302Error:
303    if (p -> Segments) _cmsFree(ContextID, p ->Segments);
304    if (p -> Evals) _cmsFree(ContextID, p -> Evals);
305    if (p ->Table16) _cmsFree(ContextID, p ->Table16);
306    _cmsFree(ContextID, p);
307    return NULL;
308}
309
310
311// Parametric Fn using floating point
312static
313cmsFloat64Number DefaultEvalParametricFn(cmsInt32Number Type, const cmsFloat64Number Params[], cmsFloat64Number R)
314{
315    cmsFloat64Number e, Val, disc;
316
317    switch (Type) {
318
319   // X = Y ^ Gamma
320    case 1:
321        if (R < 0) {
322
323            if (fabs(Params[0] - 1.0) < MATRIX_DET_TOLERANCE)
324                Val = R;
325            else
326                Val = 0;
327        }
328        else
329            Val = pow(R, Params[0]);
330        break;
331
332    // Type 1 Reversed: X = Y ^1/gamma
333    case -1:
334         if (R < 0) {
335
336            if (fabs(Params[0] - 1.0) < MATRIX_DET_TOLERANCE)
337                Val = R;
338            else
339                Val = 0;
340        }
341        else
342            Val = pow(R, 1/Params[0]);
343        break;
344
345    // CIE 122-1966
346    // Y = (aX + b)^Gamma  | X >= -b/a
347    // Y = 0               | else
348    case 2:
349        disc = -Params[2] / Params[1];
350
351        if (R >= disc ) {
352
353            e = Params[1]*R + Params[2];
354
355            if (e > 0)
356                Val = pow(e, Params[0]);
357            else
358                Val = 0;
359        }
360        else
361            Val = 0;
362        break;
363
364     // Type 2 Reversed
365     // X = (Y ^1/g  - b) / a
366     case -2:
367         if (R < 0)
368             Val = 0;
369         else
370             Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1];
371
372         if (Val < 0)
373              Val = 0;
374         break;
375
376
377    // IEC 61966-3
378    // Y = (aX + b)^Gamma | X <= -b/a
379    // Y = c              | else
380    case 3:
381        disc = -Params[2] / Params[1];
382        if (disc < 0)
383            disc = 0;
384
385        if (R >= disc) {
386
387            e = Params[1]*R + Params[2];
388
389            if (e > 0)
390                Val = pow(e, Params[0]) + Params[3];
391            else
392                Val = 0;
393        }
394        else
395            Val = Params[3];
396        break;
397
398
399    // Type 3 reversed
400    // X=((Y-c)^1/g - b)/a      | (Y>=c)
401    // X=-b/a                   | (Y<c)
402    case -3:
403        if (R >= Params[3])  {
404
405            e = R - Params[3];
406
407            if (e > 0)
408                Val = (pow(e, 1/Params[0]) - Params[2]) / Params[1];
409            else
410                Val = 0;
411        }
412        else {
413            Val = -Params[2] / Params[1];
414        }
415        break;
416
417
418    // IEC 61966-2.1 (sRGB)
419    // Y = (aX + b)^Gamma | X >= d
420    // Y = cX             | X < d
421    case 4:
422        if (R >= Params[4]) {
423
424            e = Params[1]*R + Params[2];
425
426            if (e > 0)
427                Val = pow(e, Params[0]);
428            else
429                Val = 0;
430        }
431        else
432            Val = R * Params[3];
433        break;
434
435    // Type 4 reversed
436    // X=((Y^1/g-b)/a)    | Y >= (ad+b)^g
437    // X=Y/c              | Y< (ad+b)^g
438    case -4:
439        e = Params[1] * Params[4] + Params[2];
440        if (e < 0)
441            disc = 0;
442        else
443            disc = pow(e, Params[0]);
444
445        if (R >= disc) {
446
447            Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1];
448        }
449        else {
450            Val = R / Params[3];
451        }
452        break;
453
454
455    // Y = (aX + b)^Gamma + e | X >= d
456    // Y = cX + f             | X < d
457    case 5:
458        if (R >= Params[4]) {
459
460            e = Params[1]*R + Params[2];
461
462            if (e > 0)
463                Val = pow(e, Params[0]) + Params[5];
464            else
465                Val = Params[5];
466        }
467        else
468            Val = R*Params[3] + Params[6];
469        break;
470
471
472    // Reversed type 5
473    // X=((Y-e)1/g-b)/a   | Y >=(ad+b)^g+e), cd+f
474    // X=(Y-f)/c          | else
475    case -5:
476
477        disc = Params[3] * Params[4] + Params[6];
478        if (R >= disc) {
479
480            e = R - Params[5];
481            if (e < 0)
482                Val = 0;
483            else
484                Val = (pow(e, 1.0/Params[0]) - Params[2]) / Params[1];
485        }
486        else {
487            Val = (R - Params[6]) / Params[3];
488        }
489        break;
490
491
492    // Types 6,7,8 comes from segmented curves as described in ICCSpecRevision_02_11_06_Float.pdf
493    // Type 6 is basically identical to type 5 without d
494
495    // Y = (a * X + b) ^ Gamma + c
496    case 6:
497        e = Params[1]*R + Params[2];
498
499        if (e < 0)
500            Val = Params[3];
501        else
502            Val = pow(e, Params[0]) + Params[3];
503        break;
504
505    // ((Y - c) ^1/Gamma - b) / a
506    case -6:
507        e = R - Params[3];
508        if (e < 0)
509            Val = 0;
510        else
511        Val = (pow(e, 1.0/Params[0]) - Params[2]) / Params[1];
512        break;
513
514
515    // Y = a * log (b * X^Gamma + c) + d
516    case 7:
517
518       e = Params[2] * pow(R, Params[0]) + Params[3];
519       if (e <= 0)
520           Val = Params[4];
521       else
522           Val = Params[1]*log10(e) + Params[4];
523       break;
524
525    // (Y - d) / a = log(b * X ^Gamma + c)
526    // pow(10, (Y-d) / a) = b * X ^Gamma + c
527    // pow((pow(10, (Y-d) / a) - c) / b, 1/g) = X
528    case -7:
529       Val = pow((pow(10.0, (R-Params[4]) / Params[1]) - Params[3]) / Params[2], 1.0 / Params[0]);
530       break;
531
532
533   //Y = a * b^(c*X+d) + e
534   case 8:
535       Val = (Params[0] * pow(Params[1], Params[2] * R + Params[3]) + Params[4]);
536       break;
537
538
539   // Y = (log((y-e) / a) / log(b) - d ) / c
540   // a=0, b=1, c=2, d=3, e=4,
541   case -8:
542
543       disc = R - Params[4];
544       if (disc < 0) Val = 0;
545       else
546           Val = (log(disc / Params[0]) / log(Params[1]) - Params[3]) / Params[2];
547       break;
548
549   // S-Shaped: (1 - (1-x)^1/g)^1/g
550   case 108:
551      Val = pow(1.0 - pow(1 - R, 1/Params[0]), 1/Params[0]);
552      break;
553
554    // y = (1 - (1-x)^1/g)^1/g
555    // y^g = (1 - (1-x)^1/g)
556    // 1 - y^g = (1-x)^1/g
557    // (1 - y^g)^g = 1 - x
558    // 1 - (1 - y^g)^g
559    case -108:
560        Val = 1 - pow(1 - pow(R, Params[0]), Params[0]);
561        break;
562
563    default:
564        // Unsupported parametric curve. Should never reach here
565        return 0;
566    }
567
568    return Val;
569}
570
571// Evaluate a segmented funtion for a single value. Return -1 if no valid segment found .
572// If fn type is 0, perform an interpolation on the table
573static
574cmsFloat64Number EvalSegmentedFn(const cmsToneCurve *g, cmsFloat64Number R)
575{
576    int i;
577
578    for (i = g ->nSegments-1; i >= 0 ; --i) {
579
580        // Check for domain
581        if ((R > g ->Segments[i].x0) && (R <= g ->Segments[i].x1)) {
582
583            // Type == 0 means segment is sampled
584            if (g ->Segments[i].Type == 0) {
585
586                cmsFloat32Number R1 = (cmsFloat32Number) (R - g ->Segments[i].x0) / (g ->Segments[i].x1 - g ->Segments[i].x0);
587                cmsFloat32Number Out;
588
589                // Setup the table (TODO: clean that)
590                g ->SegInterp[i]-> Table = g ->Segments[i].SampledPoints;
591
592                g ->SegInterp[i] -> Interpolation.LerpFloat(&R1, &Out, g ->SegInterp[i]);
593
594                return Out;
595            }
596            else
597                return g ->Evals[i](g->Segments[i].Type, g ->Segments[i].Params, R);
598        }
599    }
600
601    return MINUS_INF;
602}
603
604// Access to estimated low-res table
605cmsUInt32Number CMSEXPORT cmsGetToneCurveEstimatedTableEntries(const cmsToneCurve* t)
606{
607    _cmsAssert(t != NULL);
608    return t ->nEntries;
609}
610
611const cmsUInt16Number* CMSEXPORT cmsGetToneCurveEstimatedTable(const cmsToneCurve* t)
612{
613    _cmsAssert(t != NULL);
614    return t ->Table16;
615}
616
617
618// Create an empty gamma curve, by using tables. This specifies only the limited-precision part, and leaves the
619// floating point description empty.
620cmsToneCurve* CMSEXPORT cmsBuildTabulatedToneCurve16(cmsContext ContextID, cmsInt32Number nEntries, const cmsUInt16Number Values[])
621{
622    return AllocateToneCurveStruct(ContextID, nEntries, 0, NULL, Values);
623}
624
625static
626int EntriesByGamma(cmsFloat64Number Gamma)
627{
628    if (fabs(Gamma - 1.0) < 0.001) return 2;
629    return 4096;
630}
631
632
633// Create a segmented gamma, fill the table
634cmsToneCurve* CMSEXPORT cmsBuildSegmentedToneCurve(cmsContext ContextID,
635                                                   cmsInt32Number nSegments, const cmsCurveSegment Segments[])
636{
637    int i;
638    cmsFloat64Number R, Val;
639    cmsToneCurve* g;
640    int nGridPoints = 4096;
641
642    _cmsAssert(Segments != NULL);
643
644    // Optimizatin for identity curves.
645    if (nSegments == 1 && Segments[0].Type == 1) {
646
647        nGridPoints = EntriesByGamma(Segments[0].Params[0]);
648    }
649
650    g = AllocateToneCurveStruct(ContextID, nGridPoints, nSegments, Segments, NULL);
651    if (g == NULL) return NULL;
652
653    // Once we have the floating point version, we can approximate a 16 bit table of 4096 entries
654    // for performance reasons. This table would normally not be used except on 8/16 bits transforms.
655    for (i=0; i < nGridPoints; i++) {
656
657        R   = (cmsFloat64Number) i / (nGridPoints-1);
658
659        Val = EvalSegmentedFn(g, R);
660
661        // Round and saturate
662        g ->Table16[i] = _cmsQuickSaturateWord(Val * 65535.0);
663    }
664
665    return g;
666}
667
668// Use a segmented curve to store the floating point table
669cmsToneCurve* CMSEXPORT cmsBuildTabulatedToneCurveFloat(cmsContext ContextID, cmsUInt32Number nEntries, const cmsFloat32Number values[])
670{
671    cmsCurveSegment Seg[3];
672
673    // A segmented tone curve should have function segments in the first and last positions
674    // Initialize segmented curve part up to 0 to constant value = samples[0]
675    Seg[0].x0 = MINUS_INF;
676    Seg[0].x1 = 0;
677    Seg[0].Type = 6;
678
679    Seg[0].Params[0] = 1;
680    Seg[0].Params[1] = 0;
681    Seg[0].Params[2] = 0;
682    Seg[0].Params[3] = values[0];
683    Seg[0].Params[4] = 0;
684
685    // From zero to 1
686    Seg[1].x0 = 0;
687    Seg[1].x1 = 1.0;
688    Seg[1].Type = 0;
689
690    Seg[1].nGridPoints = nEntries;
691    Seg[1].SampledPoints = (cmsFloat32Number*) values;
692
693    // Final segment is constant = lastsample
694    Seg[2].x0 = 1.0;
695    Seg[2].x1 = PLUS_INF;
696    Seg[2].Type = 6;
697
698    Seg[2].Params[0] = 1;
699    Seg[2].Params[1] = 0;
700    Seg[2].Params[2] = 0;
701    Seg[2].Params[3] = values[nEntries-1];
702    Seg[2].Params[4] = 0;
703
704
705    return cmsBuildSegmentedToneCurve(ContextID, 3, Seg);
706}
707
708// Parametric curves
709//
710// Parameters goes as: Curve, a, b, c, d, e, f
711// Type is the ICC type +1
712// if type is negative, then the curve is analyticaly inverted
713cmsToneCurve* CMSEXPORT cmsBuildParametricToneCurve(cmsContext ContextID, cmsInt32Number Type, const cmsFloat64Number Params[])
714{
715    cmsCurveSegment Seg0;
716    int Pos = 0;
717    cmsUInt32Number size;
718    _cmsParametricCurvesCollection* c = GetParametricCurveByType(ContextID, Type, &Pos);
719
720    _cmsAssert(Params != NULL);
721
722    if (c == NULL) {
723        cmsSignalError(ContextID, cmsERROR_UNKNOWN_EXTENSION, "Invalid parametric curve type %d", Type);
724        return NULL;
725    }
726
727    memset(&Seg0, 0, sizeof(Seg0));
728
729    Seg0.x0   = MINUS_INF;
730    Seg0.x1   = PLUS_INF;
731    Seg0.Type = Type;
732
733    size = c->ParameterCount[Pos] * sizeof(cmsFloat64Number);
734    memmove(Seg0.Params, Params, size);
735
736    return cmsBuildSegmentedToneCurve(ContextID, 1, &Seg0);
737}
738
739
740
741// Build a gamma table based on gamma constant
742cmsToneCurve* CMSEXPORT cmsBuildGamma(cmsContext ContextID, cmsFloat64Number Gamma)
743{
744    return cmsBuildParametricToneCurve(ContextID, 1, &Gamma);
745}
746
747
748// Free all memory taken by the gamma curve
749void CMSEXPORT cmsFreeToneCurve(cmsToneCurve* Curve)
750{
751    cmsContext ContextID;
752
753	// added by Xiaochuan Liu
754	// Curve->InterpParams may be null
755    if (Curve == NULL || Curve->InterpParams == NULL) return;
756
757    ContextID = Curve ->InterpParams->ContextID;
758
759    _cmsFreeInterpParams(Curve ->InterpParams);
760	Curve ->InterpParams = NULL;
761
762    if (Curve -> Table16)
763	{
764        _cmsFree(ContextID, Curve ->Table16);
765		Curve ->Table16 = NULL;
766	}
767
768    if (Curve ->Segments) {
769
770        cmsUInt32Number i;
771
772        for (i=0; i < Curve ->nSegments; i++) {
773
774            if (Curve ->Segments[i].SampledPoints) {
775                _cmsFree(ContextID, Curve ->Segments[i].SampledPoints);
776				Curve ->Segments[i].SampledPoints = NULL;
777            }
778
779            if (Curve ->SegInterp[i] != 0)
780			{
781                _cmsFreeInterpParams(Curve->SegInterp[i]);
782				Curve->SegInterp[i] = NULL;
783			}
784        }
785
786        _cmsFree(ContextID, Curve ->Segments);
787		Curve ->Segments = NULL;
788        _cmsFree(ContextID, Curve ->SegInterp);
789		Curve ->SegInterp = NULL;
790    }
791
792    if (Curve -> Evals)
793	{
794        _cmsFree(ContextID, Curve -> Evals);
795		Curve -> Evals = NULL;
796	}
797
798    if (Curve)
799	{
800		_cmsFree(ContextID, Curve);
801		Curve = NULL;
802	}
803}
804
805// Utility function, free 3 gamma tables
806void CMSEXPORT cmsFreeToneCurveTriple(cmsToneCurve* Curve[3])
807{
808
809    _cmsAssert(Curve != NULL);
810
811    if (Curve[0] != NULL) cmsFreeToneCurve(Curve[0]);
812    if (Curve[1] != NULL) cmsFreeToneCurve(Curve[1]);
813    if (Curve[2] != NULL) cmsFreeToneCurve(Curve[2]);
814
815    Curve[0] = Curve[1] = Curve[2] = NULL;
816}
817
818
819// Duplicate a gamma table
820cmsToneCurve* CMSEXPORT cmsDupToneCurve(const cmsToneCurve* In)
821{
822	// Xiaochuan Liu
823	// fix openpdf bug(mantis id:0055683, google id:360198)
824	// the function CurveSetElemTypeFree in cmslut.c also needs to check pointer
825    if (In == NULL || In ->InterpParams == NULL) return NULL;
826
827    return  AllocateToneCurveStruct(In ->InterpParams ->ContextID, In ->nEntries, In ->nSegments, In ->Segments, In ->Table16);
828}
829
830// Joins two curves for X and Y. Curves should be monotonic.
831// We want to get
832//
833//      y = Y^-1(X(t))
834//
835cmsToneCurve* CMSEXPORT cmsJoinToneCurve(cmsContext ContextID,
836                                      const cmsToneCurve* X,
837                                      const cmsToneCurve* Y, cmsUInt32Number nResultingPoints)
838{
839    cmsToneCurve* out = NULL;
840    cmsToneCurve* Yreversed = NULL;
841    cmsFloat32Number t, x;
842    cmsFloat32Number* Res = NULL;
843    cmsUInt32Number i;
844
845
846    _cmsAssert(X != NULL);
847    _cmsAssert(Y != NULL);
848
849    Yreversed = cmsReverseToneCurveEx(nResultingPoints, Y);
850    if (Yreversed == NULL) goto Error;
851
852    Res = (cmsFloat32Number*) _cmsCalloc(ContextID, nResultingPoints, sizeof(cmsFloat32Number));
853    if (Res == NULL) goto Error;
854
855    //Iterate
856    for (i=0; i <  nResultingPoints; i++) {
857
858        t = (cmsFloat32Number) i / (nResultingPoints-1);
859        x = cmsEvalToneCurveFloat(X,  t);
860        Res[i] = cmsEvalToneCurveFloat(Yreversed, x);
861    }
862
863    // Allocate space for output
864    out = cmsBuildTabulatedToneCurveFloat(ContextID, nResultingPoints, Res);
865
866Error:
867
868    if (Res != NULL) _cmsFree(ContextID, Res);
869    if (Yreversed != NULL) cmsFreeToneCurve(Yreversed);
870
871    return out;
872}
873
874
875
876// Get the surrounding nodes. This is tricky on non-monotonic tables
877static
878int GetInterval(cmsFloat64Number In, const cmsUInt16Number LutTable[], const struct _cms_interp_struc* p)
879{
880    int i;
881    int y0, y1;
882
883    // A 1 point table is not allowed
884    if (p -> Domain[0] < 1) return -1;
885
886    // Let's see if ascending or descending.
887    if (LutTable[0] < LutTable[p ->Domain[0]]) {
888
889        // Table is overall ascending
890        for (i=p->Domain[0]-1; i >=0; --i) {
891
892            y0 = LutTable[i];
893            y1 = LutTable[i+1];
894
895            if (y0 <= y1) { // Increasing
896                if (In >= y0 && In <= y1) return i;
897            }
898            else
899                if (y1 < y0) { // Decreasing
900                    if (In >= y1 && In <= y0) return i;
901                }
902        }
903    }
904    else {
905        // Table is overall descending
906        for (i=0; i < (int) p -> Domain[0]; i++) {
907
908            y0 = LutTable[i];
909            y1 = LutTable[i+1];
910
911            if (y0 <= y1) { // Increasing
912                if (In >= y0 && In <= y1) return i;
913            }
914            else
915                if (y1 < y0) { // Decreasing
916                    if (In >= y1 && In <= y0) return i;
917                }
918        }
919    }
920
921    return -1;
922}
923
924// Reverse a gamma table
925cmsToneCurve* CMSEXPORT cmsReverseToneCurveEx(cmsInt32Number nResultSamples, const cmsToneCurve* InCurve)
926{
927    cmsToneCurve *out;
928    cmsFloat64Number a = 0, b = 0, y, x1, y1, x2, y2;
929    int i, j;
930    int Ascending;
931
932    _cmsAssert(InCurve != NULL);
933
934    // Try to reverse it analytically whatever possible
935
936    if (InCurve ->nSegments == 1 && InCurve ->Segments[0].Type > 0 &&
937        /* InCurve -> Segments[0].Type <= 5 */
938        GetParametricCurveByType(InCurve ->InterpParams->ContextID, InCurve ->Segments[0].Type, NULL) != NULL) {
939
940        return cmsBuildParametricToneCurve(InCurve ->InterpParams->ContextID,
941                                       -(InCurve -> Segments[0].Type),
942                                       InCurve -> Segments[0].Params);
943    }
944
945    // Nope, reverse the table.
946    out = cmsBuildTabulatedToneCurve16(InCurve ->InterpParams->ContextID, nResultSamples, NULL);
947    if (out == NULL)
948        return NULL;
949
950    // We want to know if this is an ascending or descending table
951    Ascending = !cmsIsToneCurveDescending(InCurve);
952
953    // Iterate across Y axis
954    for (i=0; i <  nResultSamples; i++) {
955
956        y = (cmsFloat64Number) i * 65535.0 / (nResultSamples - 1);
957
958        // Find interval in which y is within.
959        j = GetInterval(y, InCurve->Table16, InCurve->InterpParams);
960        if (j >= 0) {
961
962
963            // Get limits of interval
964            x1 = InCurve ->Table16[j];
965            x2 = InCurve ->Table16[j+1];
966
967            y1 = (cmsFloat64Number) (j * 65535.0) / (InCurve ->nEntries - 1);
968            y2 = (cmsFloat64Number) ((j+1) * 65535.0 ) / (InCurve ->nEntries - 1);
969
970            // If collapsed, then use any
971            if (x1 == x2) {
972
973                out ->Table16[i] = _cmsQuickSaturateWord(Ascending ? y2 : y1);
974                continue;
975
976            } else {
977
978                // Interpolate
979                a = (y2 - y1) / (x2 - x1);
980                b = y2 - a * x2;
981            }
982        }
983
984        out ->Table16[i] = _cmsQuickSaturateWord(a* y + b);
985    }
986
987
988    return out;
989}
990
991// Reverse a gamma table
992cmsToneCurve* CMSEXPORT cmsReverseToneCurve(const cmsToneCurve* InGamma)
993{
994    _cmsAssert(InGamma != NULL);
995
996    return cmsReverseToneCurveEx(4096, InGamma);
997}
998
999// From: Eilers, P.H.C. (1994) Smoothing and interpolation with finite
1000// differences. in: Graphic Gems IV, Heckbert, P.S. (ed.), Academic press.
1001//
1002// Smoothing and interpolation with second differences.
1003//
1004//   Input:  weights (w), data (y): vector from 1 to m.
1005//   Input:  smoothing parameter (lambda), length (m).
1006//   Output: smoothed vector (z): vector from 1 to m.
1007
1008static
1009cmsBool smooth2(cmsContext ContextID, cmsFloat32Number w[], cmsFloat32Number y[], cmsFloat32Number z[], cmsFloat32Number lambda, int m)
1010{
1011    int i, i1, i2;
1012    cmsFloat32Number *c, *d, *e;
1013    cmsBool st;
1014
1015
1016    c = (cmsFloat32Number*) _cmsCalloc(ContextID, MAX_NODES_IN_CURVE, sizeof(cmsFloat32Number));
1017    d = (cmsFloat32Number*) _cmsCalloc(ContextID, MAX_NODES_IN_CURVE, sizeof(cmsFloat32Number));
1018    e = (cmsFloat32Number*) _cmsCalloc(ContextID, MAX_NODES_IN_CURVE, sizeof(cmsFloat32Number));
1019
1020    if (c != NULL && d != NULL && e != NULL) {
1021
1022
1023    d[1] = w[1] + lambda;
1024    c[1] = -2 * lambda / d[1];
1025    e[1] = lambda /d[1];
1026    z[1] = w[1] * y[1];
1027    d[2] = w[2] + 5 * lambda - d[1] * c[1] *  c[1];
1028    c[2] = (-4 * lambda - d[1] * c[1] * e[1]) / d[2];
1029    e[2] = lambda / d[2];
1030    z[2] = w[2] * y[2] - c[1] * z[1];
1031
1032    for (i = 3; i < m - 1; i++) {
1033        i1 = i - 1; i2 = i - 2;
1034        d[i]= w[i] + 6 * lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
1035        c[i] = (-4 * lambda -d[i1] * c[i1] * e[i1])/ d[i];
1036        e[i] = lambda / d[i];
1037        z[i] = w[i] * y[i] - c[i1] * z[i1] - e[i2] * z[i2];
1038    }
1039
1040    i1 = m - 2; i2 = m - 3;
1041
1042    d[m - 1] = w[m - 1] + 5 * lambda -c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
1043    c[m - 1] = (-2 * lambda - d[i1] * c[i1] * e[i1]) / d[m - 1];
1044    z[m - 1] = w[m - 1] * y[m - 1] - c[i1] * z[i1] - e[i2] * z[i2];
1045    i1 = m - 1; i2 = m - 2;
1046
1047    d[m] = w[m] + lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
1048    z[m] = (w[m] * y[m] - c[i1] * z[i1] - e[i2] * z[i2]) / d[m];
1049    z[m - 1] = z[m - 1] / d[m - 1] - c[m - 1] * z[m];
1050
1051    for (i = m - 2; 1<= i; i--)
1052        z[i] = z[i] / d[i] - c[i] * z[i + 1] - e[i] * z[i + 2];
1053
1054      st = TRUE;
1055    }
1056    else st = FALSE;
1057
1058    if (c != NULL) _cmsFree(ContextID, c);
1059    if (d != NULL) _cmsFree(ContextID, d);
1060    if (e != NULL) _cmsFree(ContextID, e);
1061
1062    return st;
1063}
1064
1065// Smooths a curve sampled at regular intervals.
1066cmsBool  CMSEXPORT cmsSmoothToneCurve(cmsToneCurve* Tab, cmsFloat64Number lambda)
1067{
1068    cmsFloat32Number w[MAX_NODES_IN_CURVE], y[MAX_NODES_IN_CURVE], z[MAX_NODES_IN_CURVE];
1069    int i, nItems, Zeros, Poles;
1070
1071    if (Tab == NULL) return FALSE;
1072
1073    if (cmsIsToneCurveLinear(Tab)) return TRUE; // Nothing to do
1074
1075    nItems = Tab -> nEntries;
1076
1077    if (nItems >= MAX_NODES_IN_CURVE) {
1078        cmsSignalError(Tab ->InterpParams->ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: too many points.");
1079        return FALSE;
1080    }
1081
1082    memset(w, 0, nItems * sizeof(cmsFloat32Number));
1083    memset(y, 0, nItems * sizeof(cmsFloat32Number));
1084    memset(z, 0, nItems * sizeof(cmsFloat32Number));
1085
1086    for (i=0; i < nItems; i++)
1087    {
1088        y[i+1] = (cmsFloat32Number) Tab -> Table16[i];
1089        w[i+1] = 1.0;
1090    }
1091
1092    if (!smooth2(Tab ->InterpParams->ContextID, w, y, z, (cmsFloat32Number) lambda, nItems)) return FALSE;
1093
1094    // Do some reality - checking...
1095    Zeros = Poles = 0;
1096    for (i=nItems; i > 1; --i) {
1097
1098        if (z[i] == 0.) Zeros++;
1099        if (z[i] >= 65535.) Poles++;
1100        if (z[i] < z[i-1]) {
1101            cmsSignalError(Tab ->InterpParams->ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Non-Monotonic.");
1102            return FALSE;
1103        }
1104    }
1105
1106    if (Zeros > (nItems / 3)) {
1107        cmsSignalError(Tab ->InterpParams->ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Degenerated, mostly zeros.");
1108        return FALSE;
1109    }
1110    if (Poles > (nItems / 3)) {
1111        cmsSignalError(Tab ->InterpParams->ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Degenerated, mostly poles.");
1112        return FALSE;
1113    }
1114
1115    // Seems ok
1116    for (i=0; i < nItems; i++) {
1117
1118        // Clamp to cmsUInt16Number
1119        Tab -> Table16[i] = _cmsQuickSaturateWord(z[i+1]);
1120    }
1121
1122    return TRUE;
1123}
1124
1125// Is a table linear? Do not use parametric since we cannot guarantee some weird parameters resulting
1126// in a linear table. This way assures it is linear in 12 bits, which should be enought in most cases.
1127cmsBool CMSEXPORT cmsIsToneCurveLinear(const cmsToneCurve* Curve)
1128{
1129    cmsUInt32Number i;
1130    int diff;
1131
1132    _cmsAssert(Curve != NULL);
1133
1134    for (i=0; i < Curve ->nEntries; i++) {
1135
1136        diff = abs((int) Curve->Table16[i] - (int) _cmsQuantizeVal(i, Curve ->nEntries));
1137        if (diff > 0x0f)
1138            return FALSE;
1139    }
1140
1141    return TRUE;
1142}
1143
1144// Same, but for monotonicity
1145cmsBool  CMSEXPORT cmsIsToneCurveMonotonic(const cmsToneCurve* t)
1146{
1147    int n;
1148    int i, last;
1149    cmsBool lDescending;
1150
1151    _cmsAssert(t != NULL);
1152
1153    // Degenerated curves are monotonic? Ok, let's pass them
1154    n = t ->nEntries;
1155    if (n < 2) return TRUE;
1156
1157    // Curve direction
1158    lDescending = cmsIsToneCurveDescending(t);
1159
1160    if (lDescending) {
1161
1162        last = t ->Table16[0];
1163
1164        for (i = 1; i < n; i++) {
1165
1166            if (t ->Table16[i] - last > 2) // We allow some ripple
1167                return FALSE;
1168            else
1169                last = t ->Table16[i];
1170
1171        }
1172    }
1173    else {
1174
1175        last = t ->Table16[n-1];
1176
1177        for (i = n-2; i >= 0; --i) {
1178
1179            if (t ->Table16[i] - last > 2)
1180                return FALSE;
1181            else
1182                last = t ->Table16[i];
1183
1184        }
1185    }
1186
1187    return TRUE;
1188}
1189
1190// Same, but for descending tables
1191cmsBool  CMSEXPORT cmsIsToneCurveDescending(const cmsToneCurve* t)
1192{
1193    _cmsAssert(t != NULL);
1194
1195    return t ->Table16[0] > t ->Table16[t ->nEntries-1];
1196}
1197
1198
1199// Another info fn: is out gamma table multisegment?
1200cmsBool  CMSEXPORT cmsIsToneCurveMultisegment(const cmsToneCurve* t)
1201{
1202    _cmsAssert(t != NULL);
1203
1204    return t -> nSegments > 1;
1205}
1206
1207cmsInt32Number  CMSEXPORT cmsGetToneCurveParametricType(const cmsToneCurve* t)
1208{
1209    _cmsAssert(t != NULL);
1210
1211    if (t -> nSegments != 1) return 0;
1212    return t ->Segments[0].Type;
1213}
1214
1215// We need accuracy this time
1216cmsFloat32Number CMSEXPORT cmsEvalToneCurveFloat(const cmsToneCurve* Curve, cmsFloat32Number v)
1217{
1218    _cmsAssert(Curve != NULL);
1219
1220    // Check for 16 bits table. If so, this is a limited-precision tone curve
1221    if (Curve ->nSegments == 0) {
1222
1223        cmsUInt16Number In, Out;
1224
1225        In = (cmsUInt16Number) _cmsQuickSaturateWord(v * 65535.0);
1226        Out = cmsEvalToneCurve16(Curve, In);
1227
1228        return (cmsFloat32Number) (Out / 65535.0);
1229    }
1230
1231    return (cmsFloat32Number) EvalSegmentedFn(Curve, v);
1232}
1233
1234// We need xput over here
1235cmsUInt16Number CMSEXPORT cmsEvalToneCurve16(const cmsToneCurve* Curve, cmsUInt16Number v)
1236{
1237    cmsUInt16Number out;
1238
1239    _cmsAssert(Curve != NULL);
1240
1241    Curve ->InterpParams ->Interpolation.Lerp16(&v, &out, Curve ->InterpParams);
1242    return out;
1243}
1244
1245
1246// Least squares fitting.
1247// A mathematical procedure for finding the best-fitting curve to a given set of points by
1248// minimizing the sum of the squares of the offsets ("the residuals") of the points from the curve.
1249// The sum of the squares of the offsets is used instead of the offset absolute values because
1250// this allows the residuals to be treated as a continuous differentiable quantity.
1251//
1252// y = f(x) = x ^ g
1253//
1254// R  = (yi - (xi^g))
1255// R2 = (yi - (xi^g))2
1256// SUM R2 = SUM (yi - (xi^g))2
1257//
1258// dR2/dg = -2 SUM x^g log(x)(y - x^g)
1259// solving for dR2/dg = 0
1260//
1261// g = 1/n * SUM(log(y) / log(x))
1262
1263cmsFloat64Number CMSEXPORT cmsEstimateGamma(const cmsToneCurve* t, cmsFloat64Number Precision)
1264{
1265    cmsFloat64Number gamma, sum, sum2;
1266    cmsFloat64Number n, x, y, Std;
1267    cmsUInt32Number i;
1268
1269    _cmsAssert(t != NULL);
1270
1271    sum = sum2 = n = 0;
1272
1273    // Excluding endpoints
1274    for (i=1; i < (MAX_NODES_IN_CURVE-1); i++) {
1275
1276        x = (cmsFloat64Number) i / (MAX_NODES_IN_CURVE-1);
1277        y = (cmsFloat64Number) cmsEvalToneCurveFloat(t, (cmsFloat32Number) x);
1278
1279        // Avoid 7% on lower part to prevent
1280        // artifacts due to linear ramps
1281
1282        if (y > 0. && y < 1. && x > 0.07) {
1283
1284            gamma = log(y) / log(x);
1285            sum  += gamma;
1286            sum2 += gamma * gamma;
1287            n++;
1288        }
1289    }
1290
1291    // Take a look on SD to see if gamma isn't exponential at all
1292    Std = sqrt((n * sum2 - sum * sum) / (n*(n-1)));
1293
1294    if (Std > Precision)
1295        return -1.0;
1296
1297    return (sum / n);   // The mean
1298}
1299