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