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