workbench/libs/lcms2/src/cmslut.c

   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
  31 cmsStage* 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
  61 static
  62 void EvaluateIdentity(const cmsFloat32Number In[],
  63                             cmsFloat32Number Out[],
  64                       const cmsStage *mpe)
  65 {
  66     memmove(Out, In, mpe ->InputChannels * sizeof(cmsFloat32Number));
  67 }
  68
  69
  70 cmsStage* 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
  82 static
  83 void 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
  93 static
  94 void 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.
 110 cmsBool  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
 159 cmsToneCurve** _cmsStageGetPtrToCurveSet(const cmsStage* mpe)
 160 {
 161     _cmsStageToneCurvesData* Data = (_cmsStageToneCurvesData*) mpe ->Data;
 162
 163     return Data ->TheCurves;
 164 }
 165
 166 static
 167 void 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
 186 static
 187 void 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
 208 static
 209 void* 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
 233 Error:
 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
 248 cmsStage* 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
 295 cmsStage* _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
 311 static
 312 void 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
 340 static
 341 void* 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
 362 static
 363 void 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
 379 cmsStage*  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, Cols, sizeof(cmsFloat64Number));
 418         if (NewElem->Offset == NULL) {
 419            MatrixElemTypeFree(NewMPE);
 420            return NULL;
 421         }
 422
 423         for (i=0; i < Cols; 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
 440 static
 441 void 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
 450 static
 451 void 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
 466 static
 467 cmsUInt32Number 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
 487 static
 488 void* 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
 530 static
 531 void 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.
 550 cmsStage* 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
 610 cmsStage* 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
 627 cmsStage* 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
 645 cmsStage* 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
 702 static
 703 int 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
 715 cmsStage* _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
 739 cmsUInt16Number _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.
 750 cmsBool 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
 812 cmsBool 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.
 871 cmsBool 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
 901 cmsInt32Number 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
 936 static
 937 void 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.
 965 cmsStage* _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
 978 cmsStage* _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
1016 cmsStage* _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
1032 cmsStage* _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
1052 cmsStage* _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
1074 cmsStage* _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
1091 cmsStage* _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
1111 cmsStage* _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
1134 static
1135 void 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
1159 cmsStage* _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
1169 cmsStage* _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
1183 void CMSEXPORT cmsStageFree(cmsStage* mpe)
1184 {
1185     if (mpe ->FreePtr)
1186         mpe ->FreePtr(mpe);
1187
1188     _cmsFree(mpe ->ContextID, mpe);
1189 }
1190
1191
1192 cmsUInt32Number  CMSEXPORT cmsStageInputChannels(const cmsStage* mpe)
1193 {
1194     return mpe ->InputChannels;
1195 }
1196
1197 cmsUInt32Number  CMSEXPORT cmsStageOutputChannels(const cmsStage* mpe)
1198 {
1199     return mpe ->OutputChannels;
1200 }
1201
1202 cmsStageSignature CMSEXPORT cmsStageType(const cmsStage* mpe)
1203 {
1204     return mpe -> Type;
1205 }
1206
1207 void* CMSEXPORT cmsStageData(const cmsStage* mpe)
1208 {
1209     return mpe -> Data;
1210 }
1211
1212 cmsStage*  CMSEXPORT cmsStageNext(const cmsStage* mpe)
1213 {
1214     return mpe -> Next;
1215 }
1216
1217
1218 // Duplicates an MPE
1219 cmsStage* 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
1258
1259 static
1260 void BlessLUT(cmsPipeline* lut)
1261 {
1262     // We can set the input/ouput channels only if we have elements.
1263     if (lut ->Elements != NULL) {
1264
1265         cmsStage *First, *Last;
1266
1267         First  = cmsPipelineGetPtrToFirstStage(lut);
1268         Last   = cmsPipelineGetPtrToLastStage(lut);
1269
1270         if (First != NULL)lut ->InputChannels = First ->InputChannels;
1271         if (Last != NULL) lut ->OutputChannels = Last ->OutputChannels;
1272     }
1273 }
1274
1275
1276 // Default to evaluate the LUT on 16 bit-basis. Precision is retained.
1277 static
1278 void _LUTeval16(register const cmsUInt16Number In[], register cmsUInt16Number Out[],  register const void* D)
1279 {
1280     cmsPipeline* lut = (cmsPipeline*) D;
1281     cmsStage *mpe;
1282     cmsFloat32Number Storage[2][MAX_STAGE_CHANNELS];
1283     int Phase = 0, NextPhase;
1284
1285     From16ToFloat(In, &Storage[Phase][0], lut ->InputChannels);
1286
1287     for (mpe = lut ->Elements;
1288          mpe != NULL;
1289          mpe = mpe ->Next) {
1290
1291              NextPhase = Phase ^ 1;
1292              mpe ->EvalPtr(&Storage[Phase][0], &Storage[NextPhase][0], mpe);
1293              Phase = NextPhase;
1294     }
1295
1296
1297     FromFloatTo16(&Storage[Phase][0], Out, lut ->OutputChannels);
1298 }
1299
1300
1301
1302 // Does evaluate the LUT on cmsFloat32Number-basis.
1303 static
1304 void _LUTevalFloat(register const cmsFloat32Number In[], register cmsFloat32Number Out[], const void* D)
1305 {
1306     cmsPipeline* lut = (cmsPipeline*) D;
1307     cmsStage *mpe;
1308     cmsFloat32Number Storage[2][MAX_STAGE_CHANNELS];
1309     int Phase = 0, NextPhase;
1310
1311     memmove(&Storage[Phase][0], In, lut ->InputChannels  * sizeof(cmsFloat32Number));
1312
1313     for (mpe = lut ->Elements;
1314          mpe != NULL;
1315          mpe = mpe ->Next) {
1316
1317               NextPhase = Phase ^ 1;
1318               mpe ->EvalPtr(&Storage[Phase][0], &Storage[NextPhase][0], mpe);
1319               Phase = NextPhase;
1320     }
1321
1322     memmove(Out, &Storage[Phase][0], lut ->OutputChannels * sizeof(cmsFloat32Number));
1323 }
1324
1325
1326
1327
1328 // LUT Creation & Destruction
1329
1330 cmsPipeline* CMSEXPORT cmsPipelineAlloc(cmsContext ContextID, cmsUInt32Number InputChannels, cmsUInt32Number OutputChannels)
1331 {
1332        cmsPipeline* NewLUT;
1333
1334        if (InputChannels >= cmsMAXCHANNELS ||
1335            OutputChannels >= cmsMAXCHANNELS) return NULL;
1336
1337        NewLUT = (cmsPipeline*) _cmsMallocZero(ContextID, sizeof(cmsPipeline));
1338        if (NewLUT == NULL) return NULL;
1339
1340
1341        NewLUT -> InputChannels  = InputChannels;
1342        NewLUT -> OutputChannels = OutputChannels;
1343
1344        NewLUT ->Eval16Fn    = _LUTeval16;
1345        NewLUT ->EvalFloatFn = _LUTevalFloat;
1346        NewLUT ->DupDataFn   = NULL;
1347        NewLUT ->FreeDataFn  = NULL;
1348        NewLUT ->Data        = NewLUT;
1349        NewLUT ->ContextID   = ContextID;
1350
1351        BlessLUT(NewLUT);
1352
1353        return NewLUT;
1354 }
1355
1356 cmsContext CMSEXPORT cmsGetPipelineContextID(const cmsPipeline* lut)
1357 {
1358     _cmsAssert(lut != NULL);
1359     return lut ->ContextID;
1360 }
1361
1362 cmsUInt32Number CMSEXPORT cmsPipelineInputChannels(const cmsPipeline* lut)
1363 {
1364     _cmsAssert(lut != NULL);
1365     return lut ->InputChannels;
1366 }
1367
1368 cmsUInt32Number CMSEXPORT cmsPipelineOutputChannels(const cmsPipeline* lut)
1369 {
1370     _cmsAssert(lut != NULL);
1371     return lut ->OutputChannels;
1372 }
1373
1374 // Free a profile elements LUT
1375 void CMSEXPORT cmsPipelineFree(cmsPipeline* lut)
1376 {
1377     cmsStage *mpe, *Next;
1378
1379     if (lut == NULL) return;
1380
1381     for (mpe = lut ->Elements;
1382         mpe != NULL;
1383         mpe = Next) {
1384
1385             Next = mpe ->Next;
1386             cmsStageFree(mpe);
1387     }
1388
1389     if (lut ->FreeDataFn) lut ->FreeDataFn(lut ->ContextID, lut ->Data);
1390
1391     _cmsFree(lut ->ContextID, lut);
1392 }
1393
1394
1395 // Default to evaluate the LUT on 16 bit-basis.
1396 void CMSEXPORT cmsPipelineEval16(const cmsUInt16Number In[], cmsUInt16Number Out[],  const cmsPipeline* lut)
1397 {
1398     _cmsAssert(lut != NULL);
1399     lut ->Eval16Fn(In, Out, lut->Data);
1400 }
1401
1402
1403 // Does evaluate the LUT on cmsFloat32Number-basis.
1404 void CMSEXPORT cmsPipelineEvalFloat(const cmsFloat32Number In[], cmsFloat32Number Out[], const cmsPipeline* lut)
1405 {
1406     _cmsAssert(lut != NULL);
1407     lut ->EvalFloatFn(In, Out, lut);
1408 }
1409
1410
1411
1412 // Duplicates a LUT
1413 cmsPipeline* CMSEXPORT cmsPipelineDup(const cmsPipeline* lut)
1414 {
1415     cmsPipeline* NewLUT;
1416     cmsStage *NewMPE, *Anterior = NULL, *mpe;
1417     cmsBool  First = TRUE;
1418
1419     if (lut == NULL) return NULL;
1420
1421     NewLUT = cmsPipelineAlloc(lut ->ContextID, lut ->InputChannels, lut ->OutputChannels);
1422     if (NewLUT == NULL) return NULL;
1423
1424     for (mpe = lut ->Elements;
1425          mpe != NULL;
1426          mpe = mpe ->Next) {
1427
1428              NewMPE = cmsStageDup(mpe);
1429
1430              if (NewMPE == NULL) {
1431                  cmsPipelineFree(NewLUT);
1432                  return NULL;
1433              }
1434
1435              if (First) {
1436                  NewLUT ->Elements = NewMPE;
1437                  First = FALSE;
1438              }
1439              else {
1440                 Anterior ->Next = NewMPE;
1441              }
1442
1443             Anterior = NewMPE;
1444     }
1445
1446     NewLUT ->Eval16Fn    = lut ->Eval16Fn;
1447     NewLUT ->EvalFloatFn = lut ->EvalFloatFn;
1448     NewLUT ->DupDataFn   = lut ->DupDataFn;
1449     NewLUT ->FreeDataFn  = lut ->FreeDataFn;
1450
1451     if (NewLUT ->DupDataFn != NULL)
1452         NewLUT ->Data = NewLUT ->DupDataFn(lut ->ContextID, lut->Data);
1453
1454
1455     NewLUT ->SaveAs8Bits    = lut ->SaveAs8Bits;
1456
1457     BlessLUT(NewLUT);
1458     return NewLUT;
1459 }
1460
1461
1462 int CMSEXPORT cmsPipelineInsertStage(cmsPipeline* lut, cmsStageLoc loc, cmsStage* mpe)
1463 {
1464     cmsStage* Anterior = NULL, *pt;
1465
1466     if (lut == NULL || mpe == NULL)
1467         return FALSE;
1468
1469     switch (loc) {
1470
1471         case cmsAT_BEGIN:
1472             mpe ->Next = lut ->Elements;
1473             lut ->Elements = mpe;
1474             break;
1475
1476         case cmsAT_END:
1477
1478             if (lut ->Elements == NULL)
1479                 lut ->Elements = mpe;
1480             else {
1481
1482                 for (pt = lut ->Elements;
1483                      pt != NULL;
1484                      pt = pt -> Next) Anterior = pt;
1485
1486                 Anterior ->Next = mpe;
1487                 mpe ->Next = NULL;
1488             }
1489             break;
1490         default:;
1491             return FALSE;
1492     }
1493
1494     BlessLUT(lut);
1495     return TRUE;
1496 }
1497
1498 // Unlink an element and return the pointer to it
1499 void CMSEXPORT cmsPipelineUnlinkStage(cmsPipeline* lut, cmsStageLoc loc, cmsStage** mpe)
1500 {
1501     cmsStage *Anterior, *pt, *Last;
1502     cmsStage *Unlinked = NULL;
1503
1504
1505     // If empty LUT, there is nothing to remove
1506     if (lut ->Elements == NULL) {
1507         if (mpe) *mpe = NULL;
1508         return;
1509     }
1510
1511     // On depending on the strategy...
1512     switch (loc) {
1513
1514         case cmsAT_BEGIN:
1515             {
1516                 cmsStage* elem = lut ->Elements;
1517
1518                 lut ->Elements = elem -> Next;
1519                 elem ->Next = NULL;
1520                 Unlinked = elem;
1521
1522             }
1523             break;
1524
1525         case cmsAT_END:
1526             Anterior = Last = NULL;
1527             for (pt = lut ->Elements;
1528                 pt != NULL;
1529                 pt = pt -> Next) {
1530                     Anterior = Last;
1531                     Last = pt;
1532             }
1533
1534             Unlinked = Last;  // Next already points to NULL
1535
1536             // Truncate the chain
1537             if (Anterior)
1538                 Anterior ->Next = NULL;
1539             else
1540                 lut ->Elements = NULL;
1541             break;
1542         default:;
1543     }
1544
1545     if (mpe)
1546         *mpe = Unlinked;
1547     else
1548         cmsStageFree(Unlinked);
1549
1550     BlessLUT(lut);
1551 }
1552
1553
1554 // Concatenate two LUT into a new single one
1555 cmsBool  CMSEXPORT cmsPipelineCat(cmsPipeline* l1, const cmsPipeline* l2)
1556 {
1557     cmsStage* mpe;
1558
1559     // If both LUTS does not have elements, we need to inherit
1560     // the number of channels
1561     if (l1 ->Elements == NULL && l2 ->Elements == NULL) {
1562         l1 ->InputChannels  = l2 ->InputChannels;
1563         l1 ->OutputChannels = l2 ->OutputChannels;
1564     }
1565
1566     // Cat second
1567     for (mpe = l2 ->Elements;
1568          mpe != NULL;
1569          mpe = mpe ->Next) {
1570
1571             // We have to dup each element
1572             if (!cmsPipelineInsertStage(l1, cmsAT_END, cmsStageDup(mpe)))
1573                 return FALSE;
1574     }
1575
1576     BlessLUT(l1);
1577     return TRUE;
1578 }
1579
1580
1581 cmsBool CMSEXPORT cmsPipelineSetSaveAs8bitsFlag(cmsPipeline* lut, cmsBool On)
1582 {
1583     cmsBool Anterior = lut ->SaveAs8Bits;
1584
1585     lut ->SaveAs8Bits = On;
1586     return Anterior;
1587 }
1588
1589
1590 cmsStage* CMSEXPORT cmsPipelineGetPtrToFirstStage(const cmsPipeline* lut)
1591 {
1592     return lut ->Elements;
1593 }
1594
1595 cmsStage* CMSEXPORT cmsPipelineGetPtrToLastStage(const cmsPipeline* lut)
1596 {
1597     cmsStage *mpe, *Anterior = NULL;
1598
1599     for (mpe = lut ->Elements; mpe != NULL; mpe = mpe ->Next)
1600         Anterior = mpe;
1601
1602     return Anterior;
1603 }
1604
1605 cmsUInt32Number CMSEXPORT cmsPipelineStageCount(const cmsPipeline* lut)
1606 {
1607     cmsStage *mpe;
1608     cmsUInt32Number n;
1609
1610     for (n=0, mpe = lut ->Elements; mpe != NULL; mpe = mpe ->Next)
1611             n++;
1612
1613     return n;
1614 }
1615
1616 // This function may be used to set the optional evaluator and a block of private data. If private data is being used, an optional
1617 // duplicator and free functions should also be specified in order to duplicate the LUT construct. Use NULL to inhibit such functionality.
1618 void CMSEXPORT _cmsPipelineSetOptimizationParameters(cmsPipeline* Lut,
1619                                         _cmsOPTeval16Fn Eval16,
1620                                         void* PrivateData,
1621                                         _cmsFreeUserDataFn FreePrivateDataFn,
1622                                         _cmsDupUserDataFn  DupPrivateDataFn)
1623 {
1624
1625     Lut ->Eval16Fn = Eval16;
1626     Lut ->DupDataFn = DupPrivateDataFn;
1627     Lut ->FreeDataFn = FreePrivateDataFn;
1628     Lut ->Data = PrivateData;
1629 }
1630
1631
1632 // ----------------------------------------------------------- Reverse interpolation
1633 // Here's how it goes. The derivative Df(x) of the function f is the linear
1634 // transformation that best approximates f near the point x. It can be represented
1635 // by a matrix A whose entries are the partial derivatives of the components of f
1636 // with respect to all the coordinates. This is know as the Jacobian
1637 //
1638 // The best linear approximation to f is given by the matrix equation:
1639 //
1640 // y-y0 = A (x-x0)
1641 //
1642 // So, if x0 is a good "guess" for the zero of f, then solving for the zero of this
1643 // linear approximation will give a "better guess" for the zero of f. Thus let y=0,
1644 // and since y0=f(x0) one can solve the above equation for x. This leads to the
1645 // Newton's method formula:
1646 //
1647 // xn+1 = xn - A-1 f(xn)
1648 //
1649 // where xn+1 denotes the (n+1)-st guess, obtained from the n-th guess xn in the
1650 // fashion described above. Iterating this will give better and better approximations
1651 // if you have a "good enough" initial guess.
1652
1653
1654 #define JACOBIAN_EPSILON            0.001f
1655 #define INVERSION_MAX_ITERATIONS    30
1656
1657 // Increment with reflexion on boundary
1658 static
1659 void IncDelta(cmsFloat32Number *Val)
1660 {
1661     if (*Val < (1.0 - JACOBIAN_EPSILON))
1662
1663         *Val += JACOBIAN_EPSILON;
1664
1665     else
1666         *Val -= JACOBIAN_EPSILON;
1667
1668 }
1669
1670
1671
1672 // Euclidean distance between two vectors of n elements each one
1673 static
1674 cmsFloat32Number EuclideanDistance(cmsFloat32Number a[], cmsFloat32Number b[], int n)
1675 {
1676     cmsFloat32Number sum = 0;
1677     int i;
1678
1679     for (i=0; i < n; i++) {
1680         cmsFloat32Number dif = b[i] - a[i];
1681         sum +=  dif * dif;
1682     }
1683
1684     return sqrtf(sum);
1685 }
1686
1687
1688 // Evaluate a LUT in reverse direction. It only searches on 3->3 LUT. Uses Newton method
1689 //
1690 // x1 <- x - [J(x)]^-1 * f(x)
1691 //
1692 // lut: The LUT on where to do the search
1693 // Target: LabK, 3 values of Lab plus destination K which is fixed
1694 // Result: The obtained CMYK
1695 // Hint:   Location where begin the search
1696
1697 cmsBool CMSEXPORT cmsPipelineEvalReverseFloat(cmsFloat32Number Target[],
1698                                               cmsFloat32Number Result[],
1699                                               cmsFloat32Number Hint[],
1700                                               const cmsPipeline* lut)
1701 {
1702     cmsUInt32Number  i, j;
1703     cmsFloat64Number  error, LastError = 1E20;
1704     cmsFloat32Number  fx[4], x[4], xd[4], fxd[4];
1705     cmsVEC3 tmp, tmp2;
1706     cmsMAT3 Jacobian;
1707
1708     // Only 3->3 and 4->3 are supported
1709     if (lut ->InputChannels != 3 && lut ->InputChannels != 4) return FALSE;
1710     if (lut ->OutputChannels != 3) return FALSE;
1711
1712     // Take the hint as starting point if specified
1713     if (Hint == NULL) {
1714
1715         // Begin at any point, we choose 1/3 of CMY axis
1716         x[0] = x[1] = x[2] = 0.3f;
1717     }
1718     else {
1719
1720         // Only copy 3 channels from hint...
1721         for (j=0; j < 3; j++)
1722             x[j] = Hint[j];
1723     }
1724
1725     // If Lut is 4-dimensions, then grab target[3], which is fixed
1726     if (lut ->InputChannels == 4) {
1727         x[3] = Target[3];
1728     }
1729     else x[3] = 0; // To keep lint happy
1730
1731
1732     // Iterate
1733     for (i = 0; i < INVERSION_MAX_ITERATIONS; i++) {
1734
1735         // Get beginning fx
1736         cmsPipelineEvalFloat(x, fx, lut);
1737
1738         // Compute error
1739         error = EuclideanDistance(fx, Target, 3);
1740
1741         // If not convergent, return last safe value
1742         if (error >= LastError)
1743             break;
1744
1745         // Keep latest values
1746         LastError     = error;
1747         for (j=0; j < lut ->InputChannels; j++)
1748                 Result[j] = x[j];
1749
1750         // Found an exact match?
1751         if (error <= 0)
1752             break;
1753
1754         // Obtain slope (the Jacobian)
1755         for (j = 0; j < 3; j++) {
1756
1757             xd[0] = x[0];
1758             xd[1] = x[1];
1759             xd[2] = x[2];
1760             xd[3] = x[3];  // Keep fixed channel
1761
1762             IncDelta(&xd[j]);
1763
1764             cmsPipelineEvalFloat(xd, fxd, lut);
1765
1766             Jacobian.v[0].n[j] = ((fxd[0] - fx[0]) / JACOBIAN_EPSILON);
1767             Jacobian.v[1].n[j] = ((fxd[1] - fx[1]) / JACOBIAN_EPSILON);
1768             Jacobian.v[2].n[j] = ((fxd[2] - fx[2]) / JACOBIAN_EPSILON);
1769         }
1770
1771         // Solve system
1772         tmp2.n[0] = fx[0] - Target[0];
1773         tmp2.n[1] = fx[1] - Target[1];
1774         tmp2.n[2] = fx[2] - Target[2];
1775
1776         if (!_cmsMAT3solve(&tmp, &Jacobian, &tmp2))
1777             return FALSE;
1778
1779         // Move our guess
1780         x[0] -= (cmsFloat32Number) tmp.n[0];
1781         x[1] -= (cmsFloat32Number) tmp.n[1];
1782         x[2] -= (cmsFloat32Number) tmp.n[2];
1783
1784         // Some clipping....
1785         for (j=0; j < 3; j++) {
1786             if (x[j] < 0) x[j] = 0;
1787             else
1788                 if (x[j] > 1.0) x[j] = 1.0;
1789         }
1790     }
1791
1792     return TRUE;
1793 }
1794
1795