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 //---------------------------------------------------------------------------------
26 #include "lcms2_internal.h"
28 // Tone curves are powerful constructs that can contain curves specified in diverse ways.
29 // The curve is stored in segments, where each segment can be sampled or specified by parameters.
30 // a 16.bit simplification of the *whole* curve is kept for optimization purposes. For float operation,
31 // each segment is evaluated separately. Plug-ins may be used to define new parametric schemes,
32 // each plug-in may define up to MAX_TYPES_IN_LCMS_PLUGIN functions types. For defining a function,
33 // the plug-in should provide the type id, how many parameters each type has, and a pointer to
34 // a procedure that evaluates the function. In the case of reverse evaluation, the evaluator will
35 // be called with the type id as a negative value, and a sampled version of the reversed curve
38 // ----------------------------------------------------------------- Implementation
39 // Maxim number of nodes
40 #define MAX_NODES_IN_CURVE 4097
41 #define MINUS_INF (-1E22F)
42 #define PLUS_INF (+1E22F)
44 // The list of supported parametric curves
45 typedef struct _cmsParametricCurvesCollection_st
{
47 cmsUInt32Number nFunctions
; // Number of supported functions in this chunk
48 cmsInt32Number FunctionTypes
[MAX_TYPES_IN_LCMS_PLUGIN
]; // The identification types
49 cmsUInt32Number ParameterCount
[MAX_TYPES_IN_LCMS_PLUGIN
]; // Number of parameters for each function
51 cmsParametricCurveEvaluator Evaluator
; // The evaluator
53 struct _cmsParametricCurvesCollection_st
* Next
; // Next in list
55 } _cmsParametricCurvesCollection
;
57 // This is the default (built-in) evaluator
58 static cmsFloat64Number
DefaultEvalParametricFn(cmsInt32Number Type
, const cmsFloat64Number Params
[], cmsFloat64Number R
);
61 static _cmsParametricCurvesCollection DefaultCurves
= {
62 10, // # of curve types
63 { 1, 2, 3, 4, 5, 6, 7, 8, 108, 109 }, // Parametric curve ID
64 { 1, 3, 4, 5, 7, 4, 5, 5, 1, 1 }, // Parameters by type
65 DefaultEvalParametricFn
, // Evaluator
69 // Duplicates the zone of memory used by the plug-in in the new context
71 void DupPluginCurvesList(struct _cmsContext_struct
* ctx
,
72 const struct _cmsContext_struct
* src
)
74 _cmsCurvesPluginChunkType newHead
= { NULL
};
75 _cmsParametricCurvesCollection
* entry
;
76 _cmsParametricCurvesCollection
* Anterior
= NULL
;
77 _cmsCurvesPluginChunkType
* head
= (_cmsCurvesPluginChunkType
*) src
->chunks
[CurvesPlugin
];
79 _cmsAssert(head
!= NULL
);
81 // Walk the list copying all nodes
82 for (entry
= head
->ParametricCurves
;
84 entry
= entry
->Next
) {
86 _cmsParametricCurvesCollection
*newEntry
= ( _cmsParametricCurvesCollection
*) _cmsSubAllocDup(ctx
->MemPool
, entry
, sizeof(_cmsParametricCurvesCollection
));
91 // We want to keep the linked list order, so this is a little bit tricky
92 newEntry
-> Next
= NULL
;
94 Anterior
-> Next
= newEntry
;
98 if (newHead
.ParametricCurves
== NULL
)
99 newHead
.ParametricCurves
= newEntry
;
102 ctx
->chunks
[CurvesPlugin
] = _cmsSubAllocDup(ctx
->MemPool
, &newHead
, sizeof(_cmsCurvesPluginChunkType
));
105 // The allocator have to follow the chain
106 void _cmsAllocCurvesPluginChunk(struct _cmsContext_struct
* ctx
,
107 const struct _cmsContext_struct
* src
)
109 _cmsAssert(ctx
!= NULL
);
113 // Copy all linked list
114 DupPluginCurvesList(ctx
, src
);
117 static _cmsCurvesPluginChunkType CurvesPluginChunk
= { NULL
};
118 ctx
->chunks
[CurvesPlugin
] = _cmsSubAllocDup(ctx
->MemPool
, &CurvesPluginChunk
, sizeof(_cmsCurvesPluginChunkType
));
123 // The linked list head
124 _cmsCurvesPluginChunkType _cmsCurvesPluginChunk
= { NULL
};
126 // As a way to install new parametric curves
127 cmsBool
_cmsRegisterParametricCurvesPlugin(cmsContext ContextID
, cmsPluginBase
* Data
)
129 _cmsCurvesPluginChunkType
* ctx
= ( _cmsCurvesPluginChunkType
*) _cmsContextGetClientChunk(ContextID
, CurvesPlugin
);
130 cmsPluginParametricCurves
* Plugin
= (cmsPluginParametricCurves
*) Data
;
131 _cmsParametricCurvesCollection
* fl
;
135 ctx
-> ParametricCurves
= NULL
;
139 fl
= (_cmsParametricCurvesCollection
*) _cmsPluginMalloc(ContextID
, sizeof(_cmsParametricCurvesCollection
));
140 if (fl
== NULL
) return FALSE
;
142 // Copy the parameters
143 fl
->Evaluator
= Plugin
->Evaluator
;
144 fl
->nFunctions
= Plugin
->nFunctions
;
146 // Make sure no mem overwrites
147 if (fl
->nFunctions
> MAX_TYPES_IN_LCMS_PLUGIN
)
148 fl
->nFunctions
= MAX_TYPES_IN_LCMS_PLUGIN
;
151 memmove(fl
->FunctionTypes
, Plugin
->FunctionTypes
, fl
->nFunctions
* sizeof(cmsUInt32Number
));
152 memmove(fl
->ParameterCount
, Plugin
->ParameterCount
, fl
->nFunctions
* sizeof(cmsUInt32Number
));
155 fl
->Next
= ctx
->ParametricCurves
;
156 ctx
->ParametricCurves
= fl
;
163 // Search in type list, return position or -1 if not found
165 int IsInSet(int Type
, _cmsParametricCurvesCollection
* c
)
169 for (i
=0; i
< (int) c
->nFunctions
; i
++)
170 if (abs(Type
) == c
->FunctionTypes
[i
]) return i
;
176 // Search for the collection which contains a specific type
178 _cmsParametricCurvesCollection
*GetParametricCurveByType(cmsContext ContextID
, int Type
, int* index
)
180 _cmsParametricCurvesCollection
* c
;
182 _cmsCurvesPluginChunkType
* ctx
= ( _cmsCurvesPluginChunkType
*) _cmsContextGetClientChunk(ContextID
, CurvesPlugin
);
184 for (c
= ctx
->ParametricCurves
; c
!= NULL
; c
= c
->Next
) {
186 Position
= IsInSet(Type
, c
);
188 if (Position
!= -1) {
194 // If none found, revert for defaults
195 for (c
= &DefaultCurves
; c
!= NULL
; c
= c
->Next
) {
197 Position
= IsInSet(Type
, c
);
199 if (Position
!= -1) {
209 // Low level allocate, which takes care of memory details. nEntries may be zero, and in this case
210 // no optimization curve is computed. nSegments may also be zero in the inverse case, where only the
211 // optimization curve is given. Both features simultaneously is an error
213 cmsToneCurve
* AllocateToneCurveStruct(cmsContext ContextID
, cmsUInt32Number nEntries
,
214 cmsUInt32Number nSegments
, const cmsCurveSegment
* Segments
,
215 const cmsUInt16Number
* Values
)
220 // We allow huge tables, which are then restricted for smoothing operations
221 if (nEntries
> 65530) {
222 cmsSignalError(ContextID
, cmsERROR_RANGE
, "Couldn't create tone curve of more than 65530 entries");
226 if (nEntries
== 0 && nSegments
== 0) {
227 cmsSignalError(ContextID
, cmsERROR_RANGE
, "Couldn't create tone curve with zero segments and no table");
231 // Allocate all required pointers, etc.
232 p
= (cmsToneCurve
*) _cmsMallocZero(ContextID
, sizeof(cmsToneCurve
));
235 // In this case, there are no segments
236 if (nSegments
== 0) {
241 p
->Segments
= (cmsCurveSegment
*) _cmsCalloc(ContextID
, nSegments
, sizeof(cmsCurveSegment
));
242 if (p
->Segments
== NULL
) goto Error
;
244 p
->Evals
= (cmsParametricCurveEvaluator
*) _cmsCalloc(ContextID
, nSegments
, sizeof(cmsParametricCurveEvaluator
));
245 if (p
->Evals
== NULL
) goto Error
;
248 p
-> nSegments
= nSegments
;
250 // This 16-bit table contains a limited precision representation of the whole curve and is kept for
251 // increasing xput on certain operations.
256 p
->Table16
= (cmsUInt16Number
*) _cmsCalloc(ContextID
, nEntries
, sizeof(cmsUInt16Number
));
257 if (p
->Table16
== NULL
) goto Error
;
260 p
-> nEntries
= nEntries
;
262 // Initialize members if requested
263 if (Values
!= NULL
&& (nEntries
> 0)) {
265 for (i
=0; i
< nEntries
; i
++)
266 p
->Table16
[i
] = Values
[i
];
269 // Initialize the segments stuff. The evaluator for each segment is located and a pointer to it
270 // is placed in advance to maximize performance.
271 if (Segments
!= NULL
&& (nSegments
> 0)) {
273 _cmsParametricCurvesCollection
*c
;
275 p
->SegInterp
= (cmsInterpParams
**) _cmsCalloc(ContextID
, nSegments
, sizeof(cmsInterpParams
*));
276 if (p
->SegInterp
== NULL
) goto Error
;
278 for (i
=0; i
< nSegments
; i
++) {
280 // Type 0 is a special marker for table-based curves
281 if (Segments
[i
].Type
== 0)
282 p
->SegInterp
[i
] = _cmsComputeInterpParams(ContextID
, Segments
[i
].nGridPoints
, 1, 1, NULL
, CMS_LERP_FLAGS_FLOAT
);
284 memmove(&p
->Segments
[i
], &Segments
[i
], sizeof(cmsCurveSegment
));
286 if (Segments
[i
].Type
== 0 && Segments
[i
].SampledPoints
!= NULL
)
287 p
->Segments
[i
].SampledPoints
= (cmsFloat32Number
*) _cmsDupMem(ContextID
, Segments
[i
].SampledPoints
, sizeof(cmsFloat32Number
) * Segments
[i
].nGridPoints
);
289 p
->Segments
[i
].SampledPoints
= NULL
;
292 c
= GetParametricCurveByType(ContextID
, Segments
[i
].Type
, NULL
);
294 p
->Evals
[i
] = c
->Evaluator
;
298 p
->InterpParams
= _cmsComputeInterpParams(ContextID
, p
->nEntries
, 1, 1, p
->Table16
, CMS_LERP_FLAGS_16BITS
);
299 if (p
->InterpParams
!= NULL
)
303 for (i
=0; i
< nSegments
; i
++) {
304 if (p
->Segments
&& p
->Segments
[i
].SampledPoints
) _cmsFree(ContextID
, p
->Segments
[i
].SampledPoints
);
305 if (p
->SegInterp
&& p
->SegInterp
[i
]) _cmsFree(ContextID
, p
->SegInterp
[i
]);
307 if (p
-> SegInterp
) _cmsFree(ContextID
, p
-> SegInterp
);
308 if (p
-> Segments
) _cmsFree(ContextID
, p
-> Segments
);
309 if (p
-> Evals
) _cmsFree(ContextID
, p
-> Evals
);
310 if (p
->Table16
) _cmsFree(ContextID
, p
->Table16
);
311 _cmsFree(ContextID
, p
);
316 // Generates a sigmoidal function with desired steepness.
317 cmsINLINE
double sigmoid_base(double k
, double t
)
319 return (1.0 / (1.0 + exp(-k
* t
))) - 0.5;
322 cmsINLINE
double inverted_sigmoid_base(double k
, double t
)
324 return -log((1.0 / (t
+ 0.5)) - 1.0) / k
;
327 cmsINLINE
double sigmoid_factory(double k
, double t
)
329 double correction
= 0.5 / sigmoid_base(k
, 1);
331 return correction
* sigmoid_base(k
, 2.0 * t
- 1.0) + 0.5;
334 cmsINLINE
double inverse_sigmoid_factory(double k
, double t
)
336 double correction
= 0.5 / sigmoid_base(k
, 1);
338 return (inverted_sigmoid_base(k
, (t
- 0.5) / correction
) + 1.0) / 2.0;
342 // Parametric Fn using floating point
344 cmsFloat64Number
DefaultEvalParametricFn(cmsInt32Number Type
, const cmsFloat64Number Params
[], cmsFloat64Number R
)
346 cmsFloat64Number e
, Val
, disc
;
354 if (fabs(Params
[0] - 1.0) < MATRIX_DET_TOLERANCE
)
360 Val
= pow(R
, Params
[0]);
363 // Type 1 Reversed: X = Y ^1/gamma
367 if (fabs(Params
[0] - 1.0) < MATRIX_DET_TOLERANCE
)
374 if (fabs(Params
[0]) < MATRIX_DET_TOLERANCE
)
377 Val
= pow(R
, 1 / Params
[0]);
382 // Y = (aX + b)^Gamma | X >= -b/a
387 if (fabs(Params
[1]) < MATRIX_DET_TOLERANCE
)
393 disc
= -Params
[2] / Params
[1];
397 e
= Params
[1] * R
+ Params
[2];
400 Val
= pow(e
, Params
[0]);
411 // X = (Y ^1/g - b) / a
414 if (fabs(Params
[0]) < MATRIX_DET_TOLERANCE
||
415 fabs(Params
[1]) < MATRIX_DET_TOLERANCE
)
424 Val
= (pow(R
, 1.0 / Params
[0]) - Params
[2]) / Params
[1];
434 // Y = (aX + b)^Gamma + c | X <= -b/a
438 if (fabs(Params
[1]) < MATRIX_DET_TOLERANCE
)
444 disc
= -Params
[2] / Params
[1];
450 e
= Params
[1] * R
+ Params
[2];
453 Val
= pow(e
, Params
[0]) + Params
[3];
465 // X=((Y-c)^1/g - b)/a | (Y>=c)
469 if (fabs(Params
[0]) < MATRIX_DET_TOLERANCE
||
470 fabs(Params
[1]) < MATRIX_DET_TOLERANCE
)
476 if (R
>= Params
[3]) {
481 Val
= (pow(e
, 1 / Params
[0]) - Params
[2]) / Params
[1];
486 Val
= -Params
[2] / Params
[1];
493 // IEC 61966-2.1 (sRGB)
494 // Y = (aX + b)^Gamma | X >= d
497 if (R
>= Params
[4]) {
499 e
= Params
[1]*R
+ Params
[2];
502 Val
= pow(e
, Params
[0]);
511 // X=((Y^1/g-b)/a) | Y >= (ad+b)^g
512 // X=Y/c | Y< (ad+b)^g
516 e
= Params
[1] * Params
[4] + Params
[2];
520 disc
= pow(e
, Params
[0]);
524 if (fabs(Params
[0]) < MATRIX_DET_TOLERANCE
||
525 fabs(Params
[1]) < MATRIX_DET_TOLERANCE
)
530 Val
= (pow(R
, 1.0 / Params
[0]) - Params
[2]) / Params
[1];
534 if (fabs(Params
[3]) < MATRIX_DET_TOLERANCE
)
544 // Y = (aX + b)^Gamma + e | X >= d
545 // Y = cX + f | X < d
547 if (R
>= Params
[4]) {
549 e
= Params
[1]*R
+ Params
[2];
552 Val
= pow(e
, Params
[0]) + Params
[5];
557 Val
= R
*Params
[3] + Params
[6];
562 // X=((Y-e)1/g-b)/a | Y >=(ad+b)^g+e), cd+f
566 disc
= Params
[3] * Params
[4] + Params
[6];
574 if (fabs(Params
[0]) < MATRIX_DET_TOLERANCE
||
575 fabs(Params
[1]) < MATRIX_DET_TOLERANCE
)
579 Val
= (pow(e
, 1.0 / Params
[0]) - Params
[2]) / Params
[1];
583 if (fabs(Params
[3]) < MATRIX_DET_TOLERANCE
)
586 Val
= (R
- Params
[6]) / Params
[3];
593 // Types 6,7,8 comes from segmented curves as described in ICCSpecRevision_02_11_06_Float.pdf
594 // Type 6 is basically identical to type 5 without d
596 // Y = (a * X + b) ^ Gamma + c
598 e
= Params
[1]*R
+ Params
[2];
600 // On gamma 1.0, don't clamp
601 if (Params
[0] == 1.0) {
608 Val
= pow(e
, Params
[0]) + Params
[3];
612 // ((Y - c) ^1/Gamma - b) / a
615 if (fabs(Params
[0]) < MATRIX_DET_TOLERANCE
||
616 fabs(Params
[1]) < MATRIX_DET_TOLERANCE
)
626 Val
= (pow(e
, 1.0 / Params
[0]) - Params
[2]) / Params
[1];
632 // Y = a * log (b * X^Gamma + c) + d
635 e
= Params
[2] * pow(R
, Params
[0]) + Params
[3];
639 Val
= Params
[1]*log10(e
) + Params
[4];
642 // (Y - d) / a = log(b * X ^Gamma + c)
643 // pow(10, (Y-d) / a) = b * X ^Gamma + c
644 // pow((pow(10, (Y-d) / a) - c) / b, 1/g) = X
647 if (fabs(Params
[0]) < MATRIX_DET_TOLERANCE
||
648 fabs(Params
[1]) < MATRIX_DET_TOLERANCE
||
649 fabs(Params
[2]) < MATRIX_DET_TOLERANCE
)
655 Val
= pow((pow(10.0, (R
- Params
[4]) / Params
[1]) - Params
[3]) / Params
[2], 1.0 / Params
[0]);
661 //Y = a * b^(c*X+d) + e
663 Val
= (Params
[0] * pow(Params
[1], Params
[2] * R
+ Params
[3]) + Params
[4]);
667 // Y = (log((y-e) / a) / log(b) - d ) / c
668 // a=0, b=1, c=2, d=3, e=4,
671 disc
= R
- Params
[4];
672 if (disc
< 0) Val
= 0;
675 if (fabs(Params
[0]) < MATRIX_DET_TOLERANCE
||
676 fabs(Params
[2]) < MATRIX_DET_TOLERANCE
)
682 Val
= (log(disc
/ Params
[0]) / log(Params
[1]) - Params
[3]) / Params
[2];
688 // S-Shaped: (1 - (1-x)^1/g)^1/g
690 if (fabs(Params
[0]) < MATRIX_DET_TOLERANCE
)
693 Val
= pow(1.0 - pow(1 - R
, 1/Params
[0]), 1/Params
[0]);
696 // y = (1 - (1-x)^1/g)^1/g
697 // y^g = (1 - (1-x)^1/g)
698 // 1 - y^g = (1-x)^1/g
699 // (1 - y^g)^g = 1 - x
702 Val
= 1 - pow(1 - pow(R
, Params
[0]), Params
[0]);
707 Val
= sigmoid_factory(Params
[0], R
);
711 Val
= inverse_sigmoid_factory(Params
[0], R
);
715 // Unsupported parametric curve. Should never reach here
722 // Evaluate a segmented function for a single value. Return -Inf if no valid segment found .
723 // If fn type is 0, perform an interpolation on the table
725 cmsFloat64Number
EvalSegmentedFn(const cmsToneCurve
*g
, cmsFloat64Number R
)
728 cmsFloat32Number Out32
;
729 cmsFloat64Number Out
;
731 for (i
= (int) g
->nSegments
- 1; i
>= 0; --i
) {
734 if ((R
> g
->Segments
[i
].x0
) && (R
<= g
->Segments
[i
].x1
)) {
736 // Type == 0 means segment is sampled
737 if (g
->Segments
[i
].Type
== 0) {
739 cmsFloat32Number R1
= (cmsFloat32Number
)(R
- g
->Segments
[i
].x0
) / (g
->Segments
[i
].x1
- g
->Segments
[i
].x0
);
741 // Setup the table (TODO: clean that)
742 g
->SegInterp
[i
]->Table
= g
->Segments
[i
].SampledPoints
;
744 g
->SegInterp
[i
]->Interpolation
.LerpFloat(&R1
, &Out32
, g
->SegInterp
[i
]);
745 Out
= (cmsFloat64Number
) Out32
;
749 Out
= g
->Evals
[i
](g
->Segments
[i
].Type
, g
->Segments
[i
].Params
, R
);
767 // Access to estimated low-res table
768 cmsUInt32Number CMSEXPORT
cmsGetToneCurveEstimatedTableEntries(const cmsToneCurve
* t
)
770 _cmsAssert(t
!= NULL
);
774 const cmsUInt16Number
* CMSEXPORT
cmsGetToneCurveEstimatedTable(const cmsToneCurve
* t
)
776 _cmsAssert(t
!= NULL
);
781 // Create an empty gamma curve, by using tables. This specifies only the limited-precision part, and leaves the
782 // floating point description empty.
783 cmsToneCurve
* CMSEXPORT
cmsBuildTabulatedToneCurve16(cmsContext ContextID
, cmsUInt32Number nEntries
, const cmsUInt16Number Values
[])
785 return AllocateToneCurveStruct(ContextID
, nEntries
, 0, NULL
, Values
);
789 cmsUInt32Number
EntriesByGamma(cmsFloat64Number Gamma
)
791 if (fabs(Gamma
- 1.0) < 0.001) return 2;
796 // Create a segmented gamma, fill the table
797 cmsToneCurve
* CMSEXPORT
cmsBuildSegmentedToneCurve(cmsContext ContextID
,
798 cmsUInt32Number nSegments
, const cmsCurveSegment Segments
[])
801 cmsFloat64Number R
, Val
;
803 cmsUInt32Number nGridPoints
= 4096;
805 _cmsAssert(Segments
!= NULL
);
807 // Optimizatin for identity curves.
808 if (nSegments
== 1 && Segments
[0].Type
== 1) {
810 nGridPoints
= EntriesByGamma(Segments
[0].Params
[0]);
813 g
= AllocateToneCurveStruct(ContextID
, nGridPoints
, nSegments
, Segments
, NULL
);
814 if (g
== NULL
) return NULL
;
816 // Once we have the floating point version, we can approximate a 16 bit table of 4096 entries
817 // for performance reasons. This table would normally not be used except on 8/16 bits transforms.
818 for (i
= 0; i
< nGridPoints
; i
++) {
820 R
= (cmsFloat64Number
) i
/ (nGridPoints
-1);
822 Val
= EvalSegmentedFn(g
, R
);
824 // Round and saturate
825 g
->Table16
[i
] = _cmsQuickSaturateWord(Val
* 65535.0);
831 // Use a segmented curve to store the floating point table
832 cmsToneCurve
* CMSEXPORT
cmsBuildTabulatedToneCurveFloat(cmsContext ContextID
, cmsUInt32Number nEntries
, const cmsFloat32Number values
[])
834 cmsCurveSegment Seg
[3];
836 // Do some housekeeping
837 if (nEntries
== 0 || values
== NULL
)
840 // A segmented tone curve should have function segments in the first and last positions
841 // Initialize segmented curve part up to 0 to constant value = samples[0]
842 Seg
[0].x0
= MINUS_INF
;
846 Seg
[0].Params
[0] = 1;
847 Seg
[0].Params
[1] = 0;
848 Seg
[0].Params
[2] = 0;
849 Seg
[0].Params
[3] = values
[0];
850 Seg
[0].Params
[4] = 0;
857 Seg
[1].nGridPoints
= nEntries
;
858 Seg
[1].SampledPoints
= (cmsFloat32Number
*) values
;
860 // Final segment is constant = lastsample
862 Seg
[2].x1
= PLUS_INF
;
865 Seg
[2].Params
[0] = 1;
866 Seg
[2].Params
[1] = 0;
867 Seg
[2].Params
[2] = 0;
868 Seg
[2].Params
[3] = values
[nEntries
-1];
869 Seg
[2].Params
[4] = 0;
872 return cmsBuildSegmentedToneCurve(ContextID
, 3, Seg
);
877 // Parameters goes as: Curve, a, b, c, d, e, f
878 // Type is the ICC type +1
879 // if type is negative, then the curve is analytically inverted
880 cmsToneCurve
* CMSEXPORT
cmsBuildParametricToneCurve(cmsContext ContextID
, cmsInt32Number Type
, const cmsFloat64Number Params
[])
882 cmsCurveSegment Seg0
;
884 cmsUInt32Number size
;
885 _cmsParametricCurvesCollection
* c
= GetParametricCurveByType(ContextID
, Type
, &Pos
);
887 _cmsAssert(Params
!= NULL
);
890 cmsSignalError(ContextID
, cmsERROR_UNKNOWN_EXTENSION
, "Invalid parametric curve type %d", Type
);
894 memset(&Seg0
, 0, sizeof(Seg0
));
900 size
= c
->ParameterCount
[Pos
] * sizeof(cmsFloat64Number
);
901 memmove(Seg0
.Params
, Params
, size
);
903 return cmsBuildSegmentedToneCurve(ContextID
, 1, &Seg0
);
908 // Build a gamma table based on gamma constant
909 cmsToneCurve
* CMSEXPORT
cmsBuildGamma(cmsContext ContextID
, cmsFloat64Number Gamma
)
911 return cmsBuildParametricToneCurve(ContextID
, 1, &Gamma
);
915 // Free all memory taken by the gamma curve
916 void CMSEXPORT
cmsFreeToneCurve(cmsToneCurve
* Curve
)
918 cmsContext ContextID
;
920 if (Curve
== NULL
) return;
922 ContextID
= Curve
->InterpParams
->ContextID
;
924 _cmsFreeInterpParams(Curve
->InterpParams
);
926 if (Curve
-> Table16
)
927 _cmsFree(ContextID
, Curve
->Table16
);
929 if (Curve
->Segments
) {
933 for (i
=0; i
< Curve
->nSegments
; i
++) {
935 if (Curve
->Segments
[i
].SampledPoints
) {
936 _cmsFree(ContextID
, Curve
->Segments
[i
].SampledPoints
);
939 if (Curve
->SegInterp
[i
] != 0)
940 _cmsFreeInterpParams(Curve
->SegInterp
[i
]);
943 _cmsFree(ContextID
, Curve
->Segments
);
944 _cmsFree(ContextID
, Curve
->SegInterp
);
948 _cmsFree(ContextID
, Curve
-> Evals
);
950 _cmsFree(ContextID
, Curve
);
953 // Utility function, free 3 gamma tables
954 void CMSEXPORT
cmsFreeToneCurveTriple(cmsToneCurve
* Curve
[3])
957 _cmsAssert(Curve
!= NULL
);
959 if (Curve
[0] != NULL
) cmsFreeToneCurve(Curve
[0]);
960 if (Curve
[1] != NULL
) cmsFreeToneCurve(Curve
[1]);
961 if (Curve
[2] != NULL
) cmsFreeToneCurve(Curve
[2]);
963 Curve
[0] = Curve
[1] = Curve
[2] = NULL
;
967 // Duplicate a gamma table
968 cmsToneCurve
* CMSEXPORT
cmsDupToneCurve(const cmsToneCurve
* In
)
970 if (In
== NULL
) return NULL
;
972 return AllocateToneCurveStruct(In
->InterpParams
->ContextID
, In
->nEntries
, In
->nSegments
, In
->Segments
, In
->Table16
);
975 // Joins two curves for X and Y. Curves should be monotonic.
980 cmsToneCurve
* CMSEXPORT
cmsJoinToneCurve(cmsContext ContextID
,
981 const cmsToneCurve
* X
,
982 const cmsToneCurve
* Y
, cmsUInt32Number nResultingPoints
)
984 cmsToneCurve
* out
= NULL
;
985 cmsToneCurve
* Yreversed
= NULL
;
986 cmsFloat32Number t
, x
;
987 cmsFloat32Number
* Res
= NULL
;
991 _cmsAssert(X
!= NULL
);
992 _cmsAssert(Y
!= NULL
);
994 Yreversed
= cmsReverseToneCurveEx(nResultingPoints
, Y
);
995 if (Yreversed
== NULL
) goto Error
;
997 Res
= (cmsFloat32Number
*) _cmsCalloc(ContextID
, nResultingPoints
, sizeof(cmsFloat32Number
));
998 if (Res
== NULL
) goto Error
;
1001 for (i
=0; i
< nResultingPoints
; i
++) {
1003 t
= (cmsFloat32Number
) i
/ (cmsFloat32Number
)(nResultingPoints
-1);
1004 x
= cmsEvalToneCurveFloat(X
, t
);
1005 Res
[i
] = cmsEvalToneCurveFloat(Yreversed
, x
);
1008 // Allocate space for output
1009 out
= cmsBuildTabulatedToneCurveFloat(ContextID
, nResultingPoints
, Res
);
1013 if (Res
!= NULL
) _cmsFree(ContextID
, Res
);
1014 if (Yreversed
!= NULL
) cmsFreeToneCurve(Yreversed
);
1021 // Get the surrounding nodes. This is tricky on non-monotonic tables
1023 int GetInterval(cmsFloat64Number In
, const cmsUInt16Number LutTable
[], const struct _cms_interp_struc
* p
)
1028 // A 1 point table is not allowed
1029 if (p
-> Domain
[0] < 1) return -1;
1031 // Let's see if ascending or descending.
1032 if (LutTable
[0] < LutTable
[p
->Domain
[0]]) {
1034 // Table is overall ascending
1035 for (i
= (int) p
->Domain
[0] - 1; i
>= 0; --i
) {
1040 if (y0
<= y1
) { // Increasing
1041 if (In
>= y0
&& In
<= y1
) return i
;
1044 if (y1
< y0
) { // Decreasing
1045 if (In
>= y1
&& In
<= y0
) return i
;
1050 // Table is overall descending
1051 for (i
=0; i
< (int) p
-> Domain
[0]; i
++) {
1056 if (y0
<= y1
) { // Increasing
1057 if (In
>= y0
&& In
<= y1
) return i
;
1060 if (y1
< y0
) { // Decreasing
1061 if (In
>= y1
&& In
<= y0
) return i
;
1069 // Reverse a gamma table
1070 cmsToneCurve
* CMSEXPORT
cmsReverseToneCurveEx(cmsUInt32Number nResultSamples
, const cmsToneCurve
* InCurve
)
1073 cmsFloat64Number a
= 0, b
= 0, y
, x1
, y1
, x2
, y2
;
1077 _cmsAssert(InCurve
!= NULL
);
1079 // Try to reverse it analytically whatever possible
1081 if (InCurve
->nSegments
== 1 && InCurve
->Segments
[0].Type
> 0 &&
1082 /* InCurve -> Segments[0].Type <= 5 */
1083 GetParametricCurveByType(InCurve
->InterpParams
->ContextID
, InCurve
->Segments
[0].Type
, NULL
) != NULL
) {
1085 return cmsBuildParametricToneCurve(InCurve
->InterpParams
->ContextID
,
1086 -(InCurve
-> Segments
[0].Type
),
1087 InCurve
-> Segments
[0].Params
);
1090 // Nope, reverse the table.
1091 out
= cmsBuildTabulatedToneCurve16(InCurve
->InterpParams
->ContextID
, nResultSamples
, NULL
);
1095 // We want to know if this is an ascending or descending table
1096 Ascending
= !cmsIsToneCurveDescending(InCurve
);
1098 // Iterate across Y axis
1099 for (i
=0; i
< (int) nResultSamples
; i
++) {
1101 y
= (cmsFloat64Number
) i
* 65535.0 / (nResultSamples
- 1);
1103 // Find interval in which y is within.
1104 j
= GetInterval(y
, InCurve
->Table16
, InCurve
->InterpParams
);
1108 // Get limits of interval
1109 x1
= InCurve
->Table16
[j
];
1110 x2
= InCurve
->Table16
[j
+1];
1112 y1
= (cmsFloat64Number
) (j
* 65535.0) / (InCurve
->nEntries
- 1);
1113 y2
= (cmsFloat64Number
) ((j
+1) * 65535.0 ) / (InCurve
->nEntries
- 1);
1115 // If collapsed, then use any
1118 out
->Table16
[i
] = _cmsQuickSaturateWord(Ascending
? y2
: y1
);
1124 a
= (y2
- y1
) / (x2
- x1
);
1129 out
->Table16
[i
] = _cmsQuickSaturateWord(a
* y
+ b
);
1136 // Reverse a gamma table
1137 cmsToneCurve
* CMSEXPORT
cmsReverseToneCurve(const cmsToneCurve
* InGamma
)
1139 _cmsAssert(InGamma
!= NULL
);
1141 return cmsReverseToneCurveEx(4096, InGamma
);
1144 // From: Eilers, P.H.C. (1994) Smoothing and interpolation with finite
1145 // differences. in: Graphic Gems IV, Heckbert, P.S. (ed.), Academic press.
1147 // Smoothing and interpolation with second differences.
1149 // Input: weights (w), data (y): vector from 1 to m.
1150 // Input: smoothing parameter (lambda), length (m).
1151 // Output: smoothed vector (z): vector from 1 to m.
1154 cmsBool
smooth2(cmsContext ContextID
, cmsFloat32Number w
[], cmsFloat32Number y
[],
1155 cmsFloat32Number z
[], cmsFloat32Number lambda
, int m
)
1158 cmsFloat32Number
*c
, *d
, *e
;
1162 c
= (cmsFloat32Number
*) _cmsCalloc(ContextID
, MAX_NODES_IN_CURVE
, sizeof(cmsFloat32Number
));
1163 d
= (cmsFloat32Number
*) _cmsCalloc(ContextID
, MAX_NODES_IN_CURVE
, sizeof(cmsFloat32Number
));
1164 e
= (cmsFloat32Number
*) _cmsCalloc(ContextID
, MAX_NODES_IN_CURVE
, sizeof(cmsFloat32Number
));
1166 if (c
!= NULL
&& d
!= NULL
&& e
!= NULL
) {
1169 d
[1] = w
[1] + lambda
;
1170 c
[1] = -2 * lambda
/ d
[1];
1171 e
[1] = lambda
/d
[1];
1173 d
[2] = w
[2] + 5 * lambda
- d
[1] * c
[1] * c
[1];
1174 c
[2] = (-4 * lambda
- d
[1] * c
[1] * e
[1]) / d
[2];
1175 e
[2] = lambda
/ d
[2];
1176 z
[2] = w
[2] * y
[2] - c
[1] * z
[1];
1178 for (i
= 3; i
< m
- 1; i
++) {
1179 i1
= i
- 1; i2
= i
- 2;
1180 d
[i
]= w
[i
] + 6 * lambda
- c
[i1
] * c
[i1
] * d
[i1
] - e
[i2
] * e
[i2
] * d
[i2
];
1181 c
[i
] = (-4 * lambda
-d
[i1
] * c
[i1
] * e
[i1
])/ d
[i
];
1182 e
[i
] = lambda
/ d
[i
];
1183 z
[i
] = w
[i
] * y
[i
] - c
[i1
] * z
[i1
] - e
[i2
] * z
[i2
];
1186 i1
= m
- 2; i2
= m
- 3;
1188 d
[m
- 1] = w
[m
- 1] + 5 * lambda
-c
[i1
] * c
[i1
] * d
[i1
] - e
[i2
] * e
[i2
] * d
[i2
];
1189 c
[m
- 1] = (-2 * lambda
- d
[i1
] * c
[i1
] * e
[i1
]) / d
[m
- 1];
1190 z
[m
- 1] = w
[m
- 1] * y
[m
- 1] - c
[i1
] * z
[i1
] - e
[i2
] * z
[i2
];
1191 i1
= m
- 1; i2
= m
- 2;
1193 d
[m
] = w
[m
] + lambda
- c
[i1
] * c
[i1
] * d
[i1
] - e
[i2
] * e
[i2
] * d
[i2
];
1194 z
[m
] = (w
[m
] * y
[m
] - c
[i1
] * z
[i1
] - e
[i2
] * z
[i2
]) / d
[m
];
1195 z
[m
- 1] = z
[m
- 1] / d
[m
- 1] - c
[m
- 1] * z
[m
];
1197 for (i
= m
- 2; 1<= i
; i
--)
1198 z
[i
] = z
[i
] / d
[i
] - c
[i
] * z
[i
+ 1] - e
[i
] * z
[i
+ 2];
1204 if (c
!= NULL
) _cmsFree(ContextID
, c
);
1205 if (d
!= NULL
) _cmsFree(ContextID
, d
);
1206 if (e
!= NULL
) _cmsFree(ContextID
, e
);
1211 // Smooths a curve sampled at regular intervals.
1212 cmsBool CMSEXPORT
cmsSmoothToneCurve(cmsToneCurve
* Tab
, cmsFloat64Number lambda
)
1214 cmsBool SuccessStatus
= TRUE
;
1215 cmsFloat32Number
*w
, *y
, *z
;
1216 cmsUInt32Number i
, nItems
, Zeros
, Poles
;
1217 cmsBool notCheck
= FALSE
;
1219 if (Tab
!= NULL
&& Tab
->InterpParams
!= NULL
)
1221 cmsContext ContextID
= Tab
->InterpParams
->ContextID
;
1223 if (!cmsIsToneCurveLinear(Tab
)) // Only non-linear curves need smoothing
1225 nItems
= Tab
->nEntries
;
1226 if (nItems
< MAX_NODES_IN_CURVE
)
1228 // Allocate one more item than needed
1229 w
= (cmsFloat32Number
*)_cmsCalloc(ContextID
, nItems
+ 1, sizeof(cmsFloat32Number
));
1230 y
= (cmsFloat32Number
*)_cmsCalloc(ContextID
, nItems
+ 1, sizeof(cmsFloat32Number
));
1231 z
= (cmsFloat32Number
*)_cmsCalloc(ContextID
, nItems
+ 1, sizeof(cmsFloat32Number
));
1233 if (w
!= NULL
&& y
!= NULL
&& z
!= NULL
) // Ensure no memory allocation failure
1235 memset(w
, 0, (nItems
+ 1) * sizeof(cmsFloat32Number
));
1236 memset(y
, 0, (nItems
+ 1) * sizeof(cmsFloat32Number
));
1237 memset(z
, 0, (nItems
+ 1) * sizeof(cmsFloat32Number
));
1239 for (i
= 0; i
< nItems
; i
++)
1241 y
[i
+ 1] = (cmsFloat32Number
)Tab
->Table16
[i
];
1251 if (smooth2(ContextID
, w
, y
, z
, (cmsFloat32Number
)lambda
, (int)nItems
))
1253 // Do some reality - checking...
1256 for (i
= nItems
; i
> 1; --i
)
1258 if (z
[i
] == 0.) Zeros
++;
1259 if (z
[i
] >= 65535.) Poles
++;
1260 if (z
[i
] < z
[i
- 1])
1262 cmsSignalError(ContextID
, cmsERROR_RANGE
, "cmsSmoothToneCurve: Non-Monotonic.");
1263 SuccessStatus
= notCheck
;
1268 if (SuccessStatus
&& Zeros
> (nItems
/ 3))
1270 cmsSignalError(ContextID
, cmsERROR_RANGE
, "cmsSmoothToneCurve: Degenerated, mostly zeros.");
1271 SuccessStatus
= notCheck
;
1274 if (SuccessStatus
&& Poles
> (nItems
/ 3))
1276 cmsSignalError(ContextID
, cmsERROR_RANGE
, "cmsSmoothToneCurve: Degenerated, mostly poles.");
1277 SuccessStatus
= notCheck
;
1280 if (SuccessStatus
) // Seems ok
1282 for (i
= 0; i
< nItems
; i
++)
1284 // Clamp to cmsUInt16Number
1285 Tab
->Table16
[i
] = _cmsQuickSaturateWord(z
[i
+ 1]);
1289 else // Could not smooth
1291 cmsSignalError(ContextID
, cmsERROR_RANGE
, "cmsSmoothToneCurve: Function smooth2 failed.");
1292 SuccessStatus
= FALSE
;
1295 else // One or more buffers could not be allocated
1297 cmsSignalError(ContextID
, cmsERROR_RANGE
, "cmsSmoothToneCurve: Could not allocate memory.");
1298 SuccessStatus
= FALSE
;
1302 _cmsFree(ContextID
, z
);
1305 _cmsFree(ContextID
, y
);
1308 _cmsFree(ContextID
, w
);
1310 else // too many items in the table
1312 cmsSignalError(ContextID
, cmsERROR_RANGE
, "cmsSmoothToneCurve: Too many points.");
1313 SuccessStatus
= FALSE
;
1317 else // Tab parameter or Tab->InterpParams is NULL
1319 // Can't signal an error here since the ContextID is not known at this point
1320 SuccessStatus
= FALSE
;
1323 return SuccessStatus
;
1326 // Is a table linear? Do not use parametric since we cannot guarantee some weird parameters resulting
1327 // in a linear table. This way assures it is linear in 12 bits, which should be enough in most cases.
1328 cmsBool CMSEXPORT
cmsIsToneCurveLinear(const cmsToneCurve
* Curve
)
1333 _cmsAssert(Curve
!= NULL
);
1335 for (i
=0; i
< (int) Curve
->nEntries
; i
++) {
1337 diff
= abs((int) Curve
->Table16
[i
] - (int) _cmsQuantizeVal(i
, Curve
->nEntries
));
1345 // Same, but for monotonicity
1346 cmsBool CMSEXPORT
cmsIsToneCurveMonotonic(const cmsToneCurve
* t
)
1350 cmsBool lDescending
;
1352 _cmsAssert(t
!= NULL
);
1354 // Degenerated curves are monotonic? Ok, let's pass them
1356 if (n
< 2) return TRUE
;
1359 lDescending
= cmsIsToneCurveDescending(t
);
1363 last
= t
->Table16
[0];
1365 for (i
= 1; i
< (int) n
; i
++) {
1367 if (t
->Table16
[i
] - last
> 2) // We allow some ripple
1370 last
= t
->Table16
[i
];
1376 last
= t
->Table16
[n
-1];
1378 for (i
= (int) n
- 2; i
>= 0; --i
) {
1380 if (t
->Table16
[i
] - last
> 2)
1383 last
= t
->Table16
[i
];
1391 // Same, but for descending tables
1392 cmsBool CMSEXPORT
cmsIsToneCurveDescending(const cmsToneCurve
* t
)
1394 _cmsAssert(t
!= NULL
);
1396 return t
->Table16
[0] > t
->Table16
[t
->nEntries
-1];
1400 // Another info fn: is out gamma table multisegment?
1401 cmsBool CMSEXPORT
cmsIsToneCurveMultisegment(const cmsToneCurve
* t
)
1403 _cmsAssert(t
!= NULL
);
1405 return t
-> nSegments
> 1;
1408 cmsInt32Number CMSEXPORT
cmsGetToneCurveParametricType(const cmsToneCurve
* t
)
1410 _cmsAssert(t
!= NULL
);
1412 if (t
-> nSegments
!= 1) return 0;
1413 return t
->Segments
[0].Type
;
1416 // We need accuracy this time
1417 cmsFloat32Number CMSEXPORT
cmsEvalToneCurveFloat(const cmsToneCurve
* Curve
, cmsFloat32Number v
)
1419 _cmsAssert(Curve
!= NULL
);
1421 // Check for 16 bits table. If so, this is a limited-precision tone curve
1422 if (Curve
->nSegments
== 0) {
1424 cmsUInt16Number In
, Out
;
1426 In
= (cmsUInt16Number
) _cmsQuickSaturateWord(v
* 65535.0);
1427 Out
= cmsEvalToneCurve16(Curve
, In
);
1429 return (cmsFloat32Number
) (Out
/ 65535.0);
1432 return (cmsFloat32Number
) EvalSegmentedFn(Curve
, v
);
1435 // We need xput over here
1436 cmsUInt16Number CMSEXPORT
cmsEvalToneCurve16(const cmsToneCurve
* Curve
, cmsUInt16Number v
)
1438 cmsUInt16Number out
;
1440 _cmsAssert(Curve
!= NULL
);
1442 Curve
->InterpParams
->Interpolation
.Lerp16(&v
, &out
, Curve
->InterpParams
);
1447 // Least squares fitting.
1448 // A mathematical procedure for finding the best-fitting curve to a given set of points by
1449 // minimizing the sum of the squares of the offsets ("the residuals") of the points from the curve.
1450 // The sum of the squares of the offsets is used instead of the offset absolute values because
1451 // this allows the residuals to be treated as a continuous differentiable quantity.
1455 // R = (yi - (xi^g))
1456 // R2 = (yi - (xi^g))2
1457 // SUM R2 = SUM (yi - (xi^g))2
1459 // dR2/dg = -2 SUM x^g log(x)(y - x^g)
1460 // solving for dR2/dg = 0
1462 // g = 1/n * SUM(log(y) / log(x))
1464 cmsFloat64Number CMSEXPORT
cmsEstimateGamma(const cmsToneCurve
* t
, cmsFloat64Number Precision
)
1466 cmsFloat64Number gamma
, sum
, sum2
;
1467 cmsFloat64Number n
, x
, y
, Std
;
1470 _cmsAssert(t
!= NULL
);
1474 // Excluding endpoints
1475 for (i
=1; i
< (MAX_NODES_IN_CURVE
-1); i
++) {
1477 x
= (cmsFloat64Number
) i
/ (MAX_NODES_IN_CURVE
-1);
1478 y
= (cmsFloat64Number
) cmsEvalToneCurveFloat(t
, (cmsFloat32Number
) x
);
1480 // Avoid 7% on lower part to prevent
1481 // artifacts due to linear ramps
1483 if (y
> 0. && y
< 1. && x
> 0.07) {
1485 gamma
= log(y
) / log(x
);
1487 sum2
+= gamma
* gamma
;
1492 // We need enough valid samples
1493 if (n
<= 1) return -1.0;
1495 // Take a look on SD to see if gamma isn't exponential at all
1496 Std
= sqrt((n
* sum2
- sum
* sum
) / (n
*(n
-1)));
1498 if (Std
> Precision
)
1501 return (sum
/ n
); // The mean
1504 // Retrieve segments on tone curves
1506 const cmsCurveSegment
* CMSEXPORT
cmsGetToneCurveSegment(cmsInt32Number n
, const cmsToneCurve
* t
)
1508 _cmsAssert(t
!= NULL
);
1510 if (n
< 0 || n
>= (cmsInt32Number
) t
->nSegments
) return NULL
;
1511 return t
->Segments
+ n
;