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