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 #include "lcms2_internal.h"
  56 
  57 // Tone curves are powerful constructs that can contain curves specified in diverse ways.
  58 // The curve is stored in segments, where each segment can be sampled or specified by parameters.
  59 // a 16.bit simplification of the *whole* curve is kept for optimization purposes. For float operation,
  60 // each segment is evaluated separately. Plug-ins may be used to define new parametric schemes,
  61 // each plug-in may define up to MAX_TYPES_IN_LCMS_PLUGIN functions types. For defining a function,
  62 // the plug-in should provide the type id, how many parameters each type has, and a pointer to
  63 // a procedure that evaluates the function. In the case of reverse evaluation, the evaluator will
  64 // be called with the type id as a negative value, and a sampled version of the reversed curve
  65 // will be built.
  66 
  67 // ----------------------------------------------------------------- Implementation
  68 // Maxim number of nodes
  69 #define MAX_NODES_IN_CURVE   4097
  70 #define MINUS_INF            (-1E22F)
  71 #define PLUS_INF             (+1E22F)
  72 
  73 // The list of supported parametric curves
  74 typedef struct _cmsParametricCurvesCollection_st {
  75 
  76     cmsUInt32Number nFunctions;                                     // Number of supported functions in this chunk
  77     cmsInt32Number  FunctionTypes[MAX_TYPES_IN_LCMS_PLUGIN];        // The identification types
  78     cmsUInt32Number ParameterCount[MAX_TYPES_IN_LCMS_PLUGIN];       // Number of parameters for each function
  79 
  80     cmsParametricCurveEvaluator Evaluator;                          // The evaluator
  81 
  82     struct _cmsParametricCurvesCollection_st* Next; // Next in list
  83 
  84 } _cmsParametricCurvesCollection;
  85 
  86 // This is the default (built-in) evaluator
  87 static cmsFloat64Number DefaultEvalParametricFn(cmsInt32Number Type, const cmsFloat64Number Params[], cmsFloat64Number R);
  88 
  89 // The built-in list
  90 static _cmsParametricCurvesCollection DefaultCurves = {
  91     10,                                      // # of curve types
  92     { 1, 2, 3, 4, 5, 6, 7, 8, 108, 109 },    // Parametric curve ID
  93     { 1, 3, 4, 5, 7, 4, 5, 5,   1,   1 },    // Parameters by type
  94     DefaultEvalParametricFn,                 // Evaluator
  95     NULL                                     // Next in chain
  96 };
  97 
  98 // Duplicates the zone of memory used by the plug-in in the new context
  99 static
 100 void DupPluginCurvesList(struct _cmsContext_struct* ctx,
 101                                                const struct _cmsContext_struct* src)
 102 {
 103    _cmsCurvesPluginChunkType newHead = { NULL };
 104    _cmsParametricCurvesCollection*  entry;
 105    _cmsParametricCurvesCollection*  Anterior = NULL;
 106    _cmsCurvesPluginChunkType* head = (_cmsCurvesPluginChunkType*) src->chunks[CurvesPlugin];
 107 
 108     _cmsAssert(head != NULL);
 109 
 110     // Walk the list copying all nodes
 111    for (entry = head->ParametricCurves;
 112         entry != NULL;
 113         entry = entry ->Next) {
 114 
 115             _cmsParametricCurvesCollection *newEntry = ( _cmsParametricCurvesCollection *) _cmsSubAllocDup(ctx ->MemPool, entry, sizeof(_cmsParametricCurvesCollection));
 116 
 117             if (newEntry == NULL)
 118                 return;
 119 
 120             // We want to keep the linked list order, so this is a little bit tricky
 121             newEntry -> Next = NULL;
 122             if (Anterior)
 123                 Anterior -> Next = newEntry;
 124 
 125             Anterior = newEntry;
 126 
 127             if (newHead.ParametricCurves == NULL)
 128                 newHead.ParametricCurves = newEntry;
 129     }
 130 
 131   ctx ->chunks[CurvesPlugin] = _cmsSubAllocDup(ctx->MemPool, &newHead, sizeof(_cmsCurvesPluginChunkType));
 132 }
 133 
 134 // The allocator have to follow the chain
 135 void _cmsAllocCurvesPluginChunk(struct _cmsContext_struct* ctx,
 136                                 const struct _cmsContext_struct* src)
 137 {
 138     _cmsAssert(ctx != NULL);
 139 
 140     if (src != NULL) {
 141 
 142         // Copy all linked list
 143        DupPluginCurvesList(ctx, src);
 144     }
 145     else {
 146         static _cmsCurvesPluginChunkType CurvesPluginChunk = { NULL };
 147         ctx ->chunks[CurvesPlugin] = _cmsSubAllocDup(ctx ->MemPool, &CurvesPluginChunk, sizeof(_cmsCurvesPluginChunkType));
 148     }
 149 }
 150 
 151 
 152 // The linked list head
 153 _cmsCurvesPluginChunkType _cmsCurvesPluginChunk = { NULL };
 154 
 155 // As a way to install new parametric curves
 156 cmsBool _cmsRegisterParametricCurvesPlugin(cmsContext ContextID, cmsPluginBase* Data)
 157 {
 158     _cmsCurvesPluginChunkType* ctx = ( _cmsCurvesPluginChunkType*) _cmsContextGetClientChunk(ContextID, CurvesPlugin);
 159     cmsPluginParametricCurves* Plugin = (cmsPluginParametricCurves*) Data;
 160     _cmsParametricCurvesCollection* fl;
 161 
 162     if (Data == NULL) {
 163 
 164           ctx -> ParametricCurves =  NULL;
 165           return TRUE;
 166     }
 167 
 168     fl = (_cmsParametricCurvesCollection*) _cmsPluginMalloc(ContextID, sizeof(_cmsParametricCurvesCollection));
 169     if (fl == NULL) return FALSE;
 170 
 171     // Copy the parameters
 172     fl ->Evaluator  = Plugin ->Evaluator;
 173     fl ->nFunctions = Plugin ->nFunctions;
 174 
 175     // Make sure no mem overwrites
 176     if (fl ->nFunctions > MAX_TYPES_IN_LCMS_PLUGIN)
 177         fl ->nFunctions = MAX_TYPES_IN_LCMS_PLUGIN;
 178 
 179     // Copy the data
 180     memmove(fl->FunctionTypes,  Plugin ->FunctionTypes,   fl->nFunctions * sizeof(cmsUInt32Number));
 181     memmove(fl->ParameterCount, Plugin ->ParameterCount,  fl->nFunctions * sizeof(cmsUInt32Number));
 182 
 183     // Keep linked list
 184     fl ->Next = ctx->ParametricCurves;
 185     ctx->ParametricCurves = fl;
 186 
 187     // All is ok
 188     return TRUE;
 189 }
 190 
 191 
 192 // Search in type list, return position or -1 if not found
 193 static
 194 int IsInSet(int Type, _cmsParametricCurvesCollection* c)
 195 {
 196     int i;
 197 
 198     for (i=0; i < (int) c ->nFunctions; i++)
 199         if (abs(Type) == c ->FunctionTypes[i]) return i;
 200 
 201     return -1;
 202 }
 203 
 204 
 205 // Search for the collection which contains a specific type
 206 static
 207 _cmsParametricCurvesCollection *GetParametricCurveByType(cmsContext ContextID, int Type, int* index)
 208 {
 209     _cmsParametricCurvesCollection* c;
 210     int Position;
 211     _cmsCurvesPluginChunkType* ctx = ( _cmsCurvesPluginChunkType*) _cmsContextGetClientChunk(ContextID, CurvesPlugin);
 212 
 213     for (c = ctx->ParametricCurves; c != NULL; c = c ->Next) {
 214 
 215         Position = IsInSet(Type, c);
 216 
 217         if (Position != -1) {
 218             if (index != NULL)
 219                 *index = Position;
 220             return c;
 221         }
 222     }
 223     // If none found, revert for defaults
 224     for (c = &DefaultCurves; c != NULL; c = c ->Next) {
 225 
 226         Position = IsInSet(Type, c);
 227 
 228         if (Position != -1) {
 229             if (index != NULL)
 230                 *index = Position;
 231             return c;
 232         }
 233     }
 234 
 235     return NULL;
 236 }
 237 
 238 // Low level allocate, which takes care of memory details. nEntries may be zero, and in this case
 239 // no optimization curve is computed. nSegments may also be zero in the inverse case, where only the
 240 // optimization curve is given. Both features simultaneously is an error
 241 static
 242 cmsToneCurve* AllocateToneCurveStruct(cmsContext ContextID, cmsUInt32Number nEntries,
 243                                       cmsUInt32Number nSegments, const cmsCurveSegment* Segments,
 244                                       const cmsUInt16Number* Values)
 245 {
 246     cmsToneCurve* p;
 247     cmsUInt32Number i;
 248 
 249     // We allow huge tables, which are then restricted for smoothing operations
 250     if (nEntries > 65530) {
 251         cmsSignalError(ContextID, cmsERROR_RANGE, "Couldn't create tone curve of more than 65530 entries");
 252         return NULL;
 253     }
 254 
 255     if (nEntries == 0 && nSegments == 0) {
 256         cmsSignalError(ContextID, cmsERROR_RANGE, "Couldn't create tone curve with zero segments and no table");
 257         return NULL;
 258     }
 259 
 260     // Allocate all required pointers, etc.
 261     p = (cmsToneCurve*) _cmsMallocZero(ContextID, sizeof(cmsToneCurve));
 262     if (!p) return NULL;
 263 
 264     // In this case, there are no segments
 265     if (nSegments == 0) {
 266         p ->Segments = NULL;
 267         p ->Evals = NULL;
 268     }
 269     else {
 270         p ->Segments = (cmsCurveSegment*) _cmsCalloc(ContextID, nSegments, sizeof(cmsCurveSegment));
 271         if (p ->Segments == NULL) goto Error;
 272 
 273         p ->Evals    = (cmsParametricCurveEvaluator*) _cmsCalloc(ContextID, nSegments, sizeof(cmsParametricCurveEvaluator));
 274         if (p ->Evals == NULL) goto Error;
 275     }
 276 
 277     p -> nSegments = nSegments;
 278 
 279     // This 16-bit table contains a limited precision representation of the whole curve and is kept for
 280     // increasing xput on certain operations.
 281     if (nEntries == 0) {
 282         p ->Table16 = NULL;
 283     }
 284     else {
 285        p ->Table16 = (cmsUInt16Number*)  _cmsCalloc(ContextID, nEntries, sizeof(cmsUInt16Number));
 286        if (p ->Table16 == NULL) goto Error;
 287     }
 288 
 289     p -> nEntries  = nEntries;
 290 
 291     // Initialize members if requested
 292     if (Values != NULL && (nEntries > 0)) {
 293 
 294         for (i=0; i < nEntries; i++)
 295             p ->Table16[i] = Values[i];
 296     }
 297 
 298     // Initialize the segments stuff. The evaluator for each segment is located and a pointer to it
 299     // is placed in advance to maximize performance.
 300     if (Segments != NULL && (nSegments > 0)) {
 301 
 302         _cmsParametricCurvesCollection *c;
 303 
 304         p ->SegInterp = (cmsInterpParams**) _cmsCalloc(ContextID, nSegments, sizeof(cmsInterpParams*));
 305         if (p ->SegInterp == NULL) goto Error;
 306 
 307         for (i=0; i < nSegments; i++) {
 308 
 309             // Type 0 is a special marker for table-based curves
 310             if (Segments[i].Type == 0)
 311                 p ->SegInterp[i] = _cmsComputeInterpParams(ContextID, Segments[i].nGridPoints, 1, 1, NULL, CMS_LERP_FLAGS_FLOAT);
 312 
 313             memmove(&p ->Segments[i], &Segments[i], sizeof(cmsCurveSegment));
 314 
 315             if (Segments[i].Type == 0 && Segments[i].SampledPoints != NULL)
 316                 p ->Segments[i].SampledPoints = (cmsFloat32Number*) _cmsDupMem(ContextID, Segments[i].SampledPoints, sizeof(cmsFloat32Number) * Segments[i].nGridPoints);
 317             else
 318                 p ->Segments[i].SampledPoints = NULL;
 319 
 320 
 321             c = GetParametricCurveByType(ContextID, Segments[i].Type, NULL);
 322             if (c != NULL)
 323                     p ->Evals[i] = c ->Evaluator;
 324         }
 325     }
 326 
 327     p ->InterpParams = _cmsComputeInterpParams(ContextID, p ->nEntries, 1, 1, p->Table16, CMS_LERP_FLAGS_16BITS);
 328     if (p->InterpParams != NULL)
 329         return p;
 330 
 331 Error:
 332     for (i=0; i < nSegments; i++) {
 333         if (p ->Segments && p ->Segments[i].SampledPoints) _cmsFree(ContextID, p ->Segments[i].SampledPoints);
 334         if (p ->SegInterp && p ->SegInterp[i]) _cmsFree(ContextID, p ->SegInterp[i]);
 335     }
 336     if (p -> SegInterp) _cmsFree(ContextID, p -> SegInterp);
 337     if (p -> Segments) _cmsFree(ContextID, p -> Segments);
 338     if (p -> Evals) _cmsFree(ContextID, p -> Evals);
 339     if (p ->Table16) _cmsFree(ContextID, p ->Table16);
 340     _cmsFree(ContextID, p);
 341     return NULL;
 342 }
 343 
 344 
 345 // Generates a sigmoidal function with desired steepness.
 346 cmsINLINE double sigmoid_base(double k, double t)
 347 {
 348     return (1.0 / (1.0 + exp(-k * t))) - 0.5;
 349 }
 350 
 351 cmsINLINE double inverted_sigmoid_base(double k, double t)
 352 {
 353     return -log((1.0 / (t + 0.5)) - 1.0) / k;
 354 }
 355 
 356 cmsINLINE double sigmoid_factory(double k, double t)
 357 {
 358     double correction = 0.5 / sigmoid_base(k, 1);
 359 
 360     return correction * sigmoid_base(k, 2.0 * t - 1.0) + 0.5;
 361 }
 362 
 363 cmsINLINE double inverse_sigmoid_factory(double k, double t)
 364 {
 365     double correction = 0.5 / sigmoid_base(k, 1);
 366 
 367     return (inverted_sigmoid_base(k, (t - 0.5) / correction) + 1.0) / 2.0;
 368 }
 369 
 370 
 371 // Parametric Fn using floating point
 372 static
 373 cmsFloat64Number DefaultEvalParametricFn(cmsInt32Number Type, const cmsFloat64Number Params[], cmsFloat64Number R)
 374 {
 375     cmsFloat64Number e, Val, disc;
 376 
 377     switch (Type) {
 378 
 379    // X = Y ^ Gamma
 380     case 1:
 381         if (R < 0) {
 382 
 383             if (fabs(Params[0] - 1.0) < MATRIX_DET_TOLERANCE)
 384                 Val = R;
 385             else
 386                 Val = 0;
 387         }
 388         else
 389             Val = pow(R, Params[0]);
 390         break;
 391 
 392     // Type 1 Reversed: X = Y ^1/gamma
 393     case -1:
 394         if (R < 0) {
 395 
 396             if (fabs(Params[0] - 1.0) < MATRIX_DET_TOLERANCE)
 397                 Val = R;
 398             else
 399                 Val = 0;
 400         }
 401         else
 402         {
 403             if (fabs(Params[0]) < MATRIX_DET_TOLERANCE)
 404                 Val = PLUS_INF;
 405             else
 406                 Val = pow(R, 1 / Params[0]);
 407         }
 408         break;
 409 
 410     // CIE 122-1966
 411     // Y = (aX + b)^Gamma  | X >= -b/a
 412     // Y = 0               | else
 413     case 2:
 414     {
 415 
 416         if (fabs(Params[1]) < MATRIX_DET_TOLERANCE)
 417         {
 418             Val = 0;
 419         }
 420         else
 421         {
 422             disc = -Params[2] / Params[1];
 423 
 424             if (R >= disc) {
 425 
 426                 e = Params[1] * R + Params[2];
 427 
 428                 if (e > 0)
 429                     Val = pow(e, Params[0]);
 430                 else
 431                     Val = 0;
 432             }
 433             else
 434                 Val = 0;
 435         }
 436     }
 437     break;
 438 
 439      // Type 2 Reversed
 440      // X = (Y ^1/g  - b) / a
 441      case -2:
 442      {
 443          if (fabs(Params[0]) < MATRIX_DET_TOLERANCE ||
 444              fabs(Params[1]) < MATRIX_DET_TOLERANCE)
 445          {
 446              Val = 0;
 447          }
 448          else
 449          {
 450              if (R < 0)
 451                  Val = 0;
 452              else
 453                  Val = (pow(R, 1.0 / Params[0]) - Params[2]) / Params[1];
 454 
 455              if (Val < 0)
 456                  Val = 0;
 457          }
 458      }
 459      break;
 460 
 461 
 462     // IEC 61966-3
 463     // Y = (aX + b)^Gamma + c | X <= -b/a
 464     // Y = c                  | else
 465     case 3:
 466     {
 467         if (fabs(Params[1]) < MATRIX_DET_TOLERANCE)
 468         {
 469             Val = 0;
 470         }
 471         else
 472         {
 473             disc = -Params[2] / Params[1];
 474             if (disc < 0)
 475                 disc = 0;
 476 
 477             if (R >= disc) {
 478 
 479                 e = Params[1] * R + Params[2];
 480 
 481                 if (e > 0)
 482                     Val = pow(e, Params[0]) + Params[3];
 483                 else
 484                     Val = 0;
 485             }
 486             else
 487                 Val = Params[3];
 488         }
 489     }
 490     break;
 491 
 492 
 493     // Type 3 reversed
 494     // X=((Y-c)^1/g - b)/a      | (Y>=c)
 495     // X=-b/a                   | (Y<c)
 496     case -3:
 497     {
 498         if (fabs(Params[0]) < MATRIX_DET_TOLERANCE ||
 499             fabs(Params[1]) < MATRIX_DET_TOLERANCE)
 500         {
 501             Val = 0;
 502         }
 503         else
 504         {
 505             if (R >= Params[3]) {
 506 
 507                 e = R - Params[3];
 508 
 509                 if (e > 0)
 510                     Val = (pow(e, 1 / Params[0]) - Params[2]) / Params[1];
 511                 else
 512                     Val = 0;
 513             }
 514             else {
 515                 Val = -Params[2] / Params[1];
 516             }
 517         }
 518     }
 519     break;
 520 
 521 
 522     // IEC 61966-2.1 (sRGB)
 523     // Y = (aX + b)^Gamma | X >= d
 524     // Y = cX             | X < d
 525     case 4:
 526         if (R >= Params[4]) {
 527 
 528             e = Params[1]*R + Params[2];
 529 
 530             if (e > 0)
 531                 Val = pow(e, Params[0]);
 532             else
 533                 Val = 0;
 534         }
 535         else
 536             Val = R * Params[3];
 537         break;
 538 
 539     // Type 4 reversed
 540     // X=((Y^1/g-b)/a)    | Y >= (ad+b)^g
 541     // X=Y/c              | Y< (ad+b)^g
 542     case -4:
 543     {
 544 
 545         e = Params[1] * Params[4] + Params[2];
 546         if (e < 0)
 547             disc = 0;
 548         else
 549             disc = pow(e, Params[0]);
 550 
 551         if (R >= disc) {
 552 
 553             if (fabs(Params[0]) < MATRIX_DET_TOLERANCE ||
 554                 fabs(Params[1]) < MATRIX_DET_TOLERANCE)
 555 
 556                 Val = 0;
 557 
 558             else
 559                 Val = (pow(R, 1.0 / Params[0]) - Params[2]) / Params[1];
 560         }
 561         else {
 562 
 563             if (fabs(Params[3]) < MATRIX_DET_TOLERANCE)
 564                 Val = 0;
 565             else
 566                 Val = R / Params[3];
 567         }
 568 
 569     }
 570     break;
 571 
 572 
 573     // Y = (aX + b)^Gamma + e | X >= d
 574     // Y = cX + f             | X < d
 575     case 5:
 576         if (R >= Params[4]) {
 577 
 578             e = Params[1]*R + Params[2];
 579 
 580             if (e > 0)
 581                 Val = pow(e, Params[0]) + Params[5];
 582             else
 583                 Val = Params[5];
 584         }
 585         else
 586             Val = R*Params[3] + Params[6];
 587         break;
 588 
 589 
 590     // Reversed type 5
 591     // X=((Y-e)1/g-b)/a   | Y >=(ad+b)^g+e), cd+f
 592     // X=(Y-f)/c          | else
 593     case -5:
 594     {
 595         disc = Params[3] * Params[4] + Params[6];
 596         if (R >= disc) {
 597 
 598             e = R - Params[5];
 599             if (e < 0)
 600                 Val = 0;
 601             else
 602             {
 603                 if (fabs(Params[0]) < MATRIX_DET_TOLERANCE ||
 604                     fabs(Params[1]) < MATRIX_DET_TOLERANCE)
 605 
 606                     Val = 0;
 607                 else
 608                     Val = (pow(e, 1.0 / Params[0]) - Params[2]) / Params[1];
 609             }
 610         }
 611         else {
 612             if (fabs(Params[3]) < MATRIX_DET_TOLERANCE)
 613                 Val = 0;
 614             else
 615                 Val = (R - Params[6]) / Params[3];
 616         }
 617 
 618     }
 619     break;
 620 
 621 
 622     // Types 6,7,8 comes from segmented curves as described in ICCSpecRevision_02_11_06_Float.pdf
 623     // Type 6 is basically identical to type 5 without d
 624 
 625     // Y = (a * X + b) ^ Gamma + c
 626     case 6:
 627         e = Params[1]*R + Params[2];
 628 
 629         // On gamma 1.0, don't clamp
 630         if (Params[0] == 1.0) {
 631             Val = e + Params[3];
 632         }
 633         else {
 634             if (e < 0)
 635                 Val = Params[3];
 636             else
 637                 Val = pow(e, Params[0]) + Params[3];
 638         }
 639         break;
 640 
 641     // ((Y - c) ^1/Gamma - b) / a
 642     case -6:
 643     {
 644         if (fabs(Params[0]) < MATRIX_DET_TOLERANCE ||
 645             fabs(Params[1]) < MATRIX_DET_TOLERANCE)
 646         {
 647             Val = 0;
 648         }
 649         else
 650         {
 651             e = R - Params[3];
 652             if (e < 0)
 653                 Val = 0;
 654             else
 655                 Val = (pow(e, 1.0 / Params[0]) - Params[2]) / Params[1];
 656         }
 657     }
 658     break;
 659 
 660 
 661     // Y = a * log (b * X^Gamma + c) + d
 662     case 7:
 663 
 664        e = Params[2] * pow(R, Params[0]) + Params[3];
 665        if (e <= 0)
 666            Val = Params[4];
 667        else
 668            Val = Params[1]*log10(e) + Params[4];
 669        break;
 670 
 671     // (Y - d) / a = log(b * X ^Gamma + c)
 672     // pow(10, (Y-d) / a) = b * X ^Gamma + c
 673     // pow((pow(10, (Y-d) / a) - c) / b, 1/g) = X
 674     case -7:
 675     {
 676         if (fabs(Params[0]) < MATRIX_DET_TOLERANCE ||
 677             fabs(Params[1]) < MATRIX_DET_TOLERANCE ||
 678             fabs(Params[2]) < MATRIX_DET_TOLERANCE)
 679         {
 680             Val = 0;
 681         }
 682         else
 683         {
 684             Val = pow((pow(10.0, (R - Params[4]) / Params[1]) - Params[3]) / Params[2], 1.0 / Params[0]);
 685         }
 686     }
 687     break;
 688 
 689 
 690    //Y = a * b^(c*X+d) + e
 691    case 8:
 692        Val = (Params[0] * pow(Params[1], Params[2] * R + Params[3]) + Params[4]);
 693        break;
 694 
 695 
 696    // Y = (log((y-e) / a) / log(b) - d ) / c
 697    // a=0, b=1, c=2, d=3, e=4,
 698    case -8:
 699 
 700        disc = R - Params[4];
 701        if (disc < 0) Val = 0;
 702        else
 703        {
 704            if (fabs(Params[0]) < MATRIX_DET_TOLERANCE ||
 705                fabs(Params[2]) < MATRIX_DET_TOLERANCE)
 706            {
 707                Val = 0;
 708            }
 709            else
 710            {
 711                Val = (log(disc / Params[0]) / log(Params[1]) - Params[3]) / Params[2];
 712            }
 713        }
 714        break;
 715 
 716 
 717    // S-Shaped: (1 - (1-x)^1/g)^1/g
 718    case 108:
 719        if (fabs(Params[0]) < MATRIX_DET_TOLERANCE)
 720            Val = 0;
 721        else
 722            Val = pow(1.0 - pow(1 - R, 1/Params[0]), 1/Params[0]);
 723       break;
 724 
 725     // y = (1 - (1-x)^1/g)^1/g
 726     // y^g = (1 - (1-x)^1/g)
 727     // 1 - y^g = (1-x)^1/g
 728     // (1 - y^g)^g = 1 - x
 729     // 1 - (1 - y^g)^g
 730     case -108:
 731         Val = 1 - pow(1 - pow(R, Params[0]), Params[0]);
 732         break;
 733 
 734     // Sigmoidals
 735     case 109:
 736         Val = sigmoid_factory(Params[0], R);
 737         break;
 738 
 739     case -109:
 740         Val = inverse_sigmoid_factory(Params[0], R);
 741         break;
 742 
 743     default:
 744         // Unsupported parametric curve. Should never reach here
 745         return 0;
 746     }
 747 
 748     return Val;
 749 }
 750 
 751 // Evaluate a segmented function for a single value. Return -Inf if no valid segment found .
 752 // If fn type is 0, perform an interpolation on the table
 753 static
 754 cmsFloat64Number EvalSegmentedFn(const cmsToneCurve *g, cmsFloat64Number R)
 755 {
 756     int i;
 757     cmsFloat32Number Out32;
 758     cmsFloat64Number Out;
 759 
 760     for (i = (int) g->nSegments - 1; i >= 0; --i) {
 761 
 762         // Check for domain
 763         if ((R > g->Segments[i].x0) && (R <= g->Segments[i].x1)) {
 764 
 765             // Type == 0 means segment is sampled
 766             if (g->Segments[i].Type == 0) {
 767 
 768                 cmsFloat32Number R1 = (cmsFloat32Number)(R - g->Segments[i].x0) / (g->Segments[i].x1 - g->Segments[i].x0);
 769 
 770                 // Setup the table (TODO: clean that)
 771                 g->SegInterp[i]->Table = g->Segments[i].SampledPoints;
 772 
 773                 g->SegInterp[i]->Interpolation.LerpFloat(&R1, &Out32, g->SegInterp[i]);
 774                 Out = (cmsFloat64Number) Out32;
 775 
 776             }
 777             else {
 778                 Out = g->Evals[i](g->Segments[i].Type, g->Segments[i].Params, R);
 779             }
 780 
 781             if (isinf(Out))
 782                 return PLUS_INF;
 783             else
 784             {
 785                 if (isinf(-Out))
 786                     return MINUS_INF;
 787             }
 788 
 789             return Out;
 790         }
 791     }
 792 
 793     return MINUS_INF;
 794 }
 795 
 796 // Access to estimated low-res table
 797 cmsUInt32Number CMSEXPORT cmsGetToneCurveEstimatedTableEntries(const cmsToneCurve* t)
 798 {
 799     _cmsAssert(t != NULL);
 800     return t ->nEntries;
 801 }
 802 
 803 const cmsUInt16Number* CMSEXPORT cmsGetToneCurveEstimatedTable(const cmsToneCurve* t)
 804 {
 805     _cmsAssert(t != NULL);
 806     return t ->Table16;
 807 }
 808 
 809 
 810 // Create an empty gamma curve, by using tables. This specifies only the limited-precision part, and leaves the
 811 // floating point description empty.
 812 cmsToneCurve* CMSEXPORT cmsBuildTabulatedToneCurve16(cmsContext ContextID, cmsUInt32Number nEntries, const cmsUInt16Number Values[])
 813 {
 814     return AllocateToneCurveStruct(ContextID, nEntries, 0, NULL, Values);
 815 }
 816 
 817 static
 818 cmsUInt32Number EntriesByGamma(cmsFloat64Number Gamma)
 819 {
 820     if (fabs(Gamma - 1.0) < 0.001) return 2;
 821     return 4096;
 822 }
 823 
 824 
 825 // Create a segmented gamma, fill the table
 826 cmsToneCurve* CMSEXPORT cmsBuildSegmentedToneCurve(cmsContext ContextID,
 827                                                    cmsUInt32Number nSegments, const cmsCurveSegment Segments[])
 828 {
 829     cmsUInt32Number i;
 830     cmsFloat64Number R, Val;
 831     cmsToneCurve* g;
 832     cmsUInt32Number nGridPoints = 4096;
 833 
 834     _cmsAssert(Segments != NULL);
 835 
 836     // Optimizatin for identity curves.
 837     if (nSegments == 1 && Segments[0].Type == 1) {
 838 
 839         nGridPoints = EntriesByGamma(Segments[0].Params[0]);
 840     }
 841 
 842     g = AllocateToneCurveStruct(ContextID, nGridPoints, nSegments, Segments, NULL);
 843     if (g == NULL) return NULL;
 844 
 845     // Once we have the floating point version, we can approximate a 16 bit table of 4096 entries
 846     // for performance reasons. This table would normally not be used except on 8/16 bits transforms.
 847     for (i = 0; i < nGridPoints; i++) {
 848 
 849         R   = (cmsFloat64Number) i / (nGridPoints-1);
 850 
 851         Val = EvalSegmentedFn(g, R);
 852 
 853         // Round and saturate
 854         g ->Table16[i] = _cmsQuickSaturateWord(Val * 65535.0);
 855     }
 856 
 857     return g;
 858 }
 859 
 860 // Use a segmented curve to store the floating point table
 861 cmsToneCurve* CMSEXPORT cmsBuildTabulatedToneCurveFloat(cmsContext ContextID, cmsUInt32Number nEntries, const cmsFloat32Number values[])
 862 {
 863     cmsCurveSegment Seg[3];
 864 
 865     // Do some housekeeping
 866     if (nEntries == 0 || values == NULL)
 867         return NULL;
 868 
 869     // A segmented tone curve should have function segments in the first and last positions
 870     // Initialize segmented curve part up to 0 to constant value = samples[0]
 871     Seg[0].x0 = MINUS_INF;
 872     Seg[0].x1 = 0;
 873     Seg[0].Type = 6;
 874 
 875     Seg[0].Params[0] = 1;
 876     Seg[0].Params[1] = 0;
 877     Seg[0].Params[2] = 0;
 878     Seg[0].Params[3] = values[0];
 879     Seg[0].Params[4] = 0;
 880 
 881     // From zero to 1
 882     Seg[1].x0 = 0;
 883     Seg[1].x1 = 1.0;
 884     Seg[1].Type = 0;
 885 
 886     Seg[1].nGridPoints = nEntries;
 887     Seg[1].SampledPoints = (cmsFloat32Number*) values;
 888 
 889     // Final segment is constant = lastsample
 890     Seg[2].x0 = 1.0;
 891     Seg[2].x1 = PLUS_INF;
 892     Seg[2].Type = 6;
 893 
 894     Seg[2].Params[0] = 1;
 895     Seg[2].Params[1] = 0;
 896     Seg[2].Params[2] = 0;
 897     Seg[2].Params[3] = values[nEntries-1];
 898     Seg[2].Params[4] = 0;
 899 
 900 
 901     return cmsBuildSegmentedToneCurve(ContextID, 3, Seg);
 902 }
 903 
 904 // Parametric curves
 905 //
 906 // Parameters goes as: Curve, a, b, c, d, e, f
 907 // Type is the ICC type +1
 908 // if type is negative, then the curve is analytically inverted
 909 cmsToneCurve* CMSEXPORT cmsBuildParametricToneCurve(cmsContext ContextID, cmsInt32Number Type, const cmsFloat64Number Params[])
 910 {
 911     cmsCurveSegment Seg0;
 912     int Pos = 0;
 913     cmsUInt32Number size;
 914     _cmsParametricCurvesCollection* c = GetParametricCurveByType(ContextID, Type, &Pos);
 915 
 916     _cmsAssert(Params != NULL);
 917 
 918     if (c == NULL) {
 919         cmsSignalError(ContextID, cmsERROR_UNKNOWN_EXTENSION, "Invalid parametric curve type %d", Type);
 920         return NULL;
 921     }
 922 
 923     memset(&Seg0, 0, sizeof(Seg0));
 924 
 925     Seg0.x0   = MINUS_INF;
 926     Seg0.x1   = PLUS_INF;
 927     Seg0.Type = Type;
 928 
 929     size = c->ParameterCount[Pos] * sizeof(cmsFloat64Number);
 930     memmove(Seg0.Params, Params, size);
 931 
 932     return cmsBuildSegmentedToneCurve(ContextID, 1, &Seg0);
 933 }
 934 
 935 
 936 
 937 // Build a gamma table based on gamma constant
 938 cmsToneCurve* CMSEXPORT cmsBuildGamma(cmsContext ContextID, cmsFloat64Number Gamma)
 939 {
 940     return cmsBuildParametricToneCurve(ContextID, 1, &Gamma);
 941 }
 942 
 943 
 944 // Free all memory taken by the gamma curve
 945 void CMSEXPORT cmsFreeToneCurve(cmsToneCurve* Curve)
 946 {
 947     cmsContext ContextID;
 948 
 949     if (Curve == NULL) return;
 950 
 951     ContextID = Curve ->InterpParams->ContextID;
 952 
 953     _cmsFreeInterpParams(Curve ->InterpParams);
 954 
 955     if (Curve -> Table16)
 956         _cmsFree(ContextID, Curve ->Table16);
 957 
 958     if (Curve ->Segments) {
 959 
 960         cmsUInt32Number i;
 961 
 962         for (i=0; i < Curve ->nSegments; i++) {
 963 
 964             if (Curve ->Segments[i].SampledPoints) {
 965                 _cmsFree(ContextID, Curve ->Segments[i].SampledPoints);
 966             }
 967 
 968             if (Curve ->SegInterp[i] != 0)
 969                 _cmsFreeInterpParams(Curve->SegInterp[i]);
 970         }
 971 
 972         _cmsFree(ContextID, Curve ->Segments);
 973         _cmsFree(ContextID, Curve ->SegInterp);
 974     }
 975 
 976     if (Curve -> Evals)
 977         _cmsFree(ContextID, Curve -> Evals);
 978 
 979     _cmsFree(ContextID, Curve);
 980 }
 981 
 982 // Utility function, free 3 gamma tables
 983 void CMSEXPORT cmsFreeToneCurveTriple(cmsToneCurve* Curve[3])
 984 {
 985 
 986     _cmsAssert(Curve != NULL);
 987 
 988     if (Curve[0] != NULL) cmsFreeToneCurve(Curve[0]);
 989     if (Curve[1] != NULL) cmsFreeToneCurve(Curve[1]);
 990     if (Curve[2] != NULL) cmsFreeToneCurve(Curve[2]);
 991 
 992     Curve[0] = Curve[1] = Curve[2] = NULL;
 993 }
 994 
 995 
 996 // Duplicate a gamma table
 997 cmsToneCurve* CMSEXPORT cmsDupToneCurve(const cmsToneCurve* In)
 998 {
 999     if (In == NULL) return NULL;
1000 
1001     return  AllocateToneCurveStruct(In ->InterpParams ->ContextID, In ->nEntries, In ->nSegments, In ->Segments, In ->Table16);
1002 }
1003 
1004 // Joins two curves for X and Y. Curves should be monotonic.
1005 // We want to get
1006 //
1007 //      y = Y^-1(X(t))
1008 //
1009 cmsToneCurve* CMSEXPORT cmsJoinToneCurve(cmsContext ContextID,
1010                                       const cmsToneCurve* X,
1011                                       const cmsToneCurve* Y, cmsUInt32Number nResultingPoints)
1012 {
1013     cmsToneCurve* out = NULL;
1014     cmsToneCurve* Yreversed = NULL;
1015     cmsFloat32Number t, x;
1016     cmsFloat32Number* Res = NULL;
1017     cmsUInt32Number i;
1018 
1019 
1020     _cmsAssert(X != NULL);
1021     _cmsAssert(Y != NULL);
1022 
1023     Yreversed = cmsReverseToneCurveEx(nResultingPoints, Y);
1024     if (Yreversed == NULL) goto Error;
1025 
1026     Res = (cmsFloat32Number*) _cmsCalloc(ContextID, nResultingPoints, sizeof(cmsFloat32Number));
1027     if (Res == NULL) goto Error;
1028 
1029     //Iterate
1030     for (i=0; i <  nResultingPoints; i++) {
1031 
1032         t = (cmsFloat32Number) i / (cmsFloat32Number)(nResultingPoints-1);
1033         x = cmsEvalToneCurveFloat(X,  t);
1034         Res[i] = cmsEvalToneCurveFloat(Yreversed, x);
1035     }
1036 
1037     // Allocate space for output
1038     out = cmsBuildTabulatedToneCurveFloat(ContextID, nResultingPoints, Res);
1039 
1040 Error:
1041 
1042     if (Res != NULL) _cmsFree(ContextID, Res);
1043     if (Yreversed != NULL) cmsFreeToneCurve(Yreversed);
1044 
1045     return out;
1046 }
1047 
1048 
1049 
1050 // Get the surrounding nodes. This is tricky on non-monotonic tables
1051 static
1052 int GetInterval(cmsFloat64Number In, const cmsUInt16Number LutTable[], const struct _cms_interp_struc* p)
1053 {
1054     int i;
1055     int y0, y1;
1056 
1057     // A 1 point table is not allowed
1058     if (p -> Domain[0] < 1) return -1;
1059 
1060     // Let's see if ascending or descending.
1061     if (LutTable[0] < LutTable[p ->Domain[0]]) {
1062 
1063         // Table is overall ascending
1064         for (i = (int) p->Domain[0] - 1; i >= 0; --i) {
1065 
1066             y0 = LutTable[i];
1067             y1 = LutTable[i+1];
1068 
1069             if (y0 <= y1) { // Increasing
1070                 if (In >= y0 && In <= y1) return i;
1071             }
1072             else
1073                 if (y1 < y0) { // Decreasing
1074                     if (In >= y1 && In <= y0) return i;
1075                 }
1076         }
1077     }
1078     else {
1079         // Table is overall descending
1080         for (i=0; i < (int) p -> Domain[0]; i++) {
1081 
1082             y0 = LutTable[i];
1083             y1 = LutTable[i+1];
1084 
1085             if (y0 <= y1) { // Increasing
1086                 if (In >= y0 && In <= y1) return i;
1087             }
1088             else
1089                 if (y1 < y0) { // Decreasing
1090                     if (In >= y1 && In <= y0) return i;
1091                 }
1092         }
1093     }
1094 
1095     return -1;
1096 }
1097 
1098 // Reverse a gamma table
1099 cmsToneCurve* CMSEXPORT cmsReverseToneCurveEx(cmsUInt32Number nResultSamples, const cmsToneCurve* InCurve)
1100 {
1101     cmsToneCurve *out;
1102     cmsFloat64Number a = 0, b = 0, y, x1, y1, x2, y2;
1103     int i, j;
1104     int Ascending;
1105 
1106     _cmsAssert(InCurve != NULL);
1107 
1108     // Try to reverse it analytically whatever possible
1109 
1110     if (InCurve ->nSegments == 1 && InCurve ->Segments[0].Type > 0 &&
1111         /* InCurve -> Segments[0].Type <= 5 */
1112         GetParametricCurveByType(InCurve ->InterpParams->ContextID, InCurve ->Segments[0].Type, NULL) != NULL) {
1113 
1114         return cmsBuildParametricToneCurve(InCurve ->InterpParams->ContextID,
1115                                        -(InCurve -> Segments[0].Type),
1116                                        InCurve -> Segments[0].Params);
1117     }
1118 
1119     // Nope, reverse the table.
1120     out = cmsBuildTabulatedToneCurve16(InCurve ->InterpParams->ContextID, nResultSamples, NULL);
1121     if (out == NULL)
1122         return NULL;
1123 
1124     // We want to know if this is an ascending or descending table
1125     Ascending = !cmsIsToneCurveDescending(InCurve);
1126 
1127     // Iterate across Y axis
1128     for (i=0; i < (int) nResultSamples; i++) {
1129 
1130         y = (cmsFloat64Number) i * 65535.0 / (nResultSamples - 1);
1131 
1132         // Find interval in which y is within.
1133         j = GetInterval(y, InCurve->Table16, InCurve->InterpParams);
1134         if (j >= 0) {
1135 
1136 
1137             // Get limits of interval
1138             x1 = InCurve ->Table16[j];
1139             x2 = InCurve ->Table16[j+1];
1140 
1141             y1 = (cmsFloat64Number) (j * 65535.0) / (InCurve ->nEntries - 1);
1142             y2 = (cmsFloat64Number) ((j+1) * 65535.0 ) / (InCurve ->nEntries - 1);
1143 
1144             // If collapsed, then use any
1145             if (x1 == x2) {
1146 
1147                 out ->Table16[i] = _cmsQuickSaturateWord(Ascending ? y2 : y1);
1148                 continue;
1149 
1150             } else {
1151 
1152                 // Interpolate
1153                 a = (y2 - y1) / (x2 - x1);
1154                 b = y2 - a * x2;
1155             }
1156         }
1157 
1158         out ->Table16[i] = _cmsQuickSaturateWord(a* y + b);
1159     }
1160 
1161 
1162     return out;
1163 }
1164 
1165 // Reverse a gamma table
1166 cmsToneCurve* CMSEXPORT cmsReverseToneCurve(const cmsToneCurve* InGamma)
1167 {
1168     _cmsAssert(InGamma != NULL);
1169 
1170     return cmsReverseToneCurveEx(4096, InGamma);
1171 }
1172 
1173 // From: Eilers, P.H.C. (1994) Smoothing and interpolation with finite
1174 // differences. in: Graphic Gems IV, Heckbert, P.S. (ed.), Academic press.
1175 //
1176 // Smoothing and interpolation with second differences.
1177 //
1178 //   Input:  weights (w), data (y): vector from 1 to m.
1179 //   Input:  smoothing parameter (lambda), length (m).
1180 //   Output: smoothed vector (z): vector from 1 to m.
1181 
1182 static
1183 cmsBool smooth2(cmsContext ContextID, cmsFloat32Number w[], cmsFloat32Number y[],
1184                 cmsFloat32Number z[], cmsFloat32Number lambda, int m)
1185 {
1186     int i, i1, i2;
1187     cmsFloat32Number *c, *d, *e;
1188     cmsBool st;
1189 
1190 
1191     c = (cmsFloat32Number*) _cmsCalloc(ContextID, MAX_NODES_IN_CURVE, sizeof(cmsFloat32Number));
1192     d = (cmsFloat32Number*) _cmsCalloc(ContextID, MAX_NODES_IN_CURVE, sizeof(cmsFloat32Number));
1193     e = (cmsFloat32Number*) _cmsCalloc(ContextID, MAX_NODES_IN_CURVE, sizeof(cmsFloat32Number));
1194 
1195     if (c != NULL && d != NULL && e != NULL) {
1196 
1197 
1198     d[1] = w[1] + lambda;
1199     c[1] = -2 * lambda / d[1];
1200     e[1] = lambda /d[1];
1201     z[1] = w[1] * y[1];
1202     d[2] = w[2] + 5 * lambda - d[1] * c[1] *  c[1];
1203     c[2] = (-4 * lambda - d[1] * c[1] * e[1]) / d[2];
1204     e[2] = lambda / d[2];
1205     z[2] = w[2] * y[2] - c[1] * z[1];
1206 
1207     for (i = 3; i < m - 1; i++) {
1208         i1 = i - 1; i2 = i - 2;
1209         d[i]= w[i] + 6 * lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
1210         c[i] = (-4 * lambda -d[i1] * c[i1] * e[i1])/ d[i];
1211         e[i] = lambda / d[i];
1212         z[i] = w[i] * y[i] - c[i1] * z[i1] - e[i2] * z[i2];
1213     }
1214 
1215     i1 = m - 2; i2 = m - 3;
1216 
1217     d[m - 1] = w[m - 1] + 5 * lambda -c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
1218     c[m - 1] = (-2 * lambda - d[i1] * c[i1] * e[i1]) / d[m - 1];
1219     z[m - 1] = w[m - 1] * y[m - 1] - c[i1] * z[i1] - e[i2] * z[i2];
1220     i1 = m - 1; i2 = m - 2;
1221 
1222     d[m] = w[m] + lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
1223     z[m] = (w[m] * y[m] - c[i1] * z[i1] - e[i2] * z[i2]) / d[m];
1224     z[m - 1] = z[m - 1] / d[m - 1] - c[m - 1] * z[m];
1225 
1226     for (i = m - 2; 1<= i; i--)
1227         z[i] = z[i] / d[i] - c[i] * z[i + 1] - e[i] * z[i + 2];
1228 
1229       st = TRUE;
1230     }
1231     else st = FALSE;
1232 
1233     if (c != NULL) _cmsFree(ContextID, c);
1234     if (d != NULL) _cmsFree(ContextID, d);
1235     if (e != NULL) _cmsFree(ContextID, e);
1236 
1237     return st;
1238 }
1239 
1240 // Smooths a curve sampled at regular intervals.
1241 cmsBool  CMSEXPORT cmsSmoothToneCurve(cmsToneCurve* Tab, cmsFloat64Number lambda)
1242 {
1243     cmsBool SuccessStatus = TRUE;
1244     cmsFloat32Number *w, *y, *z;
1245     cmsUInt32Number i, nItems, Zeros, Poles;
1246     cmsBool notCheck = FALSE;
1247 
1248     if (Tab != NULL && Tab->InterpParams != NULL)
1249     {
1250         cmsContext ContextID = Tab->InterpParams->ContextID;
1251 
1252         if (!cmsIsToneCurveLinear(Tab)) // Only non-linear curves need smoothing
1253         {
1254             nItems = Tab->nEntries;
1255             if (nItems < MAX_NODES_IN_CURVE)
1256             {
1257                 // Allocate one more item than needed
1258                 w = (cmsFloat32Number *)_cmsCalloc(ContextID, nItems + 1, sizeof(cmsFloat32Number));
1259                 y = (cmsFloat32Number *)_cmsCalloc(ContextID, nItems + 1, sizeof(cmsFloat32Number));
1260                 z = (cmsFloat32Number *)_cmsCalloc(ContextID, nItems + 1, sizeof(cmsFloat32Number));
1261 
1262                 if (w != NULL && y != NULL && z != NULL) // Ensure no memory allocation failure
1263                 {
1264                     memset(w, 0, (nItems + 1) * sizeof(cmsFloat32Number));
1265                     memset(y, 0, (nItems + 1) * sizeof(cmsFloat32Number));
1266                     memset(z, 0, (nItems + 1) * sizeof(cmsFloat32Number));
1267 
1268                     for (i = 0; i < nItems; i++)
1269                     {
1270                         y[i + 1] = (cmsFloat32Number)Tab->Table16[i];
1271                         w[i + 1] = 1.0;
1272                     }
1273 
1274                     if (lambda < 0)
1275                     {
1276                         notCheck = TRUE;
1277                         lambda = -lambda;
1278                     }
1279 
1280                     if (smooth2(ContextID, w, y, z, (cmsFloat32Number)lambda, (int)nItems))
1281                     {
1282                         // Do some reality - checking...
1283 
1284                         Zeros = Poles = 0;
1285                         for (i = nItems; i > 1; --i)
1286                         {
1287                             if (z[i] == 0.) Zeros++;
1288                             if (z[i] >= 65535.) Poles++;
1289                             if (z[i] < z[i - 1])
1290                             {
1291                                 cmsSignalError(ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Non-Monotonic.");
1292                                 SuccessStatus = notCheck;
1293                                 break;
1294                             }
1295                         }
1296 
1297                         if (SuccessStatus && Zeros > (nItems / 3))
1298                         {
1299                             cmsSignalError(ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Degenerated, mostly zeros.");
1300                             SuccessStatus = notCheck;
1301                         }
1302 
1303                         if (SuccessStatus && Poles > (nItems / 3))
1304                         {
1305                             cmsSignalError(ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Degenerated, mostly poles.");
1306                             SuccessStatus = notCheck;
1307                         }
1308 
1309                         if (SuccessStatus) // Seems ok
1310                         {
1311                             for (i = 0; i < nItems; i++)
1312                             {
1313                                 // Clamp to cmsUInt16Number
1314                                 Tab->Table16[i] = _cmsQuickSaturateWord(z[i + 1]);
1315                             }
1316                         }
1317                     }
1318                     else // Could not smooth
1319                     {
1320                         cmsSignalError(ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Function smooth2 failed.");
1321                         SuccessStatus = FALSE;
1322                     }
1323                 }
1324                 else // One or more buffers could not be allocated
1325                 {
1326                     cmsSignalError(ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Could not allocate memory.");
1327                     SuccessStatus = FALSE;
1328                 }
1329 
1330                 if (z != NULL)
1331                     _cmsFree(ContextID, z);
1332 
1333                 if (y != NULL)
1334                     _cmsFree(ContextID, y);
1335 
1336                 if (w != NULL)
1337                     _cmsFree(ContextID, w);
1338             }
1339             else // too many items in the table
1340             {
1341                 cmsSignalError(ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Too many points.");
1342                 SuccessStatus = FALSE;
1343             }
1344         }
1345     }
1346     else // Tab parameter or Tab->InterpParams is NULL
1347     {
1348         // Can't signal an error here since the ContextID is not known at this point
1349         SuccessStatus = FALSE;
1350     }
1351 
1352     return SuccessStatus;
1353 }
1354 
1355 // Is a table linear? Do not use parametric since we cannot guarantee some weird parameters resulting
1356 // in a linear table. This way assures it is linear in 12 bits, which should be enough in most cases.
1357 cmsBool CMSEXPORT cmsIsToneCurveLinear(const cmsToneCurve* Curve)
1358 {
1359     int i;
1360     int diff;
1361 
1362     _cmsAssert(Curve != NULL);
1363 
1364     for (i=0; i < (int) Curve ->nEntries; i++) {
1365 
1366         diff = abs((int) Curve->Table16[i] - (int) _cmsQuantizeVal(i, Curve ->nEntries));
1367         if (diff > 0x0f)
1368             return FALSE;
1369     }
1370 
1371     return TRUE;
1372 }
1373 
1374 // Same, but for monotonicity
1375 cmsBool  CMSEXPORT cmsIsToneCurveMonotonic(const cmsToneCurve* t)
1376 {
1377     cmsUInt32Number n;
1378     int i, last;
1379     cmsBool lDescending;
1380 
1381     _cmsAssert(t != NULL);
1382 
1383     // Degenerated curves are monotonic? Ok, let's pass them
1384     n = t ->nEntries;
1385     if (n < 2) return TRUE;
1386 
1387     // Curve direction
1388     lDescending = cmsIsToneCurveDescending(t);
1389 
1390     if (lDescending) {
1391 
1392         last = t ->Table16[0];
1393 
1394         for (i = 1; i < (int) n; i++) {
1395 
1396             if (t ->Table16[i] - last > 2) // We allow some ripple
1397                 return FALSE;
1398             else
1399                 last = t ->Table16[i];
1400 
1401         }
1402     }
1403     else {
1404 
1405         last = t ->Table16[n-1];
1406 
1407         for (i = (int) n - 2; i >= 0; --i) {
1408 
1409             if (t ->Table16[i] - last > 2)
1410                 return FALSE;
1411             else
1412                 last = t ->Table16[i];
1413 
1414         }
1415     }
1416 
1417     return TRUE;
1418 }
1419 
1420 // Same, but for descending tables
1421 cmsBool  CMSEXPORT cmsIsToneCurveDescending(const cmsToneCurve* t)
1422 {
1423     _cmsAssert(t != NULL);
1424 
1425     return t ->Table16[0] > t ->Table16[t ->nEntries-1];
1426 }
1427 
1428 
1429 // Another info fn: is out gamma table multisegment?
1430 cmsBool  CMSEXPORT cmsIsToneCurveMultisegment(const cmsToneCurve* t)
1431 {
1432     _cmsAssert(t != NULL);
1433 
1434     return t -> nSegments > 1;
1435 }
1436 
1437 cmsInt32Number  CMSEXPORT cmsGetToneCurveParametricType(const cmsToneCurve* t)
1438 {
1439     _cmsAssert(t != NULL);
1440 
1441     if (t -> nSegments != 1) return 0;
1442     return t ->Segments[0].Type;
1443 }
1444 
1445 // We need accuracy this time
1446 cmsFloat32Number CMSEXPORT cmsEvalToneCurveFloat(const cmsToneCurve* Curve, cmsFloat32Number v)
1447 {
1448     _cmsAssert(Curve != NULL);
1449 
1450     // Check for 16 bits table. If so, this is a limited-precision tone curve
1451     if (Curve ->nSegments == 0) {
1452 
1453         cmsUInt16Number In, Out;
1454 
1455         In = (cmsUInt16Number) _cmsQuickSaturateWord(v * 65535.0);
1456         Out = cmsEvalToneCurve16(Curve, In);
1457 
1458         return (cmsFloat32Number) (Out / 65535.0);
1459     }
1460 
1461     return (cmsFloat32Number) EvalSegmentedFn(Curve, v);
1462 }
1463 
1464 // We need xput over here
1465 cmsUInt16Number CMSEXPORT cmsEvalToneCurve16(const cmsToneCurve* Curve, cmsUInt16Number v)
1466 {
1467     cmsUInt16Number out;
1468 
1469     _cmsAssert(Curve != NULL);
1470 
1471     Curve ->InterpParams ->Interpolation.Lerp16(&v, &out, Curve ->InterpParams);
1472     return out;
1473 }
1474 
1475 
1476 // Least squares fitting.
1477 // A mathematical procedure for finding the best-fitting curve to a given set of points by
1478 // minimizing the sum of the squares of the offsets ("the residuals") of the points from the curve.
1479 // The sum of the squares of the offsets is used instead of the offset absolute values because
1480 // this allows the residuals to be treated as a continuous differentiable quantity.
1481 //
1482 // y = f(x) = x ^ g
1483 //
1484 // R  = (yi - (xi^g))
1485 // R2 = (yi - (xi^g))2
1486 // SUM R2 = SUM (yi - (xi^g))2
1487 //
1488 // dR2/dg = -2 SUM x^g log(x)(y - x^g)
1489 // solving for dR2/dg = 0
1490 //
1491 // g = 1/n * SUM(log(y) / log(x))
1492 
1493 cmsFloat64Number CMSEXPORT cmsEstimateGamma(const cmsToneCurve* t, cmsFloat64Number Precision)
1494 {
1495     cmsFloat64Number gamma, sum, sum2;
1496     cmsFloat64Number n, x, y, Std;
1497     cmsUInt32Number i;
1498 
1499     _cmsAssert(t != NULL);
1500 
1501     sum = sum2 = n = 0;
1502 
1503     // Excluding endpoints
1504     for (i=1; i < (MAX_NODES_IN_CURVE-1); i++) {
1505 
1506         x = (cmsFloat64Number) i / (MAX_NODES_IN_CURVE-1);
1507         y = (cmsFloat64Number) cmsEvalToneCurveFloat(t, (cmsFloat32Number) x);
1508 
1509         // Avoid 7% on lower part to prevent
1510         // artifacts due to linear ramps
1511 
1512         if (y > 0. && y < 1. && x > 0.07) {
1513 
1514             gamma = log(y) / log(x);
1515             sum  += gamma;
1516             sum2 += gamma * gamma;
1517             n++;
1518         }
1519     }
1520 
1521     // We need enough valid samples
1522     if (n <= 1) return -1.0;
1523 
1524     // Take a look on SD to see if gamma isn't exponential at all
1525     Std = sqrt((n * sum2 - sum * sum) / (n*(n-1)));
1526 
1527     if (Std > Precision)
1528         return -1.0;
1529 
1530     return (sum / n);   // The mean
1531 }
1532 
1533 // Retrieve segments on tone curves
1534 
1535 const cmsCurveSegment* CMSEXPORT cmsGetToneCurveSegment(cmsInt32Number n, const cmsToneCurve* t)
1536 {
1537     _cmsAssert(t != NULL);
1538 
1539     if (n < 0 || n >= (cmsInt32Number) t->nSegments) return NULL;
1540     return t->Segments + n;
1541 }
1542