diff options
Diffstat (limited to 'icc/icc.c')
-rw-r--r-- | icc/icc.c | 1191 |
1 files changed, 1043 insertions, 148 deletions
@@ -7,7 +7,7 @@ * Date: 2002/04/22 * Version: 2.15 * - * Copyright 1997 - 2012 Graeme W. Gill + * Copyright 1997 - 2013 Graeme W. Gill * * This material is licensed with an "MIT" free use license:- * see the License.txt file in this directory for licensing details. @@ -51,10 +51,10 @@ #define _ICC_C_ /* Turn on implimentation code */ -#undef DEBUG_SETLUT /* Show each value being set in setting lut contents */ -#undef DEBUG_SETLUT_CLIP /* Show clipped values when setting LUT */ -#undef DEBUG_LULUT /* Show each value being looked up from lut contents */ -#undef DEBUG_LLULUT /* Debug individual lookup steps (not fully implemented) */ +#undef DEBUG_SETLUT /* [Und] Show each value being set in setting lut contents */ +#undef DEBUG_SETLUT_CLIP /* [Und] Show clipped values when setting LUT */ +#undef DEBUG_LULUT /* [Und] Show each value being looked up from lut contents */ +#undef DEBUG_LLULUT /* [Und] Debug individual lookup steps (not fully implemented) */ #include <stdio.h> #include <stdlib.h> @@ -127,6 +127,10 @@ #define DBLLL(xxx) #endif +#ifndef M_PI +# define M_PI 3.14159265358979323846 +#endif + /* =========================================================== */ /* Overflow protected unsigned int arithmatic functions. */ /* These functions saturate rather than wrapping around. */ @@ -706,6 +710,12 @@ static int write_S15Fixed16Number(double d, char *p) { return 0; } +static double round_S15Fixed16Number(double d) { + d = floor(d * 65536.0 + 0.5); /* Beware! (int)(d + 0.5) doesn't work! */ + d = d/65536.0; + return d; +} + /* Device coordinate as 8 bit value range 0.0 - 1.0 */ static double read_DCS8Number(char *p) { unsigned int rv; @@ -775,7 +785,7 @@ static void read_PCSNumber(icc *icp, icColorSpaceSignature csig, double pcs[3], if (csig == icmSigPCSData) csig = icp->header->pcs; if (csig == icSigLabData) { - if (icp->ver != 0) + if (icp->ver >= icmVersion4_1) csig = icmSigLabV4Data; else csig = icmSigLabV2Data; @@ -820,7 +830,7 @@ static int write_PCSNumber(icc *icp, icColorSpaceSignature csig, double pcs[3], if (csig == icmSigPCSData) csig = icp->header->pcs; if (csig == icSigLabData) { - if (icp->ver != 0) + if (icp->ver >= icmVersion4_1) csig = icmSigLabV4Data; else csig = icmSigLabV2Data; @@ -5374,6 +5384,257 @@ double *in /* Input array[outputChan] */ } /* ----------------------------------------------- */ +/* Tune a single interpolated value. Based on lookup_clut functions (above) */ + +/* Helper function to fine tune a single value interpolation */ +/* Return 0 on success, 1 if input clipping occured, 2 if output clipping occured */ +int icmLut_tune_value_nl( +icmLut *p, /* Pointer to Lut object */ +double *out, /* Output array[inputChan] */ +double *in /* Input array[outputChan] */ +) { + icc *icp = p->icp; + int rv = 0; + double *gp; /* Pointer to grid cube base */ + double co[MAX_CHAN]; /* Coordinate offset with the grid cell */ + double *gw, GW[1 << 8]; /* weight for each grid cube corner */ + double cout[MAX_CHAN]; /* Current output value */ + + if (p->inputChan <= 8) { + gw = GW; /* Use stack allocation */ + } else { + if ((gw = (double *) icp->al->malloc(icp->al, sat_mul((1 << p->inputChan), sizeof(double)))) == NULL) { + sprintf(icp->err,"icmLut_lookup_clut: malloc() failed"); + return icp->errc = 2; + } + } + + /* We are using an multi-linear (ie. Trilinear for 3D input) interpolation. */ + /* The implementation here uses more multiplies that some other schemes, */ + /* (for instance, see "Tri-Linear Interpolation" by Steve Hill, */ + /* Graphics Gems IV, page 521), but has less involved bookeeping, */ + /* needs less local storage for intermediate output values, does fewer */ + /* output and intermediate value reads, and fp multiplies are fast on */ + /* todays processors! */ + + /* Compute base index into grid and coordinate offsets */ + { + unsigned int e; + double clutPoints_1 = (double)(p->clutPoints-1); + int clutPoints_2 = p->clutPoints-2; + gp = p->clutTable; /* Base of grid array */ + + for (e = 0; e < p->inputChan; e++) { + unsigned int x; + double val; + val = in[e] * clutPoints_1; + if (val < 0.0) { + val = 0.0; + rv |= 1; + } else if (val > clutPoints_1) { + val = clutPoints_1; + rv |= 1; + } + x = (unsigned int)floor(val); /* Grid coordinate */ + if (x > clutPoints_2) + x = clutPoints_2; + co[e] = val - (double)x; /* 1.0 - weight */ + gp += x * p->dinc[e]; /* Add index offset for base of cube */ + } + } + /* Compute corner weights needed for interpolation */ + { + unsigned int e; + int i, g = 1; + gw[0] = 1.0; + for (e = 0; e < p->inputChan; e++) { + for (i = 0; i < g; i++) { + gw[g+i] = gw[i] * co[e]; + gw[i] *= (1.0 - co[e]); + } + g *= 2; + } + } + /* Now compute the current output value, and distribute the correction */ + { + int i; + unsigned int f; + double w, *d, ww = 0.0; + for (f = 0; f < p->outputChan; f++) + cout[f] = 0.0; + for (i = 0; i < (1 << p->inputChan); i++) { /* For all other corners of cube */ + w = gw[i]; /* Strength reduce */ + ww += w * w; /* Sum of weights squared */ + d = gp + p->dcube[i]; + for (f = 0; f < p->outputChan; f++) + cout[f] += w * d[f]; + } + + /* We distribute the correction needed in proportion to the */ + /* interpolation weighting, so the biggest correction is to the */ + /* closest vertex. */ + + for (f = 0; f < p->outputChan; f++) + cout[f] = (out[f] - cout[f])/ww; /* Amount to distribute */ + + for (i = 0; i < (1 << p->inputChan); i++) { /* For all other corners of cube */ + w = gw[i]; /* Strength reduce */ + d = gp + p->dcube[i]; + for (f = 0; f < p->outputChan; f++) { + d[f] += w * cout[f]; /* Apply correction */ + if (d[f] < 0.0) { + d[f] = 0.0; + rv |= 2; + } else if (d[f] > 1.0) { + d[f] = 1.0; + rv |= 2; + } + } + } + } + + if (gw != GW) + icp->al->free(icp->al, (void *)gw); + return rv; +} + +/* Helper function to fine tune a single value interpolation */ +/* Return 0 on success, 1 if input clipping occured, 2 if output clipping occured */ +int icmLut_tune_value_sx( +icmLut *p, /* Pointer to Lut object */ +double *out, /* Output array[inputChan] */ +double *in /* Input array[outputChan] */ +) { + int rv = 0; + double *gp; /* Pointer to grid cube base */ + double co[MAX_CHAN]; /* Coordinate offset with the grid cell */ + int si[MAX_CHAN]; /* co[] Sort index, [0] = smalest */ + double cout[MAX_CHAN]; /* Current output value */ + + /* We are using a simplex (ie. tetrahedral for 3D input) interpolation. */ + /* This method is more appropriate for XYZ/RGB/CMYK input spaces, */ + + /* Compute base index into grid and coordinate offsets */ + { + unsigned int e; + double clutPoints_1 = (double)(p->clutPoints-1); + int clutPoints_2 = p->clutPoints-2; + gp = p->clutTable; /* Base of grid array */ + + for (e = 0; e < p->inputChan; e++) { + unsigned int x; + double val; + val = in[e] * clutPoints_1; + if (val < 0.0) { + val = 0.0; + rv |= 1; + } else if (val > clutPoints_1) { + val = clutPoints_1; + rv |= 1; + } + x = (unsigned int)floor(val); /* Grid coordinate */ + if (x > clutPoints_2) + x = clutPoints_2; + co[e] = val - (double)x; /* 1.0 - weight */ + gp += x * p->dinc[e]; /* Add index offset for base of cube */ + } + } + /* Do insertion sort on coordinates, smallest to largest. */ + { + int f, vf; + unsigned int e; + double v; + for (e = 0; e < p->inputChan; e++) + si[e] = e; /* Initial unsorted indexes */ + + for (e = 1; e < p->inputChan; e++) { + f = e; + v = co[si[f]]; + vf = f; + while (f > 0 && co[si[f-1]] > v) { + si[f] = si[f-1]; + f--; + } + si[f] = vf; + } + } + /* Now compute the current output value, and distribute the correction */ + { + unsigned int e, f; + double w, ww = 0.0; /* Current vertex weight, sum of weights squared */ + double *ogp = gp; /* Pointer to grid cube base */ + + w = 1.0 - co[si[p->inputChan-1]]; /* Vertex at base of cell */ + ww += w * w; /* Sum of weights squared */ + for (f = 0; f < p->outputChan; f++) + cout[f] = w * gp[f]; + + for (e = p->inputChan-1; e > 0; e--) { /* Middle verticies */ + w = co[si[e]] - co[si[e-1]]; + ww += w * w; /* Sum of weights squared */ + gp += p->dinc[si[e]]; /* Move to top of cell in next largest dimension */ + for (f = 0; f < p->outputChan; f++) + cout[f] += w * gp[f]; + } + + w = co[si[0]]; + ww += w * w; /* Sum of weights squared */ + gp += p->dinc[si[0]]; /* Far corner from base of cell */ + for (f = 0; f < p->outputChan; f++) + cout[f] += w * gp[f]; + + /* We distribute the correction needed in proportion to the */ + /* interpolation weighting, so the biggest correction is to the */ + /* closest vertex. */ + for (f = 0; f < p->outputChan; f++) + cout[f] = (out[f] - cout[f])/ww; /* Amount to distribute */ + + gp = ogp; + w = 1.0 - co[si[p->inputChan-1]]; /* Vertex at base of cell */ + for (f = 0; f < p->outputChan; f++) { + gp[f] += w * cout[f]; /* Apply correction */ + if (gp[f] < 0.0) { + gp[f] = 0.0; + rv |= 2; + } else if (gp[f] > 1.0) { + gp[f] = 1.0; + rv |= 2; + } + } + + for (e = p->inputChan-1; e > 0; e--) { /* Middle verticies */ + w = co[si[e]] - co[si[e-1]]; + gp += p->dinc[si[e]]; /* Move to top of cell in next largest dimension */ + for (f = 0; f < p->outputChan; f++) { + gp[f] += w * cout[f]; /* Apply correction */ + if (gp[f] < 0.0) { + gp[f] = 0.0; + rv |= 2; + } else if (gp[f] > 1.0) { + gp[f] = 1.0; + rv |= 2; + } + } + } + + w = co[si[0]]; + gp += p->dinc[si[0]]; /* Far corner from base of cell */ + for (f = 0; f < p->outputChan; f++) { + gp[f] += w * cout[f]; /* Apply correction */ + if (gp[f] < 0.0) { + gp[f] = 0.0; + rv |= 2; + } else if (gp[f] > 1.0) { + gp[f] = 1.0; + rv |= 2; + } + } + } + return rv; +} + + +/* ----------------------------------------------- */ /* Pseudo - Hilbert count sequencer */ /* This multi-dimensional count sequence is a distributed */ @@ -5550,35 +5811,44 @@ static int getNormFunc( #define CLIP_MARGIN 0.005 /* Margine to allow before reporting clipping = 0.5% */ +/* NOTE that ICM_CLUT_SET_FILTER turns out to be not very useful, */ +/* as it can result in reversals. Could #ifdef out the code ?? */ + /* Helper function to set multiple Lut tables simultaneously. */ /* Note that these tables all have to be compatible in */ /* having the same configuration and resolution. */ /* Set errc and return error number in underlying icc */ -/* Set warnc if there is clipping in the output values */ -/* 1 = input table, 2 = main clut, 3 = clut midpoin, 4 = midpoint interp, 5 = output table */ +/* Set warnc if there is clipping in the output values: */ +/* 1 = input table, 2 = main clut, 3 = clut midpoint, 4 = midpoint interp, 5 = output table */ +/* Note that clutfunc in[] value has "index under", ie: */ +/* at ((int *)in)[-chan-1], and for primary grid is simply the */ +/* grid index (ie. 5,3,8), and for the center of cells grid, is */ +/* the -index-1, ie. -6,-3,-8 */ int icmSetMultiLutTables( - int ntables, /* Number of tables to be set, 1..n */ - icmLut **pp, /* Pointer to array of Lut objects */ - int flags, /* Setting flags */ - void *cbctx, /* Opaque callback context pointer value */ - icColorSpaceSignature insig, /* Input color space */ - icColorSpaceSignature outsig, /* Output color space */ + int ntables, /* Number of tables to be set, 1..n */ + icmLut **pp, /* Pointer to array of Lut objects */ + int flags, /* Setting flags */ + void *cbctx, /* Opaque callback context pointer value */ + icColorSpaceSignature insig, /* Input color space */ + icColorSpaceSignature outsig, /* Output color space */ void (*infunc)(void *cbctx, double *out, double *in), /* Input transfer function, inspace->inspace' (NULL = default) */ /* Will be called ntables times for each input grid value */ - double *inmin, double *inmax, /* Maximum range of inspace' values */ - /* (NULL = default) */ + double *inmin, double *inmax, /* Maximum range of inspace' values */ + /* (NULL = default) */ void (*clutfunc)(void *cbntx, double *out, double *in), /* inspace' -> outspace[ntables]' transfer function */ /* will be called once for each input' grid value, and */ /* ntables output values should be written consecutively */ /* to out[]. */ - double *clutmin, double *clutmax, /* Maximum range of outspace' values */ - /* (NULL = default) */ - void (*outfunc)(void *cbntx, double *out, double *in)) + double *clutmin, double *clutmax, /* Maximum range of outspace' values */ + /* (NULL = default) */ + void (*outfunc)(void *cbntx, double *out, double *in), /* Output transfer function, outspace'->outspace (NULL = deflt) */ /* Will be called ntables times on each output value */ -{ + int *apxls_gmin, int *apxls_gmax /* If not NULL, the grid indexes not to be affected */ + /* by ICM_CLUT_SET_APXLS, defaulting to 0..>clutPoints-1 */ +) { icmLut *p, *pn; /* Pointer to 0'th nd tn'th Lut object */ icc *icp; /* Pointer to common icc */ int tn; @@ -5594,6 +5864,7 @@ int icmSetMultiLutTables( double *_iv, *iv, *ivn; /* Real index value/table value */ double imin[MAX_CHAN], imax[MAX_CHAN]; double omin[MAX_CHAN], omax[MAX_CHAN]; + int def_apxls_gmin[MAX_CHAN], def_apxls_gmax[MAX_CHAN]; void (*ifromindex)(double *out, double *in); /* Index to input color space function */ void (*itoentry)(double *out, double *in); /* Input color space to entry function */ void (*ifromentry)(double *out, double *in); /* Entry to input color space function */ @@ -5806,6 +6077,17 @@ int icmSetMultiLutTables( /* Allocate space for cell center value lookup */ if (flags & ICM_CLUT_SET_APXLS) { + if (apxls_gmin == NULL) { + apxls_gmin = def_apxls_gmin; + for (e = 0; e < p->inputChan; e++) + apxls_gmin[e] = 0; + } + if (apxls_gmax == NULL) { + apxls_gmax = def_apxls_gmax; + for (e = 0; e < p->inputChan; e++) + apxls_gmax[e] = p->clutPoints-1; + } + if ((clutTable2 = (double **) icp->al->calloc(icp->al,sizeof(double *), ntables)) == NULL) { sprintf(icp->err,"icmLut_set_tables malloc of cube center array failed"); icp->al->free(icp->al, _iv); @@ -5884,7 +6166,7 @@ int icmSetMultiLutTables( ti += ii[e] * p->dinc[e]; /* Clut index */ iv[e] = ii[e]/(p->clutPoints-1.0); /* Vertex coordinates */ iv[e] = iv[e] * (imax[e] - imin[e]) + imin[e]; /* Undo expansion to 0.0 - 1.0 */ - *((int *)&iv[-((int)e)-1]) = ii[e]; /* Trick to supply grid index in iv[] */ + *((int *)&iv[-((int)e)-1]) = ii[e]; /* Trick to supply grid index in iv[] */ } if (flags & ICM_CLUT_SET_FILTER) { @@ -5941,12 +6223,13 @@ int icmSetMultiLutTables( if (clutTable2 != NULL) { for (e = 0; e < p->inputChan; e++) { - if (ii[e] >= (p->clutPoints-1)) - break; /* Don't lookup beyond last */ + if (ii[e] < apxls_gmin[e] + || ii[e] >= apxls_gmax[e]) + break; /* Don't lookup outside least squares area */ iv[e] = (ii[e] + 0.5)/(p->clutPoints-1.0); /* Vertex coordinates + 0.5 */ iv[e] = iv[e] * (imax[e] - imin[e]) + imin[e]; /* Undo expansion to 0.0 - 1.0 */ - *((int *)&iv[-((int)e)-1]) = ii[e]; /* Trick to supply grid index in iv[] */ - /* (Not this is only the base for +0.5) */ + *((int *)&iv[-((int)e)-1]) = -ii[e]-1; /* Trick to supply -ve grid index in iv[] */ + /* (Not this is only the base for +0.5 center) */ } if (e >= p->inputChan) { /* We're not on the last row */ @@ -5995,7 +6278,7 @@ int icmSetMultiLutTables( #define APXLS_DIFF_THRHESH 0.2 /* Deal with cell center value, aproximate least squares adjustment. */ /* Subtract some of the mean of the surrounding center values from each grid value. */ - /* Skip the edges so that things like the white point are not changed. */ + /* Skip the range edges so that things like the white point or Video sync are not changed. */ /* Avoid modifying the value if the difference between the */ /* interpolated value and the current value is too great, */ /* and there is the possibility of different color aliases. */ @@ -6005,9 +6288,9 @@ int icmSetMultiLutTables( int ee; double cw = 1.0/(double)(1 << p->inputChan); /* Weight for each cube corner */ - /* For each cell center point except last row */ + /* For each cell center point except last row because we access ii[e]+1 */ for (e = 0; e < p->inputChan; e++) - ii[e] = 0; /* init coords */ + ii[e] = apxls_gmin[e]; /* init coords */ /* Compute linear interpolated value from center values */ for (ee = 0; ee < p->inputChan;) { @@ -6051,14 +6334,15 @@ int icmSetMultiLutTables( } pn->clutTable[ti + f] = vv; } + DBGSL(("nix %s apxls ov %s\n",icmPiv(p->inputChan, ii), icmPdv(pn->outputChan, ivn))); } } /* Increment coord */ for (ee = 0; ee < p->inputChan; ee++) { - if (++ii[ee] < (p->clutPoints-2)) /* Don't go through upper edge */ + if (++ii[ee] < (apxls_gmax[ee]-1)) /* Stop short of upper row of clutTable2 */ break; /* No carry */ - ii[ee] = 0; + ii[ee] = apxls_gmin[ee]; } } @@ -6069,6 +6353,7 @@ int icmSetMultiLutTables( } /* Apply any smoothing in the clipped region to the resulting clutTable */ + /* !!! should avoid smoothing outside apxls_gmin[e] & apxls_gmax[e] region !!! */ if (clutTable3 != NULL) { double *clutTable1; /* Copy of current unfilted values */ FCOUNT(cc, MAX_CHAN, p->inputChan); /* Surrounding counter */ @@ -6180,7 +6465,7 @@ int icmSetMultiLutTables( free(clutTable3); } - /* Create the output table entry values */ + /* Create the 1D output table entry values */ for (tn = 0; tn < ntables; tn++) { pn = pp[tn]; for (n = 0; n < pn->outputEnt; n++) { @@ -6239,19 +6524,21 @@ int icmSetMultiLutTables( /* Set errc and return error number */ /* Set warnc if there is clipping in the output values */ static int icmLut_set_tables ( -icmLut *p, /* Pointer to Lut object */ -int flags, /* Setting flags */ -void *cbctx, /* Opaque callback context pointer value */ -icColorSpaceSignature insig, /* Input color space */ -icColorSpaceSignature outsig, /* Output color space */ +icmLut *p, /* Pointer to Lut object */ +int flags, /* Setting flags */ +void *cbctx, /* Opaque callback context pointer value */ +icColorSpaceSignature insig, /* Input color space */ +icColorSpaceSignature outsig, /* Output color space */ void (*infunc)(void *cbcntx, double *out, double *in), /* Input transfer function, inspace->inspace' (NULL = default) */ -double *inmin, double *inmax, /* Maximum range of inspace' values (NULL = default) */ +double *inmin, double *inmax, /* Maximum range of inspace' values (NULL = default) */ void (*clutfunc)(void *cbctx, double *out, double *in), /* inspace' -> outspace' transfer function */ -double *clutmin, double *clutmax, /* Maximum range of outspace' values (NULL = default) */ -void (*outfunc)(void *cbctx, double *out, double *in) - /* Output transfer function, outspace'->outspace (NULL = deflt) */ +double *clutmin, double *clutmax, /* Maximum range of outspace' values (NULL = default) */ +void (*outfunc)(void *cbctx, double *out, double *in), + /* Output transfer function, outspace'->outspace (NULL = deflt) */ +int *apxls_gmin, int *apxls_gmax /* If not NULL, the grid indexes not to be affected */ + /* by ICM_CLUT_SET_APXLS, defaulting to 0..>clutPoints-1 */ ) { struct _icmLut *pp[3]; @@ -6263,7 +6550,8 @@ void (*outfunc)(void *cbctx, double *out, double *in) inmin, inmax, clutfunc, clutmin, clutmax, - outfunc); + outfunc, + apxls_gmin, apxls_gmax); } /* - - - - - - - - - - - - - - - - */ @@ -6347,21 +6635,8 @@ static int icmLut_read( p->outputChan = read_UInt8Number(bp+9); p->clutPoints = read_UInt8Number(bp+10); - /* Sanity check */ - if (p->inputChan < 1) { - sprintf(icp->err,"icmLut_read: No input channels!"); - return icp->errc = 1; - } - - if (p->inputChan > MAX_CHAN) { - sprintf(icp->err,"icmLut_read: Can't handle > %d input channels\n",MAX_CHAN); - return icp->errc = 1; - } - - if (p->outputChan > MAX_CHAN) { - sprintf(icp->err,"icmLut_read: Can't handle > %d output channels\n",MAX_CHAN); - return icp->errc = 1; - } + if (icp->allowclutPoints256 && p->clutPoints == 0) /* Special case */ + p->clutPoints = 256; /* Read 3x3 transform matrix */ for (j = 0; j < 3; j++) { /* Rows */ @@ -6380,7 +6655,7 @@ static int icmLut_read( bp = buf+52; } - /* Sanity check dimensions. This protects against */ + /* Sanity check tag size. This protects against */ /* subsequent integer overflows involving the dimensions. */ if ((size = icmLut_get_size((icmBase *)p)) == UINT_MAX || size > len) { @@ -6389,12 +6664,15 @@ static int icmLut_read( return icp->errc = 1; } - /* Read the input tables */ - size = (p->inputChan * p->inputEnt); + /* Sanity check the dimensions and resolution values agains limits, */ + /* allocate space for them and generate internal offset tables. */ if ((rv = p->allocate((icmBase *)p)) != 0) { icp->al->free(icp->al, buf); return rv; } + + /* Read the input tables */ + size = (p->inputChan * p->inputEnt); if (p->ttype == icSigLut8Type) { for (i = 0; i < size; i++, bp += 1) p->inputTable[i] = read_DCS8Number(bp); @@ -6405,10 +6683,6 @@ static int icmLut_read( /* Read the clut table */ size = (p->outputChan * sat_pow(p->clutPoints,p->inputChan)); - if ((rv = p->allocate((icmBase *)p)) != 0) { - icp->al->free(icp->al, buf); - return rv; - } if (p->ttype == icSigLut8Type) { for (i = 0; i < size; i++, bp += 1) p->clutTable[i] = read_DCS8Number(bp); @@ -6419,10 +6693,6 @@ static int icmLut_read( /* Read the output tables */ size = (p->outputChan * p->outputEnt); - if ((rv = p->allocate((icmBase *)p)) != 0) { - icp->al->free(icp->al, buf); - return rv; - } if (p->ttype == icSigLut8Type) { for (i = 0; i < size; i++, bp += 1) p->outputTable[i] = read_DCS8Number(bp); @@ -6431,20 +6701,6 @@ static int icmLut_read( p->outputTable[i] = read_DCS16Number(bp); } - /* Private: compute dimensional increment though clut */ - /* Note that first channel varies least rapidly. */ - i = p->inputChan-1; - p->dinc[i--] = p->outputChan; - for (; i < p->inputChan; i--) - p->dinc[i] = p->dinc[i+1] * p->clutPoints; - - /* Private: compute offsets from base of cube to other corners */ - for (p->dcube[0] = 0, g = 1, j = 0; j < p->inputChan; j++) { - for (i = 0; i < g; i++) - p->dcube[g+i] = p->dcube[i] + p->dinc[j]; - g *= 2; - } - icp->al->free(icp->al, buf); return 0; } @@ -6491,10 +6747,19 @@ static int icmLut_write( icp->al->free(icp->al, buf); return icp->errc = rv; } - if ((rv = write_UInt8Number(p->clutPoints, bp+10)) != 0) { - sprintf(icp->err,"icmLut_write: write_UInt8Number() failed"); - icp->al->free(icp->al, buf); - return icp->errc = rv; + + if (icp->allowclutPoints256 && p->clutPoints == 256) { + if ((rv = write_UInt8Number(0, bp+10)) != 0) { + sprintf(icp->err,"icmLut_write: write_UInt8Number() failed"); + icp->al->free(icp->al, buf); + return icp->errc = rv; + } + } else { + if ((rv = write_UInt8Number(p->clutPoints, bp+10)) != 0) { + sprintf(icp->err,"icmLut_write: write_UInt8Number() failed"); + icp->al->free(icp->al, buf); + return icp->errc = rv; + } } write_UInt8Number(0, bp+11); /* Set padding to 0 */ @@ -6683,7 +6948,9 @@ static void icmLut_dump( } } -/* Allocate variable sized data elements */ +/* Sanity check the input & output dimensions, and */ +/* allocate variable sized data elements, and */ +/* generate internal dimension offset tables */ static int icmLut_allocate( icmBase *pp ) { @@ -6691,7 +6958,7 @@ static int icmLut_allocate( icmLut *p = (icmLut *)pp; icc *icp = p->icp; - /* Sanity check */ + /* Sanity check, so that dinc[] comp. won't fail */ if (p->inputChan < 1) { sprintf(icp->err,"icmLut_alloc: Can't handle %d input channels\n",p->inputChan); return icp->errc = 1; @@ -6805,6 +7072,8 @@ static icmBase *new_icmLut( icmLut *p; if ((p = (icmLut *) icp->al->calloc(icp->al,1,sizeof(icmLut))) == NULL) return NULL; + p->icp = icp; + p->ttype = icSigLut16Type; p->refcount = 1; p->get_size = icmLut_get_size; @@ -6825,8 +7094,7 @@ static icmBase *new_icmLut( /* Set method */ p->set_tables = icmLut_set_tables; - - p->icp = icp; + p->tune_value = icmLut_tune_value_sx; /* Default to most likely simplex */ /* Set matrix to reasonable default */ for (i = 0; i < 3; i++) @@ -7946,7 +8214,7 @@ static int icmColorantTable_allocate( if (p->count != p->_count) { unsigned int i; if (ovr_mul(p->count, sizeof(icmColorantTableVal))) { - sprintf(icp->err,"icmColorantTable_alloc: count overflow (%d of %ld bytes)",p->count,sizeof(icmColorantTableVal)); + sprintf(icp->err,"icmColorantTable_alloc: count overflow (%d of %lu bytes)",p->count,sizeof(icmColorantTableVal)); return icp->errc = 1; } if (p->data != NULL) @@ -10789,10 +11057,18 @@ static int icmHeader_read( p->cmmId = read_SInt32Number(buf + 4); /* CMM for profile */ tt = read_UInt8Number(buf + 8); /* Raw major version number */ p->majv = (tt >> 4) * 10 + (tt & 0xf); /* Integer major version number */ - icp->ver = p->majv > 3 ? 1 : 0; /* Set major version flag in icc */ - tt = read_UInt8Number(buf + 9); /* Raw minor/bug fix version numbers */ + tt = read_UInt8Number(buf + 9); /* Raw minor & bug fix version numbers */ p->minv = (tt >> 4); /* Integer minor version number */ p->bfv = (tt & 0xf); /* Integer bug fix version number */ + if (p->majv < 3) { /* Set version class */ + if (p->minv >= 4) + icp->ver = icmVersion2_4; + else if (p->minv >= 3) + icp->ver = icmVersion2_3; + else + icp->ver = icmVersionDefault; + } else + icp->ver = icmVersion4_1; p->deviceClass = (icProfileClassSignature) read_SInt32Number(buf + 12); /* Type of profile */ p->colorSpace = (icColorSpaceSignature) @@ -10819,13 +11095,13 @@ static int icmHeader_read( } p->creator = read_SInt32Number(buf + 80); /* Profile creator */ - for (tt = 0; tt < 16; tt++) - p->id[tt] = icp->ver ? read_UInt8Number(buf + 84 + tt) : 0; /* Profile ID */ + for (tt = 0; tt < 16; tt++) /* Profile ID */ + p->id[tt] = icp->ver >= icmVersion4_1 ? read_UInt8Number(buf + 84 + tt) : 0; icp->al->free(icp->al, buf); #ifndef ENABLE_V4 - if (icp->ver) { + if (icp->ver >= icmVersion4_1) { sprintf(icp->err,"icmHeader_read: ICC V4 not supported!"); return icp->errc = 1; } @@ -10874,6 +11150,7 @@ static int icmHeader_write( icp->al->free(icp->al, buf); return icp->errc = 1; } + // ~~~ Hmm. We're not checking ->ver is >= corresponding header version number ~~ tt = ((p->majv/10) << 4) + (p->majv % 10); if ((rv = write_UInt8Number(tt, buf + 8)) != 0) { /* Raw major version number */ sprintf(icp->err,"icmHeader_write: Uint8Number major version"); @@ -10951,7 +11228,7 @@ static int icmHeader_write( icp->al->free(icp->al, buf); return icp->errc = rv; } - if (doid == 0 && icp->ver) { /* ID is V4.0+ feature */ + if (doid == 0 && icp->ver >= icmVersion4_1) { /* ID is V4.0+ feature */ for (tt = 0; tt < 16; tt++) { if ((rv = write_UInt8Number(p->id[tt], buf + 84 + tt)) != 0) { /* Profile ID */ sprintf(icp->err,"icmHeader_write: UInt8Number creator"); @@ -10997,7 +11274,7 @@ static void icmHeader_dump( op->gprintf(op," Rndrng Intnt = %s\n", string_RenderingIntent(p->renderingIntent)); op->gprintf(op," Illuminant = %s\n", string_XYZNumber_and_Lab(&p->illuminant)); op->gprintf(op," Creator = %s\n", tag2str(p->creator)); /* ~~~ */ - if (p->icp->ver) { /* V4.0+ feature */ + if (p->icp->ver >= icmVersion4_1) { /* V4.0+ feature */ for (i = 0; i < 16; i++) { /* Check if ID has been set */ if (p->id[i] != 0) break; @@ -11374,10 +11651,15 @@ static int check_icc_legal( /* Found entry, so now check that all the required tags are present. */ for (j = 0; tagchecktable[i].tags[j] != icMaxEnumType; j++) { if (p->find_tag(p, tagchecktable[i].tags[j]) != 0) { /* Not present! */ +#ifdef NEVER + printf("icc_check_legal: deviceClass %s is missing required tag %s", tag2str(sig), tag2str(tagchecktable[i].tags[j])); +#endif if (tagchecktable[i].chans == -200 || tagchecktable[i].chans == -dchans) { /* But can try next table */ break; } + /* ~~99 Hmm. Should report all possible missing tags from */ + /* previous failed tables ~~~999 */ sprintf(p->err,"icc_check_legal: deviceClass %s is missing required tag %s", tag2str(sig), tag2str(tagchecktable[i].tags[j])); return p->errc = 1; @@ -11781,7 +12063,7 @@ static int icc_write_x( /* If V4.0+, Compute the MD5 id for the profile. */ /* We do this by writing to a fake icmFile */ - if (p->ver) { + if (p->ver >= icmVersion4_1) { icmMD5 *md5 = NULL; icmFile *ofp, *dfp = NULL; @@ -11801,13 +12083,13 @@ static int icc_write_x( ofp = p->fp; p->fp = dfp; - /* Dumy write the header */ + /* Dummy write the header */ if ((rv = p->header->write(p->header, 0, 1)) != 0) { p->al->free(p->al, buf); return rv; } - /* Dumy write the tag table */ + /* Dummy write the tag table */ if ( p->fp->seek(p->fp, 128) != 0 || p->fp->write(p->fp, buf, 1, len) != len) { sprintf(p->err,"icc_write: seek() or write() failed"); @@ -11815,7 +12097,7 @@ static int icc_write_x( return p->errc = 1; } - /* Dumy write all the tag element data */ + /* Dummy write all the tag element data */ /* (We invert meaning of touched here) */ for (i = 0; i < p->count; i++) { /* For all the tag element data */ if (p->data[i].objp->touched == 0) @@ -12202,8 +12484,8 @@ static icmBase *icc_read_tag_ix( } /* Read the tag element data of the first matching, and return a pointer to the object */ -/* Returns NULL if error - icc->errc will contain: */ -/* 2 if not found */ +/* Returns NULL if error - icc->errc will contain: */ +/* 2 if not found */ /* Doesn't read uknown type tags */ static icmBase *icc_read_tag( icc *p, @@ -12227,7 +12509,7 @@ static icmBase *icc_read_tag( } /* Read the tag element data of the first matching, and return a pointer to the object */ -/* Returns NULL if error. +/* Returns NULL if error. */ /* Returns an icmSigUnknownType object if the tag type isn't handled by a specific object. */ /* NOTE: we don't handle tag duplication - you'll always get the first in the file. */ static icmBase *icc_read_tag_any( @@ -12656,6 +12938,25 @@ static void Lut_Luv2Lut(double *out, double *in) { } /* - - - - - - - - - - - - - - - - */ +/* Convert YCbCr to Lut number */ +/* We are assuming full range here. foot/head scaling */ +/* should be done outside the icc profile. */ + +/* Convert Lut table index/value to YCbCr */ +static void Lut_Lut2YCbCr(double *out, double *in) { + out[0] = in[0]; /* Y */ + out[1] = in[1] - 0.5; /* Cb */ + out[2] = in[2] - 0.5; /* Cr */ +} + +/* Convert YCbCr to Lut table index/value */ +static void Lut_YCbCr2Lut(double *out, double *in) { + out[0] = in[0]; /* Y */ + out[1] = in[1] + 0.5; /* Cb */ + out[2] = in[2] + 0.5; /* Cr */ +} + +/* - - - - - - - - - - - - - - - - */ /* Default N component conversions */ static void Lut_N(double *out, double *in, int nc) { for (--nc; nc >= 0; nc--) @@ -12754,7 +13055,6 @@ static void Lut_15(double *out, double *in) { /* Function table - match conversions to color spaces. */ /* Anything not here, we don't know how to convert. */ -/* (ie. YCbCr) */ static struct { icColorSpaceSignature csig; void (*fromLut)(double *out, double *in); /* from Lut index/entry */ @@ -12769,6 +13069,7 @@ static struct { {icmSigLV2Data, Lut_Lut2LV2_16, Lut_L2LutV2_16 }, {icmSigLV4Data, Lut_Lut2LV4_16, Lut_L2LutV4_16 }, {icSigLuvData, Lut_Lut2Luv, Lut_Luv2Lut }, + {icSigYCbCrData, Lut_Lut2YCbCr, Lut_YCbCr2Lut }, {icSigYxyData, Lut_3, Lut_3 }, {icSigRgbData, Lut_3, Lut_3 }, {icSigGrayData, Lut_1, Lut_1 }, @@ -12842,7 +13143,7 @@ static int getNormFunc( if (tagType == icSigLut16Type) /* Lut16 retains legacy encoding */ csig = icmSigLabV2Data; else { /* Other tag types use version specific encoding */ - if (icp->ver) + if (icp->ver >= icmVersion4_1) csig = icmSigLabV4Data; else csig = icmSigLabV2Data; @@ -12875,7 +13176,6 @@ static int getNormFunc( /* Function table - match ranges to color spaces. */ /* Anything not here, we don't know how to convert. */ -/* (ie. YCbCr) */ /* Hmm. we're not handling Lab8 properly ?? ~~~8888 */ static struct { icColorSpaceSignature csig; @@ -12894,7 +13194,7 @@ static struct { {icmSigLV4Data, 1, { 0.0 }, { 100.0 } }, {icSigLuvData, 0, { 0.0, -128.0, -128.0 }, { 100.0, 127.0 + 255.0/256.0, 127.0 + 255.0/256.0 } }, - {icSigYCbCrData, 1, { 0.0 }, { 1.0 } }, /* ??? */ + {icSigYCbCrData, 0, { 0.0, -0.5, -0.5 }, { 1.0, 0.5, 0.5 } }, /* Full range */ {icSigYxyData, 1, { 0.0 }, { 1.0 } }, /* ??? */ {icSigRgbData, 1, { 0.0 }, { 1.0 } }, {icSigGrayData, 1, { 0.0 }, { 1.0 } }, @@ -12945,7 +13245,7 @@ static int getRange( if (tagType == icSigLut16Type) /* Lut16 retains legacy encoding */ csig = icmSigLabV2Data; else { /* Other tag types use version specific encoding */ - if (icp->ver) + if (icp->ver >= icmVersion4_1) csig = icmSigLabV4Data; else csig = icmSigLabV2Data; @@ -12979,8 +13279,8 @@ static int getRange( return 0; } -/* ============================================= */ -/* Misc. support functions. */ +/* =============================================================== */ +/* Misc. support functions. */ /* Clamp a 3 vector to be +ve */ void icmClamp3(double out[3], double in[3]) { @@ -13003,6 +13303,20 @@ void icmSub3(double out[3], double in1[3], double in2[3]) { out[2] = in1[2] - in2[2]; } +/* Divide two 3 vectors, out = in1/in2 */ +void icmDiv3(double out[3], double in1[3], double in2[3]) { + out[0] = in1[0]/in2[0]; + out[1] = in1[1]/in2[1]; + out[2] = in1[2]/in2[2]; +} + +/* Multiply two 3 vectors, out = in1 * in2 */ +void icmMul3(double out[3], double in1[3], double in2[3]) { + out[0] = in1[0] * in2[0]; + out[1] = in1[1] * in2[1]; + out[2] = in1[2] * in2[2]; +} + /* - - - - - - - - - - - - - - - - - - - - - - - - */ /* Set a 3x3 matrix to unity */ @@ -13312,6 +13626,18 @@ void icmBlend3(double out[3], double in0[3], double in1[3], double bf) { out[2] = (1.0 - bf) * in0[2] + bf * in1[2]; } +/* Clip a vector to the range 0.0 .. 1.0 */ +void icmClip3(double out[3], double in[3]) { + int j; + for (j = 0; j < 3; j++) { + out[j] = in[j]; + if (out[j] < 0.0) + out[j] = 0.0; + else if (out[j] > 1.0) + out[j] = 1.0; + } +} + /* Normalise a 3 vector to the given length. Return nz if not normalisable */ int icmNormalize3(double out[3], double in[3], double len) { double tt = sqrt(in[0] * in[0] + in[1] * in[1] + in[2] * in[2]); @@ -13686,6 +14012,74 @@ double icmPlaneDist3(double eq[4], double p[3]) { } /* - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Given 2 2D points, compute a plane equation. */ +/* The normal will be right handed given the order of the points */ +/* The plane equation will be the 2 normal components and the constant. */ +/* Return nz if any points are cooincident or co-linear */ +int icmPlaneEqn2(double eq[3], double p0[2], double p1[2]) { + double ll, v1[3]; + + /* Compute vectors along edge */ + v1[0] = p1[0] - p0[0]; + v1[1] = p1[1] - p0[1]; + + /* Normal to vector */ + eq[0] = v1[1]; + eq[1] = -v1[0]; + + /* Normalise the equation */ + ll = sqrt(eq[0] * eq[0] + eq[1] * eq[1]); + if (ll < 1e-10) { + return 1; + } + eq[0] /= ll; + eq[1] /= ll; + + /* Compute the plane equation constant */ + eq[2] = - (eq[0] * p0[0]) + - (eq[1] * p0[1]); + + return 0; +} + +/* Given a 2D point and a plane equation, return the signed */ +/* distance from the plane. The distance will be +ve if the point */ +/* is to the right of the plane formed by two points in order */ +double icmPlaneDist2(double eq[3], double p[2]) { + double rv; + + rv = eq[0] * p[0] + + eq[1] * p[1] + + eq[2]; + + return rv; +} + +/* Given two infinite 2D lines define by 4 points, compute the intersection. */ +/* Return nz if there is no intersection (lines are parallel) */ +int icmLineIntersect2(double res[2], double p1[2], double p2[2], double p3[2], double p4[2]) { + /* Compute by determinants */ + double x1y2_y1x2 = p1[0] * p2[1] - p1[1] * p2[0]; + double x3y4_y3x4 = p3[0] * p4[1] - p3[1] * p4[0]; + double x1_x2 = p1[0] - p2[0]; + double y1_y2 = p1[1] - p2[1]; + double x3_x4 = p3[0] - p4[0]; + double y3_y4 = p3[1] - p4[1]; + double num; /* Numerator */ + + num = x1_x2 * y3_y4 - y1_y2 * x3_x4; + + if (fabs(num) < 1e-10) + return 1; + + res[0] = (x1y2_y1x2 * x3_x4 - x1_x2 * x3y4_y3x4)/num; + res[1] = (x1y2_y1x2 * y3_y4 - y1_y2 * x3y4_y3x4)/num; + + return 0; +} + + +/* - - - - - - - - - - - - - - - - - - - - - - - - */ /* CIE Y (range 0 .. 1) to perceptual CIE 1976 L* (range 0 .. 100) */ double icmY2L(double val) { @@ -13783,7 +14177,7 @@ void icmLCh2Lab(double *out, double *in) { out[2] = C * sin(h); } -/* Lab to LCh */ +/* Lab to LCh (general to polar, works with Luv too) */ void icmLab2LCh(double *out, double *in) { double C, h; @@ -13888,9 +14282,6 @@ extern ICCLIB_API void icmLuv2XYZ(icmXYZNumber *w, double *out, double *in) { out[2] = Z; } -/* NOTE :- none of the following seven have been protected */ -/* against arithmmetic issues (ie. for black) */ - /* CIE XYZ to perceptual CIE 1976 UCS diagram Yu'v'*/ /* (Yu'v' is a better chromaticity space than Yxy) */ extern ICCLIB_API void icmXYZ21976UCS(double *out, double *in) { @@ -13898,8 +14289,15 @@ extern ICCLIB_API void icmXYZ21976UCS(double *out, double *in) { double den, u, v; den = (X + 15.0 * Y + 3.0 * Z); - u = (4.0 * X) / den; - v = (9.0 * Y) / den; + + if (den < 1e-9) { + Y = 0.0; + u = 4.0/19.0; + v = 9.0/19.0; + } else { + u = (4.0 * X) / den; + v = (9.0 * Y) / den; + } out[0] = Y; out[1] = u; @@ -13914,8 +14312,12 @@ extern ICCLIB_API void icm1976UCS2XYZ(double *out, double *in) { u = in[1]; v = in[2]; - X = ((9.0 * u * Y)/(4.0 * v)); - Z = -(((20.0 * v + 3.0 * u - 12.0) * Y)/(4.0 * v)); + if (v < 1e-9) { + X = Y = Z = 0.0; + } else { + X = ((9.0 * u * Y)/(4.0 * v)); + Z = -(((20.0 * v + 3.0 * u - 12.0) * Y)/(4.0 * v)); + } out[0] = X; out[1] = Y; @@ -13927,10 +14329,18 @@ extern ICCLIB_API void icm1976UCS2XYZ(double *out, double *in) { /* in computing color temperatures.) */ extern ICCLIB_API void icmXYZ21960UCS(double *out, double *in) { double X = in[0], Y = in[1], Z = in[2]; - double u, v; + double den, u, v; - u = (4.0 * X) / (X + 15.0 * Y + 3.0 * Z); - v = (6.0 * Y) / (X + 15.0 * Y + 3.0 * Z); + den = (X + 15.0 * Y + 3.0 * Z); + + if (den < 1e-9) { + Y = 0.0; + u = 4.0/19.0; + v = 6.0/19.0; + } else { + u = (4.0 * X) / den; + v = (6.0 * Y) / den; + } out[0] = Y; out[1] = u; @@ -13945,8 +14355,12 @@ extern ICCLIB_API void icm1960UCS2XYZ(double *out, double *in) { u = in[1]; v = in[2]; - X = ((3.0 * u * Y)/(2.0 * v)); - Z = -(((10.0 * v + u - 4.0) * Y)/(2.0 * v)); + if (v < 1e-9) { + X = Y = Z = 0.0; + } else { + X = ((3.0 * u * Y)/(2.0 * v)); + Z = -(((10.0 * v + u - 4.0) * Y)/(2.0 * v)); + } out[0] = X; out[1] = Y; @@ -13994,6 +14408,7 @@ extern ICCLIB_API void icm1964WUV2XYZ(icmXYZNumber *w, double *out, double *in) } /* CIE CIE1960 UCS to perceptual CIE 1964 WUV (U*V*W*) */ +/* (This is used in computing CRI) */ extern ICCLIB_API void icm1960UCS21964WUV(icmXYZNumber *w, double *out, double *in) { double W, U, V; double wucs[3]; @@ -14011,6 +14426,7 @@ extern ICCLIB_API void icm1960UCS21964WUV(icmXYZNumber *w, double *out, double * } /* - - - - - - - - - - - - - - - - - - - - - - - - */ +/* NOTE :- that these values are for the 1931 standard observer */ /* available D50 Illuminant */ icmXYZNumber icmD50 = { /* Profile illuminant - D50 */ @@ -14046,6 +14462,7 @@ double icmD65_100_ary3[3] = { /* Profile illuminant - D65, scaled to 100 */ 95.05, 100.00, 108.90 }; + /* Default black point */ icmXYZNumber icmBlack = { 0.0000, 0.0000, 0.0000 @@ -14131,7 +14548,7 @@ double icmCIE94sq(double Lab0[3], double Lab1[3]) { c12 = sqrt(c1 * c2); /* Symetric chromanance */ /* delta chromanance squared */ - dc = c2 - c1; + dc = c1 - c2; dcsq = dc * dc; } @@ -14142,8 +14559,8 @@ double icmCIE94sq(double Lab0[3], double Lab1[3]) { double sc, sh; /* Weighting factors for delta chromanance & delta hue */ - sc = 1.0 + 0.048 * c12; - sh = 1.0 + 0.014 * c12; + sc = 1.0 + 0.045 * c12; + sh = 1.0 + 0.015 * c12; return dlsq + dcsq/(sc * sc) + dhsq/(sh * sh); } } @@ -14291,7 +14708,7 @@ double icmCIE2K(double *Lab0, double *Lab1) { } /* Return the CIEDE2000 Delta E color difference measure for two XYZ values */ -extern ICCLIB_API double icmXYZCIE2K(icmXYZNumber *w, double *in0, double *in1) { +ICCLIB_API double icmXYZCIE2K(icmXYZNumber *w, double *in0, double *in1) { double lab0[3], lab1[3]; icmXYZ2Lab(w, lab0, in0); @@ -14300,12 +14717,15 @@ extern ICCLIB_API double icmXYZCIE2K(icmXYZNumber *w, double *in0, double *in1) } + /* - - - - - - - - - - - - - - - - - - - - - - - - */ /* Chromatic adaptation transform utility */ /* Return a 3x3 chromatic adaptation matrix */ /* Use icmMulBy3x3(dst, mat, src) */ +/* NOTE that to transform primaries they */ +/* must be mat[XYZ][RGB] format! */ void icmChromAdaptMatrix( - int flags, /* Use bradford, Transform given matrix flags */ + int flags, /* Use Bradford, Transform given matrix flags */ icmXYZNumber d_wp, /* Destination white point */ icmXYZNumber s_wp, /* Source white point */ double mat[3][3] /* Destination matrix */ @@ -14381,7 +14801,7 @@ int icmRGBprim2matrix( icmXYZNumber red, /* Red colorant */ icmXYZNumber green, /* Green colorant */ icmXYZNumber blue, /* Blue colorant */ - double mat[3][3] /* Destination matrix */ + double mat[3][3] /* Destination matrix[RGB][XYZ] */ ) { double tmat[3][3]; double t[3]; @@ -14414,19 +14834,73 @@ int icmRGBprim2matrix( /* Now formulate the transform matrix */ mat[0][0] = red.X * t[0]; - mat[0][1] = green.X * t[1]; - mat[0][2] = blue.X * t[2]; - mat[1][0] = red.Y * t[0]; + mat[1][0] = green.X * t[1]; + mat[2][0] = blue.X * t[2]; + mat[0][1] = red.Y * t[0]; mat[1][1] = green.Y * t[1]; - mat[1][2] = blue.Y * t[2]; - mat[2][0] = red.Z * t[0]; - mat[2][1] = green.Z * t[1]; + mat[2][1] = blue.Y * t[2]; + mat[0][2] = red.Z * t[0]; + mat[1][2] = green.Z * t[1]; mat[2][2] = blue.Z * t[2]; return 0; } /* - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Pre-round RGB device primary values to ensure that */ +/* the sum of the quantized primaries is the same as */ +/* the quantized sum. */ +void quantizeRGBprimsS15Fixed16( + double mat[3][3] /* matrix[RGB][XYZ] */ +) { + int i, j; + double sum[3]; + +// printf("D50 = %f %f %f\n",icmD50.X, icmD50.Y, icmD50.Z); + + /* Compute target sum of primary XYZ */ + for (i = 0; i < 3; i++) { + sum[i] = 0.0; + for (j = 0; j < 3; j++) + sum[i] += mat[j][i]; + } +// printf("Sum = %f %f %f\n",sum[0], sum[1], sum[2]); + + /* Pre-quantize the primary XYZ's, and then ensure that the */ + /* sum of the quantized values is the quantized sum by assigning */ + /* the rounding error to the largest component. */ + for (i = 0; i < 3; i++) { + int bix = 0; + double bval = -1e9; + + /* locate the largest and quantize each component */ + for (j = 0; j < 3; j++) { + if (fabs(mat[j][i]) > bval) { /* Locate largest */ + bix = j; + bval = fabs(mat[j][i]); + } + mat[j][i] = round_S15Fixed16Number(mat[j][i]); + } + + /* Compute the value the largest has to be */ + /* to ensure that sum of the quantized values is */ + /* equal to the quantized sum */ + for (j = 0; j < 3; j++) { + if (j == bix) + continue; + sum[i] -= mat[j][i]; + } + mat[bix][i] = round_S15Fixed16Number(sum[i]); + + /* Check the sum of the quantized values */ +// sum[i] = 0.0; +// for (j = 0; j < 3; j++) +// sum[i] += mat[j][i]; + } +// printf("Q Sum = %f %f %f\n",sum[0], sum[1], sum[2]); +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - */ /* Some PCS utility functions */ /* Clip Lab, while maintaining hue angle. */ @@ -14530,7 +15004,426 @@ int icmClipXYZ(double out[3], double in[3]) { return 1; } -/* - - - - - - - - - - - - - - - - - - - - - - - - */ +/* --------------------------------------------------------------- */ +/* Some video specific functions */ + +/* Convert Lut table index/value to YPbPr */ +/* (Same as Lut_Lut2YPbPr() ) */ +void icmLut2YPbPr(double *out, double *in) { + out[0] = in[0]; /* Y */ + out[1] = in[1] - 0.5; /* Cb */ + out[2] = in[2] - 0.5; /* Cr */ +} + +/* Convert YPbPr to Lut table index/value */ +/* (Same as Lut_YPbPr2Lut() ) */ +void icmYPbPr2Lut(double *out, double *in) { + out[0] = in[0]; /* Y */ + out[1] = in[1] + 0.5; /* Cb */ + out[2] = in[2] + 0.5; /* Cr */ +} + +/* Convert Rec601 RGB' into YPbPr, or "full range YCbCr" */ +/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */ +/* [From the Rec601 spec. ] */ +void icmRec601_RGBd_2_YPbPr(double out[3], double in[3]) { + double tt[3]; + + tt[0] = 0.299 * in[0] + 0.587 * in[1] + 0.114 * in[2]; + + tt[1] = -0.299 /1.772 * in[0] + + -0.587 /1.772 * in[1] + + (1.0-0.114)/1.772 * in[2]; + + tt[2] = (1.0-0.299)/1.402 * in[0] + + -0.587 /1.402 * in[1] + + -0.114 /1.402 * in[2]; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + +/* Convert Rec601 YPbPr to RGB' (== "full range YCbCr") */ +/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */ +/* [Inverse of above] */ +void icmRec601_YPbPr_2_RGBd(double out[3], double in[3]) { + double tt[3]; + + tt[0] = 1.000000000 * in[0] + 0.000000000 * in[1] + 1.402000000 * in[2]; + tt[1] = 1.000000000 * in[0] + -0.344136286 * in[1] + -0.714136286 * in[2]; + tt[2] = 1.000000000 * in[0] + 1.772000000 * in[1] + 0.000000000 * in[2]; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + + +/* Convert Rec709 1150/60/2:1 RGB' into YPbPr, or "full range YCbCr" */ +/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */ +/* [From the Rec709 specification] */ +void icmRec709_RGBd_2_YPbPr(double out[3], double in[3]) { + double tt[3]; + + tt[0] = 0.2126 * in[0] + 0.7152 * in[1] + 0.0722 * in[2]; + + tt[1] = 0.5389 * -0.2126 * in[0] + + 0.5389 * -0.7152 * in[1] + + 0.5389 * (1.0-0.0722) * in[2]; + + tt[2] = 0.6350 * (1.0-0.2126) * in[0] + + 0.6350 * -0.7152 * in[1] + + 0.6350 * -0.0722 * in[2]; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + +/* Convert Rec709 1150/60/2:1 YPbPr to RGB' (== "full range YCbCr") */ +/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */ +/* [Inverse of above] */ +void icmRec709_YPbPr_2_RGBd(double out[3], double in[3]) { + double tt[3]; + + tt[0] = 1.000000000 * in[0] + 0.000000000 * in[1] + 1.574803150 * in[2]; + tt[1] = 1.000000000 * in[0] + -0.187327487 * in[1] + -0.468125209 * in[2]; + tt[2] = 1.000000000 * in[0] + 1.855631843 * in[1] + 0.000000000 * in[2]; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + +/* Convert Rec709 1250/50/2:1 RGB' into YPbPr, or "full range YCbCr" */ +/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */ +/* [From the Rec709 specification] */ +void icmRec709_50_RGBd_2_YPbPr(double out[3], double in[3]) { + double tt[3]; + + tt[0] = 0.299 * in[0] + 0.587 * in[1] + 0.114 * in[2]; + + tt[1] = 0.564 * -0.299 * in[0] + + 0.564 * -0.587 * in[1] + + 0.564 * (1.0-0.114) * in[2]; + + tt[2] = 0.713 * (1.0-0.299) * in[0] + + 0.713 * -0.587 * in[1] + + 0.713 * -0.114 * in[2]; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + +/* Convert Rec709 1250/50/2:1 YPbPr to RGB' (== "full range YCbCr") */ +/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */ +/* [Inverse of above] */ +void icmRec709_50_YPbPr_2_RGBd(double out[3], double in[3]) { + double tt[3]; + + tt[0] = 1.000000000 * in[0] + 0.000000000 * in[1] + 1.402524544 * in[2]; + tt[1] = 1.000000000 * in[0] + -0.344340136 * in[1] + -0.714403473 * in[2]; + tt[2] = 1.000000000 * in[0] + 1.773049645 * in[1] + 0.000000000 * in[2]; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + + +/* Convert Rec2020 RGB' into Non-constant liminance YPbPr, or "full range YCbCr" */ +/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */ +/* [From the Rec2020 specification] */ +void icmRec2020_NCL_RGBd_2_YPbPr(double out[3], double in[3]) { + double tt[3]; + + tt[0] = 0.2627 * in[0] + 0.6780 * in[1] + 0.0593 * in[2]; + + tt[1] = 1/1.8814 * -0.2627 * in[0] + + 1/1.8814 * -0.6780 * in[1] + + 1/1.8814 * (1.0-0.0593) * in[2]; + + tt[2] = 1/1.4746 * (1.0-0.2627) * in[0] + + 1/1.4746 * -0.6780 * in[1] + + 1/1.4746 * -0.0593 * in[2]; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + +/* Convert Rec2020 Non-constant liminance YPbPr into RGB' (== "full range YCbCr") */ +/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */ +/* [Inverse of above] */ +void icmRec2020_NCL_YPbPr_2_RGBd(double out[3], double in[3]) { + double tt[3]; + + tt[0] = 1.000000000 * in[0] + 0.000000000 * in[1] + 1.474600000 * in[2]; + tt[1] = 1.000000000 * in[0] + -0.164553127 * in[1] + -0.571353127 * in[2]; + tt[2] = 1.000000000 * in[0] + 1.881400000 * in[1] + 0.000000000 * in[2]; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + +/* Convert Rec2020 RGB' into Constant liminance YPbPr, or "full range YCbCr" */ +/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */ +/* [From the Rec2020 specification] */ +void icmRec2020_CL_RGBd_2_YPbPr(double out[3], double in[3]) { + int i; + double tt[3]; + + /* Convert RGB' to RGB */ + for (i = 0; i < 3; i++) { + if (in[i] < (4.5 * 0.0181)) + tt[i] = in[i]/4.5; + else + tt[i] = pow((in[i] + 0.0993)/1.0993, 1.0/0.45); + } + + /* Y value */ + tt[0] = 0.2627 * tt[0] + 0.6780 * tt[1] + 0.0593 * tt[2]; + + /* Y' value */ + if (tt[0] < 0.0181) + tt[0] = tt[0] * 4.5; + else + tt[0] = 1.0993 * pow(tt[0], 0.45) - 0.0993; + + tt[1] = in[2] - tt[0]; + if (tt[1] <= 0.0) + tt[1] /= 1.9404; + else + tt[1] /= 1.5816; + + tt[2] = in[0] - tt[0]; + if (tt[2] <= 0.0) + tt[2] /= 1.7184; + else + tt[2] /= 0.9936; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + +/* Convert Rec2020 Constant liminance YPbPr into RGB' (== "full range YCbCr") */ +/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */ +/* [Inverse of above] */ +void icmRec2020_CL_YPbPr_2_RGBd(double out[3], double in[3]) { + int i; + double tin[3], tt[3]; + + /* Y' */ + tin[0] = in[0]; + + /* B' - Y' */ + if (in[1] <= 0.0) + tin[1] = 1.9404 * in[1]; + else + tin[1] = 1.5816 * in[1]; + + /* R' - Y' */ + if (in[2] <= 0.0) + tin[2] = 1.7184 * in[2]; + else + tin[2] = 0.9936 * in[2]; + + + /* R' */ + tt[0] = tin[2] + tin[0]; + + /* Y' */ + tt[1] = tin[0]; + + /* B' */ + tt[2] = tin[1] + tin[0]; + + /* Convert RYB' to RYB */ + for (i = 0; i < 3; i++) { + if (tt[i] < (4.5 * 0.0181)) + tin[i] = tt[i]/4.5; + else + tin[i] = pow((tt[i] + 0.0993)/1.0993, 1.0/0.45); + } + + /* G */ + tt[1] = (tin[1] - 0.2627 * tin[0] - 0.0593 * tin[2])/0.6780; + + /* G' */ + if (tt[1] < 0.0181) + tt[1] = tt[1] * 4.5; + else + tt[1] = 1.0993 * pow(tt[1], 0.45) - 0.0993; + + out[0] = tt[0]; + out[1] = tt[1]; + out[2] = tt[2]; +} + + +/* Convert Rec601/Rec709/Rec2020 YPbPr to YCbCr Video range. */ +/* input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, */ +/* output 16/255 .. 235/255, 16/255 .. 240/255, 16/255 .. 240/255 */ +void icmRecXXX_YPbPr_2_YCbCr(double out[3], double in[3]) { + out[0] = ((235.0 - 16.0) * in[0] + 16.0)/255.0; + out[1] = ((128.0 - 16.0) * 2.0 * in[1] + 128.0)/255.0; + out[2] = ((128.0 - 16.0) * 2.0 * in[2] + 128.0)/255.0; +} + +/* Convert Rec601/Rec709/Rec2020 Video YCbCr to YPbPr range. */ +/* input 16/255 .. 235/255, 16/255 .. 240/255, 16/255 .. 240/255 */ +/* output 0..1, -0.5 .. 0.5, -0.5 .. 0.5, */ +void icmRecXXX_YCbCr_2_YPbPr(double out[3], double in[3]) { + out[0] = (255.0 * in[0] - 16.0)/(235.0 - 16.0); + out[1] = (255.0 * in[1] - 128.0)/(2.0 * (128.0 - 16.0)); + out[2] = (255.0 * in[2] - 128.0)/(2.0 * (128.0 - 16.0)); +} + +/* Convert full range RGB to Video range 16..235 RGB */ +void icmRGB_2_VidRGB(double out[3], double in[3]) { + out[0] = ((235.0 - 16.0) * in[0] + 16.0)/255.0; + out[1] = ((235.0 - 16.0) * in[1] + 16.0)/255.0; + out[2] = ((235.0 - 16.0) * in[2] + 16.0)/255.0; +} + +/* Convert Video range 16..235 RGB to full range RGB */ +/* Return nz if outside RGB range */ +void icmVidRGB_2_RGB(double out[3], double in[3]) { + out[0] = (255.0 * in[0] - 16.0)/(235.0 - 16.0); + out[1] = (255.0 * in[1] - 16.0)/(235.0 - 16.0); + out[2] = (255.0 * in[2] - 16.0)/(235.0 - 16.0); +} + +/* =============================================================== */ +/* PS 3.14-2009, Digital Imaging and Communications in Medicine */ +/* (DICOM) Part 14: Grayscale Standard Display Function */ + +/* JND index value 1..1023 to L 0.05 .. 3993.404 cd/m^2 */ +static double icmDICOM_fwd_nl(double jnd) { + double a = -1.3011877; + double b = -2.5840191e-2; + double c = 8.0242636e-2; + double d = -1.0320229e-1; + double e = 1.3646699e-1; + double f = 2.8745620e-2; + double g = -2.5468404e-2; + double h = -3.1978977e-3; + double k = 1.2992634e-4; + double m = 1.3635334e-3; + double jj, num, den, rv; + + jj = jnd = log(jnd); + + num = a; + den = 1.0; + num += c * jj; + den += b * jj; + jj *= jnd; + num += e * jj; + den += d * jj; + jj *= jnd; + num += g * jj; + den += f * jj; + jj *= jnd; + num += m * jj; + den += h * jj; + jj *= jnd; + den += k * jj; + + rv = pow(10.0, num/den); + + return rv; +} + +/* JND index value 1..1023 to L 0.05 .. 3993.404 cd/m^2 */ +double icmDICOM_fwd(double jnd) { + + if (jnd < 0.5) + jnd = 0.5; + if (jnd > 1024.0) + jnd = 1024.0; + + return icmDICOM_fwd_nl(jnd); +} + +/* L 0.05 .. 3993.404 cd/m^2 to JND index value 1..1023 */ +/* This is not super accurate - typically to 0.03 .. 0.1 jne. */ +static double icmDICOM_bwd_apx(double L) { + double A = 71.498068; + double B = 94.593053; + double C = 41.912053; + double D = 9.8247004; + double E = 0.28175407; + double F = -1.1878455; + double G = -0.18014349; + double H = 0.14710899; + double I = -0.017046845; + double rv, LL; + + if (L < 0.049982) { /* == jnd 0.5 */ + return 0.5; + } + if (L > 4019.354716) /* == jnd 1024 */ + L = 4019.354716; + + LL = L = log10(L); + rv = A; + rv += B * LL; + LL *= L; + rv += C * LL; + LL *= L; + rv += D * LL; + LL *= L; + rv += E * LL; + LL *= L; + rv += F * LL; + LL *= L; + rv += G * LL; + LL *= L; + rv += H * LL; + LL *= L; + rv += I * LL; + + return rv; +} + +/* L 0.05 .. 3993.404 cd/m^2 to JND index value 1..1023 */ +/* Polish the aproximate solution twice using Newton's itteration */ +double icmDICOM_bwd(double L) { + double rv, Lc, prv, pLc, de; + int i; + + if (L < 0.045848) /* == jnd 0.5 */ + L = 0.045848; + if (L > 4019.354716) /* == jnd 1024 */ + L = 4019.354716; + + /* Approx solution */ + rv = icmDICOM_bwd_apx(L); + + /* Compute aprox derivative */ + Lc = icmDICOM_fwd_nl(rv); + + prv = rv + 0.01; + pLc = icmDICOM_fwd_nl(prv); + + do { + de = (rv - prv)/(Lc - pLc); + prv = rv; + rv -= (Lc - L) * de; + pLc = Lc; + Lc = icmDICOM_fwd_nl(rv); + } while (fabs(Lc - L) > 1e-8); + + return rv; +} + + +/* =============================================================== */ /* Object for computing RFC 1321 MD5 checksums. */ /* Derived from Colin Plumb's 1993 public domain code. */ @@ -16174,25 +17067,25 @@ double *in /* Vector of input values */ icmLut *lut = p->lut; double temp[MAX_CHAN]; - DBGLL(("icmLuLut_lookup: in = %s\n", icmPdv(p->inputChan, in))); + DBGLL(("icmLuLut_lookup: in = %s\n", icmPdv(p->lut->inputChan, in))); rv |= p->in_abs(p,temp,in); /* Possible absolute conversion */ - DBGLL(("icmLuLut_lookup: in_abs = %s\n", icmPdv(p->inputChan, temp))); + DBGLL(("icmLuLut_lookup: in_abs = %s\n", icmPdv(p->lut->inputChan, temp))); if (p->usematrix) { rv |= lut->lookup_matrix(lut,temp,temp);/* If XYZ, multiply by non-unity matrix */ - DBGLL(("icmLuLut_lookup: matrix = %s\n", icmPdv(p->inputChan, temp))); + DBGLL(("icmLuLut_lookup: matrix = %s\n", icmPdv(p->lut->inputChan, temp))); } p->in_normf(temp, temp); /* Normalize for input color space */ - DBGLL(("icmLuLut_lookup: norm = %s\n", icmPdv(p->inputChan, temp))); + DBGLL(("icmLuLut_lookup: norm = %s\n", icmPdv(p->lut->inputChan, temp))); rv |= lut->lookup_input(lut,temp,temp); /* Lookup though input tables */ - DBGLL(("icmLuLut_lookup: input = %s\n", icmPdv(p->inputChan, temp))); + DBGLL(("icmLuLut_lookup: input = %s\n", icmPdv(p->lut->inputChan, temp))); rv |= p->lookup_clut(lut,out,temp); /* Lookup though clut tables */ - DBGLL(("icmLuLut_lookup: clut = %s\n", icmPdv(p->outputChan, out))); + DBGLL(("icmLuLut_lookup: clut = %s\n", icmPdv(p->lut->outputChan, out))); rv |= lut->lookup_output(lut,out,out); /* Lookup though output tables */ - DBGLL(("icmLuLut_lookup: output = %s\n", icmPdv(p->outputChan, out))); + DBGLL(("icmLuLut_lookup: output = %s\n", icmPdv(p->lut->outputChan, out))); p->out_denormf(out,out); /* Normalize for output color space */ - DBGLL(("icmLuLut_lookup: denorm = %s\n", icmPdv(p->outputChan, out))); + DBGLL(("icmLuLut_lookup: denorm = %s\n", icmPdv(p->lut->outputChan, out))); rv |= p->out_abs(p,out,out); /* Possible absolute conversion */ - DBGLL(("icmLuLut_lookup: out_abse = %s\n", icmPdv(p->outputChan, out))); + DBGLL(("icmLuLut_lookup: out_abse = %s\n", icmPdv(p->lut->outputChan, out))); return rv; } @@ -16896,8 +17789,10 @@ icc_new_icmLuLut( if (use_sx) { p->lookup_clut = p->lut->lookup_clut_sx; + p->lut->tune_value = icmLut_tune_value_sx; } else { p->lookup_clut = p->lut->lookup_clut_nl; + p->lut->tune_value = icmLut_tune_value_nl; } } #else /* Development code */ @@ -17670,7 +18565,7 @@ icmAlloc *al /* Memory allocator */ if ((p = (icc *) al->calloc(al, 1,sizeof(icc))) == NULL) { return NULL; } - p->ver = 0; /* default is V2 profile */ + p->ver = icmVersionDefault; /* default is V2.2.0 profile */ p->al = al; /* Heap allocator */ |