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