1 //---------------------------------------------------------------------------------
3 // Little Color Management System
4 // Copyright (c) 1998-2023 Marti Maria Saguer
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:
13 // The above copyright notice and this permission notice shall be included in
14 // all copies or substantial portions of the Software.
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.
24 //---------------------------------------------------------------------------------
27 #include "lcms2_internal.h"
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
32 // Some people may want to compile as C++ with all warnings on, in this case make compiler silent
34 # if (_MSC_VER >= 1400)
35 # pragma warning( disable : 4365 )
39 // Interpolation routines by default
40 static cmsInterpFunction
DefaultInterpolatorsFactory(cmsUInt32Number nInputChannels
, cmsUInt32Number nOutputChannels
, cmsUInt32Number dwFlags
);
42 // This is the default factory
43 _cmsInterpPluginChunkType _cmsInterpPluginChunk
= { NULL
};
45 // The interpolation plug-in memory chunk allocator/dup
46 void _cmsAllocInterpPluginChunk(struct _cmsContext_struct
* ctx
, const struct _cmsContext_struct
* src
)
50 _cmsAssert(ctx
!= NULL
);
53 from
= src
->chunks
[InterpPlugin
];
56 static _cmsInterpPluginChunkType InterpPluginChunk
= { NULL
};
58 from
= &InterpPluginChunk
;
61 _cmsAssert(from
!= NULL
);
62 ctx
->chunks
[InterpPlugin
] = _cmsSubAllocDup(ctx
->MemPool
, from
, sizeof(_cmsInterpPluginChunkType
));
67 cmsBool
_cmsRegisterInterpPlugin(cmsContext ContextID
, cmsPluginBase
* Data
)
69 cmsPluginInterpolation
* Plugin
= (cmsPluginInterpolation
*) Data
;
70 _cmsInterpPluginChunkType
* ptr
= (_cmsInterpPluginChunkType
*) _cmsContextGetClientChunk(ContextID
, InterpPlugin
);
74 ptr
->Interpolators
= NULL
;
78 // Set replacement functions
79 ptr
->Interpolators
= Plugin
->InterpolatorsFactory
;
84 // Set the interpolation method
85 cmsBool
_cmsSetInterpolationRoutine(cmsContext ContextID
, cmsInterpParams
* p
)
87 _cmsInterpPluginChunkType
* ptr
= (_cmsInterpPluginChunkType
*) _cmsContextGetClientChunk(ContextID
, InterpPlugin
);
89 p
->Interpolation
.Lerp16
= NULL
;
91 // Invoke factory, possibly in the Plug-in
92 if (ptr
->Interpolators
!= NULL
)
93 p
->Interpolation
= ptr
->Interpolators(p
-> nInputs
, p
->nOutputs
, p
->dwFlags
);
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
);
100 // Check for valid interpolator (we just check one member of the union)
101 if (p
->Interpolation
.Lerp16
== NULL
) {
109 // This function precalculates as many parameters as possible to speed up the interpolation.
110 cmsInterpParams
* _cmsComputeInterpParamsEx(cmsContext ContextID
,
111 const cmsUInt32Number nSamples
[],
112 cmsUInt32Number InputChan
, cmsUInt32Number OutputChan
,
114 cmsUInt32Number dwFlags
)
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
);
125 // Creates an empty object
126 p
= (cmsInterpParams
*) _cmsMallocZero(ContextID
, sizeof(cmsInterpParams
));
127 if (p
== NULL
) return NULL
;
129 // Keep original parameters
130 p
-> dwFlags
= dwFlags
;
131 p
-> nInputs
= InputChan
;
132 p
-> nOutputs
= OutputChan
;
134 p
->ContextID
= ContextID
;
136 // Fill samples per input direction and domain (which is number of nodes minus one)
137 for (i
=0; i
< InputChan
; i
++) {
139 p
-> nSamples
[i
] = nSamples
[i
];
140 p
-> Domain
[i
] = nSamples
[i
] - 1;
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
];
149 if (!_cmsSetInterpolationRoutine(ContextID
, p
)) {
150 cmsSignalError(ContextID
, cmsERROR_UNKNOWN_EXTENSION
, "Unsupported interpolation (%d->%d channels)", InputChan
, OutputChan
);
151 _cmsFree(ContextID
, p
);
160 // This one is a wrapper on the anterior, but assuming all directions have same number of nodes
161 cmsInterpParams
* CMSEXPORT
_cmsComputeInterpParams(cmsContext ContextID
, cmsUInt32Number nSamples
,
162 cmsUInt32Number InputChan
, cmsUInt32Number OutputChan
, const void* Table
, cmsUInt32Number dwFlags
)
165 cmsUInt32Number Samples
[MAX_INPUT_DIMENSIONS
];
167 // Fill the auxiliary array
168 for (i
=0; i
< MAX_INPUT_DIMENSIONS
; i
++)
169 Samples
[i
] = nSamples
;
171 // Call the extended function
172 return _cmsComputeInterpParamsEx(ContextID
, Samples
, InputChan
, OutputChan
, Table
, dwFlags
);
176 // Free all associated memory
177 void CMSEXPORT
_cmsFreeInterpParams(cmsInterpParams
* p
)
179 if (p
!= NULL
) _cmsFree(p
->ContextID
, p
);
183 // Inline fixed point interpolation
184 cmsINLINE CMS_NO_SANITIZE cmsUInt16Number
LinearInterp(cmsS15Fixed16Number a
, cmsS15Fixed16Number l
, cmsS15Fixed16Number h
)
186 cmsUInt32Number dif
= (cmsUInt32Number
) (h
- l
) * a
+ 0x8000;
187 dif
= (dif
>> 16) + l
;
188 return (cmsUInt16Number
) (dif
);
192 // Linear interpolation (Fixed-point optimized)
194 void LinLerp1D(CMSREGISTER
const cmsUInt16Number Value
[],
195 CMSREGISTER cmsUInt16Number Output
[],
196 CMSREGISTER
const cmsInterpParams
* p
)
198 cmsUInt16Number y1
, y0
;
201 const cmsUInt16Number
* LutTable
= (cmsUInt16Number
*) p
->Table
;
203 // if last value or just one point
204 if (Value
[0] == 0xffff || p
->Domain
[0] == 0) {
206 Output
[0] = LutTable
[p
-> Domain
[0]];
210 val3
= p
->Domain
[0] * Value
[0];
211 val3
= _cmsToFixedDomain(val3
); // To fixed 15.16
213 cell0
= FIXED_TO_INT(val3
); // Cell is 16 MSB bits
214 rest
= FIXED_REST_TO_INT(val3
); // Rest is 16 LSB bits
216 y0
= LutTable
[cell0
];
217 y1
= LutTable
[cell0
+ 1];
219 Output
[0] = LinearInterp(rest
, y0
, y1
);
223 // To prevent out of bounds indexing
224 cmsINLINE cmsFloat32Number
fclamp(cmsFloat32Number v
)
226 return ((v
< 1.0e-9f
) || isnan(v
)) ? 0.0f
: (v
> 1.0f
? 1.0f
: v
);
229 // Floating-point version of 1D interpolation
231 void LinLerp1Dfloat(const cmsFloat32Number Value
[],
232 cmsFloat32Number Output
[],
233 const cmsInterpParams
* p
)
235 cmsFloat32Number y1
, y0
;
236 cmsFloat32Number val2
, rest
;
238 const cmsFloat32Number
* LutTable
= (cmsFloat32Number
*) p
->Table
;
240 val2
= fclamp(Value
[0]);
243 if (val2
== 1.0 || p
->Domain
[0] == 0) {
244 Output
[0] = LutTable
[p
-> Domain
[0]];
248 val2
*= p
->Domain
[0];
250 cell0
= (int)floor(val2
);
251 cell1
= (int)ceil(val2
);
253 // Rest is 16 LSB bits
256 y0
= LutTable
[cell0
];
257 y1
= LutTable
[cell1
];
259 Output
[0] = y0
+ (y1
- y0
) * rest
;
265 // Eval gray LUT having only one input channel
266 static CMS_NO_SANITIZE
267 void Eval1Input(CMSREGISTER
const cmsUInt16Number Input
[],
268 CMSREGISTER cmsUInt16Number Output
[],
269 CMSREGISTER
const cmsInterpParams
* p16
)
271 cmsS15Fixed16Number fk
;
272 cmsS15Fixed16Number k0
, k1
, rk
, K0
, K1
;
274 cmsUInt32Number OutChan
;
275 const cmsUInt16Number
* LutTable
= (cmsUInt16Number
*) p16
-> Table
;
279 if (Input
[0] == 0xffff || p16
->Domain
[0] == 0) {
281 cmsUInt32Number y0
= p16
->Domain
[0] * p16
->opta
[0];
283 for (OutChan
= 0; OutChan
< p16
->nOutputs
; OutChan
++) {
284 Output
[OutChan
] = LutTable
[y0
+ OutChan
];
290 v
= Input
[0] * p16
->Domain
[0];
291 fk
= _cmsToFixedDomain(v
);
293 k0
= FIXED_TO_INT(fk
);
294 rk
= (cmsUInt16Number
)FIXED_REST_TO_INT(fk
);
296 k1
= k0
+ (Input
[0] != 0xFFFFU
? 1 : 0);
298 K0
= p16
->opta
[0] * k0
;
299 K1
= p16
->opta
[0] * k1
;
301 for (OutChan
= 0; OutChan
< p16
->nOutputs
; OutChan
++) {
303 Output
[OutChan
] = LinearInterp(rk
, LutTable
[K0
+ OutChan
], LutTable
[K1
+ OutChan
]);
310 // Eval gray LUT having only one input channel
312 void Eval1InputFloat(const cmsFloat32Number Value
[],
313 cmsFloat32Number Output
[],
314 const cmsInterpParams
* p
)
316 cmsFloat32Number y1
, y0
;
317 cmsFloat32Number val2
, rest
;
319 cmsUInt32Number OutChan
;
320 const cmsFloat32Number
* LutTable
= (cmsFloat32Number
*) p
->Table
;
322 val2
= fclamp(Value
[0]);
325 if (val2
== 1.0 || p
->Domain
[0] == 0) {
327 cmsUInt32Number start
= p
->Domain
[0] * p
->opta
[0];
329 for (OutChan
= 0; OutChan
< p
->nOutputs
; OutChan
++) {
330 Output
[OutChan
] = LutTable
[start
+ OutChan
];
335 val2
*= p
->Domain
[0];
337 cell0
= (int)floor(val2
);
338 cell1
= (int)ceil(val2
);
340 // Rest is 16 LSB bits
346 for (OutChan
= 0; OutChan
< p
->nOutputs
; OutChan
++) {
348 y0
= LutTable
[cell0
+ OutChan
];
349 y1
= LutTable
[cell1
+ OutChan
];
351 Output
[OutChan
] = y0
+ (y1
- y0
) * rest
;
356 // Bilinear interpolation (16 bits) - cmsFloat32Number version
358 void BilinearInterpFloat(const cmsFloat32Number Input
[],
359 cmsFloat32Number Output
[],
360 const cmsInterpParams
* p
)
363 # define LERP(a,l,h) (cmsFloat32Number) ((l)+(((h)-(l))*(a)))
364 # define DENS(i,j) (LutTable[(i)+(j)+OutChan])
366 const cmsFloat32Number
* LutTable
= (cmsFloat32Number
*) p
->Table
;
367 cmsFloat32Number px
, py
;
370 int TotalOut
, OutChan
;
371 cmsFloat32Number fx
, fy
,
376 TotalOut
= p
-> nOutputs
;
377 px
= fclamp(Input
[0]) * p
->Domain
[0];
378 py
= fclamp(Input
[1]) * p
->Domain
[1];
380 x0
= (int) _cmsQuickFloor(px
); fx
= px
- (cmsFloat32Number
) x0
;
381 y0
= (int) _cmsQuickFloor(py
); fy
= py
- (cmsFloat32Number
) y0
;
383 X0
= p
-> opta
[1] * x0
;
384 X1
= X0
+ (fclamp(Input
[0]) >= 1.0 ? 0 : p
->opta
[1]);
386 Y0
= p
-> opta
[0] * y0
;
387 Y1
= Y0
+ (fclamp(Input
[1]) >= 1.0 ? 0 : p
->opta
[0]);
389 for (OutChan
= 0; OutChan
< TotalOut
; OutChan
++) {
396 dx0
= LERP(fx
, d00
, d10
);
397 dx1
= LERP(fx
, d01
, d11
);
399 dxy
= LERP(fy
, dx0
, dx1
);
401 Output
[OutChan
] = dxy
;
409 // Bilinear interpolation (16 bits) - optimized version
410 static CMS_NO_SANITIZE
411 void BilinearInterp16(CMSREGISTER
const cmsUInt16Number Input
[],
412 CMSREGISTER cmsUInt16Number Output
[],
413 CMSREGISTER
const cmsInterpParams
* p
)
416 #define DENS(i,j) (LutTable[(i)+(j)+OutChan])
417 #define LERP(a,l,h) (cmsUInt16Number) (l + ROUND_FIXED_TO_INT(((h-l)*a)))
419 const cmsUInt16Number
* LutTable
= (cmsUInt16Number
*) p
->Table
;
420 int OutChan
, TotalOut
;
421 cmsS15Fixed16Number fx
, fy
;
422 CMSREGISTER
int rx
, ry
;
424 CMSREGISTER
int X0
, X1
, Y0
, Y1
;
426 int d00
, d01
, d10
, d11
,
430 TotalOut
= p
-> nOutputs
;
432 fx
= _cmsToFixedDomain((int) Input
[0] * p
-> Domain
[0]);
433 x0
= FIXED_TO_INT(fx
);
434 rx
= FIXED_REST_TO_INT(fx
); // Rest in 0..1.0 domain
437 fy
= _cmsToFixedDomain((int) Input
[1] * p
-> Domain
[1]);
438 y0
= FIXED_TO_INT(fy
);
439 ry
= FIXED_REST_TO_INT(fy
);
442 X0
= p
-> opta
[1] * x0
;
443 X1
= X0
+ (Input
[0] == 0xFFFFU
? 0 : p
->opta
[1]);
445 Y0
= p
-> opta
[0] * y0
;
446 Y1
= Y0
+ (Input
[1] == 0xFFFFU
? 0 : p
->opta
[0]);
448 for (OutChan
= 0; OutChan
< TotalOut
; OutChan
++) {
455 dx0
= LERP(rx
, d00
, d10
);
456 dx1
= LERP(rx
, d01
, d11
);
458 dxy
= LERP(ry
, dx0
, dx1
);
460 Output
[OutChan
] = (cmsUInt16Number
) dxy
;
469 // Trilinear interpolation (16 bits) - cmsFloat32Number version
471 void TrilinearInterpFloat(const cmsFloat32Number Input
[],
472 cmsFloat32Number Output
[],
473 const cmsInterpParams
* p
)
476 # define LERP(a,l,h) (cmsFloat32Number) ((l)+(((h)-(l))*(a)))
477 # define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
479 const cmsFloat32Number
* LutTable
= (cmsFloat32Number
*) p
->Table
;
480 cmsFloat32Number px
, py
, pz
;
482 X0
, Y0
, Z0
, X1
, Y1
, Z1
;
483 int TotalOut
, OutChan
;
485 cmsFloat32Number fx
, fy
, fz
,
486 d000
, d001
, d010
, d011
,
487 d100
, d101
, d110
, d111
,
488 dx00
, dx01
, dx10
, dx11
,
491 TotalOut
= p
-> nOutputs
;
493 // We need some clipping here
494 px
= fclamp(Input
[0]) * p
->Domain
[0];
495 py
= fclamp(Input
[1]) * p
->Domain
[1];
496 pz
= fclamp(Input
[2]) * p
->Domain
[2];
498 x0
= (int) floor(px
); fx
= px
- (cmsFloat32Number
) x0
; // We need full floor funcionality here
499 y0
= (int) floor(py
); fy
= py
- (cmsFloat32Number
) y0
;
500 z0
= (int) floor(pz
); fz
= pz
- (cmsFloat32Number
) z0
;
502 X0
= p
-> opta
[2] * x0
;
503 X1
= X0
+ (fclamp(Input
[0]) >= 1.0 ? 0 : p
->opta
[2]);
505 Y0
= p
-> opta
[1] * y0
;
506 Y1
= Y0
+ (fclamp(Input
[1]) >= 1.0 ? 0 : p
->opta
[1]);
508 Z0
= p
-> opta
[0] * z0
;
509 Z1
= Z0
+ (fclamp(Input
[2]) >= 1.0 ? 0 : p
->opta
[0]);
511 for (OutChan
= 0; OutChan
< TotalOut
; OutChan
++) {
513 d000
= DENS(X0
, Y0
, Z0
);
514 d001
= DENS(X0
, Y0
, Z1
);
515 d010
= DENS(X0
, Y1
, Z0
);
516 d011
= DENS(X0
, Y1
, Z1
);
518 d100
= DENS(X1
, Y0
, Z0
);
519 d101
= DENS(X1
, Y0
, Z1
);
520 d110
= DENS(X1
, Y1
, Z0
);
521 d111
= DENS(X1
, Y1
, Z1
);
524 dx00
= LERP(fx
, d000
, d100
);
525 dx01
= LERP(fx
, d001
, d101
);
526 dx10
= LERP(fx
, d010
, d110
);
527 dx11
= LERP(fx
, d011
, d111
);
529 dxy0
= LERP(fy
, dx00
, dx10
);
530 dxy1
= LERP(fy
, dx01
, dx11
);
532 dxyz
= LERP(fz
, dxy0
, dxy1
);
534 Output
[OutChan
] = dxyz
;
542 // Trilinear interpolation (16 bits) - optimized version
543 static CMS_NO_SANITIZE
544 void TrilinearInterp16(CMSREGISTER
const cmsUInt16Number Input
[],
545 CMSREGISTER cmsUInt16Number Output
[],
546 CMSREGISTER
const cmsInterpParams
* p
)
549 #define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
550 #define LERP(a,l,h) (cmsUInt16Number) (l + ROUND_FIXED_TO_INT(((h-l)*a)))
552 const cmsUInt16Number
* LutTable
= (cmsUInt16Number
*) p
->Table
;
553 int OutChan
, TotalOut
;
554 cmsS15Fixed16Number fx
, fy
, fz
;
555 CMSREGISTER
int rx
, ry
, rz
;
557 CMSREGISTER
int X0
, X1
, Y0
, Y1
, Z0
, Z1
;
558 int d000
, d001
, d010
, d011
,
559 d100
, d101
, d110
, d111
,
560 dx00
, dx01
, dx10
, dx11
,
563 TotalOut
= p
-> nOutputs
;
565 fx
= _cmsToFixedDomain((int) Input
[0] * p
-> Domain
[0]);
566 x0
= FIXED_TO_INT(fx
);
567 rx
= FIXED_REST_TO_INT(fx
); // Rest in 0..1.0 domain
570 fy
= _cmsToFixedDomain((int) Input
[1] * p
-> Domain
[1]);
571 y0
= FIXED_TO_INT(fy
);
572 ry
= FIXED_REST_TO_INT(fy
);
574 fz
= _cmsToFixedDomain((int) Input
[2] * p
-> Domain
[2]);
575 z0
= FIXED_TO_INT(fz
);
576 rz
= FIXED_REST_TO_INT(fz
);
579 X0
= p
-> opta
[2] * x0
;
580 X1
= X0
+ (Input
[0] == 0xFFFFU
? 0 : p
->opta
[2]);
582 Y0
= p
-> opta
[1] * y0
;
583 Y1
= Y0
+ (Input
[1] == 0xFFFFU
? 0 : p
->opta
[1]);
585 Z0
= p
-> opta
[0] * z0
;
586 Z1
= Z0
+ (Input
[2] == 0xFFFFU
? 0 : p
->opta
[0]);
588 for (OutChan
= 0; OutChan
< TotalOut
; OutChan
++) {
590 d000
= DENS(X0
, Y0
, Z0
);
591 d001
= DENS(X0
, Y0
, Z1
);
592 d010
= DENS(X0
, Y1
, Z0
);
593 d011
= DENS(X0
, Y1
, Z1
);
595 d100
= DENS(X1
, Y0
, Z0
);
596 d101
= DENS(X1
, Y0
, Z1
);
597 d110
= DENS(X1
, Y1
, Z0
);
598 d111
= DENS(X1
, Y1
, Z1
);
601 dx00
= LERP(rx
, d000
, d100
);
602 dx01
= LERP(rx
, d001
, d101
);
603 dx10
= LERP(rx
, d010
, d110
);
604 dx11
= LERP(rx
, d011
, d111
);
606 dxy0
= LERP(ry
, dx00
, dx10
);
607 dxy1
= LERP(ry
, dx01
, dx11
);
609 dxyz
= LERP(rz
, dxy0
, dxy1
);
611 Output
[OutChan
] = (cmsUInt16Number
) dxyz
;
620 // Tetrahedral interpolation, using Sakamoto algorithm.
621 #define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
623 void TetrahedralInterpFloat(const cmsFloat32Number Input
[],
624 cmsFloat32Number Output
[],
625 const cmsInterpParams
* p
)
627 const cmsFloat32Number
* LutTable
= (cmsFloat32Number
*) p
-> Table
;
628 cmsFloat32Number px
, py
, pz
;
630 X0
, Y0
, Z0
, X1
, Y1
, Z1
;
631 cmsFloat32Number rx
, ry
, rz
;
632 cmsFloat32Number c0
, c1
=0, c2
=0, c3
=0;
633 int OutChan
, TotalOut
;
635 TotalOut
= p
-> nOutputs
;
637 // We need some clipping here
638 px
= fclamp(Input
[0]) * p
->Domain
[0];
639 py
= fclamp(Input
[1]) * p
->Domain
[1];
640 pz
= fclamp(Input
[2]) * p
->Domain
[2];
642 x0
= (int) floor(px
); rx
= (px
- (cmsFloat32Number
) x0
); // We need full floor functionality here
643 y0
= (int) floor(py
); ry
= (py
- (cmsFloat32Number
) y0
);
644 z0
= (int) floor(pz
); rz
= (pz
- (cmsFloat32Number
) z0
);
647 X0
= p
-> opta
[2] * x0
;
648 X1
= X0
+ (fclamp(Input
[0]) >= 1.0 ? 0 : p
->opta
[2]);
650 Y0
= p
-> opta
[1] * y0
;
651 Y1
= Y0
+ (fclamp(Input
[1]) >= 1.0 ? 0 : p
->opta
[1]);
653 Z0
= p
-> opta
[0] * z0
;
654 Z1
= Z0
+ (fclamp(Input
[2]) >= 1.0 ? 0 : p
->opta
[0]);
656 for (OutChan
=0; OutChan
< TotalOut
; OutChan
++) {
658 // These are the 6 Tetrahedral
660 c0
= DENS(X0
, Y0
, Z0
);
662 if (rx
>= ry
&& ry
>= rz
) {
664 c1
= DENS(X1
, Y0
, Z0
) - c0
;
665 c2
= DENS(X1
, Y1
, Z0
) - DENS(X1
, Y0
, Z0
);
666 c3
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y1
, Z0
);
670 if (rx
>= rz
&& rz
>= ry
) {
672 c1
= DENS(X1
, Y0
, Z0
) - c0
;
673 c2
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y0
, Z1
);
674 c3
= DENS(X1
, Y0
, Z1
) - DENS(X1
, Y0
, Z0
);
678 if (rz
>= rx
&& rx
>= ry
) {
680 c1
= DENS(X1
, Y0
, Z1
) - DENS(X0
, Y0
, Z1
);
681 c2
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y0
, Z1
);
682 c3
= DENS(X0
, Y0
, Z1
) - c0
;
686 if (ry
>= rx
&& rx
>= rz
) {
688 c1
= DENS(X1
, Y1
, Z0
) - DENS(X0
, Y1
, Z0
);
689 c2
= DENS(X0
, Y1
, Z0
) - c0
;
690 c3
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y1
, Z0
);
694 if (ry
>= rz
&& rz
>= rx
) {
696 c1
= DENS(X1
, Y1
, Z1
) - DENS(X0
, Y1
, Z1
);
697 c2
= DENS(X0
, Y1
, Z0
) - c0
;
698 c3
= DENS(X0
, Y1
, Z1
) - DENS(X0
, Y1
, Z0
);
702 if (rz
>= ry
&& ry
>= rx
) {
704 c1
= DENS(X1
, Y1
, Z1
) - DENS(X0
, Y1
, Z1
);
705 c2
= DENS(X0
, Y1
, Z1
) - DENS(X0
, Y0
, Z1
);
706 c3
= DENS(X0
, Y0
, Z1
) - c0
;
713 Output
[OutChan
] = c0
+ c1
* rx
+ c2
* ry
+ c3
* rz
;
720 static CMS_NO_SANITIZE
721 void TetrahedralInterp16(CMSREGISTER
const cmsUInt16Number Input
[],
722 CMSREGISTER cmsUInt16Number Output
[],
723 CMSREGISTER
const cmsInterpParams
* p
)
725 const cmsUInt16Number
* LutTable
= (cmsUInt16Number
*) p
-> Table
;
726 cmsS15Fixed16Number fx
, fy
, fz
;
727 cmsS15Fixed16Number rx
, ry
, rz
;
729 cmsS15Fixed16Number c0
, c1
, c2
, c3
, Rest
;
730 cmsUInt32Number X0
, X1
, Y0
, Y1
, Z0
, Z1
;
731 cmsUInt32Number TotalOut
= p
-> nOutputs
;
733 fx
= _cmsToFixedDomain((int) Input
[0] * p
-> Domain
[0]);
734 fy
= _cmsToFixedDomain((int) Input
[1] * p
-> Domain
[1]);
735 fz
= _cmsToFixedDomain((int) Input
[2] * p
-> Domain
[2]);
737 x0
= FIXED_TO_INT(fx
);
738 y0
= FIXED_TO_INT(fy
);
739 z0
= FIXED_TO_INT(fz
);
741 rx
= FIXED_REST_TO_INT(fx
);
742 ry
= FIXED_REST_TO_INT(fy
);
743 rz
= FIXED_REST_TO_INT(fz
);
745 X0
= p
-> opta
[2] * x0
;
746 X1
= (Input
[0] == 0xFFFFU
? 0 : p
->opta
[2]);
748 Y0
= p
-> opta
[1] * y0
;
749 Y1
= (Input
[1] == 0xFFFFU
? 0 : p
->opta
[1]);
751 Z0
= p
-> opta
[0] * z0
;
752 Z1
= (Input
[2] == 0xFFFFU
? 0 : p
->opta
[0]);
754 LutTable
+= X0
+Y0
+Z0
;
756 // Output should be computed as x = ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest))
757 // which expands as: x = (Rest + ((Rest+0x7fff)/0xFFFF) + 0x8000)>>16
758 // This can be replaced by: t = Rest+0x8001, x = (t + (t>>16))>>16
759 // at the cost of being off by one at 7fff and 17ffe.
765 for (; TotalOut
; TotalOut
--) {
773 Rest
= c1
* rx
+ c2
* ry
+ c3
* rz
+ 0x8001;
774 *Output
++ = (cmsUInt16Number
) c0
+ ((Rest
+ (Rest
>>16))>>16);
776 } else if (rz
>= rx
) {
779 for (; TotalOut
; TotalOut
--) {
787 Rest
= c1
* rx
+ c2
* ry
+ c3
* rz
+ 0x8001;
788 *Output
++ = (cmsUInt16Number
) c0
+ ((Rest
+ (Rest
>>16))>>16);
793 for (; TotalOut
; TotalOut
--) {
801 Rest
= c1
* rx
+ c2
* ry
+ c3
* rz
+ 0x8001;
802 *Output
++ = (cmsUInt16Number
) c0
+ ((Rest
+ (Rest
>>16))>>16);
809 for (; TotalOut
; TotalOut
--) {
817 Rest
= c1
* rx
+ c2
* ry
+ c3
* rz
+ 0x8001;
818 *Output
++ = (cmsUInt16Number
) c0
+ ((Rest
+ (Rest
>>16))>>16);
820 } else if (ry
>= rz
) {
823 for (; TotalOut
; TotalOut
--) {
831 Rest
= c1
* rx
+ c2
* ry
+ c3
* rz
+ 0x8001;
832 *Output
++ = (cmsUInt16Number
) c0
+ ((Rest
+ (Rest
>>16))>>16);
837 for (; TotalOut
; TotalOut
--) {
845 Rest
= c1
* rx
+ c2
* ry
+ c3
* rz
+ 0x8001;
846 *Output
++ = (cmsUInt16Number
) c0
+ ((Rest
+ (Rest
>>16))>>16);
853 #define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])
854 static CMS_NO_SANITIZE
855 void Eval4Inputs(CMSREGISTER
const cmsUInt16Number Input
[],
856 CMSREGISTER cmsUInt16Number Output
[],
857 CMSREGISTER
const cmsInterpParams
* p16
)
859 const cmsUInt16Number
* LutTable
;
860 cmsS15Fixed16Number fk
;
861 cmsS15Fixed16Number k0
, rk
;
863 cmsS15Fixed16Number fx
, fy
, fz
;
864 cmsS15Fixed16Number rx
, ry
, rz
;
866 cmsS15Fixed16Number X0
, X1
, Y0
, Y1
, Z0
, Z1
;
868 cmsS15Fixed16Number c0
, c1
, c2
, c3
, Rest
;
869 cmsUInt32Number OutChan
;
870 cmsUInt16Number Tmp1
[MAX_STAGE_CHANNELS
], Tmp2
[MAX_STAGE_CHANNELS
];
873 fk
= _cmsToFixedDomain((int) Input
[0] * p16
-> Domain
[0]);
874 fx
= _cmsToFixedDomain((int) Input
[1] * p16
-> Domain
[1]);
875 fy
= _cmsToFixedDomain((int) Input
[2] * p16
-> Domain
[2]);
876 fz
= _cmsToFixedDomain((int) Input
[3] * p16
-> Domain
[3]);
878 k0
= FIXED_TO_INT(fk
);
879 x0
= FIXED_TO_INT(fx
);
880 y0
= FIXED_TO_INT(fy
);
881 z0
= FIXED_TO_INT(fz
);
883 rk
= FIXED_REST_TO_INT(fk
);
884 rx
= FIXED_REST_TO_INT(fx
);
885 ry
= FIXED_REST_TO_INT(fy
);
886 rz
= FIXED_REST_TO_INT(fz
);
888 K0
= p16
-> opta
[3] * k0
;
889 K1
= K0
+ (Input
[0] == 0xFFFFU
? 0 : p16
->opta
[3]);
891 X0
= p16
-> opta
[2] * x0
;
892 X1
= X0
+ (Input
[1] == 0xFFFFU
? 0 : p16
->opta
[2]);
894 Y0
= p16
-> opta
[1] * y0
;
895 Y1
= Y0
+ (Input
[2] == 0xFFFFU
? 0 : p16
->opta
[1]);
897 Z0
= p16
-> opta
[0] * z0
;
898 Z1
= Z0
+ (Input
[3] == 0xFFFFU
? 0 : p16
->opta
[0]);
900 LutTable
= (cmsUInt16Number
*) p16
-> Table
;
903 for (OutChan
=0; OutChan
< p16
-> nOutputs
; OutChan
++) {
905 c0
= DENS(X0
, Y0
, Z0
);
907 if (rx
>= ry
&& ry
>= rz
) {
909 c1
= DENS(X1
, Y0
, Z0
) - c0
;
910 c2
= DENS(X1
, Y1
, Z0
) - DENS(X1
, Y0
, Z0
);
911 c3
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y1
, Z0
);
915 if (rx
>= rz
&& rz
>= ry
) {
917 c1
= DENS(X1
, Y0
, Z0
) - c0
;
918 c2
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y0
, Z1
);
919 c3
= DENS(X1
, Y0
, Z1
) - DENS(X1
, Y0
, Z0
);
923 if (rz
>= rx
&& rx
>= ry
) {
925 c1
= DENS(X1
, Y0
, Z1
) - DENS(X0
, Y0
, Z1
);
926 c2
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y0
, Z1
);
927 c3
= DENS(X0
, Y0
, Z1
) - c0
;
931 if (ry
>= rx
&& rx
>= rz
) {
933 c1
= DENS(X1
, Y1
, Z0
) - DENS(X0
, Y1
, Z0
);
934 c2
= DENS(X0
, Y1
, Z0
) - c0
;
935 c3
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y1
, Z0
);
939 if (ry
>= rz
&& rz
>= rx
) {
941 c1
= DENS(X1
, Y1
, Z1
) - DENS(X0
, Y1
, Z1
);
942 c2
= DENS(X0
, Y1
, Z0
) - c0
;
943 c3
= DENS(X0
, Y1
, Z1
) - DENS(X0
, Y1
, Z0
);
947 if (rz
>= ry
&& ry
>= rx
) {
949 c1
= DENS(X1
, Y1
, Z1
) - DENS(X0
, Y1
, Z1
);
950 c2
= DENS(X0
, Y1
, Z1
) - DENS(X0
, Y0
, Z1
);
951 c3
= DENS(X0
, Y0
, Z1
) - c0
;
958 Rest
= c1
* rx
+ c2
* ry
+ c3
* rz
;
960 Tmp1
[OutChan
] = (cmsUInt16Number
)(c0
+ ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest
)));
964 LutTable
= (cmsUInt16Number
*) p16
-> Table
;
967 for (OutChan
=0; OutChan
< p16
-> nOutputs
; OutChan
++) {
969 c0
= DENS(X0
, Y0
, Z0
);
971 if (rx
>= ry
&& ry
>= rz
) {
973 c1
= DENS(X1
, Y0
, Z0
) - c0
;
974 c2
= DENS(X1
, Y1
, Z0
) - DENS(X1
, Y0
, Z0
);
975 c3
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y1
, Z0
);
979 if (rx
>= rz
&& rz
>= ry
) {
981 c1
= DENS(X1
, Y0
, Z0
) - c0
;
982 c2
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y0
, Z1
);
983 c3
= DENS(X1
, Y0
, Z1
) - DENS(X1
, Y0
, Z0
);
987 if (rz
>= rx
&& rx
>= ry
) {
989 c1
= DENS(X1
, Y0
, Z1
) - DENS(X0
, Y0
, Z1
);
990 c2
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y0
, Z1
);
991 c3
= DENS(X0
, Y0
, Z1
) - c0
;
995 if (ry
>= rx
&& rx
>= rz
) {
997 c1
= DENS(X1
, Y1
, Z0
) - DENS(X0
, Y1
, Z0
);
998 c2
= DENS(X0
, Y1
, Z0
) - c0
;
999 c3
= DENS(X1
, Y1
, Z1
) - DENS(X1
, Y1
, Z0
);
1003 if (ry
>= rz
&& rz
>= rx
) {
1005 c1
= DENS(X1
, Y1
, Z1
) - DENS(X0
, Y1
, Z1
);
1006 c2
= DENS(X0
, Y1
, Z0
) - c0
;
1007 c3
= DENS(X0
, Y1
, Z1
) - DENS(X0
, Y1
, Z0
);
1011 if (rz
>= ry
&& ry
>= rx
) {
1013 c1
= DENS(X1
, Y1
, Z1
) - DENS(X0
, Y1
, Z1
);
1014 c2
= DENS(X0
, Y1
, Z1
) - DENS(X0
, Y0
, Z1
);
1015 c3
= DENS(X0
, Y0
, Z1
) - c0
;
1022 Rest
= c1
* rx
+ c2
* ry
+ c3
* rz
;
1024 Tmp2
[OutChan
] = (cmsUInt16Number
) (c0
+ ROUND_FIXED_TO_INT(_cmsToFixedDomain(Rest
)));
1029 for (i
=0; i
< p16
-> nOutputs
; i
++) {
1030 Output
[i
] = LinearInterp(rk
, Tmp1
[i
], Tmp2
[i
]);
1036 // For more that 3 inputs (i.e., CMYK)
1037 // evaluate two 3-dimensional interpolations and then linearly interpolate between them.
1039 void Eval4InputsFloat(const cmsFloat32Number Input
[],
1040 cmsFloat32Number Output
[],
1041 const cmsInterpParams
* p
)
1043 const cmsFloat32Number
* LutTable
= (cmsFloat32Number
*) p
-> Table
;
1044 cmsFloat32Number rest
;
1045 cmsFloat32Number pk
;
1047 const cmsFloat32Number
* T
;
1049 cmsFloat32Number Tmp1
[MAX_STAGE_CHANNELS
], Tmp2
[MAX_STAGE_CHANNELS
];
1052 pk
= fclamp(Input
[0]) * p
->Domain
[0];
1053 k0
= _cmsQuickFloor(pk
);
1054 rest
= pk
- (cmsFloat32Number
) k0
;
1056 K0
= p
-> opta
[3] * k0
;
1057 K1
= K0
+ (fclamp(Input
[0]) >= 1.0 ? 0 : p
->opta
[3]);
1060 memmove(&p1
.Domain
[0], &p
->Domain
[1], 3*sizeof(cmsUInt32Number
));
1065 TetrahedralInterpFloat(Input
+ 1, Tmp1
, &p1
);
1069 TetrahedralInterpFloat(Input
+ 1, Tmp2
, &p1
);
1071 for (i
=0; i
< p
-> nOutputs
; i
++)
1073 cmsFloat32Number y0
= Tmp1
[i
];
1074 cmsFloat32Number y1
= Tmp2
[i
];
1076 Output
[i
] = y0
+ (y1
- y0
) * rest
;
1080 #define EVAL_FNS(N,NM) static CMS_NO_SANITIZE \
1081 void Eval##N##Inputs(CMSREGISTER const cmsUInt16Number Input[], CMSREGISTER cmsUInt16Number Output[], CMSREGISTER const cmsInterpParams* p16)\
1083 const cmsUInt16Number* LutTable = (cmsUInt16Number*) p16 -> Table;\
1084 cmsS15Fixed16Number fk;\
1085 cmsS15Fixed16Number k0, rk;\
1087 const cmsUInt16Number* T;\
1089 cmsUInt16Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\
1090 cmsInterpParams p1;\
1092 fk = _cmsToFixedDomain((cmsS15Fixed16Number) Input[0] * p16 -> Domain[0]);\
1093 k0 = FIXED_TO_INT(fk);\
1094 rk = FIXED_REST_TO_INT(fk);\
1096 K0 = p16 -> opta[NM] * k0;\
1097 K1 = p16 -> opta[NM] * (k0 + (Input[0] != 0xFFFFU ? 1 : 0));\
1100 memmove(&p1.Domain[0], &p16 ->Domain[1], NM*sizeof(cmsUInt32Number));\
1105 Eval##NM##Inputs(Input + 1, Tmp1, &p1);\
1110 Eval##NM##Inputs(Input + 1, Tmp2, &p1);\
1112 for (i=0; i < p16 -> nOutputs; i++) {\
1114 Output[i] = LinearInterp(rk, Tmp1[i], Tmp2[i]);\
1118 static void Eval##N##InputsFloat(const cmsFloat32Number Input[], \
1119 cmsFloat32Number Output[],\
1120 const cmsInterpParams * p)\
1122 const cmsFloat32Number* LutTable = (cmsFloat32Number*) p -> Table;\
1123 cmsFloat32Number rest;\
1124 cmsFloat32Number pk;\
1126 const cmsFloat32Number* T;\
1128 cmsFloat32Number Tmp1[MAX_STAGE_CHANNELS], Tmp2[MAX_STAGE_CHANNELS];\
1129 cmsInterpParams p1;\
1131 pk = fclamp(Input[0]) * p->Domain[0];\
1132 k0 = _cmsQuickFloor(pk);\
1133 rest = pk - (cmsFloat32Number) k0;\
1135 K0 = p -> opta[NM] * k0;\
1136 K1 = K0 + (fclamp(Input[0]) >= 1.0 ? 0 : p->opta[NM]);\
1139 memmove(&p1.Domain[0], &p ->Domain[1], NM*sizeof(cmsUInt32Number));\
1144 Eval##NM##InputsFloat(Input + 1, Tmp1, &p1);\
1149 Eval##NM##InputsFloat(Input + 1, Tmp2, &p1);\
1151 for (i=0; i < p -> nOutputs; i++) {\
1153 cmsFloat32Number y0 = Tmp1[i];\
1154 cmsFloat32Number y1 = Tmp2[i];\
1156 Output[i] = y0 + (y1 - y0) * rest;\
1162 * Thanks to Carles Llopis for the templating idea
1177 // The default factory
1179 cmsInterpFunction
DefaultInterpolatorsFactory(cmsUInt32Number nInputChannels
, cmsUInt32Number nOutputChannels
, cmsUInt32Number dwFlags
)
1182 cmsInterpFunction Interpolation
;
1183 cmsBool IsFloat
= (dwFlags
& CMS_LERP_FLAGS_FLOAT
);
1184 cmsBool IsTrilinear
= (dwFlags
& CMS_LERP_FLAGS_TRILINEAR
);
1186 memset(&Interpolation
, 0, sizeof(Interpolation
));
1189 if (nInputChannels
>= 4 && nOutputChannels
>= MAX_STAGE_CHANNELS
)
1190 return Interpolation
;
1192 switch (nInputChannels
) {
1194 case 1: // Gray LUT / linear
1196 if (nOutputChannels
== 1) {
1199 Interpolation
.LerpFloat
= LinLerp1Dfloat
;
1201 Interpolation
.Lerp16
= LinLerp1D
;
1207 Interpolation
.LerpFloat
= Eval1InputFloat
;
1209 Interpolation
.Lerp16
= Eval1Input
;
1215 Interpolation
.LerpFloat
= BilinearInterpFloat
;
1217 Interpolation
.Lerp16
= BilinearInterp16
;
1220 case 3: // RGB et al
1225 Interpolation
.LerpFloat
= TrilinearInterpFloat
;
1227 Interpolation
.Lerp16
= TrilinearInterp16
;
1232 Interpolation
.LerpFloat
= TetrahedralInterpFloat
;
1235 Interpolation
.Lerp16
= TetrahedralInterp16
;
1243 Interpolation
.LerpFloat
= Eval4InputsFloat
;
1245 Interpolation
.Lerp16
= Eval4Inputs
;
1250 Interpolation
.LerpFloat
= Eval5InputsFloat
;
1252 Interpolation
.Lerp16
= Eval5Inputs
;
1257 Interpolation
.LerpFloat
= Eval6InputsFloat
;
1259 Interpolation
.Lerp16
= Eval6Inputs
;
1264 Interpolation
.LerpFloat
= Eval7InputsFloat
;
1266 Interpolation
.Lerp16
= Eval7Inputs
;
1271 Interpolation
.LerpFloat
= Eval8InputsFloat
;
1273 Interpolation
.Lerp16
= Eval8Inputs
;
1278 Interpolation
.LerpFloat
= Eval9InputsFloat
;
1280 Interpolation
.Lerp16
= Eval9Inputs
;
1285 Interpolation
.LerpFloat
= Eval10InputsFloat
;
1287 Interpolation
.Lerp16
= Eval10Inputs
;
1292 Interpolation
.LerpFloat
= Eval11InputsFloat
;
1294 Interpolation
.Lerp16
= Eval11Inputs
;
1299 Interpolation
.LerpFloat
= Eval12InputsFloat
;
1301 Interpolation
.Lerp16
= Eval12Inputs
;
1306 Interpolation
.LerpFloat
= Eval13InputsFloat
;
1308 Interpolation
.Lerp16
= Eval13Inputs
;
1313 Interpolation
.LerpFloat
= Eval14InputsFloat
;
1315 Interpolation
.Lerp16
= Eval14Inputs
;
1320 Interpolation
.LerpFloat
= Eval15InputsFloat
;
1322 Interpolation
.Lerp16
= Eval15Inputs
;
1326 Interpolation
.Lerp16
= NULL
;
1329 return Interpolation
;