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