Imported Upstream version 2.4
[platform/upstream/lcms2.git] / src / cmsintrp.c
index 5b9a09a..7cd7ce9 100644 (file)
-//---------------------------------------------------------------------------------
-//
-//  Little Color Management System
-//  Copyright (c) 1998-2010 Marti Maria Saguer
-//
-// Permission is hereby granted, free of charge, to any person obtaining 
-// a copy of this software and associated documentation files (the "Software"), 
-// to deal in the Software without restriction, including without limitation 
-// the rights to use, copy, modify, merge, publish, distribute, sublicense, 
-// and/or sell copies of the Software, and to permit persons to whom the Software 
-// is furnished to do so, subject to the following conditions:
-//
-// The above copyright notice and this permission notice shall be included in 
-// all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
-// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO 
-// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 
-// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 
-// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION 
-// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 
-// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-//
-//---------------------------------------------------------------------------------
-//
-
-#include "lcms2_internal.h"
-
-// This module incorporates several interpolation routines, for 1 to 8 channels on input and
-// up to 65535 channels on output. The user may change those by using the interpolation plug-in
-
-// Interpolation routines by default
-static cmsInterpFunction DefaultInterpolatorsFactory(cmsUInt32Number nInputChannels, cmsUInt32Number nOutputChannels, cmsUInt32Number dwFlags);
-
-// This is the default factory
-static cmsInterpFnFactory Interpolators = DefaultInterpolatorsFactory;
-
-
-// Main plug-in entry
-cmsBool  _cmsRegisterInterpPlugin(cmsPluginBase* Data)
-{
-    cmsPluginInterpolation* Plugin = (cmsPluginInterpolation*) Data;
-
-    if (Data == NULL) {
-    
-        Interpolators = DefaultInterpolatorsFactory;
-        return TRUE;
-    }
-
-    // Set replacement functions
-    Interpolators = Plugin ->InterpolatorsFactory;  
-    return TRUE;
-}
-
-
-// Set the interpolation method
-
-cmsBool _cmsSetInterpolationRoutine(cmsInterpParams* p)
-{
-    // Invoke factory, possibly in the Plug-in
-    p ->Interpolation = Interpolators(p -> nInputs, p ->nOutputs, p ->dwFlags);
-
-    // If unsupported by the plug-in, go for the LittleCMS default. 
-    // If happens only if an extern plug-in is being used
-    if (p ->Interpolation.Lerp16 == NULL)
-        p ->Interpolation = DefaultInterpolatorsFactory(p ->nInputs, p ->nOutputs, p ->dwFlags);
-
-    // Check for valid interpolator (we just check one member of the union)
-    if (p ->Interpolation.Lerp16 == NULL) {    
-            return FALSE;
-    }
-    return TRUE;
-}
-
-
-// This function precalculates as many parameters as possible to speed up the interpolation.
-cmsInterpParams* _cmsComputeInterpParamsEx(cmsContext ContextID,
-                                           const cmsUInt32Number nSamples[], 
-                                           int InputChan, int OutputChan, 
-                                           const void *Table,
-                                           cmsUInt32Number dwFlags)
-{        
-    cmsInterpParams* p;
-    int i;
-    
-    // Check for maximum inputs
-    if (InputChan > MAX_INPUT_DIMENSIONS) {
-             cmsSignalError(ContextID, cmsERROR_RANGE, "Too many input channels (%d channels, max=%d)", InputChan, MAX_INPUT_DIMENSIONS);            
-            return NULL;
-    }
-
-    // Creates an empty object
-    p = (cmsInterpParams*) _cmsMallocZero(ContextID, sizeof(cmsInterpParams));
-    if (p == NULL) return NULL;
-
-    // Keep original parameters
-    p -> dwFlags  = dwFlags;
-    p -> nInputs  = InputChan;
-    p -> nOutputs = OutputChan;
-    p ->Table     = Table;
-    p ->ContextID  = ContextID;
-
-    // Fill samples per input direction and domain (which is number of nodes minus one)
-    for (i=0; i < InputChan; i++) {
-
-        p -> nSamples[i] = nSamples[i];
-        p -> Domain[i]   = nSamples[i] - 1;     
-    }
-
-    // Compute factors to apply to each component to index the grid array
-    p -> opta[0] = p -> nOutputs;  
-    for (i=1; i < InputChan; i++)
-        p ->opta[i] = p ->opta[i-1] * nSamples[InputChan-i];
-
-
-    if (!_cmsSetInterpolationRoutine(p)) {
-         cmsSignalError(ContextID, cmsERROR_UNKNOWN_EXTENSION, "Unsupported interpolation (%d->%d channels)", InputChan, OutputChan);            
-        _cmsFree(ContextID, p);
-        return NULL;
-    }
-
-    // All seems ok
-    return p;
-}
-
-
-// This one is a wrapper on the anterior, but assuming all directions have same number of nodes
-cmsInterpParams* _cmsComputeInterpParams(cmsContext ContextID, int nSamples, int InputChan, int OutputChan, const void* Table, cmsUInt32Number dwFlags)
-{
-    int i;
-    cmsUInt32Number Samples[MAX_INPUT_DIMENSIONS];
-
-    // Fill the auxiliar array
-    for (i=0; i < MAX_INPUT_DIMENSIONS; i++) 
-        Samples[i] = nSamples;
-
-    // Call the extended function
-    return _cmsComputeInterpParamsEx(ContextID, Samples, InputChan, OutputChan, Table, dwFlags);
-}
-
-
-// Free all associated memory
-void _cmsFreeInterpParams(cmsInterpParams* p)
-{
-    if (p != NULL) _cmsFree(p ->ContextID, p);
-}
-
-
-// Inline fixed point interpolation 
-cmsINLINE cmsUInt16Number LinearInterp(cmsS15Fixed16Number a, cmsS15Fixed16Number l, cmsS15Fixed16Number h)
-{
-    cmsUInt32Number dif = (cmsUInt32Number) (h - l) * a + 0x8000;           
-    dif = (dif >> 16) + l;        
-    return (cmsUInt16Number) (dif);
-}
-
-
-//  Linear interpolation (Fixed-point optimized)
-static
-void LinLerp1D(register const cmsUInt16Number Value[], 
-               register cmsUInt16Number Output[], 
-               register const cmsInterpParams* p)
-{
-    cmsUInt16Number y1, y0; 
-    int cell0, rest;
-    int val3;
-    const cmsUInt16Number* LutTable = (cmsUInt16Number*) p ->Table;
-
-    // if last value... 
-    if (Value[0] == 0xffff) {
-
-        Output[0] = LutTable[p -> Domain[0]];
-        return;
-    }
-
-    val3 = p -> Domain[0] * Value[0];
-    val3 = _cmsToFixedDomain(val3);    // To fixed 15.16
-
-    cell0 = FIXED_TO_INT(val3);             // Cell is 16 MSB bits
-    rest  = FIXED_REST_TO_INT(val3);        // Rest is 16 LSB bits
-
-    y0 = LutTable[cell0];
-    y1 = LutTable[cell0+1];
-
-    
-    Output[0] = LinearInterp(rest, y0, y1);
-}
-
-
-// Floating-point version of 1D interpolation
-static
-void LinLerp1Dfloat(const cmsFloat32Number Value[], 
-                    cmsFloat32Number Output[],  
-                    const cmsInterpParams* p)
-{
-       cmsFloat32Number y1, y0;
-       cmsFloat32Number val2, rest;
-       int cell0, cell1;
-       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table;
-
-       // if last value...
-       if (Value[0] == 1.0) {
-           Output[0] = LutTable[p -> Domain[0]];
-           return;
-       }
-
-       val2 = p -> Domain[0] * Value[0];
-
-       cell0 = (int) floor(val2);
-       cell1 = (int) ceil(val2);
-
-       // Rest is 16 LSB bits
-       rest = val2 - cell0;
-
-       y0 = LutTable[cell0] ;
-       y1 = LutTable[cell1] ;
-
-       Output[0] = y0 + (y1 - y0) * rest;     
-}
-
-
-
-// Eval gray LUT having only one input channel 
-static
-void Eval1Input(register const cmsUInt16Number Input[], 
-                register cmsUInt16Number Output[], 
-                register const cmsInterpParams* p16)
-{
-       cmsS15Fixed16Number fk;
-       cmsS15Fixed16Number k0, k1, rk, K0, K1;
-       int v;
-       cmsUInt32Number OutChan;
-       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;
-
-       v = Input[0] * p16 -> Domain[0];
-       fk = _cmsToFixedDomain(v);
-
-       k0 = FIXED_TO_INT(fk);
-       rk = (cmsUInt16Number) FIXED_REST_TO_INT(fk);
-
-       k1 = k0 + (Input[0] != 0xFFFFU ? 1 : 0);
-
-       K0 = p16 -> opta[0] * k0;
-       K1 = p16 -> opta[0] * k1;
-
-       for (OutChan=0; OutChan < p16->nOutputs; OutChan++) {
-
-           Output[OutChan] = LinearInterp(rk, LutTable[K0+OutChan], LutTable[K1+OutChan]);
-       }
-}
-
-
-
-// Eval gray LUT having only one input channel 
-static
-void Eval1InputFloat(const cmsFloat32Number Value[], 
-                     cmsFloat32Number Output[], 
-                     const cmsInterpParams* p)
-{
-    cmsFloat32Number y1, y0;
-    cmsFloat32Number val2, rest;
-    int cell0, cell1;
-    cmsUInt32Number OutChan;
-    const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table; 
-
-        // if last value...
-       if (Value[0] == 1.0) {
-           Output[0] = LutTable[p -> Domain[0]];
-           return;
-       }
-
-       val2 = p -> Domain[0] * Value[0];
-
-       cell0 = (int) floor(val2);
-       cell1 = (int) ceil(val2);
-
-       // Rest is 16 LSB bits
-       rest = val2 - cell0;
-       
-       cell0 *= p -> opta[0];
-       cell1 *= p -> opta[0];
-
-       for (OutChan=0; OutChan < p->nOutputs; OutChan++) {
-
-            y0 = LutTable[cell0 + OutChan] ;
-            y1 = LutTable[cell1 + OutChan] ;
-
-            Output[OutChan] = y0 + (y1 - y0) * rest;              
-       }
-}
-
-// Bilinear interpolation (16 bits) - cmsFloat32Number version
-static
-void BilinearInterpFloat(const cmsFloat32Number Input[], 
-                         cmsFloat32Number Output[],                          
-                         const cmsInterpParams* p)
-
-{
-#   define LERP(a,l,h)    (cmsFloat32Number) ((l)+(((h)-(l))*(a)))
-#   define DENS(i,j)      (LutTable[(i)+(j)+OutChan])
-
-    const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table; 
-    cmsFloat32Number      px, py;
-    int        x0, y0,
-               X0, Y0, X1, Y1;
-    int        TotalOut, OutChan;
-    cmsFloat32Number      fx, fy,
-        d00, d01, d10, d11,
-        dx0, dx1,
-        dxy;
-    
-    TotalOut   = p -> nOutputs;
-    px = Input[0] * p->Domain[0];
-    py = Input[1] * p->Domain[1];
-
-    x0 = (int) _cmsQuickFloor(px); fx = px - (cmsFloat32Number) x0;
-    y0 = (int) _cmsQuickFloor(py); fy = py - (cmsFloat32Number) y0;
-    
-    X0 = p -> opta[1] * x0;
-    X1 = X0 + (Input[0] >= 1.0 ? 0 : p->opta[1]);
-
-    Y0 = p -> opta[0] * y0;
-    Y1 = Y0 + (Input[1] >= 1.0 ? 0 : p->opta[0]);
-   
-    for (OutChan = 0; OutChan < TotalOut; OutChan++) {
-       
-        d00 = DENS(X0, Y0);
-        d01 = DENS(X0, Y1);
-        d10 = DENS(X1, Y0);
-        d11 = DENS(X1, Y1);
-
-        dx0 = LERP(fx, d00, d10);
-        dx1 = LERP(fx, d01, d11);
-
-        dxy = LERP(fy, dx0, dx1);
-
-        Output[OutChan] = dxy;
-    }
-
-
-#   undef LERP
-#   undef DENS
-}
-
-// Bilinear interpolation (16 bits) - optimized version
-static
-void BilinearInterp16(register const cmsUInt16Number Input[], 
-                      register cmsUInt16Number Output[],
-                      register const cmsInterpParams* p)
-
-{
-#define DENS(i,j) (LutTable[(i)+(j)+OutChan])
-#define LERP(a,l,h)     (cmsUInt16Number) (l + ROUND_FIXED_TO_INT(((h-l)*a)))
-
-           const cmsUInt16Number* LutTable = (cmsUInt16Number*) p ->Table; 
-           int        OutChan, TotalOut;
-           cmsS15Fixed16Number    fx, fy;
-  register int        rx, ry;
-           int        x0, y0;
-  register int        X0, X1, Y0, Y1;
-           int        d00, d01, d10, d11,
-                      dx0, dx1,
-                      dxy;
-
-    TotalOut   = p -> nOutputs;
-
-    fx = _cmsToFixedDomain((int) Input[0] * p -> Domain[0]);
-    x0  = FIXED_TO_INT(fx);
-    rx  = FIXED_REST_TO_INT(fx);    // Rest in 0..1.0 domain
-
-
-    fy = _cmsToFixedDomain((int) Input[1] * p -> Domain[1]);
-    y0  = FIXED_TO_INT(fy);
-    ry  = FIXED_REST_TO_INT(fy);
-
-
-    X0 = p -> opta[1] * x0;
-    X1 = X0 + (Input[0] == 0xFFFFU ? 0 : p->opta[1]);
-
-    Y0 = p -> opta[0] * y0;
-    Y1 = Y0 + (Input[1] == 0xFFFFU ? 0 : p->opta[0]);
-   
-    for (OutChan = 0; OutChan < TotalOut; OutChan++) {
-
-        d00 = DENS(X0, Y0);
-        d01 = DENS(X0, Y1);
-        d10 = DENS(X1, Y0);
-        d11 = DENS(X1, Y1);
-
-        dx0 = LERP(rx, d00, d10);
-        dx1 = LERP(rx, d01, d11);
-
-        dxy = LERP(ry, dx0, dx1);
-
-        Output[OutChan] = (cmsUInt16Number) dxy;
-    }
-
-
-#   undef LERP
-#   undef DENS
-}
-
-
-// Trilinear interpolation (16 bits) - cmsFloat32Number version
-static
-void TrilinearInterpFloat(const cmsFloat32Number Input[], 
-                          cmsFloat32Number Output[],                          
-                          const cmsInterpParams* p)
-
-{
-#   define LERP(a,l,h)      (cmsFloat32Number) ((l)+(((h)-(l))*(a)))
-#   define DENS(i,j,k)      (LutTable[(i)+(j)+(k)+OutChan])
-
-    const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table; 
-    cmsFloat32Number      px, py, pz;
-    int        x0, y0, z0,
-               X0, Y0, Z0, X1, Y1, Z1;
-    int        TotalOut, OutChan;
-    cmsFloat32Number      fx, fy, fz,
-        d000, d001, d010, d011,
-        d100, d101, d110, d111,
-        dx00, dx01, dx10, dx11,
-        dxy0, dxy1, dxyz;
-    
-    TotalOut   = p -> nOutputs;
-    
-    // We need some clipping here
-    px = Input[0];
-    py = Input[1];
-    pz = Input[2];
-
-    if (px < 0) px = 0; 
-    if (px > 1) px = 1;
-    if (py < 0) py = 0; 
-    if (py > 1) py = 1;
-    if (pz < 0) pz = 0; 
-    if (pz > 1) pz = 1;
-
-    px *= p->Domain[0];
-    py *= p->Domain[1];
-    pz *= p->Domain[2];
-
-    x0 = (int) _cmsQuickFloor(px); fx = px - (cmsFloat32Number) x0;
-    y0 = (int) _cmsQuickFloor(py); fy = py - (cmsFloat32Number) y0;
-    z0 = (int) _cmsQuickFloor(pz); fz = pz - (cmsFloat32Number) z0;
-    
-    X0 = p -> opta[2] * x0;
-    X1 = X0 + (Input[0] >= 1.0 ? 0 : p->opta[2]);
-
-    Y0 = p -> opta[1] * y0;
-    Y1 = Y0 + (Input[1] >= 1.0 ? 0 : p->opta[1]);
-   
-    Z0 = p -> opta[0] * z0;
-    Z1 = Z0 + (Input[2] >= 1.0 ? 0 : p->opta[0]);
-   
-    for (OutChan = 0; OutChan < TotalOut; OutChan++) {
-
-        d000 = DENS(X0, Y0, Z0);
-        d001 = DENS(X0, Y0, Z1);
-        d010 = DENS(X0, Y1, Z0);
-        d011 = DENS(X0, Y1, Z1);
-
-        d100 = DENS(X1, Y0, Z0);
-        d101 = DENS(X1, Y0, Z1);
-        d110 = DENS(X1, Y1, Z0);
-        d111 = DENS(X1, Y1, Z1);
-
-
-        dx00 = LERP(fx, d000, d100);
-        dx01 = LERP(fx, d001, d101);
-        dx10 = LERP(fx, d010, d110);
-        dx11 = LERP(fx, d011, d111);
-
-        dxy0 = LERP(fy, dx00, dx10);
-        dxy1 = LERP(fy, dx01, dx11);
-
-        dxyz = LERP(fz, dxy0, dxy1);
-
-        Output[OutChan] = dxyz;
-    }
-
-
-#   undef LERP
-#   undef DENS
-}
-
-// Trilinear interpolation (16 bits) - optimized version
-static
-void TrilinearInterp16(register const cmsUInt16Number Input[], 
-                       register cmsUInt16Number Output[],
-                       register const cmsInterpParams* p)
-
-{
-#define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
-#define LERP(a,l,h)     (cmsUInt16Number) (l + ROUND_FIXED_TO_INT(((h-l)*a)))
-
-           const cmsUInt16Number* LutTable = (cmsUInt16Number*) p ->Table; 
-           int        OutChan, TotalOut;
-           cmsS15Fixed16Number    fx, fy, fz;
-  register int        rx, ry, rz;
-           int        x0, y0, z0;
-  register int        X0, X1, Y0, Y1, Z0, Z1;
-           int        d000, d001, d010, d011,
-                      d100, d101, d110, d111,
-                      dx00, dx01, dx10, dx11,
-                      dxy0, dxy1, dxyz;
-
-    TotalOut   = p -> nOutputs;
-
-    fx = _cmsToFixedDomain((int) Input[0] * p -> Domain[0]);
-    x0  = FIXED_TO_INT(fx);
-    rx  = FIXED_REST_TO_INT(fx);    // Rest in 0..1.0 domain
-
-
-    fy = _cmsToFixedDomain((int) Input[1] * p -> Domain[1]);
-    y0  = FIXED_TO_INT(fy);
-    ry  = FIXED_REST_TO_INT(fy);
-
-    fz = _cmsToFixedDomain((int) Input[2] * p -> Domain[2]);
-    z0 = FIXED_TO_INT(fz);
-    rz = FIXED_REST_TO_INT(fz);
-
-
-    X0 = p -> opta[2] * x0;
-    X1 = X0 + (Input[0] == 0xFFFFU ? 0 : p->opta[2]);
-
-    Y0 = p -> opta[1] * y0;
-    Y1 = Y0 + (Input[1] == 0xFFFFU ? 0 : p->opta[1]);
-   
-    Z0 = p -> opta[0] * z0;
-    Z1 = Z0 + (Input[2] == 0xFFFFU ? 0 : p->opta[0]);
-
-    for (OutChan = 0; OutChan < TotalOut; OutChan++) {
-
-        d000 = DENS(X0, Y0, Z0);
-        d001 = DENS(X0, Y0, Z1);
-        d010 = DENS(X0, Y1, Z0);
-        d011 = DENS(X0, Y1, Z1);
-
-        d100 = DENS(X1, Y0, Z0);
-        d101 = DENS(X1, Y0, Z1);
-        d110 = DENS(X1, Y1, Z0);
-        d111 = DENS(X1, Y1, Z1);
-
-
-        dx00 = LERP(rx, d000, d100);
-        dx01 = LERP(rx, d001, d101);
-        dx10 = LERP(rx, d010, d110);
-        dx11 = LERP(rx, d011, d111);
-
-        dxy0 = LERP(ry, dx00, dx10);
-        dxy1 = LERP(ry, dx01, dx11);
-
-        dxyz = LERP(rz, dxy0, dxy1);
-
-        Output[OutChan] = (cmsUInt16Number) dxyz;
-    }
-
-
-#   undef LERP
-#   undef DENS
-}
-
-
-// Tetrahedral interpolation, using Sakamoto algorithm. 
-#define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
-static
-void TetrahedralInterpFloat(const cmsFloat32Number Input[], 
-                            cmsFloat32Number Output[],                          
-                            const cmsInterpParams* p)
-{
-    const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table; 
-    cmsFloat32Number     px, py, pz;
-    int        x0, y0, z0,
-               X0, Y0, Z0, X1, Y1, Z1;
-    cmsFloat32Number     rx, ry, rz;
-    cmsFloat32Number     c0, c1=0, c2=0, c3=0;
-    int                  OutChan, TotalOut;
-
-    TotalOut   = p -> nOutputs;
-
-    // We need some clipping here
-    px = Input[0];
-    py = Input[1];
-    pz = Input[2];
-
-    if (px < 0) px = 0; 
-    if (px > 1) px = 1;
-    if (py < 0) py = 0; 
-    if (py > 1) py = 1;
-    if (pz < 0) pz = 0; 
-    if (pz > 1) pz = 1;
-
-    px *= p->Domain[0];
-    py *= p->Domain[1];
-    pz *= p->Domain[2];
-
-    x0 = (int) _cmsQuickFloor(px); rx = (px - (cmsFloat32Number) x0);
-    y0 = (int) _cmsQuickFloor(py); ry = (py - (cmsFloat32Number) y0);
-    z0 = (int) _cmsQuickFloor(pz); rz = (pz - (cmsFloat32Number) z0);
-
-
-    X0 = p -> opta[2] * x0;
-    X1 = X0 + (Input[0] >= 1.0 ? 0 : p->opta[2]);
-
-    Y0 = p -> opta[1] * y0;
-    Y1 = Y0 + (Input[1] >= 1.0 ? 0 : p->opta[1]);
-   
-    Z0 = p -> opta[0] * z0;
-    Z1 = Z0 + (Input[2] >= 1.0 ? 0 : p->opta[0]);
-   
-    for (OutChan=0; OutChan < TotalOut; OutChan++) {
-
-       // These are the 6 Tetrahedral
-
-        c0 = DENS(X0, Y0, Z0);
-
-        if (rx >= ry && ry >= rz) {
-
-            c1 = DENS(X1, Y0, Z0) - c0;
-            c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);
-            c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
-
-        }
-        else
-            if (rx >= rz && rz >= ry) {            
-
-                c1 = DENS(X1, Y0, Z0) - c0;
-                c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
-                c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);
-
-            }
-            else
-                if (rz >= rx && rx >= ry) {
-
-                    c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);
-                    c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
-                    c3 = DENS(X0, Y0, Z1) - c0;                            
-
-                }
-                else
-                    if (ry >= rx && rx >= rz) {
-
-                        c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);
-                        c2 = DENS(X0, Y1, Z0) - c0;
-                        c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
-
-                    }
-                    else
-                        if (ry >= rz && rz >= rx) {
-
-                            c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
-                            c2 = DENS(X0, Y1, Z0) - c0;
-                            c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);
-
-                        }
-                        else
-                            if (rz >= ry && ry >= rx) {             
-
-                                c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
-                                c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);
-                                c3 = DENS(X0, Y0, Z1) - c0;
-
-                            }
-                            else  {
-                                c1 = c2 = c3 = 0;                               
-                            }
-
-       Output[OutChan] = c0 + c1 * rx + c2 * ry + c3 * rz;
-       }
-
-}
-
-#undef DENS
-
-
-
-#define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
-
-static
-void TetrahedralInterp16(register const cmsUInt16Number Input[],
-                         register cmsUInt16Number Output[],
-                         register const cmsInterpParams* p)
-{
-    const cmsUInt16Number* LutTable = (cmsUInt16Number*) p -> Table;
-    cmsS15Fixed16Number    fx, fy, fz;
-    cmsS15Fixed16Number    rx, ry, rz;
-    int                    x0, y0, z0;
-    cmsS15Fixed16Number    c0, c1, c2, c3, Rest;       
-    cmsUInt32Number        OutChan;
-    cmsS15Fixed16Number    X0, X1, Y0, Y1, Z0, Z1;
-    cmsUInt32Number        TotalOut = p -> nOutputs;
-    
-
-    fx  = _cmsToFixedDomain((int) Input[0] * p -> Domain[0]);
-    fy  = _cmsToFixedDomain((int) Input[1] * p -> Domain[1]);
-    fz  = _cmsToFixedDomain((int) Input[2] * p -> Domain[2]);
-
-    x0  = FIXED_TO_INT(fx);
-    y0  = FIXED_TO_INT(fy); 
-    z0  = FIXED_TO_INT(fz);
-
-    rx  = FIXED_REST_TO_INT(fx);   
-    ry  = FIXED_REST_TO_INT(fy);      
-    rz  = FIXED_REST_TO_INT(fz);
-
-    X0 = p -> opta[2] * x0;
-    X1 = X0 + (Input[0] == 0xFFFFU ? 0 : p->opta[2]);
-
-    Y0 = p -> opta[1] * y0;
-    Y1 = Y0 + (Input[1] == 0xFFFFU ? 0 : p->opta[1]);
-
-    Z0 = p -> opta[0] * z0;
-    Z1 = Z0 + (Input[2] == 0xFFFFU ? 0 : p->opta[0]);
-
-    // These are the 6 Tetrahedral
-    for (OutChan=0; OutChan < TotalOut; OutChan++) {
-
-        c0 = DENS(X0, Y0, Z0);
-
-        if (rx >= ry && ry >= rz) {
-
-            c1 = DENS(X1, Y0, Z0) - c0;
-            c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);
-            c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
-
-        }
-        else
-            if (rx >= rz && rz >= ry) {            
-
-                c1 = DENS(X1, Y0, Z0) - c0;
-                c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
-                c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);
-
-            }
-            else
-                if (rz >= rx && rx >= ry) {
-
-                    c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);
-                    c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
-                    c3 = DENS(X0, Y0, Z1) - c0;                            
-
-                }
-                else
-                    if (ry >= rx && rx >= rz) {
-
-                        c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);
-                        c2 = DENS(X0, Y1, Z0) - c0;
-                        c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
-
-                    }
-                    else
-                        if (ry >= rz && rz >= rx) {
-
-                            c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
-                            c2 = DENS(X0, Y1, Z0) - c0;
-                            c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);
-
-                        }
-                        else
-                            if (rz >= ry && ry >= rx) {             
-
-                                c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
-                                c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);
-                                c3 = DENS(X0, Y0, Z1) - c0;
-
-                            }
-                            else  {
-                                c1 = c2 = c3 = 0;                               
-                            }
-
-                            Rest = c1 * rx + c2 * ry + c3 * rz;                
-
-                            Output[OutChan] = (cmsUInt16Number) c0 + ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest));
-    }
-
-}
-#undef DENS
-
-
-#define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
-static
-void Eval4Inputs(register const cmsUInt16Number Input[], 
-                     register cmsUInt16Number Output[], 
-                     register const cmsInterpParams* p16)
-{                        
-    const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table; 
-    cmsS15Fixed16Number fk;
-    cmsS15Fixed16Number k0, rk;
-    int K0, K1;
-    cmsS15Fixed16Number    fx, fy, fz;
-    cmsS15Fixed16Number    rx, ry, rz;
-    int                    x0, y0, z0;
-    cmsS15Fixed16Number    X0, X1, Y0, Y1, Z0, Z1;
-    cmsUInt32Number i;
-    cmsS15Fixed16Number    c0, c1, c2, c3, Rest;       
-    cmsUInt32Number        OutChan;
-    cmsUInt16Number        Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-
-
-    fk  = _cmsToFixedDomain((int) Input[0] * p16 -> Domain[0]);
-    fx  = _cmsToFixedDomain((int) Input[1] * p16 -> Domain[1]);
-    fy  = _cmsToFixedDomain((int) Input[2] * p16 -> Domain[2]);
-    fz  = _cmsToFixedDomain((int) Input[3] * p16 -> Domain[3]);
-
-    k0  = FIXED_TO_INT(fk);     
-    x0  = FIXED_TO_INT(fx);
-    y0  = FIXED_TO_INT(fy); 
-    z0  = FIXED_TO_INT(fz);
-
-    rk  = FIXED_REST_TO_INT(fk);
-    rx  = FIXED_REST_TO_INT(fx);   
-    ry  = FIXED_REST_TO_INT(fy);      
-    rz  = FIXED_REST_TO_INT(fz);
-
-    K0 = p16 -> opta[3] * k0;
-    K1 = K0 + (Input[0] == 0xFFFFU ? 0 : p16->opta[3]);
-
-    X0 = p16 -> opta[2] * x0;
-    X1 = X0 + (Input[1] == 0xFFFFU ? 0 : p16->opta[2]);
-
-    Y0 = p16 -> opta[1] * y0;
-    Y1 = Y0 + (Input[2] == 0xFFFFU ? 0 : p16->opta[1]);
-
-    Z0 = p16 -> opta[0] * z0;
-    Z1 = Z0 + (Input[3] == 0xFFFFU ? 0 : p16->opta[0]);
-
-    LutTable = (cmsUInt16Number*) p16 -> Table;
-    LutTable += K0;
-
-    for (OutChan=0; OutChan < p16 -> nOutputs; OutChan++) {
-
-        c0 = DENS(X0, Y0, Z0);
-
-        if (rx >= ry && ry >= rz) {
-
-            c1 = DENS(X1, Y0, Z0) - c0;
-            c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);
-            c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
-
-        }
-        else
-            if (rx >= rz && rz >= ry) {            
-
-                c1 = DENS(X1, Y0, Z0) - c0;
-                c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
-                c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);
-
-            }
-            else
-                if (rz >= rx && rx >= ry) {
-
-                    c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);
-                    c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
-                    c3 = DENS(X0, Y0, Z1) - c0;                            
-
-                }
-                else
-                    if (ry >= rx && rx >= rz) {
-
-                        c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);
-                        c2 = DENS(X0, Y1, Z0) - c0;
-                        c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
-
-                    }
-                    else
-                        if (ry >= rz && rz >= rx) {
-
-                            c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
-                            c2 = DENS(X0, Y1, Z0) - c0;
-                            c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);
-
-                        }
-                        else
-                            if (rz >= ry && ry >= rx) {             
-
-                                c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
-                                c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);
-                                c3 = DENS(X0, Y0, Z1) - c0;
-
-                            }
-                            else  {
-                                c1 = c2 = c3 = 0;                               
-                            }
-
-                            Rest = c1 * rx + c2 * ry + c3 * rz;                
-
-                            Tmp1[OutChan] = (cmsUInt16Number) c0 + ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest));
-    }
-
-
-    LutTable = (cmsUInt16Number*) p16 -> Table;
-    LutTable += K1;
-
-    for (OutChan=0; OutChan < p16 -> nOutputs; OutChan++) {
-
-        c0 = DENS(X0, Y0, Z0);
-
-        if (rx >= ry && ry >= rz) {
-
-            c1 = DENS(X1, Y0, Z0) - c0;
-            c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);
-            c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
-
-        }
-        else
-            if (rx >= rz && rz >= ry) {            
-
-                c1 = DENS(X1, Y0, Z0) - c0;
-                c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
-                c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);
-
-            }
-            else
-                if (rz >= rx && rx >= ry) {
-
-                    c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);
-                    c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
-                    c3 = DENS(X0, Y0, Z1) - c0;                            
-
-                }
-                else
-                    if (ry >= rx && rx >= rz) {
-
-                        c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);
-                        c2 = DENS(X0, Y1, Z0) - c0;
-                        c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
-
-                    }
-                    else
-                        if (ry >= rz && rz >= rx) {
-
-                            c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
-                            c2 = DENS(X0, Y1, Z0) - c0;
-                            c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);
-
-                        }
-                        else
-                            if (rz >= ry && ry >= rx) {             
-
-                                c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
-                                c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);
-                                c3 = DENS(X0, Y0, Z1) - c0;
-
-                            }
-                            else  {
-                                c1 = c2 = c3 = 0;                               
-                            }
-
-                            Rest = c1 * rx + c2 * ry + c3 * rz;                
-
-                            Tmp2[OutChan] = (cmsUInt16Number) c0 + ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest));
-    }
-
-
-
-    for (i=0; i < p16 -> nOutputs; i++) {
-        Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);              
-    }
-}
-#undef DENS
-
-
-// For more that 3 inputs (i.e., CMYK)
-// evaluate two 3-dimensional interpolations and then linearly interpolate between them.
-
-
-static
-void Eval4InputsFloat(const cmsFloat32Number Input[], 
-                      cmsFloat32Number Output[], 
-                      const cmsInterpParams* p)
-{                             
-       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table; 
-       cmsFloat32Number rest;
-       cmsFloat32Number pk;
-       int k0, K0, K1;
-       const cmsFloat32Number* T;
-       cmsUInt32Number i;
-       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-       cmsInterpParams p1;
-
-       
-       pk = Input[0] * p->Domain[0];    
-       k0 = _cmsQuickFloor(pk); 
-       rest = pk - (cmsFloat32Number) k0;
-
-       K0 = p -> opta[3] * k0;
-       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[3]);
-
-       p1 = *p;
-       memmove(&p1.Domain[0], &p ->Domain[1], 3*sizeof(cmsUInt32Number));
-       
-       T = LutTable + K0;
-       p1.Table = T;
-
-       TetrahedralInterpFloat(Input + 1,  Tmp1, &p1);
-
-       T = LutTable + K1;
-       p1.Table = T;
-       TetrahedralInterpFloat(Input + 1,  Tmp2, &p1);
-
-       for (i=0; i < p -> nOutputs; i++)
-       {
-              cmsFloat32Number y0 = Tmp1[i];
-              cmsFloat32Number y1 = Tmp2[i];
-
-              Output[i] = y0 + (y1 - y0) * rest;               
-       }
-}
-
-
-static
-void Eval5Inputs(register const cmsUInt16Number Input[], 
-                 register cmsUInt16Number Output[], 
-
-                 register const cmsInterpParams* p16)
-{       
-       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;  
-       cmsS15Fixed16Number fk;
-       cmsS15Fixed16Number k0, rk;
-       int K0, K1;
-       const cmsUInt16Number* T;
-       cmsUInt32Number i;
-       cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-       cmsInterpParams p1;
-
-       
-       fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);
-       k0 = FIXED_TO_INT(fk);
-       rk = FIXED_REST_TO_INT(fk);
-
-       K0 = p16 -> opta[4] * k0;
-       K1 = p16 -> opta[4] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));
-
-       p1 = *p16;
-       memmove(&p1.Domain[0], &p16 ->Domain[1], 4*sizeof(cmsUInt32Number));
-
-       T = LutTable + K0;
-       p1.Table = T;
-
-       Eval4Inputs(Input + 1, Tmp1, &p1);
-
-       T = LutTable + K1;
-       p1.Table = T;
-
-       Eval4Inputs(Input + 1, Tmp2, &p1);
-
-       for (i=0; i < p16 -> nOutputs; i++) {
-
-              Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);                         
-       }
-
-}
-
-
-static
-void Eval5InputsFloat(const cmsFloat32Number Input[], 
-                      cmsFloat32Number Output[], 
-                      const cmsInterpParams* p)
-{                             
-       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table; 
-       cmsFloat32Number rest;
-       cmsFloat32Number pk;
-       int k0, K0, K1;
-       const cmsFloat32Number* T;
-       cmsUInt32Number i;
-       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-       cmsInterpParams p1;
-
-       pk = Input[0] * p->Domain[0];    
-       k0 = _cmsQuickFloor(pk); 
-       rest = pk - (cmsFloat32Number) k0;
-
-       K0 = p -> opta[4] * k0;
-       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[4]);
-
-       p1 = *p;
-       memmove(&p1.Domain[0], &p ->Domain[1], 4*sizeof(cmsUInt32Number));
-
-       T = LutTable + K0;
-       p1.Table = T;
-       
-       Eval4InputsFloat(Input + 1,  Tmp1, &p1);
-      
-       T = LutTable + K1;
-       p1.Table = T;
-
-       Eval4InputsFloat(Input + 1,  Tmp2, &p1);
-      
-       for (i=0; i < p -> nOutputs; i++) {
-
-              cmsFloat32Number y0 = Tmp1[i];
-              cmsFloat32Number y1 = Tmp2[i];
-
-              Output[i] = y0 + (y1 - y0) * rest;             
-       }
-}
-
-
-
-static
-void Eval6Inputs(register const cmsUInt16Number Input[], 
-                 register cmsUInt16Number Output[], 
-                 register const cmsInterpParams* p16)
-{       
-       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table; 
-       cmsS15Fixed16Number fk;
-       cmsS15Fixed16Number k0, rk;
-       int K0, K1;
-       const cmsUInt16Number* T;
-       cmsUInt32Number i;
-       cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-       cmsInterpParams p1;
-       
-       fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);
-       k0 = FIXED_TO_INT(fk);
-       rk = FIXED_REST_TO_INT(fk);
-
-       K0 = p16 -> opta[5] * k0;
-       K1 = p16 -> opta[5] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));
-
-       p1 = *p16;
-       memmove(&p1.Domain[0], &p16 ->Domain[1], 5*sizeof(cmsUInt32Number));
-
-       T = LutTable + K0;
-       p1.Table = T;
-
-       Eval5Inputs(Input + 1, Tmp1, &p1);
-
-       T = LutTable + K1;
-       p1.Table = T;
-
-       Eval5Inputs(Input + 1, Tmp2, &p1);
-
-       for (i=0; i < p16 -> nOutputs; i++) {
-
-              Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);
-       }
-
-}
-
-
-static
-void Eval6InputsFloat(const cmsFloat32Number Input[], 
-                      cmsFloat32Number Output[], 
-                      const cmsInterpParams* p)
-{                             
-       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table; 
-       cmsFloat32Number rest;
-       cmsFloat32Number pk;
-       int k0, K0, K1;
-       const cmsFloat32Number* T;
-       cmsUInt32Number i;
-       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-       cmsInterpParams p1;
-
-       pk = Input[0] * p->Domain[0];    
-       k0 = _cmsQuickFloor(pk); 
-       rest = pk - (cmsFloat32Number) k0;
-
-       K0 = p -> opta[5] * k0;
-       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[5]);
-
-       p1 = *p;
-       memmove(&p1.Domain[0], &p ->Domain[1], 5*sizeof(cmsUInt32Number));
-       
-       T = LutTable + K0;
-       p1.Table = T;
-  
-       Eval5InputsFloat(Input + 1,  Tmp1, &p1);
-      
-       T = LutTable + K1;
-       p1.Table = T;
-
-       Eval5InputsFloat(Input + 1,  Tmp2, &p1);
-      
-       for (i=0; i < p -> nOutputs; i++) {
-
-              cmsFloat32Number y0 = Tmp1[i];
-              cmsFloat32Number y1 = Tmp2[i];
-
-              Output[i] = y0 + (y1 - y0) * rest;                   
-       }
-}
-
-
-static
-void Eval7Inputs(register const cmsUInt16Number Input[], 
-                 register cmsUInt16Number Output[], 
-                 register const cmsInterpParams* p16)
-{       
-       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table; 
-       cmsS15Fixed16Number fk;
-       cmsS15Fixed16Number k0, rk;
-       int K0, K1;
-       const cmsUInt16Number* T;
-       cmsUInt32Number i;
-       cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-       cmsInterpParams p1;
-
-       
-       fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);
-       k0 = FIXED_TO_INT(fk);
-       rk = FIXED_REST_TO_INT(fk);
-
-       K0 = p16 -> opta[6] * k0;
-       K1 = p16 -> opta[6] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));
-
-       p1 = *p16;
-       memmove(&p1.Domain[0], &p16 ->Domain[1], 6*sizeof(cmsUInt32Number));
-       
-       T = LutTable + K0;
-       p1.Table = T;
-
-       Eval6Inputs(Input + 1, Tmp1, &p1);
-
-       T = LutTable + K1;
-       p1.Table = T;
-
-       Eval6Inputs(Input + 1, Tmp2, &p1);
-
-       for (i=0; i < p16 -> nOutputs; i++) {
-              Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);
-       }
-}
-
-
-static
-void Eval7InputsFloat(const cmsFloat32Number Input[], 
-                      cmsFloat32Number Output[], 
-                      const cmsInterpParams* p)
-{                             
-       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table; 
-       cmsFloat32Number rest;
-       cmsFloat32Number pk;
-       int k0, K0, K1;
-       const cmsFloat32Number* T;
-       cmsUInt32Number i;
-       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-       cmsInterpParams p1;
-
-       pk = Input[0] * p->Domain[0];    
-       k0 = _cmsQuickFloor(pk); 
-       rest = pk - (cmsFloat32Number) k0;
-
-       K0 = p -> opta[6] * k0;
-       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[6]);
-
-       p1 = *p;
-       memmove(&p1.Domain[0], &p ->Domain[1], 6*sizeof(cmsUInt32Number));
-       
-       T = LutTable + K0;
-       p1.Table = T;
-
-       Eval6InputsFloat(Input + 1,  Tmp1, &p1);
-      
-       T = LutTable + K1;
-       p1.Table = T;
-       
-       Eval6InputsFloat(Input + 1,  Tmp2, &p1);
-      
-       
-       for (i=0; i < p -> nOutputs; i++) {
-
-              cmsFloat32Number y0 = Tmp1[i];
-              cmsFloat32Number y1 = Tmp2[i];
-
-              Output[i] = y0 + (y1 - y0) * rest;     
-              
-       }
-}
-
-static
-void Eval8Inputs(register const cmsUInt16Number Input[], 
-                 register cmsUInt16Number Output[], 
-                 register const cmsInterpParams* p16)
-{       
-       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table; 
-       cmsS15Fixed16Number fk;
-       cmsS15Fixed16Number k0, rk;
-       int K0, K1;
-       const cmsUInt16Number* T;
-       cmsUInt32Number i;
-       cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-       cmsInterpParams p1;
-       
-       fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);
-       k0 = FIXED_TO_INT(fk);
-       rk = FIXED_REST_TO_INT(fk);
-
-       K0 = p16 -> opta[7] * k0;
-       K1 = p16 -> opta[7] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));
-
-       p1 = *p16;
-       memmove(&p1.Domain[0], &p16 ->Domain[1], 7*sizeof(cmsUInt32Number));
-       
-       T = LutTable + K0;
-       p1.Table = T;
-
-       Eval7Inputs(Input + 1, Tmp1, &p1);
-
-       T = LutTable + K1;
-       p1.Table = T;
-       Eval7Inputs(Input + 1, Tmp2, &p1);
-
-       for (i=0; i < p16 -> nOutputs; i++) {
-              Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);
-       }
-}
-
-
-
-static
-void Eval8InputsFloat(const cmsFloat32Number Input[], 
-                      cmsFloat32Number Output[], 
-                      const cmsInterpParams* p)
-{                             
-       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table; 
-       cmsFloat32Number rest;
-       cmsFloat32Number pk;
-       int k0, K0, K1;
-       const cmsFloat32Number* T;
-       cmsUInt32Number i;
-       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];
-       cmsInterpParams p1;
-       
-       pk = Input[0] * p->Domain[0];    
-       k0 = _cmsQuickFloor(pk); 
-       rest = pk - (cmsFloat32Number) k0;
-
-       K0 = p -> opta[7] * k0;
-       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[7]);
-
-       p1 = *p;
-       memmove(&p1.Domain[0], &p ->Domain[1], 7*sizeof(cmsUInt32Number));
-
-       T = LutTable + K0;
-       p1.Table = T;
-
-       Eval7InputsFloat(Input + 1,  Tmp1, &p1);
-      
-       T = LutTable + K1;
-       p1.Table = T;
-
-       Eval7InputsFloat(Input + 1,  Tmp2, &p1);
-      
-
-       for (i=0; i < p -> nOutputs; i++) {
-
-              cmsFloat32Number y0 = Tmp1[i];
-              cmsFloat32Number y1 = Tmp2[i];
-
-              Output[i] = y0 + (y1 - y0) * rest;                   
-       }
-}
-
-// The default factory
-static 
-cmsInterpFunction DefaultInterpolatorsFactory(cmsUInt32Number nInputChannels, cmsUInt32Number nOutputChannels, cmsUInt32Number dwFlags)
-{
-
-    cmsInterpFunction Interpolation;
-    cmsBool  IsFloat     = (dwFlags & CMS_LERP_FLAGS_FLOAT);
-    cmsBool  IsTrilinear = (dwFlags & CMS_LERP_FLAGS_TRILINEAR);
-
-    memset(&Interpolation, 0, sizeof(Interpolation));
-
-    // Safety check
-    if (nInputChannels >= 4 && nOutputChannels >= MAX_STAGE_CHANNELS) 
-        return Interpolation;
-
-    switch (nInputChannels) {
-
-           case 1: // Gray LUT / linear
-
-               if (nOutputChannels == 1) {
-
-                   if (IsFloat) 
-                       Interpolation.LerpFloat = LinLerp1Dfloat;
-                   else
-                       Interpolation.Lerp16 = LinLerp1D;
-
-               }
-               else {
-
-                   if (IsFloat)
-                       Interpolation.LerpFloat = Eval1InputFloat;               
-                   else
-                       Interpolation.Lerp16 = Eval1Input;               
-               }
-               break;
-
-           case 2: // Duotone
-               if (IsFloat) 
-                      Interpolation.LerpFloat =  BilinearInterpFloat;
-               else
-                      Interpolation.Lerp16    =  BilinearInterp16;
-               break;
-
-           case 3:  // RGB et al   
-
-               if (IsTrilinear) {
-
-                   if (IsFloat) 
-                       Interpolation.LerpFloat = TrilinearInterpFloat;
-                   else
-                       Interpolation.Lerp16 = TrilinearInterp16;
-               }
-               else {
-
-                   if (IsFloat) 
-                       Interpolation.LerpFloat = TetrahedralInterpFloat;
-                   else {
-
-                       Interpolation.Lerp16 = TetrahedralInterp16;
-                   }
-               }
-               break;
-
-           case 4:  // CMYK lut             
-
-               if (IsFloat) 
-                   Interpolation.LerpFloat =  Eval4InputsFloat;
-               else
-                   Interpolation.Lerp16    =  Eval4Inputs;
-               break;
-
-           case 5: // 5 Inks
-               if (IsFloat) 
-                   Interpolation.LerpFloat =  Eval5InputsFloat;
-               else
-                   Interpolation.Lerp16    =  Eval5Inputs;
-               break;              
-
-           case 6: // 6 Inks
-               if (IsFloat) 
-                   Interpolation.LerpFloat =  Eval6InputsFloat;
-               else
-                   Interpolation.Lerp16    =  Eval6Inputs;
-               break;
-
-           case 7: // 7 inks
-               if (IsFloat) 
-                   Interpolation.LerpFloat =  Eval7InputsFloat;
-               else
-                   Interpolation.Lerp16    =  Eval7Inputs;
-               break;
-
-           case 8: // 8 inks
-               if (IsFloat) 
-                   Interpolation.LerpFloat =  Eval8InputsFloat;
-               else
-                   Interpolation.Lerp16    =  Eval8Inputs;
-               break;
-
-               break;
-
-           default:
-               Interpolation.Lerp16 = NULL;
-    }
-
-    return Interpolation;
-}
+//---------------------------------------------------------------------------------\r
+//\r
+//  Little Color Management System\r
+//  Copyright (c) 1998-2012 Marti Maria Saguer\r
+//\r
+// Permission is hereby granted, free of charge, to any person obtaining\r
+// a copy of this software and associated documentation files (the "Software"),\r
+// to deal in the Software without restriction, including without limitation\r
+// the rights to use, copy, modify, merge, publish, distribute, sublicense,\r
+// and/or sell copies of the Software, and to permit persons to whom the Software\r
+// is furnished to do so, subject to the following conditions:\r
+//\r
+// The above copyright notice and this permission notice shall be included in\r
+// all copies or substantial portions of the Software.\r
+//\r
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,\r
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO\r
+// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND\r
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE\r
+// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION\r
+// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION\r
+// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\r
+//\r
+//---------------------------------------------------------------------------------\r
+//\r
+\r
+#include "lcms2_internal.h"\r
+\r
+// This module incorporates several interpolation routines, for 1 to 8 channels on input and\r
+// up to 65535 channels on output. The user may change those by using the interpolation plug-in\r
+\r
+// Interpolation routines by default\r
+static cmsInterpFunction DefaultInterpolatorsFactory(cmsUInt32Number nInputChannels, cmsUInt32Number nOutputChannels, cmsUInt32Number dwFlags);\r
+\r
+// This is the default factory\r
+static cmsInterpFnFactory Interpolators = DefaultInterpolatorsFactory;\r
+\r
+\r
+// Main plug-in entry\r
+cmsBool  _cmsRegisterInterpPlugin(cmsPluginBase* Data)\r
+{\r
+    cmsPluginInterpolation* Plugin = (cmsPluginInterpolation*) Data;\r
+\r
+    if (Data == NULL) {\r
+\r
+        Interpolators = DefaultInterpolatorsFactory;\r
+        return TRUE;\r
+    }\r
+\r
+    // Set replacement functions\r
+    Interpolators = Plugin ->InterpolatorsFactory;\r
+    return TRUE;\r
+}\r
+\r
+\r
+// Set the interpolation method\r
+\r
+cmsBool _cmsSetInterpolationRoutine(cmsInterpParams* p)\r
+{\r
+    // Invoke factory, possibly in the Plug-in\r
+    p ->Interpolation = Interpolators(p -> nInputs, p ->nOutputs, p ->dwFlags);\r
+\r
+    // If unsupported by the plug-in, go for the LittleCMS default.\r
+    // If happens only if an extern plug-in is being used\r
+    if (p ->Interpolation.Lerp16 == NULL)\r
+        p ->Interpolation = DefaultInterpolatorsFactory(p ->nInputs, p ->nOutputs, p ->dwFlags);\r
+\r
+    // Check for valid interpolator (we just check one member of the union)\r
+    if (p ->Interpolation.Lerp16 == NULL) {\r
+            return FALSE;\r
+    }\r
+    return TRUE;\r
+}\r
+\r
+\r
+// This function precalculates as many parameters as possible to speed up the interpolation.\r
+cmsInterpParams* _cmsComputeInterpParamsEx(cmsContext ContextID,\r
+                                           const cmsUInt32Number nSamples[],\r
+                                           int InputChan, int OutputChan,\r
+                                           const void *Table,\r
+                                           cmsUInt32Number dwFlags)\r
+{\r
+    cmsInterpParams* p;\r
+    int i;\r
+\r
+    // Check for maximum inputs\r
+    if (InputChan > MAX_INPUT_DIMENSIONS) {\r
+             cmsSignalError(ContextID, cmsERROR_RANGE, "Too many input channels (%d channels, max=%d)", InputChan, MAX_INPUT_DIMENSIONS);\r
+            return NULL;\r
+    }\r
+\r
+    // Creates an empty object\r
+    p = (cmsInterpParams*) _cmsMallocZero(ContextID, sizeof(cmsInterpParams));\r
+    if (p == NULL) return NULL;\r
+\r
+    // Keep original parameters\r
+    p -> dwFlags  = dwFlags;\r
+    p -> nInputs  = InputChan;\r
+    p -> nOutputs = OutputChan;\r
+    p ->Table     = Table;\r
+    p ->ContextID  = ContextID;\r
+\r
+    // Fill samples per input direction and domain (which is number of nodes minus one)\r
+    for (i=0; i < InputChan; i++) {\r
+\r
+        p -> nSamples[i] = nSamples[i];\r
+        p -> Domain[i]   = nSamples[i] - 1;\r
+    }\r
+\r
+    // Compute factors to apply to each component to index the grid array\r
+    p -> opta[0] = p -> nOutputs;\r
+    for (i=1; i < InputChan; i++)\r
+        p ->opta[i] = p ->opta[i-1] * nSamples[InputChan-i];\r
+\r
+\r
+    if (!_cmsSetInterpolationRoutine(p)) {\r
+         cmsSignalError(ContextID, cmsERROR_UNKNOWN_EXTENSION, "Unsupported interpolation (%d->%d channels)", InputChan, OutputChan);\r
+        _cmsFree(ContextID, p);\r
+        return NULL;\r
+    }\r
+\r
+    // All seems ok\r
+    return p;\r
+}\r
+\r
+\r
+// This one is a wrapper on the anterior, but assuming all directions have same number of nodes\r
+cmsInterpParams* _cmsComputeInterpParams(cmsContext ContextID, int nSamples, int InputChan, int OutputChan, const void* Table, cmsUInt32Number dwFlags)\r
+{\r
+    int i;\r
+    cmsUInt32Number Samples[MAX_INPUT_DIMENSIONS];\r
+\r
+    // Fill the auxiliar array\r
+    for (i=0; i < MAX_INPUT_DIMENSIONS; i++)\r
+        Samples[i] = nSamples;\r
+\r
+    // Call the extended function\r
+    return _cmsComputeInterpParamsEx(ContextID, Samples, InputChan, OutputChan, Table, dwFlags);\r
+}\r
+\r
+\r
+// Free all associated memory\r
+void _cmsFreeInterpParams(cmsInterpParams* p)\r
+{\r
+    if (p != NULL) _cmsFree(p ->ContextID, p);\r
+}\r
+\r
+\r
+// Inline fixed point interpolation\r
+cmsINLINE cmsUInt16Number LinearInterp(cmsS15Fixed16Number a, cmsS15Fixed16Number l, cmsS15Fixed16Number h)\r
+{\r
+    cmsUInt32Number dif = (cmsUInt32Number) (h - l) * a + 0x8000;\r
+    dif = (dif >> 16) + l;\r
+    return (cmsUInt16Number) (dif);\r
+}\r
+\r
+\r
+//  Linear interpolation (Fixed-point optimized)\r
+static\r
+void LinLerp1D(register const cmsUInt16Number Value[],\r
+               register cmsUInt16Number Output[],\r
+               register const cmsInterpParams* p)\r
+{\r
+    cmsUInt16Number y1, y0;\r
+    int cell0, rest;\r
+    int val3;\r
+    const cmsUInt16Number* LutTable = (cmsUInt16Number*) p ->Table;\r
+\r
+    // if last value...\r
+    if (Value[0] == 0xffff) {\r
+\r
+        Output[0] = LutTable[p -> Domain[0]];\r
+        return;\r
+    }\r
+\r
+    val3 = p -> Domain[0] * Value[0];\r
+    val3 = _cmsToFixedDomain(val3);    // To fixed 15.16\r
+\r
+    cell0 = FIXED_TO_INT(val3);             // Cell is 16 MSB bits\r
+    rest  = FIXED_REST_TO_INT(val3);        // Rest is 16 LSB bits\r
+\r
+    y0 = LutTable[cell0];\r
+    y1 = LutTable[cell0+1];\r
+\r
+\r
+    Output[0] = LinearInterp(rest, y0, y1);\r
+}\r
+\r
+\r
+// Floating-point version of 1D interpolation\r
+static\r
+void LinLerp1Dfloat(const cmsFloat32Number Value[],\r
+                    cmsFloat32Number Output[],\r
+                    const cmsInterpParams* p)\r
+{\r
+       cmsFloat32Number y1, y0;\r
+       cmsFloat32Number val2, rest;\r
+       int cell0, cell1;\r
+       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table;\r
+\r
+       // if last value...\r
+       if (Value[0] == 1.0) {\r
+           Output[0] = LutTable[p -> Domain[0]];\r
+           return;\r
+       }\r
+\r
+       val2 = p -> Domain[0] * Value[0];\r
+\r
+       cell0 = (int) floor(val2);\r
+       cell1 = (int) ceil(val2);\r
+\r
+       // Rest is 16 LSB bits\r
+       rest = val2 - cell0;\r
+\r
+       y0 = LutTable[cell0] ;\r
+       y1 = LutTable[cell1] ;\r
+\r
+       Output[0] = y0 + (y1 - y0) * rest;\r
+}\r
+\r
+\r
+\r
+// Eval gray LUT having only one input channel\r
+static\r
+void Eval1Input(register const cmsUInt16Number Input[],\r
+                register cmsUInt16Number Output[],\r
+                register const cmsInterpParams* p16)\r
+{\r
+       cmsS15Fixed16Number fk;\r
+       cmsS15Fixed16Number k0, k1, rk, K0, K1;\r
+       int v;\r
+       cmsUInt32Number OutChan;\r
+       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;\r
+\r
+       v = Input[0] * p16 -> Domain[0];\r
+       fk = _cmsToFixedDomain(v);\r
+\r
+       k0 = FIXED_TO_INT(fk);\r
+       rk = (cmsUInt16Number) FIXED_REST_TO_INT(fk);\r
+\r
+       k1 = k0 + (Input[0] != 0xFFFFU ? 1 : 0);\r
+\r
+       K0 = p16 -> opta[0] * k0;\r
+       K1 = p16 -> opta[0] * k1;\r
+\r
+       for (OutChan=0; OutChan < p16->nOutputs; OutChan++) {\r
+\r
+           Output[OutChan] = LinearInterp(rk, LutTable[K0+OutChan], LutTable[K1+OutChan]);\r
+       }\r
+}\r
+\r
+\r
+\r
+// Eval gray LUT having only one input channel\r
+static\r
+void Eval1InputFloat(const cmsFloat32Number Value[],\r
+                     cmsFloat32Number Output[],\r
+                     const cmsInterpParams* p)\r
+{\r
+    cmsFloat32Number y1, y0;\r
+    cmsFloat32Number val2, rest;\r
+    int cell0, cell1;\r
+    cmsUInt32Number OutChan;\r
+    const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table;\r
+\r
+        // if last value...\r
+       if (Value[0] == 1.0) {\r
+           Output[0] = LutTable[p -> Domain[0]];\r
+           return;\r
+       }\r
+\r
+       val2 = p -> Domain[0] * Value[0];\r
+\r
+       cell0 = (int) floor(val2);\r
+       cell1 = (int) ceil(val2);\r
+\r
+       // Rest is 16 LSB bits\r
+       rest = val2 - cell0;\r
+\r
+       cell0 *= p -> opta[0];\r
+       cell1 *= p -> opta[0];\r
+\r
+       for (OutChan=0; OutChan < p->nOutputs; OutChan++) {\r
+\r
+            y0 = LutTable[cell0 + OutChan] ;\r
+            y1 = LutTable[cell1 + OutChan] ;\r
+\r
+            Output[OutChan] = y0 + (y1 - y0) * rest;\r
+       }\r
+}\r
+\r
+// Bilinear interpolation (16 bits) - cmsFloat32Number version\r
+static\r
+void BilinearInterpFloat(const cmsFloat32Number Input[],\r
+                         cmsFloat32Number Output[],\r
+                         const cmsInterpParams* p)\r
+\r
+{\r
+#   define LERP(a,l,h)    (cmsFloat32Number) ((l)+(((h)-(l))*(a)))\r
+#   define DENS(i,j)      (LutTable[(i)+(j)+OutChan])\r
+\r
+    const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table;\r
+    cmsFloat32Number      px, py;\r
+    int        x0, y0,\r
+               X0, Y0, X1, Y1;\r
+    int        TotalOut, OutChan;\r
+    cmsFloat32Number      fx, fy,\r
+        d00, d01, d10, d11,\r
+        dx0, dx1,\r
+        dxy;\r
+\r
+    TotalOut   = p -> nOutputs;\r
+    px = Input[0] * p->Domain[0];\r
+    py = Input[1] * p->Domain[1];\r
+\r
+    x0 = (int) _cmsQuickFloor(px); fx = px - (cmsFloat32Number) x0;\r
+    y0 = (int) _cmsQuickFloor(py); fy = py - (cmsFloat32Number) y0;\r
+\r
+    X0 = p -> opta[1] * x0;\r
+    X1 = X0 + (Input[0] >= 1.0 ? 0 : p->opta[1]);\r
+\r
+    Y0 = p -> opta[0] * y0;\r
+    Y1 = Y0 + (Input[1] >= 1.0 ? 0 : p->opta[0]);\r
+\r
+    for (OutChan = 0; OutChan < TotalOut; OutChan++) {\r
+\r
+        d00 = DENS(X0, Y0);\r
+        d01 = DENS(X0, Y1);\r
+        d10 = DENS(X1, Y0);\r
+        d11 = DENS(X1, Y1);\r
+\r
+        dx0 = LERP(fx, d00, d10);\r
+        dx1 = LERP(fx, d01, d11);\r
+\r
+        dxy = LERP(fy, dx0, dx1);\r
+\r
+        Output[OutChan] = dxy;\r
+    }\r
+\r
+\r
+#   undef LERP\r
+#   undef DENS\r
+}\r
+\r
+// Bilinear interpolation (16 bits) - optimized version\r
+static\r
+void BilinearInterp16(register const cmsUInt16Number Input[],\r
+                      register cmsUInt16Number Output[],\r
+                      register const cmsInterpParams* p)\r
+\r
+{\r
+#define DENS(i,j) (LutTable[(i)+(j)+OutChan])\r
+#define LERP(a,l,h)     (cmsUInt16Number) (l + ROUND_FIXED_TO_INT(((h-l)*a)))\r
+\r
+           const cmsUInt16Number* LutTable = (cmsUInt16Number*) p ->Table;\r
+           int        OutChan, TotalOut;\r
+           cmsS15Fixed16Number    fx, fy;\r
+  register int        rx, ry;\r
+           int        x0, y0;\r
+  register int        X0, X1, Y0, Y1;\r
+           int        d00, d01, d10, d11,\r
+                      dx0, dx1,\r
+                      dxy;\r
+\r
+    TotalOut   = p -> nOutputs;\r
+\r
+    fx = _cmsToFixedDomain((int) Input[0] * p -> Domain[0]);\r
+    x0  = FIXED_TO_INT(fx);\r
+    rx  = FIXED_REST_TO_INT(fx);    // Rest in 0..1.0 domain\r
+\r
+\r
+    fy = _cmsToFixedDomain((int) Input[1] * p -> Domain[1]);\r
+    y0  = FIXED_TO_INT(fy);\r
+    ry  = FIXED_REST_TO_INT(fy);\r
+\r
+\r
+    X0 = p -> opta[1] * x0;\r
+    X1 = X0 + (Input[0] == 0xFFFFU ? 0 : p->opta[1]);\r
+\r
+    Y0 = p -> opta[0] * y0;\r
+    Y1 = Y0 + (Input[1] == 0xFFFFU ? 0 : p->opta[0]);\r
+\r
+    for (OutChan = 0; OutChan < TotalOut; OutChan++) {\r
+\r
+        d00 = DENS(X0, Y0);\r
+        d01 = DENS(X0, Y1);\r
+        d10 = DENS(X1, Y0);\r
+        d11 = DENS(X1, Y1);\r
+\r
+        dx0 = LERP(rx, d00, d10);\r
+        dx1 = LERP(rx, d01, d11);\r
+\r
+        dxy = LERP(ry, dx0, dx1);\r
+\r
+        Output[OutChan] = (cmsUInt16Number) dxy;\r
+    }\r
+\r
+\r
+#   undef LERP\r
+#   undef DENS\r
+}\r
+\r
+\r
+// Trilinear interpolation (16 bits) - cmsFloat32Number version\r
+static\r
+void TrilinearInterpFloat(const cmsFloat32Number Input[],\r
+                          cmsFloat32Number Output[],\r
+                          const cmsInterpParams* p)\r
+\r
+{\r
+#   define LERP(a,l,h)      (cmsFloat32Number) ((l)+(((h)-(l))*(a)))\r
+#   define DENS(i,j,k)      (LutTable[(i)+(j)+(k)+OutChan])\r
+\r
+    const cmsFloat32Number* LutTable = (cmsFloat32Number*) p ->Table;\r
+    cmsFloat32Number      px, py, pz;\r
+    int        x0, y0, z0,\r
+               X0, Y0, Z0, X1, Y1, Z1;\r
+    int        TotalOut, OutChan;\r
+    cmsFloat32Number      fx, fy, fz,\r
+        d000, d001, d010, d011,\r
+        d100, d101, d110, d111,\r
+        dx00, dx01, dx10, dx11,\r
+        dxy0, dxy1, dxyz;\r
+\r
+    TotalOut   = p -> nOutputs;\r
+\r
+    // We need some clipping here\r
+    px = Input[0];\r
+    py = Input[1];\r
+    pz = Input[2];\r
+\r
+    if (px < 0) px = 0;\r
+    if (px > 1) px = 1;\r
+    if (py < 0) py = 0;\r
+    if (py > 1) py = 1;\r
+    if (pz < 0) pz = 0;\r
+    if (pz > 1) pz = 1;\r
+\r
+    px *= p->Domain[0];\r
+    py *= p->Domain[1];\r
+    pz *= p->Domain[2];\r
+\r
+    x0 = (int) _cmsQuickFloor(px); fx = px - (cmsFloat32Number) x0;\r
+    y0 = (int) _cmsQuickFloor(py); fy = py - (cmsFloat32Number) y0;\r
+    z0 = (int) _cmsQuickFloor(pz); fz = pz - (cmsFloat32Number) z0;\r
+\r
+    X0 = p -> opta[2] * x0;\r
+    X1 = X0 + (Input[0] >= 1.0 ? 0 : p->opta[2]);\r
+\r
+    Y0 = p -> opta[1] * y0;\r
+    Y1 = Y0 + (Input[1] >= 1.0 ? 0 : p->opta[1]);\r
+\r
+    Z0 = p -> opta[0] * z0;\r
+    Z1 = Z0 + (Input[2] >= 1.0 ? 0 : p->opta[0]);\r
+\r
+    for (OutChan = 0; OutChan < TotalOut; OutChan++) {\r
+\r
+        d000 = DENS(X0, Y0, Z0);\r
+        d001 = DENS(X0, Y0, Z1);\r
+        d010 = DENS(X0, Y1, Z0);\r
+        d011 = DENS(X0, Y1, Z1);\r
+\r
+        d100 = DENS(X1, Y0, Z0);\r
+        d101 = DENS(X1, Y0, Z1);\r
+        d110 = DENS(X1, Y1, Z0);\r
+        d111 = DENS(X1, Y1, Z1);\r
+\r
+\r
+        dx00 = LERP(fx, d000, d100);\r
+        dx01 = LERP(fx, d001, d101);\r
+        dx10 = LERP(fx, d010, d110);\r
+        dx11 = LERP(fx, d011, d111);\r
+\r
+        dxy0 = LERP(fy, dx00, dx10);\r
+        dxy1 = LERP(fy, dx01, dx11);\r
+\r
+        dxyz = LERP(fz, dxy0, dxy1);\r
+\r
+        Output[OutChan] = dxyz;\r
+    }\r
+\r
+\r
+#   undef LERP\r
+#   undef DENS\r
+}\r
+\r
+// Trilinear interpolation (16 bits) - optimized version\r
+static\r
+void TrilinearInterp16(register const cmsUInt16Number Input[],\r
+                       register cmsUInt16Number Output[],\r
+                       register const cmsInterpParams* p)\r
+\r
+{\r
+#define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])\r
+#define LERP(a,l,h)     (cmsUInt16Number) (l + ROUND_FIXED_TO_INT(((h-l)*a)))\r
+\r
+           const cmsUInt16Number* LutTable = (cmsUInt16Number*) p ->Table;\r
+           int        OutChan, TotalOut;\r
+           cmsS15Fixed16Number    fx, fy, fz;\r
+  register int        rx, ry, rz;\r
+           int        x0, y0, z0;\r
+  register int        X0, X1, Y0, Y1, Z0, Z1;\r
+           int        d000, d001, d010, d011,\r
+                      d100, d101, d110, d111,\r
+                      dx00, dx01, dx10, dx11,\r
+                      dxy0, dxy1, dxyz;\r
+\r
+    TotalOut   = p -> nOutputs;\r
+\r
+    fx = _cmsToFixedDomain((int) Input[0] * p -> Domain[0]);\r
+    x0  = FIXED_TO_INT(fx);\r
+    rx  = FIXED_REST_TO_INT(fx);    // Rest in 0..1.0 domain\r
+\r
+\r
+    fy = _cmsToFixedDomain((int) Input[1] * p -> Domain[1]);\r
+    y0  = FIXED_TO_INT(fy);\r
+    ry  = FIXED_REST_TO_INT(fy);\r
+\r
+    fz = _cmsToFixedDomain((int) Input[2] * p -> Domain[2]);\r
+    z0 = FIXED_TO_INT(fz);\r
+    rz = FIXED_REST_TO_INT(fz);\r
+\r
+\r
+    X0 = p -> opta[2] * x0;\r
+    X1 = X0 + (Input[0] == 0xFFFFU ? 0 : p->opta[2]);\r
+\r
+    Y0 = p -> opta[1] * y0;\r
+    Y1 = Y0 + (Input[1] == 0xFFFFU ? 0 : p->opta[1]);\r
+\r
+    Z0 = p -> opta[0] * z0;\r
+    Z1 = Z0 + (Input[2] == 0xFFFFU ? 0 : p->opta[0]);\r
+\r
+    for (OutChan = 0; OutChan < TotalOut; OutChan++) {\r
+\r
+        d000 = DENS(X0, Y0, Z0);\r
+        d001 = DENS(X0, Y0, Z1);\r
+        d010 = DENS(X0, Y1, Z0);\r
+        d011 = DENS(X0, Y1, Z1);\r
+\r
+        d100 = DENS(X1, Y0, Z0);\r
+        d101 = DENS(X1, Y0, Z1);\r
+        d110 = DENS(X1, Y1, Z0);\r
+        d111 = DENS(X1, Y1, Z1);\r
+\r
+\r
+        dx00 = LERP(rx, d000, d100);\r
+        dx01 = LERP(rx, d001, d101);\r
+        dx10 = LERP(rx, d010, d110);\r
+        dx11 = LERP(rx, d011, d111);\r
+\r
+        dxy0 = LERP(ry, dx00, dx10);\r
+        dxy1 = LERP(ry, dx01, dx11);\r
+\r
+        dxyz = LERP(rz, dxy0, dxy1);\r
+\r
+        Output[OutChan] = (cmsUInt16Number) dxyz;\r
+    }\r
+\r
+\r
+#   undef LERP\r
+#   undef DENS\r
+}\r
+\r
+\r
+// Tetrahedral interpolation, using Sakamoto algorithm.\r
+#define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])\r
+static\r
+void TetrahedralInterpFloat(const cmsFloat32Number Input[],\r
+                            cmsFloat32Number Output[],\r
+                            const cmsInterpParams* p)\r
+{\r
+    const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;\r
+    cmsFloat32Number     px, py, pz;\r
+    int        x0, y0, z0,\r
+               X0, Y0, Z0, X1, Y1, Z1;\r
+    cmsFloat32Number     rx, ry, rz;\r
+    cmsFloat32Number     c0, c1=0, c2=0, c3=0;\r
+    int                  OutChan, TotalOut;\r
+\r
+    TotalOut   = p -> nOutputs;\r
+\r
+    // We need some clipping here\r
+    px = Input[0];\r
+    py = Input[1];\r
+    pz = Input[2];\r
+\r
+    if (px < 0) px = 0;\r
+    if (px > 1) px = 1;\r
+    if (py < 0) py = 0;\r
+    if (py > 1) py = 1;\r
+    if (pz < 0) pz = 0;\r
+    if (pz > 1) pz = 1;\r
+\r
+    px *= p->Domain[0];\r
+    py *= p->Domain[1];\r
+    pz *= p->Domain[2];\r
+\r
+    x0 = (int) _cmsQuickFloor(px); rx = (px - (cmsFloat32Number) x0);\r
+    y0 = (int) _cmsQuickFloor(py); ry = (py - (cmsFloat32Number) y0);\r
+    z0 = (int) _cmsQuickFloor(pz); rz = (pz - (cmsFloat32Number) z0);\r
+\r
+\r
+    X0 = p -> opta[2] * x0;\r
+    X1 = X0 + (Input[0] >= 1.0 ? 0 : p->opta[2]);\r
+\r
+    Y0 = p -> opta[1] * y0;\r
+    Y1 = Y0 + (Input[1] >= 1.0 ? 0 : p->opta[1]);\r
+\r
+    Z0 = p -> opta[0] * z0;\r
+    Z1 = Z0 + (Input[2] >= 1.0 ? 0 : p->opta[0]);\r
+\r
+    for (OutChan=0; OutChan < TotalOut; OutChan++) {\r
+\r
+       // These are the 6 Tetrahedral\r
+\r
+        c0 = DENS(X0, Y0, Z0);\r
+\r
+        if (rx >= ry && ry >= rz) {\r
+\r
+            c1 = DENS(X1, Y0, Z0) - c0;\r
+            c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);\r
+            c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);\r
+\r
+        }\r
+        else\r
+            if (rx >= rz && rz >= ry) {\r
+\r
+                c1 = DENS(X1, Y0, Z0) - c0;\r
+                c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);\r
+                c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);\r
+\r
+            }\r
+            else\r
+                if (rz >= rx && rx >= ry) {\r
+\r
+                    c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);\r
+                    c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);\r
+                    c3 = DENS(X0, Y0, Z1) - c0;\r
+\r
+                }\r
+                else\r
+                    if (ry >= rx && rx >= rz) {\r
+\r
+                        c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);\r
+                        c2 = DENS(X0, Y1, Z0) - c0;\r
+                        c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);\r
+\r
+                    }\r
+                    else\r
+                        if (ry >= rz && rz >= rx) {\r
+\r
+                            c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);\r
+                            c2 = DENS(X0, Y1, Z0) - c0;\r
+                            c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);\r
+\r
+                        }\r
+                        else\r
+                            if (rz >= ry && ry >= rx) {\r
+\r
+                                c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);\r
+                                c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);\r
+                                c3 = DENS(X0, Y0, Z1) - c0;\r
+\r
+                            }\r
+                            else  {\r
+                                c1 = c2 = c3 = 0;\r
+                            }\r
+\r
+       Output[OutChan] = c0 + c1 * rx + c2 * ry + c3 * rz;\r
+       }\r
+\r
+}\r
+\r
+#undef DENS\r
+\r
+\r
+\r
+\r
+static\r
+void TetrahedralInterp16(register const cmsUInt16Number Input[],\r
+                         register cmsUInt16Number Output[],\r
+                         register const cmsInterpParams* p)\r
+{\r
+    const cmsUInt16Number* LutTable = (cmsUInt16Number*) p -> Table;\r
+    cmsS15Fixed16Number fx, fy, fz;\r
+    cmsS15Fixed16Number rx, ry, rz;\r
+    int x0, y0, z0;\r
+    cmsS15Fixed16Number c0, c1, c2, c3, Rest;\r
+    cmsS15Fixed16Number X0, X1, Y0, Y1, Z0, Z1;\r
+    cmsUInt32Number TotalOut = p -> nOutputs;\r
+\r
+    fx = _cmsToFixedDomain((int) Input[0] * p -> Domain[0]);\r
+    fy = _cmsToFixedDomain((int) Input[1] * p -> Domain[1]);\r
+    fz = _cmsToFixedDomain((int) Input[2] * p -> Domain[2]);\r
+\r
+    x0 = FIXED_TO_INT(fx);\r
+    y0 = FIXED_TO_INT(fy);\r
+    z0 = FIXED_TO_INT(fz);\r
+\r
+    rx = FIXED_REST_TO_INT(fx);\r
+    ry = FIXED_REST_TO_INT(fy);\r
+    rz = FIXED_REST_TO_INT(fz);\r
+\r
+    X0 = p -> opta[2] * x0;\r
+    X1 = (Input[0] == 0xFFFFU ? 0 : p->opta[2]);\r
+\r
+    Y0 = p -> opta[1] * y0;\r
+    Y1 = (Input[1] == 0xFFFFU ? 0 : p->opta[1]);\r
+\r
+    Z0 = p -> opta[0] * z0;\r
+    Z1 = (Input[2] == 0xFFFFU ? 0 : p->opta[0]);\r
+\r
+    LutTable = &LutTable[X0+Y0+Z0];\r
+\r
+    // Output should be computed as x = ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest))\r
+    // which expands as: x = (Rest + ((Rest+0x7fff)/0xFFFF) + 0x8000)>>16\r
+    // This can be replaced by: t = Rest+0x8001, x = (t + (t>>16))>>16\r
+    // at the cost of being off by one at 7fff and 17ffe.\r
+\r
+    if (rx >= ry) {\r
+        if (ry >= rz) {\r
+            Y1 += X1;\r
+            Z1 += Y1;\r
+            for (; TotalOut; TotalOut--) {\r
+                c1 = LutTable[X1];\r
+                c2 = LutTable[Y1];\r
+                c3 = LutTable[Z1];\r
+                c0 = *LutTable++;\r
+                c3 -= c2;\r
+                c2 -= c1;\r
+                c1 -= c0;\r
+                Rest = c1 * rx + c2 * ry + c3 * rz + 0x8001;\r
+                *Output++ = (cmsUInt16Number) c0 + ((Rest + (Rest>>16))>>16);\r
+            }\r
+        } else if (rz >= rx) {\r
+            X1 += Z1;\r
+            Y1 += X1;\r
+            for (; TotalOut; TotalOut--) {\r
+                c1 = LutTable[X1];\r
+                c2 = LutTable[Y1];\r
+                c3 = LutTable[Z1];\r
+                c0 = *LutTable++;\r
+                c2 -= c1;\r
+                c1 -= c3;\r
+                c3 -= c0;\r
+                Rest = c1 * rx + c2 * ry + c3 * rz + 0x8001;\r
+                *Output++ = (cmsUInt16Number) c0 + ((Rest + (Rest>>16))>>16);\r
+            }\r
+        } else {\r
+            Z1 += X1;\r
+            Y1 += Z1;\r
+            for (; TotalOut; TotalOut--) {\r
+                c1 = LutTable[X1];\r
+                c2 = LutTable[Y1];\r
+                c3 = LutTable[Z1];\r
+                c0 = *LutTable++;\r
+                c2 -= c3;\r
+                c3 -= c1;\r
+                c1 -= c0;\r
+                Rest = c1 * rx + c2 * ry + c3 * rz + 0x8001;\r
+                *Output++ = (cmsUInt16Number) c0 + ((Rest + (Rest>>16))>>16);\r
+            }\r
+        }\r
+    } else {\r
+        if (rx >= rz) {\r
+            X1 += Y1;\r
+            Z1 += X1;\r
+            for (; TotalOut; TotalOut--) {\r
+                c1 = LutTable[X1];\r
+                c2 = LutTable[Y1];\r
+                c3 = LutTable[Z1];\r
+                c0 = *LutTable++;\r
+                c3 -= c1;\r
+                c1 -= c2;\r
+                c2 -= c0;\r
+                Rest = c1 * rx + c2 * ry + c3 * rz + 0x8001;\r
+                *Output++ = (cmsUInt16Number) c0 + ((Rest + (Rest>>16))>>16);\r
+            }\r
+        } else if (ry >= rz) {\r
+            Z1 += Y1;\r
+            X1 += Z1;\r
+            for (; TotalOut; TotalOut--) {\r
+                c1 = LutTable[X1];\r
+                c2 = LutTable[Y1];\r
+                c3 = LutTable[Z1];\r
+                c0 = *LutTable++;\r
+                c1 -= c3;\r
+                c3 -= c2;\r
+                c2 -= c0;\r
+                Rest = c1 * rx + c2 * ry + c3 * rz + 0x8001;\r
+                *Output++ = (cmsUInt16Number) c0 + ((Rest + (Rest>>16))>>16);\r
+            }\r
+        } else {\r
+            Y1 += Z1;\r
+            X1 += Y1;\r
+            for (; TotalOut; TotalOut--) {\r
+                c1 = LutTable[X1];\r
+                c2 = LutTable[Y1];\r
+                c3 = LutTable[Z1];\r
+                c0 = *LutTable++;\r
+                c1 -= c2;\r
+                c2 -= c3;\r
+                c3 -= c0;\r
+                Rest = c1 * rx + c2 * ry + c3 * rz + 0x8001;\r
+                *Output++ = (cmsUInt16Number) c0 + ((Rest + (Rest>>16))>>16);\r
+            }\r
+        }\r
+    }\r
+}\r
+\r
+\r
+#define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])\r
+static\r
+void Eval4Inputs(register const cmsUInt16Number Input[],\r
+                     register cmsUInt16Number Output[],\r
+                     register const cmsInterpParams* p16)\r
+{\r
+    const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;\r
+    cmsS15Fixed16Number fk;\r
+    cmsS15Fixed16Number k0, rk;\r
+    int K0, K1;\r
+    cmsS15Fixed16Number    fx, fy, fz;\r
+    cmsS15Fixed16Number    rx, ry, rz;\r
+    int                    x0, y0, z0;\r
+    cmsS15Fixed16Number    X0, X1, Y0, Y1, Z0, Z1;\r
+    cmsUInt32Number i;\r
+    cmsS15Fixed16Number    c0, c1, c2, c3, Rest;\r
+    cmsUInt32Number        OutChan;\r
+    cmsUInt16Number        Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+\r
+\r
+    fk  = _cmsToFixedDomain((int) Input[0] * p16 -> Domain[0]);\r
+    fx  = _cmsToFixedDomain((int) Input[1] * p16 -> Domain[1]);\r
+    fy  = _cmsToFixedDomain((int) Input[2] * p16 -> Domain[2]);\r
+    fz  = _cmsToFixedDomain((int) Input[3] * p16 -> Domain[3]);\r
+\r
+    k0  = FIXED_TO_INT(fk);\r
+    x0  = FIXED_TO_INT(fx);\r
+    y0  = FIXED_TO_INT(fy);\r
+    z0  = FIXED_TO_INT(fz);\r
+\r
+    rk  = FIXED_REST_TO_INT(fk);\r
+    rx  = FIXED_REST_TO_INT(fx);\r
+    ry  = FIXED_REST_TO_INT(fy);\r
+    rz  = FIXED_REST_TO_INT(fz);\r
+\r
+    K0 = p16 -> opta[3] * k0;\r
+    K1 = K0 + (Input[0] == 0xFFFFU ? 0 : p16->opta[3]);\r
+\r
+    X0 = p16 -> opta[2] * x0;\r
+    X1 = X0 + (Input[1] == 0xFFFFU ? 0 : p16->opta[2]);\r
+\r
+    Y0 = p16 -> opta[1] * y0;\r
+    Y1 = Y0 + (Input[2] == 0xFFFFU ? 0 : p16->opta[1]);\r
+\r
+    Z0 = p16 -> opta[0] * z0;\r
+    Z1 = Z0 + (Input[3] == 0xFFFFU ? 0 : p16->opta[0]);\r
+\r
+    LutTable = (cmsUInt16Number*) p16 -> Table;\r
+    LutTable += K0;\r
+\r
+    for (OutChan=0; OutChan < p16 -> nOutputs; OutChan++) {\r
+\r
+        c0 = DENS(X0, Y0, Z0);\r
+\r
+        if (rx >= ry && ry >= rz) {\r
+\r
+            c1 = DENS(X1, Y0, Z0) - c0;\r
+            c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);\r
+            c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);\r
+\r
+        }\r
+        else\r
+            if (rx >= rz && rz >= ry) {\r
+\r
+                c1 = DENS(X1, Y0, Z0) - c0;\r
+                c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);\r
+                c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);\r
+\r
+            }\r
+            else\r
+                if (rz >= rx && rx >= ry) {\r
+\r
+                    c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);\r
+                    c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);\r
+                    c3 = DENS(X0, Y0, Z1) - c0;\r
+\r
+                }\r
+                else\r
+                    if (ry >= rx && rx >= rz) {\r
+\r
+                        c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);\r
+                        c2 = DENS(X0, Y1, Z0) - c0;\r
+                        c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);\r
+\r
+                    }\r
+                    else\r
+                        if (ry >= rz && rz >= rx) {\r
+\r
+                            c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);\r
+                            c2 = DENS(X0, Y1, Z0) - c0;\r
+                            c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);\r
+\r
+                        }\r
+                        else\r
+                            if (rz >= ry && ry >= rx) {\r
+\r
+                                c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);\r
+                                c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);\r
+                                c3 = DENS(X0, Y0, Z1) - c0;\r
+\r
+                            }\r
+                            else  {\r
+                                c1 = c2 = c3 = 0;\r
+                            }\r
+\r
+                            Rest = c1 * rx + c2 * ry + c3 * rz;\r
+\r
+                            Tmp1[OutChan] = (cmsUInt16Number) c0 + ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest));\r
+    }\r
+\r
+\r
+    LutTable = (cmsUInt16Number*) p16 -> Table;\r
+    LutTable += K1;\r
+\r
+    for (OutChan=0; OutChan < p16 -> nOutputs; OutChan++) {\r
+\r
+        c0 = DENS(X0, Y0, Z0);\r
+\r
+        if (rx >= ry && ry >= rz) {\r
+\r
+            c1 = DENS(X1, Y0, Z0) - c0;\r
+            c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);\r
+            c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);\r
+\r
+        }\r
+        else\r
+            if (rx >= rz && rz >= ry) {\r
+\r
+                c1 = DENS(X1, Y0, Z0) - c0;\r
+                c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);\r
+                c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);\r
+\r
+            }\r
+            else\r
+                if (rz >= rx && rx >= ry) {\r
+\r
+                    c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);\r
+                    c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);\r
+                    c3 = DENS(X0, Y0, Z1) - c0;\r
+\r
+                }\r
+                else\r
+                    if (ry >= rx && rx >= rz) {\r
+\r
+                        c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);\r
+                        c2 = DENS(X0, Y1, Z0) - c0;\r
+                        c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);\r
+\r
+                    }\r
+                    else\r
+                        if (ry >= rz && rz >= rx) {\r
+\r
+                            c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);\r
+                            c2 = DENS(X0, Y1, Z0) - c0;\r
+                            c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);\r
+\r
+                        }\r
+                        else\r
+                            if (rz >= ry && ry >= rx) {\r
+\r
+                                c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);\r
+                                c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);\r
+                                c3 = DENS(X0, Y0, Z1) - c0;\r
+\r
+                            }\r
+                            else  {\r
+                                c1 = c2 = c3 = 0;\r
+                            }\r
+\r
+                            Rest = c1 * rx + c2 * ry + c3 * rz;\r
+\r
+                            Tmp2[OutChan] = (cmsUInt16Number) c0 + ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest));\r
+    }\r
+\r
+\r
+\r
+    for (i=0; i < p16 -> nOutputs; i++) {\r
+        Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);\r
+    }\r
+}\r
+#undef DENS\r
+\r
+\r
+// For more that 3 inputs (i.e., CMYK)\r
+// evaluate two 3-dimensional interpolations and then linearly interpolate between them.\r
+\r
+\r
+static\r
+void Eval4InputsFloat(const cmsFloat32Number Input[],\r
+                      cmsFloat32Number Output[],\r
+                      const cmsInterpParams* p)\r
+{\r
+       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;\r
+       cmsFloat32Number rest;\r
+       cmsFloat32Number pk;\r
+       int k0, K0, K1;\r
+       const cmsFloat32Number* T;\r
+       cmsUInt32Number i;\r
+       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+       cmsInterpParams p1;\r
+\r
+\r
+       pk = Input[0] * p->Domain[0];\r
+       k0 = _cmsQuickFloor(pk);\r
+       rest = pk - (cmsFloat32Number) k0;\r
+\r
+       K0 = p -> opta[3] * k0;\r
+       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[3]);\r
+\r
+       p1 = *p;\r
+       memmove(&p1.Domain[0], &p ->Domain[1], 3*sizeof(cmsUInt32Number));\r
+\r
+       T = LutTable + K0;\r
+       p1.Table = T;\r
+\r
+       TetrahedralInterpFloat(Input + 1,  Tmp1, &p1);\r
+\r
+       T = LutTable + K1;\r
+       p1.Table = T;\r
+       TetrahedralInterpFloat(Input + 1,  Tmp2, &p1);\r
+\r
+       for (i=0; i < p -> nOutputs; i++)\r
+       {\r
+              cmsFloat32Number y0 = Tmp1[i];\r
+              cmsFloat32Number y1 = Tmp2[i];\r
+\r
+              Output[i] = y0 + (y1 - y0) * rest;\r
+       }\r
+}\r
+\r
+\r
+static\r
+void Eval5Inputs(register const cmsUInt16Number Input[],\r
+                 register cmsUInt16Number Output[],\r
+\r
+                 register const cmsInterpParams* p16)\r
+{\r
+       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;\r
+       cmsS15Fixed16Number fk;\r
+       cmsS15Fixed16Number k0, rk;\r
+       int K0, K1;\r
+       const cmsUInt16Number* T;\r
+       cmsUInt32Number i;\r
+       cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+       cmsInterpParams p1;\r
+\r
+\r
+       fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);\r
+       k0 = FIXED_TO_INT(fk);\r
+       rk = FIXED_REST_TO_INT(fk);\r
+\r
+       K0 = p16 -> opta[4] * k0;\r
+       K1 = p16 -> opta[4] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));\r
+\r
+       p1 = *p16;\r
+       memmove(&p1.Domain[0], &p16 ->Domain[1], 4*sizeof(cmsUInt32Number));\r
+\r
+       T = LutTable + K0;\r
+       p1.Table = T;\r
+\r
+       Eval4Inputs(Input + 1, Tmp1, &p1);\r
+\r
+       T = LutTable + K1;\r
+       p1.Table = T;\r
+\r
+       Eval4Inputs(Input + 1, Tmp2, &p1);\r
+\r
+       for (i=0; i < p16 -> nOutputs; i++) {\r
+\r
+              Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);\r
+       }\r
+\r
+}\r
+\r
+\r
+static\r
+void Eval5InputsFloat(const cmsFloat32Number Input[],\r
+                      cmsFloat32Number Output[],\r
+                      const cmsInterpParams* p)\r
+{\r
+       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;\r
+       cmsFloat32Number rest;\r
+       cmsFloat32Number pk;\r
+       int k0, K0, K1;\r
+       const cmsFloat32Number* T;\r
+       cmsUInt32Number i;\r
+       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+       cmsInterpParams p1;\r
+\r
+       pk = Input[0] * p->Domain[0];\r
+       k0 = _cmsQuickFloor(pk);\r
+       rest = pk - (cmsFloat32Number) k0;\r
+\r
+       K0 = p -> opta[4] * k0;\r
+       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[4]);\r
+\r
+       p1 = *p;\r
+       memmove(&p1.Domain[0], &p ->Domain[1], 4*sizeof(cmsUInt32Number));\r
+\r
+       T = LutTable + K0;\r
+       p1.Table = T;\r
+\r
+       Eval4InputsFloat(Input + 1,  Tmp1, &p1);\r
+\r
+       T = LutTable + K1;\r
+       p1.Table = T;\r
+\r
+       Eval4InputsFloat(Input + 1,  Tmp2, &p1);\r
+\r
+       for (i=0; i < p -> nOutputs; i++) {\r
+\r
+              cmsFloat32Number y0 = Tmp1[i];\r
+              cmsFloat32Number y1 = Tmp2[i];\r
+\r
+              Output[i] = y0 + (y1 - y0) * rest;\r
+       }\r
+}\r
+\r
+\r
+\r
+static\r
+void Eval6Inputs(register const cmsUInt16Number Input[],\r
+                 register cmsUInt16Number Output[],\r
+                 register const cmsInterpParams* p16)\r
+{\r
+       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;\r
+       cmsS15Fixed16Number fk;\r
+       cmsS15Fixed16Number k0, rk;\r
+       int K0, K1;\r
+       const cmsUInt16Number* T;\r
+       cmsUInt32Number i;\r
+       cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+       cmsInterpParams p1;\r
+\r
+       fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);\r
+       k0 = FIXED_TO_INT(fk);\r
+       rk = FIXED_REST_TO_INT(fk);\r
+\r
+       K0 = p16 -> opta[5] * k0;\r
+       K1 = p16 -> opta[5] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));\r
+\r
+       p1 = *p16;\r
+       memmove(&p1.Domain[0], &p16 ->Domain[1], 5*sizeof(cmsUInt32Number));\r
+\r
+       T = LutTable + K0;\r
+       p1.Table = T;\r
+\r
+       Eval5Inputs(Input + 1, Tmp1, &p1);\r
+\r
+       T = LutTable + K1;\r
+       p1.Table = T;\r
+\r
+       Eval5Inputs(Input + 1, Tmp2, &p1);\r
+\r
+       for (i=0; i < p16 -> nOutputs; i++) {\r
+\r
+              Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);\r
+       }\r
+\r
+}\r
+\r
+\r
+static\r
+void Eval6InputsFloat(const cmsFloat32Number Input[],\r
+                      cmsFloat32Number Output[],\r
+                      const cmsInterpParams* p)\r
+{\r
+       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;\r
+       cmsFloat32Number rest;\r
+       cmsFloat32Number pk;\r
+       int k0, K0, K1;\r
+       const cmsFloat32Number* T;\r
+       cmsUInt32Number i;\r
+       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+       cmsInterpParams p1;\r
+\r
+       pk = Input[0] * p->Domain[0];\r
+       k0 = _cmsQuickFloor(pk);\r
+       rest = pk - (cmsFloat32Number) k0;\r
+\r
+       K0 = p -> opta[5] * k0;\r
+       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[5]);\r
+\r
+       p1 = *p;\r
+       memmove(&p1.Domain[0], &p ->Domain[1], 5*sizeof(cmsUInt32Number));\r
+\r
+       T = LutTable + K0;\r
+       p1.Table = T;\r
+\r
+       Eval5InputsFloat(Input + 1,  Tmp1, &p1);\r
+\r
+       T = LutTable + K1;\r
+       p1.Table = T;\r
+\r
+       Eval5InputsFloat(Input + 1,  Tmp2, &p1);\r
+\r
+       for (i=0; i < p -> nOutputs; i++) {\r
+\r
+              cmsFloat32Number y0 = Tmp1[i];\r
+              cmsFloat32Number y1 = Tmp2[i];\r
+\r
+              Output[i] = y0 + (y1 - y0) * rest;\r
+       }\r
+}\r
+\r
+\r
+static\r
+void Eval7Inputs(register const cmsUInt16Number Input[],\r
+                 register cmsUInt16Number Output[],\r
+                 register const cmsInterpParams* p16)\r
+{\r
+       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;\r
+       cmsS15Fixed16Number fk;\r
+       cmsS15Fixed16Number k0, rk;\r
+       int K0, K1;\r
+       const cmsUInt16Number* T;\r
+       cmsUInt32Number i;\r
+       cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+       cmsInterpParams p1;\r
+\r
+\r
+       fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);\r
+       k0 = FIXED_TO_INT(fk);\r
+       rk = FIXED_REST_TO_INT(fk);\r
+\r
+       K0 = p16 -> opta[6] * k0;\r
+       K1 = p16 -> opta[6] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));\r
+\r
+       p1 = *p16;\r
+       memmove(&p1.Domain[0], &p16 ->Domain[1], 6*sizeof(cmsUInt32Number));\r
+\r
+       T = LutTable + K0;\r
+       p1.Table = T;\r
+\r
+       Eval6Inputs(Input + 1, Tmp1, &p1);\r
+\r
+       T = LutTable + K1;\r
+       p1.Table = T;\r
+\r
+       Eval6Inputs(Input + 1, Tmp2, &p1);\r
+\r
+       for (i=0; i < p16 -> nOutputs; i++) {\r
+              Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);\r
+       }\r
+}\r
+\r
+\r
+static\r
+void Eval7InputsFloat(const cmsFloat32Number Input[],\r
+                      cmsFloat32Number Output[],\r
+                      const cmsInterpParams* p)\r
+{\r
+       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;\r
+       cmsFloat32Number rest;\r
+       cmsFloat32Number pk;\r
+       int k0, K0, K1;\r
+       const cmsFloat32Number* T;\r
+       cmsUInt32Number i;\r
+       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+       cmsInterpParams p1;\r
+\r
+       pk = Input[0] * p->Domain[0];\r
+       k0 = _cmsQuickFloor(pk);\r
+       rest = pk - (cmsFloat32Number) k0;\r
+\r
+       K0 = p -> opta[6] * k0;\r
+       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[6]);\r
+\r
+       p1 = *p;\r
+       memmove(&p1.Domain[0], &p ->Domain[1], 6*sizeof(cmsUInt32Number));\r
+\r
+       T = LutTable + K0;\r
+       p1.Table = T;\r
+\r
+       Eval6InputsFloat(Input + 1,  Tmp1, &p1);\r
+\r
+       T = LutTable + K1;\r
+       p1.Table = T;\r
+\r
+       Eval6InputsFloat(Input + 1,  Tmp2, &p1);\r
+\r
+\r
+       for (i=0; i < p -> nOutputs; i++) {\r
+\r
+              cmsFloat32Number y0 = Tmp1[i];\r
+              cmsFloat32Number y1 = Tmp2[i];\r
+\r
+              Output[i] = y0 + (y1 - y0) * rest;\r
+\r
+       }\r
+}\r
+\r
+static\r
+void Eval8Inputs(register const cmsUInt16Number Input[],\r
+                 register cmsUInt16Number Output[],\r
+                 register const cmsInterpParams* p16)\r
+{\r
+       const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;\r
+       cmsS15Fixed16Number fk;\r
+       cmsS15Fixed16Number k0, rk;\r
+       int K0, K1;\r
+       const cmsUInt16Number* T;\r
+       cmsUInt32Number i;\r
+       cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+       cmsInterpParams p1;\r
+\r
+       fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);\r
+       k0 = FIXED_TO_INT(fk);\r
+       rk = FIXED_REST_TO_INT(fk);\r
+\r
+       K0 = p16 -> opta[7] * k0;\r
+       K1 = p16 -> opta[7] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));\r
+\r
+       p1 = *p16;\r
+       memmove(&p1.Domain[0], &p16 ->Domain[1], 7*sizeof(cmsUInt32Number));\r
+\r
+       T = LutTable + K0;\r
+       p1.Table = T;\r
+\r
+       Eval7Inputs(Input + 1, Tmp1, &p1);\r
+\r
+       T = LutTable + K1;\r
+       p1.Table = T;\r
+       Eval7Inputs(Input + 1, Tmp2, &p1);\r
+\r
+       for (i=0; i < p16 -> nOutputs; i++) {\r
+              Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);\r
+       }\r
+}\r
+\r
+\r
+\r
+static\r
+void Eval8InputsFloat(const cmsFloat32Number Input[],\r
+                      cmsFloat32Number Output[],\r
+                      const cmsInterpParams* p)\r
+{\r
+       const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;\r
+       cmsFloat32Number rest;\r
+       cmsFloat32Number pk;\r
+       int k0, K0, K1;\r
+       const cmsFloat32Number* T;\r
+       cmsUInt32Number i;\r
+       cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\r
+       cmsInterpParams p1;\r
+\r
+       pk = Input[0] * p->Domain[0];\r
+       k0 = _cmsQuickFloor(pk);\r
+       rest = pk - (cmsFloat32Number) k0;\r
+\r
+       K0 = p -> opta[7] * k0;\r
+       K1 = K0 + (Input[0] >= 1.0 ? 0 : p->opta[7]);\r
+\r
+       p1 = *p;\r
+       memmove(&p1.Domain[0], &p ->Domain[1], 7*sizeof(cmsUInt32Number));\r
+\r
+       T = LutTable + K0;\r
+       p1.Table = T;\r
+\r
+       Eval7InputsFloat(Input + 1,  Tmp1, &p1);\r
+\r
+       T = LutTable + K1;\r
+       p1.Table = T;\r
+\r
+       Eval7InputsFloat(Input + 1,  Tmp2, &p1);\r
+\r
+\r
+       for (i=0; i < p -> nOutputs; i++) {\r
+\r
+              cmsFloat32Number y0 = Tmp1[i];\r
+              cmsFloat32Number y1 = Tmp2[i];\r
+\r
+              Output[i] = y0 + (y1 - y0) * rest;\r
+       }\r
+}\r
+\r
+// The default factory\r
+static\r
+cmsInterpFunction DefaultInterpolatorsFactory(cmsUInt32Number nInputChannels, cmsUInt32Number nOutputChannels, cmsUInt32Number dwFlags)\r
+{\r
+\r
+    cmsInterpFunction Interpolation;\r
+    cmsBool  IsFloat     = (dwFlags & CMS_LERP_FLAGS_FLOAT);\r
+    cmsBool  IsTrilinear = (dwFlags & CMS_LERP_FLAGS_TRILINEAR);\r
+\r
+    memset(&Interpolation, 0, sizeof(Interpolation));\r
+\r
+    // Safety check\r
+    if (nInputChannels >= 4 && nOutputChannels >= MAX_STAGE_CHANNELS)\r
+        return Interpolation;\r
+\r
+    switch (nInputChannels) {\r
+\r
+           case 1: // Gray LUT / linear\r
+\r
+               if (nOutputChannels == 1) {\r
+\r
+                   if (IsFloat)\r
+                       Interpolation.LerpFloat = LinLerp1Dfloat;\r
+                   else\r
+                       Interpolation.Lerp16 = LinLerp1D;\r
+\r
+               }\r
+               else {\r
+\r
+                   if (IsFloat)\r
+                       Interpolation.LerpFloat = Eval1InputFloat;\r
+                   else\r
+                       Interpolation.Lerp16 = Eval1Input;\r
+               }\r
+               break;\r
+\r
+           case 2: // Duotone\r
+               if (IsFloat)\r
+                      Interpolation.LerpFloat =  BilinearInterpFloat;\r
+               else\r
+                      Interpolation.Lerp16    =  BilinearInterp16;\r
+               break;\r
+\r
+           case 3:  // RGB et al\r
+\r
+               if (IsTrilinear) {\r
+\r
+                   if (IsFloat)\r
+                       Interpolation.LerpFloat = TrilinearInterpFloat;\r
+                   else\r
+                       Interpolation.Lerp16 = TrilinearInterp16;\r
+               }\r
+               else {\r
+\r
+                   if (IsFloat)\r
+                       Interpolation.LerpFloat = TetrahedralInterpFloat;\r
+                   else {\r
+\r
+                       Interpolation.Lerp16 = TetrahedralInterp16;\r
+                   }\r
+               }\r
+               break;\r
+\r
+           case 4:  // CMYK lut\r
+\r
+               if (IsFloat)\r
+                   Interpolation.LerpFloat =  Eval4InputsFloat;\r
+               else\r
+                   Interpolation.Lerp16    =  Eval4Inputs;\r
+               break;\r
+\r
+           case 5: // 5 Inks\r
+               if (IsFloat)\r
+                   Interpolation.LerpFloat =  Eval5InputsFloat;\r
+               else\r
+                   Interpolation.Lerp16    =  Eval5Inputs;\r
+               break;\r
+\r
+           case 6: // 6 Inks\r
+               if (IsFloat)\r
+                   Interpolation.LerpFloat =  Eval6InputsFloat;\r
+               else\r
+                   Interpolation.Lerp16    =  Eval6Inputs;\r
+               break;\r
+\r
+           case 7: // 7 inks\r
+               if (IsFloat)\r
+                   Interpolation.LerpFloat =  Eval7InputsFloat;\r
+               else\r
+                   Interpolation.Lerp16    =  Eval7Inputs;\r
+               break;\r
+\r
+           case 8: // 8 inks\r
+               if (IsFloat)\r
+                   Interpolation.LerpFloat =  Eval8InputsFloat;\r
+               else\r
+                   Interpolation.Lerp16    =  Eval8Inputs;\r
+               break;\r
+\r
+               break;\r
+\r
+           default:\r
+               Interpolation.Lerp16 = NULL;\r
+    }\r
+\r
+    return Interpolation;\r
+}\r