1//---------------------------------------------------------------------------------
2//
3//  Little Color Management System
4//  Copyright (c) 1998-2012 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
30// Allocates an empty multi profile element
31cmsStage* CMSEXPORT _cmsStageAllocPlaceholder(cmsContext ContextID,
32                                cmsStageSignature Type,
33                                cmsUInt32Number InputChannels,
34                                cmsUInt32Number OutputChannels,
35                                _cmsStageEvalFn     EvalPtr,
36                                _cmsStageDupElemFn  DupElemPtr,
37                                _cmsStageFreeElemFn FreePtr,
38                                void*             Data)
39{
40    cmsStage* ph = (cmsStage*) _cmsMallocZero(ContextID, sizeof(cmsStage));
41
42    if (ph == NULL) return NULL;
43
44
45    ph ->ContextID = ContextID;
46
47    ph ->Type       = Type;
48    ph ->Implements = Type;   // By default, no clue on what is implementing
49
50    ph ->InputChannels  = InputChannels;
51    ph ->OutputChannels = OutputChannels;
52    ph ->EvalPtr        = EvalPtr;
53    ph ->DupElemPtr     = DupElemPtr;
54    ph ->FreePtr        = FreePtr;
55    ph ->Data           = Data;
56
57    return ph;
58}
59
60
61static
62void EvaluateIdentity(const cmsFloat32Number In[],
63                            cmsFloat32Number Out[],
64                      const cmsStage *mpe)
65{
66    memmove(Out, In, mpe ->InputChannels * sizeof(cmsFloat32Number));
67}
68
69
70cmsStage* CMSEXPORT cmsStageAllocIdentity(cmsContext ContextID, cmsUInt32Number nChannels)
71{
72    return _cmsStageAllocPlaceholder(ContextID,
73                                   cmsSigIdentityElemType,
74                                   nChannels, nChannels,
75                                   EvaluateIdentity,
76                                   NULL,
77                                   NULL,
78                                   NULL);
79 }
80
81// Conversion functions. From floating point to 16 bits
82static
83void FromFloatTo16(const cmsFloat32Number In[], cmsUInt16Number Out[], cmsUInt32Number n)
84{
85    cmsUInt32Number i;
86
87    for (i=0; i < n; i++) {
88        Out[i] = _cmsQuickSaturateWord(In[i] * 65535.0);
89    }
90}
91
92// From 16 bits to floating point
93static
94void From16ToFloat(const cmsUInt16Number In[], cmsFloat32Number Out[], cmsUInt32Number n)
95{
96    cmsUInt32Number i;
97
98    for (i=0; i < n; i++) {
99        Out[i] = (cmsFloat32Number) In[i] / 65535.0F;
100    }
101}
102
103
104// This function is quite useful to analyze the structure of a LUT and retrieve the MPE elements
105// that conform the LUT. It should be called with the LUT, the number of expected elements and
106// then a list of expected types followed with a list of cmsFloat64Number pointers to MPE elements. If
107// the function founds a match with current pipeline, it fills the pointers and returns TRUE
108// if not, returns FALSE without touching anything. Setting pointers to NULL does bypass
109// the storage process.
110cmsBool  CMSEXPORT cmsPipelineCheckAndRetreiveStages(const cmsPipeline* Lut, cmsUInt32Number n, ...)
111{
112    va_list args;
113    cmsUInt32Number i;
114    cmsStage* mpe;
115    cmsStageSignature Type;
116    void** ElemPtr;
117
118    // Make sure same number of elements
119    if (cmsPipelineStageCount(Lut) != n) return FALSE;
120
121    va_start(args, n);
122
123    // Iterate across asked types
124    mpe = Lut ->Elements;
125    for (i=0; i < n; i++) {
126
127        // Get asked type
128        Type  = (cmsStageSignature)va_arg(args, cmsStageSignature);
129        if (mpe ->Type != Type) {
130
131            va_end(args);       // Mismatch. We are done.
132            return FALSE;
133        }
134        mpe = mpe ->Next;
135    }
136
137    // Found a combination, fill pointers if not NULL
138    mpe = Lut ->Elements;
139    for (i=0; i < n; i++) {
140
141        ElemPtr = va_arg(args, void**);
142        if (ElemPtr != NULL)
143            *ElemPtr = mpe;
144
145        mpe = mpe ->Next;
146    }
147
148    va_end(args);
149    return TRUE;
150}
151
152// Below there are implementations for several types of elements. Each type may be implemented by a
153// evaluation function, a duplication function, a function to free resources and a constructor.
154
155// *************************************************************************************************
156// Type cmsSigCurveSetElemType (curves)
157// *************************************************************************************************
158
159cmsToneCurve** _cmsStageGetPtrToCurveSet(const cmsStage* mpe)
160{
161    _cmsStageToneCurvesData* Data = (_cmsStageToneCurvesData*) mpe ->Data;
162
163    return Data ->TheCurves;
164}
165
166static
167void EvaluateCurves(const cmsFloat32Number In[],
168                    cmsFloat32Number Out[],
169                    const cmsStage *mpe)
170{
171    _cmsStageToneCurvesData* Data;
172    cmsUInt32Number i;
173
174    _cmsAssert(mpe != NULL);
175
176    Data = (_cmsStageToneCurvesData*) mpe ->Data;
177    if (Data == NULL) return;
178
179    if (Data ->TheCurves == NULL) return;
180
181    for (i=0; i < Data ->nCurves; i++) {
182        Out[i] = cmsEvalToneCurveFloat(Data ->TheCurves[i], In[i]);
183    }
184}
185
186static
187void CurveSetElemTypeFree(cmsStage* mpe)
188{
189    _cmsStageToneCurvesData* Data;
190    cmsUInt32Number i;
191
192    _cmsAssert(mpe != NULL);
193
194    Data = (_cmsStageToneCurvesData*) mpe ->Data;
195    if (Data == NULL) return;
196
197    if (Data ->TheCurves != NULL) {
198        for (i=0; i < Data ->nCurves; i++) {
199            if (Data ->TheCurves[i] != NULL)
200                cmsFreeToneCurve(Data ->TheCurves[i]);
201        }
202    }
203    _cmsFree(mpe ->ContextID, Data ->TheCurves);
204    _cmsFree(mpe ->ContextID, Data);
205}
206
207
208static
209void* CurveSetDup(cmsStage* mpe)
210{
211    _cmsStageToneCurvesData* Data = (_cmsStageToneCurvesData*) mpe ->Data;
212    _cmsStageToneCurvesData* NewElem;
213    cmsUInt32Number i;
214
215    NewElem = (_cmsStageToneCurvesData*) _cmsMallocZero(mpe ->ContextID, sizeof(_cmsStageToneCurvesData));
216    if (NewElem == NULL) return NULL;
217
218    NewElem ->nCurves   = Data ->nCurves;
219    NewElem ->TheCurves = (cmsToneCurve**) _cmsCalloc(mpe ->ContextID, NewElem ->nCurves, sizeof(cmsToneCurve*));
220
221    if (NewElem ->TheCurves == NULL) goto Error;
222
223    for (i=0; i < NewElem ->nCurves; i++) {
224
225        // Duplicate each curve. It may fail.
226        NewElem ->TheCurves[i] = cmsDupToneCurve(Data ->TheCurves[i]);
227        if (NewElem ->TheCurves[i] == NULL) goto Error;
228
229
230    }
231    return (void*) NewElem;
232
233Error:
234
235    if (NewElem ->TheCurves != NULL) {
236        for (i=0; i < NewElem ->nCurves; i++) {
237            if (NewElem ->TheCurves[i])
238                cmsFreeToneCurve(NewElem ->TheCurves[i]);
239        }
240    }
241    _cmsFree(mpe ->ContextID, NewElem ->TheCurves);
242    _cmsFree(mpe ->ContextID, NewElem);
243    return NULL;
244}
245
246
247// Curves == NULL forces identity curves
248cmsStage* CMSEXPORT cmsStageAllocToneCurves(cmsContext ContextID, cmsUInt32Number nChannels, cmsToneCurve* const Curves[])
249{
250    cmsUInt32Number i;
251    _cmsStageToneCurvesData* NewElem;
252    cmsStage* NewMPE;
253
254
255    NewMPE = _cmsStageAllocPlaceholder(ContextID, cmsSigCurveSetElemType, nChannels, nChannels,
256                                     EvaluateCurves, CurveSetDup, CurveSetElemTypeFree, NULL );
257    if (NewMPE == NULL) return NULL;
258
259    NewElem = (_cmsStageToneCurvesData*) _cmsMallocZero(ContextID, sizeof(_cmsStageToneCurvesData));
260    if (NewElem == NULL) {
261        cmsStageFree(NewMPE);
262        return NULL;
263    }
264
265    NewMPE ->Data  = (void*) NewElem;
266
267    NewElem ->nCurves   = nChannels;
268    NewElem ->TheCurves = (cmsToneCurve**) _cmsCalloc(ContextID, nChannels, sizeof(cmsToneCurve*));
269    if (NewElem ->TheCurves == NULL) {
270        cmsStageFree(NewMPE);
271        return NULL;
272    }
273
274    for (i=0; i < nChannels; i++) {
275
276        if (Curves == NULL) {
277            NewElem ->TheCurves[i] = cmsBuildGamma(ContextID, 1.0);
278        }
279        else {
280            NewElem ->TheCurves[i] = cmsDupToneCurve(Curves[i]);
281        }
282
283        if (NewElem ->TheCurves[i] == NULL) {
284            cmsStageFree(NewMPE);
285            return NULL;
286        }
287
288    }
289
290   return NewMPE;
291}
292
293
294// Create a bunch of identity curves
295cmsStage* _cmsStageAllocIdentityCurves(cmsContext ContextID, int nChannels)
296{
297    cmsStage* mpe = cmsStageAllocToneCurves(ContextID, nChannels, NULL);
298
299    if (mpe == NULL) return NULL;
300    mpe ->Implements = cmsSigIdentityElemType;
301    return mpe;
302}
303
304
305// *************************************************************************************************
306// Type cmsSigMatrixElemType (Matrices)
307// *************************************************************************************************
308
309
310// Special care should be taken here because precision loss. A temporary cmsFloat64Number buffer is being used
311static
312void EvaluateMatrix(const cmsFloat32Number In[],
313                    cmsFloat32Number Out[],
314                    const cmsStage *mpe)
315{
316    cmsUInt32Number i, j;
317    _cmsStageMatrixData* Data = (_cmsStageMatrixData*) mpe ->Data;
318    cmsFloat64Number Tmp;
319
320    // Input is already in 0..1.0 notation
321    for (i=0; i < mpe ->OutputChannels; i++) {
322
323        Tmp = 0;
324        for (j=0; j < mpe->InputChannels; j++) {
325            Tmp += In[j] * Data->Double[i*mpe->InputChannels + j];
326        }
327
328        if (Data ->Offset != NULL)
329            Tmp += Data->Offset[i];
330
331        Out[i] = (cmsFloat32Number) Tmp;
332    }
333
334
335    // Output in 0..1.0 domain
336}
337
338
339// Duplicate a yet-existing matrix element
340static
341void* MatrixElemDup(cmsStage* mpe)
342{
343    _cmsStageMatrixData* Data = (_cmsStageMatrixData*) mpe ->Data;
344    _cmsStageMatrixData* NewElem;
345    cmsUInt32Number sz;
346
347    NewElem = (_cmsStageMatrixData*) _cmsMallocZero(mpe ->ContextID, sizeof(_cmsStageMatrixData));
348    if (NewElem == NULL) return NULL;
349
350    sz = mpe ->InputChannels * mpe ->OutputChannels;
351
352    NewElem ->Double = (cmsFloat64Number*) _cmsDupMem(mpe ->ContextID, Data ->Double, sz * sizeof(cmsFloat64Number)) ;
353
354    if (Data ->Offset)
355        NewElem ->Offset = (cmsFloat64Number*) _cmsDupMem(mpe ->ContextID,
356                                                Data ->Offset, mpe -> OutputChannels * sizeof(cmsFloat64Number)) ;
357
358    return (void*) NewElem;
359}
360
361
362static
363void MatrixElemTypeFree(cmsStage* mpe)
364{
365    _cmsStageMatrixData* Data = (_cmsStageMatrixData*) mpe ->Data;
366    if (Data == NULL)
367        return;
368    if (Data ->Double)
369        _cmsFree(mpe ->ContextID, Data ->Double);
370
371    if (Data ->Offset)
372        _cmsFree(mpe ->ContextID, Data ->Offset);
373
374    _cmsFree(mpe ->ContextID, mpe ->Data);
375}
376
377
378
379cmsStage*  CMSEXPORT cmsStageAllocMatrix(cmsContext ContextID, cmsUInt32Number Rows, cmsUInt32Number Cols,
380                                     const cmsFloat64Number* Matrix, const cmsFloat64Number* Offset)
381{
382    cmsUInt32Number i, n;
383    _cmsStageMatrixData* NewElem;
384    cmsStage* NewMPE;
385
386    n = Rows * Cols;
387
388    // Check for overflow
389    if (n == 0) return NULL;
390    if (n >= UINT_MAX / Cols) return NULL;
391    if (n >= UINT_MAX / Rows) return NULL;
392    if (n < Rows || n < Cols) return NULL;
393
394    NewMPE = _cmsStageAllocPlaceholder(ContextID, cmsSigMatrixElemType, Cols, Rows,
395                                     EvaluateMatrix, MatrixElemDup, MatrixElemTypeFree, NULL );
396    if (NewMPE == NULL) return NULL;
397
398
399    NewElem = (_cmsStageMatrixData*) _cmsMallocZero(ContextID, sizeof(_cmsStageMatrixData));
400    if (NewElem == NULL) return NULL;
401
402
403    NewElem ->Double = (cmsFloat64Number*) _cmsCalloc(ContextID, n, sizeof(cmsFloat64Number));
404
405    if (NewElem->Double == NULL) {
406        MatrixElemTypeFree(NewMPE);
407        return NULL;
408    }
409
410    for (i=0; i < n; i++) {
411        NewElem ->Double[i] = Matrix[i];
412    }
413
414
415    if (Offset != NULL) {
416
417        NewElem ->Offset = (cmsFloat64Number*) _cmsCalloc(ContextID, Rows, sizeof(cmsFloat64Number));
418        if (NewElem->Offset == NULL) {
419           MatrixElemTypeFree(NewMPE);
420           return NULL;
421        }
422
423        for (i=0; i < Rows; i++) {
424                NewElem ->Offset[i] = Offset[i];
425        }
426
427    }
428
429    NewMPE ->Data  = (void*) NewElem;
430    return NewMPE;
431}
432
433
434// *************************************************************************************************
435// Type cmsSigCLutElemType
436// *************************************************************************************************
437
438
439// Evaluate in true floating point
440static
441void EvaluateCLUTfloat(const cmsFloat32Number In[], cmsFloat32Number Out[], const cmsStage *mpe)
442{
443    _cmsStageCLutData* Data = (_cmsStageCLutData*) mpe ->Data;
444
445    Data -> Params ->Interpolation.LerpFloat(In, Out, Data->Params);
446}
447
448
449// Convert to 16 bits, evaluate, and back to floating point
450static
451void EvaluateCLUTfloatIn16(const cmsFloat32Number In[], cmsFloat32Number Out[], const cmsStage *mpe)
452{
453    _cmsStageCLutData* Data = (_cmsStageCLutData*) mpe ->Data;
454    cmsUInt16Number In16[MAX_STAGE_CHANNELS], Out16[MAX_STAGE_CHANNELS];
455
456    _cmsAssert(mpe ->InputChannels  <= MAX_STAGE_CHANNELS);
457    _cmsAssert(mpe ->OutputChannels <= MAX_STAGE_CHANNELS);
458
459    FromFloatTo16(In, In16, mpe ->InputChannels);
460    Data -> Params ->Interpolation.Lerp16(In16, Out16, Data->Params);
461    From16ToFloat(Out16, Out,  mpe ->OutputChannels);
462}
463
464
465// Given an hypercube of b dimensions, with Dims[] number of nodes by dimension, calculate the total amount of nodes
466static
467cmsUInt32Number CubeSize(const cmsUInt32Number Dims[], cmsUInt32Number b)
468{
469    cmsUInt32Number rv, dim;
470
471    _cmsAssert(Dims != NULL);
472
473    for (rv = 1; b > 0; b--) {
474
475        dim = Dims[b-1];
476        if (dim == 0) return 0;  // Error
477
478        rv *= dim;
479
480        // Check for overflow
481        if (rv > UINT_MAX / dim) return 0;
482    }
483
484    return rv;
485}
486
487static
488void* CLUTElemDup(cmsStage* mpe)
489{
490    _cmsStageCLutData* Data = (_cmsStageCLutData*) mpe ->Data;
491    _cmsStageCLutData* NewElem;
492
493
494    NewElem = (_cmsStageCLutData*) _cmsMallocZero(mpe ->ContextID, sizeof(_cmsStageCLutData));
495    if (NewElem == NULL) return NULL;
496
497    NewElem ->nEntries       = Data ->nEntries;
498    NewElem ->HasFloatValues = Data ->HasFloatValues;
499
500    if (Data ->Tab.T) {
501
502        if (Data ->HasFloatValues) {
503            NewElem ->Tab.TFloat = (cmsFloat32Number*) _cmsDupMem(mpe ->ContextID, Data ->Tab.TFloat, Data ->nEntries * sizeof (cmsFloat32Number));
504            if (NewElem ->Tab.TFloat == NULL)
505                goto Error;
506        } else {
507            NewElem ->Tab.T = (cmsUInt16Number*) _cmsDupMem(mpe ->ContextID, Data ->Tab.T, Data ->nEntries * sizeof (cmsUInt16Number));
508            if (NewElem ->Tab.TFloat == NULL)
509                goto Error;
510        }
511    }
512
513    NewElem ->Params   = _cmsComputeInterpParamsEx(mpe ->ContextID,
514                                                   Data ->Params ->nSamples,
515                                                   Data ->Params ->nInputs,
516                                                   Data ->Params ->nOutputs,
517                                                   NewElem ->Tab.T,
518                                                   Data ->Params ->dwFlags);
519    if (NewElem->Params != NULL)
520        return (void*) NewElem;
521 Error:
522    if (NewElem->Tab.T)
523        // This works for both types
524        _cmsFree(mpe ->ContextID, NewElem -> Tab.T);
525    _cmsFree(mpe ->ContextID, NewElem);
526    return NULL;
527}
528
529
530static
531void CLutElemTypeFree(cmsStage* mpe)
532{
533
534    _cmsStageCLutData* Data = (_cmsStageCLutData*) mpe ->Data;
535
536    // Already empty
537    if (Data == NULL) return;
538
539    // This works for both types
540    if (Data -> Tab.T)
541        _cmsFree(mpe ->ContextID, Data -> Tab.T);
542
543    _cmsFreeInterpParams(Data ->Params);
544    _cmsFree(mpe ->ContextID, mpe ->Data);
545}
546
547
548// Allocates a 16-bit multidimensional CLUT. This is evaluated at 16-bit precision. Table may have different
549// granularity on each dimension.
550cmsStage* CMSEXPORT cmsStageAllocCLut16bitGranular(cmsContext ContextID,
551                                         const cmsUInt32Number clutPoints[],
552                                         cmsUInt32Number inputChan,
553                                         cmsUInt32Number outputChan,
554                                         const cmsUInt16Number* Table)
555{
556    cmsUInt32Number i, n;
557    _cmsStageCLutData* NewElem;
558    cmsStage* NewMPE;
559
560    _cmsAssert(clutPoints != NULL);
561
562    if (inputChan > MAX_INPUT_DIMENSIONS) {
563        cmsSignalError(ContextID, cmsERROR_RANGE, "Too many input channels (%d channels, max=%d)", inputChan, MAX_INPUT_DIMENSIONS);
564        return NULL;
565    }
566
567    NewMPE = _cmsStageAllocPlaceholder(ContextID, cmsSigCLutElemType, inputChan, outputChan,
568                                     EvaluateCLUTfloatIn16, CLUTElemDup, CLutElemTypeFree, NULL );
569
570    if (NewMPE == NULL) return NULL;
571
572    NewElem = (_cmsStageCLutData*) _cmsMallocZero(ContextID, sizeof(_cmsStageCLutData));
573    if (NewElem == NULL) {
574        cmsStageFree(NewMPE);
575        return NULL;
576    }
577
578    NewMPE ->Data  = (void*) NewElem;
579
580    NewElem -> nEntries = n = outputChan * CubeSize(clutPoints, inputChan);
581    NewElem -> HasFloatValues = FALSE;
582
583    if (n == 0) {
584        cmsStageFree(NewMPE);
585        return NULL;
586    }
587
588
589    NewElem ->Tab.T  = (cmsUInt16Number*) _cmsCalloc(ContextID, n, sizeof(cmsUInt16Number));
590    if (NewElem ->Tab.T == NULL) {
591        cmsStageFree(NewMPE);
592        return NULL;
593    }
594
595    if (Table != NULL) {
596        for (i=0; i < n; i++) {
597            NewElem ->Tab.T[i] = Table[i];
598        }
599    }
600
601    NewElem ->Params = _cmsComputeInterpParamsEx(ContextID, clutPoints, inputChan, outputChan, NewElem ->Tab.T, CMS_LERP_FLAGS_16BITS);
602    if (NewElem ->Params == NULL) {
603        cmsStageFree(NewMPE);
604        return NULL;
605    }
606
607    return NewMPE;
608}
609
610cmsStage* CMSEXPORT cmsStageAllocCLut16bit(cmsContext ContextID,
611                                    cmsUInt32Number nGridPoints,
612                                    cmsUInt32Number inputChan,
613                                    cmsUInt32Number outputChan,
614                                    const cmsUInt16Number* Table)
615{
616    cmsUInt32Number Dimensions[MAX_INPUT_DIMENSIONS];
617    int i;
618
619   // Our resulting LUT would be same gridpoints on all dimensions
620    for (i=0; i < MAX_INPUT_DIMENSIONS; i++)
621        Dimensions[i] = nGridPoints;
622
623    return cmsStageAllocCLut16bitGranular(ContextID, Dimensions, inputChan, outputChan, Table);
624}
625
626
627cmsStage* CMSEXPORT cmsStageAllocCLutFloat(cmsContext ContextID,
628                                       cmsUInt32Number nGridPoints,
629                                       cmsUInt32Number inputChan,
630                                       cmsUInt32Number outputChan,
631                                       const cmsFloat32Number* Table)
632{
633   cmsUInt32Number Dimensions[MAX_INPUT_DIMENSIONS];
634   int i;
635
636    // Our resulting LUT would be same gridpoints on all dimensions
637    for (i=0; i < MAX_INPUT_DIMENSIONS; i++)
638        Dimensions[i] = nGridPoints;
639
640    return cmsStageAllocCLutFloatGranular(ContextID, Dimensions, inputChan, outputChan, Table);
641}
642
643
644
645cmsStage* CMSEXPORT cmsStageAllocCLutFloatGranular(cmsContext ContextID, const cmsUInt32Number clutPoints[], cmsUInt32Number inputChan, cmsUInt32Number outputChan, const cmsFloat32Number* Table)
646{
647    cmsUInt32Number i, n;
648    _cmsStageCLutData* NewElem;
649    cmsStage* NewMPE;
650
651    _cmsAssert(clutPoints != NULL);
652
653    if (inputChan > MAX_INPUT_DIMENSIONS) {
654        cmsSignalError(ContextID, cmsERROR_RANGE, "Too many input channels (%d channels, max=%d)", inputChan, MAX_INPUT_DIMENSIONS);
655        return NULL;
656    }
657
658    NewMPE = _cmsStageAllocPlaceholder(ContextID, cmsSigCLutElemType, inputChan, outputChan,
659                                             EvaluateCLUTfloat, CLUTElemDup, CLutElemTypeFree, NULL);
660    if (NewMPE == NULL) return NULL;
661
662
663    NewElem = (_cmsStageCLutData*) _cmsMallocZero(ContextID, sizeof(_cmsStageCLutData));
664    if (NewElem == NULL) {
665        cmsStageFree(NewMPE);
666        return NULL;
667    }
668
669    NewMPE ->Data  = (void*) NewElem;
670
671    // There is a potential integer overflow on conputing n and nEntries.
672    NewElem -> nEntries = n = outputChan * CubeSize(clutPoints, inputChan);
673    NewElem -> HasFloatValues = TRUE;
674
675    if (n == 0) {
676        cmsStageFree(NewMPE);
677        return NULL;
678    }
679
680    NewElem ->Tab.TFloat  = (cmsFloat32Number*) _cmsCalloc(ContextID, n, sizeof(cmsFloat32Number));
681    if (NewElem ->Tab.TFloat == NULL) {
682        cmsStageFree(NewMPE);
683        return NULL;
684    }
685
686    if (Table != NULL) {
687        for (i=0; i < n; i++) {
688            NewElem ->Tab.TFloat[i] = Table[i];
689        }
690    }
691
692    NewElem ->Params = _cmsComputeInterpParamsEx(ContextID, clutPoints,  inputChan, outputChan, NewElem ->Tab.TFloat, CMS_LERP_FLAGS_FLOAT);
693    if (NewElem ->Params == NULL) {
694        cmsStageFree(NewMPE);
695        return NULL;
696    }
697
698    return NewMPE;
699}
700
701
702static
703int IdentitySampler(register const cmsUInt16Number In[], register cmsUInt16Number Out[], register void * Cargo)
704{
705    int nChan = *(int*) Cargo;
706    int i;
707
708    for (i=0; i < nChan; i++)
709        Out[i] = In[i];
710
711    return 1;
712}
713
714// Creates an MPE that just copies input to output
715cmsStage* _cmsStageAllocIdentityCLut(cmsContext ContextID, int nChan)
716{
717    cmsUInt32Number Dimensions[MAX_INPUT_DIMENSIONS];
718    cmsStage* mpe ;
719    int i;
720
721    for (i=0; i < MAX_INPUT_DIMENSIONS; i++)
722        Dimensions[i] = 2;
723
724    mpe = cmsStageAllocCLut16bitGranular(ContextID, Dimensions, nChan, nChan, NULL);
725    if (mpe == NULL) return NULL;
726
727    if (!cmsStageSampleCLut16bit(mpe, IdentitySampler, &nChan, 0)) {
728        cmsStageFree(mpe);
729        return NULL;
730    }
731
732    mpe ->Implements = cmsSigIdentityElemType;
733    return mpe;
734}
735
736
737
738// Quantize a value 0 <= i < MaxSamples to 0..0xffff
739cmsUInt16Number _cmsQuantizeVal(cmsFloat64Number i, int MaxSamples)
740{
741    cmsFloat64Number x;
742
743    x = ((cmsFloat64Number) i * 65535.) / (cmsFloat64Number) (MaxSamples - 1);
744    return _cmsQuickSaturateWord(x);
745}
746
747
748// This routine does a sweep on whole input space, and calls its callback
749// function on knots. returns TRUE if all ok, FALSE otherwise.
750cmsBool CMSEXPORT cmsStageSampleCLut16bit(cmsStage* mpe, cmsSAMPLER16 Sampler, void * Cargo, cmsUInt32Number dwFlags)
751{
752    int i, t, nTotalPoints, index, rest;
753    int nInputs, nOutputs;
754    cmsUInt32Number* nSamples;
755    cmsUInt16Number In[MAX_INPUT_DIMENSIONS+1], Out[MAX_STAGE_CHANNELS];
756    _cmsStageCLutData* clut;
757
758    if (mpe == NULL) return FALSE;
759
760    clut = (_cmsStageCLutData*) mpe->Data;
761
762    if (clut == NULL) return FALSE;
763
764    nSamples = clut->Params ->nSamples;
765    nInputs  = clut->Params ->nInputs;
766    nOutputs = clut->Params ->nOutputs;
767
768    if (nInputs <= 0) return FALSE;
769    if (nOutputs <= 0) return FALSE;
770    if (nInputs > MAX_INPUT_DIMENSIONS) return FALSE;
771    if (nOutputs >= MAX_STAGE_CHANNELS) return FALSE;
772
773    nTotalPoints = CubeSize(nSamples, nInputs);
774    if (nTotalPoints == 0) return FALSE;
775
776    index = 0;
777    for (i = 0; i < nTotalPoints; i++) {
778
779        rest = i;
780        for (t = nInputs-1; t >=0; --t) {
781
782            cmsUInt32Number  Colorant = rest % nSamples[t];
783
784            rest /= nSamples[t];
785
786            In[t] = _cmsQuantizeVal(Colorant, nSamples[t]);
787        }
788
789        if (clut ->Tab.T != NULL) {
790            for (t=0; t < nOutputs; t++)
791                Out[t] = clut->Tab.T[index + t];
792        }
793
794        if (!Sampler(In, Out, Cargo))
795            return FALSE;
796
797        if (!(dwFlags & SAMPLER_INSPECT)) {
798
799            if (clut ->Tab.T != NULL) {
800                for (t=0; t < nOutputs; t++)
801                    clut->Tab.T[index + t] = Out[t];
802            }
803        }
804
805        index += nOutputs;
806    }
807
808    return TRUE;
809}
810
811// Same as anterior, but for floting point
812cmsBool CMSEXPORT cmsStageSampleCLutFloat(cmsStage* mpe, cmsSAMPLERFLOAT Sampler, void * Cargo, cmsUInt32Number dwFlags)
813{
814    int i, t, nTotalPoints, index, rest;
815    int nInputs, nOutputs;
816    cmsUInt32Number* nSamples;
817    cmsFloat32Number In[MAX_INPUT_DIMENSIONS+1], Out[MAX_STAGE_CHANNELS];
818    _cmsStageCLutData* clut = (_cmsStageCLutData*) mpe->Data;
819
820    nSamples = clut->Params ->nSamples;
821    nInputs  = clut->Params ->nInputs;
822    nOutputs = clut->Params ->nOutputs;
823
824    if (nInputs <= 0) return FALSE;
825    if (nOutputs <= 0) return FALSE;
826    if (nInputs  > MAX_INPUT_DIMENSIONS) return FALSE;
827    if (nOutputs >= MAX_STAGE_CHANNELS) return FALSE;
828
829    nTotalPoints = CubeSize(nSamples, nInputs);
830    if (nTotalPoints == 0) return FALSE;
831
832    index = 0;
833    for (i = 0; i < nTotalPoints; i++) {
834
835        rest = i;
836        for (t = nInputs-1; t >=0; --t) {
837
838            cmsUInt32Number  Colorant = rest % nSamples[t];
839
840            rest /= nSamples[t];
841
842            In[t] =  (cmsFloat32Number) (_cmsQuantizeVal(Colorant, nSamples[t]) / 65535.0);
843        }
844
845        if (clut ->Tab.TFloat != NULL) {
846            for (t=0; t < nOutputs; t++)
847                Out[t] = clut->Tab.TFloat[index + t];
848        }
849
850        if (!Sampler(In, Out, Cargo))
851            return FALSE;
852
853        if (!(dwFlags & SAMPLER_INSPECT)) {
854
855            if (clut ->Tab.TFloat != NULL) {
856                for (t=0; t < nOutputs; t++)
857                    clut->Tab.TFloat[index + t] = Out[t];
858            }
859        }
860
861        index += nOutputs;
862    }
863
864    return TRUE;
865}
866
867
868
869// This routine does a sweep on whole input space, and calls its callback
870// function on knots. returns TRUE if all ok, FALSE otherwise.
871cmsBool CMSEXPORT cmsSliceSpace16(cmsUInt32Number nInputs, const cmsUInt32Number clutPoints[],
872                                         cmsSAMPLER16 Sampler, void * Cargo)
873{
874    int i, t, nTotalPoints, rest;
875    cmsUInt16Number In[cmsMAXCHANNELS];
876
877    if (nInputs >= cmsMAXCHANNELS) return FALSE;
878
879    nTotalPoints = CubeSize(clutPoints, nInputs);
880    if (nTotalPoints == 0) return FALSE;
881
882    for (i = 0; i < nTotalPoints; i++) {
883
884        rest = i;
885        for (t = nInputs-1; t >=0; --t) {
886
887            cmsUInt32Number  Colorant = rest % clutPoints[t];
888
889            rest /= clutPoints[t];
890            In[t] = _cmsQuantizeVal(Colorant, clutPoints[t]);
891
892        }
893
894        if (!Sampler(In, NULL, Cargo))
895            return FALSE;
896    }
897
898    return TRUE;
899}
900
901cmsInt32Number CMSEXPORT cmsSliceSpaceFloat(cmsUInt32Number nInputs, const cmsUInt32Number clutPoints[],
902                                            cmsSAMPLERFLOAT Sampler, void * Cargo)
903{
904    int i, t, nTotalPoints, rest;
905    cmsFloat32Number In[cmsMAXCHANNELS];
906
907    if (nInputs >= cmsMAXCHANNELS) return FALSE;
908
909    nTotalPoints = CubeSize(clutPoints, nInputs);
910    if (nTotalPoints == 0) return FALSE;
911
912    for (i = 0; i < nTotalPoints; i++) {
913
914        rest = i;
915        for (t = nInputs-1; t >=0; --t) {
916
917            cmsUInt32Number  Colorant = rest % clutPoints[t];
918
919            rest /= clutPoints[t];
920            In[t] =  (cmsFloat32Number) (_cmsQuantizeVal(Colorant, clutPoints[t]) / 65535.0);
921
922        }
923
924        if (!Sampler(In, NULL, Cargo))
925            return FALSE;
926    }
927
928    return TRUE;
929}
930
931// ********************************************************************************
932// Type cmsSigLab2XYZElemType
933// ********************************************************************************
934
935
936static
937void EvaluateLab2XYZ(const cmsFloat32Number In[],
938                     cmsFloat32Number Out[],
939                     const cmsStage *mpe)
940{
941    cmsCIELab Lab;
942    cmsCIEXYZ XYZ;
943    const cmsFloat64Number XYZadj = MAX_ENCODEABLE_XYZ;
944
945    // V4 rules
946    Lab.L = In[0] * 100.0;
947    Lab.a = In[1] * 255.0 - 128.0;
948    Lab.b = In[2] * 255.0 - 128.0;
949
950    cmsLab2XYZ(NULL, &XYZ, &Lab);
951
952    // From XYZ, range 0..19997 to 0..1.0, note that 1.99997 comes from 0xffff
953    // encoded as 1.15 fixed point, so 1 + (32767.0 / 32768.0)
954
955    Out[0] = (cmsFloat32Number) ((cmsFloat64Number) XYZ.X / XYZadj);
956    Out[1] = (cmsFloat32Number) ((cmsFloat64Number) XYZ.Y / XYZadj);
957    Out[2] = (cmsFloat32Number) ((cmsFloat64Number) XYZ.Z / XYZadj);
958    return;
959
960    cmsUNUSED_PARAMETER(mpe);
961}
962
963
964// No dup or free routines needed, as the structure has no pointers in it.
965cmsStage* _cmsStageAllocLab2XYZ(cmsContext ContextID)
966{
967    return _cmsStageAllocPlaceholder(ContextID, cmsSigLab2XYZElemType, 3, 3, EvaluateLab2XYZ, NULL, NULL, NULL);
968}
969
970// ********************************************************************************
971
972// v2 L=100 is supposed to be placed on 0xFF00. There is no reasonable
973// number of gridpoints that would make exact match. However, a prelinearization
974// of 258 entries, would map 0xFF00 exactly on entry 257, and this is good to avoid scum dot.
975// Almost all what we need but unfortunately, the rest of entries should be scaled by
976// (255*257/256) and this is not exact.
977
978cmsStage* _cmsStageAllocLabV2ToV4curves(cmsContext ContextID)
979{
980    cmsStage* mpe;
981    cmsToneCurve* LabTable[3];
982    int i, j;
983
984    LabTable[0] = cmsBuildTabulatedToneCurve16(ContextID, 258, NULL);
985    LabTable[1] = cmsBuildTabulatedToneCurve16(ContextID, 258, NULL);
986    LabTable[2] = cmsBuildTabulatedToneCurve16(ContextID, 258, NULL);
987
988    for (j=0; j < 3; j++) {
989
990        if (LabTable[j] == NULL) {
991            cmsFreeToneCurveTriple(LabTable);
992            return NULL;
993        }
994
995        // We need to map * (0xffff / 0xff00), thats same as (257 / 256)
996        // So we can use 258-entry tables to do the trick (i / 257) * (255 * 257) * (257 / 256);
997        for (i=0; i < 257; i++)  {
998
999            LabTable[j]->Table16[i] = (cmsUInt16Number) ((i * 0xffff + 0x80) >> 8);
1000        }
1001
1002        LabTable[j] ->Table16[257] = 0xffff;
1003    }
1004
1005    mpe = cmsStageAllocToneCurves(ContextID, 3, LabTable);
1006    cmsFreeToneCurveTriple(LabTable);
1007
1008    if (mpe == NULL) return NULL;
1009    mpe ->Implements = cmsSigLabV2toV4;
1010    return mpe;
1011}
1012
1013// ********************************************************************************
1014
1015// Matrix-based conversion, which is more accurate, but slower and cannot properly be saved in devicelink profiles
1016cmsStage* _cmsStageAllocLabV2ToV4(cmsContext ContextID)
1017{
1018    static const cmsFloat64Number V2ToV4[] = { 65535.0/65280.0, 0, 0,
1019                                     0, 65535.0/65280.0, 0,
1020                                     0, 0, 65535.0/65280.0
1021                                     };
1022
1023    cmsStage *mpe = cmsStageAllocMatrix(ContextID, 3, 3, V2ToV4, NULL);
1024
1025    if (mpe == NULL) return mpe;
1026    mpe ->Implements = cmsSigLabV2toV4;
1027    return mpe;
1028}
1029
1030
1031// Reverse direction
1032cmsStage* _cmsStageAllocLabV4ToV2(cmsContext ContextID)
1033{
1034    static const cmsFloat64Number V4ToV2[] = { 65280.0/65535.0, 0, 0,
1035                                     0, 65280.0/65535.0, 0,
1036                                     0, 0, 65280.0/65535.0
1037                                     };
1038
1039     cmsStage *mpe = cmsStageAllocMatrix(ContextID, 3, 3, V4ToV2, NULL);
1040
1041    if (mpe == NULL) return mpe;
1042    mpe ->Implements = cmsSigLabV4toV2;
1043    return mpe;
1044}
1045
1046
1047// To Lab to float. Note that the MPE gives numbers in normal Lab range
1048// and we need 0..1.0 range for the formatters
1049// L* : 0...100 => 0...1.0  (L* / 100)
1050// ab* : -128..+127 to 0..1  ((ab* + 128) / 255)
1051
1052cmsStage* _cmsStageNormalizeFromLabFloat(cmsContext ContextID)
1053{
1054    static const cmsFloat64Number a1[] = {
1055        1.0/100.0, 0, 0,
1056        0, 1.0/255.0, 0,
1057        0, 0, 1.0/255.0
1058    };
1059
1060    static const cmsFloat64Number o1[] = {
1061        0,
1062        128.0/255.0,
1063        128.0/255.0
1064    };
1065
1066    cmsStage *mpe = cmsStageAllocMatrix(ContextID, 3, 3, a1, o1);
1067
1068    if (mpe == NULL) return mpe;
1069    mpe ->Implements = cmsSigLab2FloatPCS;
1070    return mpe;
1071}
1072
1073// Fom XYZ to floating point PCS
1074cmsStage* _cmsStageNormalizeFromXyzFloat(cmsContext ContextID)
1075{
1076#define n (32768.0/65535.0)
1077    static const cmsFloat64Number a1[] = {
1078        n, 0, 0,
1079        0, n, 0,
1080        0, 0, n
1081    };
1082#undef n
1083
1084    cmsStage *mpe =  cmsStageAllocMatrix(ContextID, 3, 3, a1, NULL);
1085
1086    if (mpe == NULL) return mpe;
1087    mpe ->Implements = cmsSigXYZ2FloatPCS;
1088    return mpe;
1089}
1090
1091cmsStage* _cmsStageNormalizeToLabFloat(cmsContext ContextID)
1092{
1093    static const cmsFloat64Number a1[] = {
1094        100.0, 0, 0,
1095        0, 255.0, 0,
1096        0, 0, 255.0
1097    };
1098
1099    static const cmsFloat64Number o1[] = {
1100        0,
1101        -128.0,
1102        -128.0
1103    };
1104
1105    cmsStage *mpe =  cmsStageAllocMatrix(ContextID, 3, 3, a1, o1);
1106    if (mpe == NULL) return mpe;
1107    mpe ->Implements = cmsSigFloatPCS2Lab;
1108    return mpe;
1109}
1110
1111cmsStage* _cmsStageNormalizeToXyzFloat(cmsContext ContextID)
1112{
1113#define n (65535.0/32768.0)
1114
1115    static const cmsFloat64Number a1[] = {
1116        n, 0, 0,
1117        0, n, 0,
1118        0, 0, n
1119    };
1120#undef n
1121
1122    cmsStage *mpe = cmsStageAllocMatrix(ContextID, 3, 3, a1, NULL);
1123    if (mpe == NULL) return mpe;
1124    mpe ->Implements = cmsSigFloatPCS2XYZ;
1125    return mpe;
1126}
1127
1128
1129
1130// ********************************************************************************
1131// Type cmsSigXYZ2LabElemType
1132// ********************************************************************************
1133
1134static
1135void EvaluateXYZ2Lab(const cmsFloat32Number In[], cmsFloat32Number Out[], const cmsStage *mpe)
1136{
1137    cmsCIELab Lab;
1138    cmsCIEXYZ XYZ;
1139    const cmsFloat64Number XYZadj = MAX_ENCODEABLE_XYZ;
1140
1141    // From 0..1.0 to XYZ
1142
1143    XYZ.X = In[0] * XYZadj;
1144    XYZ.Y = In[1] * XYZadj;
1145    XYZ.Z = In[2] * XYZadj;
1146
1147    cmsXYZ2Lab(NULL, &Lab, &XYZ);
1148
1149    // From V4 Lab to 0..1.0
1150
1151    Out[0] = (cmsFloat32Number) (Lab.L / 100.0);
1152    Out[1] = (cmsFloat32Number) ((Lab.a + 128.0) / 255.0);
1153    Out[2] = (cmsFloat32Number) ((Lab.b + 128.0) / 255.0);
1154    return;
1155
1156    cmsUNUSED_PARAMETER(mpe);
1157}
1158
1159cmsStage* _cmsStageAllocXYZ2Lab(cmsContext ContextID)
1160{
1161    return _cmsStageAllocPlaceholder(ContextID, cmsSigXYZ2LabElemType, 3, 3, EvaluateXYZ2Lab, NULL, NULL, NULL);
1162
1163}
1164
1165// ********************************************************************************
1166
1167// For v4, S-Shaped curves are placed in a/b axis to increase resolution near gray
1168
1169cmsStage* _cmsStageAllocLabPrelin(cmsContext ContextID)
1170{
1171    cmsToneCurve* LabTable[3];
1172    cmsFloat64Number Params[1] =  {2.4} ;
1173
1174    LabTable[0] = cmsBuildGamma(ContextID, 1.0);
1175    LabTable[1] = cmsBuildParametricToneCurve(ContextID, 108, Params);
1176    LabTable[2] = cmsBuildParametricToneCurve(ContextID, 108, Params);
1177
1178    return cmsStageAllocToneCurves(ContextID, 3, LabTable);
1179}
1180
1181
1182// Free a single MPE
1183void CMSEXPORT cmsStageFree(cmsStage* mpe)
1184{
1185    if (mpe ->FreePtr)
1186        mpe ->FreePtr(mpe);
1187
1188    _cmsFree(mpe ->ContextID, mpe);
1189}
1190
1191
1192cmsUInt32Number  CMSEXPORT cmsStageInputChannels(const cmsStage* mpe)
1193{
1194    return mpe ->InputChannels;
1195}
1196
1197cmsUInt32Number  CMSEXPORT cmsStageOutputChannels(const cmsStage* mpe)
1198{
1199    return mpe ->OutputChannels;
1200}
1201
1202cmsStageSignature CMSEXPORT cmsStageType(const cmsStage* mpe)
1203{
1204    return mpe -> Type;
1205}
1206
1207void* CMSEXPORT cmsStageData(const cmsStage* mpe)
1208{
1209    return mpe -> Data;
1210}
1211
1212cmsStage*  CMSEXPORT cmsStageNext(const cmsStage* mpe)
1213{
1214    return mpe -> Next;
1215}
1216
1217
1218// Duplicates an MPE
1219cmsStage* CMSEXPORT cmsStageDup(cmsStage* mpe)
1220{
1221    cmsStage* NewMPE;
1222
1223    if (mpe == NULL) return NULL;
1224    NewMPE = _cmsStageAllocPlaceholder(mpe ->ContextID,
1225                                     mpe ->Type,
1226                                     mpe ->InputChannels,
1227                                     mpe ->OutputChannels,
1228                                     mpe ->EvalPtr,
1229                                     mpe ->DupElemPtr,
1230                                     mpe ->FreePtr,
1231                                     NULL);
1232    if (NewMPE == NULL) return NULL;
1233
1234    NewMPE ->Implements = mpe ->Implements;
1235
1236    if (mpe ->DupElemPtr) {
1237
1238        NewMPE ->Data = mpe ->DupElemPtr(mpe);
1239
1240        if (NewMPE->Data == NULL) {
1241
1242            cmsStageFree(NewMPE);
1243            return NULL;
1244        }
1245
1246    } else {
1247
1248        NewMPE ->Data       = NULL;
1249    }
1250
1251    return NewMPE;
1252}
1253
1254
1255// ***********************************************************************************************************
1256
1257// This function sets up the channel count
1258static
1259cmsBool BlessLUT(cmsPipeline* lut)
1260{
1261    // We can set the input/ouput channels only if we have elements.
1262    if (lut ->Elements != NULL) {
1263
1264        cmsStage* prev;
1265        cmsStage* next;
1266        cmsStage* First;
1267        cmsStage* Last;
1268
1269        First  = cmsPipelineGetPtrToFirstStage(lut);
1270        Last   = cmsPipelineGetPtrToLastStage(lut);
1271
1272        if (First == NULL || Last == NULL) return FALSE;
1273
1274        lut->InputChannels = First->InputChannels;
1275        lut->OutputChannels = Last->OutputChannels;
1276
1277        // Check chain consistency
1278        prev = First;
1279        next = prev->Next;
1280
1281        while (next != NULL)
1282        {
1283            if (next->InputChannels != prev->OutputChannels)
1284                return FALSE;
1285
1286            next = next->Next;
1287            prev = prev->Next;
1288        }
1289    }
1290    return TRUE;
1291}
1292
1293
1294// Default to evaluate the LUT on 16 bit-basis. Precision is retained.
1295static
1296void _LUTeval16(register const cmsUInt16Number In[], register cmsUInt16Number Out[],  register const void* D)
1297{
1298    cmsPipeline* lut = (cmsPipeline*) D;
1299    cmsStage *mpe;
1300    cmsFloat32Number Storage[2][MAX_STAGE_CHANNELS] = {0.0f};
1301    int Phase = 0, NextPhase;
1302
1303    From16ToFloat(In, &Storage[Phase][0], lut ->InputChannels);
1304
1305    for (mpe = lut ->Elements;
1306         mpe != NULL;
1307         mpe = mpe ->Next) {
1308
1309             NextPhase = Phase ^ 1;
1310             mpe ->EvalPtr(&Storage[Phase][0], &Storage[NextPhase][0], mpe);
1311             Phase = NextPhase;
1312    }
1313
1314
1315    FromFloatTo16(&Storage[Phase][0], Out, lut ->OutputChannels);
1316}
1317
1318
1319
1320// Does evaluate the LUT on cmsFloat32Number-basis.
1321static
1322void _LUTevalFloat(register const cmsFloat32Number In[], register cmsFloat32Number Out[], const void* D)
1323{
1324    cmsPipeline* lut = (cmsPipeline*) D;
1325    cmsStage *mpe;
1326    cmsFloat32Number Storage[2][MAX_STAGE_CHANNELS] = {0.0f};
1327    int Phase = 0, NextPhase;
1328
1329    memmove(&Storage[Phase][0], In, lut ->InputChannels  * sizeof(cmsFloat32Number));
1330
1331    for (mpe = lut ->Elements;
1332         mpe != NULL;
1333         mpe = mpe ->Next) {
1334
1335              NextPhase = Phase ^ 1;
1336              mpe ->EvalPtr(&Storage[Phase][0], &Storage[NextPhase][0], mpe);
1337              Phase = NextPhase;
1338    }
1339
1340    memmove(Out, &Storage[Phase][0], lut ->OutputChannels * sizeof(cmsFloat32Number));
1341}
1342
1343
1344
1345
1346// LUT Creation & Destruction
1347
1348cmsPipeline* CMSEXPORT cmsPipelineAlloc(cmsContext ContextID, cmsUInt32Number InputChannels, cmsUInt32Number OutputChannels)
1349{
1350       cmsPipeline* NewLUT;
1351
1352       // A value of zero in channels is allowed as placeholder
1353       if (InputChannels >= cmsMAXCHANNELS ||
1354           OutputChannels >= cmsMAXCHANNELS) return NULL;
1355
1356       NewLUT = (cmsPipeline*) _cmsMallocZero(ContextID, sizeof(cmsPipeline));
1357       if (NewLUT == NULL) return NULL;
1358
1359
1360       NewLUT -> InputChannels  = InputChannels;
1361       NewLUT -> OutputChannels = OutputChannels;
1362
1363       NewLUT ->Eval16Fn    = _LUTeval16;
1364       NewLUT ->EvalFloatFn = _LUTevalFloat;
1365       NewLUT ->DupDataFn   = NULL;
1366       NewLUT ->FreeDataFn  = NULL;
1367       NewLUT ->Data        = NewLUT;
1368       NewLUT ->ContextID   = ContextID;
1369
1370       if (!BlessLUT(NewLUT))
1371       {
1372           _cmsFree(ContextID, NewLUT);
1373           return NULL;
1374       }
1375
1376       return NewLUT;
1377}
1378
1379cmsContext CMSEXPORT cmsGetPipelineContextID(const cmsPipeline* lut)
1380{
1381    _cmsAssert(lut != NULL);
1382    return lut ->ContextID;
1383}
1384
1385cmsUInt32Number CMSEXPORT cmsPipelineInputChannels(const cmsPipeline* lut)
1386{
1387    _cmsAssert(lut != NULL);
1388    return lut ->InputChannels;
1389}
1390
1391cmsUInt32Number CMSEXPORT cmsPipelineOutputChannels(const cmsPipeline* lut)
1392{
1393    _cmsAssert(lut != NULL);
1394    return lut ->OutputChannels;
1395}
1396
1397// Free a profile elements LUT
1398void CMSEXPORT cmsPipelineFree(cmsPipeline* lut)
1399{
1400    cmsStage *mpe, *Next;
1401
1402    if (lut == NULL) return;
1403
1404    for (mpe = lut ->Elements;
1405        mpe != NULL;
1406        mpe = Next) {
1407
1408            Next = mpe ->Next;
1409            cmsStageFree(mpe);
1410    }
1411
1412    if (lut ->FreeDataFn) lut ->FreeDataFn(lut ->ContextID, lut ->Data);
1413
1414    _cmsFree(lut ->ContextID, lut);
1415}
1416
1417
1418// Default to evaluate the LUT on 16 bit-basis.
1419void CMSEXPORT cmsPipelineEval16(const cmsUInt16Number In[], cmsUInt16Number Out[],  const cmsPipeline* lut)
1420{
1421    _cmsAssert(lut != NULL);
1422    lut ->Eval16Fn(In, Out, lut->Data);
1423}
1424
1425
1426// Does evaluate the LUT on cmsFloat32Number-basis.
1427void CMSEXPORT cmsPipelineEvalFloat(const cmsFloat32Number In[], cmsFloat32Number Out[], const cmsPipeline* lut)
1428{
1429    _cmsAssert(lut != NULL);
1430    lut ->EvalFloatFn(In, Out, lut);
1431}
1432
1433
1434
1435// Duplicates a LUT
1436cmsPipeline* CMSEXPORT cmsPipelineDup(const cmsPipeline* lut)
1437{
1438    cmsPipeline* NewLUT;
1439    cmsStage *NewMPE, *Anterior = NULL, *mpe;
1440    cmsBool  First = TRUE;
1441
1442    if (lut == NULL) return NULL;
1443
1444    NewLUT = cmsPipelineAlloc(lut ->ContextID, lut ->InputChannels, lut ->OutputChannels);
1445    if (NewLUT == NULL) return NULL;
1446
1447    for (mpe = lut ->Elements;
1448         mpe != NULL;
1449         mpe = mpe ->Next) {
1450
1451             NewMPE = cmsStageDup(mpe);
1452
1453             if (NewMPE == NULL) {
1454                 cmsPipelineFree(NewLUT);
1455                 return NULL;
1456             }
1457
1458             if (First) {
1459                 NewLUT ->Elements = NewMPE;
1460                 First = FALSE;
1461             }
1462             else {
1463                Anterior ->Next = NewMPE;
1464             }
1465
1466            Anterior = NewMPE;
1467    }
1468
1469    NewLUT ->Eval16Fn    = lut ->Eval16Fn;
1470    NewLUT ->EvalFloatFn = lut ->EvalFloatFn;
1471    NewLUT ->DupDataFn   = lut ->DupDataFn;
1472    NewLUT ->FreeDataFn  = lut ->FreeDataFn;
1473
1474    if (NewLUT ->DupDataFn != NULL)
1475        NewLUT ->Data = NewLUT ->DupDataFn(lut ->ContextID, lut->Data);
1476
1477
1478    NewLUT ->SaveAs8Bits    = lut ->SaveAs8Bits;
1479
1480    if (!BlessLUT(NewLUT))
1481    {
1482        _cmsFree(lut->ContextID, NewLUT);
1483        return NULL;
1484    }
1485
1486    return NewLUT;
1487}
1488
1489
1490int CMSEXPORT cmsPipelineInsertStage(cmsPipeline* lut, cmsStageLoc loc, cmsStage* mpe)
1491{
1492    cmsStage* Anterior = NULL, *pt;
1493
1494    if (lut == NULL || mpe == NULL)
1495        return FALSE;
1496
1497    switch (loc) {
1498
1499        case cmsAT_BEGIN:
1500            mpe ->Next = lut ->Elements;
1501            lut ->Elements = mpe;
1502            break;
1503
1504        case cmsAT_END:
1505
1506            if (lut ->Elements == NULL)
1507                lut ->Elements = mpe;
1508            else {
1509
1510                for (pt = lut ->Elements;
1511                     pt != NULL;
1512                     pt = pt -> Next) Anterior = pt;
1513
1514                Anterior ->Next = mpe;
1515                mpe ->Next = NULL;
1516            }
1517            break;
1518        default:;
1519            return FALSE;
1520    }
1521
1522    return BlessLUT(lut);
1523}
1524
1525// Unlink an element and return the pointer to it
1526void CMSEXPORT cmsPipelineUnlinkStage(cmsPipeline* lut, cmsStageLoc loc, cmsStage** mpe)
1527{
1528    cmsStage *Anterior, *pt, *Last;
1529    cmsStage *Unlinked = NULL;
1530
1531
1532    // If empty LUT, there is nothing to remove
1533    if (lut ->Elements == NULL) {
1534        if (mpe) *mpe = NULL;
1535        return;
1536    }
1537
1538    // On depending on the strategy...
1539    switch (loc) {
1540
1541        case cmsAT_BEGIN:
1542            {
1543                cmsStage* elem = lut ->Elements;
1544
1545                lut ->Elements = elem -> Next;
1546                elem ->Next = NULL;
1547                Unlinked = elem;
1548
1549            }
1550            break;
1551
1552        case cmsAT_END:
1553            Anterior = Last = NULL;
1554            for (pt = lut ->Elements;
1555                pt != NULL;
1556                pt = pt -> Next) {
1557                    Anterior = Last;
1558                    Last = pt;
1559            }
1560
1561            Unlinked = Last;  // Next already points to NULL
1562
1563            // Truncate the chain
1564            if (Anterior)
1565                Anterior ->Next = NULL;
1566            else
1567                lut ->Elements = NULL;
1568            break;
1569        default:;
1570    }
1571
1572    if (mpe)
1573        *mpe = Unlinked;
1574    else
1575        cmsStageFree(Unlinked);
1576
1577    // May fail, but we ignore it
1578    BlessLUT(lut);
1579}
1580
1581
1582// Concatenate two LUT into a new single one
1583cmsBool  CMSEXPORT cmsPipelineCat(cmsPipeline* l1, const cmsPipeline* l2)
1584{
1585    cmsStage* mpe;
1586
1587    // If both LUTS does not have elements, we need to inherit
1588    // the number of channels
1589    if (l1 ->Elements == NULL && l2 ->Elements == NULL) {
1590        l1 ->InputChannels  = l2 ->InputChannels;
1591        l1 ->OutputChannels = l2 ->OutputChannels;
1592    }
1593
1594    // Cat second
1595    for (mpe = l2 ->Elements;
1596         mpe != NULL;
1597         mpe = mpe ->Next) {
1598
1599            // We have to dup each element
1600            if (!cmsPipelineInsertStage(l1, cmsAT_END, cmsStageDup(mpe)))
1601                return FALSE;
1602    }
1603
1604    return BlessLUT(l1);
1605}
1606
1607
1608cmsBool CMSEXPORT cmsPipelineSetSaveAs8bitsFlag(cmsPipeline* lut, cmsBool On)
1609{
1610    cmsBool Anterior = lut ->SaveAs8Bits;
1611
1612    lut ->SaveAs8Bits = On;
1613    return Anterior;
1614}
1615
1616
1617cmsStage* CMSEXPORT cmsPipelineGetPtrToFirstStage(const cmsPipeline* lut)
1618{
1619    return lut ->Elements;
1620}
1621
1622cmsStage* CMSEXPORT cmsPipelineGetPtrToLastStage(const cmsPipeline* lut)
1623{
1624    cmsStage *mpe, *Anterior = NULL;
1625
1626    for (mpe = lut ->Elements; mpe != NULL; mpe = mpe ->Next)
1627        Anterior = mpe;
1628
1629    return Anterior;
1630}
1631
1632cmsUInt32Number CMSEXPORT cmsPipelineStageCount(const cmsPipeline* lut)
1633{
1634    cmsStage *mpe;
1635    cmsUInt32Number n;
1636
1637    for (n=0, mpe = lut ->Elements; mpe != NULL; mpe = mpe ->Next)
1638            n++;
1639
1640    return n;
1641}
1642
1643// This function may be used to set the optional evaluator and a block of private data. If private data is being used, an optional
1644// duplicator and free functions should also be specified in order to duplicate the LUT construct. Use NULL to inhibit such functionality.
1645void CMSEXPORT _cmsPipelineSetOptimizationParameters(cmsPipeline* Lut,
1646                                        _cmsOPTeval16Fn Eval16,
1647                                        void* PrivateData,
1648                                        _cmsFreeUserDataFn FreePrivateDataFn,
1649                                        _cmsDupUserDataFn  DupPrivateDataFn)
1650{
1651
1652    Lut ->Eval16Fn = Eval16;
1653    Lut ->DupDataFn = DupPrivateDataFn;
1654    Lut ->FreeDataFn = FreePrivateDataFn;
1655    Lut ->Data = PrivateData;
1656}
1657
1658
1659// ----------------------------------------------------------- Reverse interpolation
1660// Here's how it goes. The derivative Df(x) of the function f is the linear
1661// transformation that best approximates f near the point x. It can be represented
1662// by a matrix A whose entries are the partial derivatives of the components of f
1663// with respect to all the coordinates. This is know as the Jacobian
1664//
1665// The best linear approximation to f is given by the matrix equation:
1666//
1667// y-y0 = A (x-x0)
1668//
1669// So, if x0 is a good "guess" for the zero of f, then solving for the zero of this
1670// linear approximation will give a "better guess" for the zero of f. Thus let y=0,
1671// and since y0=f(x0) one can solve the above equation for x. This leads to the
1672// Newton's method formula:
1673//
1674// xn+1 = xn - A-1 f(xn)
1675//
1676// where xn+1 denotes the (n+1)-st guess, obtained from the n-th guess xn in the
1677// fashion described above. Iterating this will give better and better approximations
1678// if you have a "good enough" initial guess.
1679
1680
1681#define JACOBIAN_EPSILON            0.001f
1682#define INVERSION_MAX_ITERATIONS    30
1683
1684// Increment with reflexion on boundary
1685static
1686void IncDelta(cmsFloat32Number *Val)
1687{
1688    if (*Val < (1.0 - JACOBIAN_EPSILON))
1689
1690        *Val += JACOBIAN_EPSILON;
1691
1692    else
1693        *Val -= JACOBIAN_EPSILON;
1694
1695}
1696
1697
1698
1699// Euclidean distance between two vectors of n elements each one
1700static
1701cmsFloat32Number EuclideanDistance(cmsFloat32Number a[], cmsFloat32Number b[], int n)
1702{
1703    cmsFloat32Number sum = 0;
1704    int i;
1705
1706    for (i=0; i < n; i++) {
1707        cmsFloat32Number dif = b[i] - a[i];
1708        sum +=  dif * dif;
1709    }
1710
1711    return sqrtf(sum);
1712}
1713
1714
1715// Evaluate a LUT in reverse direction. It only searches on 3->3 LUT. Uses Newton method
1716//
1717// x1 <- x - [J(x)]^-1 * f(x)
1718//
1719// lut: The LUT on where to do the search
1720// Target: LabK, 3 values of Lab plus destination K which is fixed
1721// Result: The obtained CMYK
1722// Hint:   Location where begin the search
1723
1724cmsBool CMSEXPORT cmsPipelineEvalReverseFloat(cmsFloat32Number Target[],
1725                                              cmsFloat32Number Result[],
1726                                              cmsFloat32Number Hint[],
1727                                              const cmsPipeline* lut)
1728{
1729    cmsUInt32Number  i, j;
1730    cmsFloat64Number  error, LastError = 1E20;
1731    cmsFloat32Number  fx[4], x[4], xd[4], fxd[4];
1732    cmsVEC3 tmp, tmp2;
1733    cmsMAT3 Jacobian;
1734
1735    // Only 3->3 and 4->3 are supported
1736    if (lut ->InputChannels != 3 && lut ->InputChannels != 4) return FALSE;
1737    if (lut ->OutputChannels != 3) return FALSE;
1738
1739    // Take the hint as starting point if specified
1740    if (Hint == NULL) {
1741
1742        // Begin at any point, we choose 1/3 of CMY axis
1743        x[0] = x[1] = x[2] = 0.3f;
1744    }
1745    else {
1746
1747        // Only copy 3 channels from hint...
1748        for (j=0; j < 3; j++)
1749            x[j] = Hint[j];
1750    }
1751
1752    // If Lut is 4-dimensions, then grab target[3], which is fixed
1753    if (lut ->InputChannels == 4) {
1754        x[3] = Target[3];
1755    }
1756    else x[3] = 0; // To keep lint happy
1757
1758
1759    // Iterate
1760    for (i = 0; i < INVERSION_MAX_ITERATIONS; i++) {
1761
1762        // Get beginning fx
1763        cmsPipelineEvalFloat(x, fx, lut);
1764
1765        // Compute error
1766        error = EuclideanDistance(fx, Target, 3);
1767
1768        // If not convergent, return last safe value
1769        if (error >= LastError)
1770            break;
1771
1772        // Keep latest values
1773        LastError     = error;
1774        for (j=0; j < lut ->InputChannels; j++)
1775                Result[j] = x[j];
1776
1777        // Found an exact match?
1778        if (error <= 0)
1779            break;
1780
1781        // Obtain slope (the Jacobian)
1782        for (j = 0; j < 3; j++) {
1783
1784            xd[0] = x[0];
1785            xd[1] = x[1];
1786            xd[2] = x[2];
1787            xd[3] = x[3];  // Keep fixed channel
1788
1789            IncDelta(&xd[j]);
1790
1791            cmsPipelineEvalFloat(xd, fxd, lut);
1792
1793            Jacobian.v[0].n[j] = ((fxd[0] - fx[0]) / JACOBIAN_EPSILON);
1794            Jacobian.v[1].n[j] = ((fxd[1] - fx[1]) / JACOBIAN_EPSILON);
1795            Jacobian.v[2].n[j] = ((fxd[2] - fx[2]) / JACOBIAN_EPSILON);
1796        }
1797
1798        // Solve system
1799        tmp2.n[0] = fx[0] - Target[0];
1800        tmp2.n[1] = fx[1] - Target[1];
1801        tmp2.n[2] = fx[2] - Target[2];
1802
1803        if (!_cmsMAT3solve(&tmp, &Jacobian, &tmp2))
1804            return FALSE;
1805
1806        // Move our guess
1807        x[0] -= (cmsFloat32Number) tmp.n[0];
1808        x[1] -= (cmsFloat32Number) tmp.n[1];
1809        x[2] -= (cmsFloat32Number) tmp.n[2];
1810
1811        // Some clipping....
1812        for (j=0; j < 3; j++) {
1813            if (x[j] < 0) x[j] = 0;
1814            else
1815                if (x[j] > 1.0) x[j] = 1.0;
1816        }
1817    }
1818
1819    return TRUE;
1820}
1821