src/cmsintrp.c

   1 //---------------------------------------------------------------------------------
   2 //
   3 //  Little Color Management System
   4 //  Copyright (c) 1998-2010 Marti Maria Saguer
   5 //
   6 // Permission is hereby granted, free of charge, to any person obtaining
   7 // a copy of this software and associated documentation files (the "Software"),
   8 // to deal in the Software without restriction, including without limitation
   9 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  10 // and/or sell copies of the Software, and to permit persons to whom the Software
  11 // is furnished to do so, subject to the following conditions:
  12 //
  13 // The above copyright notice and this permission notice shall be included in
  14 // all copies or substantial portions of the Software.
  15 //
  16 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
  18 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  20 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  21 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  22 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  23 //
  24 //---------------------------------------------------------------------------------
  25 //
  26
  27 #include "lcms2_internal.h"
  28
  29 // This module incorporates several interpolation routines, for 1 to 8 channels on input and
  30 // up to 65535 channels on output. The user may change those by using the interpolation plug-in
  31
  32 // Interpolation routines by default
  33 static cmsInterpFunction DefaultInterpolatorsFactory(cmsUInt32Number nInputChannels, cmsUInt32Number nOutputChannels, cmsUInt32Number dwFlags);
  34
  35 // This is the default factory
  36 static cmsInterpFnFactory Interpolators = DefaultInterpolatorsFactory;
  37
  38
  39 // Main plug-in entry
  40 cmsBool  _cmsRegisterInterpPlugin(cmsPluginBase* Data)
  41 {
  42     cmsPluginInterpolation* Plugin = (cmsPluginInterpolation*) Data;
  43
  44     if (Data == NULL) {
  45
  46         Interpolators = DefaultInterpolatorsFactory;
  47         return TRUE;
  48     }
  49
  50     // Set replacement functions
  51     Interpolators = Plugin ->InterpolatorsFactory;
  52     return TRUE;
  53 }
  54
  55
  56 // Set the interpolation method
  57
  58 cmsBool _cmsSetInterpolationRoutine(cmsInterpParams* p)
  59 {
  60     // Invoke factory, possibly in the Plug-in
  61     p ->Interpolation = Interpolators(p -> nInputs, p ->nOutputs, p ->dwFlags);
  62
  63     // If unsupported by the plug-in, go for the LittleCMS default.
  64     // If happens only if an extern plug-in is being used
  65     if (p ->Interpolation.Lerp16 == NULL)
  66         p ->Interpolation = DefaultInterpolatorsFactory(p ->nInputs, p ->nOutputs, p ->dwFlags);
  67
  68     // Check for valid interpolator (we just check one member of the union)
  69     if (p ->Interpolation.Lerp16 == NULL) {
  70             return FALSE;
  71     }
  72     return TRUE;
  73 }
  74
  75
  76 // This function precalculates as many parameters as possible to speed up the interpolation.
  77 cmsInterpParams* _cmsComputeInterpParamsEx(cmsContext ContextID,
  78                                            const cmsUInt32Number nSamples[],
  79                                            int InputChan, int OutputChan,
  80                                            const void *Table,
  81                                            cmsUInt32Number dwFlags)
  82 {
  83     cmsInterpParams* p;
  84     int i;
  85
  86     // Check for maximum inputs
  87     if (InputChan > MAX_INPUT_DIMENSIONS) {
  88              cmsSignalError(ContextID, cmsERROR_RANGE, "Too many input channels (%d channels, max=%d)", InputChan, MAX_INPUT_DIMENSIONS);
  89             return NULL;
  90     }
  91
  92     // Creates an empty object
  93     p = (cmsInterpParams*) _cmsMallocZero(ContextID, sizeof(cmsInterpParams));
  94     if (p == NULL) return NULL;
  95
  96     // Keep original parameters
  97     p -> dwFlags  = dwFlags;
  98     p -> nInputs  = InputChan;
  99     p -> nOutputs = OutputChan;
 100     p ->Table     = Table;
 101     p ->ContextID  = ContextID;
 102
 103     // Fill samples per input direction and domain (which is number of nodes minus one)
 104     for (i=0; i < InputChan; i++) {
 105
 106         p -> nSamples[i] = nSamples[i];
 107         p -> Domain[i]   = nSamples[i] - 1;
 108     }
 109
 110     // Compute factors to apply to each component to index the grid array
 111     p -> opta[0] = p -> nOutputs;
 112     for (i=1; i < InputChan; i++)
 113         p ->opta[i] = p ->opta[i-1] * nSamples[InputChan-i];
 114
 115
 116     if (!_cmsSetInterpolationRoutine(p)) {
 117          cmsSignalError(ContextID, cmsERROR_UNKNOWN_EXTENSION, "Unsupported interpolation (%d->%d channels)", InputChan, OutputChan);
 118         _cmsFree(ContextID, p);
 119         return NULL;
 120     }
 121
 122     // All seems ok
 123     return p;
 124 }
 125
 126
 127 // This one is a wrapper on the anterior, but assuming all directions have same number of nodes
 128 cmsInterpParams* _cmsComputeInterpParams(cmsContext ContextID, int nSamples, int InputChan, int OutputChan, const void* Table, cmsUInt32Number dwFlags)
 129 {
 130     int i;
 131     cmsUInt32Number Samples[MAX_INPUT_DIMENSIONS];
 132
 133     // Fill the auxiliar array
 134     for (i=0; i < MAX_INPUT_DIMENSIONS; i++)
 135         Samples[i] = nSamples;
 136
 137     // Call the extended function
 138     return _cmsComputeInterpParamsEx(ContextID, Samples, InputChan, OutputChan, Table, dwFlags);
 139 }
 140
 141
 142 // Free all associated memory
 143 void _cmsFreeInterpParams(cmsInterpParams* p)
 144 {
 145     if (p != NULL) _cmsFree(p ->ContextID, p);
 146 }
 147
 148
 149 // Inline fixed point interpolation
 150 cmsINLINE cmsUInt16Number LinearInterp(cmsS15Fixed16Number a, cmsS15Fixed16Number l, cmsS15Fixed16Number h)
 151 {
 152     cmsUInt32Number dif = (cmsUInt32Number) (h - l) * a + 0x8000;
 153     dif = (dif >> 16) + l;
 154     return (cmsUInt16Number) (dif);
 155 }
 156
 157
 158 //  Linear interpolation (Fixed-point optimized)
 159 static
 160 void LinLerp1D(register const cmsUInt16Number Value[],
 161                register cmsUInt16Number Output[],
 162                register const cmsInterpParams* p)
 163 {
 164     cmsUInt16Number y1, y0;
 165     int cell0, rest;
 166     int val3;
 167     const cmsUInt16Number* LutTable = (cmsUInt16Number*) p ->Table;
 168
 169     // if last value...
 170     if (Value[0] == 0xffff) {
 171
 172         Output[0] = LutTable[p -> Domain[0]];
 173         return;
 174     }
 175
 176     val3 = p -> Domain[0] * Value[0];
 177     val3 = _cmsToFixedDomain(val3);    // To fixed 15.16
 178
 179     cell0 = FIXED_TO_INT(val3);             // Cell is 16 MSB bits
 180     rest  = FIXED_REST_TO_INT(val3);        // Rest is 16 LSB bits
 181
 182     y0 = LutTable[cell0];
 183     y1 = LutTable[cell0+1];
 184
 185
 186     Output[0] = LinearInterp(rest, y0, y1);
 187 }
 188
 189
 190 // Floating-point version of 1D interpolation
 191 static
 192 void LinLerp1Dfloat(const cmsFloat32Number Value[],
 193                     cmsFloat32Number Output[],
 194                     const cmsInterpParams* p)
 195 {
 196        cmsFloat32Number y1, y0;
 197        cmsFloat32Number val2, rest;
 198        int cell0, cell1;
 199        const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table;
 200
 201        // if last value...
 202        if (Value[0] == 1.0) {
 203            Output[0] = LutTable[p -> Domain[0]];
 204            return;
 205        }
 206
 207        val2 = p -> Domain[0] * Value[0];
 208
 209        cell0 = (int) floor(val2);
 210        cell1 = (int) ceil(val2);
 211
 212        // Rest is 16 LSB bits
 213        rest = val2 - cell0;
 214
 215        y0 = LutTable[cell0] ;
 216        y1 = LutTable[cell1] ;
 217
 218        Output[0] = y0 + (y1 - y0) * rest;
 219 }
 220
 221
 222
 223 // Eval gray LUT having only one input channel
 224 static
 225 void Eval1Input(register const cmsUInt16Number Input[],
 226                 register cmsUInt16Number Output[],
 227                 register const cmsInterpParams* p16)
 228 {
 229        cmsS15Fixed16Number fk;
 230        cmsS15Fixed16Number k0, k1, rk, K0, K1;
 231        int v;
 232        cmsUInt32Number OutChan;
 233        const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;
 234
 235        v = Input[0] * p16 -> Domain[0];
 236        fk = _cmsToFixedDomain(v);
 237
 238        k0 = FIXED_TO_INT(fk);
 239        rk = (cmsUInt16Number) FIXED_REST_TO_INT(fk);
 240
 241        k1 = k0 + (Input[0] != 0xFFFFU ? 1 : 0);
 242
 243        K0 = p16 -> opta[0] * k0;
 244        K1 = p16 -> opta[0] * k1;
 245
 246        for (OutChan=0; OutChan < p16->nOutputs; OutChan++) {
 247
 248            Output[OutChan] = LinearInterp(rk, LutTable[K0+OutChan], LutTable[K1+OutChan]);
 249        }
 250 }
 251
 252
 253
 254 // Eval gray LUT having only one input channel
 255 static
 256 void Eval1InputFloat(const cmsFloat32Number Value[],
 257                      cmsFloat32Number Output[],
 258                      const cmsInterpParams* p)
 259 {
 260     cmsFloat32Number y1, y0;
 261     cmsFloat32Number val2, rest;
 262     int cell0, cell1;
 263     cmsUInt32Number OutChan;
 264     const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table;
 265
 266         // if last value...
 267        if (Value[0] == 1.0) {
 268            Output[0] = LutTable[p -> Domain[0]];
 269            return;
 270        }
 271
 272        val2 = p -> Domain[0] * Value[0];
 273
 274        cell0 = (int) floor(val2);
 275        cell1 = (int) ceil(val2);
 276
 277        // Rest is 16 LSB bits
 278        rest = val2 - cell0;
 279
 280        cell0 *= p -> opta[0];
 281        cell1 *= p -> opta[0];
 282
 283        for (OutChan=0; OutChan < p->nOutputs; OutChan++) {
 284
 285             y0 = LutTable[cell0 + OutChan] ;
 286             y1 = LutTable[cell1 + OutChan] ;
 287
 288             Output[OutChan] = y0 + (y1 - y0) * rest;
 289        }
 290 }
 291
 292 // Bilinear interpolation (16 bits) - cmsFloat32Number version
 293 static
 294 void BilinearInterpFloat(const cmsFloat32Number Input[],
 295                          cmsFloat32Number Output[],
 296                          const cmsInterpParams* p)
 297
 298 {
 299 #   define LERP(a,l,h)    (cmsFloat32Number) ((l)+(((h)-(l))*(a)))
 300 #   define DENS(i,j)      (LutTable[(i)+(j)+OutChan])
 301
 302     const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table;
 303     cmsFloat32Number      px, py;
 304     int        x0, y0,
 305                X0, Y0, X1, Y1;
 306     int        TotalOut, OutChan;
 307     cmsFloat32Number      fx, fy,
 308         d00, d01, d10, d11,
 309         dx0, dx1,
 310         dxy;
 311
 312     TotalOut   = p -> nOutputs;
 313     px = Input[0] * p->Domain[0];
 314     py = Input[1] * p->Domain[1];
 315
 316     x0 = (int) _cmsQuickFloor(px); fx = px - (cmsFloat32Number) x0;
 317     y0 = (int) _cmsQuickFloor(py); fy = py - (cmsFloat32Number) y0;
 318
 319     X0 = p -> opta[1] * x0;
 320     X1 = X0 + (Input[0] >= 1.0 ? 0 : p->opta[1]);
 321
 322     Y0 = p -> opta[0] * y0;
 323     Y1 = Y0 + (Input[1] >= 1.0 ? 0 : p->opta[0]);
 324
 325     for (OutChan = 0; OutChan < TotalOut; OutChan++) {
 326
 327         d00 = DENS(X0, Y0);
 328         d01 = DENS(X0, Y1);
 329         d10 = DENS(X1, Y0);
 330         d11 = DENS(X1, Y1);
 331
 332         dx0 = LERP(fx, d00, d10);
 333         dx1 = LERP(fx, d01, d11);
 334
 335         dxy = LERP(fy, dx0, dx1);
 336
 337         Output[OutChan] = dxy;
 338     }
 339
 340
 341 #   undef LERP
 342 #   undef DENS
 343 }
 344
 345 // Bilinear interpolation (16 bits) - optimized version
 346 static
 347 void BilinearInterp16(register const cmsUInt16Number Input[],
 348                       register cmsUInt16Number Output[],
 349                       register const cmsInterpParams* p)
 350
 351 {
 352 #define DENS(i,j) (LutTable[(i)+(j)+OutChan])
 353 #define LERP(a,l,h)     (cmsUInt16Number) (l + ROUND_FIXED_TO_INT(((h-l)*a)))
 354
 355            const cmsUInt16Number* LutTable = (cmsUInt16Number*) p ->Table;
 356            int        OutChan, TotalOut;
 357            cmsS15Fixed16Number    fx, fy;
 358   register int        rx, ry;
 359            int        x0, y0;
 360   register int        X0, X1, Y0, Y1;
 361            int        d00, d01, d10, d11,
 362                       dx0, dx1,
 363                       dxy;
 364
 365     TotalOut   = p -> nOutputs;
 366
 367     fx = _cmsToFixedDomain((int) Input[0] * p -> Domain[0]);
 368     x0  = FIXED_TO_INT(fx);
 369     rx  = FIXED_REST_TO_INT(fx);    // Rest in 0..1.0 domain
 370
 371
 372     fy = _cmsToFixedDomain((int) Input[1] * p -> Domain[1]);
 373     y0  = FIXED_TO_INT(fy);
 374     ry  = FIXED_REST_TO_INT(fy);
 375
 376
 377     X0 = p -> opta[1] * x0;
 378     X1 = X0 + (Input[0] == 0xFFFFU ? 0 : p->opta[1]);
 379
 380     Y0 = p -> opta[0] * y0;
 381     Y1 = Y0 + (Input[1] == 0xFFFFU ? 0 : p->opta[0]);
 382
 383     for (OutChan = 0; OutChan < TotalOut; OutChan++) {
 384
 385         d00 = DENS(X0, Y0);
 386         d01 = DENS(X0, Y1);
 387         d10 = DENS(X1, Y0);
 388         d11 = DENS(X1, Y1);
 389
 390         dx0 = LERP(rx, d00, d10);
 391         dx1 = LERP(rx, d01, d11);
 392
 393         dxy = LERP(ry, dx0, dx1);
 394
 395         Output[OutChan] = (cmsUInt16Number) dxy;
 396     }
 397
 398
 399 #   undef LERP
 400 #   undef DENS
 401 }
 402
 403
 404 // Trilinear interpolation (16 bits) - cmsFloat32Number version
 405 static
 406 void TrilinearInterpFloat(const cmsFloat32Number Input[],
 407                           cmsFloat32Number Output[],
 408                           const cmsInterpParams* p)
 409
 410 {
 411 #   define LERP(a,l,h)      (cmsFloat32Number) ((l)+(((h)-(l))*(a)))
 412 #   define DENS(i,j,k)      (LutTable[(i)+(j)+(k)+OutChan])
 413
 414     const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table;
 415     cmsFloat32Number      px, py, pz;
 416     int        x0, y0, z0,
 417                X0, Y0, Z0, X1, Y1, Z1;
 418     int        TotalOut, OutChan;
 419     cmsFloat32Number      fx, fy, fz,
 420         d000, d001, d010, d011,
 421         d100, d101, d110, d111,
 422         dx00, dx01, dx10, dx11,
 423         dxy0, dxy1, dxyz;
 424
 425     TotalOut   = p -> nOutputs;
 426
 427     // We need some clipping here
 428     px = Input[0];
 429     py = Input[1];
 430     pz = Input[2];
 431
 432     if (px < 0) px = 0;
 433     if (px > 1) px = 1;
 434     if (py < 0) py = 0;
 435     if (py > 1) py = 1;
 436     if (pz < 0) pz = 0;
 437     if (pz > 1) pz = 1;
 438
 439     px *= p->Domain[0];
 440     py *= p->Domain[1];
 441     pz *= p->Domain[2];
 442
 443     x0 = (int) _cmsQuickFloor(px); fx = px - (cmsFloat32Number) x0;
 444     y0 = (int) _cmsQuickFloor(py); fy = py - (cmsFloat32Number) y0;
 445     z0 = (int) _cmsQuickFloor(pz); fz = pz - (cmsFloat32Number) z0;
 446
 447     X0 = p -> opta[2] * x0;
 448     X1 = X0 + (Input[0] >= 1.0 ? 0 : p->opta[2]);
 449
 450     Y0 = p -> opta[1] * y0;
 451     Y1 = Y0 + (Input[1] >= 1.0 ? 0 : p->opta[1]);
 452
 453     Z0 = p -> opta[0] * z0;
 454     Z1 = Z0 + (Input[2] >= 1.0 ? 0 : p->opta[0]);
 455
 456     for (OutChan = 0; OutChan < TotalOut; OutChan++) {
 457
 458         d000 = DENS(X0, Y0, Z0);
 459         d001 = DENS(X0, Y0, Z1);
 460         d010 = DENS(X0, Y1, Z0);
 461         d011 = DENS(X0, Y1, Z1);
 462
 463         d100 = DENS(X1, Y0, Z0);
 464         d101 = DENS(X1, Y0, Z1);
 465         d110 = DENS(X1, Y1, Z0);
 466         d111 = DENS(X1, Y1, Z1);
 467
 468
 469         dx00 = LERP(fx, d000, d100);
 470         dx01 = LERP(fx, d001, d101);
 471         dx10 = LERP(fx, d010, d110);
 472         dx11 = LERP(fx, d011, d111);
 473
 474         dxy0 = LERP(fy, dx00, dx10);
 475         dxy1 = LERP(fy, dx01, dx11);
 476
 477         dxyz = LERP(fz, dxy0, dxy1);
 478
 479         Output[OutChan] = dxyz;
 480     }
 481
 482
 483 #   undef LERP
 484 #   undef DENS
 485 }
 486
 487 // Trilinear interpolation (16 bits) - optimized version
 488 static
 489 void TrilinearInterp16(register const cmsUInt16Number Input[],
 490                        register cmsUInt16Number Output[],
 491                        register const cmsInterpParams* p)
 492
 493 {
 494 #define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
 495 #define LERP(a,l,h)     (cmsUInt16Number) (l + ROUND_FIXED_TO_INT(((h-l)*a)))
 496
 497            const cmsUInt16Number* LutTable = (cmsUInt16Number*) p ->Table;
 498            int        OutChan, TotalOut;
 499            cmsS15Fixed16Number    fx, fy, fz;
 500   register int        rx, ry, rz;
 501            int        x0, y0, z0;
 502   register int        X0, X1, Y0, Y1, Z0, Z1;
 503            int        d000, d001, d010, d011,
 504                       d100, d101, d110, d111,
 505                       dx00, dx01, dx10, dx11,
 506                       dxy0, dxy1, dxyz;
 507
 508     TotalOut   = p -> nOutputs;
 509
 510     fx = _cmsToFixedDomain((int) Input[0] * p -> Domain[0]);
 511     x0  = FIXED_TO_INT(fx);
 512     rx  = FIXED_REST_TO_INT(fx);    // Rest in 0..1.0 domain
 513
 514
 515     fy = _cmsToFixedDomain((int) Input[1] * p -> Domain[1]);
 516     y0  = FIXED_TO_INT(fy);
 517     ry  = FIXED_REST_TO_INT(fy);
 518
 519     fz = _cmsToFixedDomain((int) Input[2] * p -> Domain[2]);
 520     z0 = FIXED_TO_INT(fz);
 521     rz = FIXED_REST_TO_INT(fz);
 522
 523
 524     X0 = p -> opta[2] * x0;
 525     X1 = X0 + (Input[0] == 0xFFFFU ? 0 : p->opta[2]);
 526
 527     Y0 = p -> opta[1] * y0;
 528     Y1 = Y0 + (Input[1] == 0xFFFFU ? 0 : p->opta[1]);
 529
 530     Z0 = p -> opta[0] * z0;
 531     Z1 = Z0 + (Input[2] == 0xFFFFU ? 0 : p->opta[0]);
 532
 533     for (OutChan = 0; OutChan < TotalOut; OutChan++) {
 534
 535         d000 = DENS(X0, Y0, Z0);
 536         d001 = DENS(X0, Y0, Z1);
 537         d010 = DENS(X0, Y1, Z0);
 538         d011 = DENS(X0, Y1, Z1);
 539
 540         d100 = DENS(X1, Y0, Z0);
 541         d101 = DENS(X1, Y0, Z1);
 542         d110 = DENS(X1, Y1, Z0);
 543         d111 = DENS(X1, Y1, Z1);
 544
 545
 546         dx00 = LERP(rx, d000, d100);
 547         dx01 = LERP(rx, d001, d101);
 548         dx10 = LERP(rx, d010, d110);
 549         dx11 = LERP(rx, d011, d111);
 550
 551         dxy0 = LERP(ry, dx00, dx10);
 552         dxy1 = LERP(ry, dx01, dx11);
 553
 554         dxyz = LERP(rz, dxy0, dxy1);
 555
 556         Output[OutChan] = (cmsUInt16Number) dxyz;
 557     }
 558
 559
 560 #   undef LERP
 561 #   undef DENS
 562 }
 563
 564
 565 // Tetrahedral interpolation, using Sakamoto algorithm.
 566 #define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
 567 static
 568 void TetrahedralInterpFloat(const cmsFloat32Number Input[],
 569                             cmsFloat32Number Output[],
 570                             const cmsInterpParams* p)
 571 {
 572     const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;
 573     cmsFloat32Number     px, py, pz;
 574     int        x0, y0, z0,
 575                X0, Y0, Z0, X1, Y1, Z1;
 576     cmsFloat32Number     rx, ry, rz;
 577     cmsFloat32Number     c0, c1=0, c2=0, c3=0;
 578     int                  OutChan, TotalOut;
 579
 580     TotalOut   = p -> nOutputs;
 581
 582     // We need some clipping here
 583     px = Input[0];
 584     py = Input[1];
 585     pz = Input[2];
 586
 587     if (px < 0) px = 0;
 588     if (px > 1) px = 1;
 589     if (py < 0) py = 0;
 590     if (py > 1) py = 1;
 591     if (pz < 0) pz = 0;
 592     if (pz > 1) pz = 1;
 593
 594     px *= p->Domain[0];
 595     py *= p->Domain[1];
 596     pz *= p->Domain[2];
 597
 598     x0 = (int) _cmsQuickFloor(px); rx = (px - (cmsFloat32Number) x0);
 599     y0 = (int) _cmsQuickFloor(py); ry = (py - (cmsFloat32Number) y0);
 600     z0 = (int) _cmsQuickFloor(pz); rz = (pz - (cmsFloat32Number) z0);
 601
 602
 603     X0 = p -> opta[2] * x0;
 604     X1 = X0 + (Input[0] >= 1.0 ? 0 : p->opta[2]);
 605
 606     Y0 = p -> opta[1] * y0;
 607     Y1 = Y0 + (Input[1] >= 1.0 ? 0 : p->opta[1]);
 608
 609     Z0 = p -> opta[0] * z0;
 610     Z1 = Z0 + (Input[2] >= 1.0 ? 0 : p->opta[0]);
 611
 612     for (OutChan=0; OutChan < TotalOut; OutChan++) {
 613
 614        // These are the 6 Tetrahedral
 615
 616         c0 = DENS(X0, Y0, Z0);
 617
 618         if (rx >= ry && ry >= rz) {
 619
 620             c1 = DENS(X1, Y0, Z0) - c0;
 621             c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);
 622             c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
 623
 624         }
 625         else
 626             if (rx >= rz && rz >= ry) {
 627
 628                 c1 = DENS(X1, Y0, Z0) - c0;
 629                 c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
 630                 c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);
 631
 632             }
 633             else
 634                 if (rz >= rx && rx >= ry) {
 635
 636                     c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);
 637                     c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
 638                     c3 = DENS(X0, Y0, Z1) - c0;
 639
 640                 }
 641                 else
 642                     if (ry >= rx && rx >= rz) {
 643
 644                         c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);
 645                         c2 = DENS(X0, Y1, Z0) - c0;
 646                         c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
 647
 648                     }
 649                     else
 650                         if (ry >= rz && rz >= rx) {
 651
 652                             c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
 653                             c2 = DENS(X0, Y1, Z0) - c0;
 654                             c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);
 655
 656                         }
 657                         else
 658                             if (rz >= ry && ry >= rx) {
 659
 660                                 c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
 661                                 c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);
 662                                 c3 = DENS(X0, Y0, Z1) - c0;
 663
 664                             }
 665                             else  {
 666                                 c1 = c2 = c3 = 0;
 667                             }
 668
 669        Output[OutChan] = c0 + c1 * rx + c2 * ry + c3 * rz;
 670        }
 671
 672 }
 673
 674 #undef DENS
 675
 676
 677
 678 #define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
 679
 680 static
 681 void TetrahedralInterp16(register const cmsUInt16Number Input[],
 682                          register cmsUInt16Number Output[],
 683                          register const cmsInterpParams* p)
 684 {
 685     const cmsUInt16Number* LutTable = (cmsUInt16Number*) p -> Table;
 686     cmsS15Fixed16Number    fx, fy, fz;
 687     cmsS15Fixed16Number    rx, ry, rz;
 688     int                    x0, y0, z0;
 689     cmsS15Fixed16Number    c0, c1, c2, c3, Rest;
 690     cmsUInt32Number        OutChan;
 691     cmsS15Fixed16Number    X0, X1, Y0, Y1, Z0, Z1;
 692     cmsUInt32Number        TotalOut = p -> nOutputs;
 693
 694
 695     fx  = _cmsToFixedDomain((int) Input[0] * p -> Domain[0]);
 696     fy  = _cmsToFixedDomain((int) Input[1] * p -> Domain[1]);
 697     fz  = _cmsToFixedDomain((int) Input[2] * p -> Domain[2]);
 698
 699     x0  = FIXED_TO_INT(fx);
 700     y0  = FIXED_TO_INT(fy);
 701     z0  = FIXED_TO_INT(fz);
 702
 703     rx  = FIXED_REST_TO_INT(fx);
 704     ry  = FIXED_REST_TO_INT(fy);
 705     rz  = FIXED_REST_TO_INT(fz);
 706
 707     X0 = p -> opta[2] * x0;
 708     X1 = X0 + (Input[0] == 0xFFFFU ? 0 : p->opta[2]);
 709
 710     Y0 = p -> opta[1] * y0;
 711     Y1 = Y0 + (Input[1] == 0xFFFFU ? 0 : p->opta[1]);
 712
 713     Z0 = p -> opta[0] * z0;
 714     Z1 = Z0 + (Input[2] == 0xFFFFU ? 0 : p->opta[0]);
 715
 716     // These are the 6 Tetrahedral
 717     for (OutChan=0; OutChan < TotalOut; OutChan++) {
 718
 719         c0 = DENS(X0, Y0, Z0);
 720
 721         if (rx >= ry && ry >= rz) {
 722
 723             c1 = DENS(X1, Y0, Z0) - c0;
 724             c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);
 725             c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
 726
 727         }
 728         else
 729             if (rx >= rz && rz >= ry) {
 730
 731                 c1 = DENS(X1, Y0, Z0) - c0;
 732                 c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
 733                 c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);
 734
 735             }
 736             else
 737                 if (rz >= rx && rx >= ry) {
 738
 739                     c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);
 740                     c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
 741                     c3 = DENS(X0, Y0, Z1) - c0;
 742
 743                 }
 744                 else
 745                     if (ry >= rx && rx >= rz) {
 746
 747                         c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);
 748                         c2 = DENS(X0, Y1, Z0) - c0;
 749                         c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
 750
 751                     }
 752                     else
 753                         if (ry >= rz && rz >= rx) {
 754
 755                             c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
 756                             c2 = DENS(X0, Y1, Z0) - c0;
 757                             c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);
 758
 759                         }
 760                         else
 761                             if (rz >= ry && ry >= rx) {
 762
 763                                 c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
 764                                 c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);
 765                                 c3 = DENS(X0, Y0, Z1) - c0;
 766
 767                             }
 768                             else  {
 769                                 c1 = c2 = c3 = 0;
 770                             }
 771
 772                             Rest = c1 * rx + c2 * ry + c3 * rz;
 773
 774                             Output[OutChan] = (cmsUInt16Number) c0 + ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest));
 775     }
 776
 777 }
 778 #undef DENS
 779
 780
 781 #define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
 782 static
 783 void Eval4Inputs(register const cmsUInt16Number Input[],
 784                      register cmsUInt16Number Output[],
 785                      register const cmsInterpParams* p16)
 786 {
 787     const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;
 788     cmsS15Fixed16Number fk;
 789     cmsS15Fixed16Number k0, rk;
 790     int K0, K1;
 791     cmsS15Fixed16Number    fx, fy, fz;
 792     cmsS15Fixed16Number    rx, ry, rz;
 793     int                    x0, y0, z0;
 794     cmsS15Fixed16Number    X0, X1, Y0, Y1, Z0, Z1;
 795     cmsUInt32Number i;
 796     cmsS15Fixed16Number    c0, c1, c2, c3, Rest;
 797     cmsUInt32Number        OutChan;
 798     cmsUInt16Number        Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
 799
 800
 801     fk  = _cmsToFixedDomain((int) Input[0] * p16 -> Domain[0]);
 802     fx  = _cmsToFixedDomain((int) Input[1] * p16 -> Domain[1]);
 803     fy  = _cmsToFixedDomain((int) Input[2] * p16 -> Domain[2]);
 804     fz  = _cmsToFixedDomain((int) Input[3] * p16 -> Domain[3]);
 805
 806     k0  = FIXED_TO_INT(fk);
 807     x0  = FIXED_TO_INT(fx);
 808     y0  = FIXED_TO_INT(fy);
 809     z0  = FIXED_TO_INT(fz);
 810
 811     rk  = FIXED_REST_TO_INT(fk);
 812     rx  = FIXED_REST_TO_INT(fx);
 813     ry  = FIXED_REST_TO_INT(fy);
 814     rz  = FIXED_REST_TO_INT(fz);
 815
 816     K0 = p16 -> opta[3] * k0;
 817     K1 = K0 + (Input[0] == 0xFFFFU ? 0 : p16->opta[3]);
 818
 819     X0 = p16 -> opta[2] * x0;
 820     X1 = X0 + (Input[1] == 0xFFFFU ? 0 : p16->opta[2]);
 821
 822     Y0 = p16 -> opta[1] * y0;
 823     Y1 = Y0 + (Input[2] == 0xFFFFU ? 0 : p16->opta[1]);
 824
 825     Z0 = p16 -> opta[0] * z0;
 826     Z1 = Z0 + (Input[3] == 0xFFFFU ? 0 : p16->opta[0]);
 827
 828     LutTable = (cmsUInt16Number*) p16 -> Table;
 829     LutTable += K0;
 830
 831     for (OutChan=0; OutChan < p16 -> nOutputs; OutChan++) {
 832
 833         c0 = DENS(X0, Y0, Z0);
 834
 835         if (rx >= ry && ry >= rz) {
 836
 837             c1 = DENS(X1, Y0, Z0) - c0;
 838             c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);
 839             c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
 840
 841         }
 842         else
 843             if (rx >= rz && rz >= ry) {
 844
 845                 c1 = DENS(X1, Y0, Z0) - c0;
 846                 c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
 847                 c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);
 848
 849             }
 850             else
 851                 if (rz >= rx && rx >= ry) {
 852
 853                     c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);
 854                     c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
 855                     c3 = DENS(X0, Y0, Z1) - c0;
 856
 857                 }
 858                 else
 859                     if (ry >= rx && rx >= rz) {
 860
 861                         c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);
 862                         c2 = DENS(X0, Y1, Z0) - c0;
 863                         c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
 864
 865                     }
 866                     else
 867                         if (ry >= rz && rz >= rx) {
 868
 869                             c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
 870                             c2 = DENS(X0, Y1, Z0) - c0;
 871                             c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);
 872
 873                         }
 874                         else
 875                             if (rz >= ry && ry >= rx) {
 876
 877                                 c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
 878                                 c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);
 879                                 c3 = DENS(X0, Y0, Z1) - c0;
 880
 881                             }
 882                             else  {
 883                                 c1 = c2 = c3 = 0;
 884                             }
 885
 886                             Rest = c1 * rx + c2 * ry + c3 * rz;
 887
 888                             Tmp1[OutChan] = (cmsUInt16Number) c0 + ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest));
 889     }
 890
 891
 892     LutTable = (cmsUInt16Number*) p16 -> Table;
 893     LutTable += K1;
 894
 895     for (OutChan=0; OutChan < p16 -> nOutputs; OutChan++) {
 896
 897         c0 = DENS(X0, Y0, Z0);
 898
 899         if (rx >= ry && ry >= rz) {
 900
 901             c1 = DENS(X1, Y0, Z0) - c0;
 902             c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);
 903             c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
 904
 905         }
 906         else
 907             if (rx >= rz && rz >= ry) {
 908
 909                 c1 = DENS(X1, Y0, Z0) - c0;
 910                 c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
 911                 c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);
 912
 913             }
 914             else
 915                 if (rz >= rx && rx >= ry) {
 916
 917                     c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);
 918                     c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
 919                     c3 = DENS(X0, Y0, Z1) - c0;
 920
 921                 }
 922                 else
 923                     if (ry >= rx && rx >= rz) {
 924
 925                         c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);
 926                         c2 = DENS(X0, Y1, Z0) - c0;
 927                         c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
 928
 929                     }
 930                     else
 931                         if (ry >= rz && rz >= rx) {
 932
 933                             c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
 934                             c2 = DENS(X0, Y1, Z0) - c0;
 935                             c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);
 936
 937                         }
 938                         else
 939                             if (rz >= ry && ry >= rx) {
 940
 941                                 c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
 942                                 c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);
 943                                 c3 = DENS(X0, Y0, Z1) - c0;
 944
 945                             }
 946                             else  {
 947                                 c1 = c2 = c3 = 0;
 948                             }
 949
 950                             Rest = c1 * rx + c2 * ry + c3 * rz;
 951
 952                             Tmp2[OutChan] = (cmsUInt16Number) c0 + ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest));
 953     }
 954
 955
 956
 957     for (i=0; i < p16 -> nOutputs; i++) {
 958         Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);
 959     }
 960 }
 961 #undef DENS
 962
 963
 964 // For more that 3 inputs (i.e., CMYK)
 965 // evaluate two 3-dimensional interpolations and then linearly interpolate between them.
 966
 967
 968 static
 969 void Eval4InputsFloat(const cmsFloat32Number Input[],
 970                       cmsFloat32Number Output[],
 971                       const cmsInterpParams* p)
 972 {
 973        const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;
 974        cmsFloat32Number rest;
 975        cmsFloat32Number pk;
 976        int k0, K0, K1;
 977        const cmsFloat32Number* T;
 978        cmsUInt32Number i;
 979        cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
 980        cmsInterpParams p1;
 981
 982
 983        pk = Input[0] * p->Domain[0];
 984        k0 = _cmsQuickFloor(pk);
 985        rest = pk - (cmsFloat32Number) k0;
 986
 987        K0 = p -> opta[3] * k0;
 988        K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[3]);
 989
 990        p1 = *p;
 991        memmove(&p1.Domain[0], &p ->Domain[1], 3*sizeof(cmsUInt32Number));
 992
 993        T = LutTable + K0;
 994        p1.Table = T;
 995
 996        TetrahedralInterpFloat(Input + 1,  Tmp1, &p1);
 997
 998        T = LutTable + K1;
 999        p1.Table = T;
