modules/lcms/src/cmsgamma.c

   1 //
   2 //  Little cms
   3 //  Copyright (C) 1998-2007 Marti Maria
   4 //
   5 // Permission is hereby granted, free of charge, to any person obtaining
   6 // a copy of this software and associated documentation files (the "Software"),
   7 // to deal in the Software without restriction, including without limitation
   8 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
   9 // and/or sell copies of the Software, and to permit persons to whom the Software
  10 // is furnished to do so, subject to the following conditions:
  11 //
  12 // The above copyright notice and this permission notice shall be included in
  13 // all copies or substantial portions of the Software.
  14 //
  15 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  16 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
  17 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  18 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  19 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  20 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  21 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  22
  23
  24 #include "lcms.h"
  25
  26 // Gamma handling.
  27
  28 LPGAMMATABLE LCMSEXPORT cmsAllocGamma(int nEntries);
  29 void         LCMSEXPORT cmsFreeGamma(LPGAMMATABLE Gamma);
  30 void         LCMSEXPORT cmsFreeGammaTriple(LPGAMMATABLE Gamma[3]);
  31 LPGAMMATABLE LCMSEXPORT cmsBuildGamma(int nEntries, double Gamma);
  32 LPGAMMATABLE LCMSEXPORT cmsDupGamma(LPGAMMATABLE Src);
  33 LPGAMMATABLE LCMSEXPORT cmsReverseGamma(int nResultSamples, LPGAMMATABLE InGamma);
  34 LPGAMMATABLE LCMSEXPORT cmsBuildParametricGamma(int nEntries, int Type, double Params[]);
  35 LPGAMMATABLE LCMSEXPORT cmsJoinGamma(LPGAMMATABLE InGamma, LPGAMMATABLE OutGamma);
  36 LPGAMMATABLE LCMSEXPORT cmsJoinGammaEx(LPGAMMATABLE InGamma, LPGAMMATABLE OutGamma, int nPoints);
  37 LCMSBOOL         LCMSEXPORT cmsSmoothGamma(LPGAMMATABLE Tab, double lambda);
  38
  39 LCMSBOOL         cdecl _cmsSmoothEndpoints(LPWORD Table, int nPoints);
  40
  41
  42 // Sampled curves
  43
  44 LPSAMPLEDCURVE cdecl cmsAllocSampledCurve(int nItems);
  45 void           cdecl cmsFreeSampledCurve(LPSAMPLEDCURVE p);
  46 void           cdecl cmsEndpointsOfSampledCurve(LPSAMPLEDCURVE p, double* Min, double* Max);
  47 void           cdecl cmsClampSampledCurve(LPSAMPLEDCURVE p, double Min, double Max);
  48 LCMSBOOL       cdecl cmsSmoothSampledCurve(LPSAMPLEDCURVE Tab, double SmoothingLambda);
  49 void           cdecl cmsRescaleSampledCurve(LPSAMPLEDCURVE p, double Min, double Max, int nPoints);
  50
  51 LPSAMPLEDCURVE cdecl cmsJoinSampledCurves(LPSAMPLEDCURVE X, LPSAMPLEDCURVE Y, int nResultingPoints);
  52
  53 double LCMSEXPORT cmsEstimateGamma(LPGAMMATABLE t);
  54 double LCMSEXPORT cmsEstimateGammaEx(LPWORD GammaTable, int nEntries, double Thereshold);
  55
  56 // ----------------------------------------------------------------------------------------
  57
  58
  59 #define MAX_KNOTS   4096
  60 typedef float vec[MAX_KNOTS+1];
  61
  62
  63 // Ciclic-redundant-check for assuring table is a true representation of parametric curve
  64
  65 // The usual polynomial, which is used for AAL5, FDDI, and probably Ethernet
  66 #define QUOTIENT 0x04c11db7
  67
  68 static
  69 unsigned int Crc32(unsigned int result, LPVOID ptr, int len)
  70 {
  71     int          i,j;
  72     BYTE         octet;
  73     LPBYTE       data = (LPBYTE) ptr;
  74
  75     for (i=0; i < len; i++) {
  76
  77         octet = *data++;
  78
  79         for (j=0; j < 8; j++) {
  80
  81             if (result & 0x80000000) {
  82
  83                 result = (result << 1) ^ QUOTIENT ^ (octet >> 7);
  84             }
  85             else
  86             {
  87                 result = (result << 1) ^ (octet >> 7);
  88             }
  89             octet <<= 1;
  90         }
  91     }
  92
  93     return result;
  94 }
  95
  96 // Get CRC of gamma table
  97
  98 unsigned int _cmsCrc32OfGammaTable(LPGAMMATABLE Table)
  99 {
 100     unsigned int crc = ~0U;
 101
 102     crc = Crc32(crc, &Table -> Seed.Type,  sizeof(int));
 103     crc = Crc32(crc, Table ->Seed.Params,  sizeof(double)*10);
 104     crc = Crc32(crc, &Table ->nEntries,    sizeof(int));
 105     crc = Crc32(crc, Table ->GammaTable,   sizeof(WORD) * Table -> nEntries);
 106
 107     return ~crc;
 108
 109 }
 110
 111
 112 LPGAMMATABLE LCMSEXPORT cmsAllocGamma(int nEntries)
 113 {
 114        LPGAMMATABLE p;
 115        size_t size;
 116
 117        if (nEntries > 65530 || nEntries < 0) {
 118                 cmsSignalError(LCMS_ERRC_ABORTED, "Couldn't create gammatable of more than 65530 entries");
 119                 return NULL;
 120        }
 121
 122        size = sizeof(GAMMATABLE) + (sizeof(WORD) * (nEntries-1));
 123
 124        p = (LPGAMMATABLE) _cmsMalloc(size);
 125        if (!p) return NULL;
 126
 127        ZeroMemory(p, size);
 128
 129        p -> Seed.Type     = 0;
 130        p -> nEntries = nEntries;
 131
 132        return p;
 133 }
 134
 135 void LCMSEXPORT cmsFreeGamma(LPGAMMATABLE Gamma)
 136 {
 137        if (Gamma)  _cmsFree(Gamma);
 138 }
 139
 140
 141
 142 void LCMSEXPORT cmsFreeGammaTriple(LPGAMMATABLE Gamma[3])
 143 {
 144     cmsFreeGamma(Gamma[0]);
 145     cmsFreeGamma(Gamma[1]);
 146     cmsFreeGamma(Gamma[2]);
 147     Gamma[0] = Gamma[1] = Gamma[2] = NULL;
 148 }
 149
 150
 151
 152 // Duplicate a gamma table
 153
 154 LPGAMMATABLE  LCMSEXPORT cmsDupGamma(LPGAMMATABLE In)
 155 {
 156        LPGAMMATABLE Ptr;
 157        size_t size;
 158
 159        Ptr = cmsAllocGamma(In -> nEntries);
 160        if (Ptr == NULL) return NULL;
 161
 162        size = sizeof(GAMMATABLE) + (sizeof(WORD) * (In -> nEntries-1));
 163
 164        CopyMemory(Ptr, In, size);
 165        return Ptr;
 166 }
 167
 168
 169 // Handle gamma using interpolation tables. The resulting curves can become
 170 // very stange, but are pleasent to eye.
 171
 172 LPGAMMATABLE LCMSEXPORT cmsJoinGamma(LPGAMMATABLE InGamma,
 173                           LPGAMMATABLE OutGamma)
 174 {
 175        register int i;
 176        L16PARAMS L16In, L16Out;
 177        LPWORD InPtr, OutPtr;
 178        LPGAMMATABLE p;
 179
 180        p = cmsAllocGamma(256);
 181        if (!p) return NULL;
 182
 183        cmsCalcL16Params(InGamma -> nEntries, &L16In);
 184        InPtr  = InGamma -> GammaTable;
 185
 186        cmsCalcL16Params(OutGamma -> nEntries, &L16Out);
 187        OutPtr = OutGamma-> GammaTable;
 188
 189        for (i=0; i < 256; i++)
 190        {
 191               WORD wValIn, wValOut;
 192
 193               wValIn  = cmsLinearInterpLUT16(RGB_8_TO_16(i), InPtr, &L16In);
 194               wValOut = cmsReverseLinearInterpLUT16(wValIn, OutPtr, &L16Out);
 195
 196               p -> GammaTable[i] = wValOut;
 197        }
 198
 199        return p;
 200 }
 201
 202
 203
 204 // New method, using smoothed parametric curves. This works FAR better.
 205 // We want to get
 206 //
 207 //      y = f(g^-1(x))      ; f = ingamma, g = outgamma
 208 //
 209 // And this can be parametrized as
 210 //
 211 //      y = f(t)
 212 //      x = g(t)
 213
 214
 215 LPGAMMATABLE LCMSEXPORT cmsJoinGammaEx(LPGAMMATABLE InGamma,
 216                                        LPGAMMATABLE OutGamma, int nPoints)
 217 {
 218
 219     LPSAMPLEDCURVE x, y, r;
 220     LPGAMMATABLE res;
 221
 222     x = cmsConvertGammaToSampledCurve(InGamma,  nPoints);
 223     y = cmsConvertGammaToSampledCurve(OutGamma, nPoints);
 224     r = cmsJoinSampledCurves(y, x, nPoints);
 225
 226     // Does clean "hair"
 227     cmsSmoothSampledCurve(r, 0.001);
 228
 229     cmsClampSampledCurve(r, 0.0, 65535.0);
 230
 231     cmsFreeSampledCurve(x);
 232     cmsFreeSampledCurve(y);
 233
 234     res = cmsConvertSampledCurveToGamma(r, 65535.0);
 235     cmsFreeSampledCurve(r);
 236
 237     return res;
 238 }
 239
 240
 241
 242 // Reverse a gamma table
 243
 244 LPGAMMATABLE LCMSEXPORT cmsReverseGamma(int nResultSamples, LPGAMMATABLE InGamma)
 245 {
 246        register int i;
 247        L16PARAMS L16In;
 248        LPWORD InPtr;
 249        LPGAMMATABLE p;
 250
 251        // Try to reverse it analytically whatever possible
 252        if (InGamma -> Seed.Type > 0 && InGamma -> Seed.Type <= 5 &&
 253                         _cmsCrc32OfGammaTable(InGamma) == InGamma -> Seed.Crc32) {
 254
 255                 return cmsBuildParametricGamma(nResultSamples, -(InGamma -> Seed.Type), InGamma ->Seed.Params);
 256        }
 257
 258
 259        // Nope, reverse the table
 260        p = cmsAllocGamma(nResultSamples);
 261        if (!p) return NULL;
 262
 263        cmsCalcL16Params(InGamma -> nEntries, &L16In);
 264        InPtr  = InGamma -> GammaTable;
 265
 266        for (i=0; i < nResultSamples; i++)
 267        {
 268               WORD wValIn, wValOut;
 269
 270               wValIn = _cmsQuantizeVal(i, nResultSamples);
 271               wValOut = cmsReverseLinearInterpLUT16(wValIn, InPtr, &L16In);
 272               p -> GammaTable[i] = wValOut;
 273        }
 274
 275
 276        return p;
 277 }
 278
 279
 280
 281 // Parametric curves
 282 //
 283 // Parameters goes as: Gamma, a, b, c, d, e, f
 284 // Type is the ICC type +1
 285 // if type is negative, then the curve is analyticaly inverted
 286
 287 LPGAMMATABLE LCMSEXPORT cmsBuildParametricGamma(int nEntries, int Type, double Params[])
 288 {
 289         LPGAMMATABLE Table;
 290         double R, Val, dval, e;
 291         int i;
 292         int ParamsByType[] = { 0, 1, 3, 4, 5, 7 };
 293
 294         Table = cmsAllocGamma(nEntries);
 295         if (NULL == Table) return NULL;
 296
 297         Table -> Seed.Type = Type;
 298
 299         CopyMemory(Table ->Seed.Params, Params, ParamsByType[abs(Type)] * sizeof(double));
 300
 301
 302         for (i=0; i < nEntries; i++) {
 303
 304                 R   = (double) i / (nEntries-1);
 305
 306                 switch (Type) {
 307
 308                 // X = Y ^ Gamma
 309                 case 1:
 310                       Val = pow(R, Params[0]);
 311                       break;
 312
 313                 // Type 1 Reversed: X = Y ^1/gamma
 314                 case -1:
 315                       Val = pow(R, 1/Params[0]);
 316                       break;
 317
 318                 // CIE 122-1966
 319                 // Y = (aX + b)^Gamma  | X >= -b/a
 320                 // Y = 0               | else
 321                 case 2:
 322                     if (R >= -Params[2] / Params[1]) {
 323
 324                               e = Params[1]*R + Params[2];
 325
 326                               if (e > 0)
 327                                 Val = pow(e, Params[0]);
 328                               else
 329                                 Val = 0;
 330                     }
 331                     else
 332                               Val = 0;
 333                       break;
 334
 335                 // Type 2 Reversed
 336                 // X = (Y ^1/g  - b) / a
 337                 case -2:
 338
 339                     Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1];
 340                     if (Val < 0)
 341                             Val = 0;
 342                     break;
 343
 344
 345                 // IEC 61966-3
 346                 // Y = (aX + b)^Gamma | X <= -b/a
 347                 // Y = c              | else
 348                 case 3:
 349                     if (R >= -Params[2] / Params[1]) {
 350
 351                       e = Params[1]*R + Params[2];
 352                       Val = pow(e, Params[0]) + Params[3];
 353                     }
 354                     else
 355                       Val = Params[3];
 356                     break;
 357
 358
 359                 // Type 3 reversed
 360                 // X=((Y-c)^1/g - b)/a      | (Y>=c)
 361                 // X=-b/a                   | (Y<c)
 362
 363                 case -3:
 364                     if (R >= Params[3])  {
 365                         e = R - Params[3];
 366                         Val = (pow(e, 1/Params[0]) - Params[2]) / Params[1];
 367                         if (Val < 0) Val = 0;
 368                     }
 369                     else {
 370                         Val = -Params[2] / Params[1];
 371                     }
 372                     break;
 373
 374
 375                 // IEC 61966-2.1 (sRGB)
 376                 // Y = (aX + b)^Gamma | X >= d
 377                 // Y = cX             | X < d
 378                 case 4:
 379                     if (R >= Params[4]) {
 380
 381                               e = Params[1]*R + Params[2];
 382                               if (e > 0)
 383                                 Val = pow(e, Params[0]);
 384                               else
 385                                 Val = 0;
 386                     }
 387                       else
 388                               Val = R * Params[3];
 389                       break;
 390
 391                 // Type 4 reversed
 392                 // X=((Y^1/g-b)/a)    | Y >= (ad+b)^g
 393                 // X=Y/c              | Y< (ad+b)^g
 394
 395                 case -4:
 396                     if (R >= pow(Params[1] * Params[4] + Params[2], Params[0])) {
 397
 398                         Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1];
 399                     }
 400                     else {
 401                         Val = R / Params[3];
 402                     }
 403                     break;
 404
 405
 406
 407                 // Y = (aX + b)^Gamma + e | X <= d
 408                 // Y = cX + f             | else
 409                 case 5:
 410                     if (R >= Params[4]) {
 411
 412                         e = Params[1]*R + Params[2];
 413                         Val = pow(e, Params[0]) + Params[5];
 414                     }
 415                     else
 416                         Val = R*Params[3] + Params[6];
 417                     break;
 418
 419
 420                 // Reversed type 5
 421                 // X=((Y-e)1/g-b)/a   | Y >=(ad+b)^g+e)
 422                 // X=(Y-f)/c          | else
 423                 case -5:
 424
 425                 if (R >= pow(Params[1] * Params[4], Params[0]) + Params[5]) {
 426
 427                     Val = pow(R - Params[5], 1/Params[0]) - Params[2] / Params[1];
 428                 }
 429                 else {
 430                     Val = (R - Params[6]) / Params[3];
 431                 }
 432                 break;
 433
 434                 default:
 435                         cmsSignalError(LCMS_ERRC_ABORTED, "Unsupported parametric curve type=%d", abs(Type)-1);
 436                         cmsFreeGamma(Table);
 437                         return NULL;
 438                 }
 439
 440
 441         // Saturate
 442
 443         dval = Val * 65535.0 + .5;
 444         if (dval > 65535.) dval = 65535.0;
 445         if (dval < 0) dval = 0;
 446
 447         Table->GammaTable[i] = (WORD) floor(dval);
 448         }
 449
 450         Table -> Seed.Crc32 = _cmsCrc32OfGammaTable(Table);
 451
 452         return Table;
 453 }
 454
 455 // Build a gamma table based on gamma constant
 456
 457 LPGAMMATABLE LCMSEXPORT cmsBuildGamma(int nEntries, double Gamma)
 458 {
 459     return cmsBuildParametricGamma(nEntries, 1, &Gamma);
 460 }
 461
 462
 463
 464 // From: Eilers, P.H.C. (1994) Smoothing and interpolation with finite
 465 // differences. in: Graphic Gems IV, Heckbert, P.S. (ed.), Academic press.
 466 //
 467 // Smoothing and interpolation with second differences.
 468 //
 469 //   Input:  weights (w), data (y): vector from 1 to m.
 470 //   Input:  smoothing parameter (lambda), length (m).
 471 //   Output: smoothed vector (z): vector from 1 to m.
 472
 473
 474 static
 475 void smooth2(vec w, vec y, vec z, float lambda, int m)
 476 {
 477   int i, i1, i2;
 478   vec c, d, e;
 479   d[1] = w[1] + lambda;
 480   c[1] = -2 * lambda / d[1];
 481   e[1] = lambda /d[1];
 482   z[1] = w[1] * y[1];
 483   d[2] = w[2] + 5 * lambda - d[1] * c[1] *  c[1];
 484   c[2] = (-4 * lambda - d[1] * c[1] * e[1]) / d[2];
 485   e[2] = lambda / d[2];
 486   z[2] = w[2] * y[2] - c[1] * z[1];
 487   for (i = 3; i < m - 1; i++) {
 488     i1 = i - 1; i2 = i - 2;
 489     d[i]= w[i] + 6 * lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
 490     c[i] = (-4 * lambda -d[i1] * c[i1] * e[i1])/ d[i];
 491     e[i] = lambda / d[i];
 492     z[i] = w[i] * y[i] - c[i1] * z[i1] - e[i2] * z[i2];
 493   }
 494   i1 = m - 2; i2 = m - 3;
 495   d[m - 1] = w[m - 1] + 5 * lambda -c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
 496   c[m - 1] = (-2 * lambda - d[i1] * c[i1] * e[i1]) / d[m - 1];
 497   z[m - 1] = w[m - 1] * y[m - 1] - c[i1] * z[i1] - e[i2] * z[i2];
 498   i1 = m - 1; i2 = m - 2;
 499   d[m] = w[m] + lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2];
 500   z[m] = (w[m] * y[m] - c[i1] * z[i1] - e[i2] * z[i2]) / d[m];
 501   z[m - 1] = z[m - 1] / d[m - 1] - c[m - 1] * z[m];
 502   for (i = m - 2; 1<= i; i--)
 503      z[i] = z[i] / d[i] - c[i] * z[i + 1] - e[i] * z[i + 2];
 504 }
 505
 506
 507
 508 // Smooths a curve sampled at regular intervals
 509
 510 LCMSBOOL LCMSEXPORT cmsSmoothGamma(LPGAMMATABLE Tab, double lambda)
 511
 512 {
 513     vec w, y, z;
 514     int i, nItems, Zeros, Poles;
 515
 516
 517     if (cmsIsLinear(Tab->GammaTable, Tab->nEntries)) return FALSE; // Nothing to do
 518
 519     nItems = Tab -> nEntries;
 520
 521     if (nItems > MAX_KNOTS) {
 522                 cmsSignalError(LCMS_ERRC_ABORTED, "cmsSmoothGamma: too many points.");
 523                 return FALSE;
 524                 }
 525
 526     ZeroMemory(w, nItems * sizeof(float));
 527     ZeroMemory(y, nItems * sizeof(float));
 528     ZeroMemory(z, nItems * sizeof(float));
 529
 530     for (i=0; i < nItems; i++)
 531     {
 532         y[i+1] = (float) Tab -> GammaTable[i];
 533         w[i+1] = 1.0;
 534     }
 535
 536     smooth2(w, y, z, (float) lambda, nItems);
 537
 538     // Do some reality - checking...
 539     Zeros = Poles = 0;
 540     for (i=nItems; i > 1; --i) {
 541
 542             if (z[i] == 0.) Zeros++;
 543             if (z[i] >= 65535.) Poles++;
 544             if (z[i] < z[i-1]) return FALSE; // Non-Monotonic
 545     }
 546
 547     if (Zeros > (nItems / 3)) return FALSE;  // Degenerated, mostly zeros
 548     if (Poles > (nItems / 3)) return FALSE;  // Degenerated, mostly poles
 549
 550     // Seems ok
 551
 552     for (i=0; i < nItems; i++) {
 553
 554         // Clamp to WORD
 555
 556         float v = z[i+1];
 557
 558         if (v < 0) v = 0;
 559         if (v > 65535.) v = 65535.;
 560
 561         Tab -> GammaTable[i] = (WORD) floor(v + .5);
 562         }
 563
 564     return TRUE;
 565 }
 566
 567
 568 // Check if curve is exponential, return gamma if so.
 569
 570 double LCMSEXPORT cmsEstimateGammaEx(LPWORD GammaTable, int nEntries, double Thereshold)
 571 {
 572     double gamma, sum, sum2;
 573     double n, x, y, Std;
 574     int i;
 575
 576     sum = sum2 = n = 0;
 577
 578     // Does exclude endpoints
 579     for (i=1; i < nEntries - 1; i++) {
 580
 581             x = (double) i / (nEntries - 1);
 582             y = (double) GammaTable[i] / 65535.;
 583
 584             // Avoid 7% on lower part to prevent
 585             // artifacts due to linear ramps
 586
 587             if (y > 0. && y < 1. && x > 0.07) {
 588
 589             gamma = log(y) / log(x);
 590             sum  += gamma;
 591             sum2 += gamma * gamma;
 592             n++;
 593             }
 594
 595     }
 596
 597     // Take a look on SD to see if gamma isn't exponential at all
 598     Std = sqrt((n * sum2 - sum * sum) / (n*(n-1)));
 599
 600
 601     if (Std > Thereshold)
 602         return -1.0;
 603
 604     return (sum / n);   // The mean
 605 }
 606
 607
 608 double LCMSEXPORT cmsEstimateGamma(LPGAMMATABLE t)
 609 {
 610         return cmsEstimateGammaEx(t->GammaTable, t->nEntries, 0.7);
 611 }
 612
 613
 614 // -----------------------------------------------------------------Sampled curves
 615
 616 // Allocate a empty curve
 617
 618 LPSAMPLEDCURVE cmsAllocSampledCurve(int nItems)
 619 {
 620     LPSAMPLEDCURVE pOut;
 621
 622     pOut = (LPSAMPLEDCURVE) _cmsMalloc(sizeof(SAMPLEDCURVE));
 623     if (pOut == NULL)
 624             return NULL;
 625
 626     if((pOut->Values = (double *) _cmsMalloc(nItems * sizeof(double))) == NULL)
 627     {
 628          _cmsFree(pOut);
 629         return NULL;
 630     }
 631
 632     pOut->nItems = nItems;
 633     ZeroMemory(pOut->Values, nItems * sizeof(double));
 634
 635     return pOut;
 636 }
 637
 638
 639 void cmsFreeSampledCurve(LPSAMPLEDCURVE p)
 640 {
 641      _cmsFree((LPVOID) p -> Values);
 642      _cmsFree((LPVOID) p);
 643 }
 644
 645
 646
 647 // Does duplicate a sampled curve
 648
 649 LPSAMPLEDCURVE cmsDupSampledCurve(LPSAMPLEDCURVE p)
 650 {
 651     LPSAMPLEDCURVE out;
 652
 653     out = cmsAllocSampledCurve(p -> nItems);
 654     if (!out) return NULL;
 655
 656     CopyMemory(out ->Values, p ->Values, p->nItems * sizeof(double));
 657
 658     return out;
 659 }
 660
 661
 662 // Take min, max of curve
 663
 664 void cmsEndpointsOfSampledCurve(LPSAMPLEDCURVE p, double* Min, double* Max)
 665 {
 666         int i;
 667
 668         *Min = 65536.;
 669         *Max = 0.;
 670
 671         for (i=0; i < p -> nItems; i++) {
 672
 673                 double v = p -> Values[i];
 674
 675                 if (v < *Min)
 676                         *Min = v;
 677
 678                 if (v > *Max)
 679                         *Max = v;
 680         }
 681
 682         if (*Min < 0) *Min = 0;
 683         if (*Max > 65535.0) *Max = 65535.0;
 684 }
 685
 686 // Clamps to Min, Max
 687
 688 void cmsClampSampledCurve(LPSAMPLEDCURVE p, double Min, double Max)
 689 {
 690
 691         int i;
 692
 693         for (i=0; i < p -> nItems; i++) {
 694
 695                 double v = p -> Values[i];
 696
 697                 if (v < Min)
 698                         v = Min;
 699
 700                 if (v > Max)
 701                         v = Max;
 702
 703                 p -> Values[i] = v;
 704
 705         }
 706
 707 }
 708
 709
 710
 711 // Smooths a curve sampled at regular intervals
 712
 713 LCMSBOOL cmsSmoothSampledCurve(LPSAMPLEDCURVE Tab, double lambda)
 714 {
 715     vec w, y, z;
 716     int i, nItems;
 717
 718     nItems = Tab -> nItems;
 719
 720     if (nItems > MAX_KNOTS) {
 721                 cmsSignalError(LCMS_ERRC_ABORTED, "cmsSmoothSampledCurve: too many points.");
 722                 return FALSE;
 723                 }
 724
 725     ZeroMemory(w, nItems * sizeof(float));
 726     ZeroMemory(y, nItems * sizeof(float));
 727     ZeroMemory(z, nItems * sizeof(float));
 728
 729     for (i=0; i < nItems; i++)
 730     {
 731         float value = (float) Tab -> Values[i];
 732
 733         y[i+1] = value;
 734         w[i+1] = (float) ((value < 0.0) ?  0 : 1);
 735     }
 736
 737
 738     smooth2(w, y, z, (float) lambda, nItems);
 739
 740     for (i=0; i < nItems; i++) {
 741
 742         Tab -> Values[i] = z[i+1];;
 743      }
 744
 745     return TRUE;
 746
 747 }
 748
 749
 750 // Scale a value v, within domain Min .. Max
 751 // to a domain 0..(nPoints-1)
 752
 753 static
 754 double ScaleVal(double v, double Min, double Max, int nPoints)
 755 {
 756
 757         double a, b;
 758
 759         if (v <= Min) return 0;
 760         if (v >= Max) return (nPoints-1);
 761
 762         a = (double) (nPoints - 1) / (Max - Min);
 763         b = a * Min;
 764
 765         return (a * v) - b;
 766
 767 }
 768
 769
 770 // Does rescale a sampled curve to fit in a 0..(nPoints-1) domain
 771
 772 void cmsRescaleSampledCurve(LPSAMPLEDCURVE p, double Min, double Max, int nPoints)
 773 {
 774
 775         int i;
 776
 777         for (i=0; i < p -> nItems; i++) {
 778
 779                 double v = p -> Values[i];
 780
 781                 p -> Values[i] = ScaleVal(v, Min, Max, nPoints);
 782         }
 783
 784 }
 785
 786
 787 // Joins two sampled curves for X and Y. Curves should be sorted.
 788
 789 LPSAMPLEDCURVE cmsJoinSampledCurves(LPSAMPLEDCURVE X, LPSAMPLEDCURVE Y, int nResultingPoints)
 790 {
 791     int i, j;
 792     LPSAMPLEDCURVE out;
 793     double MinX, MinY, MaxX, MaxY;
 794     double x, y, x1, y1, x2, y2, a, b;
 795
 796     out = cmsAllocSampledCurve(nResultingPoints);
 797     if (out == NULL)
 798         return NULL;
 799
 800     if (X -> nItems != Y -> nItems) {
 801
 802         cmsSignalError(LCMS_ERRC_ABORTED, "cmsJoinSampledCurves: invalid curve.");
 803         cmsFreeSampledCurve(out);
 804         return NULL;
 805     }
 806
 807     // Get endpoints of sampled curves
 808     cmsEndpointsOfSampledCurve(X, &MinX, &MaxX);
 809     cmsEndpointsOfSampledCurve(Y, &MinY, &MaxY);
 810
 811
 812     // Set our points
 813     out ->Values[0] = MinY;
 814     for (i=1; i < nResultingPoints; i++) {
 815
 816         // Scale t to x domain
 817         x = (i * (MaxX - MinX) / (nResultingPoints-1)) + MinX;
 818
 819         // Find interval in which j is within (always up,
 820         // since fn should be monotonic at all)
 821
 822         j = 1;
 823         while ((j < X ->nItems - 1) && X ->Values[j] < x)
 824             j++;
 825
 826         // Now x is within X[j-1], X[j]
 827         x1 = X ->Values[j-1]; x2 = X ->Values[j];
 828         y1 = Y ->Values[j-1]; y2 = Y ->Values[j];
 829
 830         // Interpolate  the value
 831         a = (y1 - y2) / (x1 - x2);
 832         b = y1 - a * x1;
 833         y = a* x + b;
 834
 835         out ->Values[i] = y;
 836     }
 837
 838
 839     cmsClampSampledCurve(out, MinY, MaxY);
 840     return out;
 841 }
 842
 843
 844
 845 // Convert between curve types
 846
 847 LPGAMMATABLE cmsConvertSampledCurveToGamma(LPSAMPLEDCURVE Sampled, double Max)
 848 {
 849     LPGAMMATABLE Gamma;
 850     int i, nPoints;
 851
 852
 853     nPoints = Sampled ->nItems;
 854
 855     Gamma = cmsAllocGamma(nPoints);
 856     for (i=0; i < nPoints; i++) {
 857
 858         Gamma->GammaTable[i] = (WORD) floor(ScaleVal(Sampled ->Values[i], 0, Max, 65536) + .5);
 859     }
 860
 861     return Gamma;
 862
 863 }
 864
 865 // Inverse of anterior
 866
 867 LPSAMPLEDCURVE cmsConvertGammaToSampledCurve(LPGAMMATABLE Gamma, int nPoints)
 868 {
 869     LPSAMPLEDCURVE Sampled;
 870     L16PARAMS L16;
 871     int i;
 872     WORD wQuant, wValIn;
 873
 874     if (nPoints > 4096) {
 875
 876         cmsSignalError(LCMS_ERRC_ABORTED, "cmsConvertGammaToSampledCurve: too many points (max=4096)");
 877         return NULL;
 878     }
 879
 880     cmsCalcL16Params(Gamma -> nEntries, &L16);
 881
 882     Sampled = cmsAllocSampledCurve(nPoints);
 883     for (i=0; i < nPoints; i++) {
 884             wQuant  = _cmsQuantizeVal(i, nPoints);
 885             wValIn  = cmsLinearInterpLUT16(wQuant, Gamma ->GammaTable, &L16);
 886             Sampled ->Values[i] = (float) wValIn;
 887     }
 888
 889     return Sampled;
 890 }
 891
 892
 893
 894
 895 // Smooth endpoints (used in Black/White compensation)
 896
 897 LCMSBOOL _cmsSmoothEndpoints(LPWORD Table, int nEntries)
 898 {
 899     vec w, y, z;
 900     int i, Zeros, Poles;
 901
 902
 903
 904     if (cmsIsLinear(Table, nEntries)) return FALSE; // Nothing to do
 905
 906
 907     if (nEntries > MAX_KNOTS) {
 908                 cmsSignalError(LCMS_ERRC_ABORTED, "_cmsSmoothEndpoints: too many points.");
 909                 return FALSE;
 910                 }
 911
 912     ZeroMemory(w, nEntries * sizeof(float));
 913     ZeroMemory(y, nEntries * sizeof(float));
 914     ZeroMemory(z, nEntries * sizeof(float));
 915
 916     for (i=0; i < nEntries; i++)
 917     {
 918         y[i+1] = (float) Table[i];
 919         w[i+1] = 1.0;
 920     }
 921
 922     w[1]        = 65535.0;
 923     w[nEntries] = 65535.0;
 924
 925     smooth2(w, y, z, (float) nEntries, nEntries);
 926
 927     // Do some reality - checking...
 928     Zeros = Poles = 0;
 929     for (i=nEntries; i > 1; --i) {
 930
 931             if (z[i] == 0.) Zeros++;
 932             if (z[i] >= 65535.) Poles++;
 933             if (z[i] < z[i-1]) return FALSE; // Non-Monotonic
 934     }
 935
 936     if (Zeros > (nEntries / 3)) return FALSE;  // Degenerated, mostly zeros
 937     if (Poles > (nEntries / 3)) return FALSE;    // Degenerated, mostly poles
 938
 939     // Seems ok
 940
 941     for (i=0; i < nEntries; i++) {
 942
 943         // Clamp to WORD
 944
 945         float v = z[i+1];
 946
 947         if (v < 0) v = 0;
 948         if (v > 65535.) v = 65535.;
 949
 950         Table[i] = (WORD) floor(v + .5);
 951         }
 952
 953     return TRUE;
 954 }