thirdparty/liblcms2/src/cmsintrp.c

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