libs/lcms2/src/cmslut.c

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