diff options
Diffstat (limited to 'third_party/lcms2-2.6/src/cmsgamma.c')
| -rw-r--r-- | third_party/lcms2-2.6/src/cmsgamma.c | 1298 |
1 files changed, 1298 insertions, 0 deletions
diff --git a/third_party/lcms2-2.6/src/cmsgamma.c b/third_party/lcms2-2.6/src/cmsgamma.c new file mode 100644 index 0000000000..97aeb7cc16 --- /dev/null +++ b/third_party/lcms2-2.6/src/cmsgamma.c @@ -0,0 +1,1298 @@ +//--------------------------------------------------------------------------------- +// +// Little Color Management System +// Copyright (c) 1998-2013 Marti Maria Saguer +// +// Permission is hereby granted, free of charge, to any person obtaining +// a copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the Software +// is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO +// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +// +//--------------------------------------------------------------------------------- +// + +#include "lcms2_internal.h" + +// Tone curves are powerful constructs that can contain curves specified in diverse ways. +// The curve is stored in segments, where each segment can be sampled or specified by parameters. +// a 16.bit simplification of the *whole* curve is kept for optimization purposes. For float operation, +// each segment is evaluated separately. Plug-ins may be used to define new parametric schemes, +// each plug-in may define up to MAX_TYPES_IN_LCMS_PLUGIN functions types. For defining a function, +// the plug-in should provide the type id, how many parameters each type has, and a pointer to +// a procedure that evaluates the function. In the case of reverse evaluation, the evaluator will +// be called with the type id as a negative value, and a sampled version of the reversed curve +// will be built. + +// ----------------------------------------------------------------- Implementation +// Maxim number of nodes +#define MAX_NODES_IN_CURVE 4097 +#define MINUS_INF (-1E22F) +#define PLUS_INF (+1E22F) + +// The list of supported parametric curves +typedef struct _cmsParametricCurvesCollection_st { + + int nFunctions; // Number of supported functions in this chunk + int FunctionTypes[MAX_TYPES_IN_LCMS_PLUGIN]; // The identification types + int ParameterCount[MAX_TYPES_IN_LCMS_PLUGIN]; // Number of parameters for each function + cmsParametricCurveEvaluator Evaluator; // The evaluator + + struct _cmsParametricCurvesCollection_st* Next; // Next in list + +} _cmsParametricCurvesCollection; + +// This is the default (built-in) evaluator +static cmsFloat64Number DefaultEvalParametricFn(cmsInt32Number Type, const cmsFloat64Number Params[], cmsFloat64Number R); + +// The built-in list +static _cmsParametricCurvesCollection DefaultCurves = { + 9, // # of curve types + { 1, 2, 3, 4, 5, 6, 7, 8, 108 }, // Parametric curve ID + { 1, 3, 4, 5, 7, 4, 5, 5, 1 }, // Parameters by type + DefaultEvalParametricFn, // Evaluator + NULL // Next in chain +}; + +// Duplicates the zone of memory used by the plug-in in the new context +static +void DupPluginCurvesList(struct _cmsContext_struct* ctx, + const struct _cmsContext_struct* src) +{ + _cmsCurvesPluginChunkType newHead = { NULL }; + _cmsParametricCurvesCollection* entry; + _cmsParametricCurvesCollection* Anterior = NULL; + _cmsCurvesPluginChunkType* head = (_cmsCurvesPluginChunkType*) src->chunks[CurvesPlugin]; + + _cmsAssert(head != NULL); + + // Walk the list copying all nodes + for (entry = head->ParametricCurves; + entry != NULL; + entry = entry ->Next) { + + _cmsParametricCurvesCollection *newEntry = ( _cmsParametricCurvesCollection *) _cmsSubAllocDup(ctx ->MemPool, entry, sizeof(_cmsParametricCurvesCollection)); + + if (newEntry == NULL) + return; + + // We want to keep the linked list order, so this is a little bit tricky + newEntry -> Next = NULL; + if (Anterior) + Anterior -> Next = newEntry; + + Anterior = newEntry; + + if (newHead.ParametricCurves == NULL) + newHead.ParametricCurves = newEntry; + } + + ctx ->chunks[CurvesPlugin] = _cmsSubAllocDup(ctx->MemPool, &newHead, sizeof(_cmsCurvesPluginChunkType)); +} + +// The allocator have to follow the chain +void _cmsAllocCurvesPluginChunk(struct _cmsContext_struct* ctx, + const struct _cmsContext_struct* src) +{ + _cmsAssert(ctx != NULL); + + if (src != NULL) { + + // Copy all linked list + DupPluginCurvesList(ctx, src); + } + else { + static _cmsCurvesPluginChunkType CurvesPluginChunk = { NULL }; + ctx ->chunks[CurvesPlugin] = _cmsSubAllocDup(ctx ->MemPool, &CurvesPluginChunk, sizeof(_cmsCurvesPluginChunkType)); + } +} + + +// The linked list head +_cmsCurvesPluginChunkType _cmsCurvesPluginChunk = { NULL }; + +// As a way to install new parametric curves +cmsBool _cmsRegisterParametricCurvesPlugin(cmsContext ContextID, cmsPluginBase* Data) +{ + _cmsCurvesPluginChunkType* ctx = ( _cmsCurvesPluginChunkType*) _cmsContextGetClientChunk(ContextID, CurvesPlugin); + cmsPluginParametricCurves* Plugin = (cmsPluginParametricCurves*) Data; + _cmsParametricCurvesCollection* fl; + + if (Data == NULL) { + + ctx -> ParametricCurves = NULL; + return TRUE; + } + + fl = (_cmsParametricCurvesCollection*) _cmsPluginMalloc(ContextID, sizeof(_cmsParametricCurvesCollection)); + if (fl == NULL) return FALSE; + + // Copy the parameters + fl ->Evaluator = Plugin ->Evaluator; + fl ->nFunctions = Plugin ->nFunctions; + + // Make sure no mem overwrites + if (fl ->nFunctions > MAX_TYPES_IN_LCMS_PLUGIN) + fl ->nFunctions = MAX_TYPES_IN_LCMS_PLUGIN; + + // Copy the data + memmove(fl->FunctionTypes, Plugin ->FunctionTypes, fl->nFunctions * sizeof(cmsUInt32Number)); + memmove(fl->ParameterCount, Plugin ->ParameterCount, fl->nFunctions * sizeof(cmsUInt32Number)); + + // Keep linked list + fl ->Next = ctx->ParametricCurves; + ctx->ParametricCurves = fl; + + // All is ok + return TRUE; +} + + +// Search in type list, return position or -1 if not found +static +int IsInSet(int Type, _cmsParametricCurvesCollection* c) +{ + int i; + + for (i=0; i < c ->nFunctions; i++) + if (abs(Type) == c ->FunctionTypes[i]) return i; + + return -1; +} + + +// Search for the collection which contains a specific type +static +_cmsParametricCurvesCollection *GetParametricCurveByType(cmsContext ContextID, int Type, int* index) +{ + _cmsParametricCurvesCollection* c; + int Position; + _cmsCurvesPluginChunkType* ctx = ( _cmsCurvesPluginChunkType*) _cmsContextGetClientChunk(ContextID, CurvesPlugin); + + for (c = ctx->ParametricCurves; c != NULL; c = c ->Next) { + + Position = IsInSet(Type, c); + + if (Position != -1) { + if (index != NULL) + *index = Position; + return c; + } + } + // If none found, revert for defaults + for (c = &DefaultCurves; c != NULL; c = c ->Next) { + + Position = IsInSet(Type, c); + + if (Position != -1) { + if (index != NULL) + *index = Position; + return c; + } + } + + return NULL; +} + +// Low level allocate, which takes care of memory details. nEntries may be zero, and in this case +// no optimation curve is computed. nSegments may also be zero in the inverse case, where only the +// optimization curve is given. Both features simultaneously is an error +static +cmsToneCurve* AllocateToneCurveStruct(cmsContext ContextID, cmsInt32Number nEntries, + cmsInt32Number nSegments, const cmsCurveSegment* Segments, + const cmsUInt16Number* Values) +{ + cmsToneCurve* p; + int i; + + // We allow huge tables, which are then restricted for smoothing operations + if (nEntries > 65530 || nEntries < 0) { + cmsSignalError(ContextID, cmsERROR_RANGE, "Couldn't create tone curve of more than 65530 entries"); + return NULL; + } + + if (nEntries <= 0 && nSegments <= 0) { + cmsSignalError(ContextID, cmsERROR_RANGE, "Couldn't create tone curve with zero segments and no table"); + return NULL; + } + + // Allocate all required pointers, etc. + p = (cmsToneCurve*) _cmsMallocZero(ContextID, sizeof(cmsToneCurve)); + if (!p) return NULL; + + // In this case, there are no segments + if (nSegments <= 0) { + p ->Segments = NULL; + p ->Evals = NULL; + } + else { + p ->Segments = (cmsCurveSegment*) _cmsCalloc(ContextID, nSegments, sizeof(cmsCurveSegment)); + if (p ->Segments == NULL) goto Error; + + p ->Evals = (cmsParametricCurveEvaluator*) _cmsCalloc(ContextID, nSegments, sizeof(cmsParametricCurveEvaluator)); + if (p ->Evals == NULL) goto Error; + } + + p -> nSegments = nSegments; + + // This 16-bit table contains a limited precision representation of the whole curve and is kept for + // increasing xput on certain operations. + if (nEntries <= 0) { + p ->Table16 = NULL; + } + else { + p ->Table16 = (cmsUInt16Number*) _cmsCalloc(ContextID, nEntries, sizeof(cmsUInt16Number)); + if (p ->Table16 == NULL) goto Error; + } + + p -> nEntries = nEntries; + + // Initialize members if requested + if (Values != NULL && (nEntries > 0)) { + + for (i=0; i < nEntries; i++) + p ->Table16[i] = Values[i]; + } + + // Initialize the segments stuff. The evaluator for each segment is located and a pointer to it + // is placed in advance to maximize performance. + if (Segments != NULL && (nSegments > 0)) { + + _cmsParametricCurvesCollection *c; + + p ->SegInterp = (cmsInterpParams**) _cmsCalloc(ContextID, nSegments, sizeof(cmsInterpParams*)); + if (p ->SegInterp == NULL) goto Error; + + for (i=0; i< nSegments; i++) { + + // Type 0 is a special marker for table-based curves + if (Segments[i].Type == 0) + p ->SegInterp[i] = _cmsComputeInterpParams(ContextID, Segments[i].nGridPoints, 1, 1, NULL, CMS_LERP_FLAGS_FLOAT); + + memmove(&p ->Segments[i], &Segments[i], sizeof(cmsCurveSegment)); + + if (Segments[i].Type == 0 && Segments[i].SampledPoints != NULL) + p ->Segments[i].SampledPoints = (cmsFloat32Number*) _cmsDupMem(ContextID, Segments[i].SampledPoints, sizeof(cmsFloat32Number) * Segments[i].nGridPoints); + else + p ->Segments[i].SampledPoints = NULL; + + + c = GetParametricCurveByType(ContextID, Segments[i].Type, NULL); + if (c != NULL) + p ->Evals[i] = c ->Evaluator; + } + } + + p ->InterpParams = _cmsComputeInterpParams(ContextID, p ->nEntries, 1, 1, p->Table16, CMS_LERP_FLAGS_16BITS); + if (p->InterpParams != NULL) + return p; + +Error: + if (p -> Segments) _cmsFree(ContextID, p ->Segments); + if (p -> Evals) _cmsFree(ContextID, p -> Evals); + if (p ->Table16) _cmsFree(ContextID, p ->Table16); + _cmsFree(ContextID, p); + return NULL; +} + + +// Parametric Fn using floating point +static +cmsFloat64Number DefaultEvalParametricFn(cmsInt32Number Type, const cmsFloat64Number Params[], cmsFloat64Number R) +{ + cmsFloat64Number e, Val, disc; + + switch (Type) { + + // X = Y ^ Gamma + case 1: + if (R < 0) { + + if (fabs(Params[0] - 1.0) < MATRIX_DET_TOLERANCE) + Val = R; + else + Val = 0; + } + else + Val = pow(R, Params[0]); + break; + + // Type 1 Reversed: X = Y ^1/gamma + case -1: + if (R < 0) { + + if (fabs(Params[0] - 1.0) < MATRIX_DET_TOLERANCE) + Val = R; + else + Val = 0; + } + else + Val = pow(R, 1/Params[0]); + break; + + // CIE 122-1966 + // Y = (aX + b)^Gamma | X >= -b/a + // Y = 0 | else + case 2: + disc = -Params[2] / Params[1]; + + if (R >= disc ) { + + e = Params[1]*R + Params[2]; + + if (e > 0) + Val = pow(e, Params[0]); + else + Val = 0; + } + else + Val = 0; + break; + + // Type 2 Reversed + // X = (Y ^1/g - b) / a + case -2: + if (R < 0) + Val = 0; + else + Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1]; + + if (Val < 0) + Val = 0; + break; + + + // IEC 61966-3 + // Y = (aX + b)^Gamma | X <= -b/a + // Y = c | else + case 3: + disc = -Params[2] / Params[1]; + if (disc < 0) + disc = 0; + + if (R >= disc) { + + e = Params[1]*R + Params[2]; + + if (e > 0) + Val = pow(e, Params[0]) + Params[3]; + else + Val = 0; + } + else + Val = Params[3]; + break; + + + // Type 3 reversed + // X=((Y-c)^1/g - b)/a | (Y>=c) + // X=-b/a | (Y<c) + case -3: + if (R >= Params[3]) { + + e = R - Params[3]; + + if (e > 0) + Val = (pow(e, 1/Params[0]) - Params[2]) / Params[1]; + else + Val = 0; + } + else { + Val = -Params[2] / Params[1]; + } + break; + + + // IEC 61966-2.1 (sRGB) + // Y = (aX + b)^Gamma | X >= d + // Y = cX | X < d + case 4: + if (R >= Params[4]) { + + e = Params[1]*R + Params[2]; + + if (e > 0) + Val = pow(e, Params[0]); + else + Val = 0; + } + else + Val = R * Params[3]; + break; + + // Type 4 reversed + // X=((Y^1/g-b)/a) | Y >= (ad+b)^g + // X=Y/c | Y< (ad+b)^g + case -4: + e = Params[1] * Params[4] + Params[2]; + if (e < 0) + disc = 0; + else + disc = pow(e, Params[0]); + + if (R >= disc) { + + Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1]; + } + else { + Val = R / Params[3]; + } + break; + + + // Y = (aX + b)^Gamma + e | X >= d + // Y = cX + f | X < d + case 5: + if (R >= Params[4]) { + + e = Params[1]*R + Params[2]; + + if (e > 0) + Val = pow(e, Params[0]) + Params[5]; + else + Val = Params[5]; + } + else + Val = R*Params[3] + Params[6]; + break; + + + // Reversed type 5 + // X=((Y-e)1/g-b)/a | Y >=(ad+b)^g+e), cd+f + // X=(Y-f)/c | else + case -5: + + disc = Params[3] * Params[4] + Params[6]; + if (R >= disc) { + + e = R - Params[5]; + if (e < 0) + Val = 0; + else + Val = (pow(e, 1.0/Params[0]) - Params[2]) / Params[1]; + } + else { + Val = (R - Params[6]) / Params[3]; + } + break; + + + // Types 6,7,8 comes from segmented curves as described in ICCSpecRevision_02_11_06_Float.pdf + // Type 6 is basically identical to type 5 without d + + // Y = (a * X + b) ^ Gamma + c + case 6: + e = Params[1]*R + Params[2]; + + if (e < 0) + Val = Params[3]; + else + Val = pow(e, Params[0]) + Params[3]; + break; + + // ((Y - c) ^1/Gamma - b) / a + case -6: + e = R - Params[3]; + if (e < 0) + Val = 0; + else + Val = (pow(e, 1.0/Params[0]) - Params[2]) / Params[1]; + break; + + + // Y = a * log (b * X^Gamma + c) + d + case 7: + + e = Params[2] * pow(R, Params[0]) + Params[3]; + if (e <= 0) + Val = Params[4]; + else + Val = Params[1]*log10(e) + Params[4]; + break; + + // (Y - d) / a = log(b * X ^Gamma + c) + // pow(10, (Y-d) / a) = b * X ^Gamma + c + // pow((pow(10, (Y-d) / a) - c) / b, 1/g) = X + case -7: + Val = pow((pow(10.0, (R-Params[4]) / Params[1]) - Params[3]) / Params[2], 1.0 / Params[0]); + break; + + + //Y = a * b^(c*X+d) + e + case 8: + Val = (Params[0] * pow(Params[1], Params[2] * R + Params[3]) + Params[4]); + break; + + + // Y = (log((y-e) / a) / log(b) - d ) / c + // a=0, b=1, c=2, d=3, e=4, + case -8: + + disc = R - Params[4]; + if (disc < 0) Val = 0; + else + Val = (log(disc / Params[0]) / log(Params[1]) - Params[3]) / Params[2]; + break; + + // S-Shaped: (1 - (1-x)^1/g)^1/g + case 108: + Val = pow(1.0 - pow(1 - R, 1/Params[0]), 1/Params[0]); + break; + + // y = (1 - (1-x)^1/g)^1/g + // y^g = (1 - (1-x)^1/g) + // 1 - y^g = (1-x)^1/g + // (1 - y^g)^g = 1 - x + // 1 - (1 - y^g)^g + case -108: + Val = 1 - pow(1 - pow(R, Params[0]), Params[0]); + break; + + default: + // Unsupported parametric curve. Should never reach here + return 0; + } + + return Val; +} + +// Evaluate a segmented funtion for a single value. Return -1 if no valid segment found . +// If fn type is 0, perform an interpolation on the table +static +cmsFloat64Number EvalSegmentedFn(const cmsToneCurve *g, cmsFloat64Number R) +{ + int i; + + for (i = g ->nSegments-1; i >= 0 ; --i) { + + // Check for domain + if ((R > g ->Segments[i].x0) && (R <= g ->Segments[i].x1)) { + + // Type == 0 means segment is sampled + if (g ->Segments[i].Type == 0) { + + cmsFloat32Number R1 = (cmsFloat32Number) (R - g ->Segments[i].x0) / (g ->Segments[i].x1 - g ->Segments[i].x0); + cmsFloat32Number Out; + + // Setup the table (TODO: clean that) + g ->SegInterp[i]-> Table = g ->Segments[i].SampledPoints; + + g ->SegInterp[i] -> Interpolation.LerpFloat(&R1, &Out, g ->SegInterp[i]); + + return Out; + } + else + return g ->Evals[i](g->Segments[i].Type, g ->Segments[i].Params, R); + } + } + + return MINUS_INF; +} + +// Access to estimated low-res table +cmsUInt32Number CMSEXPORT cmsGetToneCurveEstimatedTableEntries(const cmsToneCurve* t) +{ + _cmsAssert(t != NULL); + return t ->nEntries; +} + +const cmsUInt16Number* CMSEXPORT cmsGetToneCurveEstimatedTable(const cmsToneCurve* t) +{ + _cmsAssert(t != NULL); + return t ->Table16; +} + + +// Create an empty gamma curve, by using tables. This specifies only the limited-precision part, and leaves the +// floating point description empty. +cmsToneCurve* CMSEXPORT cmsBuildTabulatedToneCurve16(cmsContext ContextID, cmsInt32Number nEntries, const cmsUInt16Number Values[]) +{ + return AllocateToneCurveStruct(ContextID, nEntries, 0, NULL, Values); +} + +static +int EntriesByGamma(cmsFloat64Number Gamma) +{ + if (fabs(Gamma - 1.0) < 0.001) return 2; + return 4096; +} + + +// Create a segmented gamma, fill the table +cmsToneCurve* CMSEXPORT cmsBuildSegmentedToneCurve(cmsContext ContextID, + cmsInt32Number nSegments, const cmsCurveSegment Segments[]) +{ + int i; + cmsFloat64Number R, Val; + cmsToneCurve* g; + int nGridPoints = 4096; + + _cmsAssert(Segments != NULL); + + // Optimizatin for identity curves. + if (nSegments == 1 && Segments[0].Type == 1) { + + nGridPoints = EntriesByGamma(Segments[0].Params[0]); + } + + g = AllocateToneCurveStruct(ContextID, nGridPoints, nSegments, Segments, NULL); + if (g == NULL) return NULL; + + // Once we have the floating point version, we can approximate a 16 bit table of 4096 entries + // for performance reasons. This table would normally not be used except on 8/16 bits transforms. + for (i=0; i < nGridPoints; i++) { + + R = (cmsFloat64Number) i / (nGridPoints-1); + + Val = EvalSegmentedFn(g, R); + + // Round and saturate + g ->Table16[i] = _cmsQuickSaturateWord(Val * 65535.0); + } + + return g; +} + +// Use a segmented curve to store the floating point table +cmsToneCurve* CMSEXPORT cmsBuildTabulatedToneCurveFloat(cmsContext ContextID, cmsUInt32Number nEntries, const cmsFloat32Number values[]) +{ + cmsCurveSegment Seg[3]; + + // A segmented tone curve should have function segments in the first and last positions + // Initialize segmented curve part up to 0 to constant value = samples[0] + Seg[0].x0 = MINUS_INF; + Seg[0].x1 = 0; + Seg[0].Type = 6; + + Seg[0].Params[0] = 1; + Seg[0].Params[1] = 0; + Seg[0].Params[2] = 0; + Seg[0].Params[3] = values[0]; + Seg[0].Params[4] = 0; + + // From zero to 1 + Seg[1].x0 = 0; + Seg[1].x1 = 1.0; + Seg[1].Type = 0; + + Seg[1].nGridPoints = nEntries; + Seg[1].SampledPoints = (cmsFloat32Number*) values; + + // Final segment is constant = lastsample + Seg[2].x0 = 1.0; + Seg[2].x1 = PLUS_INF; + Seg[2].Type = 6; + + Seg[2].Params[0] = 1; + Seg[2].Params[1] = 0; + Seg[2].Params[2] = 0; + Seg[2].Params[3] = values[nEntries-1]; + Seg[2].Params[4] = 0; + + + return cmsBuildSegmentedToneCurve(ContextID, 3, Seg); +} + +// Parametric curves +// +// Parameters goes as: Curve, a, b, c, d, e, f +// Type is the ICC type +1 +// if type is negative, then the curve is analyticaly inverted +cmsToneCurve* CMSEXPORT cmsBuildParametricToneCurve(cmsContext ContextID, cmsInt32Number Type, const cmsFloat64Number Params[]) +{ + cmsCurveSegment Seg0; + int Pos = 0; + cmsUInt32Number size; + _cmsParametricCurvesCollection* c = GetParametricCurveByType(ContextID, Type, &Pos); + + _cmsAssert(Params != NULL); + + if (c == NULL) { + cmsSignalError(ContextID, cmsERROR_UNKNOWN_EXTENSION, "Invalid parametric curve type %d", Type); + return NULL; + } + + memset(&Seg0, 0, sizeof(Seg0)); + + Seg0.x0 = MINUS_INF; + Seg0.x1 = PLUS_INF; + Seg0.Type = Type; + + size = c->ParameterCount[Pos] * sizeof(cmsFloat64Number); + memmove(Seg0.Params, Params, size); + + return cmsBuildSegmentedToneCurve(ContextID, 1, &Seg0); +} + + + +// Build a gamma table based on gamma constant +cmsToneCurve* CMSEXPORT cmsBuildGamma(cmsContext ContextID, cmsFloat64Number Gamma) +{ + return cmsBuildParametricToneCurve(ContextID, 1, &Gamma); +} + + +// Free all memory taken by the gamma curve +void CMSEXPORT cmsFreeToneCurve(cmsToneCurve* Curve) +{ + cmsContext ContextID; + + // added by Xiaochuan Liu + // Curve->InterpParams may be null + if (Curve == NULL || Curve->InterpParams == NULL) return; + + ContextID = Curve ->InterpParams->ContextID; + + _cmsFreeInterpParams(Curve ->InterpParams); + Curve ->InterpParams = NULL; + + if (Curve -> Table16) + { + _cmsFree(ContextID, Curve ->Table16); + Curve ->Table16 = NULL; + } + + if (Curve ->Segments) { + + cmsUInt32Number i; + + for (i=0; i < Curve ->nSegments; i++) { + + if (Curve ->Segments[i].SampledPoints) { + _cmsFree(ContextID, Curve ->Segments[i].SampledPoints); + Curve ->Segments[i].SampledPoints = NULL; + } + + if (Curve ->SegInterp[i] != 0) + { + _cmsFreeInterpParams(Curve->SegInterp[i]); + Curve->SegInterp[i] = NULL; + } + } + + _cmsFree(ContextID, Curve ->Segments); + Curve ->Segments = NULL; + _cmsFree(ContextID, Curve ->SegInterp); + Curve ->SegInterp = NULL; + } + + if (Curve -> Evals) + { + _cmsFree(ContextID, Curve -> Evals); + Curve -> Evals = NULL; + } + + if (Curve) + { + _cmsFree(ContextID, Curve); + Curve = NULL; + } +} + +// Utility function, free 3 gamma tables +void CMSEXPORT cmsFreeToneCurveTriple(cmsToneCurve* Curve[3]) +{ + + _cmsAssert(Curve != NULL); + + if (Curve[0] != NULL) cmsFreeToneCurve(Curve[0]); + if (Curve[1] != NULL) cmsFreeToneCurve(Curve[1]); + if (Curve[2] != NULL) cmsFreeToneCurve(Curve[2]); + + Curve[0] = Curve[1] = Curve[2] = NULL; +} + + +// Duplicate a gamma table +cmsToneCurve* CMSEXPORT cmsDupToneCurve(const cmsToneCurve* In) +{ + // Xiaochuan Liu + // fix openpdf bug(mantis id:0055683, google id:360198) + // the function CurveSetElemTypeFree in cmslut.c also needs to check pointer + if (In == NULL || In ->InterpParams == NULL || In ->Segments == NULL || In ->Table16 == NULL) return NULL; + + return AllocateToneCurveStruct(In ->InterpParams ->ContextID, In ->nEntries, In ->nSegments, In ->Segments, In ->Table16); +} + +// Joins two curves for X and Y. Curves should be monotonic. +// We want to get +// +// y = Y^-1(X(t)) +// +cmsToneCurve* CMSEXPORT cmsJoinToneCurve(cmsContext ContextID, + const cmsToneCurve* X, + const cmsToneCurve* Y, cmsUInt32Number nResultingPoints) +{ + cmsToneCurve* out = NULL; + cmsToneCurve* Yreversed = NULL; + cmsFloat32Number t, x; + cmsFloat32Number* Res = NULL; + cmsUInt32Number i; + + + _cmsAssert(X != NULL); + _cmsAssert(Y != NULL); + + Yreversed = cmsReverseToneCurveEx(nResultingPoints, Y); + if (Yreversed == NULL) goto Error; + + Res = (cmsFloat32Number*) _cmsCalloc(ContextID, nResultingPoints, sizeof(cmsFloat32Number)); + if (Res == NULL) goto Error; + + //Iterate + for (i=0; i < nResultingPoints; i++) { + + t = (cmsFloat32Number) i / (nResultingPoints-1); + x = cmsEvalToneCurveFloat(X, t); + Res[i] = cmsEvalToneCurveFloat(Yreversed, x); + } + + // Allocate space for output + out = cmsBuildTabulatedToneCurveFloat(ContextID, nResultingPoints, Res); + +Error: + + if (Res != NULL) _cmsFree(ContextID, Res); + if (Yreversed != NULL) cmsFreeToneCurve(Yreversed); + + return out; +} + + + +// Get the surrounding nodes. This is tricky on non-monotonic tables +static +int GetInterval(cmsFloat64Number In, const cmsUInt16Number LutTable[], const struct _cms_interp_struc* p) +{ + int i; + int y0, y1; + + // A 1 point table is not allowed + if (p -> Domain[0] < 1) return -1; + + // Let's see if ascending or descending. + if (LutTable[0] < LutTable[p ->Domain[0]]) { + + // Table is overall ascending + for (i=p->Domain[0]-1; i >=0; --i) { + + y0 = LutTable[i]; + y1 = LutTable[i+1]; + + if (y0 <= y1) { // Increasing + if (In >= y0 && In <= y1) return i; + } + else + if (y1 < y0) { // Decreasing + if (In >= y1 && In <= y0) return i; + } + } + } + else { + // Table is overall descending + for (i=0; i < (int) p -> Domain[0]; i++) { + + y0 = LutTable[i]; + y1 = LutTable[i+1]; + + if (y0 <= y1) { // Increasing + if (In >= y0 && In <= y1) return i; + } + else + if (y1 < y0) { // Decreasing + if (In >= y1 && In <= y0) return i; + } + } + } + + return -1; +} + +// Reverse a gamma table +cmsToneCurve* CMSEXPORT cmsReverseToneCurveEx(cmsInt32Number nResultSamples, const cmsToneCurve* InCurve) +{ + cmsToneCurve *out; + cmsFloat64Number a = 0, b = 0, y, x1, y1, x2, y2; + int i, j; + int Ascending; + + _cmsAssert(InCurve != NULL); + + // Try to reverse it analytically whatever possible + + if (InCurve ->nSegments == 1 && InCurve ->Segments[0].Type > 0 && + /* InCurve -> Segments[0].Type <= 5 */ + GetParametricCurveByType(InCurve ->InterpParams->ContextID, InCurve ->Segments[0].Type, NULL) != NULL) { + + return cmsBuildParametricToneCurve(InCurve ->InterpParams->ContextID, + -(InCurve -> Segments[0].Type), + InCurve -> Segments[0].Params); + } + + // Nope, reverse the table. + out = cmsBuildTabulatedToneCurve16(InCurve ->InterpParams->ContextID, nResultSamples, NULL); + if (out == NULL) + return NULL; + + // We want to know if this is an ascending or descending table + Ascending = !cmsIsToneCurveDescending(InCurve); + + // Iterate across Y axis + for (i=0; i < nResultSamples; i++) { + + y = (cmsFloat64Number) i * 65535.0 / (nResultSamples - 1); + + // Find interval in which y is within. + j = GetInterval(y, InCurve->Table16, InCurve->InterpParams); + if (j >= 0) { + + + // Get limits of interval + x1 = InCurve ->Table16[j]; + x2 = InCurve ->Table16[j+1]; + + y1 = (cmsFloat64Number) (j * 65535.0) / (InCurve ->nEntries - 1); + y2 = (cmsFloat64Number) ((j+1) * 65535.0 ) / (InCurve ->nEntries - 1); + + // If collapsed, then use any + if (x1 == x2) { + + out ->Table16[i] = _cmsQuickSaturateWord(Ascending ? y2 : y1); + continue; + + } else { + + // Interpolate + a = (y2 - y1) / (x2 - x1); + b = y2 - a * x2; + } + } + + out ->Table16[i] = _cmsQuickSaturateWord(a* y + b); + } + + + return out; +} + +// Reverse a gamma table +cmsToneCurve* CMSEXPORT cmsReverseToneCurve(const cmsToneCurve* InGamma) +{ + _cmsAssert(InGamma != NULL); + + return cmsReverseToneCurveEx(4096, InGamma); +} + +// From: Eilers, P.H.C. (1994) Smoothing and interpolation with finite +// differences. in: Graphic Gems IV, Heckbert, P.S. (ed.), Academic press. +// +// Smoothing and interpolation with second differences. +// +// Input: weights (w), data (y): vector from 1 to m. +// Input: smoothing parameter (lambda), length (m). +// Output: smoothed vector (z): vector from 1 to m. + +static +cmsBool smooth2(cmsContext ContextID, cmsFloat32Number w[], cmsFloat32Number y[], cmsFloat32Number z[], cmsFloat32Number lambda, int m) +{ + int i, i1, i2; + cmsFloat32Number *c, *d, *e; + cmsBool st; + + + c = (cmsFloat32Number*) _cmsCalloc(ContextID, MAX_NODES_IN_CURVE, sizeof(cmsFloat32Number)); + d = (cmsFloat32Number*) _cmsCalloc(ContextID, MAX_NODES_IN_CURVE, sizeof(cmsFloat32Number)); + e = (cmsFloat32Number*) _cmsCalloc(ContextID, MAX_NODES_IN_CURVE, sizeof(cmsFloat32Number)); + + if (c != NULL && d != NULL && e != NULL) { + + + d[1] = w[1] + lambda; + c[1] = -2 * lambda / d[1]; + e[1] = lambda /d[1]; + z[1] = w[1] * y[1]; + d[2] = w[2] + 5 * lambda - d[1] * c[1] * c[1]; + c[2] = (-4 * lambda - d[1] * c[1] * e[1]) / d[2]; + e[2] = lambda / d[2]; + z[2] = w[2] * y[2] - c[1] * z[1]; + + for (i = 3; i < m - 1; i++) { + i1 = i - 1; i2 = i - 2; + d[i]= w[i] + 6 * lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2]; + c[i] = (-4 * lambda -d[i1] * c[i1] * e[i1])/ d[i]; + e[i] = lambda / d[i]; + z[i] = w[i] * y[i] - c[i1] * z[i1] - e[i2] * z[i2]; + } + + i1 = m - 2; i2 = m - 3; + + d[m - 1] = w[m - 1] + 5 * lambda -c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2]; + c[m - 1] = (-2 * lambda - d[i1] * c[i1] * e[i1]) / d[m - 1]; + z[m - 1] = w[m - 1] * y[m - 1] - c[i1] * z[i1] - e[i2] * z[i2]; + i1 = m - 1; i2 = m - 2; + + d[m] = w[m] + lambda - c[i1] * c[i1] * d[i1] - e[i2] * e[i2] * d[i2]; + z[m] = (w[m] * y[m] - c[i1] * z[i1] - e[i2] * z[i2]) / d[m]; + z[m - 1] = z[m - 1] / d[m - 1] - c[m - 1] * z[m]; + + for (i = m - 2; 1<= i; i--) + z[i] = z[i] / d[i] - c[i] * z[i + 1] - e[i] * z[i + 2]; + + st = TRUE; + } + else st = FALSE; + + if (c != NULL) _cmsFree(ContextID, c); + if (d != NULL) _cmsFree(ContextID, d); + if (e != NULL) _cmsFree(ContextID, e); + + return st; +} + +// Smooths a curve sampled at regular intervals. +cmsBool CMSEXPORT cmsSmoothToneCurve(cmsToneCurve* Tab, cmsFloat64Number lambda) +{ + cmsFloat32Number w[MAX_NODES_IN_CURVE], y[MAX_NODES_IN_CURVE], z[MAX_NODES_IN_CURVE]; + int i, nItems, Zeros, Poles; + + if (Tab == NULL) return FALSE; + + if (cmsIsToneCurveLinear(Tab)) return TRUE; // Nothing to do + + nItems = Tab -> nEntries; + + if (nItems >= MAX_NODES_IN_CURVE) { + cmsSignalError(Tab ->InterpParams->ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: too many points."); + return FALSE; + } + + memset(w, 0, nItems * sizeof(cmsFloat32Number)); + memset(y, 0, nItems * sizeof(cmsFloat32Number)); + memset(z, 0, nItems * sizeof(cmsFloat32Number)); + + for (i=0; i < nItems; i++) + { + y[i+1] = (cmsFloat32Number) Tab -> Table16[i]; + w[i+1] = 1.0; + } + + if (!smooth2(Tab ->InterpParams->ContextID, w, y, z, (cmsFloat32Number) lambda, nItems)) return FALSE; + + // Do some reality - checking... + Zeros = Poles = 0; + for (i=nItems; i > 1; --i) { + + if (z[i] == 0.) Zeros++; + if (z[i] >= 65535.) Poles++; + if (z[i] < z[i-1]) { + cmsSignalError(Tab ->InterpParams->ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Non-Monotonic."); + return FALSE; + } + } + + if (Zeros > (nItems / 3)) { + cmsSignalError(Tab ->InterpParams->ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Degenerated, mostly zeros."); + return FALSE; + } + if (Poles > (nItems / 3)) { + cmsSignalError(Tab ->InterpParams->ContextID, cmsERROR_RANGE, "cmsSmoothToneCurve: Degenerated, mostly poles."); + return FALSE; + } + + // Seems ok + for (i=0; i < nItems; i++) { + + // Clamp to cmsUInt16Number + Tab -> Table16[i] = _cmsQuickSaturateWord(z[i+1]); + } + + return TRUE; +} + +// Is a table linear? Do not use parametric since we cannot guarantee some weird parameters resulting +// in a linear table. This way assures it is linear in 12 bits, which should be enought in most cases. +cmsBool CMSEXPORT cmsIsToneCurveLinear(const cmsToneCurve* Curve) +{ + cmsUInt32Number i; + int diff; + + _cmsAssert(Curve != NULL); + + for (i=0; i < Curve ->nEntries; i++) { + + diff = abs((int) Curve->Table16[i] - (int) _cmsQuantizeVal(i, Curve ->nEntries)); + if (diff > 0x0f) + return FALSE; + } + + return TRUE; +} + +// Same, but for monotonicity +cmsBool CMSEXPORT cmsIsToneCurveMonotonic(const cmsToneCurve* t) +{ + int n; + int i, last; + cmsBool lDescending; + + _cmsAssert(t != NULL); + + // Degenerated curves are monotonic? Ok, let's pass them + n = t ->nEntries; + if (n < 2) return TRUE; + + // Curve direction + lDescending = cmsIsToneCurveDescending(t); + + if (lDescending) { + + last = t ->Table16[0]; + + for (i = 1; i < n; i++) { + + if (t ->Table16[i] - last > 2) // We allow some ripple + return FALSE; + else + last = t ->Table16[i]; + + } + } + else { + + last = t ->Table16[n-1]; + + for (i = n-2; i >= 0; --i) { + + if (t ->Table16[i] - last > 2) + return FALSE; + else + last = t ->Table16[i]; + + } + } + + return TRUE; +} + +// Same, but for descending tables +cmsBool CMSEXPORT cmsIsToneCurveDescending(const cmsToneCurve* t) +{ + _cmsAssert(t != NULL); + + return t ->Table16[0] > t ->Table16[t ->nEntries-1]; +} + + +// Another info fn: is out gamma table multisegment? +cmsBool CMSEXPORT cmsIsToneCurveMultisegment(const cmsToneCurve* t) +{ + _cmsAssert(t != NULL); + + return t -> nSegments > 1; +} + +cmsInt32Number CMSEXPORT cmsGetToneCurveParametricType(const cmsToneCurve* t) +{ + _cmsAssert(t != NULL); + + if (t -> nSegments != 1) return 0; + return t ->Segments[0].Type; +} + +// We need accuracy this time +cmsFloat32Number CMSEXPORT cmsEvalToneCurveFloat(const cmsToneCurve* Curve, cmsFloat32Number v) +{ + _cmsAssert(Curve != NULL); + + // Check for 16 bits table. If so, this is a limited-precision tone curve + if (Curve ->nSegments == 0) { + + cmsUInt16Number In, Out; + + In = (cmsUInt16Number) _cmsQuickSaturateWord(v * 65535.0); + Out = cmsEvalToneCurve16(Curve, In); + + return (cmsFloat32Number) (Out / 65535.0); + } + + return (cmsFloat32Number) EvalSegmentedFn(Curve, v); +} + +// We need xput over here +cmsUInt16Number CMSEXPORT cmsEvalToneCurve16(const cmsToneCurve* Curve, cmsUInt16Number v) +{ + cmsUInt16Number out; + + _cmsAssert(Curve != NULL); + + Curve ->InterpParams ->Interpolation.Lerp16(&v, &out, Curve ->InterpParams); + return out; +} + + +// Least squares fitting. +// A mathematical procedure for finding the best-fitting curve to a given set of points by +// minimizing the sum of the squares of the offsets ("the residuals") of the points from the curve. +// The sum of the squares of the offsets is used instead of the offset absolute values because +// this allows the residuals to be treated as a continuous differentiable quantity. +// +// y = f(x) = x ^ g +// +// R = (yi - (xi^g)) +// R2 = (yi - (xi^g))2 +// SUM R2 = SUM (yi - (xi^g))2 +// +// dR2/dg = -2 SUM x^g log(x)(y - x^g) +// solving for dR2/dg = 0 +// +// g = 1/n * SUM(log(y) / log(x)) + +cmsFloat64Number CMSEXPORT cmsEstimateGamma(const cmsToneCurve* t, cmsFloat64Number Precision) +{ + cmsFloat64Number gamma, sum, sum2; + cmsFloat64Number n, x, y, Std; + cmsUInt32Number i; + + _cmsAssert(t != NULL); + + sum = sum2 = n = 0; + + // Excluding endpoints + for (i=1; i < (MAX_NODES_IN_CURVE-1); i++) { + + x = (cmsFloat64Number) i / (MAX_NODES_IN_CURVE-1); + y = (cmsFloat64Number) cmsEvalToneCurveFloat(t, (cmsFloat32Number) x); + + // Avoid 7% on lower part to prevent + // artifacts due to linear ramps + + if (y > 0. && y < 1. && x > 0.07) { + + gamma = log(y) / log(x); + sum += gamma; + sum2 += gamma * gamma; + n++; + } + } + + // Take a look on SD to see if gamma isn't exponential at all + Std = sqrt((n * sum2 - sum * sum) / (n*(n-1))); + + if (Std > Precision) + return -1.0; + + return (sum / n); // The mean +} |