1000        TetrahedralInterpFloat(Input + 1,  Tmp2, &p1);
1001
1002        for (i=0; i < p -> nOutputs; i++)
1003        {
1004               cmsFloat32Number y0 = Tmp1[i];
1005               cmsFloat32Number y1 = Tmp2[i];
1006
1007               Output[i] = y0 + (y1 - y0) * rest;
1008        }
1009 }
1010
1011
1012 static
1013 void Eval5Inputs(register const cmsUInt16Number Input[],
1014                  register cmsUInt16Number Output[],
1015
1016                  register const cmsInterpParams* p16)
1017 {
1018        const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;
1019        cmsS15Fixed16Number fk;
1020        cmsS15Fixed16Number k0, rk;
1021        int K0, K1;
1022        const cmsUInt16Number* T;
1023        cmsUInt32Number i;
1024        cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
1025        cmsInterpParams p1;
1026
1027
1028        fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);
1029        k0 = FIXED_TO_INT(fk);
1030        rk = FIXED_REST_TO_INT(fk);
1031
1032        K0 = p16 -> opta[4] * k0;
1033        K1 = p16 -> opta[4] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));
1034
1035        p1 = *p16;
1036        memmove(&p1.Domain[0], &p16 ->Domain[1], 4*sizeof(cmsUInt32Number));
1037
1038        T = LutTable + K0;
1039        p1.Table = T;
1040
1041        Eval4Inputs(Input + 1, Tmp1, &p1);
1042
1043        T = LutTable + K1;
1044        p1.Table = T;
1045
1046        Eval4Inputs(Input + 1, Tmp2, &p1);
1047
1048        for (i=0; i < p16 -> nOutputs; i++) {
1049
1050               Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);
1051        }
1052
1053 }
1054
1055
1056 static
1057 void Eval5InputsFloat(const cmsFloat32Number Input[],
1058                       cmsFloat32Number Output[],
1059                       const cmsInterpParams* p)
1060 {
1061        const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;
1062        cmsFloat32Number rest;
1063        cmsFloat32Number pk;
1064        int k0, K0, K1;
1065        const cmsFloat32Number* T;
1066        cmsUInt32Number i;
1067        cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
1068        cmsInterpParams p1;
1069
1070        pk = Input[0] * p->Domain[0];
1071        k0 = _cmsQuickFloor(pk);
1072        rest = pk - (cmsFloat32Number) k0;
1073
1074        K0 = p -> opta[4] * k0;
1075        K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[4]);
1076
1077        p1 = *p;
1078        memmove(&p1.Domain[0], &p ->Domain[1], 4*sizeof(cmsUInt32Number));
1079
1080        T = LutTable + K0;
1081        p1.Table = T;
1082
1083        Eval4InputsFloat(Input + 1,  Tmp1, &p1);
1084
1085        T = LutTable + K1;
1086        p1.Table = T;
1087
1088        Eval4InputsFloat(Input + 1,  Tmp2, &p1);
1089
1090        for (i=0; i < p -> nOutputs; i++) {
1091
1092               cmsFloat32Number y0 = Tmp1[i];
1093               cmsFloat32Number y1 = Tmp2[i];
1094
1095               Output[i] = y0 + (y1 - y0) * rest;
1096        }
1097 }
1098
1099
1100
1101 static
1102 void Eval6Inputs(register const cmsUInt16Number Input[],
1103                  register cmsUInt16Number Output[],
1104                  register const cmsInterpParams* p16)
1105 {
1106        const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;
1107        cmsS15Fixed16Number fk;
1108        cmsS15Fixed16Number k0, rk;
1109        int K0, K1;
1110        const cmsUInt16Number* T;
1111        cmsUInt32Number i;
1112        cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
1113        cmsInterpParams p1;
1114
1115        fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);
1116        k0 = FIXED_TO_INT(fk);
1117        rk = FIXED_REST_TO_INT(fk);
1118
1119        K0 = p16 -> opta[5] * k0;
1120        K1 = p16 -> opta[5] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));
1121
1122        p1 = *p16;
1123        memmove(&p1.Domain[0], &p16 ->Domain[1], 5*sizeof(cmsUInt32Number));
1124
1125        T = LutTable + K0;
1126        p1.Table = T;
1127
1128        Eval5Inputs(Input + 1, Tmp1, &p1);
1129
1130        T = LutTable + K1;
1131        p1.Table = T;
1132
1133        Eval5Inputs(Input + 1, Tmp2, &p1);
1134
1135        for (i=0; i < p16 -> nOutputs; i++) {
1136
1137               Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);
1138        }
1139
1140 }
1141
1142
1143 static
1144 void Eval6InputsFloat(const cmsFloat32Number Input[],
1145                       cmsFloat32Number Output[],
1146                       const cmsInterpParams* p)
1147 {
1148        const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;
1149        cmsFloat32Number rest;
1150        cmsFloat32Number pk;
1151        int k0, K0, K1;
1152        const cmsFloat32Number* T;
1153        cmsUInt32Number i;
1154        cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
1155        cmsInterpParams p1;
1156
1157        pk = Input[0] * p->Domain[0];
1158        k0 = _cmsQuickFloor(pk);
1159        rest = pk - (cmsFloat32Number) k0;
1160
1161        K0 = p -> opta[5] * k0;
1162        K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[5]);
1163
1164        p1 = *p;
1165        memmove(&p1.Domain[0], &p ->Domain[1], 5*sizeof(cmsUInt32Number));
1166
1167        T = LutTable + K0;
1168        p1.Table = T;
1169
1170        Eval5InputsFloat(Input + 1,  Tmp1, &p1);
1171
1172        T = LutTable + K1;
1173        p1.Table = T;
1174
1175        Eval5InputsFloat(Input + 1,  Tmp2, &p1);
1176
1177        for (i=0; i < p -> nOutputs; i++) {
1178
1179               cmsFloat32Number y0 = Tmp1[i];
1180               cmsFloat32Number y1 = Tmp2[i];
1181
1182               Output[i] = y0 + (y1 - y0) * rest;
1183        }
1184 }
1185
1186
1187 static
1188 void Eval7Inputs(register const cmsUInt16Number Input[],
1189                  register cmsUInt16Number Output[],
1190                  register const cmsInterpParams* p16)
1191 {
1192        const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;
1193        cmsS15Fixed16Number fk;
1194        cmsS15Fixed16Number k0, rk;
1195        int K0, K1;
1196        const cmsUInt16Number* T;
1197        cmsUInt32Number i;
1198        cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
1199        cmsInterpParams p1;
1200
1201
1202        fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);
1203        k0 = FIXED_TO_INT(fk);
1204        rk = FIXED_REST_TO_INT(fk);
1205
1206        K0 = p16 -> opta[6] * k0;
1207        K1 = p16 -> opta[6] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));
1208
1209        p1 = *p16;
1210        memmove(&p1.Domain[0], &p16 ->Domain[1], 6*sizeof(cmsUInt32Number));
1211
1212        T = LutTable + K0;
1213        p1.Table = T;
1214
1215        Eval6Inputs(Input + 1, Tmp1, &p1);
1216
1217        T = LutTable + K1;
1218        p1.Table = T;
1219
1220        Eval6Inputs(Input + 1, Tmp2, &p1);
1221
1222        for (i=0; i < p16 -> nOutputs; i++) {
1223               Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);
1224        }
1225 }
1226
1227
1228 static
1229 void Eval7InputsFloat(const cmsFloat32Number Input[],
1230                       cmsFloat32Number Output[],
1231                       const cmsInterpParams* p)
1232 {
1233        const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;
1234        cmsFloat32Number rest;
1235        cmsFloat32Number pk;
1236        int k0, K0, K1;
1237        const cmsFloat32Number* T;
1238        cmsUInt32Number i;
1239        cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
1240        cmsInterpParams p1;
1241
1242        pk = Input[0] * p->Domain[0];
1243        k0 = _cmsQuickFloor(pk);
1244        rest = pk - (cmsFloat32Number) k0;
1245
1246        K0 = p -> opta[6] * k0;
1247        K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[6]);
1248
1249        p1 = *p;
1250        memmove(&p1.Domain[0], &p ->Domain[1], 6*sizeof(cmsUInt32Number));
1251
1252        T = LutTable + K0;
1253        p1.Table = T;
1254
1255        Eval6InputsFloat(Input + 1,  Tmp1, &p1);
1256
1257        T = LutTable + K1;
1258        p1.Table = T;
1259
1260        Eval6InputsFloat(Input + 1,  Tmp2, &p1);
1261
1262
1263        for (i=0; i < p -> nOutputs; i++) {
1264
1265               cmsFloat32Number y0 = Tmp1[i];
1266               cmsFloat32Number y1 = Tmp2[i];
1267
1268               Output[i] = y0 + (y1 - y0) * rest;
1269
1270        }
1271 }
1272
1273 static
1274 void Eval8Inputs(register const cmsUInt16Number Input[],
1275                  register cmsUInt16Number Output[],
1276                  register const cmsInterpParams* p16)
1277 {
1278        const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;
1279        cmsS15Fixed16Number fk;
1280        cmsS15Fixed16Number k0, rk;
1281        int K0, K1;
1282        const cmsUInt16Number* T;
1283        cmsUInt32Number i;
1284        cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
1285        cmsInterpParams p1;
1286
1287        fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);
1288        k0 = FIXED_TO_INT(fk);
1289        rk = FIXED_REST_TO_INT(fk);
1290
1291        K0 = p16 -> opta[7] * k0;
1292        K1 = p16 -> opta[7] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));
1293
1294        p1 = *p16;
1295        memmove(&p1.Domain[0], &p16 ->Domain[1], 7*sizeof(cmsUInt32Number));
1296
1297        T = LutTable + K0;
1298        p1.Table = T;
1299
1300        Eval7Inputs(Input + 1, Tmp1, &p1);
1301
1302        T = LutTable + K1;
1303        p1.Table = T;
1304        Eval7Inputs(Input + 1, Tmp2, &p1);
1305
1306        for (i=0; i < p16 -> nOutputs; i++) {
1307               Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);
1308        }
1309 }
1310
1311
1312
1313 static
1314 void Eval8InputsFloat(const cmsFloat32Number Input[],
1315                       cmsFloat32Number Output[],
1316                       const cmsInterpParams* p)
1317 {
1318        const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;
1319        cmsFloat32Number rest;
1320        cmsFloat32Number pk;
1321        int k0, K0, K1;
1322        const cmsFloat32Number* T;
1323        cmsUInt32Number i;
1324        cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
1325        cmsInterpParams p1;
1326
1327        pk = Input[0] * p->Domain[0];
1328        k0 = _cmsQuickFloor(pk);
1329        rest = pk - (cmsFloat32Number) k0;
1330
1331        K0 = p -> opta[7] * k0;
1332        K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[7]);
1333
1334        p1 = *p;
1335        memmove(&p1.Domain[0], &p ->Domain[1], 7*sizeof(cmsUInt32Number));
1336
1337        T = LutTable + K0;
1338        p1.Table = T;
1339
1340        Eval7InputsFloat(Input + 1,  Tmp1, &p1);
1341
1342        T = LutTable + K1;
1343        p1.Table = T;
1344
1345        Eval7InputsFloat(Input + 1,  Tmp2, &p1);
1346
1347
1348        for (i=0; i < p -> nOutputs; i++) {
1349
1350               cmsFloat32Number y0 = Tmp1[i];
1351               cmsFloat32Number y1 = Tmp2[i];
1352
1353               Output[i] = y0 + (y1 - y0) * rest;
1354        }
1355 }
1356
1357 // The default factory
1358 static
1359 cmsInterpFunction DefaultInterpolatorsFactory(cmsUInt32Number nInputChannels, cmsUInt32Number nOutputChannels, cmsUInt32Number dwFlags)
1360 {
1361
1362     cmsInterpFunction Interpolation;
1363     cmsBool  IsFloat     = (dwFlags & CMS_LERP_FLAGS_FLOAT);
1364     cmsBool  IsTrilinear = (dwFlags & CMS_LERP_FLAGS_TRILINEAR);
1365
1366     memset(&Interpolation, 0, sizeof(Interpolation));
1367
1368     // Safety check
1369     if (nInputChannels >= 4 && nOutputChannels >= MAX_STAGE_CHANNELS)
1370         return Interpolation;
1371
1372     switch (nInputChannels) {
1373
1374            case 1: // Gray LUT / linear
1375
1376                if (nOutputChannels == 1) {
1377
1378                    if (IsFloat)
1379                        Interpolation.LerpFloat = LinLerp1Dfloat;
1380                    else
1381                        Interpolation.Lerp16 = LinLerp1D;
1382
1383                }
1384                else {
1385
1386                    if (IsFloat)
1387                        Interpolation.LerpFloat = Eval1InputFloat;
1388                    else
1389                        Interpolation.Lerp16 = Eval1Input;
1390                }
1391                break;
1392
1393            case 2: // Duotone
1394                if (IsFloat)
1395                       Interpolation.LerpFloat =  BilinearInterpFloat;
1396                else
1397                       Interpolation.Lerp16    =  BilinearInterp16;
1398                break;
1399
1400            case 3:  // RGB et al
1401
1402                if (IsTrilinear) {
1403
1404                    if (IsFloat)
1405                        Interpolation.LerpFloat = TrilinearInterpFloat;
1406                    else
1407                        Interpolation.Lerp16 = TrilinearInterp16;
1408                }
1409                else {
1410
1411                    if (IsFloat)
1412                        Interpolation.LerpFloat = TetrahedralInterpFloat;
1413                    else {
1414
1415                        Interpolation.Lerp16 = TetrahedralInterp16;
1416                    }
1417                }
1418                break;
1419
1420            case 4:  // CMYK lut
1421
1422                if (IsFloat)
1423                    Interpolation.LerpFloat =  Eval4InputsFloat;
1424                else
1425                    Interpolation.Lerp16    =  Eval4Inputs;
1426                break;
1427
1428            case 5: // 5 Inks
1429                if (IsFloat)
1430                    Interpolation.LerpFloat =  Eval5InputsFloat;
1431                else
1432                    Interpolation.Lerp16    =  Eval5Inputs;
1433                break;
1434
1435            case 6: // 6 Inks
1436                if (IsFloat)
1437                    Interpolation.LerpFloat =  Eval6InputsFloat;
1438                else
1439                    Interpolation.Lerp16    =  Eval6Inputs;
1440                break;
1441
1442            case 7: // 7 inks
1443                if (IsFloat)
1444                    Interpolation.LerpFloat =  Eval7InputsFloat;
1445                else
1446                    Interpolation.Lerp16    =  Eval7Inputs;
1447                break;
1448
1449            case 8: // 8 inks
1450                if (IsFloat)
1451                    Interpolation.LerpFloat =  Eval8InputsFloat;
1452                else
1453                    Interpolation.Lerp16    =  Eval8Inputs;
1454                break;
1455
1456                break;
1457
1458            default:
1459                Interpolation.Lerp16 = NULL;
1460     }
1461
1462     return Interpolation;
1463 }