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