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